Context Navigation

source: trunk/GSASIIpwd.py @ 4518

Last change on this file since 4518 was 4518, checked in by toby, 3 years ago
start on fullrmc update for 5.0; fix for rigid bodies
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 166.5 KB

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

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

Download in other formats: