Context Navigation

source: trunk/GSASIIpwd.py @ 5050

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