Context Navigation

source: trunk/GSASIIstrIO.py @ 1021

Last change on this file since 1021 was 1021, checked in by toby, 12 years ago
more digits in size
File size: 110.7 KB

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

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

Download in other formats: