Context Navigation

source: trunk/GSASIIstrIO.py @ 927

Last change on this file since 927 was 926, checked in by vondreele, 12 years ago
split GSASIIstruct.py into 3 pieces
File size: 107.2 KB

Line
1	# -- coding: utf-8 --
2	#GSASIIstrIO - structure computation IO routines
3	########### SVN repository information ###################
4	# $Date: 2013-05-17 12:13:15 -0500 (Fri, 17 May 2013) $
5	# $Author: vondreele $
6	# $Revision: 920 $
7	# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIstrIO.py $
8	# $Id: GSASIIstrIO.py 920 2013-05-17 17:13:15Z vondreele $
9	########### SVN repository information ###################
10	import sys
11	import os
12	import os.path as ospath
13	import time
14	import math
15	import copy
16	import random
17	import cPickle
18	import numpy as np
19	import numpy.ma as ma
20	import numpy.linalg as nl
21	import scipy.optimize as so
22	import GSASIIpath
23	GSASIIpath.SetVersionNumber("$Revision: 920 $")
24	import GSASIIElem as G2el
25	import GSASIIlattice as G2lat
26	import GSASIIspc as G2spc
27	import GSASIImapvars as G2mv
28	import GSASIImath as G2mth
29
30	sind = lambda x: np.sin(x*np.pi/180.)
31	cosd = lambda x: np.cos(x*np.pi/180.)
32	tand = lambda x: np.tan(x*np.pi/180.)
33	asind = lambda x: 180.*np.arcsin(x)/np.pi
34	acosd = lambda x: 180.*np.arccos(x)/np.pi
35	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
36
37	ateln2 = 8.0*math.log(2.0)
38
39	def GetControls(GPXfile):
40	''' Returns dictionary of control items found in GSASII gpx file
41	input:
42	GPXfile = .gpx full file name
43	return:
44	Controls = dictionary of control items
45	'''
46	Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
47	'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
48	fl = open(GPXfile,'rb')
49	while True:
50	try:
51	data = cPickle.load(fl)
52	except EOFError:
53	break
54	datum = data[0]
55	if datum[0] == 'Controls':
56	Controls.update(datum[1])
57	fl.close()
58	return Controls
59
60	def GetConstraints(GPXfile):
61	'''Read the constraints from the GPX file and interpret them
62	'''
63	constList = []
64	fl = open(GPXfile,'rb')
65	while True:
66	try:
67	data = cPickle.load(fl)
68	except EOFError:
69	break
70	datum = data[0]
71	if datum[0] == 'Constraints':
72	constDict = datum[1]
73	for item in constDict:
74	constList += constDict[item]
75	fl.close()
76	constDict,fixedList,ignored = ProcessConstraints(constList)
77	if ignored:
78	print ignored,'old-style Constraints were rejected'
79	return constDict,fixedList
80
81	def ProcessConstraints(constList):
82	"interpret constraints"
83	constDict = []
84	fixedList = []
85	ignored = 0
86	for item in constList:
87	if item[-1] == 'h':
88	# process a hold
89	fixedList.append('0')
90	constDict.append({item[0][1]:0.0})
91	elif item[-1] == 'f':
92	# process a new variable
93	fixedList.append(None)
94	constDict.append({})
95	for term in item[:-3]:
96	constDict[-1][term[1]] = term[0]
97	#constFlag[-1] = ['Vary']
98	elif item[-1] == 'c':
99	# process a contraint relationship
100	fixedList.append(str(item[-3]))
101	constDict.append({})
102	for term in item[:-3]:
103	constDict[-1][term[1]] = term[0]
104	#constFlag[-1] = ['VaryFree']
105	elif item[-1] == 'e':
106	# process an equivalence
107	firstmult = None
108	eqlist = []
109	for term in item[:-3]:
110	if term[0] == 0: term[0] = 1.0
111	if firstmult is None:
112	firstmult,firstvar = term
113	else:
114	eqlist.append([term[1],firstmult/term[0]])
115	G2mv.StoreEquivalence(firstvar,eqlist)
116	else:
117	ignored += 1
118	return constDict,fixedList,ignored
119
120	def CheckConstraints(GPXfile):
121	'''Load constraints and related info and return any error or warning messages'''
122	# init constraints
123	G2mv.InitVars()
124	# get variables
125	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
126	if not Phases:
127	return 'Error: No Phases!',''
128	if not Histograms:
129	return 'Error: no diffraction data',''
130	rigidbodyDict = GetRigidBodies(GPXfile)
131	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
132	rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
133	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)
134	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
135	histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
136	varyList = rbVary+phaseVary+hapVary+histVary
137	constrDict,fixedList = GetConstraints(GPXfile)
138	return G2mv.CheckConstraints(varyList,constrDict,fixedList)
139
140	def GetRestraints(GPXfile):
141	'''Read the restraints from the GPX file
142	'''
143	fl = open(GPXfile,'rb')
144	while True:
145	try:
146	data = cPickle.load(fl)
147	except EOFError:
148	break
149	datum = data[0]
150	if datum[0] == 'Restraints':
151	restraintDict = datum[1]
152	fl.close()
153	return restraintDict
154
155	def GetRigidBodies(GPXfile):
156	'''Read the rigid body models from the GPX file
157	'''
158	fl = open(GPXfile,'rb')
159	while True:
160	try:
161	data = cPickle.load(fl)
162	except EOFError:
163	break
164	datum = data[0]
165	if datum[0] == 'Rigid bodies':
166	rigidbodyDict = datum[1]
167	fl.close()
168	return rigidbodyDict
169
170	def GetFprime(controlDict,Histograms):
171	FFtables = controlDict['FFtables']
172	if not FFtables:
173	return
174	histoList = Histograms.keys()
175	histoList.sort()
176	for histogram in histoList:
177	if histogram[:4] in ['PWDR','HKLF']:
178	Histogram = Histograms[histogram]
179	hId = Histogram['hId']
180	hfx = ':%d:'%(hId)
181	keV = controlDict[hfx+'keV']
182	for El in FFtables:
183	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
184	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
185	FFtables[El][hfx+'FP'] = FP
186	FFtables[El][hfx+'FPP'] = FPP
187
188	def GetPhaseNames(GPXfile):
189	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
190	input:
191	GPXfile = gpx full file name
192	return:
193	PhaseNames = list of phase names
194	'''
195	fl = open(GPXfile,'rb')
196	PhaseNames = []
197	while True:
198	try:
199	data = cPickle.load(fl)
200	except EOFError:
201	break
202	datum = data[0]
203	if 'Phases' == datum[0]:
204	for datus in data[1:]:
205	PhaseNames.append(datus[0])
206	fl.close()
207	return PhaseNames
208
209	def GetAllPhaseData(GPXfile,PhaseName):
210	''' Returns the entire dictionary for PhaseName from GSASII gpx file
211	input:
212	GPXfile = gpx full file name
213	PhaseName = phase name
214	return:
215	phase dictionary
216	'''
217	fl = open(GPXfile,'rb')
218	General = {}
219	Atoms = []
220	while True:
221	try:
222	data = cPickle.load(fl)
223	except EOFError:
224	break
225	datum = data[0]
226	if 'Phases' == datum[0]:
227	for datus in data[1:]:
228	if datus[0] == PhaseName:
229	break
230	fl.close()
231	return datus[1]
232
233	def GetHistograms(GPXfile,hNames):
234	""" Returns a dictionary of histograms found in GSASII gpx file
235	input:
236	GPXfile = .gpx full file name
237	hNames = list of histogram names
238	return:
239	Histograms = dictionary of histograms (types = PWDR & HKLF)
240	"""
241	fl = open(GPXfile,'rb')
242	Histograms = {}
243	while True:
244	try:
245	data = cPickle.load(fl)
246	except EOFError:
247	break
248	datum = data[0]
249	hist = datum[0]
250	if hist in hNames:
251	if 'PWDR' in hist[:4]:
252	PWDRdata = {}
253	try:
254	PWDRdata.update(datum[1][0]) #weight factor
255	except ValueError:
256	PWDRdata['wtFactor'] = 1.0 #patch
257	PWDRdata['Data'] = datum[1][1] #powder data arrays
258	PWDRdata[data[2][0]] = data[2][1] #Limits
259	PWDRdata[data[3][0]] = data[3][1] #Background
260	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
261	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
262	try:
263	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
264	except IndexError:
265	PWDRdata['Reflection Lists'] = {}
266
267	Histograms[hist] = PWDRdata
268	elif 'HKLF' in hist[:4]:
269	HKLFdata = {}
270	try:
271	HKLFdata.update(datum[1][0]) #weight factor
272	except ValueError:
273	HKLFdata['wtFactor'] = 1.0 #patch
274	HKLFdata['Data'] = datum[1][1]
275	HKLFdata[data[1][0]] = data[1][1] #Instrument parameters
276	HKLFdata['Reflection Lists'] = None
277	Histograms[hist] = HKLFdata
278	fl.close()
279	return Histograms
280
281	def GetHistogramNames(GPXfile,hType):
282	""" Returns a list of histogram names found in GSASII gpx file
283	input:
284	GPXfile = .gpx full file name
285	hType = list ['PWDR','HKLF']
286	return:
287	HistogramNames = list of histogram names (types = PWDR & HKLF)
288	"""
289	fl = open(GPXfile,'rb')
290	HistogramNames = []
291	while True:
292	try:
293	data = cPickle.load(fl)
294	except EOFError:
295	break
296	datum = data[0]
297	if datum[0][:4] in hType:
298	HistogramNames.append(datum[0])
299	fl.close()
300	return HistogramNames
301
302	def GetUsedHistogramsAndPhases(GPXfile):
303	''' Returns all histograms that are found in any phase
304	and any phase that uses a histogram
305	input:
306	GPXfile = .gpx full file name
307	return:
308	Histograms = dictionary of histograms as {name:data,...}
309	Phases = dictionary of phases that use histograms
310	'''
311	phaseNames = GetPhaseNames(GPXfile)
312	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
313	allHistograms = GetHistograms(GPXfile,histoList)
314	phaseData = {}
315	for name in phaseNames:
316	phaseData[name] = GetAllPhaseData(GPXfile,name)
317	Histograms = {}
318	Phases = {}
319	for phase in phaseData:
320	Phase = phaseData[phase]
321	if Phase['Histograms']:
322	if phase not in Phases:
323	pId = phaseNames.index(phase)
324	Phase['pId'] = pId
325	Phases[phase] = Phase
326	for hist in Phase['Histograms']:
327	if 'Use' not in Phase['Histograms'][hist]: #patch
328	Phase['Histograms'][hist]['Use'] = True
329	if hist not in Histograms and Phase['Histograms'][hist]['Use']:
330	Histograms[hist] = allHistograms[hist]
331	hId = histoList.index(hist)
332	Histograms[hist]['hId'] = hId
333	return Histograms,Phases
334
335	def getBackupName(GPXfile,makeBack):
336	GPXpath,GPXname = ospath.split(GPXfile)
337	if GPXpath == '': GPXpath = '.'
338	Name = ospath.splitext(GPXname)[0]
339	files = os.listdir(GPXpath)
340	last = 0
341	for name in files:
342	name = name.split('.')
343	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
344	if makeBack:
345	last = max(last,int(name[1].strip('bak'))+1)
346	else:
347	last = max(last,int(name[1].strip('bak')))
348	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
349	return GPXback
350
351	def GPXBackup(GPXfile,makeBack=True):
352	import distutils.file_util as dfu
353	GPXback = getBackupName(GPXfile,makeBack)
354	dfu.copy_file(GPXfile,GPXback)
355	return GPXback
356
357	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
358	''' Updates gpxfile from all histograms that are found in any phase
359	and any phase that used a histogram. Also updates rigid body definitions.
360	input:
361	GPXfile = .gpx full file name
362	Histograms = dictionary of histograms as {name:data,...}
363	Phases = dictionary of phases that use histograms
364	RigidBodies = dictionary of rigid bodies
365	CovData = dictionary of refined variables, varyList, & covariance matrix
366	makeBack = True if new backup of .gpx file is to be made; else use the last one made
367	'''
368
369	GPXback = GPXBackup(GPXfile,makeBack)
370	print 'Read from file:',GPXback
371	print 'Save to file :',GPXfile
372	infile = open(GPXback,'rb')
373	outfile = open(GPXfile,'wb')
374	while True:
375	try:
376	data = cPickle.load(infile)
377	except EOFError:
378	break
379	datum = data[0]
380	# print 'read: ',datum[0]
381	if datum[0] == 'Phases':
382	for iphase in range(len(data)):
383	if data[iphase][0] in Phases:
384	phaseName = data[iphase][0]
385	data[iphase][1].update(Phases[phaseName])
386	elif datum[0] == 'Covariance':
387	data[0][1] = CovData
388	elif datum[0] == 'Rigid bodies':
389	data[0][1] = RigidBodies
390	try:
391	histogram = Histograms[datum[0]]
392	# print 'found ',datum[0]
393	data[0][1][1] = histogram['Data']
394	for datus in data[1:]:
395	# print ' read: ',datus[0]
396	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
397	datus[1] = histogram[datus[0]]
398	except KeyError:
399	pass
400
401	cPickle.dump(data,outfile,1)
402	infile.close()
403	outfile.close()
404	print 'GPX file save successful'
405
406	def SetSeqResult(GPXfile,Histograms,SeqResult):
407	GPXback = GPXBackup(GPXfile)
408	print 'Read from file:',GPXback
409	print 'Save to file :',GPXfile
410	infile = open(GPXback,'rb')
411	outfile = open(GPXfile,'wb')
412	while True:
413	try:
414	data = cPickle.load(infile)
415	except EOFError:
416	break
417	datum = data[0]
418	if datum[0] == 'Sequential results':
419	data[0][1] = SeqResult
420	try:
421	histogram = Histograms[datum[0]]
422	data[0][1][1] = histogram['Data']
423	for datus in data[1:]:
424	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
425	datus[1] = histogram[datus[0]]
426	except KeyError:
427	pass
428
429	cPickle.dump(data,outfile,1)
430	infile.close()
431	outfile.close()
432	print 'GPX file save successful'
433
434	def ShowBanner(pFile=None):
435	print >>pFile,80''
436	print >>pFile,' General Structure Analysis System-II Crystal Structure Refinement'
437	print >>pFile,' by Robert B. Von Dreele & Brian H. Toby'
438	print >>pFile,' Argonne National Laboratory(C), 2010'
439	print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,'
440	print >>pFile,' as Operator of Argonne National Laboratory.'
441	print >>pFile,' Please cite:'
442	print >>pFile,' B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)'
443
444	print >>pFile,80'','\n'
445
446	def ShowControls(Controls,pFile=None):
447	print >>pFile,' Least squares controls:'
448	print >>pFile,' Refinement type: ',Controls['deriv type']
449	if 'Hessian' in Controls['deriv type']:
450	print >>pFile,' Maximum number of cycles:',Controls['max cyc']
451	else:
452	print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
453	print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
454
455	def GetFFtable(General):
456	''' returns a dictionary of form factor data for atom types found in General
457	input:
458	General = dictionary of phase info.; includes AtomTypes
459	return:
460	FFtable = dictionary of form factor data; key is atom type
461	'''
462	atomTypes = General['AtomTypes']
463	FFtable = {}
464	for El in atomTypes:
465	FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
466	for item in FFs:
467	if item['Symbol'] == El.upper():
468	FFtable[El] = item
469	return FFtable
470
471	def GetBLtable(General):
472	''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
473	input:
474	General = dictionary of phase info.; includes AtomTypes & Isotopes
475	return:
476	BLtable = dictionary of scattering length data; key is atom type
477	'''
478	atomTypes = General['AtomTypes']
479	BLtable = {}
480	isotopes = General['Isotopes']
481	isotope = General['Isotope']
482	for El in atomTypes:
483	BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
484	return BLtable
485
486	def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
487	# if SGLaue in ['-1','2/m','mmm']:
488	# return #no Pawley symmetry required constraints
489	eqvDict = {}
490	for i,varyI in enumerate(pawleyVary):
491	eqvDict[varyI] = []
492	refI = int(varyI.split(':')[-1])
493	ih,ik,il = PawleyRef[refI][:3]
494	dspI = PawleyRef[refI][4]
495	for varyJ in pawleyVary[i+1:]:
496	refJ = int(varyJ.split(':')[-1])
497	jh,jk,jl = PawleyRef[refJ][:3]
498	dspJ = PawleyRef[refJ][4]
499	if SGLaue in ['4/m','4/mmm']:
500	isum = ih2+ik2
501	jsum = jh2+jk2
502	if abs(il) == abs(jl) and isum == jsum:
503	eqvDict[varyI].append(varyJ)
504	elif SGLaue in ['3R','3mR']:
505	isum = ih2+ik2+il**2
506	jsum = jh2+jk2jl*2
507	isum2 = ihik+ihil+ik*il
508	jsum2 = jhjk+jhjl+jk*jl
509	if isum == jsum and isum2 == jsum2:
510	eqvDict[varyI].append(varyJ)
511	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
512	isum = ih2+ik2+ih*ik
513	jsum = jh2+jk2+jh*jk
514	if abs(il) == abs(jl) and isum == jsum:
515	eqvDict[varyI].append(varyJ)
516	elif SGLaue in ['m3','m3m']:
517	isum = ih2+ik2+il**2
518	jsum = jh2+jk2+jl**2
519	if isum == jsum:
520	eqvDict[varyI].append(varyJ)
521	elif abs(dspI-dspJ)/dspI < 1.e-4:
522	eqvDict[varyI].append(varyJ)
523	for item in pawleyVary:
524	if eqvDict[item]:
525	for item2 in pawleyVary:
526	if item2 in eqvDict[item]:
527	eqvDict[item2] = []
528	G2mv.StoreEquivalence(item,eqvDict[item])
529
530	def cellVary(pfx,SGData):
531	if SGData['SGLaue'] in ['-1',]:
532	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
533	elif SGData['SGLaue'] in ['2/m',]:
534	if SGData['SGUniq'] == 'a':
535	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
536	elif SGData['SGUniq'] == 'b':
537	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
538	else:
539	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
540	elif SGData['SGLaue'] in ['mmm',]:
541	return [pfx+'A0',pfx+'A1',pfx+'A2']
542	elif SGData['SGLaue'] in ['4/m','4/mmm']:
543	return [pfx+'A0',pfx+'A2']
544	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
545	return [pfx+'A0',pfx+'A2']
546	elif SGData['SGLaue'] in ['3R', '3mR']:
547	return [pfx+'A0',pfx+'A3']
548	elif SGData['SGLaue'] in ['m3m','m3']:
549	return [pfx+'A0',]
550
551	################################################################################
552	##### Rigid Body Models
553	################################################################################
554
555	def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
556
557	def PrintResRBModel(RBModel):
558	atNames = RBModel['atNames']
559	rbRef = RBModel['rbRef']
560	rbSeq = RBModel['rbSeq']
561	print >>pFile,'Residue RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
562	'No. times used: ',RBModel['useCount']
563	print >>pFile,' At name x y z'
564	for name,xyz in zip(atNames,RBModel['rbXYZ']):
565	print >>pFile,' %8s %10.4f %10.4f %10.4f'%(name,xyz[0],xyz[1],xyz[2])
566	print >>pFile,'Orientation defined by:',atNames[rbRef[0]],' -> ',atNames[rbRef[1]], \
567	' & ',atNames[rbRef[0]],' -> ',atNames[rbRef[2]]
568	if rbSeq:
569	for i,rbseq in enumerate(rbSeq):
570	print >>pFile,'Torsion sequence ',i,' Bond: '+atNames[rbseq[0]],' - ', \
571	atNames[rbseq[1]],' riding: ',[atNames[i] for i in rbseq[3]]
572
573	def PrintVecRBModel(RBModel):
574	rbRef = RBModel['rbRef']
575	atTypes = RBModel['rbTypes']
576	print >>pFile,'Vector RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
577	'No. times used: ',RBModel['useCount']
578	for i in range(len(RBModel['VectMag'])):
579	print >>pFile,'Vector no.: ',i,' Magnitude: ', \
580	'%8.4f'%(RBModel['VectMag'][i]),' Refine? ',RBModel['VectRef'][i]
581	print >>pFile,' No. Type vx vy vz'
582	for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
583	print >>pFile,' %d %2s %10.4f %10.4f %10.4f'%(j,name,xyz[0],xyz[1],xyz[2])
584	print >>pFile,' No. Type x y z'
585	for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
586	print >>pFile,' %d %2s %10.4f %10.4f %10.4f'%(i,name,xyz[0],xyz[1],xyz[2])
587	print >>pFile,'Orientation defined by: atom ',rbRef[0],' -> atom ',rbRef[1], \
588	' & atom ',rbRef[0],' -> atom ',rbRef[2]
589	rbVary = []
590	rbDict = {}
591	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
592	if len(rbIds['Vector']):
593	for irb,item in enumerate(rbIds['Vector']):
594	if rigidbodyDict['Vector'][item]['useCount']:
595	RBmags = rigidbodyDict['Vector'][item]['VectMag']
596	RBrefs = rigidbodyDict['Vector'][item]['VectRef']
597	for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
598	pid = '::RBV;'+str(i)+':'+str(irb)
599	rbDict[pid] = mag
600	if ref:
601	rbVary.append(pid)
602	if Print:
603	print >>pFile,'\nVector rigid body model:'
604	PrintVecRBModel(rigidbodyDict['Vector'][item])
605	if len(rbIds['Residue']):
606	for item in rbIds['Residue']:
607	if rigidbodyDict['Residue'][item]['useCount']:
608	if Print:
609	print >>pFile,'\nResidue rigid body model:'
610	PrintResRBModel(rigidbodyDict['Residue'][item])
611	return rbVary,rbDict
612
613	def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
614
615	def PrintRBVectandSig(VectRB,VectSig):
616	print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
617	namstr = ' names :'
618	valstr = ' values:'
619	sigstr = ' esds :'
620	for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
621	namstr += '%12s'%('Vect '+str(i))
622	valstr += '%12.4f'%(val)
623	if sig:
624	sigstr += '%12.4f'%(sig)
625	else:
626	sigstr += 12*' '
627	print >>pFile,namstr
628	print >>pFile,valstr
629	print >>pFile,sigstr
630
631	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
632	if not RBIds['Vector']:
633	return
634	for irb,item in enumerate(RBIds['Vector']):
635	if rigidbodyDict['Vector'][item]['useCount']:
636	VectSig = []
637	RBmags = rigidbodyDict['Vector'][item]['VectMag']
638	for i,mag in enumerate(RBmags):
639	name = '::RBV;'+str(i)+':'+str(irb)
640	mag = parmDict[name]
641	if name in sigDict:
642	VectSig.append(sigDict[name])
643	PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)
644
645	################################################################################
646	##### Phase data
647	################################################################################
648
649	def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None):
650
651	def PrintFFtable(FFtable):
652	print >>pFile,'\n X-ray scattering factors:'
653	print >>pFile,' Symbol fa fb fc'
654	print >>pFile,99*'-'
655	for Ename in FFtable:
656	ffdata = FFtable[Ename]
657	fa = ffdata['fa']
658	fb = ffdata['fb']
659	print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \
660	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
661
662	def PrintBLtable(BLtable):
663	print >>pFile,'\n Neutron scattering factors:'
664	print >>pFile,' Symbol isotope mass b resonant terms'
665	print >>pFile,99*'-'
666	for Ename in BLtable:
667	bldata = BLtable[Ename]
668	isotope = bldata[0]
669	mass = bldata[1][0]
670	blen = bldata[1][1]
671	bres = []
672	if len(bldata[1]) > 2:
673	bres = bldata[1][2:]
674	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
675	for item in bres:
676	line += '%10.5g'%(item)
677	print >>pFile,line
678
679	def PrintRBObjects(resRBData,vecRBData):
680
681	def PrintRBThermals():
682	tlstr = ['11','22','33','12','13','23']
683	sstr = ['12','13','21','23','31','32','AA','BB']
684	TLS = RB['ThermalMotion'][1]
685	TLSvar = RB['ThermalMotion'][2]
686	if 'T' in RB['ThermalMotion'][0]:
687	print >>pFile,'TLS data'
688	text = ''
689	for i in range(6):
690	text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
691	print >>pFile,text
692	if 'L' in RB['ThermalMotion'][0]:
693	text = ''
694	for i in range(6,12):
695	text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
696	print >>pFile,text
697	if 'S' in RB['ThermalMotion'][0]:
698	text = ''
699	for i in range(12,20):
700	text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
701	print >>pFile,text
702	if 'U' in RB['ThermalMotion'][0]:
703	print >>pFile,'Uiso data'
704	text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])
705
706	if len(resRBData):
707	for RB in resRBData:
708	Oxyz = RB['Orig'][0]
709	Qrijk = RB['Orient'][0]
710	Angle = 2.0*acosd(Qrijk[0])
711	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
712	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
713	print >>pFile,'Orientation angle,vector:', \
714	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
715	Torsions = RB['Torsions']
716	if len(Torsions):
717	text = 'Torsions: '
718	for torsion in Torsions:
719	text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
720	print >>pFile,text
721	PrintRBThermals()
722	if len(vecRBData):
723	for RB in vecRBData:
724	Oxyz = RB['Orig'][0]
725	Qrijk = RB['Orient'][0]
726	Angle = 2.0*acosd(Qrijk[0])
727	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
728	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
729	print >>pFile,'Orientation angle,vector:', \
730	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
731	PrintRBThermals()
732
733	def PrintAtoms(General,Atoms):
734	cx,ct,cs,cia = General['AtomPtrs']
735	print >>pFile,'\n Atoms:'
736	line = ' name type refine? x y z '+ \
737	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
738	if General['Type'] == 'magnetic':
739	line += ' Mx My Mz'
740	elif General['Type'] == 'macromolecular':
741	line = ' res no residue chain'+line
742	print >>pFile,line
743	if General['Type'] == 'nuclear':
744	print >>pFile,135*'-'
745	for i,at in enumerate(Atoms):
746	line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
747	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
748	if at[cia] == 'I':
749	line += '%8.4f'%(at[cia+1])+48*' '
750	else:
751	line += 8*' '
752	for j in range(6):
753	line += '%8.4f'%(at[cia+1+j])
754	print >>pFile,line
755	elif General['Type'] == 'macromolecular':
756	print >>pFile,135*'-'
757	for i,at in enumerate(Atoms):
758	line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
759	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
760	if at[cia] == 'I':
761	line += '%8.4f'%(at[cia+1])+48*' '
762	else:
763	line += 8*' '
764	for j in range(6):
765	line += '%8.4f'%(at[cia+1+j])
766	print >>pFile,line
767
768	def PrintTexture(textureData):
769	topstr = '\n Spherical harmonics texture: Order:' + \
770	str(textureData['Order'])
771	if textureData['Order']:
772	print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
773	else:
774	print >>pFile,topstr
775	return
776	names = ['omega','chi','phi']
777	line = '\n'
778	for name in names:
779	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
780	print >>pFile,line
781	print >>pFile,'\n Texture coefficients:'
782	ptlbls = ' names :'
783	ptstr = ' values:'
784	SHcoeff = textureData['SH Coeff'][1]
785	for item in SHcoeff:
786	ptlbls += '%12s'%(item)
787	ptstr += '%12.4f'%(SHcoeff[item])
788	print >>pFile,ptlbls
789	print >>pFile,ptstr
790
791	def MakeRBParms(rbKey,phaseVary,phaseDict):
792	rbid = str(rbids.index(RB['RBId']))
793	pfxRB = pfx+'RB'+rbKey+'P'
794	pstr = ['x','y','z']
795	ostr = ['a','i','j','k']
796	for i in range(3):
797	name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
798	phaseDict[name] = RB['Orig'][0][i]
799	if RB['Orig'][1]:
800	phaseVary += [name,]
801	pfxRB = pfx+'RB'+rbKey+'O'
802	for i in range(4):
803	name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
804	phaseDict[name] = RB['Orient'][0][i]
805	if RB['Orient'][1] == 'AV' and i:
806	phaseVary += [name,]
807	elif RB['Orient'][1] == 'A' and not i:
808	phaseVary += [name,]
809
810	def MakeRBThermals(rbKey,phaseVary,phaseDict):
811	rbid = str(rbids.index(RB['RBId']))
812	tlstr = ['11','22','33','12','13','23']
813	sstr = ['12','13','21','23','31','32','AA','BB']
814	if 'T' in RB['ThermalMotion'][0]:
815	pfxRB = pfx+'RB'+rbKey+'T'
816	for i in range(6):
817	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
818	phaseDict[name] = RB['ThermalMotion'][1][i]
819	if RB['ThermalMotion'][2][i]:
820	phaseVary += [name,]
821	if 'L' in RB['ThermalMotion'][0]:
822	pfxRB = pfx+'RB'+rbKey+'L'
823	for i in range(6):
824	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
825	phaseDict[name] = RB['ThermalMotion'][1][i+6]
826	if RB['ThermalMotion'][2][i+6]:
827	phaseVary += [name,]
828	if 'S' in RB['ThermalMotion'][0]:
829	pfxRB = pfx+'RB'+rbKey+'S'
830	for i in range(8):
831	name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
832	phaseDict[name] = RB['ThermalMotion'][1][i+12]
833	if RB['ThermalMotion'][2][i+12]:
834	phaseVary += [name,]
835	if 'U' in RB['ThermalMotion'][0]:
836	name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
837	phaseDict[name] = RB['ThermalMotion'][1][0]
838	if RB['ThermalMotion'][2][0]:
839	phaseVary += [name,]
840
841	def MakeRBTorsions(rbKey,phaseVary,phaseDict):
842	rbid = str(rbids.index(RB['RBId']))
843	pfxRB = pfx+'RB'+rbKey+'Tr;'
844	for i,torsion in enumerate(RB['Torsions']):
845	name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
846	phaseDict[name] = torsion[0]
847	if torsion[1]:
848	phaseVary += [name,]
849
850	if Print:
851	print >>pFile,'\n Phases:'
852	phaseVary = []
853	phaseDict = {}
854	phaseConstr = {}
855	pawleyLookup = {}
856	FFtables = {} #scattering factors - xrays
857	BLtables = {} # neutrons
858	Natoms = {}
859	AtMults = {}
860	AtIA = {}
861	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
862	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
863	atomIndx = {}
864	for name in PhaseData:
865	General = PhaseData[name]['General']
866	pId = PhaseData[name]['pId']
867	pfx = str(pId)+'::'
868	FFtable = GetFFtable(General)
869	BLtable = GetBLtable(General)
870	FFtables.update(FFtable)
871	BLtables.update(BLtable)
872	Atoms = PhaseData[name]['Atoms']
873	AtLookup = G2mth.FillAtomLookUp(Atoms)
874	PawleyRef = PhaseData[name].get('Pawley ref',[])
875	SGData = General['SGData']
876	SGtext = G2spc.SGPrint(SGData)
877	cell = General['Cell']
878	A = G2lat.cell2A(cell[1:7])
879	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
880	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
881	if cell[0]:
882	phaseVary += cellVary(pfx,SGData)
883	resRBData = PhaseData[name]['RBModels'].get('Residue',[])
884	if resRBData:
885	rbids = rbIds['Residue'] #NB: used in the MakeRB routines
886	for iRB,RB in enumerate(resRBData):
887	MakeRBParms('R',phaseVary,phaseDict)
888	MakeRBThermals('R',phaseVary,phaseDict)
889	MakeRBTorsions('R',phaseVary,phaseDict)
890
891	vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
892	if vecRBData:
893	rbids = rbIds['Vector'] #NB: used in the MakeRB routines
894	for iRB,RB in enumerate(vecRBData):
895	MakeRBParms('V',phaseVary,phaseDict)
896	MakeRBThermals('V',phaseVary,phaseDict)
897
898	Natoms[pfx] = 0
899	if Atoms and not General.get('doPawley'):
900	cx,ct,cs,cia = General['AtomPtrs']
901	if General['Type'] in ['nuclear','macromolecular']:
902	Natoms[pfx] = len(Atoms)
903	for i,at in enumerate(Atoms):
904	atomIndx[at[-1]] = [pfx,i] #lookup table for restraints
905	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
906	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
907	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
908	pfx+'AI/A:'+str(i):at[cia],})
909	if at[cia] == 'I':
910	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
911	else:
912	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
913	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
914	if 'F' in at[ct+1]:
915	phaseVary.append(pfx+'Afrac:'+str(i))
916	if 'X' in at[ct+1]:
917	xId,xCoef = G2spc.GetCSxinel(at[cs])
918	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
919	equivs = [[],[],[]]
920	for j in range(3):
921	if xId[j] > 0:
922	phaseVary.append(names[j])
923	equivs[xId[j]-1].append([names[j],xCoef[j]])
924	for equiv in equivs:
925	if len(equiv) > 1:
926	name = equiv[0][0]
927	for eqv in equiv[1:]:
928	G2mv.StoreEquivalence(name,(eqv,))
929	if 'U' in at[ct+1]:
930	if at[9] == 'I':
931	phaseVary.append(pfx+'AUiso:'+str(i))
932	else:
933	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
934	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
935	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
936	equivs = [[],[],[],[],[],[]]
937	for j in range(6):
938	if uId[j] > 0:
939	phaseVary.append(names[j])
940	equivs[uId[j]-1].append([names[j],uCoef[j]])
941	for equiv in equivs:
942	if len(equiv) > 1:
943	name = equiv[0][0]
944	for eqv in equiv[1:]:
945	G2mv.StoreEquivalence(name,(eqv,))
946	# elif General['Type'] == 'magnetic':
947	# elif General['Type'] == 'macromolecular':
948	textureData = General['SH Texture']
949	if textureData['Order']:
950	phaseDict[pfx+'SHorder'] = textureData['Order']
951	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
952	for item in ['omega','chi','phi']:
953	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
954	if textureData['Sample '+item][0]:
955	phaseVary.append(pfx+'SH '+item)
956	for item in textureData['SH Coeff'][1]:
957	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
958	if textureData['SH Coeff'][0]:
959	phaseVary.append(pfx+item)
960
961	if Print:
962	print >>pFile,'\n Phase name: ',General['Name']
963	print >>pFile,135*'-'
964	PrintFFtable(FFtable)
965	PrintBLtable(BLtable)
966	print >>pFile,''
967	for line in SGtext: print >>pFile,line
968	PrintRBObjects(resRBData,vecRBData)
969	PrintAtoms(General,Atoms)
970	print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
971	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
972	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
973	PrintTexture(textureData)
974	if name in RestraintDict:
975	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
976	textureData,RestraintDict[name],pFile)
977
978	elif PawleyRef:
979	pawleyVary = []
980	for i,refl in enumerate(PawleyRef):
981	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
982	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
983	if refl[5]:
984	pawleyVary.append(pfx+'PWLref:'+str(i))
985	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
986	phaseVary += pawleyVary
987
988	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
989
990	def cellFill(pfx,SGData,parmDict,sigDict):
991	if SGData['SGLaue'] in ['-1',]:
992	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
993	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
994	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
995	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
996	elif SGData['SGLaue'] in ['2/m',]:
997	if SGData['SGUniq'] == 'a':
998	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
999	parmDict[pfx+'A3'],0,0]
1000	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1001	sigDict[pfx+'A3'],0,0]
1002	elif SGData['SGUniq'] == 'b':
1003	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1004	0,parmDict[pfx+'A4'],0]
1005	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1006	0,sigDict[pfx+'A4'],0]
1007	else:
1008	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1009	0,0,parmDict[pfx+'A5']]
1010	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1011	0,0,sigDict[pfx+'A5']]
1012	elif SGData['SGLaue'] in ['mmm',]:
1013	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1014	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1015	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1016	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1017	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1018	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1019	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1020	parmDict[pfx+'A0'],0,0]
1021	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1022	elif SGData['SGLaue'] in ['3R', '3mR']:
1023	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1024	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1025	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1026	elif SGData['SGLaue'] in ['m3m','m3']:
1027	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1028	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1029	return A,sigA
1030
1031	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1032	if phaseRest:
1033	Amat = G2lat.cell2AB(cell)[0]
1034	cx,ct,cs = AtPtrs[:3]
1035	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1036	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1037	['ChemComp','Sites'],['Texture','HKLs']]
1038	for name,rest in names:
1039	itemRest = phaseRest[name]
1040	if itemRest[rest] and itemRest['Use']:
1041	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1042	if name in ['Bond','Angle','Plane','Chiral']:
1043	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
1044	for indx,ops,obs,esd in itemRest[rest]:
1045	try:
1046	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1047	AtName = ''
1048	for i,Aname in enumerate(AtNames):
1049	AtName += Aname
1050	if ops[i] == '1':
1051	AtName += '-'
1052	else:
1053	AtName += '+('+ops[i]+')-'
1054	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1055	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1056	if name == 'Bond':
1057	calc = G2mth.getRestDist(XYZ,Amat)
1058	elif name == 'Angle':
1059	calc = G2mth.getRestAngle(XYZ,Amat)
1060	elif name == 'Plane':
1061	calc = G2mth.getRestPlane(XYZ,Amat)
1062	elif name == 'Chiral':
1063	calc = G2mth.getRestChiral(XYZ,Amat)
1064	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1065	except KeyError:
1066	del itemRest[rest]
1067	elif name in ['Torsion','Rama']:
1068	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
1069	coeffDict = itemRest['Coeff']
1070	for indx,ops,cofName,esd in itemRest[rest]:
1071	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1072	AtName = ''
1073	for i,Aname in enumerate(AtNames):
1074	AtName += Aname+'+('+ops[i]+')-'
1075	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1076	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1077	if name == 'Torsion':
1078	tor = G2mth.getRestTorsion(XYZ,Amat)
1079	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1080	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1081	else:
1082	phi,psi = G2mth.getRestRama(XYZ,Amat)
1083	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1084	print >>pFile,' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1])
1085	elif name == 'ChemComp':
1086	print >>pFile,' atoms mul*frac factor prod'
1087	for indx,factors,obs,esd in itemRest[rest]:
1088	try:
1089	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1090	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1091	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1092	mulfrac = mul*frac
1093	calcs = mulfracfactors
1094	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1095	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1096	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1097	except KeyError:
1098	del itemRest[rest]
1099	elif name == 'Texture' and textureData['Order']:
1100	Start = False
1101	SHCoef = textureData['SH Coeff'][1]
1102	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1103	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1104	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
1105	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1106	PH = np.array(hkl)
1107	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1108	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1109	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1110	Z = ma.masked_greater(Z,0.0)
1111	num = ma.count(Z)
1112	sum = 0
1113	if num:
1114	sum = np.sum(Z)
1115	print ' %d %d %d %d %8.3f %8.3f %8d %s %8.3f'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2)
1116
1117	def getCellEsd(pfx,SGData,A,covData):
1118	dpr = 180./np.pi
1119	rVsq = G2lat.calc_rVsq(A)
1120	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1121	cell = np.array(G2lat.Gmat2cell(g)) #real cell
1122	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1123	scos = cosd(cellst[3:6])
1124	ssin = sind(cellst[3:6])
1125	scot = scos/ssin
1126	rcos = cosd(cell[3:6])
1127	rsin = sind(cell[3:6])
1128	rcot = rcos/rsin
1129	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1130	varyList = covData['varyList']
1131	covMatrix = covData['covMatrix']
1132	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1133	Ax = np.array(A)
1134	Ax[3:] /= 2.
1135	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
1136	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
1137	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1138	Vol = 1/np.sqrt(rVsq)
1139	sigVol = Vol*3np.sqrt(srcvlsq)/2.
1140	R123 = Ax[0]Ax[1]Ax[2]
1141	dsasdg = np.zeros((3,6))
1142	dadg = np.zeros((6,6))
1143	for i0 in range(3): #0 1 2
1144	i1 = (i0+1)%3 #1 2 0
1145	i2 = (i1+1)%3 #2 0 1
1146	i3 = 5-i2 #3 5 4
1147	i4 = 5-i1 #4 3 5
1148	i5 = 5-i0 #5 4 3
1149	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
1150	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
1151	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1152	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
1153	denom = np.sqrt(denmsq)
1154	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
1155	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
1156	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
1157	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
1158	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
1159	dadg[i5][i5] = -Ax[i0]/denom
1160	for i0 in range(3):
1161	i1 = (i0+1)%3
1162	i2 = (i1+1)%3
1163	i3 = 5-i2
1164	for ij in range(6):
1165	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1166	if ij == i0:
1167	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1168	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
1169	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1170	var = np.diag(sigMat)
1171	CS = np.where(var>0.,np.sqrt(var),0.)
1172	cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
1173	return cellSig
1174
1175	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1176
1177	def PrintAtomsAndSig(General,Atoms,atomsSig):
1178	print >>pFile,'\n Atoms:'
1179	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1180	if General['Type'] == 'magnetic':
1181	line += ' Mx My Mz'
1182	elif General['Type'] == 'macromolecular':
1183	line = ' res no residue chain '+line
1184	print >>pFile,line
1185	if General['Type'] == 'nuclear':
1186	print >>pFile,135*'-'
1187	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1188	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1189	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
1190	for atyp in General['AtomTypes']: #zero composition
1191	General['NoAtoms'][atyp] = 0.
1192	for i,at in enumerate(Atoms):
1193	General['NoAtoms'][at[1]] += at[6]*float(at[8]) #new composition
1194	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1195	valstr = ' values:'
1196	sigstr = ' sig :'
1197	for ind in [3,4,5,6]:
1198	sigind = str(i)+':'+str(ind)
1199	valstr += fmt[ind]%(at[ind])
1200	if sigind in atomsSig:
1201	sigstr += fmt[ind]%(atomsSig[sigind])
1202	else:
1203	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1204	if at[9] == 'I':
1205	valstr += fmt[10]%(at[10])
1206	if str(i)+':10' in atomsSig:
1207	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1208	else:
1209	sigstr += 8*' '
1210	else:
1211	valstr += 8*' '
1212	sigstr += 8*' '
1213	for ind in [11,12,13,14,15,16]:
1214	sigind = str(i)+':'+str(ind)
1215	valstr += fmt[ind]%(at[ind])
1216	if sigind in atomsSig:
1217	sigstr += fmt[ind]%(atomsSig[sigind])
1218	else:
1219	sigstr += 8*' '
1220	print >>pFile,name
1221	print >>pFile,valstr
1222	print >>pFile,sigstr
1223
1224	def PrintRBObjPOAndSig(rbfx,rbsx):
1225	namstr = ' names :'
1226	valstr = ' values:'
1227	sigstr = ' esds :'
1228	for i,px in enumerate(['Px:','Py:','Pz:']):
1229	name = pfx+rbfx+px+rbsx
1230	namstr += '%12s'%('Pos '+px[1])
1231	valstr += '%12.5f'%(parmDict[name])
1232	if name in sigDict:
1233	sigstr += '%12.5f'%(sigDict[name])
1234	else:
1235	sigstr += 12*' '
1236	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1237	name = pfx+rbfx+po+rbsx
1238	namstr += '%12s'%('Ori '+po[1])
1239	valstr += '%12.5f'%(parmDict[name])
1240	if name in sigDict:
1241	sigstr += '%12.5f'%(sigDict[name])
1242	else:
1243	sigstr += 12*' '
1244	print >>pFile,namstr
1245	print >>pFile,valstr
1246	print >>pFile,sigstr
1247
1248	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1249	namstr = ' names :'
1250	valstr = ' values:'
1251	sigstr = ' esds :'
1252	if 'N' not in TLS:
1253	print >>pFile,' Thermal motion:'
1254	if 'T' in TLS:
1255	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1256	name = pfx+rbfx+pt+rbsx
1257	namstr += '%12s'%(pt[:3])
1258	valstr += '%12.5f'%(parmDict[name])
1259	if name in sigDict:
1260	sigstr += '%12.5f'%(sigDict[name])
1261	else:
1262	sigstr += 12*' '
1263	print >>pFile,namstr
1264	print >>pFile,valstr
1265	print >>pFile,sigstr
1266	if 'L' in TLS:
1267	namstr = ' names :'
1268	valstr = ' values:'
1269	sigstr = ' esds :'
1270	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1271	name = pfx+rbfx+pt+rbsx
1272	namstr += '%12s'%(pt[:3])
1273	valstr += '%12.3f'%(parmDict[name])
1274	if name in sigDict:
1275	sigstr += '%12.3f'%(sigDict[name])
1276	else:
1277	sigstr += 12*' '
1278	print >>pFile,namstr
1279	print >>pFile,valstr
1280	print >>pFile,sigstr
1281	if 'S' in TLS:
1282	namstr = ' names :'
1283	valstr = ' values:'
1284	sigstr = ' esds :'
1285	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1286	name = pfx+rbfx+pt+rbsx
1287	namstr += '%12s'%(pt[:3])
1288	valstr += '%12.4f'%(parmDict[name])
1289	if name in sigDict:
1290	sigstr += '%12.4f'%(sigDict[name])
1291	else:
1292	sigstr += 12*' '
1293	print >>pFile,namstr
1294	print >>pFile,valstr
1295	print >>pFile,sigstr
1296	if 'U' in TLS:
1297	name = pfx+rbfx+'U:'+rbsx
1298	namstr = ' names :'+'%12s'%('Uiso')
1299	valstr = ' values:'+'%12.5f'%(parmDict[name])
1300	if name in sigDict:
1301	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1302	else:
1303	sigstr = ' esds :'+12*' '
1304	print >>pFile,namstr
1305	print >>pFile,valstr
1306	print >>pFile,sigstr
1307
1308	def PrintRBObjTorAndSig(rbsx):
1309	namstr = ' names :'
1310	valstr = ' values:'
1311	sigstr = ' esds :'
1312	nTors = len(RBObj['Torsions'])
1313	if nTors:
1314	print >>pFile,' Torsions:'
1315	for it in range(nTors):
1316	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1317	namstr += '%12s'%('Tor'+str(it))
1318	valstr += '%12.4f'%(parmDict[name])
1319	if name in sigDict:
1320	sigstr += '%12.4f'%(sigDict[name])
1321	print >>pFile,namstr
1322	print >>pFile,valstr
1323	print >>pFile,sigstr
1324
1325	def PrintSHtextureAndSig(textureData,SHtextureSig):
1326	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1327	names = ['omega','chi','phi']
1328	namstr = ' names :'
1329	ptstr = ' values:'
1330	sigstr = ' esds :'
1331	for name in names:
1332	namstr += '%12s'%(name)
1333	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1334	if 'Sample '+name in SHtextureSig:
1335	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1336	else:
1337	sigstr += 12*' '
1338	print >>pFile,namstr
1339	print >>pFile,ptstr
1340	print >>pFile,sigstr
1341	print >>pFile,'\n Texture coefficients:'
1342	namstr = ' names :'
1343	ptstr = ' values:'
1344	sigstr = ' esds :'
1345	SHcoeff = textureData['SH Coeff'][1]
1346	for name in SHcoeff:
1347	namstr += '%12s'%(name)
1348	ptstr += '%12.3f'%(SHcoeff[name])
1349	if name in SHtextureSig:
1350	sigstr += '%12.3f'%(SHtextureSig[name])
1351	else:
1352	sigstr += 12*' '
1353	print >>pFile,namstr
1354	print >>pFile,ptstr
1355	print >>pFile,sigstr
1356
1357	print >>pFile,'\n Phases:'
1358	for phase in Phases:
1359	print >>pFile,' Result for phase: ',phase
1360	Phase = Phases[phase]
1361	General = Phase['General']
1362	SGData = General['SGData']
1363	Atoms = Phase['Atoms']
1364	AtLookup = G2mth.FillAtomLookUp(Atoms)
1365	cell = General['Cell']
1366	pId = Phase['pId']
1367	pfx = str(pId)+'::'
1368	if cell[0]:
1369	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1370	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1371	print >>pFile,' Reciprocal metric tensor: '
1372	ptfmt = "%15.9f"
1373	names = ['A11','A22','A33','A12','A13','A23']
1374	namstr = ' names :'
1375	ptstr = ' values:'
1376	sigstr = ' esds :'
1377	for name,a,siga in zip(names,A,sigA):
1378	namstr += '%15s'%(name)
1379	ptstr += ptfmt%(a)
1380	if siga:
1381	sigstr += ptfmt%(siga)
1382	else:
1383	sigstr += 15*' '
1384	print >>pFile,namstr
1385	print >>pFile,ptstr
1386	print >>pFile,sigstr
1387	cell[1:7] = G2lat.A2cell(A)
1388	cell[7] = G2lat.calc_V(A)
1389	print >>pFile,' New unit cell:'
1390	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1391	names = ['a','b','c','alpha','beta','gamma','Volume']
1392	namstr = ' names :'
1393	ptstr = ' values:'
1394	sigstr = ' esds :'
1395	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1396	namstr += '%12s'%(name)
1397	ptstr += fmt%(a)
1398	if siga:
1399	sigstr += fmt%(siga)
1400	else:
1401	sigstr += 12*' '
1402	print >>pFile,namstr
1403	print >>pFile,ptstr
1404	print >>pFile,sigstr
1405
1406	General['Mass'] = 0.
1407	if Phase['General'].get('doPawley'):
1408	pawleyRef = Phase['Pawley ref']
1409	for i,refl in enumerate(pawleyRef):
1410	key = pfx+'PWLref:'+str(i)
1411	refl[6] = parmDict[key]
1412	if key in sigDict:
1413	refl[7] = sigDict[key]
1414	else:
1415	refl[7] = 0
1416	else:
1417	VRBIds = RBIds['Vector']
1418	RRBIds = RBIds['Residue']
1419	RBModels = Phase['RBModels']
1420	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1421	jrb = VRBIds.index(RBObj['RBId'])
1422	rbsx = str(irb)+':'+str(jrb)
1423	print >>pFile,' Vector rigid body parameters:'
1424	PrintRBObjPOAndSig('RBV',rbsx)
1425	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1426	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1427	jrb = RRBIds.index(RBObj['RBId'])
1428	rbsx = str(irb)+':'+str(jrb)
1429	print >>pFile,' Residue rigid body parameters:'
1430	PrintRBObjPOAndSig('RBR',rbsx)
1431	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1432	PrintRBObjTorAndSig(rbsx)
1433	atomsSig = {}
1434	if General['Type'] == 'nuclear': #this needs macromolecular variant!
1435	for i,at in enumerate(Atoms):
1436	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1437	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1438	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1439	for ind in [3,4,5,6]:
1440	at[ind] = parmDict[names[ind]]
1441	if ind in [3,4,5]:
1442	name = names[ind].replace('A','dA')
1443	else:
1444	name = names[ind]
1445	if name in sigDict:
1446	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1447	if at[9] == 'I':
1448	at[10] = parmDict[names[10]]
1449	if names[10] in sigDict:
1450	atomsSig[str(i)+':10'] = sigDict[names[10]]
1451	else:
1452	for ind in [11,12,13,14,15,16]:
1453	at[ind] = parmDict[names[ind]]
1454	if names[ind] in sigDict:
1455	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1456	ind = General['AtomTypes'].index(at[1])
1457	General['Mass'] += General['AtomMass'][ind]at[6]at[8]
1458	PrintAtomsAndSig(General,Atoms,atomsSig)
1459
1460	textureData = General['SH Texture']
1461	if textureData['Order']:
1462	SHtextureSig = {}
1463	for name in ['omega','chi','phi']:
1464	aname = pfx+'SH '+name
1465	textureData['Sample '+name][1] = parmDict[aname]
1466	if aname in sigDict:
1467	SHtextureSig['Sample '+name] = sigDict[aname]
1468	for name in textureData['SH Coeff'][1]:
1469	aname = pfx+name
1470	textureData['SH Coeff'][1][name] = parmDict[aname]
1471	if aname in sigDict:
1472	SHtextureSig[name] = sigDict[aname]
1473	PrintSHtextureAndSig(textureData,SHtextureSig)
1474	if phase in RestraintDict:
1475	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1476	textureData,RestraintDict[phase],pFile)
1477
1478	################################################################################
1479	##### Histogram & Phase data
1480	################################################################################
1481
1482	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None):
1483
1484	def PrintSize(hapData):
1485	if hapData[0] in ['isotropic','uniaxial']:
1486	line = '\n Size model : %9s'%(hapData[0])
1487	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1488	if hapData[0] == 'uniaxial':
1489	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1490	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1491	print >>pFile,line
1492	else:
1493	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1494	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1495	Snames = ['S11','S22','S33','S12','S13','S23']
1496	ptlbls = ' names :'
1497	ptstr = ' values:'
1498	varstr = ' refine:'
1499	for i,name in enumerate(Snames):
1500	ptlbls += '%12s' % (name)
1501	ptstr += '%12.6f' % (hapData[4][i])
1502	varstr += '%12s' % (str(hapData[5][i]))
1503	print >>pFile,ptlbls
1504	print >>pFile,ptstr
1505	print >>pFile,varstr
1506
1507	def PrintMuStrain(hapData,SGData):
1508	if hapData[0] in ['isotropic','uniaxial']:
1509	line = '\n Mustrain model: %9s'%(hapData[0])
1510	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1511	if hapData[0] == 'uniaxial':
1512	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1513	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1514	print >>pFile,line
1515	else:
1516	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1517	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1518	Snames = G2spc.MustrainNames(SGData)
1519	ptlbls = ' names :'
1520	ptstr = ' values:'
1521	varstr = ' refine:'
1522	for i,name in enumerate(Snames):
1523	ptlbls += '%12s' % (name)
1524	ptstr += '%12.6f' % (hapData[4][i])
1525	varstr += '%12s' % (str(hapData[5][i]))
1526	print >>pFile,ptlbls
1527	print >>pFile,ptstr
1528	print >>pFile,varstr
1529
1530	def PrintHStrain(hapData,SGData):
1531	print >>pFile,'\n Hydrostatic/elastic strain: '
1532	Hsnames = G2spc.HStrainNames(SGData)
1533	ptlbls = ' names :'
1534	ptstr = ' values:'
1535	varstr = ' refine:'
1536	for i,name in enumerate(Hsnames):
1537	ptlbls += '%12s' % (name)
1538	ptstr += '%12.6f' % (hapData[0][i])
1539	varstr += '%12s' % (str(hapData[1][i]))
1540	print >>pFile,ptlbls
1541	print >>pFile,ptstr
1542	print >>pFile,varstr
1543
1544	def PrintSHPO(hapData):
1545	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1546	str(hapData[4])+' Refine? '+str(hapData[2])
1547	ptlbls = ' names :'
1548	ptstr = ' values:'
1549	for item in hapData[5]:
1550	ptlbls += '%12s'%(item)
1551	ptstr += '%12.3f'%(hapData[5][item])
1552	print >>pFile,ptlbls
1553	print >>pFile,ptstr
1554
1555	def PrintBabinet(hapData):
1556	print >>pFile,'\n Babinet form factor modification: '
1557	ptlbls = ' names :'
1558	ptstr = ' values:'
1559	varstr = ' refine:'
1560	for name in ['BabA','BabU']:
1561	ptlbls += '%12s' % (name)
1562	ptstr += '%12.6f' % (hapData[name][0])
1563	varstr += '%12s' % (str(hapData[name][1]))
1564	print >>pFile,ptlbls
1565	print >>pFile,ptstr
1566	print >>pFile,varstr
1567
1568
1569	hapDict = {}
1570	hapVary = []
1571	controlDict = {}
1572	poType = {}
1573	poAxes = {}
1574	spAxes = {}
1575	spType = {}
1576
1577	for phase in Phases:
1578	HistoPhase = Phases[phase]['Histograms']
1579	SGData = Phases[phase]['General']['SGData']
1580	cell = Phases[phase]['General']['Cell'][1:7]
1581	A = G2lat.cell2A(cell)
1582	pId = Phases[phase]['pId']
1583	histoList = HistoPhase.keys()
1584	histoList.sort()
1585	for histogram in histoList:
1586	try:
1587	Histogram = Histograms[histogram]
1588	except KeyError:
1589	#skip if histogram not included e.g. in a sequential refinement
1590	continue
1591	hapData = HistoPhase[histogram]
1592	hId = Histogram['hId']
1593	if 'PWDR' in histogram:
1594	limits = Histogram['Limits'][1]
1595	inst = Histogram['Instrument Parameters'][0]
1596	Zero = inst['Zero'][1]
1597	if 'C' in inst['Type'][1]:
1598	try:
1599	wave = inst['Lam'][1]
1600	except KeyError:
1601	wave = inst['Lam1'][1]
1602	dmin = wave/(2.0*sind(limits[1]/2.0))
1603	pfx = str(pId)+':'+str(hId)+':'
1604	for item in ['Scale','Extinction']:
1605	hapDict[pfx+item] = hapData[item][0]
1606	if hapData[item][1]:
1607	hapVary.append(pfx+item)
1608	names = G2spc.HStrainNames(SGData)
1609	for i,name in enumerate(names):
1610	hapDict[pfx+name] = hapData['HStrain'][0][i]
1611	if hapData['HStrain'][1][i]:
1612	hapVary.append(pfx+name)
1613	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1614	if hapData['Pref.Ori.'][0] == 'MD':
1615	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1616	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1617	if hapData['Pref.Ori.'][2]:
1618	hapVary.append(pfx+'MD')
1619	else: #'SH' spherical harmonics
1620	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1621	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1622	for item in hapData['Pref.Ori.'][5]:
1623	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1624	if hapData['Pref.Ori.'][2]:
1625	hapVary.append(pfx+item)
1626	for item in ['Mustrain','Size']:
1627	controlDict[pfx+item+'Type'] = hapData[item][0]
1628	hapDict[pfx+item+';mx'] = hapData[item][1][2]
1629	if hapData[item][2][2]:
1630	hapVary.append(pfx+item+';mx')
1631	if hapData[item][0] in ['isotropic','uniaxial']:
1632	hapDict[pfx+item+';i'] = hapData[item][1][0]
1633	if hapData[item][2][0]:
1634	hapVary.append(pfx+item+';i')
1635	if hapData[item][0] == 'uniaxial':
1636	controlDict[pfx+item+'Axis'] = hapData[item][3]
1637	hapDict[pfx+item+';a'] = hapData[item][1][1]
1638	if hapData[item][2][1]:
1639	hapVary.append(pfx+item+';a')
1640	else: #generalized for mustrain or ellipsoidal for size
1641	Nterms = len(hapData[item][4])
1642	if item == 'Mustrain':
1643	names = G2spc.MustrainNames(SGData)
1644	pwrs = []
1645	for name in names:
1646	h,k,l = name[1:]
1647	pwrs.append([int(h),int(k),int(l)])
1648	controlDict[pfx+'MuPwrs'] = pwrs
1649	for i in range(Nterms):
1650	sfx = ':'+str(i)
1651	hapDict[pfx+item+sfx] = hapData[item][4][i]
1652	if hapData[item][5][i]:
1653	hapVary.append(pfx+item+sfx)
1654	for bab in ['BabA','BabU']:
1655	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1656	if hapData['Babinet'][bab][1]:
1657	hapVary.append(pfx+bab)
1658
1659	if Print:
1660	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1661	print >>pFile,135*'-'
1662	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1663	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1664	if hapData['Pref.Ori.'][0] == 'MD':
1665	Ax = hapData['Pref.Ori.'][3]
1666	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1667	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1668	else: #'SH' for spherical harmonics
1669	PrintSHPO(hapData['Pref.Ori.'])
1670	PrintSize(hapData['Size'])
1671	PrintMuStrain(hapData['Mustrain'],SGData)
1672	PrintHStrain(hapData['HStrain'],SGData)
1673	if hapData['Babinet']['BabA'][0]:
1674	PrintBabinet(hapData['Babinet'])
1675	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1676	refList = []
1677	for h,k,l,d in HKLd:
1678	ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1679	mul *= 2 # for powder overlap of Friedel pairs
1680	if ext:
1681	continue
1682	if 'C' in inst['Type'][0]:
1683	pos = 2.0asind(wave/(2.0d))+Zero
1684	if limits[0] < pos < limits[1]:
1685	refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1686	#last item should contain form factor values by atom type
1687	else:
1688	raise ValueError
1689	Histogram['Reflection Lists'][phase] = refList
1690	elif 'HKLF' in histogram:
1691	inst = Histogram['Instrument Parameters'][0]
1692	hId = Histogram['hId']
1693	hfx = ':%d:'%(hId)
1694	for item in inst:
1695	hapDict[hfx+item] = inst[item][1]
1696	pfx = str(pId)+':'+str(hId)+':'
1697	hapDict[pfx+'Scale'] = hapData['Scale'][0]
1698	if hapData['Scale'][1]:
1699	hapVary.append(pfx+'Scale')
1700
1701	extApprox,extType,extParms = hapData['Extinction']
1702	controlDict[pfx+'EType'] = extType
1703	controlDict[pfx+'EApprox'] = extApprox
1704	controlDict[pfx+'Tbar'] = extParms['Tbar']
1705	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1706	if 'Primary' in extType:
1707	Ekey = ['Ep',]
1708	elif 'I & II' in extType:
1709	Ekey = ['Eg','Es']
1710	elif 'Secondary Type II' == extType:
1711	Ekey = ['Es',]
1712	elif 'Secondary Type I' == extType:
1713	Ekey = ['Eg',]
1714	else: #'None'
1715	Ekey = []
1716	for eKey in Ekey:
1717	hapDict[pfx+eKey] = extParms[eKey][0]
1718	if extParms[eKey][1]:
1719	hapVary.append(pfx+eKey)
1720	for bab in ['BabA','BabU']:
1721	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1722	if hapData['Babinet'][bab][1]:
1723	hapVary.append(pfx+bab)
1724	if Print:
1725	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1726	print >>pFile,135*'-'
1727	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1728	if extType != 'None':
1729	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1730	text = ' Parameters :'
1731	for eKey in Ekey:
1732	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1733	print >>pFile,text
1734	if hapData['Babinet']['BabA'][0]:
1735	PrintBabinet(hapData['Babinet'])
1736	Histogram['Reflection Lists'] = phase
1737
1738	return hapVary,hapDict,controlDict
1739
1740	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1741
1742	def PrintSizeAndSig(hapData,sizeSig):
1743	line = '\n Size model: %9s'%(hapData[0])
1744	refine = False
1745	if hapData[0] in ['isotropic','uniaxial']:
1746	line += ' equatorial:%12.4f'%(hapData[1][0])
1747	if sizeSig[0][0]:
1748	line += ', sig:%8.4f'%(sizeSig[0][0])
1749	refine = True
1750	if hapData[0] == 'uniaxial':
1751	line += ' axial:%12.4f'%(hapData[1][1])
1752	if sizeSig[0][1]:
1753	refine = True
1754	line += ', sig:%8.4f'%(sizeSig[0][1])
1755	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1756	if sizeSig[0][2]:
1757	refine = True
1758	line += ', sig:%8.4f'%(sizeSig[0][2])
1759	if refine:
1760	print >>pFile,line
1761	else:
1762	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1763	if sizeSig[0][2]:
1764	refine = True
1765	line += ', sig:%8.4f'%(sizeSig[0][2])
1766	Snames = ['S11','S22','S33','S12','S13','S23']
1767	ptlbls = ' name :'
1768	ptstr = ' value :'
1769	sigstr = ' sig :'
1770	for i,name in enumerate(Snames):
1771	ptlbls += '%12s' % (name)
1772	ptstr += '%12.6f' % (hapData[4][i])
1773	if sizeSig[1][i]:
1774	refine = True
1775	sigstr += '%12.6f' % (sizeSig[1][i])
1776	else:
1777	sigstr += 12*' '
1778	if refine:
1779	print >>pFile,line
1780	print >>pFile,ptlbls
1781	print >>pFile,ptstr
1782	print >>pFile,sigstr
1783
1784	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1785	line = '\n Mustrain model: %9s'%(hapData[0])
1786	refine = False
1787	if hapData[0] in ['isotropic','uniaxial']:
1788	line += ' equatorial:%12.1f'%(hapData[1][0])
1789	if mustrainSig[0][0]:
1790	line += ', sig:%8.1f'%(mustrainSig[0][0])
1791	refine = True
1792	if hapData[0] == 'uniaxial':
1793	line += ' axial:%12.1f'%(hapData[1][1])
1794	if mustrainSig[0][1]:
1795	line += ', sig:%8.1f'%(mustrainSig[0][1])
1796	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1797	if mustrainSig[0][2]:
1798	refine = True
1799	line += ', sig:%8.3f'%(mustrainSig[0][2])
1800	if refine:
1801	print >>pFile,line
1802	else:
1803	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1804	if mustrainSig[0][2]:
1805	refine = True
1806	line += ', sig:%8.3f'%(mustrainSig[0][2])
1807	Snames = G2spc.MustrainNames(SGData)
1808	ptlbls = ' name :'
1809	ptstr = ' value :'
1810	sigstr = ' sig :'
1811	for i,name in enumerate(Snames):
1812	ptlbls += '%12s' % (name)
1813	ptstr += '%12.6f' % (hapData[4][i])
1814	if mustrainSig[1][i]:
1815	refine = True
1816	sigstr += '%12.6f' % (mustrainSig[1][i])
1817	else:
1818	sigstr += 12*' '
1819	if refine:
1820	print >>pFile,line
1821	print >>pFile,ptlbls
1822	print >>pFile,ptstr
1823	print >>pFile,sigstr
1824
1825	def PrintHStrainAndSig(hapData,strainSig,SGData):
1826	Hsnames = G2spc.HStrainNames(SGData)
1827	ptlbls = ' name :'
1828	ptstr = ' value :'
1829	sigstr = ' sig :'
1830	refine = False
1831	for i,name in enumerate(Hsnames):
1832	ptlbls += '%12s' % (name)
1833	ptstr += '%12.6g' % (hapData[0][i])
1834	if name in strainSig:
1835	refine = True
1836	sigstr += '%12.6g' % (strainSig[name])
1837	else:
1838	sigstr += 12*' '
1839	if refine:
1840	print >>pFile,'\n Hydrostatic/elastic strain: '
1841	print >>pFile,ptlbls
1842	print >>pFile,ptstr
1843	print >>pFile,sigstr
1844
1845	def PrintSHPOAndSig(pfx,hapData,POsig):
1846	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1847	ptlbls = ' names :'
1848	ptstr = ' values:'
1849	sigstr = ' sig :'
1850	for item in hapData[5]:
1851	ptlbls += '%12s'%(item)
1852	ptstr += '%12.3f'%(hapData[5][item])
1853	if pfx+item in POsig:
1854	sigstr += '%12.3f'%(POsig[pfx+item])
1855	else:
1856	sigstr += 12*' '
1857	print >>pFile,ptlbls
1858	print >>pFile,ptstr
1859	print >>pFile,sigstr
1860
1861	def PrintExtAndSig(pfx,hapData,ScalExtSig):
1862	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
1863	text = ''
1864	for item in hapData[2]:
1865	if pfx+item in ScalExtSig:
1866	text += ' %s: '%(item)
1867	text += '%12.2e'%(hapData[2][item][0])
1868	if pfx+item in ScalExtSig:
1869	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
1870	print >>pFile,text
1871
1872	def PrintBabinetAndSig(pfx,hapData,BabSig):
1873	print >>pFile,'\n Babinet form factor modification: '
1874	ptlbls = ' names :'
1875	ptstr = ' values:'
1876	sigstr = ' sig :'
1877	for item in hapData:
1878	ptlbls += '%12s'%(item)
1879	ptstr += '%12.3f'%(hapData[item][0])
1880	if pfx+item in BabSig:
1881	sigstr += '%12.3f'%(BabSig[pfx+item])
1882	else:
1883	sigstr += 12*' '
1884	print >>pFile,ptlbls
1885	print >>pFile,ptstr
1886	print >>pFile,sigstr
1887
1888	PhFrExtPOSig = {}
1889	SizeMuStrSig = {}
1890	ScalExtSig = {}
1891	BabSig = {}
1892	wtFrSum = {}
1893	for phase in Phases:
1894	HistoPhase = Phases[phase]['Histograms']
1895	General = Phases[phase]['General']
1896	SGData = General['SGData']
1897	pId = Phases[phase]['pId']
1898	histoList = HistoPhase.keys()
1899	histoList.sort()
1900	for histogram in histoList:
1901	try:
1902	Histogram = Histograms[histogram]
1903	except KeyError:
1904	#skip if histogram not included e.g. in a sequential refinement
1905	continue
1906	hapData = HistoPhase[histogram]
1907	hId = Histogram['hId']
1908	pfx = str(pId)+':'+str(hId)+':'
1909	if hId not in wtFrSum:
1910	wtFrSum[hId] = 0.
1911	if 'PWDR' in histogram:
1912	for item in ['Scale','Extinction']:
1913	hapData[item][0] = parmDict[pfx+item]
1914	if pfx+item in sigDict:
1915	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1916	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
1917	if hapData['Pref.Ori.'][0] == 'MD':
1918	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1919	if pfx+'MD' in sigDict:
1920	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
1921	else: #'SH' spherical harmonics
1922	for item in hapData['Pref.Ori.'][5]:
1923	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1924	if pfx+item in sigDict:
1925	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1926	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1927	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
1928	pfx+'HStrain':{}})
1929	for item in ['Mustrain','Size']:
1930	hapData[item][1][2] = parmDict[pfx+item+';mx']
1931	hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
1932	if pfx+item+';mx' in sigDict:
1933	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
1934	if hapData[item][0] in ['isotropic','uniaxial']:
1935	hapData[item][1][0] = parmDict[pfx+item+';i']
1936	if item == 'Size':
1937	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1938	if pfx+item+';i' in sigDict:
1939	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
1940	if hapData[item][0] == 'uniaxial':
1941	hapData[item][1][1] = parmDict[pfx+item+';a']
1942	if item == 'Size':
1943	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
1944	if pfx+item+';a' in sigDict:
1945	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
1946	else: #generalized for mustrain or ellipsoidal for size
1947	Nterms = len(hapData[item][4])
1948	for i in range(Nterms):
1949	sfx = ':'+str(i)
1950	hapData[item][4][i] = parmDict[pfx+item+sfx]
1951	if pfx+item+sfx in sigDict:
1952	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
1953	names = G2spc.HStrainNames(SGData)
1954	for i,name in enumerate(names):
1955	hapData['HStrain'][0][i] = parmDict[pfx+name]
1956	if pfx+name in sigDict:
1957	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
1958	for name in ['BabA','BabU']:
1959	hapData['Babinet'][name][0] = parmDict[pfx+name]
1960	if pfx+name in sigDict:
1961	BabSig[pfx+name] = sigDict[pfx+name]
1962
1963	elif 'HKLF' in histogram:
1964	for item in ['Scale',]:
1965	if parmDict.get(pfx+item):
1966	hapData[item][0] = parmDict[pfx+item]
1967	if pfx+item in sigDict:
1968	ScalExtSig[pfx+item] = sigDict[pfx+item]
1969	for item in ['Ep','Eg','Es']:
1970	if parmDict.get(pfx+item):
1971	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
1972	if pfx+item in sigDict:
1973	ScalExtSig[pfx+item] = sigDict[pfx+item]
1974	for name in ['BabA','BabU']:
1975	hapData['Babinet'][name][0] = parmDict[pfx+name]
1976	if pfx+name in sigDict:
1977	BabSig[pfx+name] = sigDict[pfx+name]
1978
1979	if Print:
1980	for phase in Phases:
1981	HistoPhase = Phases[phase]['Histograms']
1982	General = Phases[phase]['General']
1983	SGData = General['SGData']
1984	pId = Phases[phase]['pId']
1985	histoList = HistoPhase.keys()
1986	histoList.sort()
1987	for histogram in histoList:
1988	try:
1989	Histogram = Histograms[histogram]
1990	except KeyError:
1991	#skip if histogram not included e.g. in a sequential refinement
1992	continue
1993	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1994	print >>pFile,130*'-'
1995	hapData = HistoPhase[histogram]
1996	hId = Histogram['hId']
1997	pfx = str(pId)+':'+str(hId)+':'
1998	if 'PWDR' in histogram:
1999	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2000	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2001
2002	if pfx+'Scale' in PhFrExtPOSig:
2003	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2004	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2005	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2006	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2007	if pfx+'Extinction' in PhFrExtPOSig:
2008	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2009	if hapData['Pref.Ori.'][0] == 'MD':
2010	if pfx+'MD' in PhFrExtPOSig:
2011	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2012	else:
2013	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2014	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2015	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2016	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2017	if len(BabSig):
2018	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2019
2020	elif 'HKLF' in histogram:
2021	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2022	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2023	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2024	if pfx+'Scale' in ScalExtSig:
2025	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2026	if hapData['Extinction'][0] != 'None':
2027	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2028	if len(BabSig):
2029	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2030
2031	################################################################################
2032	##### Histogram data
2033	################################################################################
2034
2035	def GetHistogramData(Histograms,Print=True,pFile=None):
2036
2037	def GetBackgroundParms(hId,Background):
2038	Back = Background[0]
2039	DebyePeaks = Background[1]
2040	bakType,bakFlag = Back[:2]
2041	backVals = Back[3:]
2042	backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2043	backDict = dict(zip(backNames,backVals))
2044	backVary = []
2045	if bakFlag:
2046	backVary = backNames
2047	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2048	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2049	debyeDict = {}
2050	debyeList = []
2051	for i in range(DebyePeaks['nDebye']):
2052	debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2053	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2054	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2055	debyeVary = []
2056	for item in debyeList:
2057	if item[1]:
2058	debyeVary.append(item[0])
2059	backDict.update(debyeDict)
2060	backVary += debyeVary
2061	peakDict = {}
2062	peakList = []
2063	for i in range(DebyePeaks['nPeaks']):
2064	peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
2065	':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
2066	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2067	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2068	peakVary = []
2069	for item in peakList:
2070	if item[1]:
2071	peakVary.append(item[0])
2072	backDict.update(peakDict)
2073	backVary += peakVary
2074	return bakType,backDict,backVary
2075
2076	def GetInstParms(hId,Inst):
2077	dataType = Inst['Type'][0]
2078	instDict = {}
2079	insVary = []
2080	pfx = ':'+str(hId)+':'
2081	insKeys = Inst.keys()
2082	insKeys.sort()
2083	for item in insKeys:
2084	insName = pfx+item
2085	instDict[insName] = Inst[item][1]
2086	if Inst[item][2]:
2087	insVary.append(insName)
2088	# instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2089	# instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2090	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2091	return dataType,instDict,insVary
2092
2093	def GetSampleParms(hId,Sample):
2094	sampVary = []
2095	hfx = ':'+str(hId)+':'
2096	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2097	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2098	Type = Sample['Type']
2099	if 'Bragg' in Type: #Bragg-Brentano
2100	for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering?
2101	sampDict[hfx+item] = Sample[item][0]
2102	if Sample[item][1]:
2103	sampVary.append(hfx+item)
2104	elif 'Debye' in Type: #Debye-Scherrer
2105	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2106	sampDict[hfx+item] = Sample[item][0]
2107	if Sample[item][1]:
2108	sampVary.append(hfx+item)
2109	return Type,sampDict,sampVary
2110
2111	def PrintBackground(Background):
2112	Back = Background[0]
2113	DebyePeaks = Background[1]
2114	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2115	line = ' Coefficients: '
2116	for i,back in enumerate(Back[3:]):
2117	line += '%10.3f'%(back)
2118	if i and not i%10:
2119	line += '\n'+15*' '
2120	print >>pFile,line
2121	if DebyePeaks['nDebye']:
2122	print >>pFile,'\n Debye diffuse scattering coefficients'
2123	parms = ['DebyeA','DebyeR','DebyeU']
2124	line = ' names : '
2125	for parm in parms:
2126	line += '%8s refine?'%(parm)
2127	print >>pFile,line
2128	for j,term in enumerate(DebyePeaks['debyeTerms']):
2129	line = ' term'+'%2d'%(j)+':'
2130	for i in range(3):
2131	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2132	print >>pFile,line
2133	if DebyePeaks['nPeaks']:
2134	print >>pFile,'\n Single peak coefficients'
2135	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2136	line = ' names : '
2137	for parm in parms:
2138	line += '%8s refine?'%(parm)
2139	print >>pFile,line
2140	for j,term in enumerate(DebyePeaks['peaksList']):
2141	line = ' peak'+'%2d'%(j)+':'
2142	for i in range(4):
2143	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2144	print >>pFile,line
2145
2146	def PrintInstParms(Inst):
2147	print >>pFile,'\n Instrument Parameters:'
2148	ptlbls = ' name :'
2149	ptstr = ' value :'
2150	varstr = ' refine:'
2151	insKeys = Inst.keys()
2152	insKeys.sort()
2153	for item in insKeys:
2154	if item != 'Type':
2155	ptlbls += '%12s' % (item)
2156	ptstr += '%12.6f' % (Inst[item][1])
2157	if item in ['Lam1','Lam2','Azimuth']:
2158	varstr += 12*' '
2159	else:
2160	varstr += '%12s' % (str(bool(Inst[item][2])))
2161	print >>pFile,ptlbls
2162	print >>pFile,ptstr
2163	print >>pFile,varstr
2164
2165	def PrintSampleParms(Sample):
2166	print >>pFile,'\n Sample Parameters:'
2167	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2168	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2169	ptlbls = ' name :'
2170	ptstr = ' value :'
2171	varstr = ' refine:'
2172	if 'Bragg' in Sample['Type']:
2173	for item in ['Scale','Shift','Transparency']:
2174	ptlbls += '%14s'%(item)
2175	ptstr += '%14.4f'%(Sample[item][0])
2176	varstr += '%14s'%(str(bool(Sample[item][1])))
2177
2178	elif 'Debye' in Type: #Debye-Scherrer
2179	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2180	ptlbls += '%14s'%(item)
2181	ptstr += '%14.4f'%(Sample[item][0])
2182	varstr += '%14s'%(str(bool(Sample[item][1])))
2183
2184	print >>pFile,ptlbls
2185	print >>pFile,ptstr
2186	print >>pFile,varstr
2187
2188	histDict = {}
2189	histVary = []
2190	controlDict = {}
2191	histoList = Histograms.keys()
2192	histoList.sort()
2193	for histogram in histoList:
2194	if 'PWDR' in histogram:
2195	Histogram = Histograms[histogram]
2196	hId = Histogram['hId']
2197	pfx = ':'+str(hId)+':'
2198	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2199	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2200
2201	Background = Histogram['Background']
2202	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2203	controlDict[pfx+'bakType'] = Type
2204	histDict.update(bakDict)
2205	histVary += bakVary
2206
2207	Inst = Histogram['Instrument Parameters'][0]
2208	Type,instDict,insVary = GetInstParms(hId,Inst)
2209	controlDict[pfx+'histType'] = Type
2210	if pfx+'Lam1' in instDict:
2211	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2212	else:
2213	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2214	histDict.update(instDict)
2215	histVary += insVary
2216
2217	Sample = Histogram['Sample Parameters']
2218	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2219	controlDict[pfx+'instType'] = Type
2220	histDict.update(sampDict)
2221	histVary += sampVary
2222
2223
2224	if Print:
2225	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2226	print >>pFile,135*'-'
2227	Units = {'C':' deg','T':' msec'}
2228	units = Units[controlDict[pfx+'histType'][2]]
2229	Limits = controlDict[pfx+'Limits']
2230	print >>pFile,' Instrument type: ',Sample['Type']
2231	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2232	PrintSampleParms(Sample)
2233	PrintInstParms(Inst)
2234	PrintBackground(Background)
2235	elif 'HKLF' in histogram:
2236	Histogram = Histograms[histogram]
2237	hId = Histogram['hId']
2238	pfx = ':'+str(hId)+':'
2239	controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
2240	Inst = Histogram['Instrument Parameters'][0]
2241	controlDict[pfx+'histType'] = Inst['Type'][0]
2242	histDict[pfx+'Lam'] = Inst['Lam'][1]
2243	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
2244	return histVary,histDict,controlDict
2245
2246	def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2247
2248	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2249	Back = Background[0]
2250	DebyePeaks = Background[1]
2251	lenBack = len(Back[3:])
2252	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
2253	for i in range(lenBack):
2254	Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2255	if pfx+'Back:'+str(i) in sigDict:
2256	backSig[i] = sigDict[pfx+'Back:'+str(i)]
2257	if DebyePeaks['nDebye']:
2258	for i in range(DebyePeaks['nDebye']):
2259	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2260	for j,name in enumerate(names):
2261	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2262	if name in sigDict:
2263	backSig[lenBack+3*i+j] = sigDict[name]
2264	if DebyePeaks['nPeaks']:
2265	for i in range(DebyePeaks['nPeaks']):
2266	names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
2267	pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
2268	for j,name in enumerate(names):
2269	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2270	if name in sigDict:
2271	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
2272	return backSig
2273
2274	def SetInstParms(pfx,Inst,parmDict,sigDict):
2275	instSig = {}
2276	insKeys = Inst.keys()
2277	insKeys.sort()
2278	for item in insKeys:
2279	insName = pfx+item
2280	Inst[item][1] = parmDict[insName]
2281	if insName in sigDict:
2282	instSig[item] = sigDict[insName]
2283	else:
2284	instSig[item] = 0
2285	return instSig
2286
2287	def SetSampleParms(pfx,Sample,parmDict,sigDict):
2288	if 'Bragg' in Sample['Type']: #Bragg-Brentano
2289	sampSig = [0 for i in range(3)]
2290	for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering?
2291	Sample[item][0] = parmDict[pfx+item]
2292	if pfx+item in sigDict:
2293	sampSig[i] = sigDict[pfx+item]
2294	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2295	sampSig = [0 for i in range(4)]
2296	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2297	Sample[item][0] = parmDict[pfx+item]
2298	if pfx+item in sigDict:
2299	sampSig[i] = sigDict[pfx+item]
2300	return sampSig
2301
2302	def PrintBackgroundSig(Background,backSig):
2303	Back = Background[0]
2304	DebyePeaks = Background[1]
2305	lenBack = len(Back[3:])
2306	valstr = ' value : '
2307	sigstr = ' sig : '
2308	refine = False
2309	for i,back in enumerate(Back[3:]):
2310	valstr += '%10.4g'%(back)
2311	if Back[1]:
2312	refine = True
2313	sigstr += '%10.4g'%(backSig[i])
2314	else:
2315	sigstr += 10*' '
2316	if refine:
2317	print >>pFile,'\n Background function: ',Back[0]
2318	print >>pFile,valstr
2319	print >>pFile,sigstr
2320	if DebyePeaks['nDebye']:
2321	ifAny = False
2322	ptfmt = "%12.3f"
2323	names = ' names :'
2324	ptstr = ' values:'
2325	sigstr = ' esds :'
2326	for item in sigDict:
2327	if 'Debye' in item:
2328	ifAny = True
2329	names += '%12s'%(item)
2330	ptstr += ptfmt%(parmDict[item])
2331	sigstr += ptfmt%(sigDict[item])
2332	if ifAny:
2333	print >>pFile,'\n Debye diffuse scattering coefficients'
2334	print >>pFile,names
2335	print >>pFile,ptstr
2336	print >>pFile,sigstr
2337	if DebyePeaks['nPeaks']:
2338	ifAny = False
2339	ptfmt = "%14.3f"
2340	names = ' names :'
2341	ptstr = ' values:'
2342	sigstr = ' esds :'
2343	for item in sigDict:
2344	if 'BkPk' in item:
2345	ifAny = True
2346	names += '%14s'%(item)
2347	ptstr += ptfmt%(parmDict[item])
2348	sigstr += ptfmt%(sigDict[item])
2349	if ifAny:
2350	print >>pFile,'\n Single peak coefficients'
2351	print >>pFile,names
2352	print >>pFile,ptstr
2353	print >>pFile,sigstr
2354
2355	def PrintInstParmsSig(Inst,instSig):
2356	ptlbls = ' names :'
2357	ptstr = ' value :'
2358	sigstr = ' sig :'
2359	refine = False
2360	insKeys = instSig.keys()
2361	insKeys.sort()
2362	for name in insKeys:
2363	if name not in ['Type','Lam1','Lam2','Azimuth']:
2364	ptlbls += '%12s' % (name)
2365	ptstr += '%12.6f' % (Inst[name][1])
2366	if instSig[name]:
2367	refine = True
2368	sigstr += '%12.6f' % (instSig[name])
2369	else:
2370	sigstr += 12*' '
2371	if refine:
2372	print >>pFile,'\n Instrument Parameters:'
2373	print >>pFile,ptlbls
2374	print >>pFile,ptstr
2375	print >>pFile,sigstr
2376
2377	def PrintSampleParmsSig(Sample,sampleSig):
2378	ptlbls = ' names :'
2379	ptstr = ' values:'
2380	sigstr = ' sig :'
2381	refine = False
2382	if 'Bragg' in Sample['Type']:
2383	for i,item in enumerate(['Scale','Shift','Transparency']):
2384	ptlbls += '%14s'%(item)
2385	ptstr += '%14.4f'%(Sample[item][0])
2386	if sampleSig[i]:
2387	refine = True
2388	sigstr += '%14.4f'%(sampleSig[i])
2389	else:
2390	sigstr += 14*' '
2391
2392	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2393	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2394	ptlbls += '%14s'%(item)
2395	ptstr += '%14.4f'%(Sample[item][0])
2396	if sampleSig[i]:
2397	refine = True
2398	sigstr += '%14.4f'%(sampleSig[i])
2399	else:
2400	sigstr += 14*' '
2401
2402	if refine:
2403	print >>pFile,'\n Sample Parameters:'
2404	print >>pFile,ptlbls
2405	print >>pFile,ptstr
2406	print >>pFile,sigstr
2407
2408	histoList = Histograms.keys()
2409	histoList.sort()
2410	for histogram in histoList:
2411	if 'PWDR' in histogram:
2412	Histogram = Histograms[histogram]
2413	hId = Histogram['hId']
2414	pfx = ':'+str(hId)+':'
2415	Background = Histogram['Background']
2416	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2417
2418	Inst = Histogram['Instrument Parameters'][0]
2419	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2420
2421	Sample = Histogram['Sample Parameters']
2422	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2423
2424	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2425	print >>pFile,135*'-'
2426	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
2427	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2428	if Print:
2429	print >>pFile,' Instrument type: ',Sample['Type']
2430	PrintSampleParmsSig(Sample,sampSig)
2431	PrintInstParmsSig(Inst,instSig)
2432	PrintBackgroundSig(Background,backSig)
2433

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

Download in other formats: