Context Navigation

source: trunk/GSASIIstruct.py @ 635

Last change on this file since 635 was 635, checked in by vondreele, 13 years ago
HKLF refinement "numeric" now OK; revamp output a bit for more generalized R-values.
Property svn:keywords set to `Date Author Revision URL Id`
File size: 156.9 KB

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

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

Download in other formats: