Context Navigation

source: trunk/GSASIIstrIO.py @ 1664

Last change on this file since 1664 was 1664, checked in by vondreele, 8 years ago
fix rhombohedral lattice constraints
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 131.0 KB

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

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

Download in other formats: