Context Navigation

source: trunk/GSASIIpwd.py @ 5281

Last change on this file since 5281 was 5281, checked in by toby, 19 months ago
start putting paths into PDFfit script
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 220.2 KB

Line
1	#/usr/bin/env python
2	# -- coding: utf-8 --
3	'''
4	GSASIIpwd: Powder calculations module
5	==============================================
6
7	This version hacked to provide Laue Fringe fitting.
8
9	'''
10	########### SVN repository information ###################
11	# $Date: 2022-05-18 19:16:10 +0000 (Wed, 18 May 2022) $
12	# $Author: toby $
13	# $Revision: 5281 $
14	# $URL: trunk/GSASIIpwd.py $
15	# $Id: GSASIIpwd.py 5281 2022-05-18 19:16:10Z toby $
16	########### SVN repository information ###################
17	from __future__ import division, print_function
18	import sys
19	import math
20	import time
21	import os
22	import os.path
23	import subprocess as subp
24	import datetime as dt
25	import copy
26
27	import numpy as np
28	import numpy.linalg as nl
29	import numpy.ma as ma
30	import random as rand
31	import numpy.fft as fft
32	import scipy.interpolate as si
33	import scipy.stats as st
34	import scipy.optimize as so
35	import scipy.special as sp
36	import scipy.signal as signal
37
38	import GSASIIpath
39	filversion = "$Revision: 5281 $"
40	GSASIIpath.SetVersionNumber("$Revision: 5281 $")
41	import GSASIIlattice as G2lat
42	import GSASIIspc as G2spc
43	import GSASIIElem as G2elem
44	import GSASIImath as G2mth
45	try:
46	import pypowder as pyd
47	except ImportError:
48	print ('pypowder is not available - profile calcs. not allowed')
49	try:
50	import pydiffax as pyx
51	except ImportError:
52	print ('pydiffax is not available for this platform')
53	import GSASIIfiles as G2fil
54
55	# trig functions in degrees
56	tand = lambda x: math.tan(x*math.pi/180.)
57	atand = lambda x: 180.*math.atan(x)/math.pi
58	atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
59	cosd = lambda x: math.cos(x*math.pi/180.)
60	acosd = lambda x: 180.*math.acos(x)/math.pi
61	rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
62	#numpy versions
63	npsind = lambda x: np.sin(x*np.pi/180.)
64	npasind = lambda x: 180.*np.arcsin(x)/math.pi
65	npcosd = lambda x: np.cos(x*math.pi/180.)
66	npacosd = lambda x: 180.*np.arccos(x)/math.pi
67	nptand = lambda x: np.tan(x*math.pi/180.)
68	npatand = lambda x: 180.*np.arctan(x)/np.pi
69	npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
70	npT2stl = lambda tth, wave: 2.0npsind(tth/2.0)/wave #=d
71	npT2q = lambda tth,wave: 2.0np.pinpT2stl(tth,wave) #=2pid
72	npq2T = lambda Q,wave: 2.0npasind(0.25Q*wave/np.pi)
73	ateln2 = 8.0*math.log(2.0)
74	sateln2 = np.sqrt(ateln2)
75	nxs = np.newaxis
76	is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK)
77
78	#### Powder utilities ################################################################################
79	def PhaseWtSum(G2frame,histo):
80	'''
81	Calculate sum of phase mass*phase fraction for PWDR data (exclude magnetic phases)
82
83	:param G2frame: GSASII main frame structure
84	:param str histo: histogram name
85	:returns: sum(scale*mass) for phases in histo
86	'''
87	Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
88	wtSum = 0.0
89	for phase in Phases:
90	if Phases[phase]['General']['Type'] != 'magnetic':
91	if histo in Phases[phase]['Histograms']:
92	if not Phases[phase]['Histograms'][histo]['Use']: continue
93	mass = Phases[phase]['General']['Mass']
94	phFr = Phases[phase]['Histograms'][histo]['Scale'][0]
95	wtSum += mass*phFr
96	return wtSum
97
98	#### GSASII pwdr & pdf calculation routines ################################################################################
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 - can be numpy array
229	:param Tth: 2-theta scattering angle - can be numpy array
230	which (if either) of these is "right"?
231	:return: (pola, dpdPola) - both 2-d arrays
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)**2
238	sazm = npsind(Azm)**2
239	pola = ((1.0-Pola)cazm+Polasazm)npcosd(Tth)2+(1.0-Pola)sazm+Pola*cazm
240	dpdPola = -npsind(Tth)*2(sazm-cazm)
241	return pola,dpdPola
242
243	def Oblique(ObCoeff,Tth):
244	'currently assumes detector is normal to beam'
245	if ObCoeff:
246	K = (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
247	return K
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 #Compton wavelength in A
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 KleinNishina(wave,Q):
262	hmc = 0.024262734687 #Compton wavelength in A
263	TTh = npq2T(Q,wave)
264	P = 1./(1.+(1.-npcosd(TTh)*(hmc/wave)))
265	KN = (P*3-(Pnpsind(TTh))2+P)/(1.+npcosd(TTh)2)
266	return KN
267
268	def LorchWeight(Q):
269	'needs a doc string'
270	return np.sin(np.pi(Q[-1]-Q)/(2.0Q[-1]))
271
272	def GetAsfMean(ElList,Sthl2):
273	'''Calculate various scattering factor terms for PDF calcs
274
275	:param dict ElList: element dictionary contains scattering factor coefficients, etc.
276	:param np.array Sthl2: numpy array of sin theta/lambda squared values
277	:returns: mean(f^2), mean(f)^2, mean(compton)
278	'''
279	sumNoAtoms = 0.0
280	FF = np.zeros_like(Sthl2)
281	FF2 = np.zeros_like(Sthl2)
282	CF = np.zeros_like(Sthl2)
283	for El in ElList:
284	sumNoAtoms += ElList[El]['FormulaNo']
285	for El in ElList:
286	el = ElList[El]
287	ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])2+el['fpp']2
288	cf = G2elem.ComptonFac(el,Sthl2)
289	FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
290	FF2 += ff2*el['FormulaNo']/sumNoAtoms
291	CF += cf*el['FormulaNo']/sumNoAtoms
292	return FF2,FF**2,CF
293
294	def GetNumDensity(ElList,Vol):
295	'needs a doc string'
296	sumNoAtoms = 0.0
297	for El in ElList:
298	sumNoAtoms += ElList[El]['FormulaNo']
299	return sumNoAtoms/Vol
300
301	def CalcPDF(data,inst,limits,xydata):
302	'''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
303	dict xydata; results are placed in xydata.
304	Calculation parameters are found in dicts data and inst and list limits.
305	The return value is at present an empty list.
306	'''
307	auxPlot = []
308	if 'T' in inst['Type'][0]:
309	Ibeg = 0
310	Ifin = len(xydata['Sample'][1][0])
311	else:
312	Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0])
313	Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1
314	#subtract backgrounds - if any & use PWDR limits
315	IofQ = copy.deepcopy(xydata['Sample'])
316	IofQ[1] = np.array([I[Ibeg:Ifin] for I in IofQ[1]])
317	if data['Sample Bkg.']['Name']:
318	IofQ[1][1] += xydata['Sample Bkg.'][1][1][Ibeg:Ifin]*data['Sample Bkg.']['Mult']
319	if data['Container']['Name']:
320	xycontainer = xydata['Container'][1][1]*data['Container']['Mult']
321	if data['Container Bkg.']['Name']:
322	xycontainer += xydata['Container Bkg.'][1][1][Ibeg:Ifin]*data['Container Bkg.']['Mult']
323	IofQ[1][1] += xycontainer[Ibeg:Ifin]
324	data['IofQmin'] = IofQ[1][1][-1]
325	IofQ[1][1] -= data.get('Flat Bkg',0.)
326	#get element data & absorption coeff.
327	ElList = data['ElList']
328	Tth = IofQ[1][0] #2-theta or TOF!
329	if 'X' in inst['Type'][0]:
330	Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
331	#Apply angle dependent corrections
332	MuR = Abs*data['Diam']/20.0
333	IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
334	IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
335	if data['DetType'] == 'Area detector':
336	IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
337	elif 'T' in inst['Type'][0]: #neutron TOF normalized data - needs wavelength dependent absorption
338	wave = 2.G2lat.TOF2dsp(inst,IofQ[1][0])npsind(inst['2-theta'][1]/2.)
339	Els = ElList.keys()
340	Isotope = {El:'Nat. abund.' for El in Els}
341	GD = {'AtomTypes':ElList,'Isotope':Isotope}
342	BLtables = G2elem.GetBLtable(GD)
343	FP,FPP = G2elem.BlenResTOF(Els,BLtables,wave)
344	Abs = np.zeros(len(wave))
345	for iel,El in enumerate(Els):
346	BL = BLtables[El][1]
347	SA = BL['SA']wave/1.798197+4.0np.piFPP[iel]*2 #+BL['SL'][1]?
348	SA *= ElList[El]['FormulaNo']/data['Form Vol']
349	Abs += SA
350	MuR = Abs*data['Diam']/2.
351	IofQ[1][1] /= Absorb(data['Geometry'],MuR,inst['2-theta'][1]*np.ones(len(wave)))
352	# improves look of F(Q) but no impact on G(R)
353	# bBut,aBut = signal.butter(8,.5,"lowpass")
354	# IofQ[1][1] = signal.filtfilt(bBut,aBut,IofQ[1][1])
355	XY = IofQ[1]
356	#convert to Q
357	nQpoints = 5000
358	if 'C' in inst['Type'][0]:
359	wave = G2mth.getWave(inst)
360	minQ = npT2q(Tth[0],wave)
361	maxQ = npT2q(Tth[-1],wave)
362	Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
363	dq = Qpoints[1]-Qpoints[0]
364	XY[0] = npT2q(XY[0],wave)
365	Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0]) #interpolate I(Q)
366	elif 'T' in inst['Type'][0]:
367	difC = inst['difC'][1]
368	minQ = 2.np.pidifC/Tth[-1]
369	maxQ = 2.np.pidifC/Tth[0]
370	Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
371	dq = Qpoints[1]-Qpoints[0]
372	XY[0] = 2.np.pidifC/XY[0]
373	Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][-1]) #interpolate I(Q)
374	Qdata -= np.min(Qdata)*data['BackRatio']
375
376	qLimits = data['QScaleLim']
377	maxQ = np.searchsorted(Qpoints,min(Qpoints[-1],qLimits[1]))+1
378	minQ = np.searchsorted(Qpoints,min(qLimits[0],0.90*Qpoints[-1]))
379	qLimits = [Qpoints[minQ],Qpoints[maxQ-1]]
380	newdata = []
381	if len(IofQ) < 3:
382	xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],'']
383	else:
384	xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
385	for item in xydata['IofQ'][1]:
386	newdata.append(item[:maxQ])
387	xydata['IofQ'][1] = newdata
388
389	xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
390	if 'XC' in inst['Type'][0]:
391	FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0np.pi))*2) #these are <f^2>,<f>^2,Cf
392	else: #TOF
393	CF = np.zeros(len(xydata['SofQ'][1][0]))
394	FFSq = np.ones(len(xydata['SofQ'][1][0]))
395	SqFF = np.ones(len(xydata['SofQ'][1][0]))
396	Q = xydata['SofQ'][1][0]
397	# auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
398	if 'XC' in inst['Type'][0]:
399	# CF *= KleinNishina(wave,Q)
400	ruland = Ruland(data['Ruland'],wave,Q,CF)
401	# auxPlot.append([Q,ruland,'Ruland'])
402	CF *= ruland
403	# auxPlot.append([Q,CF,'CF-Corr'])
404	scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
405	xydata['SofQ'][1][1] *= scale
406	if 'XC' in inst['Type'][0]:
407	xydata['SofQ'][1][1] -= CF
408	xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
409	scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
410	xydata['SofQ'][1][1] *= scale
411	xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
412	xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
413	if data['Lorch']:
414	xydata['FofQ'][1][1] *= LorchWeight(Q)
415	xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
416	xydata['gofr'] = copy.deepcopy(xydata['FofQ'])
417	nR = len(xydata['GofR'][1][1])
418	Rmax = GSASIIpath.GetConfigValue('PDF_Rmax',100.)
419	mul = int(round(2.np.pinR/(Rmax*qLimits[1])))
420	# mul = int(round(2.np.pinR/(data.get('Rmax',100.)*qLimits[1])))
421	R = 2.np.pinp.linspace(0,nR,nR,endpoint=True)/(mul*qLimits[1])
422	xydata['GofR'][1][0] = R
423	xydata['gofr'][1][0] = R
424	GR = -dqnp.imag(fft.fft(xydata['FofQ'][1][1],mulnR)[:nR])*data.get('GR Scale',1.0)
425	xydata['GofR'][1][1] = GR
426	gr = GR/(np.pi*R)
427	xydata['gofr'][1][1] = gr
428	numbDen = 0.
429	if 'ElList' in data:
430	numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
431	if data.get('noRing',True):
432	Rmin = data['Rmin']
433	xydata['gofr'][1][1] = np.where(R<Rmin,-4.*numbDen,xydata['gofr'][1][1])
434	xydata['GofR'][1][1] = np.where(R<Rmin,-4.Rnp.pi*numbDen,xydata['GofR'][1][1])
435	return auxPlot
436
437	def PDFPeakFit(peaks,data):
438	rs2pi = 1./np.sqrt(2*np.pi)
439
440	def MakeParms(peaks):
441	varyList = []
442	parmDict = {'slope':peaks['Background'][1][1]}
443	if peaks['Background'][2]:
444	varyList.append('slope')
445	for i,peak in enumerate(peaks['Peaks']):
446	parmDict['PDFpos;'+str(i)] = peak[0]
447	parmDict['PDFmag;'+str(i)] = peak[1]
448	parmDict['PDFsig;'+str(i)] = peak[2]
449	if 'P' in peak[3]:
450	varyList.append('PDFpos;'+str(i))
451	if 'M' in peak[3]:
452	varyList.append('PDFmag;'+str(i))
453	if 'S' in peak[3]:
454	varyList.append('PDFsig;'+str(i))
455	return parmDict,varyList
456
457	def SetParms(peaks,parmDict,varyList):
458	if 'slope' in varyList:
459	peaks['Background'][1][1] = parmDict['slope']
460	for i,peak in enumerate(peaks['Peaks']):
461	if 'PDFpos;'+str(i) in varyList:
462	peak[0] = parmDict['PDFpos;'+str(i)]
463	if 'PDFmag;'+str(i) in varyList:
464	peak[1] = parmDict['PDFmag;'+str(i)]
465	if 'PDFsig;'+str(i) in varyList:
466	peak[2] = parmDict['PDFsig;'+str(i)]
467
468
469	def CalcPDFpeaks(parmdict,Xdata):
470	Z = parmDict['slope']*Xdata
471	ipeak = 0
472	while True:
473	try:
474	pos = parmdict['PDFpos;'+str(ipeak)]
475	mag = parmdict['PDFmag;'+str(ipeak)]
476	wid = parmdict['PDFsig;'+str(ipeak)]
477	wid2 = 2.wid*2
478	Z += magrs2pinp.exp(-(Xdata-pos)**2/wid2)/wid
479	ipeak += 1
480	except KeyError: #no more peaks to process
481	return Z
482
483	def errPDFProfile(values,xdata,ydata,parmdict,varylist):
484	parmdict.update(zip(varylist,values))
485	M = CalcPDFpeaks(parmdict,xdata)-ydata
486	return M
487
488	newpeaks = copy.copy(peaks)
489	iBeg = np.searchsorted(data[1][0],newpeaks['Limits'][0])
490	iFin = np.searchsorted(data[1][0],newpeaks['Limits'][1])+1
491	X = data[1][0][iBeg:iFin]
492	Y = data[1][1][iBeg:iFin]
493	parmDict,varyList = MakeParms(peaks)
494	if not len(varyList):
495	G2fil.G2Print (' Nothing varied')
496	return newpeaks,None,None,None,None,None
497
498	Rvals = {}
499	values = np.array(Dict2Values(parmDict, varyList))
500	result = so.leastsq(errPDFProfile,values,full_output=True,ftol=0.0001,
501	args=(X,Y,parmDict,varyList))
502	chisq = np.sum(result[2]['fvec']**2)
503	Values2Dict(parmDict, varyList, result[0])
504	SetParms(peaks,parmDict,varyList)
505	Rvals['Rwp'] = np.sqrt(chisq/np.sum(Y*2))100. #to %
506	chisq = np.sum(result[2]['fvec']**2)/(len(X)-len(values)) #reduced chi^2 = M/(Nobs-Nvar)
507	sigList = list(np.sqrt(chisq*np.diag(result[1])))
508	Z = CalcPDFpeaks(parmDict,X)
509	newpeaks['calc'] = [X,Z]
510	return newpeaks,result[0],varyList,sigList,parmDict,Rvals
511
512	def MakeRDF(RDFcontrols,background,inst,pwddata):
513	auxPlot = []
514	if 'C' in inst['Type'][0] or 'B' in inst['Type'][0]:
515	Tth = pwddata[0]
516	wave = G2mth.getWave(inst)
517	minQ = npT2q(Tth[0],wave)
518	maxQ = npT2q(Tth[-1],wave)
519	powQ = npT2q(Tth,wave)
520	elif 'T' in inst['Type'][0]:
521	TOF = pwddata[0]
522	difC = inst['difC'][1]
523	minQ = 2.np.pidifC/TOF[-1]
524	maxQ = 2.np.pidifC/TOF[0]
525	powQ = 2.np.pidifC/TOF
526	piDQ = np.pi/(maxQ-minQ)
527	Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
528	if RDFcontrols['UseObsCalc'] == 'obs-calc':
529	Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
530	elif RDFcontrols['UseObsCalc'] == 'obs-back':
531	Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
532	elif RDFcontrols['UseObsCalc'] == 'calc-back':
533	Qdata = si.griddata(powQ,pwddata[3]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
534	Qdata = np.sin((Qpoints-minQ)piDQ)/piDQ
535	Qdata = 0.5np.sqrt(Qpoints) #Qbin normalization
536	dq = Qpoints[1]-Qpoints[0]
537	nR = len(Qdata)
538	R = 0.5np.pinp.linspace(0,nR,nR)/(4.*maxQ)
539	iFin = np.searchsorted(R,RDFcontrols['maxR'])+1
540	bBut,aBut = signal.butter(4,0.01)
541	Qsmooth = signal.filtfilt(bBut,aBut,Qdata)
542	# auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']])
543	# auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']])
544	DofR = dqnp.imag(fft.fft(Qsmooth,16nR)[:nR])
545	auxPlot.append([R[:iFin],DofR[:iFin],'D(R) for '+RDFcontrols['UseObsCalc']])
546	return auxPlot
547
548	# PDF optimization =============================================================
549	def OptimizePDF(data,xydata,limits,inst,showFit=True,maxCycles=25):
550	import scipy.optimize as opt
551	numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
552	Min,Init,Done = SetupPDFEval(data,xydata,limits,inst,numbDen)
553	xstart = Init()
554	bakMul = data['Sample Bkg.']['Mult']
555	if showFit:
556	rms = Min(xstart)
557	G2fil.G2Print(' Optimizing corrections to improve G(r) at low r')
558	if data['Sample Bkg.'].get('Refine',False):
559	# data['Flat Bkg'] = 0.
560	G2fil.G2Print(' start: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f})'.format(
561	data['Ruland'],data['Sample Bkg.']['Mult'],rms))
562	else:
563	G2fil.G2Print(' start: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} (RMS:{:.4f})'.format(
564	data['Flat Bkg'],data['BackRatio'],data['Ruland'],rms))
565	if data['Sample Bkg.'].get('Refine',False):
566	res = opt.minimize(Min,xstart,bounds=([0.01,1.],[1.2bakMul,0.8bakMul]),
567	method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
568	else:
569	res = opt.minimize(Min,xstart,bounds=([0.,None],[0,1],[0.01,1.]),
570	method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
571	Done(res['x'])
572	if showFit:
573	if res['success']:
574	msg = 'Converged'
575	else:
576	msg = 'Not Converged'
577	if data['Sample Bkg.'].get('Refine',False):
578	G2fil.G2Print(' end: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f}) * {} *\n'.format(
579	data['Ruland'],data['Sample Bkg.']['Mult'],res['fun'],msg))
580	else:
581	G2fil.G2Print(' end: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} RMS:{:.4f}) * {} *\n'.format(
582	data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun'],msg))
583	return res
584
585	def SetupPDFEval(data,xydata,limits,inst,numbDen):
586	Data = copy.deepcopy(data)
587	BkgMax = 1.
588	def EvalLowPDF(arg):
589	'''Objective routine -- evaluates the RMS deviations in G(r)
590	from -4(pi)#densityr for for r<Rmin
591	arguments are ['Flat Bkg','BackRatio','Ruland'] scaled so that
592	the min & max values are between 0 and 1.
593	'''
594	if Data['Sample Bkg.'].get('Refine',False):
595	R,S = arg
596	Data['Sample Bkg.']['Mult'] = S
597	else:
598	F,B,R = arg
599	Data['Flat Bkg'] = BkgMax(2.F-1.)
600	Data['BackRatio'] = B
601	Data['Ruland'] = R
602	CalcPDF(Data,inst,limits,xydata)
603	# test low r computation
604	g = xydata['GofR'][1][1]
605	r = xydata['GofR'][1][0]
606	g0 = g[r < Data['Rmin']] + 4np.pir[r < Data['Rmin']]*numbDen
607	M = sum(g0**2)/len(g0)
608	return M
609	def GetCurrentVals():
610	'''Get the current ['Flat Bkg','BackRatio','Ruland'] with scaling
611	'''
612	if data['Sample Bkg.'].get('Refine',False):
613	return [max(data['Ruland'],.05),data['Sample']['Mult']]
614	try:
615	F = 0.5+0.5*data['Flat Bkg']/BkgMax
616	except:
617	F = 0
618	return [F,data['BackRatio'],max(data['Ruland'],.05)]
619	def SetFinalVals(arg):
620	'''Set the 'Flat Bkg', 'BackRatio' & 'Ruland' values from the
621	scaled, refined values and plot corrected region of G(r)
622	'''
623	if data['Sample Bkg.'].get('Refine',False):
624	R,S = arg
625	data['Sample Bkg.']['Mult'] = S
626	else:
627	F,B,R = arg
628	data['Flat Bkg'] = BkgMax(2.F-1.)
629	data['BackRatio'] = B
630	data['Ruland'] = R
631	CalcPDF(data,inst,limits,xydata)
632	EvalLowPDF(GetCurrentVals())
633	BkgMax = max(xydata['IofQ'][1][1])/50.
634	return EvalLowPDF,GetCurrentVals,SetFinalVals
635
636	#### GSASII convolution peak fitting routines: Finger, Cox & Jephcoat model
637	def factorize(num):
638	''' Provide prime number factors for integer num
639	:returns: dictionary of prime factors (keys) & power for each (data)
640	'''
641	factors = {}
642	orig = num
643
644	# we take advantage of the fact that (i +1)2 = i2 + 2*i +1
645	i, sqi = 2, 4
646	while sqi <= num:
647	while not num%i:
648	num /= i
649	factors[i] = factors.get(i, 0) + 1
650
651	sqi += 2*i + 1
652	i += 1
653
654	if num != 1 and num != orig:
655	factors[num] = factors.get(num, 0) + 1
656
657	if factors:
658	return factors
659	else:
660	return {num:1} #a prime number!
661
662	def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
663	''' Provide list of optimal data sizes for FFT calculations
664
665	:param int nmin: minimum data size >= 1
666	:param int nmax: maximum data size > nmin
667	:param int thresh: maximum prime factor allowed
668	:Returns: list of data sizes where the maximum prime factor is < thresh
669	'''
670	plist = []
671	nmin = max(1,nmin)
672	nmax = max(nmin+1,nmax)
673	for p in range(nmin,nmax):
674	if max(list(factorize(p).keys())) < thresh:
675	plist.append(p)
676	return plist
677
678	np.seterr(divide='ignore')
679
680	# Normal distribution
681
682	# loc = mu, scale = std
683	_norm_pdf_C = 1./math.sqrt(2*math.pi)
684	class norm_gen(st.rv_continuous):
685	'needs a doc string'
686
687	def pdf(self,x,args,*kwds):
688	loc,scale=kwds['loc'],kwds['scale']
689	x = (x-loc)/scale
690	return np.exp(-x*2/2.0) _norm_pdf_C / scale
691
692	norm = norm_gen(name='norm',longname='A normal',extradoc="""
693
694	Normal distribution
695
696	The location (loc) keyword specifies the mean.
697	The scale (scale) keyword specifies the standard deviation.
698
699	normal.pdf(x) = exp(-x*2/2)/sqrt(2pi)
700	""")
701
702	## Cauchy
703
704	# median = loc
705
706	class cauchy_gen(st.rv_continuous):
707	'needs a doc string'
708
709	def pdf(self,x,args,*kwds):
710	loc,scale=kwds['loc'],kwds['scale']
711	x = (x-loc)/scale
712	return 1.0/np.pi/(1.0+x*x) / scale
713
714	cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
715
716	Cauchy distribution
717
718	cauchy.pdf(x) = 1/(pi(1+x*2))
719
720	This is the t distribution with one degree of freedom.
721	""")
722
723
724	class fcjde_gen(st.rv_continuous):
725	"""
726	Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
727	Ref: J. Appl. Cryst. (1994) 27, 892-900.
728
729	:param x: array -1 to 1
730	:param t: 2-theta position of peak
731	:param s: sum(S/L,H/L); S: sample height, H: detector opening,
732	L: sample to detector opening distance
733	:param dx: 2-theta step size in deg
734
735	:returns: for fcj.pdf
736
737	* T = x*dx+t
738	* s = S/L+H/L
739	* if x < 0::
740
741	fcj.pdf = [1/sqrt({cos(T)2/cos(t)2}-1) - 1/s]/\|cos(T)\|
742
743	* if x >= 0: fcj.pdf = 0
744
745	"""
746	def _pdf(self,x,t,s,dx):
747	T = dx*x+t
748	ax2 = abs(npcosd(T))
749	ax = ax2**2
750	bx = npcosd(t)**2
751	bx = np.where(ax>bx,bx,ax)
752	fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
753	fx = np.where(fx > 0.,fx,0.0)
754	return fx
755
756	def pdf(self,x,args,*kwds):
757	loc=kwds['loc']
758	return self._pdf(x-loc,*args)
759
760	fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
761
762	def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
763	'''Compute the Finger-Cox-Jepcoat modified Voigt function for a
764	CW powder peak by direct convolution. This version is not used.
765	'''
766	DX = xdata[1]-xdata[0]
767	widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
768	x = np.linspace(pos-fmin,pos+fmin,256)
769	dx = x[1]-x[0]
770	Norm = norm.pdf(x,loc=pos,scale=widths[0])
771	Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
772	arg = [pos,shl/57.2958,dx,]
773	FCJ = fcjde.pdf(x,*arg,loc=pos)
774	if len(np.nonzero(FCJ)[0])>5:
775	z = np.column_stack([Norm,Cauchy,FCJ]).T
776	Z = fft.fft(z)
777	Df = fft.ifft(Z.prod(axis=0)).real
778	else:
779	z = np.column_stack([Norm,Cauchy]).T
780	Z = fft.fft(z)
781	Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real
782	Df /= np.sum(Df)
783	Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
784	return intensDf(xdata)DX/dx
785
786	#### GSASII peak fitting routine: Finger, Cox & Jephcoat model
787
788	def getWidthsCW(pos,sig,gam,shl):
789	'''Compute the peak widths used for computing the range of a peak
790	for constant wavelength data. On low-angle side, 50 FWHM are used,
791	on high-angle side 75 are used, high angle side extended for axial divergence
792	(for peaks above 90 deg, these are reversed.)
793
794	:param pos: peak position; 2-theta in degrees
795	:param sig: Gaussian peak variance in centideg^2
796	:param gam: Lorentzian peak width in centidegrees
797	:param shl: axial divergence parameter (S+H)/L
798
799	:returns: widths; [Gaussian sigma, Lorentzian gamma] in degrees, and
800	low angle, high angle ends of peak; 20 FWHM & 50 FWHM from position
801	reversed for 2-theta > 90 deg.
802	'''
803	widths = [np.sqrt(sig)/100.,gam/100.]
804	fwhm = 2.355*widths[0]+widths[1]
805	fmin = 50.(fwhm+shlabs(npcosd(pos)))
806	fmax = 75.0*fwhm
807	if pos > 90:
808	fmin,fmax = [fmax,fmin]
809	return widths,fmin,fmax
810
811	def getWidthsED(pos,sig):
812	'''Compute the peak widths used for computing the range of a peak
813	for energy dispersive data. On low-energy side, 20 FWHM are used,
814	on high-energy side 20 are used
815
816	:param pos: peak position; energy in keV (not used)
817	:param sig: Gaussian peak variance in keV^2
818
819	:returns: widths; [Gaussian sigma] in keV, and
820	low angle, high angle ends of peak; 20 FWHM & 50 FWHM from position
821	'''
822	widths = [np.sqrt(sig),.001]
823	fwhm = 2.355*widths[0]
824	fmin = 5.*fwhm
825	fmax = 5.*fwhm
826	return widths,fmin,fmax
827
828	def getWidthsTOF(pos,alp,bet,sig,gam):
829	'''Compute the peak widths used for computing the range of a peak
830	for constant wavelength data. 50 FWHM are used on both sides each
831	extended by exponential coeff.
832
833	param pos: peak position; TOF in musec (not used)
834	param alp,bet: TOF peak exponential rise & decay parameters
835	param sig: Gaussian peak variance in musec^2
836	param gam: Lorentzian peak width in musec
837
838	returns: widths; [Gaussian sigma, Lornetzian gamma] in musec
839	returns: low TOF, high TOF ends of peak; 50FWHM from position
840	'''
841	widths = [np.sqrt(sig),gam]
842	fwhm = 2.355widths[0]+2.widths[1]
843	fmin = 50.fwhm(1.+1./alp)
844	fmax = 50.fwhm(1.+1./bet)
845	return widths,fmin,fmax
846
847	def getFWHM(pos,Inst):
848	'''Compute total FWHM from Thompson, Cox & Hastings (1987) , J. Appl. Cryst. 20, 79-83
849	via getgamFW(g,s).
850
851	:param pos: float peak position in deg 2-theta or tof in musec
852	:param Inst: dict instrument parameters
853
854	:returns float: total FWHM of pseudoVoigt in deg or musec
855	'''
856
857	sig = lambda Th,U,V,W: np.sqrt(max(0.001,Utand(Th)2+Vtand(Th)+W))
858	sigED = lambda E,A,B,C: np.sqrt(max(0.001,AE2+BE+C))
859	sigTOF = lambda dsp,S0,S1,S2,Sq: np.sqrt(S0+S1dsp2+S2dsp*4+Sqdsp)
860	gam = lambda Th,X,Y,Z: Z+X/cosd(Th)+Y*tand(Th)
861	gamTOF = lambda dsp,X,Y,Z: Z+Xdsp+Ydsp**2
862	alpTOF = lambda dsp,alp: alp/dsp
863	betTOF = lambda dsp,bet0,bet1,betq: bet0+bet1/dsp4+betq/dsp2
864	alpPink = lambda pos,alp0,alp1: alp0+alp1*tand(pos/2.)
865	betPink = lambda pos,bet0,bet1: bet0+bet1*tand(pos/2.)
866	if 'LF' in Inst['Type'][0]:
867	return 3
868	elif 'T' in Inst['Type'][0]:
869	dsp = pos/Inst['difC'][1]
870	alp = alpTOF(dsp,Inst['alpha'][1])
871	bet = betTOF(dsp,Inst['beta-0'][1],Inst['beta-1'][1],Inst['beta-q'][1])
872	s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
873	g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
874	return getgamFW(g,s)+np.log(2.0)(alp+bet)/(alpbet)
875	elif 'C' in Inst['Type'][0]:
876	s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
877	g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
878	return getgamFW(g,s)/100. #returns FWHM in deg
879	elif 'E' in Inst['Type'][0]:
880	s = sigED(pos,Inst['A'][1],Inst['B'][1],Inst['C'][1])
881	return 2.35482*s
882	else: #'B'
883	alp = alpPink(pos,Inst['alpha-0'][1],Inst['alpha-1'][1])
884	bet = betPink(pos,Inst['beta-0'][1],Inst['beta-1'][1])
885	s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
886	g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
887	return getgamFW(g,s)/100.+np.log(2.0)(alp+bet)/(alpbet) #returns FWHM in deg
888
889	def getgamFW(g,s):
890	'''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
891	lambda fxn needs FWHM for both Gaussian & Lorentzian components
892
893	:param g: float Lorentzian gamma = FWHM(L)
894	:param s: float Gaussian sig
895
896	:returns float: total FWHM of pseudoVoigt
897	'''
898	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.)
899	return gamFW(2.35482s,g) #sqrt(8ln2)sig = FWHM(G)
900
901	def getBackground(pfx,parmDict,bakType,dataType,xdata,fixback=None):
902	'''Computes the background from vars pulled from gpx file or tree.
903	'''
904	if 'T' in dataType:
905	q = 2.np.piparmDict[pfx+'difC']/xdata
906	elif 'E' in dataType:
907	const = 4.np.pinpsind(parmDict[pfx+'2-theta']/2.0)
908	q = const*xdata
909	else:
910	wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
911	q = npT2q(xdata,wave)
912	yb = np.zeros_like(xdata)
913	nBak = 0
914	sumBk = [0.,0.,0]
915	while True:
916	key = pfx+'Back;'+str(nBak)
917	if key in parmDict:
918	nBak += 1
919	else:
920	break
921	#empirical functions
922	if bakType in ['chebyschev','cosine','chebyschev-1']:
923	dt = xdata[-1]-xdata[0]
924	for iBak in range(nBak):
925	key = pfx+'Back;'+str(iBak)
926	if bakType == 'chebyschev':
927	ybi = parmDict[key](-1.+2.(xdata-xdata[0])/dt)**iBak
928	elif bakType == 'chebyschev-1':
929	xpos = -1.+2.*(xdata-xdata[0])/dt
930	ybi = parmDict[key]np.cos(iBaknp.arccos(xpos))
931	elif bakType == 'cosine':
932	ybi = parmDict[key]npcosd(180.xdata*iBak/xdata[-1])
933	yb += ybi
934	sumBk[0] = np.sum(yb)
935	elif bakType in ['Q^2 power series','Q^-2 power series']:
936	QT = 1.
937	yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
938	for iBak in range(nBak-1):
939	key = pfx+'Back;'+str(iBak+1)
940	if '-2' in bakType:
941	QT = (iBak+1)q**-2
942	else:
943	QT = q*2/(iBak+1)
944	yb += QT*parmDict[key]
945	sumBk[0] = np.sum(yb)
946	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
947	if nBak == 1:
948	yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
949	elif nBak == 2:
950	dX = xdata[-1]-xdata[0]
951	T2 = (xdata-xdata[0])/dX
952	T1 = 1.0-T2
953	yb = parmDict[pfx+'Back;0']T1+parmDict[pfx+'Back;1']T2
954	else:
955	xnomask = ma.getdata(xdata)
956	xmin,xmax = xnomask[0],xnomask[-1]
957	if bakType == 'lin interpolate':
958	bakPos = np.linspace(xmin,xmax,nBak,True)
959	elif bakType == 'inv interpolate':
960	bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
961	elif bakType == 'log interpolate':
962	bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
963	bakPos[0] = xmin
964	bakPos[-1] = xmax
965	bakVals = np.zeros(nBak)
966	for i in range(nBak):
967	bakVals[i] = parmDict[pfx+'Back;'+str(i)]
968	bakInt = si.interp1d(bakPos,bakVals,'linear')
969	yb = bakInt(ma.getdata(xdata))
970	sumBk[0] = np.sum(yb)
971	#Debye function
972	if pfx+'difC' in parmDict or 'E' in dataType:
973	ff = 1.
974	else:
975	try:
976	wave = parmDict[pfx+'Lam']
977	except KeyError:
978	wave = parmDict[pfx+'Lam1']
979	SQ = (q/(4.np.pi))*2
980	FF = G2elem.GetFormFactorCoeff('Si')[0]
981	ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
982	iD = 0
983	while True:
984	try:
985	dbA = parmDict[pfx+'DebyeA;'+str(iD)]
986	dbR = parmDict[pfx+'DebyeR;'+str(iD)]
987	dbU = parmDict[pfx+'DebyeU;'+str(iD)]
988	ybi = ffdbAnp.sin(qdbR)np.exp(-dbUq2)/(qdbR)
989	yb += ybi
990	sumBk[1] += np.sum(ybi)
991	iD += 1
992	except KeyError:
993	break
994	#peaks
995	iD = 0
996	while True:
997	try:
998	pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
999	pkI = max(parmDict[pfx+'BkPkint;'+str(iD)],0.1)
1000	pkS = max(parmDict[pfx+'BkPksig;'+str(iD)],0.01)
1001	pkG = max(parmDict[pfx+'BkPkgam;'+str(iD)],0.1)
1002	if 'C' in dataType:
1003	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
1004	elif 'E' in dataType:
1005	Wd,fmin,fmax = getWidthsED(pkP,pkS)
1006	else: #'T'OF
1007	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
1008	iBeg = np.searchsorted(xdata,pkP-fmin)
1009	iFin = np.searchsorted(xdata,pkP+fmax)
1010	lenX = len(xdata)
1011	if not iBeg:
1012	iFin = np.searchsorted(xdata,pkP+fmax)
1013	elif iBeg == lenX:
1014	iFin = iBeg
1015	else:
1016	iFin = np.searchsorted(xdata,pkP+fmax)
1017	if 'C' in dataType:
1018	ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])[0]
1019	elif 'T' in dataType:
1020	ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])[0]
1021	elif 'B' in dataType:
1022	ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS/100.,pkG/1.e4,xdata[iBeg:iFin])[0]
1023	elif 'E' in dataType:
1024	ybi = pkIgetPsVoigt(pkP,pkS10.*4,pkG100.,xdata[iBeg:iFin])[0]
1025	else:
1026	raise Exception('dataType of {:} should not happen!'.format(dataType))
1027	yb[iBeg:iFin] += ybi
1028	sumBk[2] += np.sum(ybi)
1029	iD += 1
1030	except KeyError:
1031	break
1032	except ValueError:
1033	G2fil.G2Print ('** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **')
1034	break
1035	if fixback is not None:
1036	yb += parmDict[pfx+'BF mult']*fixback
1037	sumBk[0] = sum(yb)
1038	return yb,sumBk
1039
1040	def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata,fixback=None):
1041	'needs a doc string'
1042	if 'T' in dataType:
1043	q = 2.np.piparmDict[hfx+'difC']/xdata
1044	elif 'E' in dataType:
1045	const = 4.np.pinpsind(parmDict[hfx+'2-theta']/2.0)
1046	q = const*xdata
1047	else:
1048	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1049	q = 2.np.pinpsind(xdata/2.)/wave
1050	nBak = 0
1051	while True:
1052	key = hfx+'Back;'+str(nBak)
1053	if key in parmDict:
1054	nBak += 1
1055	else:
1056	break
1057	dydb = np.zeros(shape=(nBak,len(xdata)))
1058	dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
1059	dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
1060	dydfb = []
1061
1062	if bakType in ['chebyschev','cosine','chebyschev-1']:
1063	dt = xdata[-1]-xdata[0]
1064	for iBak in range(nBak):
1065	if bakType == 'chebyschev':
1066	dydb[iBak] = (-1.+2.(xdata-xdata[0])/dt)*iBak
1067	elif bakType == 'chebyschev-1':
1068	xpos = -1.+2.*(xdata-xdata[0])/dt
1069	dydb[iBak] = np.cos(iBak*np.arccos(xpos))
1070	elif bakType == 'cosine':
1071	dydb[iBak] = npcosd(180.xdataiBak/xdata[-1])
1072	elif bakType in ['Q^2 power series','Q^-2 power series']:
1073	QT = 1.
1074	dydb[0] = np.ones_like(xdata)
1075	for iBak in range(nBak-1):
1076	if '-2' in bakType:
1077	QT = (iBak+1)q**-2
1078	else:
1079	QT = q*2/(iBak+1)
1080	dydb[iBak+1] = QT
1081	elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
1082	if nBak == 1:
1083	dydb[0] = np.ones_like(xdata)
1084	elif nBak == 2:
1085	dX = xdata[-1]-xdata[0]
1086	T2 = (xdata-xdata[0])/dX
1087	T1 = 1.0-T2
1088	dydb = [T1,T2]
1089	else:
1090	xnomask = ma.getdata(xdata)
1091	xmin,xmax = xnomask[0],xnomask[-1]
1092	if bakType == 'lin interpolate':
1093	bakPos = np.linspace(xmin,xmax,nBak,True)
1094	elif bakType == 'inv interpolate':
1095	bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
1096	elif bakType == 'log interpolate':
1097	bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
1098	bakPos[0] = xmin
1099	bakPos[-1] = xmax
1100	for i,pos in enumerate(bakPos):
1101	if i == 0:
1102	dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
1103	elif i == len(bakPos)-1:
1104	dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
1105	else:
1106	dydb[i] = np.where(xdata>bakPos[i],
1107	np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
1108	np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
1109	if hfx+'difC' in parmDict:
1110	ff = 1.
1111	else:
1112	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1113	q = npT2q(xdata,wave)
1114	SQ = (q/(4np.pi))*2
1115	FF = G2elem.GetFormFactorCoeff('Si')[0]
1116	ff = np.array(G2elem.ScatFac(FF,SQ)[0])np.pi*2 #needs pi^2~10. for cw data (why?)
1117	iD = 0
1118	while True:
1119	try:
1120	if hfx+'difC' in parmDict:
1121	q = 2np.piparmDict[hfx+'difC']/xdata
1122	dbA = parmDict[hfx+'DebyeA;'+str(iD)]
1123	dbR = parmDict[hfx+'DebyeR;'+str(iD)]
1124	dbU = parmDict[hfx+'DebyeU;'+str(iD)]
1125	sqr = np.sin(qdbR)/(qdbR)
1126	cqr = np.cos(q*dbR)
1127	temp = np.exp(-dbUq*2)
1128	dyddb[3iD] = ffsqr*temp
1129	dyddb[3iD+1] = ffdbAtemp(cqr-sqr)/(dbR)
1130	dyddb[3iD+2] = -ffdbAsqrtempq*2
1131	iD += 1
1132	except KeyError:
1133	break
1134	iD = 0
1135	while True:
1136	try:
1137	pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
1138	pkI = max(parmDict[hfx+'BkPkint;'+str(iD)],0.1)
1139	pkS = max(parmDict[hfx+'BkPksig;'+str(iD)],0.01)
1140	pkG = max(parmDict[hfx+'BkPkgam;'+str(iD)],0.1)
1141	if 'C' in dataType:
1142	Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
1143	elif 'E' in dataType:
1144	Wd,fmin,fmax = getWidthsED(pkP,pkS)
1145	else: #'T' or 'B'
1146	Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
1147	iBeg = np.searchsorted(xdata,pkP-fmin)
1148	iFin = np.searchsorted(xdata,pkP+fmax)
1149	lenX = len(xdata)
1150	if not iBeg:
1151	iFin = np.searchsorted(xdata,pkP+fmax)
1152	elif iBeg == lenX:
1153	iFin = iBeg
1154	else:
1155	iFin = np.searchsorted(xdata,pkP+fmax)
1156	if 'C' in dataType:
1157	Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
1158	elif 'E' in dataType:
1159	Df,dFdp,dFds,dFdg = getdPsVoigt(pkP,pkS10.4,pkG100.,xdata[iBeg:iFin])
1160	else: #'T'OF
1161	Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
1162	dydpk[4iD][iBeg:iFin] += pkIdFdp
1163	dydpk[4*iD+1][iBeg:iFin] += Df
1164	dydpk[4iD+2][iBeg:iFin] += pkIdFds
1165	dydpk[4iD+3][iBeg:iFin] += pkIdFdg
1166	iD += 1
1167	except KeyError:
1168	break
1169	except ValueError:
1170	G2fil.G2Print ('** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again **')
1171	break
1172	# fixed background from file
1173	if fixback is not None:
1174	dydfb = fixback
1175	return dydb,dyddb,dydpk,dydfb
1176
1177	#### Using old gsas fortran routines for powder peak shapes & derivatives
1178	def getFCJVoigt3(pos,sig,gam,shl,xdata):
1179	'''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
1180	CW powder peak in external Fortran routine
1181
1182	param pos: peak position in degrees
1183	param sig: Gaussian variance in centideg^2
1184	param gam: Lorentzian width in centideg
1185	param shl: axial divergence parameter (S+H)/L
1186	param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv)
1187
1188	returns: array: calculated peak function at each xdata
1189	returns: integral of peak; nominally = 1.0
1190	'''
1191	if len(xdata):
1192	cw = np.diff(xdata)
1193	cw = np.append(cw,cw[-1])
1194	Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1195	return Df,np.sum(100.Dfcw)
1196	else:
1197	return 0.,1.
1198
1199	def getdFCJVoigt3(pos,sig,gam,shl,xdata):
1200	'''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
1201	function for a CW powder peak
1202
1203	param pos: peak position in degrees
1204	param sig: Gaussian variance in centideg^2
1205	param gam: Lorentzian width in centideg
1206	param shl: axial divergence parameter (S+H)/L
1207	param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv)
1208
1209	returns: arrays: function and derivatives of pos, sig, gam, & shl
1210	'''
1211	Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1212	return Df,dFdp,dFds,dFdg,dFdsh
1213
1214	def getPsVoigt(pos,sig,gam,xdata):
1215	'''Compute the simple Pseudo-Voigt function for a
1216	small angle Bragg peak in external Fortran routine
1217
1218	param pos: peak position in degrees
1219	param sig: Gaussian variance in centideg^2
1220	param gam: Lorentzian width in centideg
1221
1222	returns: array: calculated peak function at each xdata
1223	returns: integral of peak; nominally = 1.0
1224	'''
1225
1226	cw = np.diff(xdata)
1227	cw = np.append(cw,cw[-1])
1228	Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
1229	return Df,np.sum(100.Dfcw)
1230
1231	def getdPsVoigt(pos,sig,gam,xdata):
1232	'''Compute the simple Pseudo-Voigt function derivatives for a
1233	small angle Bragg peak peak in external Fortran routine
1234
1235	param pos: peak position in degrees
1236	param sig: Gaussian variance in centideg^2
1237	param gam: Lorentzian width in centideg
1238
1239	returns: arrays: function and derivatives of pos, sig & gam
1240	NB: the pos derivative has the opposite sign compared to that in other profile functions
1241	'''
1242
1243	Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
1244	return Df,dFdp,dFds,dFdg
1245
1246	def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
1247	'''Compute the double exponential Pseudo-Voigt convolution function for a
1248	neutron TOF & CW pink peak in external Fortran routine
1249	'''
1250
1251	cw = np.diff(xdata)
1252	cw = np.append(cw,cw[-1])
1253	Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1254	return Df,np.sum(Df*cw)
1255
1256	def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
1257	'''Compute the double exponential Pseudo-Voigt convolution function derivatives for a
1258	neutron TOF & CW pink peak in external Fortran routine
1259	'''
1260
1261	Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1262	return Df,dFdp,dFda,dFdb,dFds,dFdg
1263
1264	def ellipseSize(H,Sij,GB):
1265	'''Implements r=1/sqrt(sum((1/S)*(q.v)^2) per note from Alexander Brady
1266	'''
1267
1268	HX = np.inner(H.T,GB)
1269	lenHX = np.sqrt(np.sum(HX**2))
1270	Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))
1271	R = np.inner(HX/lenHX,Rsize)*2Esize #want column length for hkl in crystal
1272	lenR = 1./np.sqrt(np.sum(R))
1273	return lenR
1274
1275	def ellipseSizeDerv(H,Sij,GB):
1276	'''Implements r=1/sqrt(sum((1/S)*(q.v)^2) derivative per note from Alexander Brady
1277	'''
1278
1279	lenR = ellipseSize(H,Sij,GB)
1280	delt = 0.001
1281	dRdS = np.zeros(6)
1282	for i in range(6):
1283	Sij[i] -= delt
1284	lenM = ellipseSize(H,Sij,GB)
1285	Sij[i] += 2.*delt
1286	lenP = ellipseSize(H,Sij,GB)
1287	Sij[i] -= delt
1288	dRdS[i] = (lenP-lenM)/(2.*delt)
1289	return lenR,dRdS
1290
1291	def getMustrain(HKL,G,SGData,muStrData):
1292	if muStrData[0] == 'isotropic':
1293	return np.ones(HKL.shape[1])*muStrData[1][0]
1294	elif muStrData[0] == 'uniaxial':
1295	H = np.array(HKL)
1296	P = np.array(muStrData[3])
1297	cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1298	Si = muStrData[1][0]
1299	Sa = muStrData[1][1]
1300	return SiSa/(np.sqrt((SicosP)*2+(SasinP)**2))
1301	else: #generalized - P.W. Stephens model
1302	H = np.array(HKL)
1303	rdsq = np.array([G2lat.calc_rDsq2(h,G) for h in H.T])
1304	Strms = np.array(G2spc.MustrainCoeff(H,SGData))
1305	Sum = np.sum(np.array(muStrData[4])[:,nxs]*Strms,axis=0)
1306	return np.sqrt(Sum)/rdsq
1307
1308	def getCrSize(HKL,G,GB,sizeData):
1309	if sizeData[0] == 'isotropic':
1310	return np.ones(HKL.shape[1])*sizeData[1][0]
1311	elif sizeData[0] == 'uniaxial':
1312	H = np.array(HKL)
1313	P = np.array(sizeData[3])
1314	cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1315	Si = sizeData[1][0]
1316	Sa = sizeData[1][1]
1317	return SiSa/(np.sqrt((SicosP)*2+(SasinP)**2))
1318	else:
1319	Sij =[sizeData[4][i] for i in range(6)]
1320	H = np.array(HKL)
1321	return 1./np.array([ellipseSize(h,Sij,GB) for h in H.T])**2
1322
1323	def getHKLpeak(dmin,SGData,A,Inst=None,nodup=False):
1324	'''
1325	Generates allowed by symmetry reflections with d >= dmin
1326	NB: GenHKLf & checkMagextc return True for extinct reflections
1327
1328	:param dmin: minimum d-spacing
1329	:param SGData: space group data obtained from SpcGroup
1330	:param A: lattice parameter terms A1-A6
1331	:param Inst: instrument parameter info
1332	:returns: HKLs: np.array hkl, etc for allowed reflections
1333
1334	'''
1335	HKL = G2lat.GenHLaue(dmin,SGData,A)
1336	HKLs = []
1337	ds = []
1338	for h,k,l,d in HKL:
1339	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1340	if ext and 'MagSpGrp' in SGData:
1341	ext = G2spc.checkMagextc([h,k,l],SGData)
1342	if not ext:
1343	if nodup and int(10000*d) in ds:
1344	continue
1345	ds.append(int(10000*d))
1346	if Inst == None:
1347	HKLs.append([h,k,l,d,0,-1])
1348	else:
1349	HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1350	return np.array(HKLs)
1351
1352	def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1353	'needs a doc string'
1354	HKLs = []
1355	vec = np.array(Vec)
1356	vstar = np.sqrt(G2lat.calc_rDsq(vec,A)) #find extra needed for -n SS reflections
1357	dvec = 1./(maxH*vstar+1./dmin)
1358	HKL = G2lat.GenHLaue(dvec,SGData,A)
1359	SSdH = [vec*h for h in range(-maxH,maxH+1)]
1360	SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1361	ifMag = False
1362	if 'MagSpGrp' in SGData:
1363	ifMag = True
1364	for h,k,l,d in HKL:
1365	ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1366	if not ext and d >= dmin:
1367	HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1368	for dH in SSdH:
1369	if dH:
1370	DH = SSdH[dH]
1371	H = [h+DH[0],k+DH[1],l+DH[2]]
1372	d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1373	if d >= dmin:
1374	HKLM = np.array([h,k,l,dH])
1375	if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1376	HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])
1377	return G2lat.sortHKLd(HKLs,True,True,True)
1378
1379	peakInstPrmMode = True
1380	'''Determines the mode used for peak fitting. When peakInstPrmMode=True peak
1381	width parameters are computed from the instrument parameters (UVW,... or
1382	alpha,... etc) unless the individual parameter is refined. This allows the
1383	instrument parameters to be refined. When peakInstPrmMode=False, the instrument
1384	parameters are not used and cannot be refined.
1385	The default is peakFitMode=True. This is changed only in
1386	:func:`setPeakInstPrmMode`, which is called from :mod:`GSASIIscriptable`
1387	or GSASIIphsGUI.OnSetPeakWidMode ('Gen unvaried widths' menu item).
1388	'''
1389
1390	def setPeakInstPrmMode(normal=True):
1391	'''Determines the mode used for peak fitting. If normal=True (default)
1392	peak width parameters are computed from the instrument parameters
1393	unless the individual parameter is refined. If normal=False,
1394	peak widths are used as supplied for each peak.
1395
1396	Note that normal=True unless this routine is called. Also,
1397	instrument parameters can only be refined with normal=True.
1398
1399	:param bool normal: setting to apply to global variable
1400	:data:`peakInstPrmMode`
1401	'''
1402	global peakInstPrmMode
1403	peakInstPrmMode = normal
1404
1405	def getPeakProfile(dataType,parmDict,xdata,fixback,varyList,bakType):
1406	'''Computes the profiles from multiple peaks for individual peak fitting
1407	for powder patterns.
1408	NB: not used for Rietveld refinement
1409	'''
1410
1411	yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0]
1412	yc = np.zeros_like(yb)
1413	if 'LF' in dataType:
1414	if 'Lam1' in parmDict.keys():
1415	lam = parmDict['Lam1']
1416	lam2 = parmDict['Lam2']
1417	Ka2 = True
1418	lamRatio = 360(lam2-lam)/(np.pilam)
1419	kRatio = parmDict['I(L2)/I(L1)']
1420	else:
1421	lam = parmDict['Lam']
1422	Ka2 = False
1423	shol = 0
1424	# loop over peaks
1425	iPeak = -1
1426	try:
1427	ncells = parmDict['ncell']
1428	clat = parmDict['clat']
1429	except KeyError: # no Laue info must be bkg fit
1430	print('Laue Fit: no params, assuming bkg fit')
1431	return yb
1432	while True:
1433	iPeak += 1
1434	try:
1435	#Qcen = 2 * np.pi * lam * (iPeak+1) / parmDict['clat']
1436	l = parmDict['l'+str(iPeak)]
1437	pos = 360 * np.arcsin(0.5 * lam * l / parmDict['clat']) / np.pi
1438	parmDict['pos'+str(iPeak)] = pos
1439	#tth = (pos-parmDict['Zero'])
1440	intens = parmDict['int'+str(iPeak)]
1441	damp = parmDict['damp'+str(iPeak)]
1442	asym = parmDict['asym'+str(iPeak)]
1443	sig = parmDict['sig'+str(iPeak)]
1444	gam = parmDict['gam'+str(iPeak)]
1445	fmin = 8 # for now make peaks 8 degrees wide
1446	fmin = min(0.9*abs(xdata[-1] - xdata[0]),fmin) # unless the data range is smaller
1447	iBeg = np.searchsorted(xdata,pos-fmin/2)
1448	iFin = np.searchsorted(xdata,pos+fmin/2)
1449	if not iBeg+iFin: # skip peak below low limit
1450	continue
1451	elif not iBeg-iFin: # got peak above high limit (peaks sorted, so we can stop)
1452	break
1453	#LF.plotme(fmin,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym)
1454	LaueFringePeakCalc(xdata,yc,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym,fmin)
1455	if Ka2:
1456	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1457	iBeg = np.searchsorted(xdata,pos2-fmin)
1458	iFin = np.searchsorted(xdata,pos2+fmin)
1459	if iBeg-iFin:
1460	LaueFringePeakCalc(xdata,yc,lam2,pos2,intens*kRatio,sig,gam,shol,ncells,clat,damp,asym,fmin)
1461	except KeyError: #no more peaks to process
1462	return yb+yc
1463	elif 'C' in dataType:
1464	shl = max(parmDict['SH/L'],0.002)
1465	Ka2 = False
1466	if 'Lam1' in parmDict.keys():
1467	Ka2 = True
1468	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1469	kRatio = parmDict['I(L2)/I(L1)']
1470	iPeak = 0
1471	while True:
1472	try:
1473	pos = parmDict['pos'+str(iPeak)]
1474	tth = (pos-parmDict['Zero'])
1475	intens = parmDict['int'+str(iPeak)]
1476	sigName = 'sig'+str(iPeak)
1477	if sigName in varyList or not peakInstPrmMode:
1478	sig = parmDict[sigName]
1479	else:
1480	sig = G2mth.getCWsig(parmDict,tth)
1481	sig = max(sig,0.001) #avoid neg sigma^2
1482	gamName = 'gam'+str(iPeak)
1483	if gamName in varyList or not peakInstPrmMode:
1484	gam = parmDict[gamName]
1485	else:
1486	gam = G2mth.getCWgam(parmDict,tth)
1487	gam = max(gam,0.001) #avoid neg gamma
1488	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1489	iBeg = np.searchsorted(xdata,pos-fmin)
1490	iFin = np.searchsorted(xdata,pos+fmin)
1491	if not iBeg+iFin: #peak below low limit
1492	iPeak += 1
1493	continue
1494	elif not iBeg-iFin: #peak above high limit
1495	return yb+yc
1496	fp = getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])[0]
1497	yc[iBeg:iFin] += intens*fp
1498	if Ka2:
1499	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1500	iBeg = np.searchsorted(xdata,pos2-fmin)
1501	iFin = np.searchsorted(xdata,pos2+fmin)
1502	if iBeg-iFin:
1503	fp2 = getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])[0]
1504	yc[iBeg:iFin] += intenskRatiofp2
1505	iPeak += 1
1506	except KeyError: #no more peaks to process
1507	return yb+yc
1508	elif 'E' in dataType:
1509	iPeak = 0
1510	dsp = 1.0 #for now - fix later
1511	while True:
1512	try:
1513	pos = parmDict['pos'+str(iPeak)]
1514	intens = parmDict['int'+str(iPeak)]
1515	sigName = 'sig'+str(iPeak)
1516	if sigName in varyList or not peakInstPrmMode:
1517	sig = parmDict[sigName]
1518	else:
1519	sig = G2mth.getEDsig(parmDict,pos)
1520	sig = max(sig,0.001) #avoid neg sigma^2
1521	Wd,fmin,fmax = getWidthsED(pos,sig)
1522	iBeg = np.searchsorted(xdata,pos-fmin)
1523	iFin = max(iBeg+3,np.searchsorted(xdata,pos+fmin))
1524	if not iBeg+iFin: #peak below low limit
1525	iPeak += 1
1526	continue
1527	elif not iBeg-iFin: #peak above high limit
1528	return yb+yc
1529	yc[iBeg:iFin] += intensgetPsVoigt(pos,sig10.**4,0.001,xdata[iBeg:iFin])[0]
1530	iPeak += 1
1531	except KeyError: #no more peaks to process
1532	return yb+yc
1533	elif 'B' in dataType:
1534	iPeak = 0
1535	dsp = 1.0 #for now - fix later
1536	while True:
1537	try:
1538	pos = parmDict['pos'+str(iPeak)]
1539	tth = (pos-parmDict['Zero'])
1540	intens = parmDict['int'+str(iPeak)]
1541	alpName = 'alp'+str(iPeak)
1542	if alpName in varyList or not peakInstPrmMode:
1543	alp = parmDict[alpName]
1544	else:
1545	alp = G2mth.getPinkalpha(parmDict,tth)
1546	alp = max(0.1,alp)
1547	betName = 'bet'+str(iPeak)
1548	if betName in varyList or not peakInstPrmMode:
1549	bet = parmDict[betName]
1550	else:
1551	bet = G2mth.getPinkbeta(parmDict,tth)
1552	bet = max(0.1,bet)
1553	sigName = 'sig'+str(iPeak)
1554	if sigName in varyList or not peakInstPrmMode:
1555	sig = parmDict[sigName]
1556	else:
1557	sig = G2mth.getCWsig(parmDict,tth)
1558	sig = max(sig,0.001) #avoid neg sigma^2
1559	gamName = 'gam'+str(iPeak)
1560	if gamName in varyList or not peakInstPrmMode:
1561	gam = parmDict[gamName]
1562	else:
1563	gam = G2mth.getCWgam(parmDict,tth)
1564	gam = max(gam,0.001) #avoid neg gamma
1565	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1566	iBeg = np.searchsorted(xdata,pos-fmin)
1567	iFin = np.searchsorted(xdata,pos+fmin)
1568	if not iBeg+iFin: #peak below low limit
1569	iPeak += 1
1570	continue
1571	elif not iBeg-iFin: #peak above high limit
1572	return yb+yc
1573	yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])[0]
1574	iPeak += 1
1575	except KeyError: #no more peaks to process
1576	return yb+yc
1577	else:
1578	Pdabc = parmDict['Pdabc']
1579	difC = parmDict['difC']
1580	iPeak = 0
1581	while True:
1582	try:
1583	pos = parmDict['pos'+str(iPeak)]
1584	tof = pos-parmDict['Zero']
1585	dsp = tof/difC
1586	intens = parmDict['int'+str(iPeak)]
1587	alpName = 'alp'+str(iPeak)
1588	if alpName in varyList or not peakInstPrmMode:
1589	alp = parmDict[alpName]
1590	else:
1591	if len(Pdabc):
1592	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1593	else:
1594	alp = G2mth.getTOFalpha(parmDict,dsp)
1595	alp = max(0.1,alp)
1596	betName = 'bet'+str(iPeak)
1597	if betName in varyList or not peakInstPrmMode:
1598	bet = parmDict[betName]
1599	else:
1600	if len(Pdabc):
1601	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1602	else:
1603	bet = G2mth.getTOFbeta(parmDict,dsp)
1604	bet = max(0.0001,bet)
1605	sigName = 'sig'+str(iPeak)
1606	if sigName in varyList or not peakInstPrmMode:
1607	sig = parmDict[sigName]
1608	else:
1609	sig = G2mth.getTOFsig(parmDict,dsp)
1610	gamName = 'gam'+str(iPeak)
1611	if gamName in varyList or not peakInstPrmMode:
1612	gam = parmDict[gamName]
1613	else:
1614	gam = G2mth.getTOFgamma(parmDict,dsp)
1615	gam = max(gam,0.001) #avoid neg gamma
1616	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1617	iBeg = np.searchsorted(xdata,pos-fmin)
1618	iFin = np.searchsorted(xdata,pos+fmax)
1619	lenX = len(xdata)
1620	if not iBeg:
1621	iFin = np.searchsorted(xdata,pos+fmax)
1622	elif iBeg == lenX:
1623	iFin = iBeg
1624	else:
1625	iFin = np.searchsorted(xdata,pos+fmax)
1626	if not iBeg+iFin: #peak below low limit
1627	iPeak += 1
1628	continue
1629	elif not iBeg-iFin: #peak above high limit
1630	return yb+yc
1631	yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])[0]
1632	iPeak += 1
1633	except KeyError: #no more peaks to process
1634	return yb+yc
1635
1636	def getPeakProfileDerv(dataType,parmDict,xdata,fixback,varyList,bakType):
1637	'''Computes the profile derivatives for a powder pattern for single peak fitting
1638
1639	return: np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1640
1641	NB: not used for Rietveld refinement
1642	'''
1643	dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1644	dMdb,dMddb,dMdpk,dMdfb = getBackgroundDerv('',parmDict,bakType,dataType,xdata,fixback)
1645	if 'Back;0' in varyList: #background derivs are in front if present
1646	dMdv[0:len(dMdb)] = dMdb
1647	names = ['DebyeA','DebyeR','DebyeU']
1648	for name in varyList:
1649	if 'Debye' in name:
1650	parm,Id = name.split(';')
1651	ip = names.index(parm)
1652	dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip]
1653	names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1654	for name in varyList:
1655	if 'BkPk' in name:
1656	parm,Id = name.split(';')
1657	ip = names.index(parm)
1658	dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
1659	if 'LF' in dataType:
1660	for i,name in enumerate(varyList):
1661	if not np.all(dMdv[i] == 0): continue
1662	deltaParmDict = parmDict.copy()
1663	delta = max(parmDict[name]/1e5,0.001)
1664	deltaParmDict[name] += delta
1665	#print('num. deriv for',name,'val',deltaParmDict[name],'delta',delta)
1666	intArrP = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType)
1667	deltaParmDict[name] -= 2*delta
1668	intArrM = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType)
1669	dMdv[i] = 0.5 * (intArrP - intArrM) / delta
1670	return dMdv
1671	if 'C' in dataType:
1672	shl = max(parmDict['SH/L'],0.002)
1673	Ka2 = False
1674	if 'Lam1' in parmDict.keys():
1675	Ka2 = True
1676	lamRatio = 360(parmDict['Lam2']-parmDict['Lam1'])/(np.piparmDict['Lam1'])
1677	kRatio = parmDict['I(L2)/I(L1)']
1678	iPeak = 0
1679	while True:
1680	try:
1681	pos = parmDict['pos'+str(iPeak)]
1682	tth = (pos-parmDict['Zero'])
1683	intens = parmDict['int'+str(iPeak)]
1684	sigName = 'sig'+str(iPeak)
1685	if sigName in varyList or not peakInstPrmMode:
1686	sig = parmDict[sigName]
1687	dsdU = dsdV = dsdW = 0
1688	else:
1689	sig = G2mth.getCWsig(parmDict,tth)
1690	dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1691	sig = max(sig,0.001) #avoid neg sigma
1692	gamName = 'gam'+str(iPeak)
1693	if gamName in varyList or not peakInstPrmMode:
1694	gam = parmDict[gamName]
1695	dgdX = dgdY = dgdZ = 0
1696	else:
1697	gam = G2mth.getCWgam(parmDict,tth)
1698	dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1699	gam = max(gam,0.001) #avoid neg gamma
1700	Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1701	iBeg = np.searchsorted(xdata,pos-fmin)
1702	iFin = max(iBeg+3,np.searchsorted(xdata,pos+fmin))
1703	if not iBeg+iFin: #peak below low limit
1704	iPeak += 1
1705	continue
1706	elif not iBeg-iFin: #peak above high limit
1707	break
1708	dMdpk = np.zeros(shape=(6,len(xdata)))
1709	dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1710	for i in range(1,5):
1711	dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1712	dMdpk[0][iBeg:iFin] += dMdipk[0]
1713	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1714	if Ka2:
1715	pos2 = pos+lamRatiotand(pos/2.0) # + 360/pi Dlam/lam * tan(th)
1716	iBeg = np.searchsorted(xdata,pos2-fmin)
1717	iFin = np.searchsorted(xdata,pos2+fmin)
1718	if iBeg-iFin:
1719	dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1720	for i in range(1,5):
1721	dMdpk[i][iBeg:iFin] += intenskRatiodMdipk2[i]
1722	dMdpk[0][iBeg:iFin] += kRatio*dMdipk2[0]
1723	dMdpk[5][iBeg:iFin] += dMdipk2[0]
1724	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1725	for parmName in ['pos','int','sig','gam']:
1726	try:
1727	idx = varyList.index(parmName+str(iPeak))
1728	dMdv[idx] = dervDict[parmName]
1729	except ValueError:
1730	pass
1731	if 'U' in varyList:
1732	dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1733	if 'V' in varyList:
1734	dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1735	if 'W' in varyList:
1736	dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1737	if 'X' in varyList:
1738	dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1739	if 'Y' in varyList:
1740	dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1741	if 'Z' in varyList:
1742	dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1743	if 'SH/L' in varyList:
1744	dMdv[varyList.index('SH/L')] += dervDict['shl'] #problem here
1745	if 'I(L2)/I(L1)' in varyList:
1746	dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1747	iPeak += 1
1748	except KeyError: #no more peaks to process
1749	break
1750	elif 'E' in dataType:
1751	iPeak = 0
1752	while True:
1753	try:
1754	pos = parmDict['pos'+str(iPeak)]
1755	intens = parmDict['int'+str(iPeak)]
1756	sigName = 'sig'+str(iPeak)
1757	if sigName in varyList or not peakInstPrmMode:
1758	sig = parmDict[sigName]
1759	dsdA = dsdB = dsdC = 0
1760	else:
1761	sig = G2mth.getEDsig(parmDict,pos)
1762	dsdA,dsdB,dsdC = G2mth.getEDsigDeriv(parmDict,pos)
1763	sig = max(sig,0.001) #avoid neg sigma
1764	Wd,fmin,fmax = getWidthsED(pos,sig)
1765	iBeg = np.searchsorted(xdata,pos-fmin)
1766	iFin = np.searchsorted(xdata,pos+fmin)
1767	if not iBeg+iFin: #peak below low limit
1768	iPeak += 1
1769	continue
1770	elif not iBeg-iFin: #peak above high limit
1771	break
1772	dMdpk = np.zeros(shape=(4,len(xdata)))
1773	dMdipk = getdPsVoigt(pos,sig10.*4,0.001,xdata[iBeg:iFin])
1774	dMdpk[0][iBeg:iFin] += dMdipk[0]
1775	for i in range(1,4):
1776	dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1777	dervDict = {'int':dMdpk[0],'pos':-dMdpk[1],'sig':dMdpk[2]10*4}
1778	for parmName in ['pos','int','sig']:
1779	try:
1780	idx = varyList.index(parmName+str(iPeak))
1781	dMdv[idx] = dervDict[parmName]
1782	except ValueError:
1783	pass
1784	if 'A' in varyList:
1785	dMdv[varyList.index('A')] += dsdA*dervDict['sig']
1786	if 'B' in varyList:
1787	dMdv[varyList.index('B')] += dsdB*dervDict['sig']
1788	if 'C' in varyList:
1789	dMdv[varyList.index('C')] += dsdC*dervDict['sig']
1790	iPeak += 1
1791	except KeyError: #no more peaks to process
1792	break
1793
1794	elif 'B' in dataType:
1795	iPeak = 0
1796	while True:
1797	try:
1798	pos = parmDict['pos'+str(iPeak)]
1799	tth = (pos-parmDict['Zero'])
1800	intens = parmDict['int'+str(iPeak)]
1801	alpName = 'alp'+str(iPeak)
1802	if alpName in varyList or not peakInstPrmMode:
1803	alp = parmDict[alpName]
1804	dada0 = dada1 = 0.0
1805	else:
1806	alp = G2mth.getPinkalpha(parmDict,tth)
1807	dada0,dada1 = G2mth.getPinkalphaDeriv(tth)
1808	alp = max(0.0001,alp)
1809	betName = 'bet'+str(iPeak)
1810	if betName in varyList or not peakInstPrmMode:
1811	bet = parmDict[betName]
1812	dbdb0 = dbdb1 = 0.0
1813	else:
1814	bet = G2mth.getPinkbeta(parmDict,tth)
1815	dbdb0,dbdb1 = G2mth.getPinkbetaDeriv(tth)
1816	bet = max(0.0001,bet)
1817	sigName = 'sig'+str(iPeak)
1818	if sigName in varyList or not peakInstPrmMode:
1819	sig = parmDict[sigName]
1820	dsdU = dsdV = dsdW = 0
1821	else:
1822	sig = G2mth.getCWsig(parmDict,tth)
1823	dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1824	sig = max(sig,0.001) #avoid neg sigma
1825	gamName = 'gam'+str(iPeak)
1826	if gamName in varyList or not peakInstPrmMode:
1827	gam = parmDict[gamName]
1828	dgdX = dgdY = dgdZ = 0
1829	else:
1830	gam = G2mth.getCWgam(parmDict,tth)
1831	dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1832	gam = max(gam,0.001) #avoid neg gamma
1833	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig/1.e4,gam/100.)
1834	iBeg = np.searchsorted(xdata,pos-fmin)
1835	iFin = np.searchsorted(xdata,pos+fmin)
1836	if not iBeg+iFin: #peak below low limit
1837	iPeak += 1
1838	continue
1839	elif not iBeg-iFin: #peak above high limit
1840	break
1841	dMdpk = np.zeros(shape=(7,len(xdata)))
1842	dMdipk = getdEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
1843	for i in range(1,6):
1844	dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1845	dMdpk[0][iBeg:iFin] += dMdipk[0]
1846	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4]/1.e4,'gam':dMdpk[5]/100.}
1847	for parmName in ['pos','int','alp','bet','sig','gam']:
1848	try:
1849	idx = varyList.index(parmName+str(iPeak))
1850	dMdv[idx] = dervDict[parmName]
1851	except ValueError:
1852	pass
1853	if 'U' in varyList:
1854	dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1855	if 'V' in varyList:
1856	dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1857	if 'W' in varyList:
1858	dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1859	if 'X' in varyList:
1860	dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1861	if 'Y' in varyList:
1862	dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1863	if 'Z' in varyList:
1864	dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1865	if 'alpha-0' in varyList:
1866	dMdv[varyList.index('alpha-0')] += dada0*dervDict['alp']
1867	if 'alpha-1' in varyList:
1868	dMdv[varyList.index('alpha-1')] += dada1*dervDict['alp']
1869	if 'beta-0' in varyList:
1870	dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1871	if 'beta-1' in varyList:
1872	dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1873	iPeak += 1
1874	except KeyError: #no more peaks to process
1875	break
1876	else:
1877	Pdabc = parmDict['Pdabc']
1878	difC = parmDict['difC']
1879	iPeak = 0
1880	while True:
1881	try:
1882	pos = parmDict['pos'+str(iPeak)]
1883	tof = pos-parmDict['Zero']
1884	dsp = tof/difC
1885	intens = parmDict['int'+str(iPeak)]
1886	alpName = 'alp'+str(iPeak)
1887	if alpName in varyList or not peakInstPrmMode:
1888	alp = parmDict[alpName]
1889	else:
1890	if len(Pdabc):
1891	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1892	dada0 = 0
1893	else:
1894	alp = G2mth.getTOFalpha(parmDict,dsp)
1895	dada0 = G2mth.getTOFalphaDeriv(dsp)
1896	betName = 'bet'+str(iPeak)
1897	if betName in varyList or not peakInstPrmMode:
1898	bet = parmDict[betName]
1899	else:
1900	if len(Pdabc):
1901	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1902	dbdb0 = dbdb1 = dbdb2 = 0
1903	else:
1904	bet = G2mth.getTOFbeta(parmDict,dsp)
1905	dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1906	sigName = 'sig'+str(iPeak)
1907	if sigName in varyList or not peakInstPrmMode:
1908	sig = parmDict[sigName]
1909	dsds0 = dsds1 = dsds2 = dsds3 = 0
1910	else:
1911	sig = G2mth.getTOFsig(parmDict,dsp)
1912	dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1913	gamName = 'gam'+str(iPeak)
1914	if gamName in varyList or not peakInstPrmMode:
1915	gam = parmDict[gamName]
1916	dsdX = dsdY = dsdZ = 0
1917	else:
1918	gam = G2mth.getTOFgamma(parmDict,dsp)
1919	dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1920	gam = max(gam,0.001) #avoid neg gamma
1921	Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1922	iBeg = np.searchsorted(xdata,pos-fmin)
1923	lenX = len(xdata)
1924	if not iBeg:
1925	iFin = np.searchsorted(xdata,pos+fmax)
1926	elif iBeg == lenX:
1927	iFin = iBeg
1928	else:
1929	iFin = np.searchsorted(xdata,pos+fmax)
1930	if not iBeg+iFin: #peak below low limit
1931	iPeak += 1
1932	continue
1933	elif not iBeg-iFin: #peak above high limit
1934	break
1935	dMdpk = np.zeros(shape=(7,len(xdata)))
1936	dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1937	for i in range(1,6):
1938	dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1939	dMdpk[0][iBeg:iFin] += dMdipk[0]
1940	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1941	for parmName in ['pos','int','alp','bet','sig','gam']:
1942	try:
1943	idx = varyList.index(parmName+str(iPeak))
1944	dMdv[idx] = dervDict[parmName]
1945	except ValueError:
1946	pass
1947	if 'alpha' in varyList:
1948	dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1949	if 'beta-0' in varyList:
1950	dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1951	if 'beta-1' in varyList:
1952	dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1953	if 'beta-q' in varyList:
1954	dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1955	if 'sig-0' in varyList:
1956	dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1957	if 'sig-1' in varyList:
1958	dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1959	if 'sig-2' in varyList:
1960	dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1961	if 'sig-q' in varyList:
1962	dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1963	if 'X' in varyList:
1964	dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1965	if 'Y' in varyList:
1966	dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1967	if 'Z' in varyList:
1968	dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1969	iPeak += 1
1970	except KeyError: #no more peaks to process
1971	break
1972	if 'BF mult' in varyList:
1973	dMdv[varyList.index('BF mult')] = fixback
1974
1975	return dMdv
1976
1977	def Dict2Values(parmdict, varylist):
1978	'''Use before call to leastsq to setup list of values for the parameters
1979	in parmdict, as selected by key in varylist'''
1980	return [parmdict[key] for key in varylist]
1981
1982	def Values2Dict(parmdict, varylist, values):
1983	''' Use after call to leastsq to update the parameter dictionary with
1984	values corresponding to keys in varylist'''
1985	parmdict.update(zip(varylist,values))
1986
1987	def SetBackgroundParms(Background):
1988	'Loads background parameters into dicts/lists to create varylist & parmdict'
1989	if len(Background) == 1: # fix up old backgrounds
1990	Background.append({'nDebye':0,'debyeTerms':[]})
1991	bakType,bakFlag = Background[0][:2]
1992	backVals = Background[0][3:]
1993	backNames = ['Back;'+str(i) for i in range(len(backVals))]
1994	Debye = Background[1] #also has background peaks stuff
1995	backDict = dict(zip(backNames,backVals))
1996	backVary = []
1997	if bakFlag:
1998	backVary = backNames
1999
2000	backDict['nDebye'] = Debye['nDebye']
2001	debyeDict = {}
2002	debyeList = []
2003	for i in range(Debye['nDebye']):
2004	debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
2005	debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
2006	debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
2007	debyeVary = []
2008	for item in debyeList:
2009	if item[1]:
2010	debyeVary.append(item[0])
2011	backDict.update(debyeDict)
2012	backVary += debyeVary
2013
2014	backDict['nPeaks'] = Debye['nPeaks']
2015	peaksDict = {}
2016	peaksList = []
2017	for i in range(Debye['nPeaks']):
2018	peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
2019	peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
2020	peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
2021	peaksVary = []
2022	for item in peaksList:
2023	if item[1]:
2024	peaksVary.append(item[0])
2025	backDict.update(peaksDict)
2026	backVary += peaksVary
2027	if 'background PWDR' in Background[1]:
2028	backDict['Back File'] = Background[1]['background PWDR'][0]
2029	backDict['BF mult'] = Background[1]['background PWDR'][1]
2030	if len(Background[1]['background PWDR']) > 2:
2031	if Background[1]['background PWDR'][2]:
2032	backVary += ['BF mult',]
2033	return bakType,backDict,backVary
2034
2035	def DoCalibInst(IndexPeaks,Inst):
2036
2037	def SetInstParms():
2038	dataType = Inst['Type'][0]
2039	insVary = []
2040	insNames = []
2041	insVals = []
2042	for parm in Inst:
2043	insNames.append(parm)
2044	insVals.append(Inst[parm][1])
2045	if parm in ['Lam','difC','difA','difB','Zero','2-theta','XE','YE','ZE']:
2046	if Inst[parm][2]:
2047	insVary.append(parm)
2048	instDict = dict(zip(insNames,insVals))
2049	return dataType,instDict,insVary
2050
2051	def GetInstParms(parmDict,Inst,varyList):
2052	for name in Inst:
2053	Inst[name][1] = parmDict[name]
2054
2055	def InstPrint(Inst,sigDict):
2056	print ('Instrument Parameters:')
2057	if 'C' in Inst['Type'][0] or 'B' in Inst['Type'][0]:
2058	ptfmt = "%12.6f"
2059	else:
2060	ptfmt = "%12.3f"
2061	ptlbls = 'names :'
2062	ptstr = 'values:'
2063	sigstr = 'esds :'
2064	for parm in Inst:
2065	if parm in ['Lam','difC','difA','difB','Zero','2-theta','XE','YE','ZE']:
2066	ptlbls += "%s" % (parm.center(12))
2067	ptstr += ptfmt % (Inst[parm][1])
2068	if parm in sigDict:
2069	sigstr += ptfmt % (sigDict[parm])
2070	else:
2071	sigstr += 12*' '
2072	print (ptlbls)
2073	print (ptstr)
2074	print (sigstr)
2075
2076	def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
2077	parmDict.update(zip(varyList,values))
2078	return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
2079
2080	peakPos = []
2081	peakDsp = []
2082	peakWt = []
2083	for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
2084	if peak[2] and peak[3] and sig > 0.:
2085	peakPos.append(peak[0])
2086	peakDsp.append(peak[-1]) #d-calc
2087	# peakWt.append(peak[-1]2/sig2) #weight by d**2
2088	peakWt.append(1./(sig*peak[-1])) #
2089	peakPos = np.array(peakPos)
2090	peakDsp = np.array(peakDsp)
2091	peakWt = np.array(peakWt)
2092	dataType,insDict,insVary = SetInstParms()
2093	parmDict = {}
2094	parmDict.update(insDict)
2095	varyList = insVary
2096	if not len(varyList):
2097	G2fil.G2Print ('** ERROR - nothing to refine! **')
2098	return False
2099	while True:
2100	begin = time.time()
2101	values = np.array(Dict2Values(parmDict, varyList))
2102	result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
2103	args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
2104	ncyc = int(result[2]['nfev']/2)
2105	runtime = time.time()-begin
2106	chisq = np.sum(result[2]['fvec']**2)
2107	Values2Dict(parmDict, varyList, result[0])
2108	GOF = chisq/(len(peakPos)-len(varyList)) #reduced chi^2
2109	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
2110	G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2111	G2fil.G2Print ('chi2 = %12.6g, reduced chi2 = %6.2f'%(chisq,GOF))
2112	try:
2113	sig = np.sqrt(np.diag(result[1])*GOF)
2114	if np.any(np.isnan(sig)):
2115	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *')
2116	break #refinement succeeded - finish up!
2117	except ValueError: #result[1] is None on singular matrix
2118	G2fil.G2Print ('** Refinement failed - singular matrix **')
2119
2120	sigDict = dict(zip(varyList,sig))
2121	GetInstParms(parmDict,Inst,varyList)
2122	InstPrint(Inst,sigDict)
2123	return True
2124
2125	def getHeaderInfo(dataType):
2126	'''Provide parameter name, label name and formatting information for the
2127	contents of the peak table and where used in DoPeakFit
2128	'''
2129	names = ['pos','int']
2130	lnames = ['position','intensity']
2131	if 'LF' in dataType:
2132	names = ['int','sig','gam','damp','asym','l']
2133	lnames = ['intensity','sigma','gamma','satellite\ndamping','satellite\nasym','00l']
2134	fmt = ["%10.2f","%10.3f","%10.3f","%10.3f","%10.3f","%4d"]
2135	elif 'C' in dataType:
2136	names += ['sig','gam']
2137	lnames += ['sigma','gamma']
2138	fmt = ["%10.5f","%10.1f","%10.3f","%10.3f"]
2139	elif 'T' in dataType:
2140	names += ['alp','bet','sig','gam']
2141	lnames += ['alpha','beta','sigma','gamma']
2142	fmt = ["%10.2f","%10.4f","%8.3f","%8.5f","%10.3f","%10.3f"]
2143	elif 'E' in dataType:
2144	names += ['sig']
2145	lnames += ['sigma']
2146	fmt = ["%10.5f","%10.1f","%8.3f"]
2147	else: # 'B'
2148	names += ['alp','bet','sig','gam']
2149	lnames += ['alpha','beta','sigma','gamma']
2150	fmt = ["%10.5f","%10.1f","%8.2f","%8.4f","%10.3f","%10.3f"]
2151	return names, fmt, lnames
2152
2153	def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],
2154	oneCycle=False,controls=None,wtFactor=1.0,dlg=None,noFit=False):
2155	'''Called to perform a peak fit, refining the selected items in the peak
2156	table as well as selected items in the background.
2157
2158	:param str FitPgm: type of fit to perform. At present this is ignored.
2159	:param list Peaks: a list of peaks. Each peak entry is a list with paired values:
2160	The number of pairs depends on the data type (see :func:`getHeaderInfo`).
2161	For CW data there are
2162	four values each followed by a refine flag where the values are: position, intensity,
2163	sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
2164	tree entry, dict item "peaks". For some types of fits, overall parameters are placed
2165	in a dict entry.
2166	:param list Background: describes the background. List with two items.
2167	Item 0 specifies a background model and coefficients. Item 1 is a dict.
2168	From the Histogram/Background tree entry.
2169	:param list Limits: min and max x-value to use
2170	:param dict Inst: Instrument parameters
2171	:param dict Inst2: more Instrument parameters
2172	:param numpy.array data: a 5xn array. data[0] is the x-values,
2173	data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
2174	calc, background and difference intensities, respectively.
2175	:param array fixback: fixed background array; same size as data[0-5]
2176	:param list prevVaryList: Used in sequential refinements to override the
2177	variable list. Defaults as an empty list.
2178	:param bool oneCycle: True if only one cycle of fitting should be performed
2179	:param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
2180	and derivType = controls['deriv type']. If None default values are used.
2181	:param float wtFactor: weight multiplier; = 1.0 by default
2182	:param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
2183	Defaults to None, which means no updates are done.
2184	:param bool noFit: When noFit is True, a refinement is not performed. Default
2185	is False.
2186
2187	'''
2188	def GetBackgroundParms(parmList,Background):
2189	iBak = 0
2190	while True:
2191	try:
2192	bakName = 'Back;'+str(iBak)
2193	Background[0][iBak+3] = parmList[bakName]
2194	iBak += 1
2195	except KeyError:
2196	break
2197	iDb = 0
2198	while True:
2199	names = ['DebyeA;','DebyeR;','DebyeU;']
2200	try:
2201	for i,name in enumerate(names):
2202	val = parmList[name+str(iDb)]
2203	Background[1]['debyeTerms'][iDb][2*i] = val
2204	iDb += 1
2205	except KeyError:
2206	break
2207	iDb = 0
2208	while True:
2209	names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
2210	try:
2211	for i,name in enumerate(names):
2212	val = parmList[name+str(iDb)]
2213	Background[1]['peaksList'][iDb][2*i] = val
2214	iDb += 1
2215	except KeyError:
2216	break
2217	if 'BF mult' in parmList:
2218	Background[1]['background PWDR'][1] = parmList['BF mult']
2219
2220	def BackgroundPrint(Background,sigDict):
2221	print ('Background coefficients for '+Background[0][0]+' function')
2222	ptfmt = "%12.5f"
2223	ptstr = 'value: '
2224	sigstr = 'esd : '
2225	for i,back in enumerate(Background[0][3:]):
2226	ptstr += ptfmt % (back)
2227	if Background[0][1]:
2228	prm = 'Back;'+str(i)
2229	if prm in sigDict:
2230	sigstr += ptfmt % (sigDict[prm])
2231	else:
2232	sigstr += " "*12
2233	if len(ptstr) > 75:
2234	print (ptstr)
2235	if Background[0][1]: print (sigstr)
2236	ptstr = 'value: '
2237	sigstr = 'esd : '
2238	if len(ptstr) > 8:
2239	print (ptstr)
2240	if Background[0][1]: print (sigstr)
2241
2242	if Background[1]['nDebye']:
2243	parms = ['DebyeA;','DebyeR;','DebyeU;']
2244	print ('Debye diffuse scattering coefficients')
2245	ptfmt = "%12.5f"
2246	print (' term DebyeA esd DebyeR esd DebyeU esd')
2247	for term in range(Background[1]['nDebye']):
2248	line = ' term %d'%(term)
2249	for ip,name in enumerate(parms):
2250	line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
2251	if name+str(term) in sigDict:
2252	line += ptfmt%(sigDict[name+str(term)])
2253	else:
2254	line += " "*12
2255	print (line)
2256	if Background[1]['nPeaks']:
2257	print ('Coefficients for Background Peaks')
2258	ptfmt = "%15.3f"
2259	for j,pl in enumerate(Background[1]['peaksList']):
2260	names = 'peak %3d:'%(j+1)
2261	ptstr = 'values :'
2262	sigstr = 'esds :'
2263	for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
2264	val = pl[2*i]
2265	prm = lbl+";"+str(j)
2266	names += '%15s'%(prm)
2267	ptstr += ptfmt%(val)
2268	if prm in sigDict:
2269	sigstr += ptfmt%(sigDict[prm])
2270	else:
2271	sigstr += " "*15
2272	print (names)
2273	print (ptstr)
2274	print (sigstr)
2275	if 'BF mult' in sigDict:
2276	print('Background file mult: %.3f(%d)'%(Background[1]['background PWDR'][1],int(1000*sigDict['BF mult'])))
2277
2278	def SetInstParms(Inst):
2279	dataType = Inst['Type'][0]
2280	insVary = []
2281	insNames = []
2282	insVals = []
2283	for parm in Inst:
2284	insNames.append(parm)
2285	insVals.append(Inst[parm][1])
2286	if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha','A','B','C',
2287	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1'] and Inst[parm][2]:
2288	insVary.append(parm)
2289	instDict = dict(zip(insNames,insVals))
2290	if 'SH/L' in instDict:
2291	instDict['SH/L'] = max(instDict['SH/L'],0.002)
2292	return dataType,instDict,insVary
2293
2294	def GetPkInstParms(parmDict,Inst,varyList):
2295	for name in Inst:
2296	Inst[name][1] = parmDict[name]
2297	iPeak = 0
2298	while True:
2299	try:
2300	sigName = 'sig'+str(iPeak)
2301	pos = parmDict['pos'+str(iPeak)]
2302	if sigName not in varyList and peakInstPrmMode:
2303	if 'T' in Inst['Type'][0]:
2304	dsp = G2lat.Pos2dsp(Inst,pos)
2305	parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
2306	if 'E' in Inst['Type'][0]:
2307	parmDict[sigName] = G2mth.getEDsig(parmDict,pos)
2308	else:
2309	parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
2310	gamName = 'gam'+str(iPeak)
2311	if gamName not in varyList and peakInstPrmMode:
2312	if 'T' in Inst['Type'][0]:
2313	dsp = G2lat.Pos2dsp(Inst,pos)
2314	parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
2315	else:
2316	parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
2317	iPeak += 1
2318	except KeyError:
2319	break
2320
2321	def InstPrint(Inst,sigDict):
2322	print ('Instrument Parameters:')
2323	ptfmt = "%12.6f"
2324	ptlbls = 'names :'
2325	ptstr = 'values:'
2326	sigstr = 'esds :'
2327	for parm in Inst:
2328	if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha','A','B','C',
2329	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1']:
2330	ptlbls += "%s" % (parm.center(12))
2331	ptstr += ptfmt % (Inst[parm][1])
2332	if parm in sigDict:
2333	sigstr += ptfmt % (sigDict[parm])
2334	else:
2335	sigstr += 12*' '
2336	print (ptlbls)
2337	print (ptstr)
2338	print (sigstr)
2339
2340	def SetPeaksParms(dataType,Peaks):
2341	'''Set the contents of peakDict from list Peaks
2342	'''
2343	peakDict = {}
2344	peakVary = []
2345	names,_,_ = getHeaderInfo(dataType)
2346	off = 0
2347	if 'LF' in dataType: off = 2
2348	for i,peak in enumerate(Peaks):
2349	if type(peak) is dict:
2350	peakDict.update(peak)
2351	continue
2352	for j,name in enumerate(names):
2353	parName = name+str(i)
2354	peakDict[parName] = peak[off+2*j]
2355	if peak[off+2*j+1]:
2356	peakVary.append(parName)
2357	return peakDict,peakVary
2358
2359	def GetPeaksParms(Inst,parmDict,Peaks,varyList):
2360	'''Put values into the Peaks list from the refinement results from inside
2361	the parmDict array
2362	'''
2363	off = 0
2364	if 'LF' in Inst['Type'][0]:
2365	off = 2
2366	if 'clat' in varyList:
2367	Peaks[-1]['clat'] = parmDict['clat']
2368	names,_,_ = getHeaderInfo(Inst['Type'][0])
2369	for i,peak in enumerate(Peaks):
2370	if type(peak) is dict:
2371	continue
2372	for j in range(len(names)):
2373	parName = names[j]+str(i)
2374	if parName in varyList or not peakInstPrmMode:
2375	peak[2*j+off] = parmDict[parName]
2376	if 'pos'+str(i) not in parmDict: continue
2377	pos = parmDict['pos'+str(i)]
2378	if 'LF' in Inst['Type'][0]: peak[0] = pos
2379	if 'difC' in Inst:
2380	dsp = pos/Inst['difC'][1]
2381	for j in range(len(names)):
2382	parName = names[j]+str(i)
2383	if peak[2*j+off + 1] or not peakInstPrmMode: continue
2384	if 'alp' in parName:
2385	if 'T' in Inst['Type'][0]:
2386	peak[2*j+off] = G2mth.getTOFalpha(parmDict,dsp)
2387	else: #'B'
2388	peak[2*j+off] = G2mth.getPinkalpha(parmDict,pos)
2389	elif 'bet' in parName:
2390	if 'T' in Inst['Type'][0]:
2391	peak[2*j+off] = G2mth.getTOFbeta(parmDict,dsp)
2392	else: #'B'
2393	peak[2*j+off] = G2mth.getPinkbeta(parmDict,pos)
2394	elif 'sig' in parName:
2395	if 'T' in Inst['Type'][0]:
2396	peak[2*j+off] = G2mth.getTOFsig(parmDict,dsp)
2397	elif 'E' in Inst['Type'][0]:
2398	peak[2*j+off] = G2mth.getEDsig(parmDict,pos)
2399	else: #'C' & 'B'
2400	peak[2*j+off] = G2mth.getCWsig(parmDict,pos)
2401	elif 'gam' in parName:
2402	if 'T' in Inst['Type'][0]:
2403	peak[2*j+off] = G2mth.getTOFgamma(parmDict,dsp)
2404	else: #'C' & 'B'
2405	peak[2*j+off] = G2mth.getCWgam(parmDict,pos)
2406
2407	def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
2408	if 'clat' in varyList:
2409	print('c = {:.6f} esd {:.6f}'.format(parmDict['clat'],sigDict['clat']))
2410	print ('Peak coefficients:')
2411	names,fmt,_ = getHeaderInfo(dataType)
2412	head = 13*' '
2413	for name in names:
2414	if name in ['alp','bet']:
2415	head += name.center(8)+'esd'.center(8)
2416	else:
2417	head += name.center(10)+'esd'.center(10)
2418	head += 'bins'.center(8)
2419	print (head)
2420	ptfmt = dict(zip(names,fmt))
2421	for i,peak in enumerate(Peaks):
2422	if type(peak) is dict:
2423	continue
2424	ptstr = ':'
2425	for j in range(len(names)):
2426	name = names[j]
2427	parName = name+str(i)
2428	ptstr += ptfmt[name] % (parmDict[parName])
2429	if parName in varyList:
2430	ptstr += ptfmt[name] % (sigDict[parName])
2431	else:
2432	if name in ['alp','bet']:
2433	ptstr += 8*' '
2434	else:
2435	ptstr += 10*' '
2436	ptstr += '%9.2f'%(ptsperFW[i])
2437	print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
2438
2439	def devPeakProfile(values,xdata,ydata,fixback, weights,dataType,parmdict,varylist,bakType,dlg):
2440	'''Computes a matrix where each row is the derivative of the calc-obs
2441	values (see :func:`errPeakProfile`) with respect to each parameter
2442	in backVary,insVary,peakVary. Used for peak fitting.
2443	'''
2444	parmdict.update(zip(varylist,values))
2445	return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,fixback,varylist,bakType)
2446
2447	def errPeakProfile(values,xdata,ydata,fixback,weights,dataType,parmdict,varylist,bakType,dlg):
2448	'''Computes a vector with the weighted calc-obs values differences
2449	for peak fitting
2450	'''
2451	parmdict.update(zip(varylist,values))
2452	M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,fixback,varylist,bakType)-ydata)
2453	Rwp = min(100.,np.sqrt(np.sum(M*2)/np.sum(weightsydata*2))100.)
2454	if dlg:
2455	dlg.Raise()
2456	GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
2457	if not GoOn:
2458	return -M #abort!!
2459	return M
2460
2461	#---- beginning of DoPeakFit ---------------------------------------------
2462	if controls:
2463	Ftol = controls['min dM/M']
2464	else:
2465	Ftol = 0.0001
2466	if oneCycle:
2467	Ftol = 1.0
2468	x,y,w,yc,yb,yd = data #these are numpy arrays - remove masks!
2469	if fixback is None:
2470	fixback = np.zeros_like(y)
2471	yc.fill(0.) #set calcd ones to zero
2472	yb.fill(0.)
2473	yd.fill(0.)
2474	xBeg = np.searchsorted(x,Limits[0])
2475	xFin = np.searchsorted(x,Limits[1])+1
2476	# find out what is varied
2477	bakType,bakDict,bakVary = SetBackgroundParms(Background)
2478	dataType,insDict,insVary = SetInstParms(Inst)
2479	peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
2480	parmDict = {}
2481	parmDict.update(bakDict)
2482	parmDict.update(insDict)
2483	parmDict.update(peakDict)
2484	parmDict['Pdabc'] = [] #dummy Pdabc
2485	parmDict.update(Inst2) #put in real one if there
2486	if prevVaryList:
2487	varyList = prevVaryList[:]
2488	else:
2489	varyList = bakVary+insVary+peakVary
2490	if 'LF' in Inst['Type'][0] and Peaks:
2491	if Peaks[-1].get('clat-ref'): varyList += ['clat']
2492	fullvaryList = varyList[:]
2493	if not peakInstPrmMode:
2494	for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1','A','B','C',
2495	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',):
2496	if v in varyList:
2497	raise Exception('Instrumental profile terms cannot be varied '+
2498	'after setPeakInstPrmMode(False) is used')
2499	if 'LF' in Inst['Type'][0]:
2500	warn = []
2501	for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1',
2502	'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',):
2503	if v in varyList:
2504	warn.append(v)
2505	del varyList[varyList.index(v)]
2506	if warn:
2507	print('Instrumental profile terms cannot be varied '+
2508	'in Laue Fringe fits:',warn)
2509
2510	while not noFit:
2511	begin = time.time()
2512	values = np.array(Dict2Values(parmDict, varyList))
2513	Rvals = {}
2514	badVary = []
2515	try:
2516	result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
2517	args=(x[xBeg:xFin],y[xBeg:xFin],fixback[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
2518	except Exception as msg:
2519	if GSASIIpath.GetConfigValue('debug'):
2520	print('peak fit failure\n',msg)
2521	import traceback
2522	print (traceback.format_exc())
2523	else:
2524	print('peak fit failure')
2525	return
2526	ncyc = int(result[2]['nfev']/2)
2527	runtime = time.time()-begin
2528	chisq = np.sum(result[2]['fvec']**2)
2529	Values2Dict(parmDict, varyList, result[0])
2530	Rvals['Rwp'] = np.sqrt(chisq/np.sum(wtFactorw[xBeg:xFin]y[xBeg:xFin]*2))100. #to %
2531	Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList)) #reduced chi^2
2532	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
2533	if ncyc:
2534	G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2535	G2fil.G2Print ('Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2536	sig = [0]*len(varyList)
2537	if len(varyList) == 0: break # if nothing was refined
2538	try:
2539	sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
2540	if np.any(np.isnan(sig)):
2541	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *')
2542	break #refinement succeeded - finish up!
2543	except ValueError: #result[1] is None on singular matrix
2544	G2fil.G2Print ('** Refinement failed - singular matrix **')
2545	Ipvt = result[2]['ipvt']
2546	for i,ipvt in enumerate(Ipvt):
2547	if not np.sum(result[2]['fjac'],axis=1)[i]:
2548	G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
2549	badVary.append(varyList[ipvt-1])
2550	del(varyList[ipvt-1])
2551	break
2552	else: # nothing removed
2553	break
2554	if dlg: dlg.Destroy()
2555	yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin],fixback[xBeg:xFin])[0]
2556	yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],fixback[xBeg:xFin],varyList,bakType)
2557	yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
2558	if noFit:
2559	GetPeaksParms(Inst,parmDict,Peaks,varyList)
2560	return
2561	sigDict = dict(zip(varyList,sig))
2562	GetBackgroundParms(parmDict,Background)
2563	if bakVary: BackgroundPrint(Background,sigDict)
2564	GetPkInstParms(parmDict,Inst,varyList)
2565	if insVary: InstPrint(Inst,sigDict)
2566	GetPeaksParms(Inst,parmDict,Peaks,varyList)
2567	binsperFWHM = []
2568	for peak in Peaks:
2569	if type(peak) is dict:
2570	continue
2571	FWHM = getFWHM(peak[0],Inst)
2572	try:
2573	xpk = x.searchsorted(peak[0])
2574	cw = x[xpk]-x[xpk-1]
2575	binsperFWHM.append(FWHM/cw)
2576	except IndexError:
2577	binsperFWHM.append(0.)
2578	if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
2579	if len(binsperFWHM):
2580	if min(binsperFWHM) < 1.:
2581	G2fil.G2Print ('* Warning: calculated peak widths are too narrow to refine profile coefficients *')
2582	if 'T' in Inst['Type'][0]:
2583	G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
2584	else:
2585	G2fil.G2Print (' Manually increase W in Instrument Parameters')
2586	elif min(binsperFWHM) < 4.:
2587	G2fil.G2Print ('* Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement *')
2588	G2fil.G2Print ('* recommended is 6-10; you have %.2f *'%(min(binsperFWHM)))
2589	return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
2590
2591	def calcIncident(Iparm,xdata):
2592	'needs a doc string'
2593
2594	def IfunAdv(Iparm,xdata):
2595	Itype = Iparm['Itype']
2596	Icoef = Iparm['Icoeff']
2597	DYI = np.ones((12,xdata.shape[0]))
2598	YI = np.ones_like(xdata)*Icoef[0]
2599
2600	x = xdata/1000. #expressions are in ms
2601	if Itype == 'Exponential':
2602	for i in [1,3,5,7,9]:
2603	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
2604	YI += Icoef[i]*Eterm
2605	DYI[i] *= Eterm
2606	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
2607	elif 'Maxwell'in Itype:
2608	Eterm = np.exp(-Icoef[2]/x**2)
2609	DYI[1] = Eterm/x**5
2610	DYI[2] = -Icoef[1]DYI[1]/x*2
2611	YI += (Icoef[1]Eterm/x*5)
2612	if 'Exponential' in Itype:
2613	for i in range(3,11,2):
2614	Eterm = np.exp(-Icoef[i+1]x*((i+1)/2))
2615	YI += Icoef[i]*Eterm
2616	DYI[i] *= Eterm
2617	DYI[i+1] = -Icoef[i]Etermx*((i+1)/2)
2618	else: #Chebyschev
2619	T = (2./x)-1.
2620	Ccof = np.ones((12,xdata.shape[0]))
2621	Ccof[1] = T
2622	for i in range(2,12):
2623	Ccof[i] = 2TCcof[i-1]-Ccof[i-2]
2624	for i in range(1,10):
2625	YI += Ccof[i]*Icoef[i+2]
2626	DYI[i+2] =Ccof[i]
2627	return YI,DYI
2628
2629	Iesd = np.array(Iparm['Iesd'])
2630	Icovar = Iparm['Icovar']
2631	YI,DYI = IfunAdv(Iparm,xdata)
2632	YI = np.where(YI>0,YI,1.)
2633	WYI = np.zeros_like(xdata)
2634	vcov = np.zeros((12,12))
2635	k = 0
2636	for i in range(12):
2637	for j in range(i,12):
2638	vcov[i][j] = Icovar[k]Iesd[i]Iesd[j]
2639	vcov[j][i] = Icovar[k]Iesd[i]Iesd[j]
2640	k += 1
2641	M = np.inner(vcov,DYI.T)
2642	WYI = np.sum(M*DYI,axis=0)
2643	WYI = np.where(WYI>0.,WYI,0.)
2644	return YI,WYI
2645
2646	#### RMCutilities ################################################################################
2647	def MakeInst(PWDdata,Name,Size,Mustrain,useSamBrd):
2648	inst = PWDdata['Instrument Parameters'][0]
2649	Xsb = 0.
2650	Ysb = 0.
2651	if 'T' in inst['Type'][1]:
2652	difC = inst['difC'][1]
2653	if useSamBrd[0]:
2654	if 'ellipsoidal' not in Size[0]: #take the isotropic term only
2655	Xsb = 1.e-4*difC/Size[1][0]
2656	if useSamBrd[1]:
2657	if 'generalized' not in Mustrain[0]: #take the isotropic term only
2658	Ysb = 1.e-6difCMustrain[1][0]
2659	prms = ['Bank',
2660	'difC','difA','Zero','2-theta',
2661	'alpha','beta-0','beta-1',
2662	'sig-0','sig-1','sig-2',
2663	'Z','X','Y']
2664	fname = Name+'.inst'
2665	fl = open(fname,'w')
2666	fl.write('1\n')
2667	fl.write('%d\n'%int(inst[prms[0]][1]))
2668	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]))
2669	fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[5]][1],inst[prms[6]][1],inst[prms[7]][1]))
2670	fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[8]][1],inst[prms[9]][1],inst[prms[10]][1]))
2671	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))
2672	fl.write('\n\n\n')
2673	fl.close()
2674	else:
2675	if useSamBrd[0]:
2676	wave = G2mth.getWave(inst)
2677	if 'ellipsoidal' not in Size[0]: #take the isotropic term only
2678	Xsb = 1.8wave/(np.piSize[1][0])
2679	if useSamBrd[1]:
2680	if 'generalized' not in Mustrain[0]: #take the isotropic term only
2681	Ysb = 0.0180*Mustrain[1][0]/np.pi
2682	prms = ['Bank',
2683	'Lam','Zero','Polariz.',
2684	'U','V','W',
2685	'X','Y']
2686	fname = Name+'.inst'
2687	fl = open(fname,'w')
2688	fl.write('1\n')
2689	fl.write('%d\n'%int(inst[prms[0]][1]))
2690	fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2691	fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2692	fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))
2693	fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2694	fl.write('\n\n\n')
2695	fl.close()
2696	return fname
2697
2698	def MakeBack(PWDdata,Name):
2699	Back = PWDdata['Background'][0]
2700	inst = PWDdata['Instrument Parameters'][0]
2701	if 'chebyschev-1' != Back[0]:
2702	return None
2703	Nback = Back[2]
2704	BackVals = Back[3:]
2705	fname = Name+'.back'
2706	fl = open(fname,'w')
2707	fl.write('%10d\n'%Nback)
2708	for val in BackVals:
2709	if 'T' in inst['Type'][1]:
2710	fl.write('%12.6g\n'%(float(val)))
2711	else:
2712	fl.write('%12.6g\n'%val)
2713	fl.close()
2714	return fname
2715
2716	def findDup(Atoms):
2717	Dup = []
2718	Fracs = []
2719	for iat1,at1 in enumerate(Atoms):
2720	if any([at1[0] in dup for dup in Dup]):
2721	continue
2722	else:
2723	Dup.append([at1[0],])
2724	Fracs.append([at1[6],])
2725	for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2726	if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2727	Dup[-1] += [at2[0],]
2728	Fracs[-1]+= [at2[6],]
2729	return Dup,Fracs
2730
2731	def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):
2732
2733	Meta = RMCPdict['metadata']
2734	Atseq = RMCPdict['atSeq']
2735	Supercell = RMCPdict['SuperCell']
2736	generalData = Phase['General']
2737	Dups,Fracs = findDup(Phase['Atoms'])
2738	Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2739	ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
2740	Sample = PWDdata['Sample Parameters']
2741	Meta['temperature'] = Sample['Temperature']
2742	Meta['pressure'] = Sample['Pressure']
2743	Cell = generalData['Cell'][1:7]
2744	Trans = np.eye(3)*np.array(Supercell)
2745	newPhase = copy.deepcopy(Phase)
2746	newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2747	newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2748	GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2749	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.
2750	newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
2751	Natm = np.core.defchararray.count(np.array(Atcodes),'+') #no. atoms in original unit cell
2752	Natm = np.count_nonzero(Natm-1)
2753	Atoms = newPhase['Atoms']
2754	reset = False
2755
2756	if ifSfracs:
2757	Natm = np.core.defchararray.count(np.array(Atcodes),'+') #no. atoms in original unit cell
2758	Natm = np.count_nonzero(Natm-1)
2759	Satoms = []
2760	for i in range(len(Atoms)//Natm):
2761	ind = i*Natm
2762	Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
2763	Natoms = []
2764	for satoms in Satoms:
2765	for idup,dup in enumerate(Dups):
2766	ldup = len(dup)
2767	natm = len(satoms)
2768	i = 0
2769	while i < natm:
2770	if satoms[i][0] in dup:
2771	atoms = satoms[i:i+ldup]
2772	try:
2773	atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2774	Natoms.append(atom)
2775	except IndexError: #what about vacancies?
2776	if 'Va' not in Atseq:
2777	reset = True
2778	Atseq.append('Va')
2779	RMCPdict['aTypes']['Va'] = 0.0
2780	atom = atoms[0]
2781	atom[1] = 'Va'
2782	Natoms.append(atom)
2783	i += ldup
2784	else:
2785	i += 1
2786	else:
2787	Natoms = Atoms
2788
2789	NAtype = np.zeros(len(Atseq))
2790	for atom in Natoms:
2791	NAtype[Atseq.index(atom[1])] += 1
2792	NAstr = ['%6d'%i for i in NAtype]
2793	Cell = newPhase['General']['Cell'][1:7]
2794	if os.path.exists(Name+'.his6f'):
2795	os.remove(Name+'.his6f')
2796	if os.path.exists(Name+'.neigh'):
2797	os.remove(Name+'.neigh')
2798	fname = Name+'.rmc6f'
2799	fl = open(fname,'w')
2800	fl.write('(Version 6f format configuration file)\n')
2801	for item in Meta:
2802	fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2803	fl.write('Atom types present: %s\n'%' '.join(Atseq))
2804	fl.write('Number of each atom type: %s\n'%''.join(NAstr))
2805	fl.write('Number of atoms: %d\n'%len(Natoms))
2806	fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2807	fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2808	Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2809	A,B = G2lat.cell2AB(Cell,True)
2810	fl.write('Lattice vectors (Ang):\n')
2811	for i in [0,1,2]:
2812	fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2813	fl.write('Atoms (fractional coordinates):\n')
2814	nat = 0
2815	for atm in Atseq:
2816	for iat,atom in enumerate(Natoms):
2817	if atom[1] == atm:
2818	nat += 1
2819	atcode = Atcodes[iat].split(':')
2820	cell = [0,0,0]
2821	if '+' in atcode[1]:
2822	cell = eval(atcode[1].split('+')[1])
2823	fl.write('%6d%4s [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(
2824	nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2825	fl.close()
2826	return fname,reset
2827
2828	def MakeBragg(PWDdata,Name,Phase):
2829	generalData = Phase['General']
2830	Vol = generalData['Cell'][7]
2831	Data = PWDdata['Data']
2832	Inst = PWDdata['Instrument Parameters'][0]
2833	Bank = int(Inst['Bank'][1])
2834	Sample = PWDdata['Sample Parameters']
2835	Scale = Sample['Scale'][0]
2836	Limits = PWDdata['Limits'][1]
2837	Ibeg = np.searchsorted(Data[0],Limits[0])
2838	Ifin = np.searchsorted(Data[0],Limits[1])+1
2839	fname = Name+'.bragg'
2840	fl = open(fname,'w')
2841	fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2842	if 'T' in Inst['Type'][0]:
2843	fl.write('%12s%12s\n'%(' TOF,ms',' I(obs)'))
2844	for i in range(Ibeg,Ifin-1):
2845	fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2846	else:
2847	fl.write('%12s%12s\n'%(' 2-theta, deg',' I(obs)'))
2848	for i in range(Ibeg,Ifin-1):
2849	fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))
2850	fl.close()
2851	return fname
2852
2853	def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2854	Meta = RMCPdict['metadata']
2855	Times = RMCPdict['runTimes']
2856	Atseq = RMCPdict['atSeq']
2857	Natoms = RMCPdict['NoAtoms']
2858	sumatms = np.sum(np.array([Natoms[iatm] for iatm in Natoms]))
2859	Isotope = RMCPdict['Isotope']
2860	Isotopes = RMCPdict['Isotopes']
2861	Atypes = RMCPdict['aTypes']
2862	if 'Va' in Atypes:
2863	Isotope['Va'] = 'Nat. Abund.'
2864	Isotopes['Va'] = {'Nat. Abund.':{'SL':[0.0,0.0]}}
2865	atPairs = RMCPdict['Pairs']
2866	Files = RMCPdict['files']
2867	BraggWt = RMCPdict['histogram'][1]
2868	inst = PWDdata['Instrument Parameters'][0]
2869	try:
2870	refList = PWDdata['Reflection Lists'][Name]['RefList']
2871	except KeyError:
2872	return 'Error - missing reflection list; you must do Refine first'
2873	dMin = refList[-1][4]
2874	gsasType = 'xray2'
2875	if 'T' in inst['Type'][1]:
2876	gsasType = 'gsas3'
2877	elif 'X' in inst['Type'][1]:
2878	XFF = G2elem.GetFFtable(Atseq)
2879	Xfl = open(Name+'.xray','w')
2880	for atm in Atseq:
2881	fa = XFF[atm]['fa']
2882	fb = XFF[atm]['fb']
2883	fc = XFF[atm]['fc']
2884	Xfl.write('%2s %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2885	atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2886	Xfl.close()
2887	lenA = len(Atseq)
2888	Pairs = []
2889	Ncoeff = []
2890	Nblen = [Isotopes[at][Isotope[at]]['SL'][0] for at in Atypes]
2891	for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2892	Pairs += pair
2893	for pair in Pairs:
2894	pair = pair.replace(' ','')
2895	at1,at2 = pair.split('-')
2896	if at1 == 'Va' or at2 == 'Va':
2897	ncoef = 0.0
2898	else:
2899	ncoef = Isotopes[at1][Isotope[at1]]['SL'][0]*Natoms[at1]/sumatms
2900	ncoef = Isotopes[at2][Isotope[at2]]['SL'][0]Natoms[at2]/sumatms
2901	if at1 != at2:
2902	ncoef *= 2.
2903	Ncoeff += [ncoef,]
2904	pairMin = [atPairs[pair] if pair in atPairs else [0.0,0.,0.] for pair in Pairs ]
2905	maxMoves = [Atypes[atm] if atm in Atypes else 0.0 for atm in Atseq ]
2906	fname = Name+'.dat'
2907	fl = open(fname,'w')
2908	fl.write(' %% Hand edit the following as needed\n')
2909	fl.write('TITLE :: '+Name+'\n')
2910	fl.write('MATERIAL :: '+Meta['material']+'\n')
2911	fl.write('PHASE :: '+Meta['phase']+'\n')
2912	fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2913	fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2914	if RMCPdict.get('useGPU',False):
2915	fl.write('GPU_ACCELERATOR :: 0\n')
2916	minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2917	minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2918	maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2919	fl.write('MINIMUM_DISTANCES :: %s Angstrom\n'%minHD)
2920	maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2921	fl.write('MAXIMUM_MOVES :: %s Angstrom\n'%maxMv)
2922	fl.write('R_SPACING :: 0.0200 Angstrom\n')
2923	fl.write('PRINT_PERIOD :: 100\n')
2924	fl.write('TIME_LIMIT :: %.2f MINUTES\n'%Times[0])
2925	fl.write('SAVE_PERIOD :: %.2f MINUTES\n'%Times[1])
2926	fl.write('\n')
2927	fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2928	fl.write('\n')
2929	fl.write('FLAGS ::\n')
2930	fl.write(' > NO_MOVEOUT\n')
2931	fl.write(' > NO_SAVE_CONFIGURATIONS\n')
2932	fl.write(' > NO_RESOLUTION_CONVOLUTION\n')
2933	fl.write('\n')
2934	fl.write('INPUT_CONFIGURATION_FORMAT :: rmc6f\n')
2935	fl.write('SAVE_CONFIGURATION_FORMAT :: rmc6f\n')
2936	fl.write('IGNORE_HISTORY_FILE ::\n')
2937	fl.write('\n')
2938	if 'T' in inst['Type'][1]:
2939	fl.write('NEUTRON_COEFFICIENTS :: '+''.join(['%9.5f'%coeff for coeff in Ncoeff])+'\n')
2940	fl.write('DISTANCE_WINDOW ::\n')
2941	fl.write(' > MNDIST :: %s\n'%minD)
2942	fl.write(' > MXDIST :: %s\n'%maxD)
2943	if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2944	fl.write('\n')
2945	fl.write('POTENTIALS ::\n')
2946	fl.write(' > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2947	fl.write(' > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2948	if len(RMCPdict['Potentials']['Stretch']):
2949	fl.write(' > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2950	for bond in RMCPdict['Potentials']['Stretch']:
2951	fl.write(' > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))
2952	if len(RMCPdict['Potentials']['Angles']):
2953	fl.write(' > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2954	for angle in RMCPdict['Potentials']['Angles']:
2955	fl.write(' > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2956	(angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2957	if RMCPdict['useBVS']:
2958	fl.write('BVS ::\n')
2959	fl.write(' > ATOM :: '+' '.join(Atseq)+'\n')
2960	fl.write(' > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2961	oxid = []
2962	for val in RMCPdict['Oxid']:
2963	if len(val) == 3:
2964	oxid.append(val[0][1:])
2965	else:
2966	oxid.append(val[0][2:])
2967	fl.write(' > OXID :: %s\n'%' '.join(oxid))
2968	fl.write(' > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2969	fl.write(' > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2970	fl.write(' > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2971	fl.write(' > SAVE :: 100000\n')
2972	fl.write(' > UPDATE :: 100000\n')
2973	if len(RMCPdict['Swap']):
2974	fl.write('\n')
2975	fl.write('SWAP_MULTI ::\n')
2976	for swap in RMCPdict['Swap']:
2977	try:
2978	at1 = Atseq.index(swap[0])
2979	at2 = Atseq.index(swap[1])
2980	except ValueError:
2981	break
2982	fl.write(' > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2983
2984	if len(RMCPdict['FxCN']):
2985	fl.write('FIXED_COORDINATION_CONSTRAINTS :: %d\n'%len(RMCPdict['FxCN']))
2986	for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2987	try:
2988	at1 = Atseq.index(fxcn[0])
2989	at2 = Atseq.index(fxcn[1])
2990	except ValueError:
2991	break
2992	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]))
2993	if len(RMCPdict['AveCN']):
2994	fl.write('AVERAGE_COORDINATION_CONSTRAINTS :: %d\n'%len(RMCPdict['AveCN']))
2995	for iav,avcn in enumerate(RMCPdict['AveCN']):
2996	try:
2997	at1 = Atseq.index(avcn[0])
2998	at2 = Atseq.index(avcn[1])
2999	except ValueError:
3000	break
3001	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]))
3002	for File in Files:
3003	if Files[File][0] and Files[File][0] != 'Select':
3004	if 'Xray' in File and 'F(Q)' in File:
3005	fqdata = open(Files[File][0],'r')
3006	lines = int(fqdata.readline()[:-1])
3007	fqdata.close()
3008	fl.write('\n')
3009	fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
3010	fl.write(' > FILENAME :: %s\n'%Files[File][0])
3011	fl.write(' > DATA_TYPE :: %s\n'%Files[File][2])
3012	fl.write(' > FIT_TYPE :: %s\n'%Files[File][2])
3013	if 'Xray' not in File:
3014	fl.write(' > START_POINT :: 1\n')
3015	fl.write(' > END_POINT :: 3000\n')
3016	fl.write(' > WEIGHT :: %.4f\n'%Files[File][1])
3017	fl.write(' > CONSTANT_OFFSET 0.000\n')
3018	fl.write(' > NO_FITTED_OFFSET\n')
3019	if RMCPdict['FitScale']:
3020	fl.write(' > FITTED_SCALE\n')
3021	else:
3022	fl.write(' > NO_FITTED_SCALE\n')
3023	if Files[File][3] !='RMC':
3024	fl.write(' > %s\n'%Files[File][3])
3025	if 'reciprocal' in File:
3026	fl.write(' > CONVOLVE ::\n')
3027	if 'Xray' in File:
3028	fl.write(' > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
3029	fl.write(' > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
3030	fl.write(' > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
3031	fl.write(' > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
3032	fl.write('\n')
3033	fl.write('BRAGG ::\n')
3034	fl.write(' > BRAGG_SHAPE :: %s\n'%gsasType)
3035	fl.write(' > RECALCUATE\n')
3036	fl.write(' > DMIN :: %.2f\n'%(dMin-0.02))
3037	fl.write(' > WEIGHT :: %10.3f\n'%BraggWt)
3038	if 'T' in inst['Type'][1]:
3039	fl.write(' > SCATTERING LENGTH :: '+''.join(['%8.4f'%blen for blen in Nblen])+'\n')
3040	fl.write('\n')
3041	fl.write('END ::\n')
3042	fl.close()
3043	return fname
3044
3045	# def FindBonds(Phase,RMCPdict):
3046	# generalData = Phase['General']
3047	# cx,ct,cs,cia = generalData['AtomPtrs']
3048	# atomData = Phase['Atoms']
3049	# Res = 'RMC'
3050	# if 'macro' in generalData['Type']:
3051	# Res = atomData[0][ct-3]
3052	# AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
3053	# Pairs = RMCPdict['Pairs'] #dict!
3054	# BondList = []
3055	# notNames = []
3056	# for FrstName in AtDict:
3057	# nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
3058	# Atyp1 = AtDict[FrstName]
3059	# if 'Va' in Atyp1:
3060	# continue
3061	# for nbr in nbrs:
3062	# Atyp2 = AtDict[nbr[0]]
3063	# if 'Va' in Atyp2:
3064	# continue
3065	# try:
3066	# bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
3067	# except KeyError:
3068	# bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
3069	# if any(bndData):
3070	# if bndData[0] <= nbr[1] <= bndData[1]:
3071	# bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
3072	# revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
3073	# if bondStr not in BondList and revbondStr not in BondList:
3074	# BondList.append(bondStr)
3075	# notNames.append(FrstName)
3076	# return Res,BondList
3077
3078	# def FindAngles(Phase,RMCPdict):
3079	# generalData = Phase['General']
3080	# Cell = generalData['Cell'][1:7]
3081	# Amat = G2lat.cell2AB(Cell)[0]
3082	# cx,ct,cs,cia = generalData['AtomPtrs']
3083	# atomData = Phase['Atoms']
3084	# AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
3085	# AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
3086	# Angles = RMCPdict['Angles']
3087	# AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
3088	# AngleList = []
3089	# for MidName in AtDict:
3090	# nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
3091	# if len(nbrs) < 2: #need 2 neighbors to make an angle
3092	# continue
3093	# Atyp2 = AtDict[MidName]
3094	# for i,nbr1 in enumerate(nbrs):
3095	# Atyp1 = AtDict[nbr1[0]]
3096	# for j,nbr3 in enumerate(nbrs[i+1:]):
3097	# Atyp3 = AtDict[nbr3[0]]
3098	# IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
3099	# try:
3100	# angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
3101	# except KeyError:
3102	# try:
3103	# angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
3104	# except KeyError:
3105	# continue
3106	# XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
3107	# calAngle = G2mth.getRestAngle(XYZ,Amat)
3108	# if angData[0] <= calAngle <= angData[1]:
3109	# angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
3110	# revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
3111	# if angStr not in AngleList and revangStr not in AngleList:
3112	# AngleList.append(angStr)
3113	# return AngleList
3114
3115	# def GetSqConvolution(XY,d):
3116
3117	# n = XY.shape[1]
3118	# snew = np.zeros(n)
3119	# dq = np.zeros(n)
3120	# sold = XY[1]
3121	# q = XY[0]
3122	# dq[1:] = np.diff(q)
3123	# dq[0] = dq[1]
3124
3125	# for j in range(n):
3126	# for i in range(n):
3127	# b = abs(q[i]-q[j])
3128	# t = q[i]+q[j]
3129	# if j == i:
3130	# snew[j] += q[i]sold[i](d-np.sin(td)/t)dq[i]
3131	# else:
3132	# snew[j] += q[i]sold[i](np.sin(bd)/b-np.sin(td)/t)*dq[i]
3133	# snew[j] /= np.pi*q[j]
3134
3135	# snew[0] = snew[1]
3136	# return snew
3137
3138	# def GetMaxSphere(pdbName):
3139	# try:
3140	# pFil = open(pdbName,'r')
3141	# except FileNotFoundError:
3142	# return None
3143	# while True:
3144	# line = pFil.readline()
3145	# if 'Boundary' in line:
3146	# line = line.split()[3:]
3147	# G = np.array([float(item) for item in line])
3148	# G = np.reshape(G,(3,3))**2
3149	# G = nl.inv(G)
3150	# pFil.close()
3151	# break
3152	# dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
3153	# return min(dspaces)
3154
3155	def findfullrmc():
3156	'''Find where fullrmc is installed. Tries the following:
3157
3158	1. Returns the Config var 'fullrmc_exec', if defined. No check
3159	is done that the interpreter has fullrmc
3160	2. The current Python interpreter if fullrmc can be imported
3161	and fullrmc is version 5+
3162	3. The path is checked for a fullrmc image as named by Bachir
3163
3164	:returns: the full path to a python executable that is assumed to
3165	have fullrmc installed or None, if it was not found.
3166	'''
3167	if GSASIIpath.GetConfigValue('fullrmc_exec') is not None and is_exe(
3168	GSASIIpath.GetConfigValue('fullrmc_exec')):
3169	return GSASIIpath.GetConfigValue('fullrmc_exec')
3170	try:
3171	import fullrmc
3172	if int(fullrmc.__version__.split('.')[0]) >= 5:
3173	return sys.executable
3174	except:
3175	pass
3176	pathlist = os.environ["PATH"].split(os.pathsep)
3177	for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(),
3178	os.path.split(sys.executable)[0]):
3179	if p not in pathlist: pathlist.insert(0,p)
3180	import glob
3181	for p in pathlist:
3182	if sys.platform == "darwin":
3183	lookfor = "fullrmcmacOSi386-64bit"
3184	elif sys.platform == "win32":
3185	lookfor = "fullrmc*.exe"
3186	else:
3187	lookfor = "fullrmc*"
3188	fl = glob.glob(lookfor)
3189	if len(fl) > 0:
3190	return os.path.abspath(sorted(fl)[0])
3191
3192	def findPDFfit():
3193	'''Find if PDFfit2 is installed (may be local to GSAS-II). Does the following:
3194	:returns: two items: (1) the full path to a python executable or None, if
3195	it was not found and (2) path(s) to the PDFfit2 location(s) as a list.
3196
3197	'''
3198	if GSASIIpath.GetConfigValue('pdffit2_exec') is not None and is_exe(
3199	GSASIIpath.GetConfigValue('pdffit2_exec')):
3200	return GSASIIpath.GetConfigValue('pdffit2_exec'),None
3201	pdffitloc = os.path.join(GSASIIpath.path2GSAS2,'PDFfit2')
3202	if not os.path.exists(pdffitloc):
3203	print('PDFfit2 not found in GSAS-II \n\t(expected in '+pdffitloc+')')
3204	return None,[]
3205	if pdffitloc not in sys.path: sys.path.append(pdffitloc)
3206	try:
3207	from diffpy.pdffit2 import PdfFit
3208	import diffpy
3209	import inspect
3210	pdffitloc = [os.path.dirname(os.path.dirname(inspect.getfile(diffpy)))]
3211	# is this the original version of diffpy (w/pdffit2.py)
3212	try:
3213	from diffpy.pdffit2 import pdffit2
3214	except ImportError:
3215	# or the GSAS-II version w/o; for this we need to find the binary's location
3216	try:
3217	import pdffit2 # added for GSAS-II to relocate binary file
3218	except ImportError:
3219	print('\nError: pdffit2 failed to load with this python\n')
3220	return None,[]
3221	except ModuleNotFoundError:
3222	print('\nGSAS-II does not have a PDFfit2 module compatible\nwith this Python interpreter\n')
3223	return None,[]
3224	pdffitloc += [os.path.dirname(inspect.getfile(pdffit2))]
3225	return sys.executable,pdffitloc
3226	except Exception as msg:
3227	print('Error importing PDFfit2:\n',msg)
3228	return None,[]
3229
3230	def GetPDFfitAtomVar(Phase,RMCPdict):
3231	''' Find dict of independent "@n" variables for PDFfit in atom constraints
3232	'''
3233	General = Phase['General']
3234	Atoms = Phase['Atoms']
3235	cx,ct,cs,cia = General['AtomPtrs']
3236	AtomVar = RMCPdict['AtomVar']
3237	varnames = []
3238	for iat,atom in enumerate(RMCPdict['AtomConstr']):
3239	for it,item in enumerate(atom):
3240	if it > 1 and item:
3241	itms = item.split('@')
3242	for itm in itms[1:]:
3243	itnum = itm[:2]
3244	varname = '@%s'%itnum
3245	varnames.append(varname)
3246	if it < 6:
3247	if varname not in AtomVar:
3248	AtomVar[varname] = 0.0 #put ISODISTORT mode displ here?
3249	else:
3250	for i in range(3):
3251	if varname not in AtomVar:
3252	AtomVar[varname] = Atoms[iat][cia+i+2]
3253	varnames = set(varnames)
3254	for name in list(AtomVar.keys()): #clear out unused parameters
3255	if name not in varnames:
3256	del AtomVar[name]
3257
3258	def MakePDFfitAtomsFile(Phase,RMCPdict):
3259	'''Make the PDFfit atoms file
3260	'''
3261	General = Phase['General']
3262	if General['SGData']['SpGrp'] != 'P 1':
3263	return 'Space group symmetry must be lowered to P 1 for PDFfit'
3264	fName = General['Name']+'-PDFfit.stru'
3265	fName = fName.replace(' ','_')
3266	if 'sequential' in RMCPdict['refinement']:
3267	fName = 'Sequential_PDFfit.stru'
3268	fatm = open(fName,'w')
3269	fatm.write('title structure of '+General['Name']+'\n')
3270	fatm.write('format pdffit\n')
3271	fatm.write('scale 1.000000\n') #fixed
3272	sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut'])
3273	fatm.write('sharp '+sharp)
3274	shape = ''
3275	if RMCPdict['shape'] == 'sphere' and RMCPdict['spdiameter'][0] > 0.:
3276	shape = ' sphere, %10.6f\n'%RMCPdict['spdiameter'][0]
3277	elif RMCPdict['stepcut'] > 0.:
3278	shape = 'stepcut, %10.6f\n'%RMCPdict['stepcut']
3279	if shape:
3280	fatm.write('shape '+shape)
3281	fatm.write('spcgr %s\n'%RMCPdict['SGData']['SpGrp'].replace(' ',''))
3282	cell = General['Cell'][1:7]
3283	fatm.write('cell %10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n'%(
3284	cell[0],cell[1],cell[2],cell[3],cell[4],cell[5]))
3285	fatm.write('dcell '+5*' 0.000000,'+' 0.000000\n')
3286	Atoms = Phase['Atoms']
3287	fatm.write('ncell %8d,%8d,%8d,%10d\n'%(1,1,1,len(Atoms)))
3288	fatm.write('atoms\n')
3289	cx,ct,cs,cia = General['AtomPtrs']
3290	for atom in Atoms:
3291	fatm.write('%4s%18.8f%18.8f%18.8f%13.4f\n'%(atom[ct][:2].ljust(2),atom[cx],atom[cx+1],atom[cx+2],atom[cx+3]))
3292	fatm.write(' '+'%18.8f%18.8f%18.8f%13.4f\n'%(0.,0.,0.,0.))
3293	fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(atom[cia+2],atom[cia+3],atom[cia+4]))
3294	fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.,))
3295	fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(atom[cia+5],atom[cia+6],atom[cia+7]))
3296	fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.))
3297	fatm.close()
3298
3299	def MakePDFfitRunFile(Phase,RMCPdict):
3300	'''Make the PDFfit python run file
3301	'''
3302
3303	def GetCellConstr(SGData):
3304	if SGData['SGLaue'] in ['m3', 'm3m']:
3305	return [1,1,1,0,0,0]
3306	elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']:
3307	return [1,1,2,0,0,0]
3308	elif SGData['SGLaue'] in ['3R','3mR']:
3309	return [1,1,1,2,2,2]
3310	elif SGData['SGLaue'] == 'mmm':
3311	return [1,2,3,0,0,0]
3312	elif SGData['SGLaue'] == '2/m':
3313	if SGData['SGUniq'] == 'a':
3314	return [1,2,3,4,0,0]
3315	elif SGData['SGUniq'] == 'b':
3316	return [1,2,3,0,4,0]
3317	elif SGData['SGUniq'] == 'c':
3318	return [1,2,3,0,0,4]
3319	else:
3320	return [1,2,3,4,5,6]
3321
3322	General = Phase['General']
3323	Cell = General['Cell'][1:7]
3324	rundata = '''#!/usr/bin/env python
3325	# -- coding: utf-8 --
3326	import sys
3327	'''
3328	PDFfit_exe,PDFfit_path = findPDFfit() # returns python loc and path(s) for pdffit
3329	if not PDFfit_exe:
3330	GSASIIpath.IPyBreak()
3331	return None
3332	for p in PDFfit_path:
3333	rundata += "sys.path.append(r'{:}')\n".format(p)
3334	rundata += 'from diffpy.pdffit2 import PdfFit\n'
3335	rundata += 'pf = PdfFit()\n'
3336	Nd = 0
3337	Np = 0
3338	parms = {}
3339	parmNames = {}
3340	if 'sequential' in RMCPdict['refinement']:
3341	Np = 3
3342	rundata += '#sequential data here\n'
3343	else:
3344	for fil in RMCPdict['files']:
3345	filNam = RMCPdict['files'][fil][0]
3346	if 'Select' in filNam:
3347	continue
3348	if 'Neutron' in fil:
3349	Nd += 1
3350	dType = 'Ndata'
3351	else:
3352	Nd += 1
3353	dType = 'Xdata'
3354	filNam = os.path.abspath(filNam)
3355	rundata += "pf.read_data(r'%s', '%s', 30.0, %.4f)\n"%(filNam,dType[0],RMCPdict[dType]['qdamp'][0])
3356	rundata += 'pf.setdata(%d)\n'%Nd
3357	rundata += 'pf.pdfrange(%d, %6.2f, %6.2f)\n'%(Nd,RMCPdict[dType]['Fitrange'][0],RMCPdict[dType]['Fitrange'][1])
3358	for item in ['dscale','qdamp','qbroad']:
3359	if RMCPdict[dType][item][1]:
3360	Np += 1
3361	rundata += 'pf.constrain(pf.%s(),"@%d")\n'%(item,Np)
3362	parms[Np] = RMCPdict[dType][item][0]
3363	parmNames[Np] = item
3364	fName = General['Name']+'-PDFfit.stru'
3365	fName = fName.replace(' ','_')
3366	if 'sequential' in RMCPdict['refinement']:
3367	fName = 'Sequential_PDFfit.stru'
3368	Np = 9
3369	rundata += "pf.read_struct(r'{:}')\n".format(os.path.abspath(fName))
3370	for item in ['delta1','delta2','sratio']:
3371	if RMCPdict[item][1]:
3372	Np += 1
3373	rundata += 'pf.constrain(pf.%s,"@%d")\n'%(item,Np)
3374	parms[Np] = RMCPdict[item][0]
3375	parmNames[Np] = item
3376	if 'sphere' in RMCPdict['shape'] and RMCPdict['spdiameter'][1]:
3377	Np += 1
3378	rundata += 'pf.constrain(pf.spdiameter,"@%d")\n'%Np
3379	parms[Np] = RMCPdict['spdiameter'][0]
3380	parmNames[Np] = 'spdiameter'
3381
3382	if RMCPdict['cellref']:
3383	cellconst = GetCellConstr(RMCPdict['SGData'])
3384	used = []
3385	cellNames = ['a','b','c','alpha','beta','gamma']
3386	for ic in range(6):
3387	if cellconst[ic]:
3388	rundata += 'pf.constrain(pf.lat(%d), "@%d")\n'%(ic+1,Np+cellconst[ic])
3389	if cellconst[ic] not in used:
3390	parms[Np+cellconst[ic]] = Cell[ic]
3391	parmNames[Np+cellconst[ic]] = cellNames[ic]
3392	used.append(cellconst[ic])
3393	#Atom constraints here -------------------------------------------------------
3394	AtomVar = RMCPdict['AtomVar']
3395	used = []
3396	for iat,atom in enumerate(RMCPdict['AtomConstr']):
3397	for it,item in enumerate(atom):
3398	names = ['pf.x(%d)'%(iat+1),'pf.y(%d)'%(iat+1),'pf.z(%d)'%(iat+1),'pf.occ(%d)'%(iat+1)]
3399	if it > 1 and item:
3400	itms = item.split('@')
3401	once = False
3402	for itm in itms[1:]:
3403	try:
3404	itnum = int(itm[:2])
3405	except ValueError:
3406	print(' * ERROR - invalid string in atom constraint %s *'%(item))
3407	return None
3408	if it < 6:
3409	if not once:
3410	rundata += 'pf.constrain(%s,"%s")\n'%(names[it-2],item)
3411	once = True
3412	if itnum not in used:
3413	parms[itnum] = AtomVar['@%d'%itnum]
3414	parmNames[itnum] = names[it-2].split('.')[1]
3415	used.append(itnum)
3416	else:
3417	uijs = ['pf.u11(%d)'%(iat+1),'pf.u22(%d)'%(iat+1),'pf.u33(%d)'%(iat+1)]
3418	for i in range(3):
3419	rundata += 'pf.constrain(%s,"%s")\n'%(uijs[i],item)
3420	if itnum not in used:
3421	parms[itnum] = AtomVar['@%d'%itnum]
3422	parmNames[itnum] = uijs[i].split('.')[1]
3423	used.append(itnum)
3424
3425	if 'sequential' in RMCPdict['refinement']:
3426	rundata += '#parameters here\n'
3427	RMCPdict['Parms'] = parms #{'n':val,...}
3428	RMCPdict['ParmNames'] = parmNames #{'n':name,...}
3429	else:
3430	# set parameter values
3431	for iprm in parms:
3432	rundata += 'pf.setpar(%d,%.6f)\n'%(iprm,parms[iprm])
3433
3434	# Save results ---------------------------------------------------------------
3435	rundata += 'pf.refine()\n'
3436	if 'sequential' in RMCPdict['refinement']:
3437	fName = 'Sequential_PDFfit'
3438	rfile = open('Seq_PDFfit_template.py','w')
3439	rundata += 'pf.save_pdf(1, "%s")\n'%(fName+'.fgr')
3440	else:
3441	fName = General['Name'].replace(' ','_')+'-PDFfit'
3442	rfile = open(fName+'.py','w')
3443	Nd = 0
3444	for file in RMCPdict['files']:
3445	if 'Select' in RMCPdict['files'][file][0]: #skip unselected
3446	continue
3447	Nd += 1
3448	rundata += 'pf.save_pdf(%d, "%s")\n'%(Nd,fName+file[0]+'.fgr')
3449
3450	rundata += 'pf.save_struct(1, "%s")\n'%(fName+'.rstr')
3451	rundata += 'pf.save_res("%s")\n'%(fName+'.res')
3452
3453	rfile.writelines(rundata)
3454	rfile.close()
3455
3456	return fName+'.py'
3457
3458	def GetSeqCell(SGData,parmDict):
3459	''' For use in processing PDFfit sequential results
3460	'''
3461	try:
3462	if SGData['SGLaue'] in ['m3', 'm3m']:
3463	cell = [parmDict['11'][0],parmDict['11'][0],parmDict['11'][0],90.,90.,90.]
3464	elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']:
3465	cell = [parmDict['11'][0],parmDict['11'][0],parmDict['12'][0],90.,90.,90.]
3466	elif SGData['SGLaue'] in ['3R','3mR']:
3467	cell = [parmDict['11'][0],parmDict['11'][0],parmDict['11'][0],
3468	parmDict['12'][0],parmDict['12'][0],parmDict['12'][0]]
3469	elif SGData['SGLaue'] == 'mmm':
3470	cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,90.,90.]
3471	elif SGData['SGLaue'] == '2/m':
3472	if SGData['SGUniq'] == 'a':
3473	cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],parmDict['14'][0],90.,90.]
3474	elif SGData['SGUniq'] == 'b':
3475	cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,parmDict['14'][0],90.]
3476	elif SGData['SGUniq'] == 'c':
3477	cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,90.,parmDict['14'][0]]
3478	else:
3479	cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],
3480	parmDict['14'][0],parmDict['15'][0],parmDict['16'][0]]
3481	return G2lat.cell2A(cell)
3482	except KeyError:
3483	return None
3484
3485	def UpdatePDFfit(Phase,RMCPdict):
3486	''' Updates various PDFfit parameters held in GSAS-II
3487	'''
3488
3489	General = Phase['General']
3490	if RMCPdict['refinement'] == 'normal':
3491	fName = General['Name']+'-PDFfit.rstr'
3492	try:
3493	rstr = open(fName.replace(' ','_'),'r')
3494	except FileNotFoundError:
3495	return [fName,'Not found - PDFfit failed']
3496	lines = rstr.readlines()
3497	rstr.close()
3498	header = [line[:-1].split(' ',1) for line in lines[:7]]
3499	resdict = dict(header)
3500	for item in ['sharp','cell']:
3501	resdict[item] = [float(val) for val in resdict[item].split(',')]
3502	General['Cell'][1:7] = resdict['cell']
3503	for inam,name in enumerate(['delta2','delta1','sratio']):
3504	RMCPdict[name][0] = float(resdict['sharp'][inam])
3505	if 'shape' in resdict:
3506	if 'sphere' in resdict['shape']:
3507	RMCPdict['spdiameter'][0] = float(resdict['shape'].split()[-1])
3508	else:
3509	RMCPdict['stepcut'][0] = float(resdict['shape'][-1])
3510	cx,ct,cs,ci = G2mth.getAtomPtrs(Phase)
3511	Atoms = Phase['Atoms']
3512	atmBeg = 0
3513	for line in lines:
3514	atmBeg += 1
3515	if 'atoms' in line:
3516	break
3517	for atom in Atoms:
3518	atstr = lines[atmBeg][:-1].split()
3519	Uiistr = lines[atmBeg+2][:-1].split()
3520	Uijstr = lines[atmBeg+4][:-1].split()
3521	atom[cx:cx+4] = [float(atstr[1]),float(atstr[2]),float(atstr[3]),float(atstr[4])]
3522	atom[ci] = 'A'
3523	atom[ci+2:ci+5] = [float(Uiistr[0]),float(Uiistr[1]),float(Uiistr[2])]
3524	atom[ci+5:ci+8] = [float(Uijstr[0]),float(Uijstr[1]),float(Uijstr[2])]
3525	atmBeg += 6
3526	fName = General['Name']+'-PDFfit.res'
3527	else:
3528	fName = 'Sequential_PDFfit.res'
3529	try:
3530	res = open(fName.replace(' ','_'),'r')
3531	except FileNotFoundError:
3532	return [fName,'Not found - PDFfit failed']
3533	lines = res.readlines()
3534	res.close()
3535	Ibeg = False
3536	resline = ''
3537	XNdata = {'Xdata':RMCPdict['Xdata'],'Ndata':RMCPdict['Ndata']}
3538	for line in lines:
3539	if 'Radiation' in line and 'X-Rays' in line:
3540	dkey = 'Xdata'
3541	if 'Radiation' in line and'Neutrons' in line:
3542	dkey = 'Ndata'
3543	if 'Qdamp' in line and '(' in line:
3544	XNdata[dkey]['qdamp'][0] = float(line.split()[4])
3545	if 'Qbroad' in line and '(' in line:
3546	XNdata[dkey]['qbroad'][0] = float(line.split()[4])
3547	if 'Scale' in line and '(' in line:
3548	XNdata[dkey]['dscale'][0] = float(line.split()[3])
3549
3550	for iline,line in enumerate(lines):
3551	if 'Refinement parameters' in line:
3552	Ibeg = True
3553	continue
3554	if Ibeg:
3555	if '---------' in line:
3556	break
3557	resline += line[:-1]
3558	for iline,line in enumerate(lines):
3559	if 'Rw - ' in line:
3560	if 'nan' in line:
3561	Rwp = 100.0
3562	else:
3563	Rwp = float(line.split(':')[1])
3564	results = resline.replace('(','').split(')')[:-1]
3565	results = ['@'+result.lstrip() for result in results]
3566	results = [item.split() for item in results]
3567	RMCPdict['Parms'] = dict([[item[0][1:-1],float(item[1])] for item in results]) #{'n':val,...}
3568	if RMCPdict['refinement'] == 'normal':
3569	fName = General['Name']+'-PDFfit.py'
3570	py = open(fName.replace(' ','_'),'r')
3571	pylines = py.readlines()
3572	py.close()
3573	py = open(fName.replace(' ','_'),'w')
3574	newpy = []
3575	for pyline in pylines:
3576	if 'setpar' in pyline:
3577	parm = pyline.split('(')[1].split(',')[0]
3578	newpy.append('pf.setpar(%s,%.5f)\n'%(parm,RMCPdict['Parms'][parm]))
3579	else:
3580	newpy.append(pyline)
3581	py.writelines(newpy)
3582	py.close()
3583	RMCPdict.update(XNdata)
3584	results = dict([[item[0][:-1],float(item[1])] for item in results if item[0][:-1] in RMCPdict['AtomVar']])
3585	RMCPdict['AtomVar'].update(results)
3586	return None
3587	else: #sequential
3588	newParms = dict([[item[0][1:-1],[float(item[1]),float(item[2])]] for item in results]) #{'n':[val,esd],...}
3589	return newParms,Rwp
3590
3591	def MakefullrmcRun(pName,Phase,RMCPdict):
3592	'''Creates a script to run fullrmc. Returns the name of the file that was
3593	created.
3594	'''
3595	BondList = {}
3596	for k in RMCPdict['Pairs']:
3597	if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
3598	BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
3599	AngleList = []
3600	for angle in RMCPdict['Angles']:
3601	if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
3602	continue
3603	for i in (0,1,2):
3604	angle[i] = angle[i].strip()
3605	AngleList.append(angle)
3606	# rmin = RMCPdict['min Contact']
3607	cell = Phase['General']['Cell'][1:7]
3608	SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
3609	cx,ct,cs,cia = Phase['General']['AtomPtrs']
3610	atomsList = []
3611	for atom in Phase['Atoms']:
3612	el = ''.join([i for i in atom[ct] if i.isalpha()])
3613	atomsList.append([el] + atom[cx:cx+4])
3614	projDir,projName = os.path.split(os.path.abspath(pName))
3615	scrname = pName+'-fullrmc.py'
3616	restart = '%s_restart.pdb'%pName
3617	Files = RMCPdict['files']
3618	rundata = ''
3619	rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
3620	rundata += '# created in '+__file__+" v"+filversion.split()[1]
3621	rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
3622	rundata += '''
3623	# fullrmc imports (all that are potentially useful)
3624	import os,glob
3625	import time
3626	import pickle
3627	import numpy as np
3628	from fullrmc.Core import Collection
3629	from fullrmc.Engine import Engine
3630	import fullrmc.Constraints.PairDistributionConstraints as fPDF
3631	from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
3632	from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
3633	from fullrmc.Constraints.BondConstraints import BondConstraint
3634	from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
3635	from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
3636	from fullrmc.Generators.Swaps import SwapPositionsGenerator
3637	# utility routines
3638	def writeHeader(ENGINE,statFP):
3639	'header for stats file'
3640	statFP.write('generated-steps, total-error, ')
3641	for c in ENGINE.constraints:
3642	statFP.write(c.constraintName)
3643	statFP.write(', ')
3644	statFP.write('\\n')
3645	statFP.flush()
3646
3647	def writeCurrentStatus(ENGINE,statFP,plotF):
3648	'line in stats file & current constraint plots'
3649	statFP.write(str(ENGINE.generated))
3650	statFP.write(', ')
3651	statFP.write(str(ENGINE.totalStandardError))
3652	statFP.write(', ')
3653	for c in ENGINE.constraints:
3654	statFP.write(str(c.standardError))
3655	statFP.write(', ')
3656	statFP.write('\\n')
3657	statFP.flush()
3658	mpl.use('agg')
3659	fp = open(plotF,'wb')
3660	for c in ENGINE.constraints:
3661	p = c.plot(show=False)
3662	p[0].canvas.draw()
3663	image = p[0].canvas.buffer_rgba()
3664	pickle.dump(c.constraintName,fp)
3665	pickle.dump(np.array(image),fp)
3666	fp.close()
3667
3668	def calcRmax(ENGINE):
3669	'from Bachir, works for non-orthorhombic cells'
3670	a,b,c = ENGINE.basisVectors
3671	lens = []
3672	ts = np.linalg.norm(np.cross(a,b))/2
3673	lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(b)] )
3674	ts = np.linalg.norm(np.cross(b,c))/2
3675	lens.extend( [ts/np.linalg.norm(b), ts/np.linalg.norm(c)] )
3676	ts = np.linalg.norm(np.cross(a,c))/2
3677	lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(c)] )
3678	return min(lens)
3679	'''
3680	rundata += '''
3681	### When True, erases an existing engine to provide a fresh start
3682	FRESH_START = {:}
3683	dirName = "{:}"
3684	prefix = "{:}"
3685	project = prefix + "-fullrmc"
3686	time0 = time.time()
3687	'''.format(RMCPdict['ReStart'][0],projDir,projName)
3688
3689	rundata += '# setup structure\n'
3690	rundata += 'cell = ' + str(cell) + '\n'
3691	rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
3692	rundata += 'atomList = ' + str(atomsList).replace('],','],\n ') + '\n'
3693	rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
3694
3695	rundata += '\n# initialize engine\n'
3696	rundata += '''
3697	engineFileName = os.path.join(dirName, project + '.rmc')
3698	projectStats = os.path.join(dirName, project + '.stats')
3699	projectPlots = os.path.join(dirName, project + '.plots')
3700	pdbFile = os.path.join(dirName, project + '_restart.pdb')
3701	# check Engine exists if so (and not FRESH_START) load it
3702	# otherwise build it
3703	ENGINE = Engine(path=None)
3704	if not ENGINE.is_engine(engineFileName) or FRESH_START:
3705	## create structure
3706	ENGINE = Engine(path=engineFileName, freshStart=True)
3707	ENGINE.build_crystal_set_pdb(symOps = SymOpList,
3708	atoms = atomList,
3709	unitcellBC = cell,
3710	supercell = supercell)
3711	'''
3712	import atmdata
3713	# rundata += ' # conversion factors (may be needed)\n'
3714	# rundata += ' sumCiBi2 = 0.\n'
3715	# for elem in Phase['General']['AtomTypes']:
3716	# rundata += ' Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
3717	# rundata += ' sumCiBi2 += (Ci{})*2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
3718	rundata += ' rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
3719	# settings that require a new Engine
3720	for File in Files:
3721	filDat = RMCPdict['files'][File]
3722	if not os.path.exists(filDat[0]): continue
3723	sfwt = 'neutronCohb'
3724	if 'Xray' in File:
3725	sfwt = 'atomicNumber'
3726	if 'G(r)' in File:
3727	rundata += ' GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3728	if filDat[3] == 0:
3729	#rundata += ''' # read and xform G(r) as defined in RMCProfile
3730	# see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
3731	#rundata += ' GR[1] = 4 np.pi * GR[0] * rho0 / sumCiBi2\n'
3732	#rundata += ' GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3733	rundata += ' # G(r) as defined in RMCProfile\n'
3734	rundata += ' GofR = fullrmc.Constraints.RadialDistributionConstraints.RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3735	elif filDat[3] == 1:
3736	rundata += ' # This is G(r) as defined in PDFFIT\n'
3737	rundata += ' GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3738	elif filDat[3] == 2:
3739	rundata += ' # This is g(r)\n'
3740	rundata += ' GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3741	else:
3742	raise ValueError('Invalid G(r) type: '+str(filDat[3]))
3743	rundata += ' ENGINE.add_constraints([GofR])\n'
3744	rundata += ' GofR.set_limits((None, calcRmax(ENGINE)))\n'
3745	elif '(Q)' in File:
3746	rundata += ' SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3747	if filDat[3] == 0:
3748	rundata += ' # F(Q) as defined in RMCProfile\n'
3749	#rundata += ' SOQ[1] *= 1 / sumCiBi2\n'
3750	if filDat[4]:
3751	rundata += ' SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3752	rundata += ' SofQ = fullrmc.Constraints.StructureFactorConstraints.NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3753	elif filDat[3] == 1:
3754	rundata += ' # S(Q) as defined in PDFFIT\n'
3755	rundata += ' SOQ[1] -= 1\n'
3756	if filDat[4]:
3757	rundata += ' SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3758	rundata += ' SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3759	else:
3760	raise ValueError('Invalid S(Q) type: '+str(filDat[3]))
3761	rundata += ' ENGINE.add_constraints([SofQ])\n'
3762	else:
3763	print('What is this?')
3764	minDists = ''
3765	if BondList:
3766	rundata += ''' B_CONSTRAINT = BondConstraint()
3767	ENGINE.add_constraints(B_CONSTRAINT)
3768	B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
3769	'''
3770	for pair in BondList:
3771	e1,e2 = pair.split('-')
3772	d1,d2 = BondList[pair]
3773	if d1 == 0: continue
3774	if d2 == 0:
3775	minDists += '("element","{}","{}",{}),'.format(e1.strip(),e2.strip(),d1)
3776	else:
3777	rundata += ' ("element","{}","{}",{},{}),\n'.format(
3778	e1.strip(),e2.strip(),d1,d2)
3779	rundata += ' ])\n'
3780	rundata += ' D_CONSTRAINT = DistanceConstraint(defaultLowerDistance={})\n'.format(RMCPdict['min Contact'])
3781	if minDists:
3782	rundata += " D_CONSTRAINT.set_pairs_definition( {'inter':[" + minDists + "]})\n"
3783	rundata += ' ENGINE.add_constraints(D_CONSTRAINT)\n'
3784
3785	if AngleList:
3786	rundata += ''' A_CONSTRAINT = BondsAngleConstraint()
3787	ENGINE.add_constraints(A_CONSTRAINT)
3788	A_CONSTRAINT.create_supercell_angles(anglesDefinition=[
3789	'''
3790	for item in AngleList:
3791	rundata += (' '+
3792	'("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
3793	rundata += ' ])\n'
3794	rundata += '''
3795	for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
3796	ENGINE.save()
3797	else:
3798	ENGINE = ENGINE.load(path=engineFileName)
3799
3800	ENGINE.set_log_file(os.path.join(dirName,prefix))
3801	'''
3802	if RMCPdict['Swaps']:
3803	rundata += '\n#set up for site swaps\n'
3804	rundata += 'aN = ENGINE.allNames\n'
3805	rundata += 'SwapGen = {}\n'
3806	for swap in RMCPdict['Swaps']:
3807	rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
3808	rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
3809	rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
3810	rundata += ' for swaps in SwapGen:\n'
3811	rundata += ' AB = swaps.split("-")\n'
3812	rundata += ' ENGINE.set_groups_as_atoms()\n'
3813	rundata += ' for g in ENGINE.groups:\n'
3814	rundata += ' if aN[g.indexes[0]]==AB[0]:\n'
3815	rundata += ' g.set_move_generator(SwapGen[swaps][0])\n'
3816	rundata += ' elif aN[g.indexes[0]]==AB[1]:\n'
3817	rundata += ' g.set_move_generator(SwapGen[swaps][1])\n'
3818	rundata += ' sProb = SwapGen[swaps][2]\n'
3819	rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
3820	# rundata += 'wtDict = {}\n'
3821	# for File in Files:
3822	# filDat = RMCPdict['files'][File]
3823	# if not os.path.exists(filDat[0]): continue
3824	# if 'Xray' in File:
3825	# sfwt = 'atomicNumber'
3826	# else:
3827	# sfwt = 'neutronCohb'
3828	# if 'G(r)' in File:
3829	# typ = 'Pair'
3830	# elif '(Q)' in File:
3831	# typ = 'Struct'
3832	# rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
3833	rundata += 'for c in ENGINE.constraints: # loop over predefined constraints\n'
3834	rundata += ' if type(c) is fPDF.PairDistributionConstraint:\n'
3835	# rundata += ' c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
3836	rundata += ' c.set_limits((None,calcRmax(ENGINE)))\n'
3837	if RMCPdict['FitScale']:
3838	rundata += ' c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3839	# rundata += ' c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
3840	if RMCPdict['FitScale']:
3841	rundata += ' elif type(c) is ReducedStructureFactorConstraint:\n'
3842	rundata += ' c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3843	# torsions difficult to implement, must be internal to cell & named with
3844	# fullrmc atom names
3845	# if len(RMCPdict['Torsions']): # Torsions currently commented out in GUI
3846	# rundata += 'for c in ENGINE.constraints: # look for Dihedral Angle Constraints\n'
3847	# rundata += ' if type(c) is DihedralAngleConstraint:\n'
3848	# rundata += ' c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
3849	# rundata += ' c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
3850	# for torsion in RMCPdict['Torsions']:
3851	# rundata += ' %s\n'%str(tuple(torsion))
3852	# rundata += ' ]})\n'
3853	rundata += '''
3854	if FRESH_START:
3855	# initialize engine with one step to get starting config energetics
3856	ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
3857	statFP = open(projectStats,'w')
3858	writeHeader(ENGINE,statFP)
3859	writeCurrentStatus(ENGINE,statFP,projectPlots)
3860	else:
3861	statFP = open(projectStats,'a')
3862
3863	# setup runs for fullrmc
3864	'''
3865	rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
3866	rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
3867	rundata += ' ENGINE.set_groups_as_atoms()\n'
3868	rundata += ' expected = ENGINE.generated+steps\n'
3869
3870	rundata += ' ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
3871	rundata += ' writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
3872	rundata += ' if ENGINE.generated != expected: break # run was stopped\n'
3873	rundata += 'statFP.close()\n'
3874	rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
3875	rfile = open(scrname,'w')
3876	rfile.writelines(rundata)
3877	rfile.close()
3878	return scrname
3879
3880	def GetRMCBonds(general,RMCPdict,Atoms,bondList):
3881	bondDist = []
3882	Cell = general['Cell'][1:7]
3883	Supercell = RMCPdict['SuperCell']
3884	Trans = np.eye(3)*np.array(Supercell)
3885	Cell = G2lat.TransformCell(Cell,Trans)[:6]
3886	Amat,Bmat = G2lat.cell2AB(Cell)
3887	indices = (-1,0,1)
3888	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3889	for bonds in bondList:
3890	Oxyz = np.array(Atoms[bonds[0]][1:])
3891	Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])
3892	Dx = np.array([Txyz-Oxyz+unit for unit in Units])
3893	Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
3894	for dx in Dx.T:
3895	bondDist.append(np.min(dx))
3896	return np.array(bondDist)
3897
3898	def GetRMCAngles(general,RMCPdict,Atoms,angleList):
3899	bondAngles = []
3900	Cell = general['Cell'][1:7]
3901	Supercell = RMCPdict['SuperCell']
3902	Trans = np.eye(3)*np.array(Supercell)
3903	Cell = G2lat.TransformCell(Cell,Trans)[:6]
3904	Amat,Bmat = G2lat.cell2AB(Cell)
3905	indices = (-1,0,1)
3906	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3907	for angle in angleList:
3908	Oxyz = np.array(Atoms[angle[0]][1:])
3909	TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])
3910	TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])
3911	DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
3912	DAx = np.sqrt(np.sum(DAxV**2,axis=2))
3913	DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
3914	DBx = np.sqrt(np.sum(DBxV**2,axis=2))
3915	iDAx = np.argmin(DAx,axis=0)
3916	iDBx = np.argmin(DBx,axis=0)
3917	for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
3918	DAv = DAxV[iA,i]/DAx[iA,i]
3919	DBv = DBxV[iB,i]/DBx[iB,i]
3920	bondAngles.append(npacosd(np.sum(DAv*DBv)))
3921	return np.array(bondAngles)
3922
3923	def ISO2PDFfit(Phase):
3924	''' Creates new phase structure to be used for PDFfit from an ISODISTORT mode displacement phase.
3925	It builds the distortion mode parameters to be used as PDFfit variables for atom displacements from
3926	the original parent positions as transformed to the child cell wiht symmetry defined from ISODISTORT.
3927
3928	:param Phase: dict GSAS-II Phase structure; must contain ISODISTORT dict. NB: not accessed otherwise
3929
3930	:returns: dict: GSAS-II Phase structure; will contain ['RMC']['PDFfit'] dict
3931	'''
3932
3933	Trans = np.eye(3)
3934	Uvec = np.zeros(3)
3935	Vvec = np.zeros(3)
3936	Phase = copy.deepcopy(Phase)
3937	Atoms = Phase['Atoms']
3938	parentXYZ = Phase['ISODISTORT']['G2parentCoords'] #starting point for mode displacements
3939	cx,ct,cs,cia = Phase['General']['AtomPtrs']
3940	for iat,atom in enumerate(Atoms):
3941	atom[cx:cx+3] = parentXYZ[iat]
3942	SGData = copy.deepcopy(Phase['General']['SGData'])
3943	SGOps = SGData['SGOps']
3944	newPhase = copy.deepcopy(Phase)
3945	newPhase['ranId'] = rand.randint(0,sys.maxsize)
3946	newPhase['General']['Name'] += '_PDFfit'
3947	newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1] #this is for filled unit cell
3948	newPhase,atCodes = G2lat.TransformPhase(Phase,newPhase,Trans,Uvec,Vvec,False)
3949	newPhase['Histograms'] = {}
3950	newPhase['Drawing'] = []
3951	Atoms = newPhase['Atoms']
3952	RMCPdict = newPhase['RMC']['PDFfit']
3953	ISOdict = newPhase['ISODISTORT']
3954	RMCPdict['AtomConstr'] = []
3955	RMCPdict['SGData'] = copy.deepcopy(SGData) #this is from the ISODISTORT child; defines PDFfit constraints
3956	Norms = ISOdict['NormList']
3957	ModeMatrix = ISOdict['Mode2VarMatrix']
3958	RMCPdict['AtomVar'] = {'@%d'%(itm+21):val for itm,val in enumerate(ISOdict['modeDispl'])}
3959	for iatm,[atom,atcode] in enumerate(zip(Atoms,atCodes)):
3960	pid,opid = [int(item) for item in atcode.split(':')]
3961	atmConstr = [atom[ct-1],atom[ct],'','','','','',atcode]
3962	for ip,pname in enumerate(['%s_d%s'%(atom[ct-1],x) for x in ['x','y','z']]):
3963	try:
3964	conStr = ''
3965	if Atoms[iatm][cx+ip]:
3966	conStr += '%.5f'%Atoms[iatm][cx+ip]
3967	pid = ISOdict['IsoVarList'].index(pname)
3968	consVec = ModeMatrix[pid]
3969	for ic,citm in enumerate(consVec): #NB: this assumes orthorhombic or lower symmetry
3970	if opid < 0:
3971	citm *= -SGOps[100-opid%100-1][0][ip][ip] #remove centering, if any
3972	else:
3973	citm *= SGOps[opid%100-1][0][ip][ip]
3974	if citm > 0.:
3975	conStr += '+%.5f@%d'%(citmNorms[ic],ic+21)
3976	elif citm < 0.:
3977	conStr += '%.5f@%d'%(citmNorms[ic],ic+21)
3978	atmConstr[ip+2] = conStr
3979	except ValueError:
3980	atmConstr[ip+2] = ''
3981	RMCPdict['AtomConstr'].append(atmConstr)
3982	return newPhase
3983
3984	def GetAtmDispList(ISOdata):
3985	atmDispList = []
3986	MT = ISOdata['Mode2VarMatrix'].T
3987	DispList = ISOdata['IsoVarList']
3988	N = len(DispList)
3989	for I in range(N):
3990	vary = []
3991	for i in range(N):
3992	if MT[I,i]:
3993	vary.append(DispList[i])
3994	atmDispList.append(vary)
3995	return atmDispList
3996
3997	#### Reflectometry calculations ################################################################################
3998	def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
3999	G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
4000
4001	class RandomDisplacementBounds(object):
4002	"""random displacement with bounds"""
4003	def __init__(self, xmin, xmax, stepsize=0.5):
4004	self.xmin = xmin
4005	self.xmax = xmax
4006	self.stepsize = stepsize
4007
4008	def __call__(self, x):
4009	"""take a random step but ensure the new position is within the bounds"""
4010	while True:
4011	# this could be done in a much more clever way, but it will work for example purposes
4012	steps = self.xmax-self.xmin
4013	xnew = x + np.random.uniform(-self.stepsizesteps, self.stepsizesteps, np.shape(x))
4014	if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4015	break
4016	return xnew
4017
4018	def GetModelParms():
4019	parmDict = {}
4020	varyList = []
4021	values = []
4022	bounds = []
4023	parmDict['dQ type'] = data['dQ type']
4024	parmDict['Res'] = data['Resolution'][0]/(100.*sateln2) #% FWHM-->decimal sig
4025	for parm in ['Scale','FltBack']:
4026	parmDict[parm] = data[parm][0]
4027	if data[parm][1]:
4028	varyList.append(parm)
4029	values.append(data[parm][0])
4030	bounds.append(Bounds[parm])
4031	parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
4032	parmDict['nLayers'] = len(parmDict['Layer Seq'])
4033	for ilay,layer in enumerate(data['Layers']):
4034	name = layer['Name']
4035	cid = str(ilay)+';'
4036	parmDict[cid+'Name'] = name
4037	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
4038	parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
4039	if layer.get(parm,[0.,False])[1]:
4040	varyList.append(cid+parm)
4041	value = layer[parm][0]
4042	values.append(value)
4043	if value:
4044	bound = [value*Bfac,value/Bfac]
4045	else:
4046	bound = [0.,10.]
4047	bounds.append(bound)
4048	if name not in ['vacuum','unit scatter']:
4049	parmDict[cid+'rho'] = Substances[name]['Scatt density']
4050	parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
4051	return parmDict,varyList,values,bounds
4052
4053	def SetModelParms():
4054	line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
4055	if 'Scale' in varyList:
4056	data['Scale'][0] = parmDict['Scale']
4057	line += ' esd: %.4g'%(sigDict['Scale'])
4058	G2fil.G2Print (line)
4059	line = ' Flat background: %15.4g'%(parmDict['FltBack'])
4060	if 'FltBack' in varyList:
4061	data['FltBack'][0] = parmDict['FltBack']
4062	line += ' esd: %15.3g'%(sigDict['FltBack'])
4063	G2fil.G2Print (line)
4064	for ilay,layer in enumerate(data['Layers']):
4065	name = layer['Name']
4066	G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
4067	cid = str(ilay)+';'
4068	line = ' '
4069	line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
4070	line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
4071	for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
4072	if parm in layer:
4073	if parm == 'Rough':
4074	layer[parm][0] = abs(parmDict[cid+parm]) #make positive
4075	else:
4076	layer[parm][0] = parmDict[cid+parm]
4077	line += ' %s: %.3f'%(parm,layer[parm][0])
4078	if cid+parm in varyList:
4079	line += ' esd: %.3g'%(sigDict[cid+parm])
4080	G2fil.G2Print (line)
4081	G2fil.G2Print (line2)
4082
4083	def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
4084	parmDict.update(zip(varyList,values))
4085	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
4086	return M
4087
4088	def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
4089	parmDict.update(zip(varyList,values))
4090	M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
4091	return np.sum(M**2)
4092
4093	def getREFD(Q,Qsig,parmDict):
4094	Ic = np.ones_like(Q)*parmDict['FltBack']
4095	Scale = parmDict['Scale']
4096	Nlayers = parmDict['nLayers']
4097	Res = parmDict['Res']
4098	depth = np.zeros(Nlayers)
4099	rho = np.zeros(Nlayers)
4100	irho = np.zeros(Nlayers)
4101	sigma = np.zeros(Nlayers)
4102	for ilay,lay in enumerate(parmDict['Layer Seq']):
4103	cid = str(lay)+';'
4104	depth[ilay] = parmDict[cid+'Thick']
4105	sigma[ilay] = parmDict[cid+'Rough']
4106	if parmDict[cid+'Name'] == u'unit scatter':
4107	rho[ilay] = parmDict[cid+'DenMul']
4108	irho[ilay] = parmDict[cid+'iDenMul']
4109	elif 'vacuum' != parmDict[cid+'Name']:
4110	rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
4111	irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
4112	if cid+'Mag SLD' in parmDict:
4113	rho[ilay] += parmDict[cid+'Mag SLD']
4114	if parmDict['dQ type'] == 'None':
4115	AB = abeles(0.5*Q,depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
4116	elif 'const' in parmDict['dQ type']:
4117	AB = SmearAbeles(0.5Q,QRes,depth,rho,irho,sigma[1:])
4118	else: #dQ/Q in data
4119	AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
4120	Ic += AB*Scale
4121	return Ic
4122
4123	def estimateT0(takestep):
4124	Mmax = -1.e-10
4125	Mmin = 1.e10
4126	for i in range(100):
4127	x0 = takestep(values)
4128	M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
4129	Mmin = min(M,Mmin)
4130	MMax = max(M,Mmax)
4131	return 1.5*(MMax-Mmin)
4132
4133	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
4134	if data.get('2% weight'):
4135	wt = 1./(0.02Io)*2
4136	Qmin = Limits[1][0]
4137	Qmax = Limits[1][1]
4138	wtFactor = ProfDict['wtFactor']
4139	Bfac = data['Toler']
4140	Ibeg = np.searchsorted(Q,Qmin)
4141	Ifin = np.searchsorted(Q,Qmax)+1 #include last point
4142	Ic[:] = 0
4143	Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
4144	'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
4145	parmDict,varyList,values,bounds = GetModelParms()
4146	Msg = 'Failed to converge'
4147	if varyList:
4148	if data['Minimizer'] == 'LMLS':
4149	result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
4150	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
4151	parmDict.update(zip(varyList,result[0]))
4152	chisq = np.sum(result[2]['fvec']**2)
4153	ncalc = result[2]['nfev']
4154	covM = result[1]
4155	newVals = result[0]
4156	elif data['Minimizer'] == 'Basin Hopping':
4157	xyrng = np.array(bounds).T
4158	take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
4159	T0 = estimateT0(take_step)
4160	G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
4161	result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
4162	interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
4163	'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
4164	chisq = result.fun
4165	ncalc = result.nfev
4166	newVals = result.x
4167	covM = []
4168	elif data['Minimizer'] == 'MC/SA Anneal':
4169	xyrng = np.array(bounds).T
4170	result = G2mth.anneal(sumREFD, values,
4171	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
4172	schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
4173	boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
4174	ranRange=0.20,autoRan=False,dlg=None)
4175	newVals = result[0]
4176	parmDict.update(zip(varyList,newVals))
4177	chisq = result[1]
4178	ncalc = result[3]
4179	covM = []
4180	G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
4181	elif data['Minimizer'] == 'L-BFGS-B':
4182	result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds, #ftol=Ftol,
4183	args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
4184	parmDict.update(zip(varyList,result['x']))
4185	chisq = result.fun
4186	ncalc = result.nfev
4187	newVals = result.x
4188	covM = []
4189	else: #nothing varied
4190	M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
4191	chisq = np.sum(M**2)
4192	ncalc = 0
4193	covM = []
4194	sig = []
4195	sigDict = {}
4196	result = []
4197	Rvals = {}
4198	Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]Io[Ibeg:Ifin]2))100. #to %
4199	Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList)) #reduced chi^2
4200	Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
4201	Ib[Ibeg:Ifin] = parmDict['FltBack']
4202	try:
4203	if not len(varyList):
4204	Msg += ' - nothing refined'
4205	raise ValueError
4206	Nans = np.isnan(newVals)
4207	if np.any(Nans):
4208	name = varyList[Nans.nonzero(True)[0]]
4209	Msg += ' Nan result for '+name+'!'
4210	raise ValueError
4211	Negs = np.less_equal(newVals,0.)
4212	if np.any(Negs):
4213	indx = Negs.nonzero()
4214	name = varyList[indx[0][0]]
4215	if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
4216	Msg += ' negative coefficient for '+name+'!'
4217	raise ValueError
4218	if len(covM):
4219	sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
4220	covMatrix = covM*Rvals['GOF']
4221	else:
4222	sig = np.zeros(len(varyList))
4223	covMatrix = []
4224	sigDict = dict(zip(varyList,sig))
4225	G2fil.G2Print (' Results of reflectometry data modelling fit:')
4226	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
4227	G2fil.G2Print ('Rwp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
4228	SetModelParms()
4229	return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
4230	except (ValueError,TypeError): #when bad LS refinement; covM missing or with nans
4231	G2fil.G2Print (Msg)
4232	return False,0,0,0,0,0,0,Msg
4233
4234	def makeSLDprofile(data,Substances):
4235
4236	sq2 = np.sqrt(2.)
4237	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
4238	Nlayers = len(laySeq)
4239	laySeq = np.array(laySeq,dtype=int)
4240	interfaces = np.zeros(Nlayers)
4241	rho = np.zeros(Nlayers)
4242	sigma = np.zeros(Nlayers)
4243	name = data['Layers'][0]['Name']
4244	thick = 0.
4245	for ilay,lay in enumerate(laySeq):
4246	layer = data['Layers'][lay]
4247	name = layer['Name']
4248	if 'Thick' in layer:
4249	thick += layer['Thick'][0]
4250	interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
4251	if 'Rough' in layer:
4252	sigma[ilay] = max(0.001,layer['Rough'][0])
4253	if name != 'vacuum':
4254	if name == 'unit scatter':
4255	rho[ilay] = np.sqrt(layer['DenMul'][0]2+layer['iDenMul'][0]2)
4256	else:
4257	rrho = Substances[name]['Scatt density']
4258	irho = Substances[name]['XImag density']
4259	rho[ilay] = np.sqrt(rrho2+irho2)*layer['DenMul'][0]
4260	if 'Mag SLD' in layer:
4261	rho[ilay] += layer['Mag SLD'][0]
4262	name = data['Layers'][-1]['Name']
4263	x = np.linspace(-0.15thick,1.15thick,1000,endpoint=True)
4264	xr = np.flipud(x)
4265	interfaces[-1] = x[-1]
4266	y = np.ones_like(x)*rho[0]
4267	iBeg = 0
4268	for ilayer in range(Nlayers-1):
4269	delt = rho[ilayer+1]-rho[ilayer]
4270	iPos = np.searchsorted(x,interfaces[ilayer])
4271	y[iBeg:] += (delt/2.)sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2sigma[ilayer+1]))
4272	iBeg = iPos
4273	return x,xr,y
4274
4275	def REFDModelFxn(Profile,Inst,Limits,Substances,data):
4276
4277	Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
4278	Qmin = Limits[1][0]
4279	Qmax = Limits[1][1]
4280	iBeg = np.searchsorted(Q,Qmin)
4281	iFin = np.searchsorted(Q,Qmax)+1 #include last point
4282	Ib[:] = data['FltBack'][0]
4283	Ic[:] = 0
4284	Scale = data['Scale'][0]
4285	if data['Layer Seq'] == []:
4286	return
4287	laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
4288	Nlayers = len(laySeq)
4289	depth = np.zeros(Nlayers)
4290	rho = np.zeros(Nlayers)
4291	irho = np.zeros(Nlayers)
4292	sigma = np.zeros(Nlayers)
4293	for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
4294	layer = data['Layers'][lay]
4295	name = layer['Name']
4296	if 'Thick' in layer: #skips first & last layers
4297	depth[ilay] = layer['Thick'][0]
4298	if 'Rough' in layer: #skips first layer
4299	sigma[ilay] = layer['Rough'][0]
4300	if 'unit scatter' == name:
4301	rho[ilay] = layer['DenMul'][0]
4302	irho[ilay] = layer['iDenMul'][0]
4303	else:
4304	rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
4305	irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
4306	if 'Mag SLD' in layer:
4307	rho[ilay] += layer['Mag SLD'][0]
4308	if data['dQ type'] == 'None':
4309	AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles
4310	elif 'const' in data['dQ type']:
4311	res = data['Resolution'][0]/(100.*sateln2)
4312	AB = SmearAbeles(0.5Q[iBeg:iFin],resQ[iBeg:iFin],depth,rho,irho,sigma[1:])
4313	else: #dQ/Q in data
4314	AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
4315	Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]