Context Navigation

source: trunk/GSASIIstrIO.py @ 1045

Last change on this file since 1045 was 1045, checked in by vondreele, 10 years ago
fix Fourier map problems - didn't like Fo<0. use ErrorDialog? for some errors.
File size: 111.1 KB

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