Context Navigation

source: trunk/GSASIIstruct.py @ 869

Last change on this file since 869 was 869, checked in by vondreele, 10 years ago
a fix to restraint refinement
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 188.5 KB

Line
1	# -- coding: utf-8 --
2	#GSASIIstructure - structure computation routines
3	########### SVN repository information ###################
4	# $Date: 2013-03-21 15:54:41 +0000 (Thu, 21 Mar 2013) $
5	# $Author: vondreele $
6	# $Revision: 869 $
7	# $URL: trunk/GSASIIstruct.py $
8	# $Id: GSASIIstruct.py 869 2013-03-21 15:54:41Z vondreele $
9	########### SVN repository information ###################
10	import sys
11	import os
12	import os.path as ospath
13	import time
14	import math
15	import copy
16	import random
17	import cPickle
18	import numpy as np
19	import numpy.ma as ma
20	import numpy.linalg as nl
21	import scipy.optimize as so
22	import GSASIIpath
23	GSASIIpath.SetVersionNumber("$Revision: 869 $")
24	import GSASIIElem as G2el
25	import GSASIIlattice as G2lat
26	import GSASIIspc as G2spc
27	import GSASIIpwd as G2pwd
28	import GSASIImapvars as G2mv
29	import GSASIImath as G2mth
30
31	sind = lambda x: np.sin(x*np.pi/180.)
32	cosd = lambda x: np.cos(x*np.pi/180.)
33	tand = lambda x: np.tan(x*np.pi/180.)
34	asind = lambda x: 180.*np.arcsin(x)/np.pi
35	acosd = lambda x: 180.*np.arccos(x)/np.pi
36	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
37
38	ateln2 = 8.0*math.log(2.0)
39	DEBUG = False
40
41	def GetControls(GPXfile):
42	''' Returns dictionary of control items found in GSASII gpx file
43	input:
44	GPXfile = .gpx full file name
45	return:
46	Controls = dictionary of control items
47	'''
48	Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
49	'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
50	fl = open(GPXfile,'rb')
51	while True:
52	try:
53	data = cPickle.load(fl)
54	except EOFError:
55	break
56	datum = data[0]
57	if datum[0] == 'Controls':
58	Controls.update(datum[1])
59	fl.close()
60	return Controls
61
62	def GetConstraints(GPXfile):
63	'''Read the constraints from the GPX file and interpret them
64	'''
65	constList = []
66	fl = open(GPXfile,'rb')
67	while True:
68	try:
69	data = cPickle.load(fl)
70	except EOFError:
71	break
72	datum = data[0]
73	if datum[0] == 'Constraints':
74	constDict = datum[1]
75	for item in constDict:
76	constList += constDict[item]
77	fl.close()
78	constDict,fixedList,ignored = ProcessConstraints(constList)
79	if ignored:
80	print ignored,'old-style Constraints were rejected'
81	return constDict,fixedList
82
83	def ProcessConstraints(constList):
84	"interpret constraints"
85	constDict = []
86	fixedList = []
87	ignored = 0
88	for item in constList:
89	if item[-1] == 'h':
90	# process a hold
91	fixedList.append('0')
92	constDict.append({item[0][1]:0.0})
93	elif item[-1] == 'f':
94	# process a new variable
95	fixedList.append(None)
96	constDict.append({})
97	for term in item[:-3]:
98	constDict[-1][term[1]] = term[0]
99	#constFlag[-1] = ['Vary']
100	elif item[-1] == 'c':
101	# process a contraint relationship
102	fixedList.append(str(item[-3]))
103	constDict.append({})
104	for term in item[:-3]:
105	constDict[-1][term[1]] = term[0]
106	#constFlag[-1] = ['VaryFree']
107	elif item[-1] == 'e':
108	# process an equivalence
109	firstmult = None
110	eqlist = []
111	for term in item[:-3]:
112	if term[0] == 0: term[0] = 1.0
113	if firstmult is None:
114	firstmult,firstvar = term
115	else:
116	eqlist.append([term[1],firstmult/term[0]])
117	G2mv.StoreEquivalence(firstvar,eqlist)
118	else:
119	ignored += 1
120	return constDict,fixedList,ignored
121
122	def CheckConstraints(GPXfile):
123	'''Load constraints and related info and return any error or warning messages'''
124	# init constraints
125	G2mv.InitVars()
126	# get variables
127	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
128	if not Phases:
129	return 'Error: No Phases!',''
130	if not Histograms:
131	return 'Error: no diffraction data',''
132	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,Print=False)
133	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
134	histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
135	varyList = phaseVary+hapVary+histVary
136	constrDict,fixedList = GetConstraints(GPXfile)
137	return G2mv.CheckConstraints(varyList,constrDict,fixedList)
138
139	def GetRestraints(GPXfile):
140	'''Read the restraints from the GPX file
141	'''
142	fl = open(GPXfile,'rb')
143	while True:
144	try:
145	data = cPickle.load(fl)
146	except EOFError:
147	break
148	datum = data[0]
149	if datum[0] == 'Restraints':
150	restraintDict = datum[1]
151	fl.close()
152	return restraintDict
153
154	def GetRigidBodies(GPXfile):
155	'''Read the rigid body models from the GPX file
156	'''
157	fl = open(GPXfile,'rb')
158	while True:
159	try:
160	data = cPickle.load(fl)
161	except EOFError:
162	break
163	datum = data[0]
164	if datum[0] == 'Rigid bodies':
165	rigidbodyDict = datum[1]
166	fl.close()
167	return rigidbodyDict
168
169	def GetPhaseNames(GPXfile):
170	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
171	input:
172	GPXfile = gpx full file name
173	return:
174	PhaseNames = list of phase names
175	'''
176	fl = open(GPXfile,'rb')
177	PhaseNames = []
178	while True:
179	try:
180	data = cPickle.load(fl)
181	except EOFError:
182	break
183	datum = data[0]
184	if 'Phases' == datum[0]:
185	for datus in data[1:]:
186	PhaseNames.append(datus[0])
187	fl.close()
188	return PhaseNames
189
190	def GetAllPhaseData(GPXfile,PhaseName):
191	''' Returns the entire dictionary for PhaseName from GSASII gpx file
192	input:
193	GPXfile = gpx full file name
194	PhaseName = phase name
195	return:
196	phase dictionary
197	'''
198	fl = open(GPXfile,'rb')
199	General = {}
200	Atoms = []
201	while True:
202	try:
203	data = cPickle.load(fl)
204	except EOFError:
205	break
206	datum = data[0]
207	if 'Phases' == datum[0]:
208	for datus in data[1:]:
209	if datus[0] == PhaseName:
210	break
211	fl.close()
212	return datus[1]
213
214	def GetHistograms(GPXfile,hNames):
215	""" Returns a dictionary of histograms found in GSASII gpx file
216	input:
217	GPXfile = .gpx full file name
218	hNames = list of histogram names
219	return:
220	Histograms = dictionary of histograms (types = PWDR & HKLF)
221	"""
222	fl = open(GPXfile,'rb')
223	Histograms = {}
224	while True:
225	try:
226	data = cPickle.load(fl)
227	except EOFError:
228	break
229	datum = data[0]
230	hist = datum[0]
231	if hist in hNames:
232	if 'PWDR' in hist[:4]:
233	PWDRdata = {}
234	try:
235	PWDRdata.update(datum[1][0]) #weight factor
236	except ValueError:
237	PWDRdata['wtFactor'] = 1.0 #patch
238	PWDRdata['Data'] = datum[1][1] #powder data arrays
239	PWDRdata[data[2][0]] = data[2][1] #Limits
240	PWDRdata[data[3][0]] = data[3][1] #Background
241	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
242	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
243	try:
244	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
245	except IndexError:
246	PWDRdata['Reflection Lists'] = {}
247
248	Histograms[hist] = PWDRdata
249	elif 'HKLF' in hist[:4]:
250	HKLFdata = {}
251	try:
252	HKLFdata.update(datum[1][0]) #weight factor
253	except ValueError:
254	HKLFdata['wtFactor'] = 1.0 #patch
255	HKLFdata['Data'] = datum[1][1]
256	HKLFdata[data[1][0]] = data[1][1] #Instrument parameters
257	HKLFdata['Reflection Lists'] = None
258	Histograms[hist] = HKLFdata
259	fl.close()
260	return Histograms
261
262	def GetHistogramNames(GPXfile,hType):
263	""" Returns a list of histogram names found in GSASII gpx file
264	input:
265	GPXfile = .gpx full file name
266	hType = list ['PWDR','HKLF']
267	return:
268	HistogramNames = list of histogram names (types = PWDR & HKLF)
269	"""
270	fl = open(GPXfile,'rb')
271	HistogramNames = []
272	while True:
273	try:
274	data = cPickle.load(fl)
275	except EOFError:
276	break
277	datum = data[0]
278	if datum[0][:4] in hType:
279	HistogramNames.append(datum[0])
280	fl.close()
281	return HistogramNames
282
283	def GetUsedHistogramsAndPhases(GPXfile):
284	''' Returns all histograms that are found in any phase
285	and any phase that uses a histogram
286	input:
287	GPXfile = .gpx full file name
288	return:
289	Histograms = dictionary of histograms as {name:data,...}
290	Phases = dictionary of phases that use histograms
291	'''
292	phaseNames = GetPhaseNames(GPXfile)
293	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
294	allHistograms = GetHistograms(GPXfile,histoList)
295	phaseData = {}
296	for name in phaseNames:
297	phaseData[name] = GetAllPhaseData(GPXfile,name)
298	Histograms = {}
299	Phases = {}
300	for phase in phaseData:
301	Phase = phaseData[phase]
302	if Phase['Histograms']:
303	if phase not in Phases:
304	pId = phaseNames.index(phase)
305	Phase['pId'] = pId
306	Phases[phase] = Phase
307	for hist in Phase['Histograms']:
308	if 'Use' not in Phase['Histograms'][hist]: #patch
309	Phase['Histograms'][hist]['Use'] = True
310	if hist not in Histograms and Phase['Histograms'][hist]['Use']:
311	Histograms[hist] = allHistograms[hist]
312	hId = histoList.index(hist)
313	Histograms[hist]['hId'] = hId
314	return Histograms,Phases
315
316	def getBackupName(GPXfile,makeBack):
317	GPXpath,GPXname = ospath.split(GPXfile)
318	if GPXpath == '': GPXpath = '.'
319	Name = ospath.splitext(GPXname)[0]
320	files = os.listdir(GPXpath)
321	last = 0
322	for name in files:
323	name = name.split('.')
324	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
325	if makeBack:
326	last = max(last,int(name[1].strip('bak'))+1)
327	else:
328	last = max(last,int(name[1].strip('bak')))
329	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
330	return GPXback
331
332	def GPXBackup(GPXfile,makeBack=True):
333	import distutils.file_util as dfu
334	GPXback = getBackupName(GPXfile,makeBack)
335	dfu.copy_file(GPXfile,GPXback)
336	return GPXback
337
338	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
339	''' Updates gpxfile from all histograms that are found in any phase
340	and any phase that used a histogram
341	input:
342	GPXfile = .gpx full file name
343	Histograms = dictionary of histograms as {name:data,...}
344	Phases = dictionary of phases that use histograms
345	CovData = dictionary of refined variables, varyList, & covariance matrix
346	makeBack = True if new backup of .gpx file is to be made; else use the last one made
347	'''
348
349	GPXback = GPXBackup(GPXfile,makeBack)
350	print 'Read from file:',GPXback
351	print 'Save to file :',GPXfile
352	infile = open(GPXback,'rb')
353	outfile = open(GPXfile,'wb')
354	while True:
355	try:
356	data = cPickle.load(infile)
357	except EOFError:
358	break
359	datum = data[0]
360	# print 'read: ',datum[0]
361	if datum[0] == 'Phases':
362	for iphase in range(len(data)):
363	if data[iphase][0] in Phases:
364	phaseName = data[iphase][0]
365	data[iphase][1].update(Phases[phaseName])
366	elif datum[0] == 'Covariance':
367	data[0][1] = CovData
368	try:
369	histogram = Histograms[datum[0]]
370	# print 'found ',datum[0]
371	data[0][1][1] = histogram['Data']
372	for datus in data[1:]:
373	# print ' read: ',datus[0]
374	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
375	datus[1] = histogram[datus[0]]
376	except KeyError:
377	pass
378
379	cPickle.dump(data,outfile,1)
380	infile.close()
381	outfile.close()
382	print 'GPX file save successful'
383
384	def SetSeqResult(GPXfile,Histograms,SeqResult):
385	GPXback = GPXBackup(GPXfile)
386	print 'Read from file:',GPXback
387	print 'Save to file :',GPXfile
388	infile = open(GPXback,'rb')
389	outfile = open(GPXfile,'wb')
390	while True:
391	try:
392	data = cPickle.load(infile)
393	except EOFError:
394	break
395	datum = data[0]
396	if datum[0] == 'Sequental results':
397	data[0][1] = SeqResult
398	try:
399	histogram = Histograms[datum[0]]
400	data[0][1][1] = histogram['Data']
401	for datus in data[1:]:
402	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
403	datus[1] = histogram[datus[0]]
404	except KeyError:
405	pass
406
407	cPickle.dump(data,outfile,1)
408	infile.close()
409	outfile.close()
410	print 'GPX file save successful'
411
412	def ShowBanner(pFile=None):
413	print >>pFile,80''
414	print >>pFile,' General Structure Analysis System-II Crystal Structure Refinement'
415	print >>pFile,' by Robert B. Von Dreele & Brian H. Toby'
416	print >>pFile,' Argonne National Laboratory(C), 2010'
417	print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,'
418	print >>pFile,' as Operator of Argonne National Laboratory.'
419	print >>pFile,80'','\n'
420
421	def ShowControls(Controls,pFile=None):
422	print >>pFile,' Least squares controls:'
423	print >>pFile,' Refinement type: ',Controls['deriv type']
424	if 'Hessian' in Controls['deriv type']:
425	print >>pFile,' Maximum number of cycles:',Controls['max cyc']
426	else:
427	print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
428	print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
429
430	def GetFFtable(General):
431	''' returns a dictionary of form factor data for atom types found in General
432	input:
433	General = dictionary of phase info.; includes AtomTypes
434	return:
435	FFtable = dictionary of form factor data; key is atom type
436	'''
437	atomTypes = General['AtomTypes']
438	FFtable = {}
439	for El in atomTypes:
440	FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
441	for item in FFs:
442	if item['Symbol'] == El.upper():
443	FFtable[El] = item
444	return FFtable
445
446	def GetBLtable(General):
447	''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
448	input:
449	General = dictionary of phase info.; includes AtomTypes & Isotopes
450	return:
451	BLtable = dictionary of scattering length data; key is atom type
452	'''
453	atomTypes = General['AtomTypes']
454	BLtable = {}
455	isotopes = General['Isotopes']
456	isotope = General['Isotope']
457	for El in atomTypes:
458	BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
459	return BLtable
460
461	def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
462	# if SGLaue in ['-1','2/m','mmm']:
463	# return #no Pawley symmetry required constraints
464	eqvDict = {}
465	for i,varyI in enumerate(pawleyVary):
466	eqvDict[varyI] = []
467	refI = int(varyI.split(':')[-1])
468	ih,ik,il = PawleyRef[refI][:3]
469	dspI = PawleyRef[refI][4]
470	for varyJ in pawleyVary[i+1:]:
471	refJ = int(varyJ.split(':')[-1])
472	jh,jk,jl = PawleyRef[refJ][:3]
473	dspJ = PawleyRef[refJ][4]
474	if SGLaue in ['4/m','4/mmm']:
475	isum = ih2+ik2
476	jsum = jh2+jk2
477	if abs(il) == abs(jl) and isum == jsum:
478	eqvDict[varyI].append(varyJ)
479	elif SGLaue in ['3R','3mR']:
480	isum = ih2+ik2+il**2
481	jsum = jh2+jk2jl*2
482	isum2 = ihik+ihil+ik*il
483	jsum2 = jhjk+jhjl+jk*jl
484	if isum == jsum and isum2 == jsum2:
485	eqvDict[varyI].append(varyJ)
486	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
487	isum = ih2+ik2+ih*ik
488	jsum = jh2+jk2+jh*jk
489	if abs(il) == abs(jl) and isum == jsum:
490	eqvDict[varyI].append(varyJ)
491	elif SGLaue in ['m3','m3m']:
492	isum = ih2+ik2+il**2
493	jsum = jh2+jk2+jl**2
494	if isum == jsum:
495	eqvDict[varyI].append(varyJ)
496	elif abs(dspI-dspJ)/dspI < 1.e-4:
497	eqvDict[varyI].append(varyJ)
498	for item in pawleyVary:
499	if eqvDict[item]:
500	for item2 in pawleyVary:
501	if item2 in eqvDict[item]:
502	eqvDict[item2] = []
503	G2mv.StoreEquivalence(item,eqvDict[item])
504
505	def cellVary(pfx,SGData):
506	if SGData['SGLaue'] in ['-1',]:
507	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
508	elif SGData['SGLaue'] in ['2/m',]:
509	if SGData['SGUniq'] == 'a':
510	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
511	elif SGData['SGUniq'] == 'b':
512	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
513	else:
514	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
515	elif SGData['SGLaue'] in ['mmm',]:
516	return [pfx+'A0',pfx+'A1',pfx+'A2']
517	elif SGData['SGLaue'] in ['4/m','4/mmm']:
518	return [pfx+'A0',pfx+'A2']
519	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
520	return [pfx+'A0',pfx+'A2']
521	elif SGData['SGLaue'] in ['3R', '3mR']:
522	return [pfx+'A0',pfx+'A3']
523	elif SGData['SGLaue'] in ['m3m','m3']:
524	return [pfx+'A0',]
525
526	################################################################################
527	##### Phase data
528	################################################################################
529
530	def GetPhaseData(PhaseData,RestraintDict={},Print=True,pFile=None):
531
532	def PrintFFtable(FFtable):
533	print >>pFile,'\n X-ray scattering factors:'
534	print >>pFile,' Symbol fa fb fc'
535	print >>pFile,99*'-'
536	for Ename in FFtable:
537	ffdata = FFtable[Ename]
538	fa = ffdata['fa']
539	fb = ffdata['fb']
540	print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \
541	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
542
543	def PrintBLtable(BLtable):
544	print >>pFile,'\n Neutron scattering factors:'
545	print >>pFile,' Symbol isotope mass b resonant terms'
546	print >>pFile,99*'-'
547	for Ename in BLtable:
548	bldata = BLtable[Ename]
549	isotope = bldata[0]
550	mass = bldata[1][0]
551	blen = bldata[1][1]
552	bres = []
553	if len(bldata[1]) > 2:
554	bres = bldata[1][2:]
555	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
556	for item in bres:
557	line += '%10.5g'%(item)
558	print >>pFile,line
559
560	#def PrintRBObjects()
561
562	def PrintAtoms(General,Atoms):
563	cx,ct,cs,cia = General['AtomPtrs']
564	print >>pFile,'\n Atoms:'
565	line = ' name type refine? x y z '+ \
566	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
567	if General['Type'] == 'magnetic':
568	line += ' Mx My Mz'
569	elif General['Type'] == 'macromolecular':
570	line = ' res no residue chain'+line
571	print >>pFile,line
572	if General['Type'] == 'nuclear':
573	print >>pFile,135*'-'
574	for i,at in enumerate(Atoms):
575	line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
576	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
577	if at[cia] == 'I':
578	line += '%8.4f'%(at[cia+1])+48*' '
579	else:
580	line += 8*' '
581	for j in range(6):
582	line += '%8.4f'%(at[cia+1+j])
583	print >>pFile,line
584	elif General['Type'] == 'macromolecular':
585	print >>pFile,135*'-'
586	for i,at in enumerate(Atoms):
587	line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
588	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
589	if at[cia] == 'I':
590	line += '%8.4f'%(at[cia+1])+48*' '
591	else:
592	line += 8*' '
593	for j in range(6):
594	line += '%8.4f'%(at[cia+1+j])
595	print >>pFile,line
596
597	def PrintTexture(textureData):
598	topstr = '\n Spherical harmonics texture: Order:' + \
599	str(textureData['Order'])
600	if textureData['Order']:
601	print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
602	else:
603	print >>pFile,topstr
604	return
605	names = ['omega','chi','phi']
606	line = '\n'
607	for name in names:
608	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
609	print >>pFile,line
610	print >>pFile,'\n Texture coefficients:'
611	ptlbls = ' names :'
612	ptstr = ' values:'
613	SHcoeff = textureData['SH Coeff'][1]
614	for item in SHcoeff:
615	ptlbls += '%12s'%(item)
616	ptstr += '%12.4f'%(SHcoeff[item])
617	print >>pFile,ptlbls
618	print >>pFile,ptstr
619
620	if Print:print >>pFile,' Phases:'
621	phaseVary = []
622	phaseDict = {}
623	phaseConstr = {}
624	pawleyLookup = {}
625	FFtables = {} #scattering factors - xrays
626	BLtables = {} # neutrons
627	Natoms = {}
628	AtMults = {}
629	AtIA = {}
630	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
631	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
632	atomIndx = {}
633	for name in PhaseData:
634	General = PhaseData[name]['General']
635	pId = PhaseData[name]['pId']
636	pfx = str(pId)+'::'
637	FFtable = GetFFtable(General)
638	BLtable = GetBLtable(General)
639	FFtables.update(FFtable)
640	BLtables.update(BLtable)
641	Atoms = PhaseData[name]['Atoms']
642	AtLookup = G2mth.FillAtomLookUp(Atoms)
643	try:
644	PawleyRef = PhaseData[name]['Pawley ref']
645	except KeyError:
646	PawleyRef = []
647	SGData = General['SGData']
648	SGtext = G2spc.SGPrint(SGData)
649	cell = General['Cell']
650	A = G2lat.cell2A(cell[1:7])
651	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
652	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
653	if cell[0]:
654	phaseVary += cellVary(pfx,SGData)
655	#rigid body model input here
656	Natoms[pfx] = 0
657	if Atoms and not General.get('doPawley'):
658	cx,ct,cs,cia = General['AtomPtrs']
659	if General['Type'] in ['nuclear','macromolecular']:
660	Natoms[pfx] = len(Atoms)
661	for i,at in enumerate(Atoms):
662	atomIndx[at[-1]] = [pfx,i] #lookup table for restraints
663	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
664	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
665	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
666	pfx+'AI/A:'+str(i):at[cia],})
667	if at[cia] == 'I':
668	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
669	else:
670	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
671	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
672	if 'F' in at[ct+1]:
673	phaseVary.append(pfx+'Afrac:'+str(i))
674	if 'X' in at[ct+1]:
675	xId,xCoef = G2spc.GetCSxinel(at[cs])
676	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
677	equivs = [[],[],[]]
678	for j in range(3):
679	if xId[j] > 0:
680	phaseVary.append(names[j])
681	equivs[xId[j]-1].append([names[j],xCoef[j]])
682	for equiv in equivs:
683	if len(equiv) > 1:
684	name = equiv[0][0]
685	for eqv in equiv[1:]:
686	G2mv.StoreEquivalence(name,(eqv,))
687	if 'U' in at[ct+1]:
688	if at[9] == 'I':
689	phaseVary.append(pfx+'AUiso:'+str(i))
690	else:
691	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
692	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
693	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
694	equivs = [[],[],[],[],[],[]]
695	for j in range(6):
696	if uId[j] > 0:
697	phaseVary.append(names[j])
698	equivs[uId[j]-1].append([names[j],uCoef[j]])
699	for equiv in equivs:
700	if len(equiv) > 1:
701	name = equiv[0][0]
702	for eqv in equiv[1:]:
703	G2mv.StoreEquivalence(name,(eqv,))
704	# elif General['Type'] == 'magnetic':
705	# elif General['Type'] == 'macromolecular':
706	textureData = General['SH Texture']
707	if textureData['Order']:
708	phaseDict[pfx+'SHorder'] = textureData['Order']
709	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
710	for item in ['omega','chi','phi']:
711	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
712	if textureData['Sample '+item][0]:
713	phaseVary.append(pfx+'SH '+item)
714	for item in textureData['SH Coeff'][1]:
715	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
716	if textureData['SH Coeff'][0]:
717	phaseVary.append(pfx+item)
718
719	if Print:
720	print >>pFile,'\n Phase name: ',General['Name']
721	print >>pFile,135*'-'
722	PrintFFtable(FFtable)
723	PrintBLtable(BLtable)
724	print >>pFile,''
725	for line in SGtext: print >>pFile,line
726	#PrintRBObjects(whatever is needed here)
727	PrintAtoms(General,Atoms)
728	print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
729	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
730	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
731	PrintTexture(textureData)
732	if name in RestraintDict:
733	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
734	textureData,RestraintDict[name],pFile)
735
736	elif PawleyRef:
737	pawleyVary = []
738	for i,refl in enumerate(PawleyRef):
739	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
740	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
741	if refl[5]:
742	pawleyVary.append(pfx+'PWLref:'+str(i))
743	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
744	phaseVary += pawleyVary
745
746	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
747
748	def cellFill(pfx,SGData,parmDict,sigDict):
749	if SGData['SGLaue'] in ['-1',]:
750	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
751	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
752	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
753	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
754	elif SGData['SGLaue'] in ['2/m',]:
755	if SGData['SGUniq'] == 'a':
756	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
757	parmDict[pfx+'A3'],0,0]
758	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
759	sigDict[pfx+'A3'],0,0]
760	elif SGData['SGUniq'] == 'b':
761	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
762	0,parmDict[pfx+'A4'],0]
763	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
764	0,sigDict[pfx+'A4'],0]
765	else:
766	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
767	0,0,parmDict[pfx+'A5']]
768	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
769	0,0,sigDict[pfx+'A5']]
770	elif SGData['SGLaue'] in ['mmm',]:
771	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
772	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
773	elif SGData['SGLaue'] in ['4/m','4/mmm']:
774	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
775	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
776	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
777	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
778	parmDict[pfx+'A0'],0,0]
779	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
780	elif SGData['SGLaue'] in ['3R', '3mR']:
781	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
782	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
783	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
784	elif SGData['SGLaue'] in ['m3m','m3']:
785	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
786	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
787	return A,sigA
788
789	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
790	if phaseRest:
791	Amat = G2lat.cell2AB(cell)[0]
792	cx,ct,cs = AtPtrs[:3]
793	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
794	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
795	['ChemComp','Sites'],['Texture','HKLs']]
796	for name,rest in names:
797	itemRest = phaseRest[name]
798	if itemRest[rest] and itemRest['Use']:
799	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
800	if name in ['Bond','Angle','Plane','Chiral']:
801	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
802	for indx,ops,obs,esd in itemRest[rest]:
803	try:
804	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
805	AtName = ''
806	for i,Aname in enumerate(AtNames):
807	AtName += Aname
808	if ops[i] == '1':
809	AtName += '-'
810	else:
811	AtName += '+('+ops[i]+')-'
812	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
813	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
814	if name == 'Bond':
815	calc = G2mth.getRestDist(XYZ,Amat)
816	elif name == 'Angle':
817	calc = G2mth.getRestAngle(XYZ,Amat)
818	elif name == 'Plane':
819	calc = G2mth.getRestPlane(XYZ,Amat)
820	elif name == 'Chiral':
821	calc = G2mth.getRestChiral(XYZ,Amat)
822	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
823	except KeyError:
824	del itemRest[rest]
825	elif name in ['Torsion','Rama']:
826	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
827	coeffDict = itemRest['Coeff']
828	for indx,ops,cofName,esd in enumerate(itemRest[rest]):
829	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
830	AtName = ''
831	for i,Aname in enumerate(AtNames):
832	AtName += Aname+'+('+ops[i]+')-'
833	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
834	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
835	if name == 'Torsion':
836	tor = G2mth.getRestTorsion(XYZ,Amat)
837	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
838	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(AtName[:-1],calc,obs,esd,(obs-calc)/esd,tor)
839	else:
840	phi,psi = G2mth.getRestRama(XYZ,Amat)
841	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
842	print >>pFile,' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s'%(AtName[:-1],calc,obs,esd,(obs-calc)/esd,phi,psi)
843	elif name == 'ChemComp':
844	print >>pFile,' atoms mul*frac factor prod'
845	for indx,factors,obs,esd in itemRest[rest]:
846	try:
847	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
848	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
849	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
850	mulfrac = mul*frac
851	calcs = mulfracfactors
852	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
853	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
854	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
855	except KeyError:
856	del itemRest[rest]
857	elif name == 'Texture' and textureData['Order']:
858	Start = False
859	SHCoef = textureData['SH Coeff'][1]
860	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
861	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
862	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
863	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
864	PH = np.array(hkl)
865	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
866	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
867	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
868	Z = ma.masked_greater(Z,0.0)
869	num = ma.count(Z)
870	sum = 0
871	if num:
872	sum = np.sum(Z)
873	print ' %d %d %d %d %8.3f %8.3f %8d %s %8.3f'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2)
874
875	def getCellEsd(pfx,SGData,A,covData):
876	dpr = 180./np.pi
877	rVsq = G2lat.calc_rVsq(A)
878	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
879	cell = np.array(G2lat.Gmat2cell(g)) #real cell
880	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
881	scos = cosd(cellst[3:6])
882	ssin = sind(cellst[3:6])
883	scot = scos/ssin
884	rcos = cosd(cell[3:6])
885	rsin = sind(cell[3:6])
886	rcot = rcos/rsin
887	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
888	varyList = covData['varyList']
889	covMatrix = covData['covMatrix']
890	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
891	Ax = np.array(A)
892	Ax[3:] /= 2.
893	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
894	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
895	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
896	Vol = 1/np.sqrt(rVsq)
897	sigVol = Vol*3np.sqrt(srcvlsq)/2.
898	R123 = Ax[0]Ax[1]Ax[2]
899	dsasdg = np.zeros((3,6))
900	dadg = np.zeros((6,6))
901	for i0 in range(3): #0 1 2
902	i1 = (i0+1)%3 #1 2 0
903	i2 = (i1+1)%3 #2 0 1
904	i3 = 5-i2 #3 5 4
905	i4 = 5-i1 #4 3 5
906	i5 = 5-i0 #5 4 3
907	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
908	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
909	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
910	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
911	denom = np.sqrt(denmsq)
912	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
913	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
914	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
915	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
916	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
917	dadg[i5][i5] = -Ax[i0]/denom
918	for i0 in range(3):
919	i1 = (i0+1)%3
920	i2 = (i1+1)%3
921	i3 = 5-i2
922	for ij in range(6):
923	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
924	if ij == i0:
925	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
926	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
927	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
928	var = np.diag(sigMat)
929	CS = np.where(var>0.,np.sqrt(var),0.)
930	cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
931	return cellSig
932
933	def SetPhaseData(parmDict,sigDict,Phases,covData,RestraintDict=None,pFile=None):
934
935	def PrintAtomsAndSig(General,Atoms,atomsSig):
936	print >>pFile,'\n Atoms:'
937	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
938	if General['Type'] == 'magnetic':
939	line += ' Mx My Mz'
940	elif General['Type'] == 'macromolecular':
941	line = ' res no residue chain '+line
942	print >>pFile,line
943	if General['Type'] == 'nuclear':
944	print >>pFile,135*'-'
945	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
946	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
947	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
948	for atyp in General['AtomTypes']: #zero composition
949	General['NoAtoms'][atyp] = 0.
950	for i,at in enumerate(Atoms):
951	General['NoAtoms'][at[1]] += at[6]*float(at[8]) #new composition
952	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
953	valstr = ' values:'
954	sigstr = ' sig :'
955	for ind in [3,4,5,6]:
956	sigind = str(i)+':'+str(ind)
957	valstr += fmt[ind]%(at[ind])
958	if sigind in atomsSig:
959	sigstr += fmt[ind]%(atomsSig[sigind])
960	else:
961	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
962	if at[9] == 'I':
963	valstr += fmt[10]%(at[10])
964	if str(i)+':10' in atomsSig:
965	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
966	else:
967	sigstr += 8*' '
968	else:
969	valstr += 8*' '
970	sigstr += 8*' '
971	for ind in [11,12,13,14,15,16]:
972	sigind = str(i)+':'+str(ind)
973	valstr += fmt[ind]%(at[ind])
974	if sigind in atomsSig:
975	sigstr += fmt[ind]%(atomsSig[sigind])
976	else:
977	sigstr += 8*' '
978	print >>pFile,name
979	print >>pFile,valstr
980	print >>pFile,sigstr
981
982	#def PrintRBObjectsAndSig()
983
984	def PrintSHtextureAndSig(textureData,SHtextureSig):
985	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
986	names = ['omega','chi','phi']
987	namstr = ' names :'
988	ptstr = ' values:'
989	sigstr = ' esds :'
990	for name in names:
991	namstr += '%12s'%(name)
992	ptstr += '%12.3f'%(textureData['Sample '+name][1])
993	if 'Sample '+name in SHtextureSig:
994	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
995	else:
996	sigstr += 12*' '
997	print >>pFile,namstr
998	print >>pFile,ptstr
999	print >>pFile,sigstr
1000	print >>pFile,'\n Texture coefficients:'
1001	namstr = ' names :'
1002	ptstr = ' values:'
1003	sigstr = ' esds :'
1004	SHcoeff = textureData['SH Coeff'][1]
1005	for name in SHcoeff:
1006	namstr += '%12s'%(name)
1007	ptstr += '%12.3f'%(SHcoeff[name])
1008	if name in SHtextureSig:
1009	sigstr += '%12.3f'%(SHtextureSig[name])
1010	else:
1011	sigstr += 12*' '
1012	print >>pFile,namstr
1013	print >>pFile,ptstr
1014	print >>pFile,sigstr
1015
1016
1017	print >>pFile,'\n Phases:'
1018	for phase in Phases:
1019	print >>pFile,' Result for phase: ',phase
1020	Phase = Phases[phase]
1021	General = Phase['General']
1022	SGData = General['SGData']
1023	Atoms = Phase['Atoms']
1024	AtLookup = G2mth.FillAtomLookUp(Atoms)
1025	cell = General['Cell']
1026	pId = Phase['pId']
1027	pfx = str(pId)+'::'
1028	if cell[0]:
1029	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1030	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1031	print >>pFile,' Reciprocal metric tensor: '
1032	ptfmt = "%15.9f"
1033	names = ['A11','A22','A33','A12','A13','A23']
1034	namstr = ' names :'
1035	ptstr = ' values:'
1036	sigstr = ' esds :'
1037	for name,a,siga in zip(names,A,sigA):
1038	namstr += '%15s'%(name)
1039	ptstr += ptfmt%(a)
1040	if siga:
1041	sigstr += ptfmt%(siga)
1042	else:
1043	sigstr += 15*' '
1044	print >>pFile,namstr
1045	print >>pFile,ptstr
1046	print >>pFile,sigstr
1047	cell[1:7] = G2lat.A2cell(A)
1048	cell[7] = G2lat.calc_V(A)
1049	print >>pFile,' New unit cell:'
1050	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1051	names = ['a','b','c','alpha','beta','gamma','Volume']
1052	namstr = ' names :'
1053	ptstr = ' values:'
1054	sigstr = ' esds :'
1055	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1056	namstr += '%12s'%(name)
1057	ptstr += fmt%(a)
1058	if siga:
1059	sigstr += fmt%(siga)
1060	else:
1061	sigstr += 12*' '
1062	print >>pFile,namstr
1063	print >>pFile,ptstr
1064	print >>pFile,sigstr
1065
1066	General['Mass'] = 0.
1067	if Phase['General'].get('doPawley'):
1068	pawleyRef = Phase['Pawley ref']
1069	for i,refl in enumerate(pawleyRef):
1070	key = pfx+'PWLref:'+str(i)
1071	refl[6] = parmDict[key]
1072	if key in sigDict:
1073	refl[7] = sigDict[key]
1074	else:
1075	refl[7] = 0
1076	else:
1077	#new RBObj parms here, mod atoms, & PrintRBObjectsAndSig
1078	atomsSig = {}
1079	if General['Type'] == 'nuclear': #this needs macromolecular variant!
1080	for i,at in enumerate(Atoms):
1081	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1082	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1083	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1084	for ind in [3,4,5,6]:
1085	at[ind] = parmDict[names[ind]]
1086	if ind in [3,4,5]:
1087	name = names[ind].replace('A','dA')
1088	else:
1089	name = names[ind]
1090	if name in sigDict:
1091	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1092	if at[9] == 'I':
1093	at[10] = parmDict[names[10]]
1094	if names[10] in sigDict:
1095	atomsSig[str(i)+':10'] = sigDict[names[10]]
1096	else:
1097	for ind in [11,12,13,14,15,16]:
1098	at[ind] = parmDict[names[ind]]
1099	if names[ind] in sigDict:
1100	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1101	ind = General['AtomTypes'].index(at[1])
1102	General['Mass'] += General['AtomMass'][ind]at[6]at[8]
1103	PrintAtomsAndSig(General,Atoms,atomsSig)
1104
1105	textureData = General['SH Texture']
1106	if textureData['Order']:
1107	SHtextureSig = {}
1108	for name in ['omega','chi','phi']:
1109	aname = pfx+'SH '+name
1110	textureData['Sample '+name][1] = parmDict[aname]
1111	if aname in sigDict:
1112	SHtextureSig['Sample '+name] = sigDict[aname]
1113	for name in textureData['SH Coeff'][1]:
1114	aname = pfx+name
1115	textureData['SH Coeff'][1][name] = parmDict[aname]
1116	if aname in sigDict:
1117	SHtextureSig[name] = sigDict[aname]
1118	PrintSHtextureAndSig(textureData,SHtextureSig)
1119	if phase in RestraintDict:
1120	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1121	textureData,RestraintDict[phase],pFile)
1122
1123	################################################################################
1124	##### Histogram & Phase data
1125	################################################################################
1126
1127	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None):
1128
1129	def PrintSize(hapData):
1130	if hapData[0] in ['isotropic','uniaxial']:
1131	line = '\n Size model : %9s'%(hapData[0])
1132	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1133	if hapData[0] == 'uniaxial':
1134	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1135	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1136	print >>pFile,line
1137	else:
1138	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1139	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1140	Snames = ['S11','S22','S33','S12','S13','S23']
1141	ptlbls = ' names :'
1142	ptstr = ' values:'
1143	varstr = ' refine:'
1144	for i,name in enumerate(Snames):
1145	ptlbls += '%12s' % (name)
1146	ptstr += '%12.6f' % (hapData[4][i])
1147	varstr += '%12s' % (str(hapData[5][i]))
1148	print >>pFile,ptlbls
1149	print >>pFile,ptstr
1150	print >>pFile,varstr
1151
1152	def PrintMuStrain(hapData,SGData):
1153	if hapData[0] in ['isotropic','uniaxial']:
1154	line = '\n Mustrain model: %9s'%(hapData[0])
1155	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1156	if hapData[0] == 'uniaxial':
1157	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1158	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1159	print >>pFile,line
1160	else:
1161	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1162	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1163	Snames = G2spc.MustrainNames(SGData)
1164	ptlbls = ' names :'
1165	ptstr = ' values:'
1166	varstr = ' refine:'
1167	for i,name in enumerate(Snames):
1168	ptlbls += '%12s' % (name)
1169	ptstr += '%12.6f' % (hapData[4][i])
1170	varstr += '%12s' % (str(hapData[5][i]))
1171	print >>pFile,ptlbls
1172	print >>pFile,ptstr
1173	print >>pFile,varstr
1174
1175	def PrintHStrain(hapData,SGData):
1176	print >>pFile,'\n Hydrostatic/elastic strain: '
1177	Hsnames = G2spc.HStrainNames(SGData)
1178	ptlbls = ' names :'
1179	ptstr = ' values:'
1180	varstr = ' refine:'
1181	for i,name in enumerate(Hsnames):
1182	ptlbls += '%12s' % (name)
1183	ptstr += '%12.6f' % (hapData[0][i])
1184	varstr += '%12s' % (str(hapData[1][i]))
1185	print >>pFile,ptlbls
1186	print >>pFile,ptstr
1187	print >>pFile,varstr
1188
1189	def PrintSHPO(hapData):
1190	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1191	str(hapData[4])+' Refine? '+str(hapData[2])
1192	ptlbls = ' names :'
1193	ptstr = ' values:'
1194	for item in hapData[5]:
1195	ptlbls += '%12s'%(item)
1196	ptstr += '%12.3f'%(hapData[5][item])
1197	print >>pFile,ptlbls
1198	print >>pFile,ptstr
1199
1200	def PrintBabinet(hapData):
1201	print >>pFile,'\n Babinet form factor modification: '
1202	ptlbls = ' names :'
1203	ptstr = ' values:'
1204	varstr = ' refine:'
1205	for name in ['BabA','BabU']:
1206	ptlbls += '%12s' % (name)
1207	ptstr += '%12.6f' % (hapData[name][0])
1208	varstr += '%12s' % (str(hapData[name][1]))
1209	print >>pFile,ptlbls
1210	print >>pFile,ptstr
1211	print >>pFile,varstr
1212
1213
1214	hapDict = {}
1215	hapVary = []
1216	controlDict = {}
1217	poType = {}
1218	poAxes = {}
1219	spAxes = {}
1220	spType = {}
1221
1222	for phase in Phases:
1223	HistoPhase = Phases[phase]['Histograms']
1224	SGData = Phases[phase]['General']['SGData']
1225	cell = Phases[phase]['General']['Cell'][1:7]
1226	A = G2lat.cell2A(cell)
1227	pId = Phases[phase]['pId']
1228	histoList = HistoPhase.keys()
1229	histoList.sort()
1230	for histogram in histoList:
1231	try:
1232	Histogram = Histograms[histogram]
1233	except KeyError:
1234	#skip if histogram not included e.g. in a sequential refinement
1235	continue
1236	hapData = HistoPhase[histogram]
1237	hId = Histogram['hId']
1238	if 'PWDR' in histogram:
1239	limits = Histogram['Limits'][1]
1240	inst = Histogram['Instrument Parameters'][0]
1241	Zero = inst['Zero'][1]
1242	if 'C' in inst['Type'][1]:
1243	try:
1244	wave = inst['Lam'][1]
1245	except KeyError:
1246	wave = inst['Lam1'][1]
1247	dmin = wave/(2.0*sind(limits[1]/2.0))
1248	pfx = str(pId)+':'+str(hId)+':'
1249	for item in ['Scale','Extinction']:
1250	hapDict[pfx+item] = hapData[item][0]
1251	if hapData[item][1]:
1252	hapVary.append(pfx+item)
1253	names = G2spc.HStrainNames(SGData)
1254	for i,name in enumerate(names):
1255	hapDict[pfx+name] = hapData['HStrain'][0][i]
1256	if hapData['HStrain'][1][i]:
1257	hapVary.append(pfx+name)
1258	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1259	if hapData['Pref.Ori.'][0] == 'MD':
1260	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1261	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1262	if hapData['Pref.Ori.'][2]:
1263	hapVary.append(pfx+'MD')
1264	else: #'SH' spherical harmonics
1265	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1266	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1267	for item in hapData['Pref.Ori.'][5]:
1268	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1269	if hapData['Pref.Ori.'][2]:
1270	hapVary.append(pfx+item)
1271	for item in ['Mustrain','Size']:
1272	controlDict[pfx+item+'Type'] = hapData[item][0]
1273	hapDict[pfx+item+';mx'] = hapData[item][1][2]
1274	if hapData[item][2][2]:
1275	hapVary.append(pfx+item+';mx')
1276	if hapData[item][0] in ['isotropic','uniaxial']:
1277	hapDict[pfx+item+';i'] = hapData[item][1][0]
1278	if hapData[item][2][0]:
1279	hapVary.append(pfx+item+';i')
1280	if hapData[item][0] == 'uniaxial':
1281	controlDict[pfx+item+'Axis'] = hapData[item][3]
1282	hapDict[pfx+item+';a'] = hapData[item][1][1]
1283	if hapData[item][2][1]:
1284	hapVary.append(pfx+item+';a')
1285	else: #generalized for mustrain or ellipsoidal for size
1286	Nterms = len(hapData[item][4])
1287	if item == 'Mustrain':
1288	names = G2spc.MustrainNames(SGData)
1289	pwrs = []
1290	for name in names:
1291	h,k,l = name[1:]
1292	pwrs.append([int(h),int(k),int(l)])
1293	controlDict[pfx+'MuPwrs'] = pwrs
1294	for i in range(Nterms):
1295	sfx = ':'+str(i)
1296	hapDict[pfx+item+sfx] = hapData[item][4][i]
1297	if hapData[item][5][i]:
1298	hapVary.append(pfx+item+sfx)
1299	for bab in ['BabA','BabU']:
1300	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1301	if hapData['Babinet'][bab][1]:
1302	hapVary.append(pfx+bab)
1303
1304	if Print:
1305	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1306	print >>pFile,135*'-'
1307	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1308	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1309	if hapData['Pref.Ori.'][0] == 'MD':
1310	Ax = hapData['Pref.Ori.'][3]
1311	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1312	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1313	else: #'SH' for spherical harmonics
1314	PrintSHPO(hapData['Pref.Ori.'])
1315	PrintSize(hapData['Size'])
1316	PrintMuStrain(hapData['Mustrain'],SGData)
1317	PrintHStrain(hapData['HStrain'],SGData)
1318	if hapData['Babinet']['BabA'][0]:
1319	PrintBabinet(hapData['Babinet'])
1320	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1321	refList = []
1322	for h,k,l,d in HKLd:
1323	ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1324	mul *= 2 # for powder overlap of Friedel pairs
1325	if ext:
1326	continue
1327	if 'C' in inst['Type'][0]:
1328	pos = 2.0asind(wave/(2.0d))+Zero
1329	if limits[0] < pos < limits[1]:
1330	refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1331	#last item should contain form factor values by atom type
1332	else:
1333	raise ValueError
1334	Histogram['Reflection Lists'][phase] = refList
1335	elif 'HKLF' in histogram:
1336	inst = Histogram['Instrument Parameters'][0]
1337	hId = Histogram['hId']
1338	hfx = ':%d:'%(hId)
1339	for item in inst:
1340	hapDict[hfx+item] = inst[item][1]
1341	pfx = str(pId)+':'+str(hId)+':'
1342	hapDict[pfx+'Scale'] = hapData['Scale'][0]
1343	if hapData['Scale'][1]:
1344	hapVary.append(pfx+'Scale')
1345
1346	extApprox,extType,extParms = hapData['Extinction']
1347	controlDict[pfx+'EType'] = extType
1348	controlDict[pfx+'EApprox'] = extApprox
1349	controlDict[pfx+'Tbar'] = extParms['Tbar']
1350	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1351	if 'Primary' in extType:
1352	Ekey = ['Ep',]
1353	elif 'I & II' in extType:
1354	Ekey = ['Eg','Es']
1355	elif 'Secondary Type II' == extType:
1356	Ekey = ['Es',]
1357	elif 'Secondary Type I' == extType:
1358	Ekey = ['Eg',]
1359	else: #'None'
1360	Ekey = []
1361	for eKey in Ekey:
1362	hapDict[pfx+eKey] = extParms[eKey][0]
1363	if extParms[eKey][1]:
1364	hapVary.append(pfx+eKey)
1365	for bab in ['BabA','BabU']:
1366	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1367	if hapData['Babinet'][bab][1]:
1368	hapVary.append(pfx+bab)
1369	if Print:
1370	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1371	print >>pFile,135*'-'
1372	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1373	if extType != 'None':
1374	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1375	text = ' Parameters :'
1376	for eKey in Ekey:
1377	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1378	print >>pFile,text
1379	if hapData['Babinet']['BabA'][0]:
1380	PrintBabinet(hapData['Babinet'])
1381	Histogram['Reflection Lists'] = phase
1382
1383	return hapVary,hapDict,controlDict
1384
1385	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1386
1387	def PrintSizeAndSig(hapData,sizeSig):
1388	line = '\n Size model: %9s'%(hapData[0])
1389	refine = False
1390	if hapData[0] in ['isotropic','uniaxial']:
1391	line += ' equatorial:%12.4f'%(hapData[1][0])
1392	if sizeSig[0][0]:
1393	line += ', sig:%8.4f'%(sizeSig[0][0])
1394	refine = True
1395	if hapData[0] == 'uniaxial':
1396	line += ' axial:%12.4f'%(hapData[1][1])
1397	if sizeSig[0][1]:
1398	refine = True
1399	line += ', sig:%8.4f'%(sizeSig[0][1])
1400	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1401	if sizeSig[0][2]:
1402	refine = True
1403	line += ', sig:%8.4f'%(sizeSig[0][2])
1404	if refine:
1405	print >>pFile,line
1406	else:
1407	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1408	if sizeSig[0][2]:
1409	refine = True
1410	line += ', sig:%8.4f'%(sizeSig[0][2])
1411	Snames = ['S11','S22','S33','S12','S13','S23']
1412	ptlbls = ' name :'
1413	ptstr = ' value :'
1414	sigstr = ' sig :'
1415	for i,name in enumerate(Snames):
1416	ptlbls += '%12s' % (name)
1417	ptstr += '%12.6f' % (hapData[4][i])
1418	if sizeSig[1][i]:
1419	refine = True
1420	sigstr += '%12.6f' % (sizeSig[1][i])
1421	else:
1422	sigstr += 12*' '
1423	if refine:
1424	print >>pFile,line
1425	print >>pFile,ptlbls
1426	print >>pFile,ptstr
1427	print >>pFile,sigstr
1428
1429	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1430	line = '\n Mustrain model: %9s'%(hapData[0])
1431	refine = False
1432	if hapData[0] in ['isotropic','uniaxial']:
1433	line += ' equatorial:%12.1f'%(hapData[1][0])
1434	if mustrainSig[0][0]:
1435	line += ', sig:%8.1f'%(mustrainSig[0][0])
1436	refine = True
1437	if hapData[0] == 'uniaxial':
1438	line += ' axial:%12.1f'%(hapData[1][1])
1439	if mustrainSig[0][1]:
1440	line += ', sig:%8.1f'%(mustrainSig[0][1])
1441	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1442	if mustrainSig[0][2]:
1443	refine = True
1444	line += ', sig:%8.3f'%(mustrainSig[0][2])
1445	if refine:
1446	print >>pFile,line
1447	else:
1448	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1449	if mustrainSig[0][2]:
1450	refine = True
1451	line += ', sig:%8.3f'%(mustrainSig[0][2])
1452	Snames = G2spc.MustrainNames(SGData)
1453	ptlbls = ' name :'
1454	ptstr = ' value :'
1455	sigstr = ' sig :'
1456	for i,name in enumerate(Snames):
1457	ptlbls += '%12s' % (name)
1458	ptstr += '%12.6f' % (hapData[4][i])
1459	if mustrainSig[1][i]:
1460	refine = True
1461	sigstr += '%12.6f' % (mustrainSig[1][i])
1462	else:
1463	sigstr += 12*' '
1464	if refine:
1465	print >>pFile,line
1466	print >>pFile,ptlbls
1467	print >>pFile,ptstr
1468	print >>pFile,sigstr
1469
1470	def PrintHStrainAndSig(hapData,strainSig,SGData):
1471	Hsnames = G2spc.HStrainNames(SGData)
1472	ptlbls = ' name :'
1473	ptstr = ' value :'
1474	sigstr = ' sig :'
1475	refine = False
1476	for i,name in enumerate(Hsnames):
1477	ptlbls += '%12s' % (name)
1478	ptstr += '%12.6g' % (hapData[0][i])
1479	if name in strainSig:
1480	refine = True
1481	sigstr += '%12.6g' % (strainSig[name])
1482	else:
1483	sigstr += 12*' '
1484	if refine:
1485	print >>pFile,'\n Hydrostatic/elastic strain: '
1486	print >>pFile,ptlbls
1487	print >>pFile,ptstr
1488	print >>pFile,sigstr
1489
1490	def PrintSHPOAndSig(pfx,hapData,POsig):
1491	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1492	ptlbls = ' names :'
1493	ptstr = ' values:'
1494	sigstr = ' sig :'
1495	for item in hapData[5]:
1496	ptlbls += '%12s'%(item)
1497	ptstr += '%12.3f'%(hapData[5][item])
1498	if pfx+item in POsig:
1499	sigstr += '%12.3f'%(POsig[pfx+item])
1500	else:
1501	sigstr += 12*' '
1502	print >>pFile,ptlbls
1503	print >>pFile,ptstr
1504	print >>pFile,sigstr
1505
1506	def PrintExtAndSig(pfx,hapData,ScalExtSig):
1507	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
1508	text = ''
1509	for item in hapData[2]:
1510	if pfx+item in ScalExtSig:
1511	text += ' %s: '%(item)
1512	text += '%12.2e'%(hapData[2][item][0])
1513	if pfx+item in ScalExtSig:
1514	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
1515	print >>pFile,text
1516
1517	def PrintBabinetAndSig(pfx,hapData,BabSig):
1518	print >>pFile,'\n Babinet form factor modification: '
1519	ptlbls = ' names :'
1520	ptstr = ' values:'
1521	sigstr = ' sig :'
1522	for item in hapData:
1523	ptlbls += '%12s'%(item)
1524	ptstr += '%12.3f'%(hapData[item][0])
1525	if pfx+item in BabSig:
1526	sigstr += '%12.3f'%(BabSig[pfx+item])
1527	else:
1528	sigstr += 12*' '
1529	print >>pFile,ptlbls
1530	print >>pFile,ptstr
1531	print >>pFile,sigstr
1532
1533	PhFrExtPOSig = {}
1534	SizeMuStrSig = {}
1535	ScalExtSig = {}
1536	BabSig = {}
1537	wtFrSum = {}
1538	for phase in Phases:
1539	HistoPhase = Phases[phase]['Histograms']
1540	General = Phases[phase]['General']
1541	SGData = General['SGData']
1542	pId = Phases[phase]['pId']
1543	histoList = HistoPhase.keys()
1544	histoList.sort()
1545	for histogram in histoList:
1546	try:
1547	Histogram = Histograms[histogram]
1548	except KeyError:
1549	#skip if histogram not included e.g. in a sequential refinement
1550	continue
1551	hapData = HistoPhase[histogram]
1552	hId = Histogram['hId']
1553	pfx = str(pId)+':'+str(hId)+':'
1554	if hId not in wtFrSum:
1555	wtFrSum[hId] = 0.
1556	if 'PWDR' in histogram:
1557	for item in ['Scale','Extinction']:
1558	hapData[item][0] = parmDict[pfx+item]
1559	if pfx+item in sigDict:
1560	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1561	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
1562	if hapData['Pref.Ori.'][0] == 'MD':
1563	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1564	if pfx+'MD' in sigDict:
1565	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
1566	else: #'SH' spherical harmonics
1567	for item in hapData['Pref.Ori.'][5]:
1568	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1569	if pfx+item in sigDict:
1570	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1571	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1572	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
1573	pfx+'HStrain':{}})
1574	for item in ['Mustrain','Size']:
1575	hapData[item][1][2] = parmDict[pfx+item+';mx']
1576	hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
1577	if pfx+item+';mx' in sigDict:
1578	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
1579	if hapData[item][0] in ['isotropic','uniaxial']:
1580	hapData[item][1][0] = parmDict[pfx+item+';i']
1581	if item == 'Size':
1582	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1583	if pfx+item+';i' in sigDict:
1584	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
1585	if hapData[item][0] == 'uniaxial':
1586	hapData[item][1][1] = parmDict[pfx+item+';a']
1587	if item == 'Size':
1588	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
1589	if pfx+item+';a' in sigDict:
1590	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
1591	else: #generalized for mustrain or ellipsoidal for size
1592	Nterms = len(hapData[item][4])
1593	for i in range(Nterms):
1594	sfx = ':'+str(i)
1595	hapData[item][4][i] = parmDict[pfx+item+sfx]
1596	if pfx+item+sfx in sigDict:
1597	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
1598	names = G2spc.HStrainNames(SGData)
1599	for i,name in enumerate(names):
1600	hapData['HStrain'][0][i] = parmDict[pfx+name]
1601	if pfx+name in sigDict:
1602	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
1603	for name in ['BabA','BabU']:
1604	hapData['Babinet'][name][0] = parmDict[pfx+name]
1605	if pfx+name in sigDict:
1606	BabSig[pfx+name] = sigDict[pfx+name]
1607
1608	elif 'HKLF' in histogram:
1609	for item in ['Scale',]:
1610	if parmDict.get(pfx+item):
1611	hapData[item][0] = parmDict[pfx+item]
1612	if pfx+item in sigDict:
1613	ScalExtSig[pfx+item] = sigDict[pfx+item]
1614	for item in ['Ep','Eg','Es']:
1615	if parmDict.get(pfx+item):
1616	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
1617	if pfx+item in sigDict:
1618	ScalExtSig[pfx+item] = sigDict[pfx+item]
1619	for name in ['BabA','BabU']:
1620	hapData['Babinet'][name][0] = parmDict[pfx+name]
1621	if pfx+name in sigDict:
1622	BabSig[pfx+name] = sigDict[pfx+name]
1623
1624	if Print:
1625	for phase in Phases:
1626	HistoPhase = Phases[phase]['Histograms']
1627	General = Phases[phase]['General']
1628	SGData = General['SGData']
1629	pId = Phases[phase]['pId']
1630	histoList = HistoPhase.keys()
1631	histoList.sort()
1632	for histogram in histoList:
1633	try:
1634	Histogram = Histograms[histogram]
1635	except KeyError:
1636	#skip if histogram not included e.g. in a sequential refinement
1637	continue
1638	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1639	print >>pFile,130*'-'
1640	hapData = HistoPhase[histogram]
1641	hId = Histogram['hId']
1642	pfx = str(pId)+':'+str(hId)+':'
1643	if 'PWDR' in histogram:
1644	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
1645	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1646
1647	if pfx+'Scale' in PhFrExtPOSig:
1648	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
1649	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
1650	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
1651	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
1652	if pfx+'Extinction' in PhFrExtPOSig:
1653	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
1654	if hapData['Pref.Ori.'][0] == 'MD':
1655	if pfx+'MD' in PhFrExtPOSig:
1656	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
1657	else:
1658	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
1659	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
1660	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
1661	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
1662	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
1663
1664	elif 'HKLF' in histogram:
1665	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
1666	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1667	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
1668	if pfx+'Scale' in ScalExtSig:
1669	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
1670	if hapData['Extinction'][0] != 'None':
1671	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
1672	if len(BabSig):
1673	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
1674
1675	################################################################################
1676	##### Histogram data
1677	################################################################################
1678
1679	def GetHistogramData(Histograms,Print=True,pFile=None):
1680
1681	def GetBackgroundParms(hId,Background):
1682	Back = Background[0]
1683	DebyePeaks = Background[1]
1684	bakType,bakFlag = Back[:2]
1685	backVals = Back[3:]
1686	backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1687	backDict = dict(zip(backNames,backVals))
1688	backVary = []
1689	if bakFlag:
1690	backVary = backNames
1691	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
1692	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
1693	debyeDict = {}
1694	debyeList = []
1695	for i in range(DebyePeaks['nDebye']):
1696	debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
1697	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
1698	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
1699	debyeVary = []
1700	for item in debyeList:
1701	if item[1]:
1702	debyeVary.append(item[0])
1703	backDict.update(debyeDict)
1704	backVary += debyeVary
1705	peakDict = {}
1706	peakList = []
1707	for i in range(DebyePeaks['nPeaks']):
1708	peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
1709	':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
1710	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
1711	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
1712	peakVary = []
1713	for item in peakList:
1714	if item[1]:
1715	peakVary.append(item[0])
1716	backDict.update(peakDict)
1717	backVary += peakVary
1718	return bakType,backDict,backVary
1719
1720	def GetInstParms(hId,Inst):
1721	dataType = Inst['Type'][0]
1722	instDict = {}
1723	insVary = []
1724	pfx = ':'+str(hId)+':'
1725	insKeys = Inst.keys()
1726	insKeys.sort()
1727	for item in insKeys:
1728	insName = pfx+item
1729	instDict[insName] = Inst[item][1]
1730	if Inst[item][2]:
1731	insVary.append(insName)
1732	# instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
1733	# instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
1734	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1735	return dataType,instDict,insVary
1736
1737	def GetSampleParms(hId,Sample):
1738	sampVary = []
1739	hfx = ':'+str(hId)+':'
1740	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1741	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1742	Type = Sample['Type']
1743	if 'Bragg' in Type: #Bragg-Brentano
1744	for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering?
1745	sampDict[hfx+item] = Sample[item][0]
1746	if Sample[item][1]:
1747	sampVary.append(hfx+item)
1748	elif 'Debye' in Type: #Debye-Scherrer
1749	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1750	sampDict[hfx+item] = Sample[item][0]
1751	if Sample[item][1]:
1752	sampVary.append(hfx+item)
1753	return Type,sampDict,sampVary
1754
1755	def PrintBackground(Background):
1756	Back = Background[0]
1757	DebyePeaks = Background[1]
1758	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
1759	line = ' Coefficients: '
1760	for i,back in enumerate(Back[3:]):
1761	line += '%10.3f'%(back)
1762	if i and not i%10:
1763	line += '\n'+15*' '
1764	print >>pFile,line
1765	if DebyePeaks['nDebye']:
1766	print >>pFile,'\n Debye diffuse scattering coefficients'
1767	parms = ['DebyeA','DebyeR','DebyeU']
1768	line = ' names : '
1769	for parm in parms:
1770	line += '%8s refine?'%(parm)
1771	print >>pFile,line
1772	for j,term in enumerate(DebyePeaks['debyeTerms']):
1773	line = ' term'+'%2d'%(j)+':'
1774	for i in range(3):
1775	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
1776	print >>pFile,line
1777	if DebyePeaks['nPeaks']:
1778	print >>pFile,'\n Single peak coefficients'
1779	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1780	line = ' names : '
1781	for parm in parms:
1782	line += '%8s refine?'%(parm)
1783	print >>pFile,line
1784	for j,term in enumerate(DebyePeaks['peaksList']):
1785	line = ' peak'+'%2d'%(j)+':'
1786	for i in range(4):
1787	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
1788	print >>pFile,line
1789
1790	def PrintInstParms(Inst):
1791	print >>pFile,'\n Instrument Parameters:'
1792	ptlbls = ' name :'
1793	ptstr = ' value :'
1794	varstr = ' refine:'
1795	insKeys = Inst.keys()
1796	insKeys.sort()
1797	for item in insKeys:
1798	if item != 'Type':
1799	ptlbls += '%12s' % (item)
1800	ptstr += '%12.6f' % (Inst[item][1])
1801	if item in ['Lam1','Lam2','Azimuth']:
1802	varstr += 12*' '
1803	else:
1804	varstr += '%12s' % (str(bool(Inst[item][2])))
1805	print >>pFile,ptlbls
1806	print >>pFile,ptstr
1807	print >>pFile,varstr
1808
1809	def PrintSampleParms(Sample):
1810	print >>pFile,'\n Sample Parameters:'
1811	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1812	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
1813	ptlbls = ' name :'
1814	ptstr = ' value :'
1815	varstr = ' refine:'
1816	if 'Bragg' in Sample['Type']:
1817	for item in ['Scale','Shift','Transparency']:
1818	ptlbls += '%14s'%(item)
1819	ptstr += '%14.4f'%(Sample[item][0])
1820	varstr += '%14s'%(str(bool(Sample[item][1])))
1821
1822	elif 'Debye' in Type: #Debye-Scherrer
1823	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1824	ptlbls += '%14s'%(item)
1825	ptstr += '%14.4f'%(Sample[item][0])
1826	varstr += '%14s'%(str(bool(Sample[item][1])))
1827
1828	print >>pFile,ptlbls
1829	print >>pFile,ptstr
1830	print >>pFile,varstr
1831
1832	histDict = {}
1833	histVary = []
1834	controlDict = {}
1835	histoList = Histograms.keys()
1836	histoList.sort()
1837	for histogram in histoList:
1838	if 'PWDR' in histogram:
1839	Histogram = Histograms[histogram]
1840	hId = Histogram['hId']
1841	pfx = ':'+str(hId)+':'
1842	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
1843	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1844
1845	Background = Histogram['Background']
1846	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1847	controlDict[pfx+'bakType'] = Type
1848	histDict.update(bakDict)
1849	histVary += bakVary
1850
1851	Inst = Histogram['Instrument Parameters'][0]
1852	Type,instDict,insVary = GetInstParms(hId,Inst)
1853	controlDict[pfx+'histType'] = Type
1854	if pfx+'Lam1' in instDict:
1855	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1856	else:
1857	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
1858	histDict.update(instDict)
1859	histVary += insVary
1860
1861	Sample = Histogram['Sample Parameters']
1862	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1863	controlDict[pfx+'instType'] = Type
1864	histDict.update(sampDict)
1865	histVary += sampVary
1866
1867
1868	if Print:
1869	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
1870	print >>pFile,135*'-'
1871	Units = {'C':' deg','T':' msec'}
1872	units = Units[controlDict[pfx+'histType'][2]]
1873	Limits = controlDict[pfx+'Limits']
1874	print >>pFile,' Instrument type: ',Sample['Type']
1875	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
1876	PrintSampleParms(Sample)
1877	PrintInstParms(Inst)
1878	PrintBackground(Background)
1879	elif 'HKLF' in histogram:
1880	Histogram = Histograms[histogram]
1881	hId = Histogram['hId']
1882	pfx = ':'+str(hId)+':'
1883	controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
1884	Inst = Histogram['Instrument Parameters'][0]
1885	controlDict[pfx+'histType'] = Inst['Type'][0]
1886	histDict[pfx+'Lam'] = Inst['Lam'][1]
1887	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
1888	return histVary,histDict,controlDict
1889
1890	def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
1891
1892	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1893	Back = Background[0]
1894	DebyePeaks = Background[1]
1895	lenBack = len(Back[3:])
1896	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
1897	for i in range(lenBack):
1898	Back[3+i] = parmDict[pfx+'Back:'+str(i)]
1899	if pfx+'Back:'+str(i) in sigDict:
1900	backSig[i] = sigDict[pfx+'Back:'+str(i)]
1901	if DebyePeaks['nDebye']:
1902	for i in range(DebyePeaks['nDebye']):
1903	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
1904	for j,name in enumerate(names):
1905	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
1906	if name in sigDict:
1907	backSig[lenBack+3*i+j] = sigDict[name]
1908	if DebyePeaks['nPeaks']:
1909	for i in range(DebyePeaks['nPeaks']):
1910	names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
1911	pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
1912	for j,name in enumerate(names):
1913	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
1914	if name in sigDict:
1915	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
1916	return backSig
1917
1918	def SetInstParms(pfx,Inst,parmDict,sigDict):
1919	instSig = {}
1920	insKeys = Inst.keys()
1921	insKeys.sort()
1922	for item in insKeys:
1923	insName = pfx+item
1924	Inst[item][1] = parmDict[insName]
1925	if insName in sigDict:
1926	instSig[item] = sigDict[insName]
1927	else:
1928	instSig[item] = 0
1929	return instSig
1930
1931	def SetSampleParms(pfx,Sample,parmDict,sigDict):
1932	if 'Bragg' in Sample['Type']: #Bragg-Brentano
1933	sampSig = [0 for i in range(3)]
1934	for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering?
1935	Sample[item][0] = parmDict[pfx+item]
1936	if pfx+item in sigDict:
1937	sampSig[i] = sigDict[pfx+item]
1938	elif 'Debye' in Sample['Type']: #Debye-Scherrer
1939	sampSig = [0 for i in range(4)]
1940	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1941	Sample[item][0] = parmDict[pfx+item]
1942	if pfx+item in sigDict:
1943	sampSig[i] = sigDict[pfx+item]
1944	return sampSig
1945
1946	def PrintBackgroundSig(Background,backSig):
1947	Back = Background[0]
1948	DebyePeaks = Background[1]
1949	lenBack = len(Back[3:])
1950	valstr = ' value : '
1951	sigstr = ' sig : '
1952	refine = False
1953	for i,back in enumerate(Back[3:]):
1954	valstr += '%10.4g'%(back)
1955	if Back[1]:
1956	refine = True
1957	sigstr += '%10.4g'%(backSig[i])
1958	else:
1959	sigstr += 10*' '
1960	if refine:
1961	print >>pFile,'\n Background function: ',Back[0]
1962	print >>pFile,valstr
1963	print >>pFile,sigstr
1964	if DebyePeaks['nDebye']:
1965	ifAny = False
1966	ptfmt = "%12.3f"
1967	names = ' names :'
1968	ptstr = ' values:'
1969	sigstr = ' esds :'
1970	for item in sigDict:
1971	if 'Debye' in item:
1972	ifAny = True
1973	names += '%12s'%(item)
1974	ptstr += ptfmt%(parmDict[item])
1975	sigstr += ptfmt%(sigDict[item])
1976	if ifAny:
1977	print >>pFile,'\n Debye diffuse scattering coefficients'
1978	print >>pFile,names
1979	print >>pFile,ptstr
1980	print >>pFile,sigstr
1981	if DebyePeaks['nPeaks']:
1982	ifAny = False
1983	ptfmt = "%14.3f"
1984	names = ' names :'
1985	ptstr = ' values:'
1986	sigstr = ' esds :'
1987	for item in sigDict:
1988	if 'BkPk' in item:
1989	ifAny = True
1990	names += '%14s'%(item)
1991	ptstr += ptfmt%(parmDict[item])
1992	sigstr += ptfmt%(sigDict[item])
1993	if ifAny:
1994	print >>pFile,'\n Single peak coefficients'
1995	print >>pFile,names
1996	print >>pFile,ptstr
1997	print >>pFile,sigstr
1998
1999	def PrintInstParmsSig(Inst,instSig):
2000	ptlbls = ' names :'
2001	ptstr = ' value :'
2002	sigstr = ' sig :'
2003	refine = False
2004	insKeys = instSig.keys()
2005	insKeys.sort()
2006	for name in insKeys:
2007	if name not in ['Type','Lam1','Lam2','Azimuth']:
2008	ptlbls += '%12s' % (name)
2009	ptstr += '%12.6f' % (Inst[name][1])
2010	if instSig[name]:
2011	refine = True
2012	sigstr += '%12.6f' % (instSig[name])
2013	else:
2014	sigstr += 12*' '
2015	if refine:
2016	print >>pFile,'\n Instrument Parameters:'
2017	print >>pFile,ptlbls
2018	print >>pFile,ptstr
2019	print >>pFile,sigstr
2020
2021	def PrintSampleParmsSig(Sample,sampleSig):
2022	ptlbls = ' names :'
2023	ptstr = ' values:'
2024	sigstr = ' sig :'
2025	refine = False
2026	if 'Bragg' in Sample['Type']:
2027	for i,item in enumerate(['Scale','Shift','Transparency']):
2028	ptlbls += '%14s'%(item)
2029	ptstr += '%14.4f'%(Sample[item][0])
2030	if sampleSig[i]:
2031	refine = True
2032	sigstr += '%14.4f'%(sampleSig[i])
2033	else:
2034	sigstr += 14*' '
2035
2036	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2037	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2038	ptlbls += '%14s'%(item)
2039	ptstr += '%14.4f'%(Sample[item][0])
2040	if sampleSig[i]:
2041	refine = True
2042	sigstr += '%14.4f'%(sampleSig[i])
2043	else:
2044	sigstr += 14*' '
2045
2046	if refine:
2047	print >>pFile,'\n Sample Parameters:'
2048	print >>pFile,ptlbls
2049	print >>pFile,ptstr
2050	print >>pFile,sigstr
2051
2052	histoList = Histograms.keys()
2053	histoList.sort()
2054	for histogram in histoList:
2055	if 'PWDR' in histogram:
2056	Histogram = Histograms[histogram]
2057	hId = Histogram['hId']
2058	pfx = ':'+str(hId)+':'
2059	Background = Histogram['Background']
2060	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2061
2062	Inst = Histogram['Instrument Parameters'][0]
2063	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2064
2065	Sample = Histogram['Sample Parameters']
2066	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2067
2068	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2069	print >>pFile,135*'-'
2070	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
2071	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2072	if Print:
2073	print >>pFile,' Instrument type: ',Sample['Type']
2074	PrintSampleParmsSig(Sample,sampSig)
2075	PrintInstParmsSig(Inst,instSig)
2076	PrintBackgroundSig(Background,backSig)
2077
2078	################################################################################
2079	##### Penalty & restraint functions
2080	################################################################################
2081
2082	def penaltyFxn(HistoPhases,parmDict,varyList):
2083	Histograms,Phases,restraintDict = HistoPhases
2084	pNames = []
2085	pVals = []
2086	pWt = []
2087	negWt = {}
2088	for phase in Phases:
2089	pId = Phases[phase]['pId']
2090	negWt[pId] = Phases[phase]['General']['Pawley neg wt']
2091	General = Phases[phase]['General']
2092	textureData = General['SH Texture']
2093	SGData = General['SGData']
2094	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
2095	cell = General['Cell'][1:7]
2096	Amat,Bmat = G2lat.cell2AB(cell)
2097	if phase not in restraintDict:
2098	continue
2099	phaseRest = restraintDict[phase]
2100	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
2101	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
2102	['ChemComp','Sites'],['Texture','HKLs']]
2103	for name,rest in names:
2104	itemRest = phaseRest[name]
2105	if itemRest[rest] and itemRest['Use']:
2106	wt = itemRest['wtFactor']
2107	if name in ['Bond','Angle','Plane','Chiral']:
2108	for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
2109	pNames.append(str(pId)+':'+name+':'+str(i))
2110	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2111	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2112	if name == 'Bond':
2113	calc = G2mth.getRestDist(XYZ,Amat)
2114	elif name == 'Angle':
2115	calc = G2mth.getRestAngle(XYZ,Amat)
2116	elif name == 'Plane':
2117	calc = G2mth.getRestPlane(XYZ,Amat)
2118	elif name == 'Chiral':
2119	calc = G2mth.getRestChiral(XYZ,Amat)
2120	pVals.append(obs-calc)
2121	pWt.append(wt/esd**2)
2122	elif name in ['Torsion','Rama']:
2123	coeffDict = itemRest['Coeff']
2124	for i,[indx,ops,cofName,esd] in enumerate(torsionList):
2125	pNames.append(str(pId)+':'+name+':'+str(i))
2126	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2127	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2128	if name == 'Torsion':
2129	tor = G2mth.getRestTorsion(XYZ,Amat)
2130	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
2131	else:
2132	phi,psi = G2mth.getRestRama(XYZ,Amat)
2133	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
2134	pVals.append(obs-calc)
2135	pWt.append(wt/esd**2)
2136	elif name == 'ChemComp':
2137	for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
2138	pNames.append(str(pId)+':'+name+':'+str(i))
2139	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
2140	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
2141	calc = np.sum(mulfracfactors)
2142	pVals.append(obs-calc)
2143	pWt.append(wt/esd**2)
2144	elif name == 'Texture':
2145	SHkeys = textureData['SH Coeff'][1].keys()
2146	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
2147	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2148	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2149	for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
2150	PH = np.array(hkl)
2151	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2152	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
2153	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
2154	Z1 = -ma.masked_greater(Z,0.0)
2155	IndZ1 = np.array(ma.nonzero(Z1))
2156	for ind in IndZ1.T:
2157	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
2158	pVals.append(Z1[ind[0]][ind[1]])
2159	pWt.append(wt/esd1**2)
2160	if ifesd2:
2161	Z2 = 1.-Z
2162	for ind in np.ndindex(grid,grid):
2163	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
2164	pVals.append(Z1[ind[0]][ind[1]])
2165	pWt.append(wt/esd2**2)
2166
2167	for item in varyList:
2168	if 'PWLref' in item and parmDict[item] < 0.:
2169	pId = int(item.split(':')[0])
2170	if negWt[pId]:
2171	pNames.append(item)
2172	pVals.append(-parmDict[item])
2173	pWt.append(negWt[pId])
2174	pVals = np.array(pVals)
2175	pWt = np.array(pWt) #should this be np.sqrt?
2176	return pNames,pVals,pWt
2177
2178	def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
2179	Histograms,Phases,restraintDict = HistoPhases
2180	pDerv = np.zeros((len(varyList),len(pVal)))
2181	for phase in Phases:
2182	if phase not in restraintDict:
2183	continue
2184	pId = Phases[phase]['pId']
2185	General = Phases[phase]['General']
2186	SGData = General['SGData']
2187	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
2188	cell = General['Cell'][1:7]
2189	Amat,Bmat = G2lat.cell2AB(cell)
2190	textureData = General['SH Texture']
2191
2192	SHkeys = textureData['SH Coeff'][1].keys()
2193	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
2194	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2195	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2196	sam = SamSym[textureData['Model']]
2197	phaseRest = restraintDict[phase]
2198	names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
2199	'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
2200	'ChemComp':'Sites','Texture':'HKLs'}
2201	lasthkl = np.array([0,0,0])
2202	for ip,pName in enumerate(pNames):
2203	pnames = pName.split(':')
2204	if pId == int(pnames[0]):
2205	name = pnames[1]
2206	if not name: #empty for Pawley restraints; pName has '::' in it
2207	pDerv[varyList.index(pName)][ip] += 1.
2208	continue
2209	id = int(pnames[2])
2210	itemRest = phaseRest[name]
2211	if name in ['Bond','Angle','Plane','Chiral']:
2212	indx,ops,obs,esd = itemRest[names[name]][id]
2213	dNames = []
2214	for ind in indx:
2215	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
2216	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2217	if name == 'Bond':
2218	deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
2219	elif name == 'Angle':
2220	deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
2221	elif name == 'Plane':
2222	deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
2223	elif name == 'Chiral':
2224	deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
2225	elif name in ['Torsion','Rama']:
2226	coffDict = itemRest['Coeff']
2227	indx,ops,cofName,esd = itemRest[names[name]][id]
2228	dNames = []
2229	for ind in indx:
2230	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
2231	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2232	if name == 'Torsion':
2233	deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
2234	else:
2235	deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
2236	elif name == 'ChemComp':
2237	indx,factors,obs,esd = itemRest[names[name]][id]
2238	dNames = []
2239	for ind in indx:
2240	dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
2241	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
2242	deriv = mul*factors
2243	elif 'Texture' in name:
2244	deriv = []
2245	dNames = []
2246	hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
2247	hkl = np.array(hkl)
2248	if np.any(lasthkl-hkl):
2249	PH = np.array(hkl)
2250	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2251	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
2252	lasthkl = copy.copy(hkl)
2253	if 'unit' in name:
2254	pass
2255	else:
2256	gam = float(pnames[3])
2257	psi = float(pnames[4])
2258	for SHname in ODFln:
2259	l,m,n = eval(SHname[1:])
2260	Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
2261	dNames += [str(pId)+'::'+SHname]
2262	deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
2263	for dName,drv in zip(dNames,deriv):
2264	try:
2265	ind = varyList.index(dName)
2266	pDerv[ind][ip] += drv
2267	except ValueError:
2268	pass
2269	return pDerv
2270
2271	################################################################################
2272	##### Function & derivative calculations
2273	################################################################################
2274
2275	def GetAtomFXU(pfx,calcControls,parmDict):
2276	Natoms = calcControls['Natoms'][pfx]
2277	Tdata = Natoms*[' ',]
2278	Mdata = np.zeros(Natoms)
2279	IAdata = Natoms*[' ',]
2280	Fdata = np.zeros(Natoms)
2281	FFdata = []
2282	BLdata = []
2283	Xdata = np.zeros((3,Natoms))
2284	dXdata = np.zeros((3,Natoms))
2285	Uisodata = np.zeros(Natoms)
2286	Uijdata = np.zeros((6,Natoms))
2287	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
2288	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
2289	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
2290	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
2291	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
2292	for iatm in range(Natoms):
2293	for key in keys:
2294	parm = pfx+key+str(iatm)
2295	if parm in parmDict:
2296	keys[key][iatm] = parmDict[parm]
2297	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
2298
2299	def getFFvalues(FFtables,SQ):
2300	FFvals = {}
2301	for El in FFtables:
2302	FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
2303	return FFvals
2304
2305	def getBLvalues(BLtables):
2306	BLvals = {}
2307	for El in BLtables:
2308	BLvals[El] = BLtables[El][1][1]
2309	return BLvals
2310
2311	def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2312	''' Compute structure factors for all h,k,l for phase
2313	input:
2314	refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2315	G: reciprocal metric tensor
2316	pfx: phase id string
2317	SGData: space group info. dictionary output from SpcGroup
2318	calcControls:
2319	ParmDict:
2320	puts result F^2 in each ref[8] in refList
2321	'''
2322	twopi = 2.0*np.pi
2323	twopisq = 2.0np.pi*2
2324	phfx = pfx.split(':')[0]+hfx
2325	ast = np.sqrt(np.diag(G))
2326	Mast = twopisq*np.multiply.outer(ast,ast)
2327	FFtables = calcControls['FFtables']
2328	BLtables = calcControls['BLtables']
2329	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2330	FF = np.zeros(len(Tdata))
2331	if 'N' in parmDict[hfx+'Type']:
2332	FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
2333	else:
2334	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2335	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2336	maxPos = len(SGData['SGOps'])
2337	Uij = np.array(G2lat.U6toUij(Uijdata))
2338	bij = Mast*Uij.T
2339	for refl in refList:
2340	fbs = np.array([0,0])
2341	H = refl[:3]
2342	SQ = 1./(2.refl[4])*2
2343	SQfactor = 4.0SQtwopisq
2344	Bab = parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor)
2345	if not len(refl[-1]): #no form factors
2346	if 'N' in parmDict[hfx+'Type']:
2347	refl[-1] = getBLvalues(BLtables)
2348	else: #'X'
2349	refl[-1] = getFFvalues(FFtables,SQ)
2350	for i,El in enumerate(Tdata):
2351	FF[i] = refl[-1][El]
2352	Uniq = refl[11]
2353	phi = refl[12]
2354	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
2355	sinp = np.sin(phase)
2356	cosp = np.cos(phase)
2357	occ = Mdata*Fdata/len(Uniq)
2358	biso = -SQfactor*Uisodata
2359	Tiso = np.where(biso<1.,np.exp(biso),1.0)
2360	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
2361	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
2362	Tcorr = Tiso*Tuij
2363	fa = np.array([(FF+FP-Bab)occcospTcorr,-FPPoccsinpTcorr])
2364	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
2365	if not SGData['SGInv']:
2366	fb = np.array([(FF+FP-Bab)occsinpTcorr,FPPocccospTcorr])
2367	fbs = np.sum(np.sum(fb,axis=1),axis=1)
2368	fasq = fas**2
2369	fbsq = fbs**2 #imaginary
2370	refl[9] = np.sum(fasq)+np.sum(fbsq)
2371	refl[10] = atan2d(fbs[0],fas[0])
2372	return refList
2373
2374	def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2375	twopi = 2.0*np.pi
2376	twopisq = 2.0np.pi*2
2377	phfx = pfx.split(':')[0]+hfx
2378	ast = np.sqrt(np.diag(G))
2379	Mast = twopisq*np.multiply.outer(ast,ast)
2380	FFtables = calcControls['FFtables']
2381	BLtables = calcControls['BLtables']
2382	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2383	FF = np.zeros(len(Tdata))
2384	if 'N' in parmDict[hfx+'Type']:
2385	FP = 0.
2386	FPP = 0.
2387	else:
2388	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2389	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2390	maxPos = len(SGData['SGOps'])
2391	Uij = np.array(G2lat.U6toUij(Uijdata))
2392	bij = Mast*Uij.T
2393	dFdvDict = {}
2394	dFdfr = np.zeros((len(refList),len(Mdata)))
2395	dFdx = np.zeros((len(refList),len(Mdata),3))
2396	dFdui = np.zeros((len(refList),len(Mdata)))
2397	dFdua = np.zeros((len(refList),len(Mdata),6))
2398	dFdbab = np.zeros((len(refList),2))
2399	for iref,refl in enumerate(refList):
2400	H = np.array(refl[:3])
2401	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
2402	SQfactor = 8.0SQnp.pi**2
2403	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2404	Bab = parmDict[phfx+'BabA']*dBabdA
2405	for i,El in enumerate(Tdata):
2406	FF[i] = refl[-1][El]
2407	Uniq = refl[11]
2408	phi = refl[12]
2409	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
2410	sinp = np.sin(phase)
2411	cosp = np.cos(phase)
2412	occ = Mdata*Fdata/len(Uniq)
2413	biso = -SQfactor*Uisodata
2414	Tiso = np.where(biso<1.,np.exp(biso),1.0)
2415	HbH = -np.inner(H,np.inner(bij,H))
2416	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
2417	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
2418	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
2419	Tcorr = Tiso*Tuij
2420	fot = (FF+FP-Bab)occTcorr
2421	fotp = FPPoccTcorr
2422	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
2423	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
2424
2425	fas = np.sum(np.sum(fa,axis=1),axis=1)
2426	fbs = np.sum(np.sum(fb,axis=1),axis=1)
2427	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
2428	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
2429	#sum below is over Uniq
2430	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
2431	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
2432	dfadui = np.sum(-SQfactor*fa,axis=2)
2433	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
2434	dfadba = np.sum(-cosp(occTcorr)[:,np.newaxis],axis=1)
2435	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for
2436	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
2437	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
2438	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
2439	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
2440	dFdbab[iref] = np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
2441	if not SGData['SGInv']:
2442	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom
2443	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
2444	dfbdui = np.sum(-SQfactor*fb,axis=2)
2445	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
2446	dfbdba = np.sum(-sinp(occTcorr)[:,np.newaxis],axis=1)
2447	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
2448	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
2449	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
2450	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
2451	dFdbab[iref] += np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
2452	#loop over atoms - each dict entry is list of derivatives for all the reflections
2453	for i in range(len(Mdata)):
2454	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2455	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2456	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2457	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2458	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2459	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2460	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2461	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2462	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2463	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2464	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2465	dFdvDict[pfx+'BabA'] = dFdbab.T[0]
2466	dFdvDict[pfx+'BabU'] = dFdbab.T[1]
2467	return dFdvDict
2468
2469	def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
2470	''' Single crystal extinction function; puts correction in ref[13] and returns
2471	corrections needed for derivatives
2472	'''
2473	ref[13] = 1.0
2474	dervCor = 1.0
2475	dervDict = {}
2476	if calcControls[phfx+'EType'] != 'None':
2477	cos2T = 1.0-0.5(parmDict[hfx+'Lam']/ref[4])*2 #cos(2theta)
2478	if 'SXC' in parmDict[hfx+'Type']:
2479	AV = 7.9406e5/parmDict[pfx+'Vol']**2
2480	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2481	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
2482	elif 'SNT' in parmDict[hfx+'Type']:
2483	AV = 1.e7/parmDict[pfx+'Vol']**2
2484	PL = 1./(4.refl[4]*2)
2485	P12 = 1.0
2486	elif 'SNC' in parmDict[hfx+'Type']:
2487	AV = 1.e7/parmDict[pfx+'Vol']**2
2488	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2489	P12 = 1.0
2490
2491	PLZ = AVP12parmDict[hfx+'Lam']*2ref[7]
2492	if 'Primary' in calcControls[phfx+'EType']:
2493	PLZ *= 1.5
2494	else:
2495	PLZ *= calcControls[phfx+'Tbar']
2496
2497	if 'Primary' in calcControls[phfx+'EType']:
2498	PSIG = parmDict[phfx+'Ep']
2499	elif 'I & II' in calcControls[phfx+'EType']:
2500	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
2501	elif 'Type II' in calcControls[phfx+'EType']:
2502	PSIG = parmDict[phfx+'Es']
2503	else: # 'Secondary Type I'
2504	PSIG = parmDict[phfx+'Eg']/PL
2505
2506	AG = 0.58+0.48cos2T+0.24cos2T**2
2507	AL = 0.025+0.285*cos2T
2508	BG = 0.02-0.025*cos2T
2509	BL = 0.15-0.2(0.75-cos2T)*2
2510	if cos2T < 0.:
2511	BL = -0.45*cos2T
2512	CG = 2.
2513	CL = 2.
2514	PF = PLZ*PSIG
2515
2516	if 'Gaussian' in calcControls[phfx+'EApprox']:
2517	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
2518	extCor = np.sqrt(PF4)
2519	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
2520	else:
2521	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
2522	extCor = np.sqrt(PF4)
2523	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
2524
2525	dervCor = (1.+PF)*PF3
2526	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2527	dervDict[phfx+'Ep'] = -ref[7]PLZPF3
2528	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2529	dervDict[phfx+'Es'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
2530	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2531	dervDict[phfx+'Eg'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
2532
2533	ref[13] = 1./extCor
2534	return dervCor,dervDict
2535
2536
2537	def Dict2Values(parmdict, varylist):
2538	'''Use before call to leastsq to setup list of values for the parameters
2539	in parmdict, as selected by key in varylist'''
2540	return [parmdict[key] for key in varylist]
2541
2542	def Values2Dict(parmdict, varylist, values):
2543	''' Use after call to leastsq to update the parameter dictionary with
2544	values corresponding to keys in varylist'''
2545	parmdict.update(zip(varylist,values))
2546
2547	def GetNewCellParms(parmDict,varyList):
2548	newCellDict = {}
2549	Anames = ['A'+str(i) for i in range(6)]
2550	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2551	for item in varyList:
2552	keys = item.split(':')
2553	if keys[2] in Ddict:
2554	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
2555	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
2556	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2557	return newCellDict # is e.g. {'0::D11':A0+D11}
2558
2559	def ApplyXYZshifts(parmDict,varyList):
2560	''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2561	input:
2562	parmDict - parameter dictionary
2563	varyList - list of variables
2564	returns:
2565	newAtomDict - dictionary of new atomic coordinate names & values;
2566	key is parameter shift name
2567	'''
2568	newAtomDict = {}
2569	for item in parmDict:
2570	if 'dA' in item:
2571	parm = ''.join(item.split('d'))
2572	parmDict[parm] += parmDict[item]
2573	newAtomDict[item] = [parm,parmDict[parm]]
2574	return newAtomDict
2575
2576	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2577	#Spherical harmonics texture
2578	IFCoup = 'Bragg' in calcControls[hfx+'instType']
2579	odfCor = 1.0
2580	H = refl[:3]
2581	cell = G2lat.Gmat2cell(g)
2582	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2583	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2584	phi,beta = G2lat.CrsAng(H,cell,SGData)
2585	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2586	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2587	for item in SHnames:
2588	L,M,N = eval(item.strip('C'))
2589	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2590	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2591	Lnorm = G2lat.Lnorm(L)
2592	odfCor += parmDict[pfx+item]LnormKcl*Ksl
2593	return odfCor
2594
2595	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2596	#Spherical harmonics texture derivatives
2597	FORPI = 4.0*np.pi
2598	IFCoup = 'Bragg' in calcControls[hfx+'instType']
2599	odfCor = 1.0
2600	dFdODF = {}
2601	dFdSA = [0,0,0]
2602	H = refl[:3]
2603	cell = G2lat.Gmat2cell(g)
2604	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2605	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2606	phi,beta = G2lat.CrsAng(H,cell,SGData)
2607	psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
2608	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2609	for item in SHnames:
2610	L,M,N = eval(item.strip('C'))
2611	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2612	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2613	Lnorm = G2lat.Lnorm(L)
2614	odfCor += parmDict[pfx+item]LnormKcl*Ksl
2615	dFdODF[pfx+item] = LnormKclKsl
2616	for i in range(3):
2617	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
2618	return odfCor,dFdODF,dFdSA
2619
2620	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2621	#sphericaql harmonics preferred orientation (cylindrical symmetry only)
2622	odfCor = 1.0
2623	H = refl[:3]
2624	cell = G2lat.Gmat2cell(g)
2625	Sangl = [0.,0.,0.]
2626	if 'Bragg' in calcControls[hfx+'instType']:
2627	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2628	IFCoup = True
2629	else:
2630	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2631	IFCoup = False
2632	phi,beta = G2lat.CrsAng(H,cell,SGData)
2633	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2634	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2635	for item in SHnames:
2636	L,N = eval(item.strip('C'))
2637	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
2638	odfCor += parmDict[phfx+item]LnormKcsl
2639	return odfCor
2640
2641	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2642	#spherical harmonics preferred orientation derivatives (cylindrical symmetry only)
2643	FORPI = 12.5663706143592
2644	odfCor = 1.0
2645	dFdODF = {}
2646	H = refl[:3]
2647	cell = G2lat.Gmat2cell(g)
2648	Sangl = [0.,0.,0.]
2649	if 'Bragg' in calcControls[hfx+'instType']:
2650	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2651	IFCoup = True
2652	else:
2653	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2654	IFCoup = False
2655	phi,beta = G2lat.CrsAng(H,cell,SGData)
2656	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2657	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2658	for item in SHnames:
2659	L,N = eval(item.strip('C'))
2660	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
2661	odfCor += parmDict[phfx+item]LnormKcsl
2662	dFdODF[phfx+item] = Kcsl*Lnorm
2663	return odfCor,dFdODF
2664
2665	def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2666	POcorr = 1.0
2667	if calcControls[phfx+'poType'] == 'MD':
2668	MD = parmDict[phfx+'MD']
2669	if MD != 1.0:
2670	MDAxis = calcControls[phfx+'MDAxis']
2671	sumMD = 0
2672	for H in refl[11]:
2673	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2674	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
2675	sumMD += A**3
2676	POcorr = sumMD/len(refl[11])
2677	else: #spherical harmonics
2678	if calcControls[phfx+'SHord']:
2679	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2680	return POcorr
2681
2682	def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2683	POcorr = 1.0
2684	POderv = {}
2685	if calcControls[phfx+'poType'] == 'MD':
2686	MD = parmDict[phfx+'MD']
2687	MDAxis = calcControls[phfx+'MDAxis']
2688	sumMD = 0
2689	sumdMD = 0
2690	for H in refl[11]:
2691	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2692	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
2693	sumMD += A**3
2694	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
2695	POcorr = sumMD/len(refl[11])
2696	POderv[phfx+'MD'] = sumdMD/len(refl[11])
2697	else: #spherical harmonics
2698	if calcControls[phfx+'SHord']:
2699	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2700	return POcorr,POderv
2701
2702	def GetAbsorb(refl,hfx,calcControls,parmDict):
2703	if 'Debye' in calcControls[hfx+'instType']:
2704	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
2705	else:
2706	return 1.0
2707
2708	def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
2709	if 'Debye' in calcControls[hfx+'instType']:
2710	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
2711	else:
2712	return 0.0
2713
2714	def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2715	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
2716	if 'X' in parmDict[hfx+'Type']:
2717	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
2718	Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2719	if pfx+'SHorder' in parmDict:
2720	Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2721	Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
2722	refl[13] = Icorr
2723
2724	def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2725	dIdsh = 1./parmDict[hfx+'Scale']
2726	dIdsp = 1./parmDict[phfx+'Scale']
2727	if 'X' in parmDict[hfx+'Type']:
2728	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
2729	dIdPola /= pola
2730	else: #'N'
2731	dIdPola = 0.0
2732	POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2733	for iPO in dIdPO:
2734	dIdPO[iPO] /= POcorr
2735	dFdODF = {}
2736	dFdSA = [0,0,0]
2737	if pfx+'SHorder' in parmDict:
2738	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2739	for iSH in dFdODF:
2740	dFdODF[iSH] /= odfCor
2741	for i in range(3):
2742	dFdSA[i] /= odfCor
2743	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
2744	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
2745
2746	def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
2747	costh = cosd(refl[5]/2.)
2748	#crystallite size
2749	if calcControls[phfx+'SizeType'] == 'isotropic':
2750	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
2751	elif calcControls[phfx+'SizeType'] == 'uniaxial':
2752	H = np.array(refl[:3])
2753	P = np.array(calcControls[phfx+'SizeAxis'])
2754	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2755	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
2756	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
2757	else: #ellipsoidal crystallites
2758	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2759	H = np.array(refl[:3])
2760	lenR = G2pwd.ellipseSize(H,Sij,GB)
2761	Sgam = 1.8wave/(np.picosth*lenR)
2762	#microstrain
2763	if calcControls[phfx+'MustrainType'] == 'isotropic':
2764	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
2765	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2766	H = np.array(refl[:3])
2767	P = np.array(calcControls[phfx+'MustrainAxis'])
2768	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2769	Si = parmDict[phfx+'Mustrain;i']
2770	Sa = parmDict[phfx+'Mustrain;a']
2771	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
2772	else: #generalized - P.W. Stephens model
2773	pwrs = calcControls[phfx+'MuPwrs']
2774	sum = 0
2775	for i,pwr in enumerate(pwrs):
2776	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
2777	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
2778	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
2779	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
2780	sig /= ateln2
2781	return sig,gam
2782
2783	def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
2784	gamDict = {}
2785	sigDict = {}
2786	costh = cosd(refl[5]/2.)
2787	tanth = tand(refl[5]/2.)
2788	#crystallite size derivatives
2789	if calcControls[phfx+'SizeType'] == 'isotropic':
2790	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
2791	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
2792	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
2793	elif calcControls[phfx+'SizeType'] == 'uniaxial':
2794	H = np.array(refl[:3])
2795	P = np.array(calcControls[phfx+'SizeAxis'])
2796	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2797	Si = parmDict[phfx+'Size;i']
2798	Sa = parmDict[phfx+'Size;a']
2799	gami = (1.8wave/np.pi)/(SiSa)
2800	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
2801	Sgam = gami*sqtrm
2802	gam = Sgam/costh
2803	dsi = (gamiSicosP*2/(sqtrmcosth)-gam/Si)
2804	dsa = (gamiSasinP*2/(sqtrmcosth)-gam/Sa)
2805	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2806	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2807	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
2808	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
2809	else: #ellipsoidal crystallites
2810	const = 1.8wave/(np.picosth)
2811	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2812	H = np.array(refl[:3])
2813	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2814	Sgam = 1.8wave/(np.picosth*lenR)
2815	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2816	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
2817	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
2818	gamDict[phfx+'Size;mx'] = Sgam
2819	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
2820
2821	#microstrain derivatives
2822	if calcControls[phfx+'MustrainType'] == 'isotropic':
2823	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
2824	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
2825	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
2826	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2827	H = np.array(refl[:3])
2828	P = np.array(calcControls[phfx+'MustrainAxis'])
2829	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2830	Si = parmDict[phfx+'Mustrain;i']
2831	Sa = parmDict[phfx+'Mustrain;a']
2832	gami = 0.018SiSa*tanth/np.pi
2833	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
2834	Mgam = gami/sqtrm
2835	dsi = -gamiSicosP2/sqtrm3
2836	dsa = -gamiSasinP2/sqtrm3
2837	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2838	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2839	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
2840	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
2841	else: #generalized - P.W. Stephens model
2842	pwrs = calcControls[phfx+'MuPwrs']
2843	const = 0.018refl[4]2tanth
2844	sum = 0
2845	for i,pwr in enumerate(pwrs):
2846	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
2847	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
2848	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain;mx']
2849	sigDict[phfx+'Mustrain:'+str(i)] = \
2850	2.constterm(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
2851	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
2852	for i in range(len(pwrs)):
2853	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
2854	gamDict[phfx+'Mustrain;mx'] = Mgam
2855	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2856	return sigDict,gamDict
2857
2858	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
2859	h,k,l = refl[:3]
2860	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
2861	d = np.sqrt(dsq)
2862
2863	refl[4] = d
2864	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
2865	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
2866	if 'Bragg' in calcControls[hfx+'instType']:
2867	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2868	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
2869	else: #Debye-Scherrer - simple but maybe not right
2870	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2871	return pos
2872
2873	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
2874	dpr = 180./np.pi
2875	h,k,l = refl[:3]
2876	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2877	dst = np.sqrt(dstsq)
2878	pos = refl[5]-parmDict[hfx+'Zero']
2879	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
2880	dpdw = const*dst
2881	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
2882	dpdA = constwave/(2.0*dst)
2883	dpdZ = 1.0
2884	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
2885	if 'Bragg' in calcControls[hfx+'instType']:
2886	dpdSh = -4.constcosd(pos/2.0)
2887	dpdTr = -constsind(pos)100.0
2888	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2889	else: #Debye-Scherrer - simple but maybe not right
2890	dpdXd = -const*cosd(pos)
2891	dpdYd = -const*sind(pos)
2892	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2893
2894	def GetHStrainShift(refl,SGData,phfx,parmDict):
2895	laue = SGData['SGLaue']
2896	uniq = SGData['SGUniq']
2897	h,k,l = refl[:3]
2898	if laue in ['m3','m3m']:
2899	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
2900	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
2901	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2902	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
2903	elif laue in ['3R','3mR']:
2904	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
2905	elif laue in ['4/m','4/mmm']:
2906	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
2907	elif laue in ['mmm']:
2908	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2909	elif laue in ['2/m']:
2910	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2911	if uniq == 'a':
2912	Dij += parmDict[phfx+'D23']kl
2913	elif uniq == 'b':
2914	Dij += parmDict[phfx+'D13']hl
2915	elif uniq == 'c':
2916	Dij += parmDict[phfx+'D12']hk
2917	else:
2918	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
2919	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
2920	return -Dijrefl[4]2tand(refl[5]/2.0)
2921
2922	def GetHStrainShiftDerv(refl,SGData,phfx):
2923	laue = SGData['SGLaue']
2924	uniq = SGData['SGUniq']
2925	h,k,l = refl[:3]
2926	if laue in ['m3','m3m']:
2927	dDijDict = {phfx+'D11':h2+k2+l**2,
2928	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
2929	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2930	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
2931	elif laue in ['3R','3mR']:
2932	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
2933	elif laue in ['4/m','4/mmm']:
2934	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
2935	elif laue in ['mmm']:
2936	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2937	elif laue in ['2/m']:
2938	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2939	if uniq == 'a':
2940	dDijDict[phfx+'D23'] = k*l
2941	elif uniq == 'b':
2942	dDijDict[phfx+'D13'] = h*l
2943	elif uniq == 'c':
2944	dDijDict[phfx+'D12'] = h*k
2945	names.append()
2946	else:
2947	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
2948	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
2949	for item in dDijDict:
2950	dDijDict[item] = -refl[4]2tand(refl[5]/2.0)
2951	return dDijDict
2952
2953	def GetFprime(controlDict,Histograms):
2954	FFtables = controlDict['FFtables']
2955	if not FFtables:
2956	return
2957	histoList = Histograms.keys()
2958	histoList.sort()
2959	for histogram in histoList:
2960	if histogram[:4] in ['PWDR','HKLF']:
2961	Histogram = Histograms[histogram]
2962	hId = Histogram['hId']
2963	hfx = ':%d:'%(hId)
2964	keV = controlDict[hfx+'keV']
2965	for El in FFtables:
2966	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2967	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2968	FFtables[El][hfx+'FP'] = FP
2969	FFtables[El][hfx+'FPP'] = FPP
2970
2971	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2972	histoList = Histograms.keys()
2973	histoList.sort()
2974	for histogram in histoList:
2975	if 'PWDR' in histogram[:4]:
2976	Histogram = Histograms[histogram]
2977	hId = Histogram['hId']
2978	hfx = ':%d:'%(hId)
2979	Limits = calcControls[hfx+'Limits']
2980	shl = max(parmDict[hfx+'SH/L'],0.0005)
2981	Ka2 = False
2982	kRatio = 0.0
2983	if hfx+'Lam1' in parmDict.keys():
2984	Ka2 = True
2985	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2986	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2987	x,y,w,yc,yb,yd = Histogram['Data']
2988	xB = np.searchsorted(x,Limits[0])
2989	xF = np.searchsorted(x,Limits[1])
2990	ymb = np.array(y-yb)
2991	ymb = np.where(ymb,ymb,1.0)
2992	ycmb = np.array(yc-yb)
2993	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
2994	refLists = Histogram['Reflection Lists']
2995	for phase in refLists:
2996	Phase = Phases[phase]
2997	pId = Phase['pId']
2998	phfx = '%d:%d:'%(pId,hId)
2999	refList = refLists[phase]
3000	sumFo = 0.0
3001	sumdF = 0.0
3002	sumFosq = 0.0
3003	sumdFsq = 0.0
3004	for refl in refList:
3005	if 'C' in calcControls[hfx+'histType']:
3006	yp = np.zeros_like(yb)
3007	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3008	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
3009	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
3010	iFin2 = iFin
3011	yp[iBeg:iFin] = refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
3012	if Ka2:
3013	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
3014	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
3015	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
3016	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
3017	if not iBeg2+iFin2: #peak below low limit - skip peak
3018	continue
3019	elif not iBeg2-iFin2: #peak above high limit - done
3020	break
3021	yp[iBeg2:iFin2] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here
3022	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[13](1.+kRatio)),0.0))
3023	elif 'T' in calcControls[hfx+'histType']:
3024	print 'TOF Undefined at present'
3025	raise Exception #no TOF yet
3026	Fo = np.sqrt(np.abs(refl[8]))
3027	Fc = np.sqrt(np.abs(refl[9]))
3028	sumFo += Fo
3029	sumFosq += refl[8]**2
3030	sumdF += np.abs(Fo-Fc)
3031	sumdFsq += (refl[8]-refl[9])**2
3032	Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3033	Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3034	Histogram[phfx+'Nref'] = len(refList)
3035
3036	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3037
3038	def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
3039	U = parmDict[hfx+'U']
3040	V = parmDict[hfx+'V']
3041	W = parmDict[hfx+'W']
3042	X = parmDict[hfx+'X']
3043	Y = parmDict[hfx+'Y']
3044	tanPos = tand(refl[5]/2.0)
3045	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
3046	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
3047	sig = max(0.001,sig)
3048	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
3049	gam = max(0.001,gam)
3050	return sig,gam
3051
3052	hId = Histogram['hId']
3053	hfx = ':%d:'%(hId)
3054	bakType = calcControls[hfx+'bakType']
3055	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
3056	yc = np.zeros_like(yb)
3057
3058	if 'C' in calcControls[hfx+'histType']:
3059	shl = max(parmDict[hfx+'SH/L'],0.002)
3060	Ka2 = False
3061	if hfx+'Lam1' in parmDict.keys():
3062	wave = parmDict[hfx+'Lam1']
3063	Ka2 = True
3064	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
3065	kRatio = parmDict[hfx+'I(L2)/I(L1)']
3066	else:
3067	wave = parmDict[hfx+'Lam']
3068	else:
3069	print 'TOF Undefined at present'
3070	raise ValueError
3071	for phase in Histogram['Reflection Lists']:
3072	refList = Histogram['Reflection Lists'][phase]
3073	Phase = Phases[phase]
3074	pId = Phase['pId']
3075	pfx = '%d::'%(pId)
3076	phfx = '%d:%d:'%(pId,hId)
3077	hfx = ':%d:'%(hId)
3078	SGData = Phase['General']['SGData']
3079	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
3080	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3081	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
3082	Vst = np.sqrt(nl.det(G)) #V*
3083	if not Phase['General'].get('doPawley'):
3084	refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
3085	for refl in refList:
3086	if 'C' in calcControls[hfx+'histType']:
3087	h,k,l = refl[:3]
3088	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
3089	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
3090	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
3091	refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam
3092	GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13]
3093	refl[13] = VstLorenz
3094	if Phase['General'].get('doPawley'):
3095	try:
3096	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3097	refl[9] = parmDict[pInd]
3098	except KeyError:
3099	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
3100	continue
3101	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3102	iBeg = np.searchsorted(x,refl[5]-fmin)
3103	iFin = np.searchsorted(x,refl[5]+fmax)
3104	if not iBeg+iFin: #peak below low limit - skip peak
3105	continue
3106	elif not iBeg-iFin: #peak above high limit - done
3107	break
3108	yc[iBeg:iFin] += refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
3109	if Ka2:
3110	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
3111	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
3112	iBeg = np.searchsorted(x,pos2-fmin)
3113	iFin = np.searchsorted(x,pos2+fmax)
3114	if not iBeg+iFin: #peak below low limit - skip peak
3115	continue
3116	elif not iBeg-iFin: #peak above high limit - done
3117	return yc,yb
3118	yc[iBeg:iFin] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here
3119	elif 'T' in calcControls[hfx+'histType']:
3120	print 'TOF Undefined at present'
3121	raise Exception #no TOF yet
3122	return yc,yb
3123
3124	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3125
3126	def cellVaryDerv(pfx,SGData,dpdA):
3127	if SGData['SGLaue'] in ['-1',]:
3128	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
3129	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
3130	elif SGData['SGLaue'] in ['2/m',]:
3131	if SGData['SGUniq'] == 'a':
3132	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
3133	elif SGData['SGUniq'] == 'b':
3134	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
3135	else:
3136	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
3137	elif SGData['SGLaue'] in ['mmm',]:
3138	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
3139	elif SGData['SGLaue'] in ['4/m','4/mmm']:
3140	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
3141	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
3142	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
3143	elif SGData['SGLaue'] in ['3R', '3mR']:
3144	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
3145	elif SGData['SGLaue'] in ['m3m','m3']:
3146	return [[pfx+'A0',dpdA[0]]]
3147	# create a list of dependent variables and set up a dictionary to hold their derivatives
3148	dependentVars = G2mv.GetDependentVars()
3149	depDerivDict = {}
3150	for j in dependentVars:
3151	depDerivDict[j] = np.zeros(shape=(len(x)))
3152	#print 'dependent vars',dependentVars
3153	lenX = len(x)
3154	hId = Histogram['hId']
3155	hfx = ':%d:'%(hId)
3156	bakType = calcControls[hfx+'bakType']
3157	dMdv = np.zeros(shape=(len(varylist),len(x)))
3158	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
3159	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
3160	bBpos =varylist.index(hfx+'Back:0')
3161	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
3162	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
3163	for name in varylist:
3164	if 'Debye' in name:
3165	id = int(name.split(':')[-1])
3166	parm = name[:int(name.rindex(':'))]
3167	ip = names.index(parm)
3168	dMdv[varylist.index(name)] = dMddb[3*id+ip]
3169	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
3170	for name in varylist:
3171	if 'BkPk' in name:
3172	id = int(name.split(':')[-1])
3173	parm = name[:int(name.rindex(':'))]
3174	ip = names.index(parm)
3175	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
3176	cw = np.diff(x)
3177	cw = np.append(cw,cw[-1])
3178	if 'C' in calcControls[hfx+'histType']:
3179	shl = max(parmDict[hfx+'SH/L'],0.002)
3180	Ka2 = False
3181	if hfx+'Lam1' in parmDict.keys():
3182	wave = parmDict[hfx+'Lam1']
3183	Ka2 = True
3184	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
3185	kRatio = parmDict[hfx+'I(L2)/I(L1)']
3186	else:
3187	wave = parmDict[hfx+'Lam']
3188	else:
3189	print 'TOF Undefined at present'
3190	raise ValueError
3191	for phase in Histogram['Reflection Lists']:
3192	refList = Histogram['Reflection Lists'][phase]
3193	Phase = Phases[phase]
3194	SGData = Phase['General']['SGData']
3195	pId = Phase['pId']
3196	pfx = '%d::'%(pId)
3197	phfx = '%d:%d:'%(pId,hId)
3198	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
3199	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3200	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
3201	if not Phase['General'].get('doPawley'):
3202	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
3203	for iref,refl in enumerate(refList):
3204	if 'C' in calcControls[hfx+'histType']: #CW powder
3205	h,k,l = refl[:3]
3206	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3207	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3208	iBeg = np.searchsorted(x,refl[5]-fmin)
3209	iFin = np.searchsorted(x,refl[5]+fmax)
3210	if not iBeg+iFin: #peak below low limit - skip peak
3211	continue
3212	elif not iBeg-iFin: #peak above high limit - done
3213	break
3214	pos = refl[5]
3215	tanth = tand(pos/2.0)
3216	costh = cosd(pos/2.0)
3217	lenBF = iFin-iBeg
3218	dMdpk = np.zeros(shape=(6,lenBF))
3219	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
3220	for i in range(5):
3221	dMdpk[i] += 100.cw[iBeg:iFin]refl[13]refl[9]dMdipk[i]
3222	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
3223	if Ka2:
3224	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
3225	iBeg2 = np.searchsorted(x,pos2-fmin)
3226	iFin2 = np.searchsorted(x,pos2+fmax)
3227	if iBeg2-iFin2:
3228	lenBF2 = iFin2-iBeg2
3229	dMdpk2 = np.zeros(shape=(6,lenBF2))
3230	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
3231	for i in range(5):
3232	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[13]refl[9]kRatio*dMdipk2[i]
3233	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[13]*dMdipk2[0]
3234	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
3235	if Phase['General'].get('doPawley'):
3236	dMdpw = np.zeros(len(x))
3237	try:
3238	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3239	idx = varylist.index(pIdx)
3240	dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
3241	if Ka2:
3242	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
3243	dMdv[idx] = dMdpw
3244	except: # ValueError:
3245	pass
3246	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
3247	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
3248	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
3249	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
3250	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
3251	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
3252	hfx+'DisplaceY':[dpdY,'pos'],hfx+'Absorption':[dFdAb,'int'],}
3253	for name in names:
3254	item = names[name]
3255	if name in varylist:
3256	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
3257	if Ka2:
3258	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
3259	elif name in dependentVars:
3260	if Ka2:
3261	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
3262	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
3263	for iPO in dIdPO:
3264	if iPO in varylist:
3265	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
3266	if Ka2:
3267	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
3268	elif iPO in dependentVars:
3269	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
3270	if Ka2:
3271	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
3272	for i,name in enumerate(['omega','chi','phi']):
3273	aname = pfx+'SH '+name
3274	if aname in varylist:
3275	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
3276	if Ka2:
3277	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
3278	elif aname in dependentVars:
3279	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
3280	if Ka2:
3281	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
3282	for iSH in dFdODF:
3283	if iSH in varylist:
3284	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
3285	if Ka2:
3286	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
3287	elif iSH in dependentVars:
3288	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
3289	if Ka2:
3290	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
3291	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
3292	for name,dpdA in cellDervNames:
3293	if name in varylist:
3294	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
3295	if Ka2:
3296	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
3297	elif name in dependentVars:
3298	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
3299	if Ka2:
3300	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
3301	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
3302	for name in dDijDict:
3303	if name in varylist:
3304	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
3305	if Ka2:
3306	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
3307	elif name in dependentVars:
3308	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
3309	if Ka2:
3310	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
3311	sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
3312	for name in gamDict:
3313	if name in varylist:
3314	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
3315	if Ka2:
3316	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
3317	elif name in dependentVars:
3318	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
3319	if Ka2:
3320	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
3321	for name in sigDict:
3322	if name in varylist:
3323	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
3324	if Ka2:
3325	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
3326	elif name in dependentVars:
3327	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
3328	if Ka2:
3329	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
3330	for name in ['BabA','BabU']:
3331	if phfx+name in varylist:
3332	dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]dervDict['int']cw[iBeg:iFin]
3333	if Ka2:
3334	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]dervDict2['int']cw[iBeg2:iFin2]
3335	elif phfx+name in dependentVars:
3336	depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]dervDict['int']cw[iBeg:iFin]
3337	if Ka2:
3338	depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]dervDict2['int']cw[iBeg2:iFin2]
3339	elif 'T' in calcControls[hfx+'histType']:
3340	print 'TOF Undefined at present'
3341	raise Exception #no TOF yet
3342	#do atom derivatives - for F,X & U so far
3343	corr = dervDict['int']/refl[9]
3344	if Ka2:
3345	corr2 = dervDict2['int']/refl[9]
3346	for name in varylist+dependentVars:
3347	try:
3348	aname = name.split(pfx)[1][:2]
3349	if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
3350	except IndexError:
3351	continue
3352	if name in varylist:
3353	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
3354	if Ka2:
3355	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
3356	elif name in dependentVars:
3357	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
3358	if Ka2:
3359	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
3360	# now process derivatives in constraints
3361	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
3362	return dMdv
3363
3364	def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
3365	parmdict.update(zip(varylist,values))
3366	G2mv.Dict2Map(parmdict,varylist)
3367	Histograms,Phases,restraintDict = HistoPhases
3368	nvar = len(varylist)
3369	dMdv = np.empty(0)
3370	histoList = Histograms.keys()
3371	histoList.sort()
3372	for histogram in histoList:
3373	if 'PWDR' in histogram[:4]:
3374	Histogram = Histograms[histogram]
3375	hId = Histogram['hId']
3376	hfx = ':%d:'%(hId)
3377	wtFactor = calcControls[hfx+'wtFactor']
3378	Limits = calcControls[hfx+'Limits']
3379	x,y,w,yc,yb,yd = Histogram['Data']
3380	W = wtFactor*w
3381	xB = np.searchsorted(x,Limits[0])
3382	xF = np.searchsorted(x,Limits[1])
3383	dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
3384	varylist,Histogram,Phases,calcControls,pawleyLookup)
3385	elif 'HKLF' in histogram[:4]:
3386	Histogram = Histograms[histogram]
3387	nobs = Histogram['Nobs']
3388	phase = Histogram['Reflection Lists']
3389	Phase = Phases[phase]
3390	hId = Histogram['hId']
3391	hfx = ':%d:'%(hId)
3392	wtFactor = calcControls[hfx+'wtFactor']
3393	pfx = '%d::'%(Phase['pId'])
3394	phfx = '%d:%d:'%(Phase['pId'],hId)
3395	SGData = Phase['General']['SGData']
3396	A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3397	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3398	refList = Histogram['Data']
3399	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3400	dMdvh = np.zeros((len(varylist),len(refList)))
3401	for iref,ref in enumerate(refList):
3402	if ref[6] > 0:
3403	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
3404	if calcControls['F**2']:
3405	if ref[5]/ref[6] >= calcControls['minF/sig']:
3406	w = wtFactor/ref[6]
3407	for j,var in enumerate(varylist):
3408	if var in dFdvDict:
3409	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor
3410	if phfx+'Scale' in varylist:
3411	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
3412	else:
3413	Fo = np.sqrt(ref[5])
3414	Fc = np.sqrt(ref[7])
3415	sig = ref[6]/(2.0*Fo)
3416	if Fo/sig >= calcControls['minF/sig']:
3417	w = wtFactor/sig
3418	for j,var in enumerate(varylist):
3419	if var in dFdvDict:
3420	dMdvh[j][iref] = wdFdvDict[var][iref]np.sqrt(dervCor)
3421	if phfx+'Scale' in varylist:
3422	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]np.sqrt(dervCor)
3423	for item in ['Ep','Es','Eg']:
3424	if phfx+item in varylist:
3425	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]
3426	else:
3427	continue #skip non-histogram entries
3428	if len(dMdv):
3429	dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
3430	else:
3431	dMdv = dMdvh
3432
3433	pNames,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
3434	if np.any(pVals):
3435	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmdict,varylist)
3436	dMdv = np.concatenate((dMdv.T,dpdv.T)).T
3437
3438	return dMdv
3439
3440	def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
3441	parmdict.update(zip(varylist,values))
3442	G2mv.Dict2Map(parmdict,varylist)
3443	Histograms,Phases,restraintDict = HistoPhases
3444	nvar = len(varylist)
3445	Hess = np.empty(0)
3446	histoList = Histograms.keys()
3447	histoList.sort()
3448	for histogram in histoList:
3449	if 'PWDR' in histogram[:4]:
3450	Histogram = Histograms[histogram]
3451	hId = Histogram['hId']
3452	hfx = ':%d:'%(hId)
3453	wtFactor = calcControls[hfx+'wtFactor']
3454	Limits = calcControls[hfx+'Limits']
3455	x,y,w,yc,yb,yd = Histogram['Data']
3456	W = wtFactor*w
3457	dy = y-yc
3458	xB = np.searchsorted(x,Limits[0])
3459	xF = np.searchsorted(x,Limits[1])
3460	dMdvh = getPowderProfileDerv(parmdict,x[xB:xF],
3461	varylist,Histogram,Phases,calcControls,pawleyLookup)
3462	Wt = np.sqrt(W[xB:xF])[np.newaxis,:]
3463	Dy = dy[xB:xF][np.newaxis,:]
3464	dMdvh *= Wt
3465	if dlg:
3466	dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3467	if len(Hess):
3468	Hess += np.inner(dMdvh,dMdvh)
3469	dMdvh = WtDy
3470	Vec += np.sum(dMdvh,axis=1)
3471	else:
3472	Hess = np.inner(dMdvh,dMdvh)
3473	dMdvh = WtDy
3474	Vec = np.sum(dMdvh,axis=1)
3475	elif 'HKLF' in histogram[:4]:
3476	Histogram = Histograms[histogram]
3477	nobs = Histogram['Nobs']
3478	phase = Histogram['Reflection Lists']
3479	Phase = Phases[phase]
3480	hId = Histogram['hId']
3481	hfx = ':%d:'%(hId)
3482	wtFactor = calcControls[hfx+'wtFactor']
3483	pfx = '%d::'%(Phase['pId'])
3484	phfx = '%d:%d:'%(Phase['pId'],hId)
3485	SGData = Phase['General']['SGData']
3486	A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3487	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3488	refList = Histogram['Data']
3489	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3490	dMdvh = np.zeros((len(varylist),len(refList)))
3491	wdf = np.zeros(len(refList))
3492	for iref,ref in enumerate(refList):
3493	if ref[6] > 0:
3494	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
3495	if calcControls['F**2']:
3496	if ref[5]/ref[6] >= calcControls['minF/sig']:
3497	w = wtFactor/ref[6]
3498	wdf[iref] = w*(ref[5]-ref[7])
3499	for j,var in enumerate(varylist):
3500	if var in dFdvDict:
3501	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor
3502	if phfx+'Scale' in varylist:
3503	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
3504	else:
3505	if ref[5] > 0.:
3506	Fo = np.sqrt(ref[5])
3507	Fc = np.sqrt(ref[7])
3508	sig = ref[6]/(2.0*Fo)
3509	w = wtFactor/sig
3510	wdf[iref] = w*(Fo-Fc)
3511	if Fo/sig >= calcControls['minF/sig']:
3512	for j,var in enumerate(varylist):
3513	if var in dFdvDict:
3514	dMdvh[j][iref] = wdFdvDict[var][iref]np.sqrt(dervCor)
3515	if phfx+'Scale' in varylist:
3516	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]np.sqrt(dervCor)
3517	for item in ['Ep','Es','Eg']:
3518	if phfx+item in varylist:
3519	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]
3520	if dlg:
3521	dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3522	if len(Hess):
3523	Vec += np.sum(dMdvh*wdf,axis=1)
3524	Hess += np.inner(dMdvh,dMdvh)
3525	else:
3526	Vec = np.sum(dMdvh*wdf,axis=1)
3527	Hess = np.inner(dMdvh,dMdvh)
3528	else:
3529	continue #skip non-histogram entries
3530	pNames,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
3531	if np.any(pVals):
3532	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmdict,varylist)
3533	Vec += np.sum(dpdvpWtpVals,axis=1)
3534	Hess += np.inner(dpdv*pWt,dpdv)
3535	return Vec,Hess
3536
3537	def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
3538	parmdict.update(zip(varylist,values))
3539	Values2Dict(parmdict, varylist, values)
3540	G2mv.Dict2Map(parmdict,varylist)
3541	Histograms,Phases,restraintDict = HistoPhases
3542	M = np.empty(0)
3543	SumwYo = 0
3544	Nobs = 0
3545	histoList = Histograms.keys()
3546	histoList.sort()
3547	for histogram in histoList:
3548	if 'PWDR' in histogram[:4]:
3549	Histogram = Histograms[histogram]
3550	hId = Histogram['hId']
3551	hfx = ':%d:'%(hId)
3552	wtFactor = calcControls[hfx+'wtFactor']
3553	Limits = calcControls[hfx+'Limits']
3554	x,y,w,yc,yb,yd = Histogram['Data']
3555	W = wtFactor*w
3556	yc *= 0.0 #zero full calcd profiles
3557	yb *= 0.0
3558	yd *= 0.0
3559	xB = np.searchsorted(x,Limits[0])
3560	xF = np.searchsorted(x,Limits[1])
3561	Histogram['Nobs'] = xF-xB
3562	Nobs += Histogram['Nobs']
3563	Histogram['sumwYo'] = np.sum(W[xB:xF]y[xB:xF]*2)
3564	SumwYo += Histogram['sumwYo']
3565	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
3566	varylist,Histogram,Phases,calcControls,pawleyLookup)
3567	yc[xB:xF] += yb[xB:xF]
3568	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
3569	Histogram['sumwYd'] = np.sum(np.sqrt(W[xB:xF])*(yd[xB:xF]))
3570	wdy = -np.sqrt(W[xB:xF])*(yd[xB:xF])
3571	Histogram['wR'] = min(100.,np.sqrt(np.sum(wdy*2)/Histogram['sumwYo'])100.)
3572	if dlg:
3573	dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3574	M = np.concatenate((M,wdy))
3575	#end of PWDR processing
3576	elif 'HKLF' in histogram[:4]:
3577	Histogram = Histograms[histogram]
3578	phase = Histogram['Reflection Lists']
3579	Phase = Phases[phase]
3580	hId = Histogram['hId']
3581	hfx = ':%d:'%(hId)
3582	wtFactor = calcControls[hfx+'wtFactor']
3583	pfx = '%d::'%(Phase['pId'])
3584	phfx = '%d:%d:'%(Phase['pId'],hId)
3585	SGData = Phase['General']['SGData']
3586	A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3587	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3588	refList = Histogram['Data']
3589	refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3590	df = np.zeros(len(refList))
3591	sumwYo = 0
3592	sumFo = 0
3593	sumFo2 = 0
3594	sumdF = 0
3595	sumdF2 = 0
3596	nobs = 0
3597	for i,ref in enumerate(refList):
3598	if ref[6] > 0:
3599	SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
3600	ref[7] = parmdict[phfx+'Scale']*ref[9]
3601	ref[7] *= ref[13]
3602	ref[8] = ref[5]/parmdict[phfx+'Scale']
3603	if calcControls['F**2']:
3604	if ref[5]/ref[6] >= calcControls['minF/sig']:
3605	sumFo2 += ref[5]
3606	Fo = np.sqrt(ref[5])
3607	sumFo += Fo
3608	sumFo2 += ref[5]
3609	sumdF += abs(Fo-np.sqrt(ref[7]))
3610	sumdF2 += abs(ref[5]-ref[7])
3611	nobs += 1
3612	df[i] = -np.sqrt(wtFactor)*(ref[5]-ref[7])/ref[6]
3613	sumwYo += wtFactor(ref[5]/ref[6])*2
3614	else:
3615	Fo = np.sqrt(ref[5])
3616	Fc = np.sqrt(ref[7])
3617	sig = ref[6]/(2.0*Fo)
3618	if Fo/sig >= calcControls['minF/sig']:
3619	sumFo += Fo
3620	sumFo2 += ref[5]
3621	sumdF += abs(Fo-Fc)
3622	sumdF2 += abs(ref[5]-ref[7])
3623	nobs += 1
3624	df[i] = -np.sqrt(wtFactor)*(Fo-Fc)/sig
3625	sumwYo += wtFactor(Fo/sig)*2
3626	Histogram['Nobs'] = nobs
3627	Histogram['sumwYo'] = sumwYo
3628	SumwYo += sumwYo
3629	Histogram['wR'] = min(100.,np.sqrt(np.sum(df*2)/Histogram['sumwYo'])100.)
3630	Histogram[phfx+'Rf'] = 100.*sumdF/sumFo
3631	Histogram[phfx+'Rf^2'] = 100.*sumdF2/sumFo2
3632	Histogram[phfx+'Nref'] = nobs
3633	Nobs += nobs
3634	if dlg:
3635	dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3636	M = np.concatenate((M,df))
3637	# end of HKLF processing
3638	Histograms['sumwYo'] = SumwYo
3639	Histograms['Nobs'] = Nobs
3640	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
3641	if dlg:
3642	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
3643	if not GoOn:
3644	parmdict['saved values'] = values
3645	dlg.Destroy()
3646	raise Exception #Abort!!
3647	pDict,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
3648	if np.any(pVals):
3649	pSum = np.sum(pWtpVals*2)
3650	print 'Penalty function: %.3f on %d terms'%(pSum,len(pVals))
3651	Nobs += len(pVals)
3652	M = np.concatenate((M,np.sqrt(pWt)*pVals))
3653	return M
3654
3655	def Refine(GPXfile,dlg):
3656	import pytexture as ptx
3657	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
3658
3659	printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
3660	ShowBanner(printFile)
3661	varyList = []
3662	parmDict = {}
3663	G2mv.InitVars()
3664	Controls = GetControls(GPXfile)
3665	ShowControls(Controls,printFile)
3666	calcControls = {}
3667	calcControls.update(Controls)
3668	constrDict,fixedList = GetConstraints(GPXfile)
3669	restraintDict = GetRestraints(GPXfile)
3670	rigidbodyDict = GetRigidBodies(GPXfile)
3671	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3672	if not Phases:
3673	print ' * ERROR - you have no phases! *'
3674	print ' * Refine aborted *'
3675	raise Exception
3676	if not Histograms:
3677	print ' * ERROR - you have no data to refine with! *'
3678	print ' * Refine aborted *'
3679	raise Exception
3680	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,restraintDict,pFile=printFile)
3681	calcControls['atomIndx'] = atomIndx
3682	calcControls['Natoms'] = Natoms
3683	calcControls['FFtables'] = FFtables
3684	calcControls['BLtables'] = BLtables
3685	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
3686	calcControls.update(controlDict)
3687	histVary,histDict,controlDict = GetHistogramData(Histograms,pFile=printFile)
3688	calcControls.update(controlDict)
3689	varyList = phaseVary+hapVary+histVary
3690	parmDict.update(phaseDict)
3691	parmDict.update(hapDict)
3692	parmDict.update(histDict)
3693	GetFprime(calcControls,Histograms)
3694	# do constraint processing
3695	try:
3696	groups,parmlist = G2mv.GroupConstraints(constrDict)
3697	G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3698	except:
3699	print ' * ERROR - your constraints are internally inconsistent *'
3700	# traceback for debug
3701	#print 'varyList',varyList
3702	#print 'constrDict',constrDict
3703	#print 'fixedList',fixedList
3704	#import traceback
3705	#print traceback.format_exc()
3706	raise Exception(' * Refine aborted *')
3707	# # check to see which generated parameters are fully varied
3708	# msg = G2mv.SetVaryFlags(varyList)
3709	# if msg:
3710	# print ' * ERROR - you have not set the refine flags for constraints consistently! *'
3711	# print msg
3712	# raise Exception(' * Refine aborted *')
3713	#print G2mv.VarRemapShow(varyList)
3714	G2mv.Map2Dict(parmDict,varyList)
3715	Rvals = {}
3716	while True:
3717	begin = time.time()
3718	values = np.array(Dict2Values(parmDict, varyList))
3719	Ftol = Controls['min dM/M']
3720	Factor = Controls['shift factor']
3721	maxCyc = Controls['max cyc']
3722	if 'Jacobian' in Controls['deriv type']:
3723	result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3724	ftol=Ftol,col_deriv=True,factor=Factor,
3725	args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3726	ncyc = int(result[2]['nfev']/2)
3727	elif 'Hessian' in Controls['deriv type']:
3728	result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3729	args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3730	ncyc = result[2]['num cyc']+1
3731	Rvals['lamMax'] = result[2]['lamMax']
3732	else: #'numeric'
3733	result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3734	args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3735	ncyc = int(result[2]['nfev']/len(varyList))
3736	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
3737	# for item in table: print item,table[item] #useful debug - are things shifting?
3738	runtime = time.time()-begin
3739	Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3740	Values2Dict(parmDict, varyList, result[0])
3741	G2mv.Dict2Map(parmDict,varyList)
3742
3743	Rvals['Nobs'] = Histograms['Nobs']
3744	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100. #to %
3745	Rvals['GOF'] = Rvals['chisq']/(Histograms['Nobs']-len(varyList))
3746	print >>printFile,'\n Refinement results:'
3747	print >>printFile,135*'-'
3748	print >>printFile,' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
3749	print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3750	print >>printFile,' wR = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3751	try:
3752	covMatrix = result[1]*Rvals['GOF']
3753	sig = np.sqrt(np.diag(covMatrix))
3754	if np.any(np.isnan(sig)):
3755	print '* Least squares aborted - some invalid esds possible *'
3756	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
3757	# for item in table: print item,table[item] #useful debug - are things shifting?
3758	break #refinement succeeded - finish up!
3759	except TypeError: #result[1] is None on singular matrix
3760	print '** Refinement failed - singular matrix **'
3761	if 'Hessian' in Controls['deriv type']:
3762	num = len(varyList)-1
3763	for i,val in enumerate(np.flipud(result[2]['psing'])):
3764	if val:
3765	print 'Removing parameter: ',varyList[num-i]
3766	del(varyList[num-i])
3767	else:
3768	Ipvt = result[2]['ipvt']
3769	for i,ipvt in enumerate(Ipvt):
3770	if not np.sum(result[2]['fjac'],axis=1)[i]:
3771	print 'Removing parameter: ',varyList[ipvt-1]
3772	del(varyList[ipvt-1])
3773	break
3774
3775	# print 'dependentParmList: ',G2mv.dependentParmList
3776	# print 'arrayList: ',G2mv.arrayList
3777	# print 'invarrayList: ',G2mv.invarrayList
3778	# print 'indParmList: ',G2mv.indParmList
3779	# print 'fixedDict: ',G2mv.fixedDict
3780	# print 'test1'
3781	GetFobsSq(Histograms,Phases,parmDict,calcControls)
3782	# print 'test2'
3783	sigDict = dict(zip(varyList,sig))
3784	newCellDict = GetNewCellParms(parmDict,varyList)
3785	newAtomDict = ApplyXYZshifts(parmDict,varyList)
3786	covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3787	'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3788	# add the uncertainties into the esd dictionary (sigDict)
3789	sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
3790	G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
3791	SetPhaseData(parmDict,sigDict,Phases,covData,restraintDict,printFile)
3792	SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,pFile=printFile)
3793	SetHistogramData(parmDict,sigDict,Histograms,pFile=printFile)
3794	SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
3795	printFile.close()
3796	print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
3797	print ' *** Refinement successful ***'
3798
3799	#for testing purposes!!!
3800	if DEBUG:
3801	import cPickle
3802	fl = open('structTestdata.dat','wb')
3803	cPickle.dump(parmDict,fl,1)
3804	cPickle.dump(varyList,fl,1)
3805	for histogram in Histograms:
3806	if 'PWDR' in histogram[:4]:
3807	Histogram = Histograms[histogram]
3808	cPickle.dump(Histogram,fl,1)
3809	cPickle.dump(Phases,fl,1)
3810	cPickle.dump(calcControls,fl,1)
3811	cPickle.dump(pawleyLookup,fl,1)
3812	fl.close()
3813
3814	if dlg:
3815	return Rvals['Rwp']
3816
3817	def SeqRefine(GPXfile,dlg):
3818	import pytexture as ptx
3819	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
3820
3821	printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
3822	print ' Sequential Refinement'
3823	ShowBanner(printFile)
3824	G2mv.InitVars()
3825	Controls = GetControls(GPXfile)
3826	ShowControls(Controls,printFile)
3827	constrDict,fixedList = GetConstraints(GPXfile)
3828	restraintDict = GetRestraints(GPXfile)
3829	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3830	if not Phases:
3831	print ' * ERROR - you have no histograms to refine! *'
3832	print ' * Refine aborted *'
3833	raise Exception
3834	if not Histograms:
3835	print ' * ERROR - you have no data to refine with! *'
3836	print ' * Refine aborted *'
3837	raise Exception
3838	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,restraintDict,False,printFile)
3839	for item in phaseVary:
3840	if '::A0' in item:
3841	print '** WARNING - lattice parameters should not be refined in a sequential refinement **'
3842	print '** instead use the Dij parameters for each powder histogram **'
3843	if 'Seq Data' in Controls:
3844	histNames = Controls['Seq Data']
3845	else:
3846	histNames = GetHistogramNames(GPXfile,['PWDR',])
3847	if 'Reverse Seq' in Controls:
3848	if Controls['Reverse Seq']:
3849	histNames.reverse()
3850	SeqResult = {'histNames':histNames}
3851	makeBack = True
3852	for ihst,histogram in enumerate(histNames):
3853	ifPrint = False
3854	if dlg:
3855	dlg.SetTitle('Residual for histogram '+str(ihst))
3856	calcControls = {}
3857	calcControls['atomIndx'] = atomIndx
3858	calcControls['Natoms'] = Natoms
3859	calcControls['FFtables'] = FFtables
3860	calcControls['BLtables'] = BLtables
3861	varyList = []
3862	parmDict = {}
3863	Histo = {histogram:Histograms[histogram],}
3864	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
3865	calcControls.update(controlDict)
3866	histVary,histDict,controlDict = GetHistogramData(Histo,False)
3867	calcControls.update(controlDict)
3868	varyList = phaseVary+hapVary+histVary
3869	if not ihst:
3870	saveVaryList = varyList[:]
3871	for i,item in enumerate(saveVaryList):
3872	items = item.split(':')
3873	if items[1]:
3874	items[1] = ''
3875	item = ':'.join(items)
3876	saveVaryList[i] = item
3877	SeqResult['varyList'] = saveVaryList
3878	else:
3879	newVaryList = varyList[:]
3880	for i,item in enumerate(newVaryList):
3881	items = item.split(':')
3882	if items[1]:
3883	items[1] = ''
3884	item = ':'.join(items)
3885	newVaryList[i] = item
3886	if newVaryList != SeqResult['varyList']:
3887	print newVaryList
3888	print SeqResult['varyList']
3889	print '**** ERROR - variable list for this histogram does not match previous'
3890	raise Exception
3891	parmDict.update(phaseDict)
3892	parmDict.update(hapDict)
3893	parmDict.update(histDict)
3894	GetFprime(calcControls,Histo)
3895	# do constraint processing
3896	try:
3897	groups,parmlist = G2mv.GroupConstraints(constrDict)
3898	G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3899	except:
3900	print ' * ERROR - your constraints are internally inconsistent *'
3901	raise Exception(' * Refine aborted *')
3902	# check to see which generated parameters are fully varied
3903	# msg = G2mv.SetVaryFlags(varyList)
3904	# if msg:
3905	# print ' * ERROR - you have not set the refine flags for constraints consistently! *'
3906	# print msg
3907	# raise Exception(' * Refine aborted *')
3908	#print G2mv.VarRemapShow(varyList)
3909	G2mv.Map2Dict(parmDict,varyList)
3910	Rvals = {}
3911	while True:
3912	begin = time.time()
3913	values = np.array(Dict2Values(parmDict,varyList))
3914	Ftol = Controls['min dM/M']
3915	Factor = Controls['shift factor']
3916	maxCyc = Controls['max cyc']
3917
3918	if 'Jacobian' in Controls['deriv type']:
3919	result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3920	ftol=Ftol,col_deriv=True,factor=Factor,
3921	args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3922	ncyc = int(result[2]['nfev']/2)
3923	elif 'Hessian' in Controls['deriv type']:
3924	result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3925	args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3926	ncyc = result[2]['num cyc']+1
3927	else: #'numeric'
3928	result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3929	args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3930	ncyc = int(result[2]['nfev']/len(varyList))
3931
3932	runtime = time.time()-begin
3933	Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3934	Values2Dict(parmDict, varyList, result[0])
3935	G2mv.Dict2Map(parmDict,varyList)
3936
3937	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histo['sumwYo'])*100. #to %
3938	Rvals['GOF'] = Rvals['Rwp']/(Histo['Nobs']-len(varyList))
3939	Rvals['Nobs'] = Histo['Nobs']
3940	print >>printFile,'\n Refinement results for histogram: v'+histogram
3941	print >>printFile,135*'-'
3942	print >>printFile,' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
3943	print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3944	print >>printFile,' wRp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3945	try:
3946	covMatrix = result[1]*Rvals['GOF']
3947	sig = np.sqrt(np.diag(covMatrix))
3948	if np.any(np.isnan(sig)):
3949	print '* Least squares aborted - some invalid esds possible *'
3950	ifPrint = True
3951	break #refinement succeeded - finish up!
3952	except TypeError: #result[1] is None on singular matrix
3953	print '** Refinement failed - singular matrix **'
3954	if 'Hessian' in Controls['deriv type']:
3955	num = len(varyList)-1
3956	for i,val in enumerate(np.flipud(result[2]['psing'])):
3957	if val:
3958	�