Context Navigation

source: trunk/GSASIIstrIO.py @ 2534

Last change on this file since 2534 was 2520, checked in by vondreele, 9 years ago
remove unused & commented Srr fctr & deriv routines Mag moment derivatives semireasonable but not close transform nucl - mag ok if no off diag terms in transform matrix
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 145.3 KB

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

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

Download in other formats: