Context Navigation

source: trunk/GSASIIstruct.py @ 874

Last change on this file since 874 was 874, checked in by vondreele, 10 years ago
add citation to banners
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 188.6 KB

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