Context Navigation

source: trunk/GSASIIpwd.py @ 4821

Last change on this file since 4821 was 4821, checked in by toby, 2 years ago
new peakfit mode in scripting: allow peak width terms to remain unchanged (neither refined nor calc from inst. params.); not accesible from GUI.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 171.2 KB

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

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

Download in other formats: