Context Navigation

source: trunk/GSASIIpwd.py @ 4523

Last change on this file since 4523 was 4523, checked in by toby, 3 years ago
more fullrmc work, but not done
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 168.9 KB

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

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

Download in other formats: