Context Navigation

source: trunk/GSASIIpwd.py @ 4800

Last change on this file since 4800 was 4800, checked in by toby, 2 years ago
finish up docs cleanups
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 169.9 KB

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

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

Download in other formats: