Context Navigation

source: trunk/GSASIIpwd.py @ 4483

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

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

Download in other formats: