Context Navigation

source: trunk/GSASIIstruct.py @ 743

Last change on this file since 743 was 743, checked in by vondreele, 11 years ago
fix PO refinement output & calcs for SH
Property svn:keywords set to `Date Author Revision URL Id`
File size: 168.2 KB

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

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

Download in other formats: