Context Navigation

source: trunk/GSASIIpwd.py @ 4427

Last change on this file since 4427 was 4427, checked in by vondreele, 3 years ago
fullrmc fixes
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 163.7 KB

Line
1	#/usr/bin/env python
2	# -- coding: utf-8 --
3	'''
4	GSASII powder calculation module
5	==================================
6
7	'''
8	########### SVN repository information ###################
9	# $Date: 2020-05-21 00:33:57 +0000 (Thu, 21 May 2020) $
10	# $Author: vondreele $
11	# $Revision: 4427 $
12	# $URL: trunk/GSASIIpwd.py $
13	# $Id: GSASIIpwd.py 4427 2020-05-21 00:33:57Z vondreele $
14	########### SVN repository information ###################
15	from __future__ import division, print_function
16	import sys
17	import math
18	import time
19	import os
20	import os.path
21	import subprocess as subp
22	import copy
23
24	import numpy as np
25	import numpy.linalg as nl
26	import numpy.ma as ma
27	import random as rand
28	import numpy.fft as fft
29	import scipy.interpolate as si
30	import scipy.stats as st
31	import scipy.optimize as so
32	import scipy.special as sp
33
34	import GSASIIpath
35	GSASIIpath.SetVersionNumber("$Revision: 4427 $")
36	import GSASIIlattice as G2lat
37	import GSASIIspc as G2spc
38	import GSASIIElem as G2elem
39	import GSASIImath as G2mth
40	try:
41	import pypowder as pyd
42	except ImportError:
43	print ('pypowder is not available - profile calcs. not allowed')
44	try:
45	import pydiffax as pyx
46	except ImportError:
47	print ('pydiffax is not available for this platform')
48	import GSASIIfiles as G2fil
49
50
51	# trig functions in degrees
52	tand = lambda x: math.tan(x*math.pi/180.)
53	atand = lambda x: 180.*math.atan(x)/math.pi
54	atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
55	cosd = lambda x: math.cos(x*math.pi/180.)
56	acosd = lambda x: 180.*math.acos(x)/math.pi
57	rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
58	#numpy versions
59	npsind = lambda x: np.sin(x*np.pi/180.)
60	npasind = lambda x: 180.*np.arcsin(x)/math.pi
61	npcosd = lambda x: np.cos(x*math.pi/180.)
62	npacosd = lambda x: 180.*np.arccos(x)/math.pi
63	nptand = lambda x: np.tan(x*math.pi/180.)
64	npatand = lambda x: 180.*np.arctan(x)/np.pi
65	npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
66	npT2stl = lambda tth, wave: 2.0npsind(tth/2.0)/wave #=d
67	npT2q = lambda tth,wave: 2.0np.pinpT2stl(tth,wave) #=2pid
68	ateln2 = 8.0*math.log(2.0)
69	sateln2 = np.sqrt(ateln2)
70	nxs = np.newaxis
71
72	################################################################################
73	#### Powder utilities
74	################################################################################
75
76	def PhaseWtSum(G2frame,histo):
77	'''
78	Calculate sum of phase mass*phase fraction for PWDR data (exclude magnetic phases)
79
80	:param G2frame: GSASII main frame structure
81	:param str histo: histogram name
82	:returns: sum(scale*mass) for phases in histo
83	'''
84	Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
85	wtSum = 0.0
86	for phase in Phases:
87	if Phases[phase]['General']['Type'] != 'magnetic':
88	if histo in Phases[phase]['Histograms']:
89	if not Phases[phase]['Histograms'][histo]['Use']: continue
90	mass = Phases[phase]['General']['Mass']
91	phFr = Phases[phase]['Histograms'][histo]['Scale'][0]
92	wtSum += mass*phFr
93	return wtSum
94
95	################################################################################
96	#### GSASII pwdr & pdf calculation routines
97	################################################################################
98
99	def Transmission(Geometry,Abs,Diam):
100	'''
101	Calculate sample transmission
102
103	:param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
104	:param float Abs: absorption coeff in cm-1
105	:param float Diam: sample thickness/diameter in mm
106	'''
107	if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
108	MuR = Abs*Diam/20.0
109	if MuR <= 3.0:
110	T0 = 16/(3.*math.pi)
111	T1 = -0.045780
112	T2 = -0.02489
113	T3 = 0.003045
114	T = -T0MuR-T1MuR*2-T2MuR*3-T3MuR**4
115	if T < -20.:
116	return 2.06e-9
117	else:
118	return math.exp(T)
119	else:
120	T1 = 1.433902
121	T2 = 0.013869+0.337894
122	T3 = 1.933433+1.163198
123	T4 = 0.044365-0.04259
124	T = (T1-T4)/(1.0+T2(MuR-3.0))*T3+T4
125	return T/100.
126	elif 'plate' in Geometry:
127	MuR = Abs*Diam/10.
128	return math.exp(-MuR)
129	elif 'Bragg' in Geometry:
130	return 0.0
131
132	def SurfaceRough(SRA,SRB,Tth):
133	''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
134	:param float SRA: Suortti surface roughness parameter
135	:param float SRB: Suortti surface roughness parameter
136	:param float Tth: 2-theta(deg) - can be numpy array
137
138	'''
139	sth = npsind(Tth/2.)
140	T1 = np.exp(-SRB/sth)
141	T2 = SRA+(1.-SRA)*np.exp(-SRB)
142	return (SRA+(1.-SRA)*T1)/T2
143
144	def SurfaceRoughDerv(SRA,SRB,Tth):
145	''' Suortti surface roughness correction derivatives
146	:param float SRA: Suortti surface roughness parameter (dimensionless)
147	:param float SRB: Suortti surface roughness parameter (dimensionless)
148	:param float Tth: 2-theta(deg) - can be numpy array
149	:return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
150	'''
151	sth = npsind(Tth/2.)
152	T1 = np.exp(-SRB/sth)
153	T2 = SRA+(1.-SRA)*np.exp(-SRB)
154	Trans = (SRA+(1.-SRA)*T1)/T2
155	dydSRA = ((1.-T1)T2-(1.-np.exp(-SRB))Trans)/T2**2
156	dydSRB = ((SRA-1.)T1T2/sth-Trans(SRA-T2))/T2*2
157	return [dydSRA,dydSRB]
158
159	def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
160	'''Calculate sample absorption
161	:param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
162	:param float MuR: absorption coeff * sample thickness/2 or radius
163	:param Tth: 2-theta scattering angle - can be numpy array
164	:param float Phi: flat plate tilt angle - future
165	:param float Psi: flat plate tilt axis - future
166	'''
167
168	def muRunder3(MuR,Sth2):
169	T0 = 16.0/(3.*np.pi)
170	T1 = (25.99978-0.01911Sth20.25)np.exp(-0.024551*Sth2)+ \
171	0.109561*np.sqrt(Sth2)-26.04556
172	T2 = -0.02489-0.39499Sth2+1.219077Sth2**1.5- \
173	1.31268Sth22+0.871081Sth2*2.5-0.2327Sth2**3
174	T3 = 0.003045+0.018167Sth2-0.03305Sth2**2
175	Trns = -T0MuR-T1MuR*2-T2MuR*3-T3MuR**4
176	return np.exp(Trns)
177
178	def muRover3(MuR,Sth2):
179	T1 = 1.433902+11.07504Sth2-8.77629Sth2*Sth2+ \
180	10.02088Sth23-3.36778Sth2**4
181	T2 = (0.013869-0.01249Sth2)np.exp(3.27094*Sth2)+ \
182	(0.337894+13.77317Sth2)/(1.0+11.53544Sth2)**1.555039
183	T3 = 1.933433/(1.0+23.12967Sth2)*1.686715- \
184	0.13576*np.sqrt(Sth2)+1.163198
185	T4 = 0.044365-0.04259/(1.0+0.41051Sth2)*148.4202
186	Trns = (T1-T4)/(1.0+T2(MuR-3.0))*T3+T4
187	return Trns/100.
188
189	Sth2 = npsind(Tth/2.0)**2
190	if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
191	if 'array' in str(type(MuR)):
192	MuRSTh2 = np.vstack((MuR,Sth2))
193	AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1]))
194	return AbsCr
195	else:
196	if MuR <= 3.0:
197	return muRunder3(MuR,Sth2)
198	else:
199	return muRover3(MuR,Sth2)
200	elif 'Bragg' in Geometry:
201	return 1.0
202	elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
203	# and only defined for 2theta < 90
204	MuT = 2.*MuR
205	T1 = np.exp(-MuT)
206	T2 = np.exp(-MuT/npcosd(Tth))
207	Tb = MuT-MuT/npcosd(Tth)
208	return (T2-T1)/Tb
209	elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
210	MuT = 2.*MuR
211	cth = npcosd(Tth/2.0)
212	return np.exp(-MuT/cth)/cth
213
214	def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
215	'needs a doc string'
216	dA = 0.001
217	AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
218	if MuR:
219	AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
220	return (AbsP-AbsM)/(2.0*dA)
221	else:
222	return (AbsP-1.)/dA
223
224	def Polarization(Pola,Tth,Azm=0.0):
225	""" Calculate angle dependent x-ray polarization correction (not scaled correctly!)
226
227	:param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
228	:param Azm: azimuthal angle e.g. 0.0 in plane of polarization
229	:param Tth: 2-theta scattering angle - can be numpy array
230	which (if either) of these is "right"?
231	:return: (pola, dpdPola)
232	* pola = ((1-Pola)npcosd(Azm)2+Polanpsind(Azm)*2)npcosd(Tth)**2+ \
233	(1-Pola)npsind(Azm)2+Polanpcosd(Azm)**2
234	* dpdPola: derivative needed for least squares
235
236	"""
237	cazm = npcosd(Azm)
238	sazm = npsind(Azm)
239	pola = ((1.0-Pola)cazm2+Polanpsind(Azm)*2)npcosd(Tth)**2+ \
240	(1.0-Pola)sazm2+Polacazm**2
241	dpdPola = -npsind(Tth)*2(sazm2-cazm2)
242	return pola,dpdPola
243
244	def Oblique(ObCoeff,Tth):
245	'currently assumes detector is normal to beam'
246	if ObCoeff:
247	return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
248	else:
249	return 1.0
250
251	def Ruland(RulCoff,wave,Q,Compton):
252	'needs a doc string'
253	C = 2.9978e8
254	D = 1.5e-3
255	hmc = 0.024262734687
256	sinth2 = (Qwave/(4.0np.pi))**2
257	dlam = (wave*2)Compton*Q/C
258	dlam_c = 2.0hmcsinth2-Dwave*2
259	return 1.0/((1.0+dlam/RulCoff)(1.0+(np.pidlam_c/(dlam+RulCoff))**2))
260
261	def LorchWeight(Q):
262	'needs a doc string'
263	return np.sin(np.pi(Q[-1]-Q)/(2.0Q[-1]))
264
265	def GetAsfMean(ElList,Sthl2):
266	'''Calculate various scattering factor terms for PDF calcs
267
268	:param dict ElList: element dictionary contains scattering factor coefficients, etc.
269	:param np.array Sthl2: numpy array of sin theta/lambda squared values
270	:returns: mean(f^2), mean(f)^2, mean(compton)
271	'''
272	sumNoAtoms = 0.0
273	FF = np.zeros_like(Sthl2)
274	FF2 = np.zeros_like(Sthl2)
275	CF = np.zeros_like(Sthl2)
276	for El in ElList:
277	sumNoAtoms += ElList[El]['FormulaNo']
278	for El in ElList:
279	el = ElList[El]
280	ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])2+el['fpp']2
281	cf = G2elem.ComptonFac(el,Sthl2)
282	FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
283	FF2 += ff2*el['FormulaNo']/sumNoAtoms
284	CF += cf*el['FormulaNo']/sumNoAtoms
285	return FF2,FF**2,CF
286
287	def GetNumDensity(ElList,Vol):
288	'needs a doc string'
289	sumNoAtoms = 0.0
290	for El in ElList:
291	sumNoAtoms += ElList[El]['FormulaNo']
292	return sumNoAtoms/Vol
293
294	def CalcPDF(data,inst,limits,xydata):
295	'''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
296	dict xydata; results are placed in xydata.
297	Calculation parameters are found in dicts data and inst and list limits.
298	The return value is at present an empty list.
299	'''
300	auxPlot = []
301	if 'T' in inst['Type'][0]:
302	Ibeg = 0
303	Ifin = len(xydata['Sample'][1][0])
304	else:
305	Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0])
306	Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1
307	#subtract backgrounds - if any & use PWDR limits
308	IofQ = copy.deepcopy(xydata['Sample'])
309	IofQ[1] = np.array(IofQ[1])[:,Ibeg:Ifin]
310	if data['Sample Bkg.']['Name']:
311	IofQ[1][1] += xydata['Sample Bkg.'][1][1][Ibeg:Ifin]*data['Sample Bkg.']['Mult']
312	if data['Container']['Name']:
313	xycontainer = xydata['Container'][1][1]*data['Container']['Mult']
314	if data['Container Bkg.']['Name']:
315	xycontainer += xydata['Container Bkg.'][1][1][Ibeg:Ifin]*data['Container Bkg.']['Mult']
316	IofQ[1][1] += xycontainer[Ibeg:Ifin]
317	data['IofQmin'] = IofQ[1][1][-1]
318	IofQ[1][1] -= data.get('Flat Bkg',0.)
319	#get element data & absorption coeff.
320	ElList = data['ElList']
321	Tth = IofQ[1][0] #2-theta or TOF!
322	if 'X' in inst['Type'][0]:
323	Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
324	#Apply angle dependent corrections
325	MuR = Abs*data['Diam']/20.0
326	IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
327	IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
328	if data['DetType'] == 'Image plate':
329	IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
330	elif 'T' in inst['Type'][0]: #neutron TOF normalized data - needs wavelength dependent absorption
331	wave = 2.G2lat.TOF2dsp(inst,IofQ[1][0])npsind(inst['2-theta'][1]/2.)
332	Els = ElList.keys()
333	Isotope = {El:'Nat. abund.' for El in Els}
334	GD = {'AtomTypes':ElList,'Isotope':Isotope}
335	BLtables = G2elem.GetBLtable(GD)
336	FP,FPP = G2elem.BlenResTOF(Els,BLtables,wave)
337	Abs = np.zeros(len(wave))
338	for iel,El in enumerate(Els):
339	BL = BLtables[El][1]
340	SA = BL['SA']wave/1.798197+4.0np.piFPP[iel]*2 #+BL['SL'][1]?
341	SA *= ElList[El]['FormulaNo']/data['Form Vol']
342	Abs += SA
343	MuR = Abs*data['Diam']/2.
344	IofQ[1][1] /= Absorb(data['Geometry'],MuR,inst['2-theta'][1]*np.ones(len(wave)))
345	XY = IofQ[1]
346	#convert to Q
347	# nQpoints = len(XY[0]) #points for Q interpolation
348	nQpoints = 5000
349	if 'C' in inst['Type'][0]:
350	wave = G2mth.getWave(inst)
351	minQ = npT2q(Tth[0],wave)
352	maxQ = npT2q(Tth[-1],wave)
353	Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
354	dq = Qpoints[1]-Qpoints[0]
355	XY[0] = npT2q(XY[0],wave)
356	Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0]) #interpolate I(Q)
357	elif 'T' in inst['Type'][0]:
358	difC = inst['difC'][1]
359	minQ = 2.np.pidifC/Tth[-1]
360	maxQ = 2.np.pidifC/Tth[0]
361	Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
362	dq = Qpoints[1]-Qpoints[0]
363	XY[0] = 2.np.pidifC/XY[0]
364	Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][-1]) #interpolate I(Q)
365	Qdata -= np.min(Qdata)*data['BackRatio']
366
367	qLimits = data['QScaleLim']
368	maxQ = np.searchsorted(Qpoints,min(Qpoints[-1],qLimits[1]))+1
369	minQ = np.searchsorted(Qpoints,min(qLimits[0],0.90*Qpoints[-1]))
370	qLimits = [Qpoints[minQ],Qpoints[maxQ-1]]
371	newdata = []
372	if len(IofQ) < 3:
373	xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],'']
374	else:
375	xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
376	for item in xydata['IofQ'][1]:
377	newdata.append(item[:maxQ])
378	xydata['IofQ'][1] = newdata
379
380	xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
381	if 'XC' in inst['Type'][0]:
382	FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0np.pi))*2) #these are <f^2>,<f>^2,Cf
383	else: #TOF
384	CF = np.zeros(len(xydata['SofQ'][1][0]))
385	FFSq = np.ones(len(xydata['SofQ'][1][0]))
386	SqFF = np.ones(len(xydata['SofQ'][1][0]))
387	Q = xydata['SofQ'][1][0]
388	# auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
389	if 'XC' in inst['Type'][0]:
390	ruland = Ruland(data['Ruland'],wave,Q,CF)
391	# auxPlot.append([Q,ruland,'Ruland'])
392	CF *= ruland
393	# auxPlot.append([Q,CF,'CF-Corr'])
394	scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
395	xydata['SofQ'][1][1] *= scale
396	if 'XC' in inst['Type'][0]:
397	xydata['SofQ'][1][1] -= CF
398	xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
399	scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
400	xydata['SofQ'][1][1] *= scale
401	xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
402	xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
403	if data['Lorch']:
404	xydata['FofQ'][1][1] *= LorchWeight(Q)
405	xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
406	xydata['gofr'] = copy.deepcopy(xydata['FofQ'])
407	nR = len(xydata['GofR'][1][1])
408	Rmax = GSASIIpath.GetConfigValue('PDF_Rmax',100.)
409	mul = int(round(2.np.pinR/(Rmax*qLimits[1])))
410	# mul = int(round(2.np.pinR/(data.get('Rmax',100.)*qLimits[1])))
411	R = 2.np.pinp.linspace(0,nR,nR,endpoint=True)/(mul*qLimits[1])
412	xydata['GofR'][1][0] = R
413	xydata['gofr'][1][0] = R
414	GR = -dqnp.imag(fft.fft(xydata['FofQ'][1][1],mulnR)[:nR])
415	xydata['GofR'][1][1] = GR
416	gr = GR/(np.pi*R)
417	xydata['gofr'][1][1] = gr
418	numbDen = 0.
419	if 'ElList' in data:
420	numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
421	if data.get('noRing',True):
422	Rmin = data['Rmin']
423	xydata['gofr'][1][1] = np.where(R<Rmin,-4.*numbDen,xydata['gofr'][1][1])
424	xydata['GofR'][1][1] = np.where(R<Rmin,-4.Rnp.pi*numbDen,xydata['GofR'][1][1])
425	return auxPlot
426
427	def PDFPeakFit(peaks,data):
428	rs2pi = 1./np.sqrt(2*np.pi)
429
430	def MakeParms(peaks):
431	varyList = []
432	parmDict = {'slope':peaks['Background'][1][1]}
433	if peaks['Background'][2]:
434	varyList.append('slope')
435	for i,peak in enumerate(peaks['Peaks']):
436	parmDict['PDFpos;'+str(i)] = peak[0]
437	parmDict['PDFmag;'+str(i)] = peak[1]
438	parmDict['PDFsig;'+str(i)] = peak[2]
439	if 'P' in peak[3]:
440	varyList.append('PDFpos;'+str(i))
441	if 'M' in peak[3]:
442	varyList.append('PDFmag;'+str(i))
443	if 'S' in peak[3]:
444	varyList.append('PDFsig;'+str(i))
445	return parmDict,varyList
446
447	def SetParms(peaks,parmDict,varyList):
448	if 'slope' in varyList:
449	peaks['Background'][1][1] = parmDict['slope']
450	for i,peak in enumerate(peaks['Peaks']):
451	if 'PDFpos;'+str(i) in varyList:
452	peak[0] = parmDict['PDFpos;'+str(i)]
453	if 'PDFmag;'+str(i) in varyList:
454	peak[1] = parmDict['PDFmag;'+str(i)]
455	if 'PDFsig;'+str(i) in varyList:
456	peak[2] = parmDict['PDFsig;'+str(i)]
457
458
459	def CalcPDFpeaks(parmdict,Xdata):
460	Z = parmDict['slope']*Xdata
461	ipeak = 0
462	while True:
463	try:
464	pos = parmdict['PDFpos;'+str(ipeak)]
465	mag = parmdict['PDFmag;'+str(ipeak)]
466	wid = parmdict['PDFsig;'+str(ipeak)]
467	wid2 = 2.wid*2
468	Z += magrs2pinp.exp(-(Xdata-pos)**2/wid2)/wid
469	ipeak += 1
470	except KeyError: #no more peaks to process
471	return Z
472
473	def errPDFProfile(values,xdata,ydata,parmdict,varylist):
474	parmdict.update(zip(varylist,values))
475	M = CalcPDFpeaks(parmdict,xdata)-ydata
476	return M
477
478	newpeaks = copy.copy(peaks)
479	iBeg = np.searchsorted(data[1][0],newpeaks['Limits'][0])
480	iFin = np.searchsorted(data[1][0],newpeaks['Limits'][1])+1
481	X = data[1][0][iBeg:iFin]
482	Y = data[1][1][iBeg:iFin]
483	parmDict,varyList = MakeParms(peaks)
484	if not len(varyList):
485	G2fil.G2Print (' Nothing varied')
486	return newpeaks,None,None,None,None,None
487
488	Rvals = {}
489	values = np.array(Dict2Values(parmDict, varyList))
490	result = so.leastsq(errPDFProfile,values,full_output=True,ftol=0.0001,
491	args=(X,Y,parmDict,varyList))
492	chisq = np.sum(result[2]['fvec']**2)
493	Values2Dict(parmDict, varyList, result[0])
494	SetParms(peaks,parmDict,varyList)
495	Rvals['Rwp'] = np.sqrt(chisq/np.sum(Y*2))100. #to %
496	chisq = np.sum(result[2]['fvec']**2)/(len(X)-len(values)) #reduced chi^2 = M/(Nobs-Nvar)
497	sigList = list(np.sqrt(chisq*np.diag(result[1])))
498	Z = CalcPDFpeaks(parmDict,X)
499	newpeaks['calc'] = [X,Z]
500	return newpeaks,result[0],varyList,sigList,parmDict,Rvals
501
502	def MakeRDF(RDFcontrols,background,inst,pwddata):
503	import scipy.signal as signal
504	auxPlot = []
505	if 'C' in inst['Type'][0]:
506	Tth = pwddata[0]
507	wave = G2mth.getWave(inst)
508	minQ = npT2q(Tth[0],wave)
509	maxQ = npT2q(Tth[-1],wave)
510	powQ = npT2q(Tth,wave)
511	elif 'T' in inst['Type'][0]:
512	TOF = pwddata[0]
513	difC = inst['difC'][1]
514	minQ = 2.np.pidifC/TOF[-1]
515	maxQ = 2.np.pidifC/TOF[0]
516	powQ = 2.np.pidifC/TOF
517	piDQ = np.pi/(maxQ-minQ)
518	Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
519	if RDFcontrols['UseObsCalc'] == 'obs-calc':
520	Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
521	elif RDFcontrols['UseObsCalc'] == 'obs-back':
522	Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
523	elif RDFcontrols['UseObsCalc'] == 'calc-back':
524	Qdata = si.griddata(powQ,pwddata[3]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
525	Qdata = np.sin((Qpoints-minQ)piDQ)/piDQ
526	Qdata = 0.5np.sqrt(Qpoints) #Qbin normalization
527	# GSASIIpath.IPyBreak()
528	dq = Qpoints[1]-Qpoints[0]
529	nR = len(Qdata)
530	R = 0.5np.pinp.linspace(0,nR,nR)/(4.*maxQ)
531	iFin = np.searchsorted(R,RDFcontrols['maxR'])+1
532	bBut,aBut = signal.butter(4,0.01)
533	Qsmooth = signal.filtfilt(bBut,aBut,Qdata)
534	# auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']])
535	# auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']])
536	DofR = dqnp.imag(fft.fft(Qsmooth,16nR)[:nR])
537	# DofR = dqnp.imag(ft.fft(Qsmooth,16nR)[:nR])
538	auxPlot.append([R[:iFin],DofR[:iFin],'D(R) for '+RDFcontrols['UseObsCalc']])
539	return auxPlot
540
541	# PDF optimization =============================================================
542	def OptimizePDF(data,xydata,limits,inst,showFit=True,maxCycles=5):
543	import scipy.optimize as opt
544	numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
545	Min,Init,Done = SetupPDFEval(data,xydata,limits,inst,numbDen)
546	xstart = Init()
547	bakMul = data['Sample Bkg.']['Mult']
548	if showFit:
549	rms = Min(xstart)
550	G2fil.G2Print(' Optimizing corrections to improve G(r) at low r')
551	if data['Sample Bkg.'].get('Refine',False):
552	# data['Flat Bkg'] = 0.
553	G2fil.G2Print(' start: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f})'.format(
554	data['Ruland'],data['Sample Bkg.']['Mult'],rms))
555	else:
556	G2fil.G2Print(' start: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} (RMS:{:.4f})'.format(
557	data['Flat Bkg'],data['BackRatio'],data['Ruland'],rms))
558	if data['Sample Bkg.'].get('Refine',False):
559	res = opt.minimize(Min,xstart,bounds=([0.01,1],[1.2bakMul,0.8bakMul]),
560	method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
561	else:
562	res = opt.minimize(Min,xstart,bounds=([0,None],[0,1],[0.01,1]),
563	method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
564	Done(res['x'])
565	if showFit:
566	if res['success']:
567	msg = 'Converged'
568	else:
569	msg = 'Not Converged'
570	if data['Sample Bkg.'].get('Refine',False):
571	G2fil.G2Print(' end: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f}) * {} *\n'.format(
572	data['Ruland'],data['Sample Bkg.']['Mult'],res['fun'],msg))
573	else:
574	G2fil.G2Print(' end: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f}) * {} *\n'.format(
575	data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun'],msg))
576	return res
577
578	def SetupPDFEval(data,xydata,limits,inst,numbDen):
579	Data = copy.deepcopy(data)
580	BkgMax = 1.
581	def EvalLowPDF(arg):
582	'''Objective routine -- evaluates the RMS deviations in G(r)
583	from -4(pi)#densityr for for r<Rmin
584	arguments are ['Flat Bkg','BackRatio','Ruland'] scaled so that
585	the min & max values are between 0 and 1.
586	'''
587	if Data['Sample Bkg.'].get('Refine',False):
588	R,S = arg
589	Data['Sample Bkg.']['Mult'] = S
590	else:
591	F,B,R = arg
592	Data['Flat Bkg'] = F*BkgMax
593	Data['BackRatio'] = B
594	Data['Ruland'] = R/10.
595	CalcPDF(Data,inst,limits,xydata)
596	# test low r computation
597	g = xydata['GofR'][1][1]
598	r = xydata['GofR'][1][0]
599	g0 = g[r < Data['Rmin']] + 4np.pir[r < Data['Rmin']]*numbDen
600	M = sum(g0**2)/len(g0)
601	return M
602	def GetCurrentVals():
603	'''Get the current ['Flat Bkg','BackRatio','Ruland'] with scaling
604	'''
605	if data['Sample Bkg.'].get('Refine',False):
606	return [max(10*data['Ruland'],.05),data['Sample']['Mult']]
607	try:
608	F = data['Flat Bkg']/BkgMax
609	except:
610	F = 0
611	return [F,data['BackRatio'],max(10*data['Ruland'],.05)]
612	def SetFinalVals(arg):
613	'''Set the 'Flat Bkg', 'BackRatio' & 'Ruland' values from the
614	scaled, refined values and plot corrected region of G(r)
615	'''
616	if data['Sample Bkg.'].get('Refine',False):
617	R,S = arg
618	data['Sample Bkg.']['Mult'] = S
619	else:
620	F,B,R = arg
621	data['Flat Bkg'] = F*BkgMax
622	data['BackRatio'] = B
623	data['Ruland'] = R/10.
624	CalcPDF(data,inst,limits,xydata)
625	EvalLowPDF(GetCurrentVals())
626	BkgMax = max(xydata['IofQ'][1][1])/50.
627	return EvalLowPDF,GetCurrentVals,SetFinalVals
628
629	################################################################################
630	#### GSASII peak fitting routines: Finger, Cox & Jephcoat model
631	################################################################################
632
633	def factorize(num):
634	''' Provide prime number factors for integer num
635	:returns: dictionary of prime factors (keys) & power for each (data)
636	'''
637	factors = {}
638	orig = num
639
640	# we take advantage of the fact that (i +1)2 = i2 + 2*i +1
641	i, sqi = 2, 4
642	while sqi <= num:
643	while not num%i:
644	num /= i
645	factors[i] = factors.get(i, 0) + 1
646
647	sqi += 2*i + 1
648	i += 1
649
650	if num != 1 and num != orig:
651	factors[num] = factors.get(num, 0) + 1
652
653	if factors:
654	return factors
655	else:
656	return {num:1} #a prime number!
657
658	def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
659	''' Provide list of optimal data sizes for FFT calculations
660
661	:param int nmin: minimum data size >= 1
662	:param int nmax: maximum data size > nmin
663	:param int thresh: maximum prime factor allowed
664	:Returns: list of data sizes where the maximum prime factor is < thresh
665	'''
666	plist = []
667	nmin = max(1,nmin)
668	nmax = max(nmin+1,nmax)
669	for p in range(nmin,nmax):
670	if max(list(factorize(p).keys())) < thresh:
671	plist.append(p)
672	return plist
673
674	np.seterr(divide='ignore')
675
676	# Normal distribution
677
678	# loc = mu, scale = std
679	_norm_pdf_C = 1./math.sqrt(2*math.pi)
680	class norm_gen(st.rv_continuous):
681	'needs a doc string'
682
683	def pdf(self,x,args,*kwds):
684	loc,scale=kwds['loc'],kwds['scale']
685	x = (x-loc)/scale
686	return np.exp(-x*2/2.0) _norm_pdf_C / scale
687
688	norm = norm_gen(name='norm',longname='A normal',extradoc="""
689
690	Normal distribution
691
692	The location (loc) keyword specifies the mean.
693	The scale (scale) keyword specifies the standard deviation.
694
695	normal.pdf(x) = exp(-x*2/2)/sqrt(2pi)
696	""")
697
698	## Cauchy
699
700	# median = loc
701
702	class cauchy_gen(st.rv_continuous):
703	'needs a doc string'
704
705	def pdf(self,x,args,*kwds):
706	loc,scale=kwds['loc'],kwds['scale']
707	x = (x-loc)/scale
708	return 1.0/np.pi/(1.0+x*x) / scale
709
710	cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
711
712	Cauchy distribution
713
714	cauchy.pdf(x) = 1/(pi(1+x*2))
715
716	This is the t distribution with one degree of freedom.
717	""")
718
719
720	#GSASII peak fitting routine: Finger, Cox & Jephcoat model
721
722
723	class fcjde_gen(st.rv_continuous):
724	"""
725	Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
726	Ref: J. Appl. Cryst. (1994) 27, 892-900.
727
728	:param x: array -1 to 1
729	:param t: 2-theta position of peak
730	:param s: sum(S/L,H/L); S: sample height, H: detector opening,
731	L: sample to detector opening distance
732	:param dx: 2-theta step size in deg
733
734	:returns: for fcj.pdf
735
736	* T = x*dx+t
737	* s = S/L+H/L
738	* if x < 0::
739
740	fcj.pdf = [1/sqrt({cos(T)2/cos(t)2}-1) - 1/s]/\|cos(T)\|
741
742	* if x >= 0: fcj.pdf = 0
743	"""
744	def _pdf(self,x,t,s,dx):
745	T = dx*x+t
746	ax2 = abs(npcosd(T))
747	ax = ax2**2
748	bx = npcosd(t)**2
749	bx = np.where(ax>bx,bx,ax)
750	fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
751	fx = np.where(fx > 0.,fx,0.0)
752	return fx
753
754	def pdf(self,x,args,*kwds):
755	loc=kwds['loc']
756	return self._pdf(x-loc,*args)
757
758	fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
759
760	def getWidthsCW(pos,sig,gam,shl):
761	'''Compute the peak widths used for computing the range of a peak
762	for constant wavelength data. On low-angle side, 50 FWHM are used,
763	on high-angle side 75 are used, low angle side extended for axial divergence
764	(for peaks above 90 deg, these are reversed.)
765	'''
766	widths = [np.sqrt(sig)/100.,gam/100.]
767	fwhm = 2.355*widths[0]+widths[1]
768	fmin = 50.(fwhm+shlabs(npcosd(pos)))
769	fmax = 75.0*fwhm
770	if pos > 90:
771	fmin,fmax = [fmax,fmin]
772	return widths,fmin,fmax
773
774	def getWidthsTOF(pos,alp,bet,sig,gam):
775	'''Compute the peak widths used for computing the range of a peak
776	for constant wavelength data. 50 FWHM are used on both sides each
777	extended by exponential coeff.
778	'''
779	widths = [np.sqrt(sig),gam]
780	fwhm = 2.355widths[0]+2.widths[1]
781	fmin = 50.fwhm(1.+1./alp)
782	fmax = 50.fwhm(1.+1./bet)
783	return widths,fmin,fmax
784
785	def getFWHM(pos,Inst):
786	'''Compute total FWHM from Thompson, Cox & Hastings (1987) , J. Appl. Cryst. 20, 79-83
787	via getgamFW(g,s).
788
789	:param pos: float peak position in deg 2-theta or tof in musec
790	:param Inst: dict instrument parameters
791
792	:returns float: total FWHM of pseudoVoigt in deg or musec
793	'''
794
795	sig = lambda Th,U,V,W: np.sqrt(max(0.001,Utand(Th)2+Vtand(Th)+W))
796	sigTOF = lambda dsp,S0,S1,S2,Sq: np.sqrt(S0+S1dsp2+S2dsp*4+Sqdsp)
797	gam = lambda Th,X,Y,Z: Z+X/cosd(Th)+Y*tand(Th)
798	gamTOF = lambda dsp,X,Y,Z: Z+Xdsp+Ydsp**2
799	alpTOF = lambda dsp,alp: alp/dsp
800	betTOF = lambda dsp,bet0,bet1,betq: bet0+bet1/dsp4+betq/dsp2
801	if 'C' in Inst['Type'][0]:
802	s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
803	g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
804	return getgamFW(g,s)/100. #returns FWHM in deg
805	else:
806	dsp = pos/Inst['difC'][0]
807	alp = alpTOF(dsp,Inst['alpha'][0])
808	bet = betTOF(dsp,Inst['beta-0'][0],Inst['beta-1'][0],Inst['beta-q'][0])
809	s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
810	g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
811	return getgamFW(g,s)+np.log(2.0)(alp+bet)/(alpbet)
812
813	def getgamFW(g,s):
814	'''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
815	lambda fxn needs FWHM for both Gaussian & Lorentzian components
816
817	:param g: float Lorentzian gamma = FWHM(L)
818	:param s: float Gaussian sig
819
820	:returns float: total FWHM of pseudoVoigt
821	'''
822	gamFW = lambda s,g: np.exp(np.log(s*5+2.69269s*4g+2.42843s3g*2+4.47163s*2g*3+0.07842sg4+g*5)/5.)
823	return gamFW(2.35482s,g) #sqrt(8ln2)sig = FWHM(G)
824
825	def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
826	'''Compute the Finger-Cox-Jepcoat modified Voigt function for a
827	CW powder peak by direct convolution. This version is not used.
828	'''
829	DX = xdata[1]-xdata[0]
830	widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
831	x = np.linspace(pos-fmin,pos+fmin,256)
832	dx = x[1]-x[0]
833	Norm = norm.pdf(x,loc=pos,scale=widths[0])
834	Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
835	arg = [pos,shl/57.2958,dx,]
836	FCJ = fcjde.pdf(x,*arg,loc=pos)
837	if len(np.nonzero(FCJ)[0])>5:
838	z = np.column_stack([Norm,Cauchy,FCJ]).T
839	Z = fft.fft(z)
840	Df = fft.ifft(Z.prod(axis=0)).real
841	else:
842	z = np.column_stack([Norm,Cauchy]).T
843	Z = fft.fft(z)
844	Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real
845	Df /= np.sum(Df)
846	Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
847	return intensDf(xdata)DX/dx
848
849	def getBackground(pfx,parmDict,bakType,dataType,xdata,fixedBkg={}):
850	'''Computes the background from vars pulled from gpx file or tree.
851	'''
852	if 'T' in dataType:
853	q = 2.np.piparmDict[pfx+'difC']/xdata
854	elif 'C' in dataType:
855	wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
856	q = npT2q(xdata,wave)
857	yb = np.zeros_like(xdata)
858	nBak = 0
859	cw = np.diff(xdata)
860	cw = np.append(cw,cw[-1])
861	sumBk = [0.,0.,0]
862	while True:
863	key = pfx+'Back;'+str(nBak)
864	if key in parmDict:
865	nBak += 1
866	else:
867	break
868	#empirical functions
869	if bakType in ['chebyschev','cosine','chebyschev-1']:
870	dt = xdata[-1]-xdata[0]
871	for iBak in range(nBak):
872	key = pfx+'Back;'+str(iBak)
873	if bakType == 'chebyschev':
874	ybi = parmDict[key](-1.+2.(xdata-xdata[0])/dt)**iBak
875	elif bakType == 'chebyschev-1':
876	xpos = -1.+2.*(xdata-xdata[0])/dt
877	ybi = parmDict[key]np.cos(iBaknp.arccos(xpos))
878	elif bakType == 'cosine':
879	ybi = parmDict[key]npcosd(180.xdata*iBak/xdata[-1])
880	yb += ybi
881	sumBk[0] = np.sum(yb)
882	elif bakType in ['Q^2 power series','Q^-2 power series']:
883	QT = 1.
884	yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
885	for iBak in range(nBak-1):
886	key = pfx+'Back;'+str(iBak+1)
887	if '-2' in bakType:
888	QT = (iBak+1)q**-2
889	else:
890	QT = q*2/(iBak+1)
891	yb += QT*parmDict[key]
892	sumBk[0] = np.sum(yb)
893	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
894	if nBak == 1:
895	yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
896	elif nBak == 2:
897	dX = xdata[-1]-xdata[0]
898	T2 = (xdata-xdata[0])/dX
899	T1 = 1.0-T2
900	yb = parmDict[pfx+'Back;0']T1+parmDict[pfx+'Back;1']T2
901	else:
902	xnomask = ma.getdata(xdata)
903	xmin,xmax = xnomask[0],xnomask[-1]
904	if bakType == 'lin interpolate':
905	bakPos = np.linspace(xmin,xmax,nBak,True)
906	elif bakType == 'inv interpolate':
907	bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
908	elif bakType == 'log interpolate':
909	bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
910	bakPos[0] = xmin
911	bakPos[-1] = xmax
912	bakVals = np.zeros(nBak)
913	for i in range(nBak):
914	bakVals[i] = parmDict[pfx+'Back;'+str(i)]
915	bakInt = si.interp1d(bakPos,bakVals,'linear')
916	yb = bakInt(ma.getdata(xdata))
917	sumBk[0] = np.sum(yb)
918	#Debye function
919	if pfx+'difC' in parmDict:
920	ff = 1.
921	else:
922	try:
923	wave = parmDict[pfx+'Lam']
924	except KeyError:
925	wave = parmDict[pfx+'Lam1']
926	SQ = (q/(4.np.pi))*2
927	FF = G2elem.GetFormFactorCoeff('Si')[0]
928	ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
929	iD = 0
930	while True:
931	try:
932	dbA = parmDict[pfx+'DebyeA;'+str(iD)]
933	dbR = parmDict[pfx+'DebyeR;'+str(iD)]
934	dbU = parmDict[pfx+'DebyeU;'+str(iD)]
935	ybi = ffdbAnp.sin(qdbR)np.exp(-dbUq2)/(qdbR)
936	yb += ybi
937	sumBk[1] += np.sum(ybi)
938	iD += 1
939	except KeyError:
940	break
941	#peaks
942	iD = 0
943	while True:
944	try:
945	pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
946	pkI = max(parmDict[pfx+'BkPkint;'+str(iD)],0.1)
947	pkS = max(parmDict[pfx+'BkPksig;'+str(iD)],1.)
948	pkG = max(parmDict[pfx+'BkPkgam;'+str(iD)],0.1)
949	if 'C' in dataType:
950	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
951	else: #'T'OF
952	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
953	iBeg = np.searchsorted(xdata,pkP-fmin)
954	iFin = np.searchsorted(xdata,pkP+fmax)
955	lenX = len(xdata)
956	if not iBeg:
957	iFin = np.searchsorted(xdata,pkP+fmax)
958	elif iBeg == lenX:
959	iFin = iBeg
960	else:
961	iFin = np.searchsorted(xdata,pkP+fmax)
962	if 'C' in dataType:
963	ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
964	yb[iBeg:iFin] += ybi
965	else: #'T'OF
966	ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
967	yb[iBeg:iFin] += ybi
968	sumBk[2] += np.sum(ybi)
969	iD += 1
970	except KeyError:
971	break
972	except ValueError:
973	G2fil.G2Print ('** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **')
974	break
975	# fixed background from file
976	if len(fixedBkg) >= 3:
977	mult = fixedBkg.get('_fixedMult',0.0)
978	if len(fixedBkg.get('_fixedValues',[])) != len(yb):
979	G2fil.G2Print('Lengths of backgrounds do not agree: yb={}, fixed={}'.format(
980	len(yb),len(fixedBkg.get('_fixedValues',[]))))
981	elif mult:
982	yb -= mult*fixedBkg.get('_fixedValues',[]) # N.B. mult is negative
983	sumBk[0] = sum(yb)
984	return yb,sumBk
985
986	def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
987	'needs a doc string'
988	if 'T' in dataType:
989	q = 2.np.piparmDict[hfx+'difC']/xdata
990	elif 'C' in dataType:
991	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
992	q = 2.np.pinpsind(xdata/2.)/wave
993	nBak = 0
994	while True:
995	key = hfx+'Back;'+str(nBak)
996	if key in parmDict:
997	nBak += 1
998	else:
999	break
1000	dydb = np.zeros(shape=(nBak,len(xdata)))
1001	dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
1002	dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
1003	cw = np.diff(xdata)
1004	cw = np.append(cw,cw[-1])
1005
1006	if bakType in ['chebyschev','cosine','chebyschev-1']:
1007	dt = xdata[-1]-xdata[0]
1008	for iBak in range(nBak):
1009	if bakType == 'chebyschev':
1010	dydb[iBak] = (-1.+2.(xdata-xdata[0])/dt)*iBak
1011	elif bakType == 'chebyschev-1':
1012	xpos = -1.+2.*(xdata-xdata[0])/dt
1013	dydb[iBak] = np.cos(iBak*np.arccos(xpos))
1014	elif bakType == 'cosine':
1015	dydb[iBak] = npcosd(180.xdataiBak/xdata[-1])
1016	elif bakType in ['Q^2 power series','Q^-2 power series']:
1017	QT = 1.
1018	dydb[0] = np.ones_like(xdata)
1019	for iBak in range(nBak-1):
1020	if '-2' in bakType:
1021	QT = (iBak+1)q**-2
1022	else:
1023	QT = q*2/(iBak+1)
1024	dydb[iBak+1] = QT
1025	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
1026	if nBak == 1:
1027	dydb[0] = np.ones_like(xdata)
1028	elif nBak == 2:
1029	dX = xdata[-1]-xdata[0]
1030	T2 = (xdata-xdata[0])/dX
1031	T1 = 1.0-T2
1032	dydb = [T1,T2]
1033	else:
1034	xnomask = ma.getdata(xdata)
1035	xmin,xmax = xnomask[0],xnomask[-1]
1036	if bakType == 'lin interpolate':
1037	bakPos = np.linspace(xmin,xmax,nBak,True)
1038	elif bakType == 'inv interpolate':
1039	bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
1040	elif bakType == 'log interpolate':
1041	bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
1042	bakPos[0] = xmin
1043	bakPos[-1] = xmax
1044	for i,pos in enumerate(bakPos):
1045	if i == 0:
1046	dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
1047	elif i == len(bakPos)-1:
1048	dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
1049	else:
1050	dydb[i] = np.where(xdata>bakPos[i],
1051	np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
1052	np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
1053	if hfx+'difC' in parmDict:
1054	ff = 1.
1055	else:
1056	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1057	q = npT2q(xdata,wave)
1058	SQ = (q/(4np.pi))*2
1059	FF = G2elem.GetFormFactorCoeff('Si')[0]
1060	ff = np.array(G2elem.ScatFac(FF,SQ)[0])np.pi*2 #needs pi^2~10. for cw data (why?)
1061	iD = 0
1062	while True:
1063	try:
1064	if hfx+'difC' in parmDict:
1065	q = 2np.piparmDict[hfx+'difC']/xdata
1066	dbA = parmDict[hfx+'DebyeA;'+str(iD)]
1067	dbR = parmDict[hfx+'DebyeR;'+str(iD)]
1068	dbU = parmDict[hfx+'DebyeU;'+str(iD)]
1069	sqr = np.sin(qdbR)/(qdbR)
1070	cqr = np.cos(q*dbR)
1071	temp = np.exp(-dbUq*2)
1072	dyddb[3iD] = ffsqr*temp
1073	dyddb[3iD+1] = ffdbAtemp(cqr-sqr)/(dbR)
1074	dyddb[3iD+2] = -ffdbAsqrtempq*2
1075	iD += 1
1076	except KeyError:
1077	break
1078	iD = 0
1079	while True:
1080	try:
1081	pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
1082	pkI = max(parmDict[hfx+'BkPkint;'+str(iD)],0.1)
1083	pkS = max(parmDict[hfx+'BkPksig;'+str(iD)],1.0)
1084	pkG = max(parmDict[hfx+'BkPkgam;'+str(iD)],0.1)
1085	if 'C' in dataType:
1086	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
1087	else: #'T'OF
1088	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
1089	iBeg = np.searchsorted(xdata,pkP-fmin)
1090	iFin = np.searchsorted(xdata,pkP+fmax)
1091	lenX = len(xdata)
1092	if not iBeg:
1093	iFin = np.searchsorted(xdata,pkP+fmax)
1094	elif iBeg == lenX:
1095	iFin = iBeg
1096	else:
1097	iFin = np.searchsorted(xdata,pkP+fmax)
1098	if 'C' in dataType:
1099	Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
1100	dydpk[4iD][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdp
1101	dydpk[4iD+1][iBeg:iFin] += 100.cw[iBeg:iFin]*Df
1102	dydpk[4iD+2][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFds
1103	dydpk[4iD+3][iBeg:iFin] += 100.cw[iBeg:iFin]pkIdFdg
1104	else: #'T'OF
1105	Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
1106	dydpk[4iD][iBeg:iFin] += pkIdFdp
1107	dydpk[4*iD+1][iBeg:iFin] += Df
1108	dydpk[4iD+2][iBeg:iFin] += pkIdFds
1109	dydpk[4iD+3][iBeg:iFin] += pkIdFdg
1110	iD += 1
1111	except KeyError:
1112	break
1113	except ValueError:
1114	G2fil.G2Print ('** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **')
1115	break
1116	return dydb,dyddb,dydpk
1117
1118	#use old fortran routine
1119	def getFCJVoigt3(pos,sig,gam,shl,xdata):
1120	'''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
1121	CW powder peak in external Fortran routine
1122	'''
1123	Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1124	# Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
1125	Df /= np.sum(Df)
1126	return Df
1127
1128	def getdFCJVoigt3(pos,sig,gam,shl,xdata):
1129	'''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
1130	function for a CW powder peak
1131	'''
1132	Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1133	# Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
1134	return Df,dFdp,dFds,dFdg,dFdsh
1135
1136	def getPsVoigt(pos,sig,gam,xdata):
1137	'needs a doc string'
1138
1139	Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
1140	Df /= np.sum(Df)
1141	return Df
1142
1143	def getdPsVoigt(pos,sig,gam,xdata):
1144	'needs a doc string'
1145
1146	Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
1147	return Df,dFdp,dFds,dFdg
1148
1149	def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
1150	'needs a doc string'
1151	Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1152	Df /= np.sum(Df)
1153	return Df
1154
1155	def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
1156	'needs a doc string'
1157	Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1158	return Df,dFdp,dFda,dFdb,dFds,dFdg
1159
1160	def ellipseSize(H,Sij,GB):
1161	'Implements r=1/sqrt(sum((1/S)*(q.v)^2) per note from Alexander Brady'
1162	HX = np.inner(H.T,GB)
1163	lenHX = np.sqrt(np.sum(HX**2))
1164	Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))
1165	R = np.inner(HX/lenHX,Rsize)*2Esize #want column length for hkl in crystal
1166	lenR = 1./np.sqrt(np.sum(R))
1167	return lenR
1168
1169	def ellipseSizeDerv(H,Sij,GB):
1170	'needs a doc string'
1171	lenR = ellipseSize(H,Sij,GB)
1172	delt = 0.001
1173	dRdS = np.zeros(6)
1174	for i in range(6):
1175	Sij[i] -= delt
1176	lenM = ellipseSize(H,Sij,GB)
1177	Sij[i] += 2.*delt
1178	lenP = ellipseSize(H,Sij,GB)
1179	Sij[i] -= delt
1180	dRdS[i] = (lenP-lenM)/(2.*delt)
1181	return lenR,dRdS
1182
1183	def getMustrain(HKL,G,SGData,muStrData):
1184	if muStrData[0] == 'isotropic':
1185	return np.ones(HKL.shape[1])*muStrData[1][0]
1186	elif muStrData[0] == 'uniaxial':
1187	H = np.array(HKL)
1188	P = np.array(muStrData[3])
1189	cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1190	Si = muStrData[1][0]
1191	Sa = muStrData[1][1]
1192	return SiSa/(np.sqrt((SicosP)*2+(SasinP)**2))
1193	else: #generalized - P.W. Stephens model
1194	H = np.array(HKL)
1195	rdsq = np.array([G2lat.calc_rDsq2(h,G) for h in H.T])
1196	Strms = np.array(G2spc.MustrainCoeff(H,SGData))
1197	Sum = np.sum(np.array(muStrData[4])[:,nxs]*Strms,axis=0)
1198	return np.sqrt(Sum)/rdsq
1199
1200	def getCrSize(HKL,G,GB,sizeData):
1201	if sizeData[0] == 'isotropic':
1202	return np.ones(HKL.shape[1])*sizeData[1][0]
1203	elif sizeData[0] == 'uniaxial':
1204	H = np.array(HKL)
1205	P = np.array(sizeData[3])
1206	cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1207	Si = sizeData[1][0]
1208	Sa = sizeData[1][1]
1209	return SiSa/(np.sqrt((SicosP)*2+(SasinP)**2))
1210	else:
1211	Sij =[sizeData[4][i] for i in range(6)]
1212	H = np.array(HKL)
1213	return 1./np.array([ellipseSize(h,Sij,GB) for h in H.T])**2
1214
1215	def getHKLpeak(dmin,SGData,A,Inst=None,nodup=False):
1216	'''
1217	Generates allowed by symmetry reflections with d >= dmin
1218	NB: GenHKLf & checkMagextc return True for extinct reflections
1219
1220	:param dmin: minimum d-spacing
1221	:param SGData: space group data obtained from SpcGroup
1222	:param A: lattice parameter terms A1-A6
1223	:param Inst: instrument parameter info
1224	:returns: HKLs: np.array hkl, etc for allowed reflections
1225
1226	'''
1227	HKL = G2lat.GenHLaue(dmin,SGData,A)
1228	HKLs = []
1229	ds = []
1230	for h,k,l,d in HKL:
1231	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1232	if ext and 'MagSpGrp' in SGData:
1233	ext = G2spc.checkMagextc([h,k,l],SGData)
1234	if not ext:
1235	if nodup and int(10000*d) in ds:
1236	continue
1237	ds.append(int(10000*d))
1238	if Inst == None:
1239	HKLs.append([h,k,l,d,0,-1])
1240	else:
1241	HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1242	return np.array(HKLs)
1243
1244	def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1245	'needs a doc string'
1246	HKLs = []
1247	vec = np.array(Vec)
1248	vstar = np.sqrt(G2lat.calc_rDsq(vec,A)) #find extra needed for -n SS reflections
1249	dvec = 1./(maxH*vstar+1./dmin)
1250	HKL = G2lat.GenHLaue(dvec,SGData,A)
1251	SSdH = [vec*h for h in range(-maxH,maxH+1)]
1252	SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1253	ifMag = False
1254	if 'MagSpGrp' in SGData:
1255	ifMag = True
1256	for h,k,l,d in HKL:
1257	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1258	if not ext and d >= dmin:
1259	HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1260	for dH in SSdH:
1261	if dH:
1262	DH = SSdH[dH]
1263	H = [h+DH[0],k+DH[1],l+DH[2]]
1264	d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1265	if d >= dmin:
1266	HKLM = np.array([h,k,l,dH])
1267	if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1268	HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])
1269	return G2lat.sortHKLd(HKLs,True,True,True)
1270
1271	def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
1272	'Computes the profile for a powder pattern'
1273
1274	yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
1275	yc = np.zeros_like(yb)
1276	cw = np.diff(xdata)
1277	cw = np.append(cw,cw[-1])
1278	if 'C' in dataType:
1279	shl = max(parmDict['SH/L'],0.002)
1280	Ka2 = False
1281	if 'Lam1' in parmDict.keys():
1282	Ka2 = True
1283	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1284	kRatio = parmDict['I(L2)/I(L1)']
1285	iPeak = 0
1286	while True:
1287	try:
1288	pos = parmDict['pos'+str(iPeak)]
1289	tth = (pos-parmDict['Zero'])
1290	intens = parmDict['int'+str(iPeak)]
1291	sigName = 'sig'+str(iPeak)
1292	if sigName in varyList:
1293	sig = parmDict[sigName]
1294	else:
1295	sig = G2mth.getCWsig(parmDict,tth)
1296	sig = max(sig,0.001) #avoid neg sigma^2
1297	gamName = 'gam'+str(iPeak)
1298	if gamName in varyList:
1299	gam = parmDict[gamName]
1300	else:
1301	gam = G2mth.getCWgam(parmDict,tth)
1302	gam = max(gam,0.001) #avoid neg gamma
1303	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1304	iBeg = np.searchsorted(xdata,pos-fmin)
1305	iFin = np.searchsorted(xdata,pos+fmin)
1306	if not iBeg+iFin: #peak below low limit
1307	iPeak += 1
1308	continue
1309	elif not iBeg-iFin: #peak above high limit
1310	return yb+yc
1311	yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1312	if Ka2:
1313	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1314	iBeg = np.searchsorted(xdata,pos2-fmin)
1315	iFin = np.searchsorted(xdata,pos2+fmin)
1316	if iBeg-iFin:
1317	yc[iBeg:iFin] += intenskRatiogetFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1318	iPeak += 1
1319	except KeyError: #no more peaks to process
1320	return yb+yc
1321	else:
1322	Pdabc = parmDict['Pdabc']
1323	difC = parmDict['difC']
1324	iPeak = 0
1325	while True:
1326	try:
1327	pos = parmDict['pos'+str(iPeak)]
1328	tof = pos-parmDict['Zero']
1329	dsp = tof/difC
1330	intens = parmDict['int'+str(iPeak)]
1331	alpName = 'alp'+str(iPeak)
1332	if alpName in varyList:
1333	alp = parmDict[alpName]
1334	else:
1335	if len(Pdabc):
1336	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1337	else:
1338	alp = G2mth.getTOFalpha(parmDict,dsp)
1339	alp = max(0.1,alp)
1340	betName = 'bet'+str(iPeak)
1341	if betName in varyList:
1342	bet = parmDict[betName]
1343	else:
1344	if len(Pdabc):
1345	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1346	else:
1347	bet = G2mth.getTOFbeta(parmDict,dsp)
1348	bet = max(0.0001,bet)
1349	sigName = 'sig'+str(iPeak)
1350	if sigName in varyList:
1351	sig = parmDict[sigName]
1352	else:
1353	sig = G2mth.getTOFsig(parmDict,dsp)
1354	gamName = 'gam'+str(iPeak)
1355	if gamName in varyList:
1356	gam = parmDict[gamName]
1357	else:
1358	gam = G2mth.getTOFgamma(parmDict,dsp)
1359	gam = max(gam,0.001) #avoid neg gamma
1360	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1361	iBeg = np.searchsorted(xdata,pos-fmin)
1362	iFin = np.searchsorted(xdata,pos+fmax)
1363	lenX = len(xdata)
1364	if not iBeg:
1365	iFin = np.searchsorted(xdata,pos+fmax)
1366	elif iBeg == lenX:
1367	iFin = iBeg
1368	else:
1369	iFin = np.searchsorted(xdata,pos+fmax)
1370	if not iBeg+iFin: #peak below low limit
1371	iPeak += 1
1372	continue
1373	elif not iBeg-iFin: #peak above high limit
1374	return yb+yc
1375	yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1376	iPeak += 1
1377	except KeyError: #no more peaks to process
1378	return yb+yc
1379
1380	def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
1381	'needs a doc string'
1382	# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1383	dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1384	dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
1385	if 'Back;0' in varyList: #background derivs are in front if present
1386	dMdv[0:len(dMdb)] = dMdb
1387	names = ['DebyeA','DebyeR','DebyeU']
1388	for name in varyList:
1389	if 'Debye' in name:
1390	parm,Id = name.split(';')
1391	ip = names.index(parm)
1392	dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip]
1393	names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1394	for name in varyList:
1395	if 'BkPk' in name:
1396	parm,Id = name.split(';')
1397	ip = names.index(parm)
1398	dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
1399	cw = np.diff(xdata)
1400	cw = np.append(cw,cw[-1])
1401	if 'C' in dataType:
1402	shl = max(parmDict['SH/L'],0.002)
1403	Ka2 = False
1404	if 'Lam1' in parmDict.keys():
1405	Ka2 = True
1406	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1407	kRatio = parmDict['I(L2)/I(L1)']
1408	iPeak = 0
1409	while True:
1410	try:
1411	pos = parmDict['pos'+str(iPeak)]
1412	tth = (pos-parmDict['Zero'])
1413	intens = parmDict['int'+str(iPeak)]
1414	sigName = 'sig'+str(iPeak)
1415	if sigName in varyList:
1416	sig = parmDict[sigName]
1417	dsdU = dsdV = dsdW = 0
1418	else:
1419	sig = G2mth.getCWsig(parmDict,tth)
1420	dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1421	sig = max(sig,0.001) #avoid neg sigma
1422	gamName = 'gam'+str(iPeak)
1423	if gamName in varyList:
1424	gam = parmDict[gamName]
1425	dgdX = dgdY = dgdZ = 0
1426	else:
1427	gam = G2mth.getCWgam(parmDict,tth)
1428	dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1429	gam = max(gam,0.001) #avoid neg gamma
1430	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1431	iBeg = np.searchsorted(xdata,pos-fmin)
1432	iFin = np.searchsorted(xdata,pos+fmin)
1433	if not iBeg+iFin: #peak below low limit
1434	iPeak += 1
1435	continue
1436	elif not iBeg-iFin: #peak above high limit
1437	break
1438	dMdpk = np.zeros(shape=(6,len(xdata)))
1439	dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1440	for i in range(1,5):
1441	dMdpk[i][iBeg:iFin] += 100.cw[iBeg:iFin]intens*dMdipk[i]
1442	dMdpk[0][iBeg:iFin] += 100.cw[iBeg:iFin]dMdipk[0]
1443	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1444	if Ka2:
1445	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1446	iBeg = np.searchsorted(xdata,pos2-fmin)
1447	iFin = np.searchsorted(xdata,pos2+fmin)
1448	if iBeg-iFin:
1449	dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1450	for i in range(1,5):
1451	dMdpk[i][iBeg:iFin] += 100.cw[iBeg:iFin]intenskRatiodMdipk2[i]
1452	dMdpk[0][iBeg:iFin] += 100.cw[iBeg:iFin]kRatio*dMdipk2[0]
1453	dMdpk[5][iBeg:iFin] += 100.cw[iBeg:iFin]dMdipk2[0]
1454	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1455	for parmName in ['pos','int','sig','gam']:
1456	try:
1457	idx = varyList.index(parmName+str(iPeak))
1458	dMdv[idx] = dervDict[parmName]
1459	except ValueError:
1460	pass
1461	if 'U' in varyList:
1462	dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1463	if 'V' in varyList:
1464	dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1465	if 'W' in varyList:
1466	dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1467	if 'X' in varyList:
1468	dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1469	if 'Y' in varyList:
1470	dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1471	if 'Z' in varyList:
1472	dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1473	if 'SH/L' in varyList:
1474	dMdv[varyList.index('SH/L')] += dervDict['shl'] #problem here
1475	if 'I(L2)/I(L1)' in varyList:
1476	dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1477	iPeak += 1
1478	except KeyError: #no more peaks to process
1479	break
1480	else:
1481	Pdabc = parmDict['Pdabc']
1482	difC = parmDict['difC']
1483	iPeak = 0
1484	while True:
1485	try:
1486	pos = parmDict['pos'+str(iPeak)]
1487	tof = pos-parmDict['Zero']
1488	dsp = tof/difC
1489	intens = parmDict['int'+str(iPeak)]
1490	alpName = 'alp'+str(iPeak)
1491	if alpName in varyList:
1492	alp = parmDict[alpName]
1493	else:
1494	if len(Pdabc):
1495	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1496	dada0 = 0
1497	else:
1498	alp = G2mth.getTOFalpha(parmDict,dsp)
1499	dada0 = G2mth.getTOFalphaDeriv(dsp)
1500	betName = 'bet'+str(iPeak)
1501	if betName in varyList:
1502	bet = parmDict[betName]
1503	else:
1504	if len(Pdabc):
1505	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1506	dbdb0 = dbdb1 = dbdb2 = 0
1507	else:
1508	bet = G2mth.getTOFbeta(parmDict,dsp)
1509	dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1510	sigName = 'sig'+str(iPeak)
1511	if sigName in varyList:
1512	sig = parmDict[sigName]
1513	dsds0 = dsds1 = dsds2 = dsds3 = 0
1514	else:
1515	sig = G2mth.getTOFsig(parmDict,dsp)
1516	dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1517	gamName = 'gam'+str(iPeak)
1518	if gamName in varyList:
1519	gam = parmDict[gamName]
1520	dsdX = dsdY = dsdZ = 0
1521	else:
1522	gam = G2mth.getTOFgamma(parmDict,dsp)
1523	dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1524	gam = max(gam,0.001) #avoid neg gamma
1525	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1526	iBeg = np.searchsorted(xdata,pos-fmin)
1527	lenX = len(xdata)
1528	if not iBeg:
1529	iFin = np.searchsorted(xdata,pos+fmax)
1530	elif iBeg == lenX:
1531	iFin = iBeg
1532	else:
1533	iFin = np.searchsorted(xdata,pos+fmax)
1534	if not iBeg+iFin: #peak below low limit
1535	iPeak += 1
1536	continue
1537	elif not iBeg-iFin: #peak above high limit
1538	break
1539	dMdpk = np.zeros(shape=(7,len(xdata)))
1540	dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1541	for i in range(1,6):
1542	dMdpk[i][iBeg:iFin] += intenscw[iBeg:iFin]dMdipk[i]
1543	dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1544	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1545	for parmName in ['pos','int','alp','bet','sig','gam']:
1546	try:
1547	idx = varyList.index(parmName+str(iPeak))
1548	dMdv[idx] = dervDict[parmName]
1549	except ValueError:
1550	pass
1551	if 'alpha' in varyList:
1552	dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1553	if 'beta-0' in varyList:
1554	dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1555	if 'beta-1' in varyList:
1556	dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1557	if 'beta-q' in varyList:
1558	dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1559	if 'sig-0' in varyList:
1560	dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1561	if 'sig-1' in varyList:
1562	dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1563	if 'sig-2' in varyList:
1564	dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1565	if 'sig-q' in varyList:
1566	dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1567	if 'X' in varyList:
1568	dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1569	if 'Y' in varyList:
1570	dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1571	if 'Z' in varyList:
1572	dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1573	iPeak += 1
1574	except KeyError: #no more peaks to process
1575	break
1576	return dMdv
1577
1578	def Dict2Values(parmdict, varylist):
1579	'''Use before call to leastsq to setup list of values for the parameters
1580	in parmdict, as selected by key in varylist'''
1581	return [parmdict[key] for key in varylist]
1582
1583	def Values2Dict(parmdict, varylist, values):
1584	''' Use after call to leastsq to update the parameter dictionary with
1585	values corresponding to keys in varylist'''
1586	parmdict.update(zip(varylist,values))
1587
1588	def SetBackgroundParms(Background):
1589	'Loads background parameters into dicts/lists to create varylist & parmdict'
1590	if len(Background) == 1: # fix up old backgrounds
1591	Background.append({'nDebye':0,'debyeTerms':[]})
1592	bakType,bakFlag = Background[0][:2]
1593	backVals = Background[0][3:]
1594	backNames = ['Back;'+str(i) for i in range(len(backVals))]
1595	Debye = Background[1] #also has background peaks stuff
1596	backDict = dict(zip(backNames,backVals))
1597	backVary = []
1598	if bakFlag:
1599	backVary = backNames
1600
1601	backDict['nDebye'] = Debye['nDebye']
1602	debyeDict = {}
1603	debyeList = []
1604	for i in range(Debye['nDebye']):
1605	debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1606	debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1607	debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1608	debyeVary = []
1609	for item in debyeList:
1610	if item[1]:
1611	debyeVary.append(item[0])
1612	backDict.update(debyeDict)
1613	backVary += debyeVary
1614
1615	backDict['nPeaks'] = Debye['nPeaks']
1616	peaksDict = {}
1617	peaksList = []
1618	for i in range(Debye['nPeaks']):
1619	peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1620	peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1621	peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1622	peaksVary = []
1623	for item in peaksList:
1624	if item[1]:
1625	peaksVary.append(item[0])
1626	backDict.update(peaksDict)
1627	backVary += peaksVary
1628	return bakType,backDict,backVary
1629
1630	def DoCalibInst(IndexPeaks,Inst):
1631
1632	def SetInstParms():
1633	dataType = Inst['Type'][0]
1634	insVary = []
1635	insNames = []
1636	insVals = []
1637	for parm in Inst:
1638	insNames.append(parm)
1639	insVals.append(Inst[parm][1])
1640	if parm in ['Lam','difC','difA','difB','Zero',]:
1641	if Inst[parm][2]:
1642	insVary.append(parm)
1643	instDict = dict(zip(insNames,insVals))
1644	return dataType,instDict,insVary
1645
1646	def GetInstParms(parmDict,Inst,varyList):
1647	for name in Inst:
1648	Inst[name][1] = parmDict[name]
1649
1650	def InstPrint(Inst,sigDict):
1651	print ('Instrument Parameters:')
1652	if 'C' in Inst['Type'][0]:
1653	ptfmt = "%12.6f"
1654	else:
1655	ptfmt = "%12.3f"
1656	ptlbls = 'names :'
1657	ptstr = 'values:'
1658	sigstr = 'esds :'
1659	for parm in Inst:
1660	if parm in ['Lam','difC','difA','difB','Zero',]:
1661	ptlbls += "%s" % (parm.center(12))
1662	ptstr += ptfmt % (Inst[parm][1])
1663	if parm in sigDict:
1664	sigstr += ptfmt % (sigDict[parm])
1665	else:
1666	sigstr += 12*' '
1667	print (ptlbls)
1668	print (ptstr)
1669	print (sigstr)
1670
1671	def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1672	parmDict.update(zip(varyList,values))
1673	return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1674
1675	peakPos = []
1676	peakDsp = []
1677	peakWt = []
1678	for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1679	if peak[2] and peak[3] and sig > 0.:
1680	peakPos.append(peak[0])
1681	peakDsp.append(peak[-1]) #d-calc
1682	# peakWt.append(peak[-1]2/sig2) #weight by d**2
1683	peakWt.append(1./(sig*peak[-1])) #
1684	peakPos = np.array(peakPos)
1685	peakDsp = np.array(peakDsp)
1686	peakWt = np.array(peakWt)
1687	dataType,insDict,insVary = SetInstParms()
1688	parmDict = {}
1689	parmDict.update(insDict)
1690	varyList = insVary
1691	if not len(varyList):
1692	G2fil.G2Print ('** ERROR - nothing to refine! **')
1693	return False
1694	while True:
1695	begin = time.time()
1696	values = np.array(Dict2Values(parmDict, varyList))
1697	result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1698	args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1699	ncyc = int(result[2]['nfev']/2)
1700	runtime = time.time()-begin
1701	chisq = np.sum(result[2]['fvec']**2)
1702	Values2Dict(parmDict, varyList, result[0])
1703	GOF = chisq/(len(peakPos)-len(varyList)) #reduced chi^2
1704	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
1705	G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
1706	G2fil.G2Print ('chi2 = %12.6g, reduced chi2 = %6.2f'%(chisq,GOF))
1707	try:
1708	sig = np.sqrt(np.diag(result[1])*GOF)
1709	if np.any(np.isnan(sig)):
1710	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *')
1711	break #refinement succeeded - finish up!
1712	except ValueError: #result[1] is None on singular matrix
1713	G2fil.G2Print ('** Refinement failed - singular matrix **')
1714
1715	sigDict = dict(zip(varyList,sig))
1716	GetInstParms(parmDict,Inst,varyList)
1717	InstPrint(Inst,sigDict)
1718	return True
1719
1720	def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1721	'''Called to perform a peak fit, refining the selected items in the peak
1722	table as well as selected items in the background.
1723
1724	:param str FitPgm: type of fit to perform. At present this is ignored.
1725	:param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
1726	four values followed by a refine flag where the values are: position, intensity,
1727	sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
1728	tree entry, dict item "peaks"
1729	:param list Background: describes the background. List with two items.
1730	Item 0 specifies a background model and coefficients. Item 1 is a dict.
1731	From the Histogram/Background tree entry.
1732	:param list Limits: min and max x-value to use
1733	:param dict Inst: Instrument parameters
1734	:param dict Inst2: more Instrument parameters
1735	:param numpy.array data: a 5xn array. data[0] is the x-values,
1736	data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
1737	calc, background and difference intensities, respectively.
1738	:param array fixback: fixed background values
1739	:param list prevVaryList: Used in sequential refinements to override the
1740	variable list. Defaults as an empty list.
1741	:param bool oneCycle: True if only one cycle of fitting should be performed
1742	:param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
1743	and derivType = controls['deriv type']. If None default values are used.
1744	:param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
1745	Defaults to None, which means no updates are done.
1746	'''
1747	def GetBackgroundParms(parmList,Background):
1748	iBak = 0
1749	while True:
1750	try:
1751	bakName = 'Back;'+str(iBak)
1752	Background[0][iBak+3] = parmList[bakName]
1753	iBak += 1
1754	except KeyError:
1755	break
1756	iDb = 0
1757	while True:
1758	names = ['DebyeA;','DebyeR;','DebyeU;']
1759	try:
1760	for i,name in enumerate(names):
1761	val = parmList[name+str(iDb)]
1762	Background[1]['debyeTerms'][iDb][2*i] = val
1763	iDb += 1
1764	except KeyError:
1765	break
1766	iDb = 0
1767	while True:
1768	names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1769	try:
1770	for i,name in enumerate(names):
1771	val = parmList[name+str(iDb)]
1772	Background[1]['peaksList'][iDb][2*i] = val
1773	iDb += 1
1774	except KeyError:
1775	break
1776
1777	def BackgroundPrint(Background,sigDict):
1778	print ('Background coefficients for '+Background[0][0]+' function')
1779	ptfmt = "%12.5f"
1780	ptstr = 'value: '
1781	sigstr = 'esd : '
1782	for i,back in enumerate(Background[0][3:]):
1783	ptstr += ptfmt % (back)
1784	if Background[0][1]:
1785	prm = 'Back;'+str(i)
1786	if prm in sigDict:
1787	sigstr += ptfmt % (sigDict[prm])
1788	else:
1789	sigstr += " "*12
1790	if len(ptstr) > 75:
1791	print (ptstr)
1792	if Background[0][1]: print (sigstr)
1793	ptstr = 'value: '
1794	sigstr = 'esd : '
1795	if len(ptstr) > 8:
1796	print (ptstr)
1797	if Background[0][1]: print (sigstr)
1798
1799	if Background[1]['nDebye']:
1800	parms = ['DebyeA;','DebyeR;','DebyeU;']
1801	print ('Debye diffuse scattering coefficients')
1802	ptfmt = "%12.5f"
1803	print (' term DebyeA esd DebyeR esd DebyeU esd')
1804	for term in range(Background[1]['nDebye']):
1805	line = ' term %d'%(term)
1806	for ip,name in enumerate(parms):
1807	line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1808	if name+str(term) in sigDict:
1809	line += ptfmt%(sigDict[name+str(term)])
1810	else:
1811	line += " "*12
1812	print (line)
1813	if Background[1]['nPeaks']:
1814	print ('Coefficients for Background Peaks')
1815	ptfmt = "%15.3f"
1816	for j,pl in enumerate(Background[1]['peaksList']):
1817	names = 'peak %3d:'%(j+1)
1818	ptstr = 'values :'
1819	sigstr = 'esds :'
1820	for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1821	val = pl[2*i]
1822	prm = lbl+";"+str(j)
1823	names += '%15s'%(prm)
1824	ptstr += ptfmt%(val)
1825	if prm in sigDict:
1826	sigstr += ptfmt%(sigDict[prm])
1827	else:
1828	sigstr += " "*15
1829	print (names)
1830	print (ptstr)
1831	print (sigstr)
1832
1833	def SetInstParms(Inst):
1834	dataType = Inst['Type'][0]
1835	insVary = []
1836	insNames = []
1837	insVals = []
1838	for parm in Inst:
1839	insNames.append(parm)
1840	insVals.append(Inst[parm][1])
1841	if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1842	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1843	insVary.append(parm)
1844	instDict = dict(zip(insNames,insVals))
1845	# instDict['X'] = max(instDict['X'],0.01)
1846	# instDict['Y'] = max(instDict['Y'],0.01)
1847	if 'SH/L' in instDict:
1848	instDict['SH/L'] = max(instDict['SH/L'],0.002)
1849	return dataType,instDict,insVary
1850
1851	def GetInstParms(parmDict,Inst,varyList,Peaks):
1852	for name in Inst:
1853	Inst[name][1] = parmDict[name]
1854	iPeak = 0
1855	while True:
1856	try:
1857	sigName = 'sig'+str(iPeak)
1858	pos = parmDict['pos'+str(iPeak)]
1859	if sigName not in varyList:
1860	if 'C' in Inst['Type'][0]:
1861	parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1862	else:
1863	dsp = G2lat.Pos2dsp(Inst,pos)
1864	parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1865	gamName = 'gam'+str(iPeak)
1866	if gamName not in varyList:
1867	if 'C' in Inst['Type'][0]:
1868	parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1869	else:
1870	dsp = G2lat.Pos2dsp(Inst,pos)
1871	parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1872	iPeak += 1
1873	except KeyError:
1874	break
1875
1876	def InstPrint(Inst,sigDict):
1877	print ('Instrument Parameters:')
1878	ptfmt = "%12.6f"
1879	ptlbls = 'names :'
1880	ptstr = 'values:'
1881	sigstr = 'esds :'
1882	for parm in Inst:
1883	if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1884	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1885	ptlbls += "%s" % (parm.center(12))
1886	ptstr += ptfmt % (Inst[parm][1])
1887	if parm in sigDict:
1888	sigstr += ptfmt % (sigDict[parm])
1889	else:
1890	sigstr += 12*' '
1891	print (ptlbls)
1892	print (ptstr)
1893	print (sigstr)
1894
1895	def SetPeaksParms(dataType,Peaks):
1896	peakNames = []
1897	peakVary = []
1898	peakVals = []
1899	if 'C' in dataType:
1900	names = ['pos','int','sig','gam']
1901	else:
1902	names = ['pos','int','alp','bet','sig','gam']
1903	for i,peak in enumerate(Peaks):
1904	for j,name in enumerate(names):
1905	peakVals.append(peak[2*j])
1906	parName = name+str(i)
1907	peakNames.append(parName)
1908	if peak[2*j+1]:
1909	peakVary.append(parName)
1910	return dict(zip(peakNames,peakVals)),peakVary
1911
1912	def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1913	if 'C' in Inst['Type'][0]:
1914	names = ['pos','int','sig','gam']
1915	else: #'T'
1916	names = ['pos','int','alp','bet','sig','gam']
1917	for i,peak in enumerate(Peaks):
1918	pos = parmDict['pos'+str(i)]
1919	if 'difC' in Inst:
1920	dsp = pos/Inst['difC'][1]
1921	for j in range(len(names)):
1922	parName = names[j]+str(i)
1923	if parName in varyList:
1924	peak[2*j] = parmDict[parName]
1925	elif 'alpha' in parName:
1926	peak[2*j] = parmDict['alpha']/dsp
1927	elif 'beta' in parName:
1928	peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1929	elif 'sig' in parName:
1930	if 'C' in Inst['Type'][0]:
1931	peak[2*j] = G2mth.getCWsig(parmDict,pos)
1932	else:
1933	peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1934	elif 'gam' in parName:
1935	if 'C' in Inst['Type'][0]:
1936	peak[2*j] = G2mth.getCWgam(parmDict,pos)
1937	else:
1938	peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1939
1940	def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
1941	print ('Peak coefficients:')
1942	if 'C' in dataType:
1943	names = ['pos','int','sig','gam']
1944	else: #'T'
1945	names = ['pos','int','alp','bet','sig','gam']
1946	head = 13*' '
1947	for name in names:
1948	if name in ['alp','bet']:
1949	head += name.center(8)+'esd'.center(8)
1950	else:
1951	head += name.center(10)+'esd'.center(10)
1952	head += 'bins'.center(8)
1953	print (head)
1954	if 'C' in dataType:
1955	ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1956	else:
1957	ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1958	for i,peak in enumerate(Peaks):
1959	ptstr = ':'
1960	for j in range(len(names)):
1961	name = names[j]
1962	parName = name+str(i)
1963	ptstr += ptfmt[name] % (parmDict[parName])
1964	if parName in varyList:
1965	ptstr += ptfmt[name] % (sigDict[parName])
1966	else:
1967	if name in ['alp','bet']:
1968	ptstr += 8*' '
1969	else:
1970	ptstr += 10*' '
1971	ptstr += '%9.2f'%(ptsperFW[i])
1972	print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
1973
1974	def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1975	parmdict.update(zip(varylist,values))
1976	return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1977
1978	def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):
1979	parmdict.update(zip(varylist,values))
1980	M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1981	Rwp = min(100.,np.sqrt(np.sum(M*2)/np.sum(weightsydata*2))100.)
1982	if dlg:
1983	dlg.Raise()
1984	GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1985	if not GoOn:
1986	return -M #abort!!
1987	return M
1988
1989	if controls:
1990	Ftol = controls['min dM/M']
1991	else:
1992	Ftol = 0.0001
1993	if oneCycle:
1994	Ftol = 1.0
1995	x,y,w,yc,yb,yd = data #these are numpy arrays - remove masks!
1996	if fixback is None:
1997	fixback = np.zeros_like(y)
1998	yc *= 0. #set calcd ones to zero
1999	yb *= 0.
2000	yd *= 0.
2001	cw = x[1:]-x[:-1]
2002	xBeg = np.searchsorted(x,Limits[0])
2003	xFin = np.searchsorted(x,Limits[1])+1
2004	bakType,bakDict,bakVary = SetBackgroundParms(Background)
2005	dataType,insDict,insVary = SetInstParms(Inst)
2006	peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
2007	parmDict = {}
2008	parmDict.update(bakDict)
2009	parmDict.update(insDict)
2010	parmDict.update(peakDict)
2011	parmDict['Pdabc'] = [] #dummy Pdabc
2012	parmDict.update(Inst2) #put in real one if there
2013	if prevVaryList:
2014	varyList = prevVaryList[:]
2015	else:
2016	varyList = bakVary+insVary+peakVary
2017	fullvaryList = varyList[:]
2018	while True:
2019	begin = time.time()
2020	values = np.array(Dict2Values(parmDict, varyList))
2021	Rvals = {}
2022	badVary = []
2023	result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
2024	args=(x[xBeg:xFin],(y+fixback)[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
2025	ncyc = int(result[2]['nfev']/2)
2026	runtime = time.time()-begin
2027	chisq = np.sum(result[2]['fvec']**2)
2028	Values2Dict(parmDict, varyList, result[0])
2029	Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin](y+fixback)[xBeg:xFin]2))100. #to %
2030	Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList)) #reduced chi^2
2031	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
2032	if ncyc:
2033	G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2034	G2fil.G2Print ('Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2035	sig = [0]*len(varyList)
2036	if len(varyList) == 0: break # if nothing was refined
2037	try:
2038	sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
2039	if np.any(np.isnan(sig)):
2040	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *')
2041	break #refinement succeeded - finish up!
2042	except ValueError: #result[1] is None on singular matrix
2043	G2fil.G2Print ('** Refinement failed - singular matrix **')
2044	Ipvt = result[2]['ipvt']
2045	for i,ipvt in enumerate(Ipvt):
2046	if not np.sum(result[2]['fjac'],axis=1)[i]:
2047	G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
2048	badVary.append(varyList[ipvt-1])
2049	del(varyList[ipvt-1])
2050	break
2051	else: # nothing removed
2052	break
2053	if dlg: dlg.Destroy()
2054	sigDict = dict(zip(varyList,sig))
2055	yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]-fixback[xBeg:xFin]
2056	yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)-fixback[xBeg:xFin]
2057	yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
2058	GetBackgroundParms(parmDict,Background)
2059	if bakVary: BackgroundPrint(Background,sigDict)
2060	GetInstParms(parmDict,Inst,varyList,Peaks)
2061	if insVary: InstPrint(Inst,sigDict)
2062	GetPeaksParms(Inst,parmDict,Peaks,varyList)
2063	binsperFWHM = []
2064	for peak in Peaks:
2065	FWHM = getFWHM(peak[0],Inst)
2066	try:
2067	binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
2068	except IndexError:
2069	binsperFWHM.append(0.)
2070	if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
2071	if len(binsperFWHM):
2072	if min(binsperFWHM) < 1.:
2073	G2fil.G2Print ('* Warning: calculated peak widths are too narrow to refine profile coefficients *')
2074	if 'T' in Inst['Type'][0]:
2075	G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
2076	else:
2077	G2fil.G2Print (' Manually increase W in Instrument Parameters')
2078	elif min(binsperFWHM) < 4.:
2079	G2fil.G2Print ('* Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement *')
2080	G2fil.G2Print ('* recommended is 6-10; you have %.2f *'%(min(binsperFWHM)))
2081	return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
2082
2083	def calcIncident(Iparm,xdata):
2084	'needs a doc string'
2085
2086	def IfunAdv(Iparm,xdata):
2087	Itype = Iparm['Itype']
2088	Icoef = Iparm['Icoeff']
2089	DYI = np.ones((12,xdata.shape[0]))
2090	YI = np.ones_like(xdata)*Icoef[0]
2091
2092	x = xdata/1000. #expressions are in ms
2093	if Itype == 'Exponential':
2094	for i in [1,3,5,7,9]:
2095	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
2096	YI += Icoef[i]*Eterm
2097	DYI[i] *= Eterm
2098	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
2099	elif 'Maxwell'in Itype:
2100	Eterm = np.exp(-Icoef[2]/x**2)
2101	DYI[1] = Eterm/x**5
2102	DYI[2] = -Icoef[1]DYI[1]/x*2
2103	YI += (Icoef[1]Eterm/x*5)
2104	if 'Exponential' in Itype:
2105	for i in range(3,11,2):
2106	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
2107	YI += Icoef[i]*Eterm
2108	DYI[i] *= Eterm
2109	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
2110	else: #Chebyschev
2111	T = (2./x)-1.
2112	Ccof = np.ones((12,xdata.shape[0]))
2113	Ccof[1] = T
2114	for i in range(2,12):
2115	Ccof[i] = 2TCcof[i-1]-Ccof[i-2]
2116	for i in range(1,10):
2117	YI += Ccof[i]*Icoef[i+2]
2118	DYI[i+2] =Ccof[i]
2119	return YI,DYI
2120
2121	Iesd = np.array(Iparm['Iesd'])
2122	Icovar = Iparm['Icovar']
2123	YI,DYI = IfunAdv(Iparm,xdata)
2124	YI = np.where(YI>0,YI,1.)
2125	WYI = np.zeros_like(xdata)
2126	vcov = np.zeros((12,12))
2127	k = 0
2128	for i in range(12):
2129	for j in range(i,12):
2130	vcov[i][j] = Icovar[k]Iesd[i]Iesd[j]
2131	vcov[j][i] = Icovar[k]Iesd[i]Iesd[j]
2132	k += 1
2133	M = np.inner(vcov,DYI.T)
2134	WYI = np.sum(M*DYI,axis=0)
2135	WYI = np.where(WYI>0.,WYI,0.)
2136	return YI,WYI
2137
2138	################################################################################
2139	#### RMCutilities
2140	################################################################################
2141
2142	def MakeInst(PWDdata,Name,Size,Mustrain,useSamBrd):
2143	inst = PWDdata['Instrument Parameters'][0]
2144	Xsb = 0.
2145	Ysb = 0.
2146	if 'T' in inst['Type'][1]:
2147	difC = inst['difC'][1]
2148	if useSamBrd[0]:
2149	if 'ellipsoidal' not in Size[0]: #take the isotropic term only
2150	Xsb = 1.e-4*difC/Size[1][0]
2151	if useSamBrd[1]:
2152	if 'generalized' not in Mustrain[0]: #take the isotropic term only
2153	Ysb = 1.e-6difCMustrain[1][0]
2154	prms = ['Bank',
2155	'difC','difA','Zero','2-theta',
2156	'alpha','beta-0','beta-1',
2157	'sig-0','sig-1','sig-2',
2158	'Z','X','Y']
2159	fname = Name+'.inst'
2160	fl = open(fname,'w')
2161	fl.write('1\n')
2162	fl.write('%d\n'%int(inst[prms[0]][1]))
2163	fl.write('%19.11f%19.11f%19.11f%19.11f\n'%(inst[prms[1]][1],inst[prms[2]][1],inst[prms[3]][1],inst[prms[4]][1]))
2164	fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[5]][1],inst[prms[6]][1],inst[prms[7]][1]))
2165	fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[8]][1],inst[prms[9]][1],inst[prms[10]][1]))
2166	fl.write('%12.6e%14.6e%14.6e%14.6e%14.6e\n'%(inst[prms[11]][1],inst[prms[12]][1]+Ysb,inst[prms[13]][1]+Xsb,0.0,0.0))
2167	fl.close()
2168	else:
2169	if useSamBrd[0]:
2170	wave = G2mth.getWave(inst)
2171	if 'ellipsoidal' not in Size[0]: #take the isotropic term only
2172	Xsb = 1.8wave/(np.piSize[1][0])
2173	if useSamBrd[1]:
2174	if 'generalized' not in Mustrain[0]: #take the isotropic term only
2175	Ysb = 0.0180*Mustrain[1][0]/np.pi
2176	prms = ['Bank',
2177	'Lam','Zero','Polariz.',
2178	'U','V','W',
2179	'X','Y']
2180	fname = Name+'.inst'
2181	fl = open(fname,'w')
2182	fl.write('1\n')
2183	fl.write('%d\n'%int(inst[prms[0]][1]))
2184	fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2185	fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2186	fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))
2187	fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2188	fl.close()
2189	return fname
2190
2191	def MakeBack(PWDdata,Name):
2192	Back = PWDdata['Background'][0]
2193	inst = PWDdata['Instrument Parameters'][0]
2194	if 'chebyschev-1' != Back[0]:
2195	return None
2196	Nback = Back[2]
2197	BackVals = Back[3:]
2198	fname = Name+'.back'
2199	fl = open(fname,'w')
2200	fl.write('%10d\n'%Nback)
2201	for val in BackVals:
2202	if 'T' in inst['Type'][1]:
2203	fl.write('%12.6g\n'%(float(val)))
2204	else:
2205	fl.write('%12.6g\n'%val)
2206	fl.close()
2207	return fname
2208
2209	def findDup(Atoms):
2210	Dup = []
2211	Fracs = []
2212	for iat1,at1 in enumerate(Atoms):
2213	if any([at1[0] in dup for dup in Dup]):
2214	continue
2215	else:
2216	Dup.append([at1[0],])
2217	Fracs.append([at1[6],])
2218	for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2219	if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2220	Dup[-1] += [at2[0],]
2221	Fracs[-1]+= [at2[6],]
2222	return Dup,Fracs
2223
2224	def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):
2225
2226	Meta = RMCPdict['metadata']
2227	Atseq = RMCPdict['atSeq']
2228	Supercell = RMCPdict['SuperCell']
2229	generalData = Phase['General']
2230	Dups,Fracs = findDup(Phase['Atoms'])
2231	Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2232	Sample = PWDdata['Sample Parameters']
2233	Meta['temperature'] = Sample['Temperature']
2234	Meta['pressure'] = Sample['Pressure']
2235	Cell = generalData['Cell'][1:7]
2236	Trans = np.eye(3)*np.array(Supercell)
2237	newPhase = copy.deepcopy(Phase)
2238	newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2239	newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2240	GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2241	RMCPdict['Rmax'] = np.min(np.sqrt(np.array([1./G2lat.calc_rDsq2(H,GB) for H in [[1,0,0],[0,1,0],[0,0,1]]])))/2.
2242	newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False)
2243	Natm = np.core.defchararray.count(np.array(Atcodes),'+') #no. atoms in original unit cell
2244	Natm = np.count_nonzero(Natm-1)
2245	Atoms = newPhase['Atoms']
2246	Satoms = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms,5),4),3)
2247	Datoms = [[atom for atom in Satoms if atom[0] in dup] for dup in Dups]
2248	Natoms = []
2249	reset = False
2250	for idup,dup in enumerate(Dups):
2251	ldup = len(dup)
2252	datoms = Datoms[idup]
2253	natm = len(datoms)
2254	i = 0
2255	while i < natm:
2256	atoms = datoms[i:i+ldup]
2257	try:
2258	atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2259	Natoms.append(atom)
2260	except IndexError: #what about vacancies?
2261	if 'Va' not in Atseq:
2262	reset = True
2263	Atseq.append('Va')
2264	RMCPdict['aTypes']['Va'] = 0.0
2265	atom = atoms[0]
2266	atom[1] = 'Va'
2267	Natoms.append(atom)
2268	i += ldup
2269	NAtype = np.zeros(len(Atseq))
2270	for atom in Natoms:
2271	NAtype[Atseq.index(atom[1])] += 1
2272	NAstr = ['%6d'%i for i in NAtype]
2273	Cell = newPhase['General']['Cell'][1:7]
2274	if os.path.exists(Name+'.his6f'):
2275	os.remove(Name+'.his6f')
2276	if os.path.exists(Name+'.neigh'):
2277	os.remove(Name+'.neigh')
2278	fname = Name+'.rmc6f'
2279	fl = open(fname,'w')
2280	fl.write('(Version 6f format configuration file)\n')
2281	for item in Meta:
2282	fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2283	fl.write('Atom types present: %s\n'%' '.join(Atseq))
2284	fl.write('Number of each atom type: %s\n'%''.join(NAstr))
2285	fl.write('Number of atoms: %d\n'%len(Natoms))
2286	fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2287	fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2288	Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2289	A,B = G2lat.cell2AB(Cell,True)
2290	fl.write('Lattice vectors (Ang):\n')
2291	for i in [0,1,2]:
2292	fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2293	fl.write('Atoms (fractional coordinates):\n')
2294	nat = 0
2295	for atm in Atseq:
2296	for iat,atom in enumerate(Natoms):
2297	if atom[1] == atm:
2298	nat += 1
2299	atcode = Atcodes[iat].split(':')
2300	cell = [0,0,0]
2301	if '+' in atcode[1]:
2302	cell = eval(atcode[1].split('+')[1])
2303	fl.write('%6d%4s [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(
2304	nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2305	fl.close()
2306	return fname,reset
2307
2308	def MakeBragg(PWDdata,Name,Phase):
2309	generalData = Phase['General']
2310	Vol = generalData['Cell'][7]
2311	Data = PWDdata['Data']
2312	Inst = PWDdata['Instrument Parameters'][0]
2313	Bank = int(Inst['Bank'][1])
2314	Sample = PWDdata['Sample Parameters']
2315	Scale = Sample['Scale'][0]
2316	if 'X' in Inst['Type'][1]:
2317	Scale *= 2.
2318	Limits = PWDdata['Limits'][1]
2319	Ibeg = np.searchsorted(Data[0],Limits[0])
2320	Ifin = np.searchsorted(Data[0],Limits[1])+1
2321	fname = Name+'.bragg'
2322	fl = open(fname,'w')
2323	fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2324	if 'T' in Inst['Type'][0]:
2325	fl.write('%12s%12s\n'%(' TOF,ms',' I(obs)'))
2326	for i in range(Ibeg,Ifin-1):
2327	fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2328	else:
2329	fl.write('%12s%12s\n'%(' 2-theta, deg',' I(obs)'))
2330	for i in range(Ibeg,Ifin-1):
2331	fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))
2332	fl.close()
2333	return fname
2334
2335	def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2336	Meta = RMCPdict['metadata']
2337	Times = RMCPdict['runTimes']
2338	Atseq = RMCPdict['atSeq']
2339	Atypes = RMCPdict['aTypes']
2340	atPairs = RMCPdict['Pairs']
2341	Files = RMCPdict['files']
2342	BraggWt = RMCPdict['histogram'][1]
2343	inst = PWDdata['Instrument Parameters'][0]
2344	refList = PWDdata['Reflection Lists'][Name]['RefList']
2345	dMin = refList[-1][4]
2346	gsasType = 'xray2'
2347	if 'T' in inst['Type'][1]:
2348	gsasType = 'gsas3'
2349	elif 'X' in inst['Type'][1]:
2350	XFF = G2elem.GetFFtable(Atseq)
2351	Xfl = open(Name+'.xray','w')
2352	for atm in Atseq:
2353	fa = XFF[atm]['fa']
2354	fb = XFF[atm]['fb']
2355	fc = XFF[atm]['fc']
2356	Xfl.write('%2s %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2357	atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2358	Xfl.close()
2359	lenA = len(Atseq)
2360	Pairs = []
2361	for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2362	Pairs += pair
2363	pairMin = [atPairs[pair]for pair in Pairs if pair in atPairs]
2364	maxMoves = [Atypes[atm] for atm in Atseq if atm in Atypes]
2365	fname = Name+'.dat'
2366	fl = open(fname,'w')
2367	fl.write(' %% Hand edit the following as needed\n')
2368	fl.write('TITLE :: '+Name+'\n')
2369	fl.write('MATERIAL :: '+Meta['material']+'\n')
2370	fl.write('PHASE :: '+Meta['phase']+'\n')
2371	fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2372	fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2373	minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2374	minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2375	maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2376	fl.write('MINIMUM_DISTANCES :: %s Angstrom\n'%minHD)
2377	maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2378	fl.write('MAXIMUM_MOVES :: %s Angstrom\n'%maxMv)
2379	fl.write('R_SPACING :: 0.0200 Angstrom\n')
2380	fl.write('PRINT_PERIOD :: 100\n')
2381	fl.write('TIME_LIMIT :: %.2f MINUTES\n'%Times[0])
2382	fl.write('SAVE_PERIOD :: %.2f MINUTES\n'%Times[1])
2383	fl.write('\n')
2384	fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2385	fl.write('\n')
2386	fl.write('FLAGS ::\n')
2387	fl.write(' > NO_MOVEOUT\n')
2388	fl.write(' > NO_SAVE_CONFIGURATIONS\n')
2389	fl.write(' > NO_RESOLUTION_CONVOLUTION\n')
2390	fl.write('\n')
2391	fl.write('INPUT_CONFIGURATION_FORMAT :: rmc6f\n')
2392	fl.write('SAVE_CONFIGURATION_FORMAT :: rmc6f\n')
2393	fl.write('IGNORE_HISTORY_FILE ::\n')
2394	fl.write('\n')
2395	fl.write('DISTANCE_WINDOW ::\n')
2396	fl.write(' > MNDIST :: %s\n'%minD)
2397	fl.write(' > MXDIST :: %s\n'%maxD)
2398	if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2399	fl.write('\n')
2400	fl.write('POTENTIALS ::\n')
2401	fl.write(' > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2402	fl.write(' > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2403	if len(RMCPdict['Potentials']['Stretch']):
2404	fl.write(' > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2405	for bond in RMCPdict['Potentials']['Stretch']:
2406	fl.write(' > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))
2407	if len(RMCPdict['Potentials']['Angles']):
2408	fl.write(' > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2409	for angle in RMCPdict['Potentials']['Angles']:
2410	fl.write(' > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2411	(angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2412	if RMCPdict['useBVS']:
2413	fl.write('BVS ::\n')
2414	fl.write(' > ATOM :: '+' '.join(Atseq)+'\n')
2415	fl.write(' > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2416	oxid = []
2417	for val in RMCPdict['Oxid']:
2418	if len(val) == 3:
2419	oxid.append(val[0][1:])
2420	else:
2421	oxid.append(val[0][2:])
2422	fl.write(' > OXID :: %s\n'%' '.join(oxid))
2423	fl.write(' > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2424	fl.write(' > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2425	fl.write(' > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2426	fl.write(' > SAVE :: 100000\n')
2427	fl.write(' > UPDATE :: 100000\n')
2428	if len(RMCPdict['Swap']):
2429	fl.write('\n')
2430	fl.write('SWAP_MULTI ::\n')
2431	for swap in RMCPdict['Swap']:
2432	try:
2433	at1 = Atseq.index(swap[0])
2434	at2 = Atseq.index(swap[1])
2435	except ValueError:
2436	break
2437	fl.write(' > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2438
2439	if len(RMCPdict['FxCN']):
2440	fl.write('FIXED_COORDINATION_CONSTRAINTS :: %d\n'%len(RMCPdict['FxCN']))
2441	for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2442	try:
2443	at1 = Atseq.index(fxcn[0])
2444	at2 = Atseq.index(fxcn[1])
2445	except ValueError:
2446	break
2447	fl.write(' > CSTR%d :: %d %d %.2f %.2f %.2f %.2f %.6f\n'%(ifx+1,at1+1,at2+1,fxcn[2],fxcn[3],fxcn[4],fxcn[5],fxcn[6]))
2448	if len(RMCPdict['AveCN']):
2449	fl.write('AVERAGE_COORDINATION_CONSTRAINTS :: %d\n'%len(RMCPdict['AveCN']))
2450	for iav,avcn in enumerate(RMCPdict['AveCN']):
2451	try:
2452	at1 = Atseq.index(avcn[0])
2453	at2 = Atseq.index(avcn[1])
2454	except ValueError:
2455	break
2456	fl.write(' > CAVSTR%d :: %d %d %.2f %.2f %.2f %.6f\n'%(iav+1,at1+1,at2+1,avcn[2],avcn[3],avcn[4],avcn[5]))
2457	for File in Files:
2458	if Files[File][0] and Files[File][0] != 'Select':
2459	if 'Xray' in File and 'F(Q)' in File:
2460	fqdata = open(Files[File][0],'r')
2461	lines = int(fqdata.readline()[:-1])
2462	fl.write('\n')
2463	fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2464	fl.write(' > FILENAME :: %s\n'%Files[File][0])
2465	fl.write(' > DATA_TYPE :: %s\n'%Files[File][2])
2466	fl.write(' > FIT_TYPE :: %s\n'%Files[File][2])
2467	if 'Xray' not in File:
2468	fl.write(' > START_POINT :: 1\n')
2469	fl.write(' > END_POINT :: 3000\n')
2470	fl.write(' > WEIGHT :: %.4f\n'%Files[File][1])
2471	fl.write(' > CONSTANT_OFFSET 0.000\n')
2472	fl.write(' > NO_FITTED_OFFSET\n')
2473	if RMCPdict['FitScale']:
2474	fl.write(' > FITTED_SCALE\n')
2475	else:
2476	fl.write(' > NO_FITTED_SCALE\n')
2477	if Files[File][3] !='RMC':
2478	fl.write(' > %s\n'%Files[File][3])
2479	if 'reciprocal' in File:
2480	fl.write(' > CONVOLVE ::\n')
2481	if 'Xray' in File:
2482	fl.write(' > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2483	fl.write(' > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2484	fl.write(' > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2485	fl.write(' > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
2486	fl.write('\n')
2487	fl.write('BRAGG ::\n')
2488	fl.write(' > BRAGG_SHAPE :: %s\n'%gsasType)
2489	fl.write(' > RECALCUATE\n')
2490	fl.write(' > DMIN :: %.2f\n'%(dMin-0.02))
2491	fl.write(' > WEIGHT :: %10.3f\n'%BraggWt)
2492	fl.write('\n')
2493	fl.write('END ::\n')
2494	fl.close()
2495	return fname
2496
2497	def FindBonds(Phase,RMCPdict):
2498	generalData = Phase['General']
2499	cx,ct,cs,cia = generalData['AtomPtrs']
2500	atomData = Phase['Atoms']
2501	Res = 'RMC'
2502	if 'macro' in generalData['Type']:
2503	Res = atomData[0][ct-3]
2504	AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2505	Pairs = RMCPdict['Pairs'] #dict!
2506	BondList = []
2507	notNames = []
2508	for FrstName in AtDict:
2509	nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
2510	Atyp1 = AtDict[FrstName]
2511	if 'Va' in Atyp1:
2512	continue
2513	for nbr in nbrs:
2514	Atyp2 = AtDict[nbr[0]]
2515	if 'Va' in Atyp2:
2516	continue
2517	try:
2518	bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
2519	except KeyError:
2520	bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
2521	if any(bndData):
2522	if bndData[0] <= nbr[1] <= bndData[1]:
2523	bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
2524	revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
2525	if bondStr not in BondList and revbondStr not in BondList:
2526	BondList.append(bondStr)
2527	notNames.append(FrstName)
2528	return Res,BondList
2529
2530	def FindAngles(Phase,RMCPdict):
2531	generalData = Phase['General']
2532	Cell = generalData['Cell'][1:7]
2533	Amat = G2lat.cell2AB(Cell)[0]
2534	cx,ct,cs,cia = generalData['AtomPtrs']
2535	atomData = Phase['Atoms']
2536	AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
2537	AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2538	Angles = RMCPdict['Angles']
2539	AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
2540	AngleList = []
2541	for MidName in AtDict:
2542	nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
2543	if len(nbrs) < 2: #need 2 neighbors to make an angle
2544	continue
2545	Atyp2 = AtDict[MidName]
2546	for i,nbr1 in enumerate(nbrs):
2547	Atyp1 = AtDict[nbr1[0]]
2548	for j,nbr3 in enumerate(nbrs[i+1:]):
2549	Atyp3 = AtDict[nbr3[0]]
2550	IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
2551	try:
2552	angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
2553	except KeyError:
2554	try:
2555	angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
2556	except KeyError:
2557	continue
2558	XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
2559	calAngle = G2mth.getRestAngle(XYZ,Amat)
2560	if angData[0] < calAngle < angData[1]:
2561	AngleList.append(str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n')
2562	return AngleList
2563
2564	def GetSqConvolution(XY,d):
2565
2566	n = XY.shape[1]
2567	snew = np.zeros(n)
2568	dq = np.zeros(n)
2569	sold = XY[1]
2570	q = XY[0]
2571	dq[1:] = np.diff(q)
2572	dq[0] = dq[1]
2573
2574	for j in range(n):
2575	for i in range(n):
2576	b = abs(q[i]-q[j])
2577	t = q[i]+q[j]
2578	if j == i:
2579	snew[j] += q[i]sold[i](d-np.sin(td)/t)dq[i]
2580	else:
2581	snew[j] += q[i]sold[i](np.sin(bd)/b-np.sin(td)/t)*dq[i]
2582	snew[j] /= np.pi*q[j]
2583
2584	snew[0] = snew[1]
2585	return snew
2586
2587	def GetMaxSphere(pdbName):
2588	try:
2589	pFil = open(pdbName,'r')
2590	except FileNotFoundError:
2591	return None
2592	while True:
2593	line = pFil.readline()
2594	if 'Boundary' in line:
2595	line = line.split()[3:]
2596	G = np.array([float(item) for item in line])
2597	G = np.reshape(G,(3,3))**2
2598	G = nl.inv(G)
2599	pFil.close()
2600	break
2601	dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
2602	return min(dspaces)
2603
2604	def MakefullrmcRun(pName,Phase,RMCPdict):
2605	Res,BondList = FindBonds(Phase,RMCPdict)
2606	AngleList = FindAngles(Phase,RMCPdict)
2607	rmin = RMCPdict['min Contact']
2608	rmax = GetMaxSphere(RMCPdict['atomPDB'])
2609	if rmax is None:
2610	return None
2611	rname = pName+'-run.py'
2612	restart = '%s_restart.pdb'%pName
2613	Files = RMCPdict['files']
2614	wtDict = {}
2615	rundata = ''
2616	rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%rname
2617	rundata += '''
2618	# fullrmc imports (all that are potentially useful)
2619	import numpy as np
2620	import time
2621	from fullrmc.sincConvolution import sincConvolution
2622	from fullrmc.Globals import LOGGER
2623	from fullrmc.Engine import Engine
2624	from fullrmc.Constraints.PairDistributionConstraints import PairDistributionConstraint
2625	from fullrmc.Constraints.PairCorrelationConstraints import PairCorrelationConstraint
2626	from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint
2627	from fullrmc.Constraints.DistanceConstraints import InterMolecularDistanceConstraint
2628	from fullrmc.Constraints.BondConstraints import BondConstraint
2629	from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
2630	from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
2631	from fullrmc.Core.MoveGenerator import MoveGeneratorCollector
2632	from fullrmc.Generators.Swaps import SwapPositionsGenerator
2633	from fullrmc.Core.GroupSelector import RecursiveGroupSelector
2634	from fullrmc.Selectors.RandomSelectors import RandomSelector
2635	from fullrmc.Selectors.OrderedSelectors import DefinedOrderSelector
2636	from fullrmc.Generators.Translations import TranslationGenerator, TranslationAlongSymmetryAxisGenerator
2637	from fullrmc.Generators.Agitations import DistanceAgitationGenerator, AngleAgitationGenerator
2638	from fullrmc.Generators.Rotations import RotationGenerator, RotationAboutAxisGenerator
2639	from fullrmc.Core.Collection import get_principal_axis
2640	from fullrmc.debugStuff import *
2641	InvokeDebugOpts()
2642	time0 = time.time()
2643	SwapGen = {}
2644	# engine setup\n'''
2645	rundata += 'LOGGER.set_log_file_basename("%s")\n'%pName
2646	rundata += 'engineFileName = "%s.rmc"\n'%pName
2647	rundata += 'ENGINE = Engine(path=None)\n'
2648	rundata += 'if not ENGINE.is_engine(engineFileName):\n'
2649	rundata += '# create engine & set atomic (pdb) model\n'
2650	rundata += ' ENGINE = Engine(path=engineFileName)\n'
2651	rundata += ' ENGINE.set_pdb("%s")\n'%RMCPdict['atomPDB']
2652	rundata += '# create experimental constraints must restart to change these\n'
2653	for File in Files:
2654	filDat = RMCPdict['files'][File]
2655	if filDat[0] != 'Select':
2656	sfwt = 'neutronCohb'
2657	if 'Xray' in File:
2658	sfwt = 'atomicNumber'
2659	if 'G(r)' in File:
2660	rundata += ' RGR = np.loadtxt("%s").T\n'%filDat[0]
2661	if filDat[3]:
2662	rundata += ' RGR[1] *= RGR[0]\n'
2663	rundata += ' GofR = PairDistributionConstraint(experimentalData=RGR.T, weighting="%s")\n'%sfwt
2664	rundata += ' ENGINE.add_constraints([GofR])\n'
2665	wtDict['Pair-'+sfwt] = filDat[1]
2666	else:
2667	rundata += ' SOQ = np.loadtxt("%s").T\n'%filDat[0]
2668	if filDat[3]:
2669	rundata += ' SOQ[1] = sincConvolution(SOQ,%.3f)\n'%rmax
2670	rundata += ' FofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
2671	rundata += ' ENGINE.add_constraints([FofQ])\n'
2672	wtDict['Struct-'+sfwt] = filDat[1]
2673	rundata += ' ENGINE.add_constraints(InterMolecularDistanceConstraint())\n'
2674	if RMCPdict['byMolec']:
2675	if len(BondList):
2676	rundata += ' B_CONSTRAINT = BondConstraint()\n'
2677	rundata += ' ENGINE.add_constraints(B_CONSTRAINT)\n'
2678	if len(AngleList):
2679	rundata += ' A_CONSTRAINT = BondsAngleConstraint()\n'
2680	rundata += ' ENGINE.add_constraints(A_CONSTRAINT)\n'
2681	if len(RMCPdict['Torsions']):
2682	rundata += ' T_CONSTRAINT = DihedralAngleConstraint()\n'
2683	rundata += ' ENGINE.add_constraints(T_CONSTRAINT)\n'
2684	rundata += ' ENGINE.save()\n'
2685	rundata += 'else:\n'
2686	rundata += ' ENGINE = ENGINE.load(path=engineFileName)\n'
2687	rundata += '#fill & change constraints - can be done without restart\n'
2688	rundata += 'wtDict = %s\n'%str(wtDict)
2689	rundata += 'Constraints = ENGINE.constraints\n'
2690	rundata += 'for constraint in Constraints:\n'
2691	rundata += ' strcons = str(type(constraint))\n'
2692	rundata += ' if "InterMolecular" in strcons:\n'
2693	rundata += ' constraint.set_default_distance(%f)\n'%RMCPdict['min Contact']
2694	rundata += ' elif "PairDistribution" in strcons:\n'
2695	rundata += ' constraint.set_variance_squared(wtDict["Pair-"+constraint.weighting])\n'
2696	rundata += ' constraint.set_limits((%.3f,%.3f))\n'%(rmin,rmax)
2697	if RMCPdict['FitScale']:
2698	rundata += ' constraint.set_adjust_scale_factor((10, 0.01, 100.))\n'
2699	rundata += ' elif "StructureFactor" in strcons:\n'
2700	rundata += ' constraint.set_variance_squared(wtDict["Struct-"+constraint.weighting])\n'
2701	if RMCPdict['FitScale']:
2702	rundata += ' constraint.set_adjust_scale_factor((10, 0.01, 100.))\n'
2703	if RMCPdict['byMolec']:
2704	if len(BondList):
2705	rundata += ' elif "BondConstraint" in strcons:\n'
2706	rundata += ' constraint.set_variance_squared(%f)\n'%RMCPdict['Bond Weight']
2707	rundata += ' constraint.create_bonds_by_definition(bondsDefinition={"%s":[\n'%Res
2708	for bond in BondList:
2709	rundata += ' %s'%bond
2710	rundata += ' ]})\n'
2711	if len(AngleList):
2712	rundata += ' elif "BondsAngleConstraint" in strcons:\n'
2713	rundata += ' constraint.set_variance_squared(%f)\n'%RMCPdict['Angle Weight']
2714	rundata += ' constraint.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
2715	for angle in AngleList:
2716	rundata += ' %s'%angle
2717	rundata += ' ]})\n'
2718	if len(RMCPdict['Torsions']):
2719	rundata += ' elif "DihedralAngleConstraint" in strcons:\n'
2720	rundata += ' constraint.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
2721	rundata += ' constraint.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
2722	for torsion in RMCPdict['Torsions']:
2723	rundata += ' %s\n'%str(tuple(torsion))
2724	rundata += ' ]})\n'
2725	if len(RMCPdict['Swaps']):
2726	rundata += ' allNames = ENGINE.allNames\n'
2727	for swap in RMCPdict['Swaps']:
2728	rundata += ' SwapA = [[idx] for idx in range(len(allNames)) if allNames[idx]=="%s"]\n'%swap[0]
2729	rundata += ' SwapB = [[idx] for idx in range(len(allNames)) if allNames[idx]=="%s"]\n'%swap[1]
2730	rundata += ' SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
2731	rundata += 'ENGINE.save()\n'
2732	rundata += '#setup runs for fullrmc\n'
2733
2734	rundata += 'for _ in range(%d):\n'%RMCPdict['Cycles']
2735	rundata += ' ENGINE.set_groups_as_atoms()\n'
2736	rundata += ' ENGINE.run(restartPdb="%s",numberOfSteps=10000, saveFrequency=1000)\n'%restart
2737	if len(RMCPdict['Swaps']):
2738	rundata += ' for swaps in SwapGen:\n'
2739	rundata += ' AB = swaps.split("-")\n'
2740	rundata += ' ENGINE.set_groups_as_atoms()\n'
2741	rundata += ' for g in ENGINE.groups:\n'
2742	rundata += ' if allNames[g.indexes[0]]==AB[0]:\n'
2743	rundata += ' g.set_move_generator(SwapGen[swaps][0])\n'
2744	rundata += ' elif allNames[g.indexes[0]]==AB[1]:\n'
2745	rundata += ' g.set_move_generator(SwapGen[swaps][1])\n'
2746	rundata += ' sProb = SwapGen[swaps][2]\n'
2747	rundata += ' ENGINE.run(restartPdb="%s",numberOfSteps=10000*sProb, saveFrequency=1000)\n'%restart
2748	rundata += ' ENGINE.set_groups_as_atoms()\n'
2749	rundata += ' ENGINE.run(restartPdb="%s",numberOfSteps=10000*(1.-sProb), saveFrequency=1000)\n'%restart
2750	rundata += 'ENGINE.close()\n'
2751	rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
2752	rfile = open(rname,'w')
2753	rfile.writelines(rundata)
2754	rfile.close()
2755
2756	return rname
2757
2758
2759	def MakefullrmcPDB(Name,Phase,RMCPdict):
2760	generalData = Phase['General']
2761	Atseq = RMCPdict['atSeq']
2762	Dups,Fracs = findDup(Phase['Atoms'])
2763	Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2764	ifSfracs = np.any(np.array(Sfracs)-1.)
2765	Supercell = RMCPdict['SuperCell']
2766	Cell = generalData['Cell'][1:7]
2767	Trans = np.eye(3)*np.array(Supercell)
2768	newPhase = copy.deepcopy(Phase)
2769	newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2770	newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
2771	newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=False)
2772	Atoms = newPhase['Atoms']
2773
2774	if ifSfracs:
2775	Natm = np.core.defchararray.count(np.array(Atcodes),'+') #no. atoms in original unit cell
2776	Natm = np.count_nonzero(Natm-1)
2777	Satoms = []
2778	for i in range(len(Atoms)//Natm):
2779	ind = i*Natm
2780	Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
2781	Natoms = []
2782	for satoms in Satoms:
2783	for idup,dup in enumerate(Dups):
2784	ldup = len(dup)
2785	natm = len(satoms)
2786	i = 0
2787	while i < natm:
2788	if satoms[i][0] in dup:
2789	atoms = satoms[i:i+ldup]
2790	try:
2791	atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2792	Natoms.append(atom)
2793	except IndexError: #what about vacancies?
2794	if 'Va' not in Atseq:
2795	Atseq.append('Va')
2796	RMCPdict['aTypes']['Va'] = 0.0
2797	atom = atoms[0]
2798	atom[1] = 'Va'
2799	Natoms.append(atom)
2800	i += ldup
2801	else:
2802	i += 1
2803	else:
2804	Natoms = Atoms
2805
2806	XYZ = np.array([atom[3:6] for atom in Natoms]).T
2807	XYZptp = np.array([ma.ptp(XYZ[0]),ma.ptp(XYZ[1]),ma.ptp(XYZ[2])])/2.
2808	Cell = newPhase['General']['Cell'][1:7]
2809	A,B = G2lat. cell2AB(Cell)
2810	fname = Name+'_cbb.pdb'
2811	fl = open(fname,'w')
2812	fl.write('REMARK Boundary Conditions:%6.2f 0.0 0.0 0.0%7.2f 0.0 0.0 0.0%7.2f\n'%(
2813	Cell[0],Cell[1],Cell[2]))
2814	fl.write('ORIGX1 1.000000 0.000000 0.000000 0.00000\n')
2815	fl.write('ORIGX2 0.000000 1.000000 0.000000 0.00000\n')
2816	fl.write('ORIGX3 0.000000 0.000000 1.000000 0.00000\n')
2817	fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1 1\n'%(
2818	Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2819
2820	Natm = np.core.defchararray.count(np.array(Atcodes),'+')
2821	Natm = np.count_nonzero(Natm-1)
2822	nat = 0
2823	if RMCPdict['byMolec']:
2824	NPM = RMCPdict['Natoms']
2825	for iat,atom in enumerate(Natoms):
2826	XYZ = np.inner(A,np.array(atom[3:6])-XYZptp) #shift origin to middle & make Cartesian;residue = 'RMC'
2827	fl.write('ATOM %5d %-4s RMC%6d%12.3f%8.3f%8.3f 1.00 0.00 %2s\n'%(
2828	1+nat%NPM,atom[0],1+nat//NPM,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
2829	nat += 1
2830	else:
2831	for atm in Atseq:
2832	for iat,atom in enumerate(Natoms):
2833	if atom[1] == atm:
2834	XYZ = np.inner(A,np.array(atom[3:6])-XYZptp) #shift origin to middle & make Cartesian
2835	fl.write('ATOM %5d %-4s RMC%6d%12.3f%8.3f%8.3f 1.00 0.00 %2s\n'%(
2836	1+nat,atom[0],1+nat,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
2837	nat += 1
2838	fl.close()
2839	return fname
2840
2841	def MakePdparse(RMCPdict):
2842	fname = 'make_pdb.py'
2843	outName = RMCPdict['moleculePdb'].split('.')
2844	outName[0] += '_rbb'
2845	outName = '.'.join(outName)
2846	RMCPdict['atomPDB'] = outName #might be empty if pdbparser run fails
2847	fl = open(fname,'w')
2848	fl.write('from pdbparser.pdbparser import pdbparser\n')
2849	fl.write('from pdbparser.Utilities.Construct import AmorphousSystem\n')
2850	fl.write("pdb = pdbparser('%s')\n"%RMCPdict['moleculePdb'])
2851	boxstr= 'boxsize=%s'%str(RMCPdict['Box'])
2852	recstr = 'recursionLimit=%d'%RMCPdict['maxRecursion']
2853	denstr = 'density=%.3f'%RMCPdict['targetDensity']
2854	fl.write('pdb = AmorphousSystem(pdb,%s,%s,%s,\n'%(boxstr,recstr,denstr))
2855	fl.write(' priorities={"boxSize":True, "insertionNumber":False, "density":True}).construct().get_pdb()\n')
2856	fl.write('pdb.export_pdb("%s")\n'%outName)
2857	fl.close
2858	return fname
2859
2860	def GetRMCBonds(general,RMCPdict,Atoms,bondList):
2861	bondDist = []
2862	Cell = general['Cell'][1:7]
2863	Supercell = RMCPdict['SuperCell']
2864	Trans = np.eye(3)*np.array(Supercell)
2865	Cell = G2lat.TransformCell(Cell,Trans)[:6]
2866	Amat,Bmat = G2lat.cell2AB(Cell)
2867	indices = (-1,0,1)
2868	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2869	for bonds in bondList:
2870	Oxyz = np.array(Atoms[bonds[0]][1:])
2871	Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])
2872	Dx = np.array([Txyz-Oxyz+unit for unit in Units])
2873	Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
2874	for dx in Dx.T:
2875	bondDist.append(np.min(dx))
2876	return np.array(bondDist)
2877
2878	def GetRMCAngles(general,RMCPdict,Atoms,angleList):
2879	bondAngles = []
2880	Cell = general['Cell'][1:7]
2881	Supercell = RMCPdict['SuperCell']
2882	Trans = np.eye(3)*np.array(Supercell)
2883	Cell = G2lat.TransformCell(Cell,Trans)[:6]
2884	Amat,Bmat = G2lat.cell2AB(Cell)
2885	indices = (-1,0,1)
2886	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2887	for angle in angleList:
2888	Oxyz = np.array(Atoms[angle[0]][1:])
2889	TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])
2890	TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])
2891	DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
2892	DAx = np.sqrt(np.sum(DAxV**2,axis=2))
2893	DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
2894	DBx = np.sqrt(np.sum(DBxV**2,axis=2))
2895	iDAx = np.argmin(DAx,axis=0)
2896	iDBx = np.argmin(DBx,axis=0)
2897	for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
2898	DAv = DAxV[iA,i]/DAx[iA,i]
2899	DBv = DBxV[iB,i]/DBx[iB,i]
2900	bondAngles.append(npacosd(np.sum(DAv*DBv)))
2901	return np.array(bondAngles)
2902
2903	################################################################################
2904	#### Reflectometry calculations
2905	################################################################################
2906
2907	def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
2908	G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
2909
2910	class RandomDisplacementBounds(object):
2911	"""random displacement with bounds"""
2912	def __init__(self, xmin, xmax, stepsize=0.5):
2913	self.xmin = xmin
2914	self.xmax = xmax
2915	self.stepsize = stepsize
2916
2917	def __call__(self, x):
2918	"""take a random step but ensure the new position is within the bounds"""
2919	while True:
2920	# this could be done in a much more clever way, but it will work for example purposes
2921	steps = self.xmax-self.xmin
2922	xnew = x + np.random.uniform(-self.stepsizesteps, self.stepsizesteps, np.shape(x))
2923	if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
2924	break
2925	return xnew
2926
2927	def GetModelParms():
2928	parmDict = {}
2929	varyList = []
2930	values = []
2931	bounds = []
2932	parmDict['dQ type'] = data['dQ type']
2933	parmDict['Res'] = data['Resolution'][0]/(100.*sateln2) #% FWHM-->decimal sig
2934	for parm in ['Scale','FltBack']:
2935	parmDict[parm] = data[parm][0]
2936	if data[parm][1]:
2937	varyList.append(parm)
2938	values.append(data[parm][0])
2939	bounds.append(Bounds[parm])
2940	parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
2941	parmDict['nLayers'] = len(parmDict['Layer Seq'])
2942	for ilay,layer in enumerate(data['Layers']):
2943	name = layer['Name']
2944	cid = str(ilay)+';'
2945	parmDict[cid+'Name'] = name
2946	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2947	parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
2948	if layer.get(parm,[0.,False])[1]:
2949	varyList.append(cid+parm)
2950	value = layer[parm][0]
2951	values.append(value)
2952	if value:
2953	bound = [value*Bfac,value/Bfac]
2954	else:
2955	bound = [0.,10.]
2956	bounds.append(bound)
2957	if name not in ['vacuum','unit scatter']:
2958	parmDict[cid+'rho'] = Substances[name]['Scatt density']
2959	parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
2960	return parmDict,varyList,values,bounds
2961
2962	def SetModelParms():
2963	line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
2964	if 'Scale' in varyList:
2965	data['Scale'][0] = parmDict['Scale']
2966	line += ' esd: %.4g'%(sigDict['Scale'])
2967	G2fil.G2Print (line)
2968	line = ' Flat background: %15.4g'%(parmDict['FltBack'])
2969	if 'FltBack' in varyList:
2970	data['FltBack'][0] = parmDict['FltBack']
2971	line += ' esd: %15.3g'%(sigDict['FltBack'])
2972	G2fil.G2Print (line)
2973	for ilay,layer in enumerate(data['Layers']):
2974	name = layer['Name']
2975	G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
2976	cid = str(ilay)+';'
2977	line = ' '
2978	line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
2979	line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
2980	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
2981	if parm in layer:
2982	if parm == 'Rough':
2983	layer[parm][0] = abs(parmDict[cid+parm]) #make positive
2984	else:
2985	layer[parm][0] = parmDict[cid+parm]
2986	line += ' %s: %.3f'%(parm,layer[parm][0])
2987	if cid+parm in varyList:
2988	line += ' esd: %.3g'%(sigDict[cid+parm])
2989	G2fil.G2Print (line)
2990	G2fil.G2Print (line2)
2991
2992	def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2993	parmDict.update(zip(varyList,values))
2994	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
2995	return M
2996
2997	def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
2998	parmDict.update(zip(varyList,values))
2999	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3000	return np.sum(M**2)
3001
3002	def getREFD(Q,Qsig,parmDict):
3003	Ic = np.ones_like(Q)*parmDict['FltBack']
3004	Scale = parmDict['Scale']
3005	Nlayers = parmDict['nLayers']
3006	Res = parmDict['Res']
3007	depth = np.zeros(Nlayers)
3008	rho = np.zeros(Nlayers)
3009	irho = np.zeros(Nlayers)
3010	sigma = np.zeros(Nlayers)
3011	for ilay,lay in enumerate(parmDict['Layer Seq']):
3012	cid = str(lay)+';'
3013	depth[ilay] = parmDict[cid+'Thick']
3014	sigma[ilay] = parmDict[cid+'Rough']
3015	if parmDict[cid+'Name'] == u'unit scatter':
3016	rho[ilay] = parmDict[cid+'DenMul']
3017	irho[ilay] = parmDict[cid+'iDenMul']
3018	elif 'vacuum' != parmDict[cid+'Name']:
3019	rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
3020	irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
3021	if cid+'Mag SLD' in parmDict:
3022	rho[ilay] += parmDict[cid+'Mag SLD']
3023	if parmDict['dQ type'] == 'None':
3024	AB = abeles(0.5*Q,depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
3025	elif 'const' in parmDict['dQ type']:
3026	AB = SmearAbeles(0.5Q,QRes,depth,rho,irho,sigma[1:])
3027	else: #dQ/Q in data
3028	AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
3029	Ic += AB*Scale
3030	return Ic
3031
3032	def estimateT0(takestep):
3033	Mmax = -1.e-10
3034	Mmin = 1.e10
3035	for i in range(100):
3036	x0 = takestep(values)
3037	M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3038	Mmin = min(M,Mmin)
3039	MMax = max(M,Mmax)
3040	return 1.5*(MMax-Mmin)
3041
3042	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3043	if data.get('2% weight'):
3044	wt = 1./(0.02Io)*2
3045	Qmin = Limits[1][0]
3046	Qmax = Limits[1][1]
3047	wtFactor = ProfDict['wtFactor']
3048	Bfac = data['Toler']
3049	Ibeg = np.searchsorted(Q,Qmin)
3050	Ifin = np.searchsorted(Q,Qmax)+1 #include last point
3051	Ic[:] = 0
3052	Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
3053	'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
3054	parmDict,varyList,values,bounds = GetModelParms()
3055	Msg = 'Failed to converge'
3056	if varyList:
3057	if data['Minimizer'] == 'LMLS':
3058	result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
3059	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3060	parmDict.update(zip(varyList,result[0]))
3061	chisq = np.sum(result[2]['fvec']**2)
3062	ncalc = result[2]['nfev']
3063	covM = result[1]
3064	newVals = result[0]
3065	elif data['Minimizer'] == 'Basin Hopping':
3066	xyrng = np.array(bounds).T
3067	take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
3068	T0 = estimateT0(take_step)
3069	G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
3070	result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
3071	interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
3072	'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
3073	chisq = result.fun
3074	ncalc = result.nfev
3075	newVals = result.x
3076	covM = []
3077	elif data['Minimizer'] == 'MC/SA Anneal':
3078	xyrng = np.array(bounds).T
3079	result = G2mth.anneal(sumREFD, values,
3080	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
3081	schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
3082	boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
3083	ranRange=0.20,autoRan=False,dlg=None)
3084	newVals = result[0]
3085	parmDict.update(zip(varyList,newVals))
3086	chisq = result[1]
3087	ncalc = result[3]
3088	covM = []
3089	G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
3090	elif data['Minimizer'] == 'L-BFGS-B':
3091	result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds, #ftol=Ftol,
3092	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3093	parmDict.update(zip(varyList,result['x']))
3094	chisq = result.fun
3095	ncalc = result.nfev
3096	newVals = result.x
3097	covM = []
3098	else: #nothing varied
3099	M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3100	chisq = np.sum(M**2)
3101	ncalc = 0
3102	covM = []
3103	sig = []
3104	sigDict = {}
3105	result = []
3106	Rvals = {}
3107	Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]Io[Ibeg:Ifin]2))100. #to %
3108	Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList)) #reduced chi^2
3109	Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
3110	Ib[Ibeg:Ifin] = parmDict['FltBack']
3111	try:
3112	if not len(varyList):
3113	Msg += ' - nothing refined'
3114	raise ValueError
3115	Nans = np.isnan(newVals)
3116	if np.any(Nans):
3117	name = varyList[Nans.nonzero(True)[0]]
3118	Msg += ' Nan result for '+name+'!'
3119	raise ValueError
3120	Negs = np.less_equal(newVals,0.)
3121	if np.any(Negs):
3122	indx = Negs.nonzero()
3123	name = varyList[indx[0][0]]
3124	if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
3125	Msg += ' negative coefficient for '+name+'!'
3126	raise ValueError
3127	if len(covM):
3128	sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
3129	covMatrix = covM*Rvals['GOF']
3130	else:
3131	sig = np.zeros(len(varyList))
3132	covMatrix = []
3133	sigDict = dict(zip(varyList,sig))
3134	G2fil.G2Print (' Results of reflectometry data modelling fit:')
3135	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
3136	G2fil.G2Print ('Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
3137	SetModelParms()
3138	return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
3139	except (ValueError,TypeError): #when bad LS refinement; covM missing or with nans
3140	G2fil.G2Print (Msg)
3141	return False,0,0,0,0,0,0,Msg
3142
3143	def makeSLDprofile(data,Substances):
3144
3145	sq2 = np.sqrt(2.)
3146	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3147	Nlayers = len(laySeq)
3148	laySeq = np.array(laySeq,dtype=int)
3149	interfaces = np.zeros(Nlayers)
3150	rho = np.zeros(Nlayers)
3151	sigma = np.zeros(Nlayers)
3152	name = data['Layers'][0]['Name']
3153	thick = 0.
3154	for ilay,lay in enumerate(laySeq):
3155	layer = data['Layers'][lay]
3156	name = layer['Name']
3157	if 'Thick' in layer:
3158	thick += layer['Thick'][0]
3159	interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
3160	if 'Rough' in layer:
3161	sigma[ilay] = max(0.001,layer['Rough'][0])
3162	if name != 'vacuum':
3163	if name == 'unit scatter':
3164	rho[ilay] = np.sqrt(layer['DenMul'][0]2+layer['iDenMul'][0]2)
3165	else:
3166	rrho = Substances[name]['Scatt density']
3167	irho = Substances[name]['XImag density']
3168	rho[ilay] = np.sqrt(rrho2+irho2)*layer['DenMul'][0]
3169	if 'Mag SLD' in layer:
3170	rho[ilay] += layer['Mag SLD'][0]
3171	name = data['Layers'][-1]['Name']
3172	x = np.linspace(-0.15thick,1.15thick,1000,endpoint=True)
3173	xr = np.flipud(x)
3174	interfaces[-1] = x[-1]
3175	y = np.ones_like(x)*rho[0]
3176	iBeg = 0
3177	for ilayer in range(Nlayers-1):
3178	delt = rho[ilayer+1]-rho[ilayer]
3179	iPos = np.searchsorted(x,interfaces[ilayer])
3180	y[iBeg:] += (delt/2.)sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2sigma[ilayer+1]))
3181	iBeg = iPos
3182	return x,xr,y
3183
3184	def REFDModelFxn(Profile,Inst,Limits,Substances,data):
3185
3186	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3187	Qmin = Limits[1][0]
3188	Qmax = Limits[1][1]
3189	iBeg = np.searchsorted(Q,Qmin)
3190	iFin = np.searchsorted(Q,Qmax)+1 #include last point
3191	Ib[:] = data['FltBack'][0]
3192	Ic[:] = 0
3193	Scale = data['Scale'][0]
3194	if data['Layer Seq'] == []:
3195	return
3196	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3197	Nlayers = len(laySeq)
3198	depth = np.zeros(Nlayers)
3199	rho = np.zeros(Nlayers)
3200	irho = np.zeros(Nlayers)
3201	sigma = np.zeros(Nlayers)
3202	for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
3203	layer = data['Layers'][lay]
3204	name = layer['Name']
3205	if 'Thick' in layer: #skips first & last layers
3206	depth[ilay] = layer['Thick'][0]
3207	if 'Rough' in layer: #skips first layer
3208	sigma[ilay] = layer['Rough'][0]
3209	if 'unit scatter' == name:
3210	rho[ilay] = layer['DenMul'][0]
3211	irho[ilay] = layer['iDenMul'][0]
3212	else:
3213	rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
3214	irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
3215	if 'Mag SLD' in layer:
3216	rho[ilay] += layer['Mag SLD'][0]
3217	if data['dQ type'] == 'None':
3218	AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
3219	elif 'const' in data['dQ type']:
3220	res = data['Resolution'][0]/(100.*sateln2)
3221	AB = SmearAbeles(0.5Q[iBeg:iFin],resQ[iBeg:iFin],depth,rho,irho,sigma[1:])
3222	else: #dQ/Q in data
3223	AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
3224	Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
3225
3226	def abeles(kz, depth, rho, irho=0, sigma=0):
3227	"""
3228	Optical matrix form of the reflectivity calculation.
3229	O.S. Heavens, Optical Properties of Thin Solid Films
3230
3231	Reflectometry as a function of kz for a set of slabs.
3232
3233	:param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
3234	This is :math:`\\tfrac12 Q_z`.
3235	:param depth: float[m] (Ang).
3236	thickness of each layer. The thickness of the incident medium
3237	and substrate are ignored.
3238	:param rho: float[n,k] (1e-6/Ang^2)
3239	Real scattering length density for each layer for each kz
3240	:param irho: float[n,k] (1e-6/Ang^2)
3241	Imaginary scattering length density for each layer for each kz
3242	Note: absorption cross section mu = 2 irho/lambda for neutrons
3243	:param sigma: float[m-1] (Ang)
3244	interfacial roughness. This is the roughness between a layer
3245	and the previous layer. The sigma array should have m-1 entries.
3246
3247	Slabs are ordered with the surface SLD at index 0 and substrate at
3248	index -1, or reversed if kz < 0.
3249	"""
3250	def calc(kz, depth, rho, irho, sigma):
3251	if len(kz) == 0: return kz
3252
3253	# Complex index of refraction is relative to the incident medium.
3254	# We can get the same effect using kz_rel^2 = kz^2 + 4pirho_o
3255	# in place of kz^2, and ignoring rho_o
3256	kz_sq = kz*2 + 4e-6np.pi*rho[:,0]
3257	k = kz
3258
3259	# According to Heavens, the initial matrix should be [ 1 F; F 1],
3260	# which we do by setting B=I and M0 to [1 F; F 1]. An extra matrix
3261	# multiply versus some coding convenience.
3262	B11 = 1
3263	B22 = 1
3264	B21 = 0
3265	B12 = 0
3266	for i in range(0, len(depth)-1):
3267	k_next = np.sqrt(kz_sq - 4e-6np.pi(rho[:,i+1] + 1j*irho[:,i+1]))
3268	F = (k - k_next) / (k + k_next)
3269	F = np.exp(-2kk_nextsigma[i]**2)
3270	#print "==== layer",i
3271	#print "kz:", kz
3272	#print "k:", k
3273	#print "k_next:",k_next
3274	#print "F:",F
3275	#print "rho:",rho[:,i+1]
3276	#print "irho:",irho[:,i+1]
3277	#print "d:",depth[i],"sigma:",sigma[i]
3278	M11 = np.exp(1jkdepth[i]) if i>0 else 1
3279	M22 = np.exp(-1jkdepth[i]) if i>0 else 1
3280	M21 = F*M11
3281	M12 = F*M22
3282	C1 = B11M11 + B21M12
3283	C2 = B11M21 + B21M22
3284	B11 = C1
3285	B21 = C2
3286	C1 = B12M11 + B22M12
3287	C2 = B12M21 + B22M22
3288	B12 = C1
3289	B22 = C2
3290	k = k_next
3291
3292	r = B12/B11
3293	return np.absolute(r)**2
3294
3295	if np.isscalar(kz): kz = np.asarray([kz], 'd')
3296
3297	m = len(depth)
3298
3299	# Make everything into arrays
3300	depth = np.asarray(depth,'d')
3301	rho = np.asarray(rho,'d')
3302	irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
3303	sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
3304
3305	# Repeat rho,irho columns as needed
3306	if len(rho.shape) == 1:
3307	rho = rho[None,:]
3308	irho = irho[None,:]
3309
3310	return calc(kz, depth, rho, irho, sigma)
3311
3312	def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
3313	y = abeles(kz, depth, rho, irho, sigma)
3314	s = dq/2.
3315	y += 0.1354(abeles(kz+2s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
3316	y += 0.24935(abeles(kz-5s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma))
3317	y += 0.4111(abeles(kz-4s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma))
3318	y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
3319	y += 0.80074(abeles(kz-2s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
3320	y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
3321	y *= 0.137023
3322	return y
3323
3324	def makeRefdFFT(Limits,Profile):
3325	G2fil.G2Print ('make fft')
3326	Q,Io = Profile[:2]
3327	Qmin = Limits[1][0]
3328	Qmax = Limits[1][1]
3329	iBeg = np.searchsorted(Q,Qmin)
3330	iFin = np.searchsorted(Q,Qmax)+1 #include last point
3331	Qf = np.linspace(0.,Q[iFin-1],5000)
3332	QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
3333	If = QI(Qf)Qf*4
3334	R = np.linspace(0.,5000.,5000)
3335	F = fft.rfft(If)
3336	return R,F
3337
3338
3339	################################################################################
3340	#### Stacking fault simulation codes
3341	################################################################################
3342
3343	def GetStackParms(Layers):
3344
3345	Parms = []
3346	#cell parms
3347	if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
3348	Parms.append('cellA')
3349	Parms.append('cellC')
3350	else:
3351	Parms.append('cellA')
3352	Parms.append('cellB')
3353	Parms.append('cellC')
3354	if Layers['Laue'] != 'mmm':
3355	Parms.append('cellG')
3356	#Transition parms
3357	for iY in range(len(Layers['Layers'])):
3358	for iX in range(len(Layers['Layers'])):
3359	Parms.append('TransP;%d;%d'%(iY,iX))
3360	Parms.append('TransX;%d;%d'%(iY,iX))
3361	Parms.append('TransY;%d;%d'%(iY,iX))
3362	Parms.append('TransZ;%d;%d'%(iY,iX))
3363	return Parms
3364
3365	def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
3366	'''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
3367
3368	:param dict Layers: dict with following items
3369
3370	::
3371
3372	{'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
3373	'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
3374	'Layers':[],'Stacking':[],'Transitions':[]}
3375
3376	:param str ctrls: controls string to be written on DIFFaX controls.dif file
3377	:param float scale: scale factor
3378	:param dict background: background parameters
3379	:param list limits: min/max 2-theta to be calculated
3380	:param dict inst: instrument parameters dictionary
3381	:param list profile: powder pattern data
3382
3383	Note that parameters all updated in place
3384	'''
3385	import atmdata
3386	path = sys.path
3387	for name in path:
3388	if 'bin' in name:
3389	DIFFaX = name+'/DIFFaX.exe'
3390	G2fil.G2Print (' Execute '+DIFFaX)
3391	break
3392	# make form factor file that DIFFaX wants - atom types are GSASII style
3393	sf = open('data.sfc','w')
3394	sf.write('GSASII special form factor file for DIFFaX\n\n')
3395	atTypes = list(Layers['AtInfo'].keys())
3396	if 'H' not in atTypes:
3397	atTypes.insert(0,'H')
3398	for atType in atTypes:
3399	if atType == 'H':
3400	blen = -.3741
3401	else:
3402	blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
3403	Adat = atmdata.XrayFF[atType]
3404	text = '%4s'%(atType.ljust(4))
3405	for i in range(4):
3406	text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
3407	text += '%11.6f%11.6f'%(Adat['fc'],blen)
3408	text += '%3d\n'%(Adat['Z'])
3409	sf.write(text)
3410	sf.close()
3411	#make DIFFaX control.dif file - future use GUI to set some of these flags
3412	cf = open('control.dif','w')
3413	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
3414	x0 = profile[0]
3415	iBeg = np.searchsorted(x0,limits[0])
3416	iFin = np.searchsorted(x0,limits[1])+1
3417	if iFin-iBeg > 20000:
3418	iFin = iBeg+20000
3419	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3420	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3421	cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
3422	else:
3423	cf.write('GSASII-DIFFaX.dat\n'+ctrls)
3424	inst = {'Type':['XSC','XSC',]}
3425	cf.close()
3426	#make DIFFaX data file
3427	df = open('GSASII-DIFFaX.dat','w')
3428	df.write('INSTRUMENTAL\n')
3429	if 'X' in inst['Type'][0]:
3430	df.write('X-RAY\n')
3431	elif 'N' in inst['Type'][0]:
3432	df.write('NEUTRON\n')
3433	if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n':
3434	df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
3435	U = ateln2*inst['U'][1]/10000.
3436	V = ateln2*inst['V'][1]/10000.
3437	W = ateln2*inst['W'][1]/10000.
3438	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
3439	HW = np.sqrt(np.mean(HWHM))
3440	# df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
3441	if 'Mean' in Layers['selInst']:
3442	df.write('GAUSSIAN %.6f TRIM\n'%(HW)) #fast option - might not really matter
3443	elif 'Gaussian' in Layers['selInst']:
3444	df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W)) #slow - make a GUI option?
3445	else:
3446	df.write('None\n')
3447	else:
3448	df.write('0.10\nNone\n')
3449	df.write('STRUCTURAL\n')
3450	a,b,c = Layers['Cell'][1:4]
3451	gam = Layers['Cell'][6]
3452	df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
3453	laue = Layers['Laue']
3454	if laue == '2/m(ab)':
3455	laue = '2/m(1)'
3456	elif laue == '2/m(c)':
3457	laue = '2/m(2)'
3458	if 'unknown' in Layers['Laue']:
3459	df.write('%s %.3f\n'%(laue,Layers['Toler']))
3460	else:
3461	df.write('%s\n'%(laue))
3462	df.write('%d\n'%(len(Layers['Layers'])))
3463	if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
3464	df.write('%.1f %.1f\n'%(Layers['Width'][0][0]10000.,Layers['Width'][0][0]10000.)) #mum to A
3465	layerNames = []
3466	for layer in Layers['Layers']:
3467	layerNames.append(layer['Name'])
3468	for il,layer in enumerate(Layers['Layers']):
3469	if layer['SameAs']:
3470	df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
3471	continue
3472	df.write('LAYER %d\n'%(il+1))
3473	if '-1' in layer['Symm']:
3474	df.write('CENTROSYMMETRIC\n')
3475	else:
3476	df.write('NONE\n')
3477	for ia,atom in enumerate(layer['Atoms']):
3478	[name,atype,x,y,z,frac,Uiso] = atom
3479	if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
3480	frac /= 2.
3481	df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
3482	df.write('STACKING\n')
3483	df.write('%s\n'%(Layers['Stacking'][0]))
3484	if 'recursive' in Layers['Stacking'][0]:
3485	df.write('%s\n'%Layers['Stacking'][1])
3486	else:
3487	if 'list' in Layers['Stacking'][1]:
3488	Slen = len(Layers['Stacking'][2])
3489	iB = 0
3490	iF = 0
3491	while True:
3492	iF += 68
3493	if iF >= Slen:
3494	break
3495	iF = min(iF,Slen)
3496	df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
3497	iB = iF
3498	else:
3499	df.write('%s\n'%Layers['Stacking'][1])
3500	df.write('TRANSITIONS\n')
3501	for iY in range(len(Layers['Layers'])):
3502	sumPx = 0.
3503	for iX in range(len(Layers['Layers'])):
3504	p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
3505	p = round(p,3)
3506	df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
3507	sumPx += p
3508	if sumPx != 1.0: #this has to be picky since DIFFaX is.
3509	G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
3510	df.close()
3511	os.remove('data.sfc')
3512	os.remove('control.dif')
3513	os.remove('GSASII-DIFFaX.dat')
3514	return
3515	df.close()
3516	time0 = time.time()
3517	try:
3518	subp.call(DIFFaX)
3519	except OSError:
3520	G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
3521	G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
3522	if os.path.exists('GSASII-DIFFaX.spc'):
3523	Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
3524	iFin = iBeg+Xpat.shape[1]
3525	bakType,backDict,backVary = SetBackgroundParms(background)
3526	backDict['Lam1'] = G2mth.getWave(inst)
3527	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
3528	profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
3529	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
3530	rv = st.poisson(profile[3][iBeg:iFin])
3531	profile[1][iBeg:iFin] = rv.rvs()
3532	Z = np.ones_like(profile[3][iBeg:iFin])
3533	Z[1::2] *= -1
3534	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3535	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3536	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3537	#cleanup files..
3538	os.remove('GSASII-DIFFaX.spc')
3539	elif os.path.exists('GSASII-DIFFaX.sadp'):
3540	Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
3541	Sadp = np.reshape(Sadp,(256,-1))
3542	Layers['Sadp']['Img'] = Sadp
3543	os.remove('GSASII-DIFFaX.sadp')
3544	os.remove('data.sfc')
3545	os.remove('control.dif')
3546	os.remove('GSASII-DIFFaX.dat')
3547
3548	def SetPWDRscan(inst,limits,profile):
3549
3550	wave = G2mth.getMeanWave(inst)
3551	x0 = profile[0]
3552	iBeg = np.searchsorted(x0,limits[0])
3553	iFin = np.searchsorted(x0,limits[1])
3554	if iFin-iBeg > 20000:
3555	iFin = iBeg+20000
3556	Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
3557	pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
3558	return iFin-iBeg
3559
3560	def SetStackingSF(Layers,debug):
3561	# Load scattering factors into DIFFaX arrays
3562	import atmdata
3563	atTypes = Layers['AtInfo'].keys()
3564	aTypes = []
3565	for atype in atTypes:
3566	aTypes.append('%4s'%(atype.ljust(4)))
3567	SFdat = []
3568	for atType in atTypes:
3569	Adat = atmdata.XrayFF[atType]
3570	SF = np.zeros(9)
3571	SF[:8:2] = Adat['fa']
3572	SF[1:8:2] = Adat['fb']
3573	SF[8] = Adat['fc']
3574	SFdat.append(SF)
3575	SFdat = np.array(SFdat)
3576	pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
3577
3578	def SetStackingClay(Layers,Type):
3579	# Controls
3580	rand.seed()
3581	ranSeed = rand.randint(1,2**16-1)
3582	try:
3583	laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
3584	'6/m','6/mmm'].index(Layers['Laue'])+1
3585	except ValueError: #for 'unknown'
3586	laueId = -1
3587	if 'SADP' in Type:
3588	planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
3589	lmax = int(Layers['Sadp']['Lmax'])
3590	else:
3591	planeId = 0
3592	lmax = 0
3593	# Sequences
3594	StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
3595	try:
3596	StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
3597	except ValueError:
3598	StkParm = -1
3599	if StkParm == 2: #list
3600	StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
3601	Nstk = len(StkSeq)
3602	else:
3603	Nstk = 1
3604	StkSeq = [0,]
3605	if StkParm == -1:
3606	StkParm = int(Layers['Stacking'][1])
3607	Wdth = Layers['Width'][0]
3608	mult = 1
3609	controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
3610	LaueSym = Layers['Laue'].ljust(12)
3611	pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
3612	return laueId,controls
3613
3614	def SetCellAtoms(Layers):
3615	Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
3616	# atoms in layers
3617	atTypes = list(Layers['AtInfo'].keys())
3618	AtomXOU = []
3619	AtomTp = []
3620	LayerSymm = []
3621	LayerNum = []
3622	layerNames = []
3623	Natm = 0
3624	Nuniq = 0
3625	for layer in Layers['Layers']:
3626	layerNames.append(layer['Name'])
3627	for il,layer in enumerate(Layers['Layers']):
3628	if layer['SameAs']:
3629	LayerNum.append(layerNames.index(layer['SameAs'])+1)
3630	continue
3631	else:
3632	LayerNum.append(il+1)
3633	Nuniq += 1
3634	if '-1' in layer['Symm']:
3635	LayerSymm.append(1)
3636	else:
3637	LayerSymm.append(0)
3638	for ia,atom in enumerate(layer['Atoms']):
3639	[name,atype,x,y,z,frac,Uiso] = atom
3640	Natm += 1
3641	AtomTp.append('%4s'%(atype.ljust(4)))
3642	Ta = atTypes.index(atype)+1
3643	AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
3644	AtomXOU = np.array(AtomXOU)
3645	Nlayers = len(layerNames)
3646	pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
3647	return Nlayers
3648
3649	def SetStackingTrans(Layers,Nlayers):
3650	# Transitions
3651	TransX = []
3652	TransP = []
3653	for Ytrans in Layers['Transitions']:
3654	TransP.append([trans[0] for trans in Ytrans]) #get just the numbers
3655	TransX.append([trans[1:4] for trans in Ytrans]) #get just the numbers
3656	TransP = np.array(TransP,dtype='float').T
3657	TransX = np.array(TransX,dtype='float')
3658	# GSASIIpath.IPyBreak()
3659	pyx.pygettrans(Nlayers,TransP,TransX)
3660
3661	def CalcStackingPWDR(Layers,scale,background,limits,inst,profile,debug):
3662	# Scattering factors
3663	SetStackingSF(Layers,debug)
3664	# Controls & sequences
3665	laueId,controls = SetStackingClay(Layers,'PWDR')
3666	# cell & atoms
3667	Nlayers = SetCellAtoms(Layers)
3668	Volume = Layers['Cell'][7]
3669	# Transitions
3670	SetStackingTrans(Layers,Nlayers)
3671	# PWDR scan
3672	Nsteps = SetPWDRscan(inst,limits,profile)
3673	# result as Spec
3674	x0 = profile[0]
3675	profile[3] = np.zeros(len(profile[0]))
3676	profile[4] = np.zeros(len(profile[0]))
3677	profile[5] = np.zeros(len(profile[0]))
3678	iBeg = np.searchsorted(x0,limits[0])
3679	iFin = np.searchsorted(x0,limits[1])+1
3680	if iFin-iBeg > 20000:
3681	iFin = iBeg+20000
3682	Nspec = 20001
3683	spec = np.zeros(Nspec,dtype='double')
3684	time0 = time.time()
3685	pyx.pygetspc(controls,Nspec,spec)
3686	G2fil.G2Print (' GETSPC time = %.2fs'%(time.time()-time0))
3687	time0 = time.time()
3688	U = ateln2*inst['U'][1]/10000.
3689	V = ateln2*inst['V'][1]/10000.
3690	W = ateln2*inst['W'][1]/10000.
3691	HWHM = Unptand(x0[iBeg:iFin]/2.)2+Vnptand(x0[iBeg:iFin]/2.)+W
3692	HW = np.sqrt(np.mean(HWHM))
3693	BrdSpec = np.zeros(Nsteps)
3694	if 'Mean' in Layers['selInst']:
3695	pyx.pyprofile(U,V,W,HW,1,Nsteps,BrdSpec)
3696	elif 'Gaussian' in Layers['selInst']:
3697	pyx.pyprofile(U,V,W,HW,4,Nsteps,BrdSpec)
3698	else:
3699	BrdSpec = spec[:Nsteps]
3700	BrdSpec /= Volume
3701	iFin = iBeg+Nsteps
3702	bakType,backDict,backVary = SetBackgroundParms(background)
3703	backDict['Lam1'] = G2mth.getWave(inst)
3704	profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]
3705	profile[3][iBeg:iFin] = BrdSpec*scale+profile[4][iBeg:iFin]
3706	if not np.any(profile[1]): #fill dummy data x,y,w,yc,yb,yd
3707	try:
3708	rv = st.poisson(profile[3][iBeg:iFin])
3709	profile[1][iBeg:iFin] = rv.rvs()
3710	except ValueError:
3711	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]
3712	Z = np.ones_like(profile[3][iBeg:iFin])
3713	Z[1::2] *= -1
3714	profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
3715	profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
3716	profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
3717	G2fil.G2Print (' Broadening time = %.2fs'%(time.time()-time0))
3718
3719	def CalcStackingSADP(Layers,debug):
3720
3721	# Scattering factors
3722	SetStackingSF(Layers,debug)
3723	# Controls & sequences
3724	laueId,controls = SetStackingClay(Layers,'SADP')
3725	# cell & atoms
3726	Nlayers = SetCellAtoms(Layers)
3727	# Transitions
3728	SetStackingTrans(Layers,Nlayers)
3729	# result as Sadp
3730	Nspec = 20001
3731	spec = np.zeros(Nspec,dtype='double')
3732	time0 = time.time()
3733	hkLim,Incr,Nblk = pyx.pygetsadp(controls,Nspec,spec)
3734	Sapd = np.zeros((256,256))
3735	iB = 0
3736	for i in range(hkLim):
3737	iF = iB+Nblk
3738	p1 = 127+int(i*Incr)
3739	p2 = 128-int(i*Incr)
3740	if Nblk == 128:
3741	if i:
3742	Sapd[128:,p1] = spec[iB:iF]
3743	Sapd[:128,p1] = spec[iF:iB:-1]
3744	Sapd[128:,p2] = spec[iB:iF]
3745	Sapd[:128,p2] = spec[iF:iB:-1]
3746	else:
3747	if i:
3748	Sapd[:,p1] = spec[iB:iF]
3749	Sapd[:,p2] = spec[iB:iF]
3750	iB += Nblk
3751	Layers['Sadp']['Img'] = Sapd
3752	G2fil.G2Print (' GETSAD time = %.2fs'%(time.time()-time0))
3753
3754	###############################################################################
3755	#### Maximum Entropy Method - Dysnomia
3756	###############################################################################
3757
3758	def makePRFfile(data,MEMtype):
3759	''' makes Dysnomia .prf control file from Dysnomia GUI controls
3760
3761	:param dict data: GSAS-II phase data
3762	:param int MEMtype: 1 for neutron data with negative scattering lengths
3763	0 otherwise
3764	:returns str: name of Dysnomia control file
3765	'''
3766
3767	generalData = data['General']
3768	pName = generalData['Name'].replace(' ','_')
3769	DysData = data['Dysnomia']
3770	prfName = pName+'.prf'
3771	prf = open(prfName,'w')
3772	prf.write('$PREFERENCES\n')
3773	prf.write(pName+'.mem\n') #or .fos?
3774	prf.write(pName+'.out\n')
3775	prf.write(pName+'.pgrid\n')
3776	prf.write(pName+'.fba\n')
3777	prf.write(pName+'_eps.raw\n')
3778	prf.write('%d\n'%MEMtype)
3779	if DysData['DenStart'] == 'uniform':
3780	prf.write('0\n')
3781	else:
3782	prf.write('1\n')
3783	if DysData['Optimize'] == 'ZSPA':
3784	prf.write('0\n')
3785	else:
3786	prf.write('1\n')
3787	prf.write('1\n')
3788	if DysData['Lagrange'][0] == 'user':
3789	prf.write('0\n')
3790	else:
3791	prf.write('1\n')
3792	prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
3793	prf.write('%.3f\n'%DysData['Lagrange'][2])
3794	prf.write('%.2f\n'%DysData['E_factor'])
3795	prf.write('1\n')
3796	prf.write('0\n')
3797	prf.write('%d\n'%DysData['Ncyc'])
3798	prf.write('1\n')
3799	prf.write('1 0 0 0 0 0 0 0\n')
3800	if DysData['prior'] == 'uniform':
3801	prf.write('0\n')
3802	else:
3803	prf.write('1\n')
3804	prf.close()
3805	return prfName
3806
3807	def makeMEMfile(data,reflData,MEMtype,DYSNOMIA):
3808	''' make Dysnomia .mem file of reflection data, etc.
3809
3810	:param dict data: GSAS-II phase data
3811	:param list reflData: GSAS-II reflection data
3812	:param int MEMtype: 1 for neutron data with negative scattering lengths
3813	0 otherwise
3814	:param str DYSNOMIA: path to dysnomia.exe
3815	'''
3816
3817	DysData = data['Dysnomia']
3818	generalData = data['General']
3819	cell = generalData['Cell'][1:7]
3820	A = G2lat.cell2A(cell)
3821	SGData = generalData['SGData']
3822	pName = generalData['Name'].replace(' ','_')
3823	memName = pName+'.mem'
3824	Map = generalData['Map']
3825	Type = Map['Type']
3826	UseList = Map['RefList']
3827	mem = open(memName,'w')
3828	mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
3829	a,b,c,alp,bet,gam = cell
3830	mem.write('%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
3831	mem.write(' 0.0000000 0.0000000 -1 0 0 0 P\n') #dummy PO stuff
3832	SGSym = generalData['SGData']['SpGrp']
3833	try:
3834	SGId = G2spc.spgbyNum.index(SGSym)
3835	except ValueError:
3836	return False
3837	org = 1
3838	if SGSym in G2spc.spg2origins:
3839	org = 2
3840	mapsize = Map['rho'].shape
3841	sumZ = 0.
3842	sumpos = 0.
3843	sumneg = 0.
3844	mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
3845	for atm in generalData['NoAtoms']:
3846	Nat = generalData['NoAtoms'][atm]
3847	AtInfo = G2elem.GetAtomInfo(atm)
3848	sumZ += Nat*AtInfo['Z']
3849	isotope = generalData['Isotope'][atm]
3850	blen = generalData['Isotopes'][atm][isotope]['SL'][0]
3851	if blen < 0.:
3852	sumneg += blen*Nat
3853	else:
3854	sumpos += blen*Nat
3855	if 'X' in Type:
3856	mem.write('%10.2f 0.001\n'%sumZ)
3857	elif 'N' in Type and MEMtype:
3858	mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
3859	else:
3860	mem.write('%10.3f 0.001\n'%sumpos)
3861
3862	dmin = DysData['MEMdmin']
3863	TOFlam = 2.0dminnpsind(80.0)
3864	refSet = G2lat.GenHLaue(dmin,SGData,A) #list of h,k,l,d
3865	refDict = {'%d %d %d'%(ref[0],ref[1],ref[2]):ref for ref in refSet}
3866
3867	refs = []
3868	prevpos = 0.
3869	for ref in reflData:
3870	if ref[3] < 0:
3871	continue
3872	if 'T' in Type:
3873	h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
3874	s = np.sqrt(max(sig,0.0001)) #var -> sig in deg
3875	FWHM = getgamFW(gam,s)
3876	if dsp < dmin:
3877	continue
3878	theta = npasind(TOFlam/(2.*dsp))
3879	FWHM *= nptand(theta)/pos
3880	pos = 2.*theta
3881	else:
3882	h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
3883	g = gam/100. #centideg -> deg
3884	s = np.sqrt(max(sig,0.0001))/100. #var -> sig in deg
3885	FWHM = getgamFW(g,s)
3886	delt = pos-prevpos
3887	refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
3888	prevpos = pos
3889
3890	ovlp = DysData['overlap']
3891	refs1 = []
3892	refs2 = []
3893	nref2 = 0
3894	iref = 0
3895	Nref = len(refs)
3896	start = False
3897	while iref < Nref-1:
3898	if refs[iref+1][-1] < ovlp*refs[iref][5]:
3899	if refs[iref][-1] > ovlp*refs[iref][5]:
3900	refs2.append([])
3901	start = True
3902	if nref2 == len(refs2):
3903	refs2.append([])
3904	refs2[nref2].append(refs[iref])
3905	else:
3906	if start:
3907	refs2[nref2].append(refs[iref])
3908	start = False
3909	nref2 += 1
3910	else:
3911	refs1.append(refs[iref])
3912	iref += 1
3913	if start:
3914	refs2[nref2].append(refs[iref])
3915	else:
3916	refs1.append(refs[iref])
3917
3918	mem.write('%5d\n'%len(refs1))
3919	for ref in refs1:
3920	h,k,l = ref[:3]
3921	hkl = '%d %d %d'%(h,k,l)
3922	if hkl in refDict:
3923	del refDict[hkl]
3924	Fobs = np.sqrt(ref[6])
3925	mem.write('%5d%5d%5d%10.3f%10.3f%10.3f\n'%(h,k,l,Fobsnpcosd(ref[7]),Fobsnpsind(ref[7]),max(0.01*Fobs,0.1)))
3926	while True and nref2:
3927	if not len(refs2[-1]):
3928	del refs2[-1]
3929	else:
3930	break
3931	mem.write('%5d\n'%len(refs2))
3932	for iref2,ref2 in enumerate(refs2):
3933	mem.write('#%5d\n'%iref2)
3934	mem.write('%5d\n'%len(ref2))
3935	Gsum = 0.
3936	Msum = 0
3937	for ref in ref2:
3938	Gsum += ref[6]*ref[3]
3939	Msum += ref[3]
3940	G = np.sqrt(Gsum/Msum)
3941	h,k,l = ref2[0][:3]
3942	hkl = '%d %d %d'%(h,k,l)
3943	if hkl in refDict:
3944	del refDict[hkl]
3945	mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
3946	for ref in ref2[1:]:
3947	h,k,l,m = ref[:4]
3948	mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
3949	hkl = '%d %d %d'%(h,k,l)
3950	if hkl in refDict:
3951	del refDict[hkl]
3952	if len(refDict):
3953	mem.write('%d\n'%len(refDict))
3954	for hkl in list(refDict.keys()):
3955	h,k,l = refDict[hkl][:3]
3956	mem.write('%5d%5d%5d\n'%(h,k,l))
3957	else:
3958	mem.write('0\n')
3959	mem.close()
3960	return True
3961
3962	def MEMupdateReflData(prfName,data,reflData):
3963	''' Update reflection data with new Fosq, phase result from Dysnomia
3964
3965	:param str prfName: phase.mem file name
3966	:param list reflData: GSAS-II reflection data
3967	'''
3968
3969	generalData = data['General']
3970	Map = generalData['Map']
3971	Type = Map['Type']
3972	cell = generalData['Cell'][1:7]
3973	A = G2lat.cell2A(cell)
3974	reflDict = {}
3975	newRefs = []
3976	for iref,ref in enumerate(reflData):
3977	if ref[3] > 0:
3978	newRefs.append(ref)
3979	reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
3980	fbaName = os.path.splitext(prfName)[0]+'.fba'
3981	if os.path.isfile(fbaName):
3982	fba = open(fbaName,'r')
3983	else:
3984	return False
3985	fba.readline()
3986	Nref = int(fba.readline()[:-1])
3987	fbalines = fba.readlines()
3988	for line in fbalines[:Nref]:
3989	info = line.split()
3990	h = int(info[0])
3991	k = int(info[1])
3992	l = int(info[2])
3993	FoR = float(info[3])
3994	FoI = float(info[4])
3995	Fosq = FoR2+FoI2
3996	phase = npatan2d(FoI,FoR)
3997	try:
3998	refId = reflDict[hash('%5d%5d%5d'%(h,k,l))]
3999	except KeyError: #added reflections at end skipped
4000	d = float(1/np.sqrt(G2lat.calc_rDsq([h,k,l],A)))
4001	if 'T' in Type:
4002	newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0,1.0,1.0,1.0])
4003	else:
4004	newRefs.append([h,k,l,-1,d,0.,0.01,1.0,Fosq,Fosq,phase,1.0,1.0,1.0,1.0])
4005	continue
4006	newRefs[refId][8] = Fosq
4007	newRefs[refId][10] = phase
4008	newRefs = np.array(newRefs)
4009	return True,newRefs
4010
4011	#### testing data
4012	NeedTestData = True
4013	def TestData():
4014	'needs a doc string'
4015	# global NeedTestData
4016	global bakType
4017	bakType = 'chebyschev'
4018	global xdata
4019	xdata = np.linspace(4.0,40.0,36000)
4020	global parmDict0
4021	parmDict0 = {
4022	'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
4023	'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
4024	'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
4025	'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
4026	'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'Z':0.0,'SH/L':0.002,
4027	'Back0':5.384,'Back1':-0.015,'Back2':.004,
4028	}
4029	global parmDict1
4030	parmDict1 = {
4031	'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
4032	'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
4033	'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
4034	'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
4035	'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
4036	'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
4037	'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'Z':0.0,'SH/L':0.002,
4038	'Back0':36.897,'Back1':-0.508,'Back2':.006,
4039	'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
4040	}
4041	global parmDict2
4042	parmDict2 = {
4043	'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
4044	'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'Z':0.0,'SH/L':0.02,
4045	'Back0':5.,'Back1':-0.02,'Back2':.004,
4046	# 'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
4047	}
4048	global varyList
4049	varyList = []
4050
4051	def test0():
4052	if NeedTestData: TestData()
4053	gplot = plotter.add('FCJ-Voigt, 11BM').gca()
4054	gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])
4055	gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
4056	fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
4057	fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])
4058	fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
4059
4060	def test1():
4061	if NeedTestData: TestData()
4062	time0 = time.time()
4063	for i in range(100):
4064	getPeakProfile(parmDict1,xdata,varyList,bakType)
4065	G2fil.G2Print ('100+6*Ka1-2 peaks=1200 peaks %.2f'%time.time()-time0)
4066
4067	def test2(name,delt):
4068	if NeedTestData: TestData()
4069	varyList = [name,]
4070	xdata = np.linspace(5.6,5.8,400)
4071	hplot = plotter.add('derivatives test for '+name).gca()
4072	hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
4073	y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
4074	parmDict2[name] += delt
4075	y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
4076	hplot.plot(xdata,(y1-y0)/delt,'r+')
4077
4078	def test3(name,delt):
4079	if NeedTestData: TestData()
4080	names = ['pos','sig','gam','shl']
4081	idx = names.index(name)
4082	myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
4083	xdata = np.linspace(5.6,5.8,800)
4084	dx = xdata[1]-xdata[0]
4085	hplot = plotter.add('derivatives test for '+name).gca()
4086	hplot.plot(xdata,100.dxgetdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
4087	y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
4088	myDict[name] += delt
4089	y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
4090	hplot.plot(xdata,(y1-y0)/delt,'r+')
4091
4092	if __name__ == '__main__':
4093	import GSASIItestplot as plot
4094	global plotter
4095	plotter = plot.PlotNotebook()
4096	# test0()
4097	# for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)']:
4098	for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
4099	['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['Z',0.0001],['SH/L',0.00005]]:
4100	test2(name,shft)
4101	for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
4102	test3(name,shft)
4103	G2fil.G2Print ("OK")
4104	plotter.StartEventLoop()

Note: See TracBrowser for help on using the repository browser.

Download in other formats: