Context Navigation

source: trunk/GSASIIstrIO.py @ 2481

Last change on this file since 2481 was 2481, checked in by vondreele, 7 years ago
force Transform to delete nonmagnetic atoms if phase made magnetic & add 'mag' to new phase name fix TOF cosine background function magnetic structure refinement works (in numeric mode only) magnetic structure factor calculation correct
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 145.1 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2016-09-30 15:33:12 +0000 (Fri, 30 Sep 2016) $
4	# $Author: vondreele $
5	# $Revision: 2481 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 2481 2016-09-30 15:33:12Z 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: 2481 $")
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 = G2spc.SytSym(at[cx:cx+3],SGData)
1212	if SpnFlp[Nop] > 0.: #black use p
1213	mId,mCoef = G2spc.GetCSpqinel(SytSym)[0]
1214	else: #red use q
1215	mId,mCoef = G2spc.GetCSpqinel(SytSym)[1]
1216	names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)]
1217	equivs = [[],[],[]]
1218	for j in range(3):
1219	if mId[j] > 0:
1220	phaseVary.append(names[j])
1221	equivs[mId[j]-1].append([names[j],mCoef[j]])
1222	for equiv in equivs:
1223	if len(equiv) > 1:
1224	name = equiv[0][0]
1225	coef = equiv[0][1]
1226	for eqv in equiv[1:]:
1227	eqv[1] /= coef
1228	G2mv.StoreEquivalence(name,(eqv,))
1229	if General.get('Modulated',False):
1230	AtomSS = at[-1]['SS1']
1231	waveType = AtomSS['waveType']
1232	phaseDict[pfx+'waveType:'+str(i)] = waveType
1233	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1234	Waves = AtomSS[Stype]
1235	nx = 0
1236	for iw,wave in enumerate(Waves):
1237	if not iw:
1238	if waveType in ['ZigZag','Block']:
1239	nx = 1
1240	CSI = G2spc.GetSSfxuinel(waveType,1,at[cx:cx+3],SGData,SSGData)
1241	else:
1242	CSI = G2spc.GetSSfxuinel('Fourier',iw+1-nx,at[cx:cx+3],SGData,SSGData)
1243	uId,uCoef = CSI[Stype]
1244	stiw = str(i)+':'+str(iw)
1245	if Stype == 'Spos':
1246	if waveType in ['ZigZag','Block',] and not iw:
1247	names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1248	equivs = [[],[], [],[],[]]
1249	else:
1250	names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1251	pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1252	equivs = [[],[],[], [],[],[]]
1253	elif Stype == 'Sadp':
1254	names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1255	pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1256	pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1257	pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1258	equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1259	elif Stype == 'Sfrac':
1260	equivs = [[],[]]
1261	if 'Crenel' in waveType and not iw:
1262	names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1263	else:
1264	names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1265	elif Stype == 'Smag':
1266	equivs = [[],[],[], [],[],[]]
1267	names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1268	pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1269	phaseDict.update(dict(zip(names,wave[0])))
1270	if wave[1]: #what do we do here for multiple terms in modulation constraints?
1271	for j in range(len(equivs)):
1272	if uId[j][0] > 0:
1273	phaseVary.append(names[j])
1274	equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1275	for equiv in equivs:
1276	if len(equiv) > 1:
1277	name = equiv[0][0]
1278	coef = equiv[0][1]
1279	for eqv in equiv[1:]:
1280	eqv[1] /= coef
1281	G2mv.StoreEquivalence(name,equiv[1:])
1282	maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1283	textureData = General['SH Texture']
1284	if textureData['Order'] and not seqRef:
1285	phaseDict[pfx+'SHorder'] = textureData['Order']
1286	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1287	for item in ['omega','chi','phi']:
1288	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1289	if textureData['Sample '+item][0]:
1290	phaseVary.append(pfx+'SH '+item)
1291	for item in textureData['SH Coeff'][1]:
1292	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1293	if textureData['SH Coeff'][0]:
1294	phaseVary.append(pfx+item)
1295
1296	if Print:
1297	print >>pFile,'\n Phase name: ',General['Name']
1298	print >>pFile,135*'-'
1299	PrintFFtable(FFtable)
1300	PrintBLtable(BLtable)
1301	if General['Type'] == 'magnetic':
1302	PrintMFtable(MFtable)
1303	print >>pFile,''
1304	#how do we print magnetic symmetry table? TBD
1305	if len(SSGtext): #if superstructure
1306	for line in SSGtext: print >>pFile,line
1307	if len(SSGtable):
1308	for item in SSGtable:
1309	line = ' %s '%(item)
1310	print >>pFile,line
1311	else:
1312	print >>pFile,' ( 1) %s'%(SSGtable[0])
1313	else:
1314	for line in SGtext: print >>pFile,line
1315	if len(SGtable):
1316	for item in SGtable:
1317	line = ' %s '%(item)
1318	print >>pFile,line
1319	else:
1320	print >>pFile,' ( 1) %s'%(SGtable[0])
1321	PrintRBObjects(resRBData,vecRBData)
1322	PrintAtoms(General,Atoms)
1323	if General['Type'] == 'magnetic':
1324	PrintMoments(General,Atoms)
1325	if General.get('Modulated',False):
1326	PrintWaves(General,Atoms)
1327	print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1328	' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1329	' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1330	if len(SSGtext): #if superstructure
1331	print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]), \
1332	' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1333	if not seqRef:
1334	PrintTexture(textureData)
1335	if name in RestraintDict:
1336	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1337	textureData,RestraintDict[name],pFile)
1338
1339	elif PawleyRef:
1340	if Print:
1341	print >>pFile,'\n Phase name: ',General['Name']
1342	print >>pFile,135*'-'
1343	print >>pFile,''
1344	if len(SSGtext): #if superstructure
1345	for line in SSGtext: print >>pFile,line
1346	if len(SSGtable):
1347	for item in SSGtable:
1348	line = ' %s '%(item)
1349	print >>pFile,line
1350	else:
1351	print >>pFile,' ( 1) %s'%(SSGtable[0])
1352	else:
1353	for line in SGtext: print >>pFile,line
1354	if len(SGtable):
1355	for item in SGtable:
1356	line = ' %s '%(item)
1357	print >>pFile,line
1358	else:
1359	print >>pFile,' ( 1) %s'%(SGtable[0])
1360	print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1361	' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1362	' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1363	if len(SSGtext): #if superstructure
1364	print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]), \
1365	' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1366	pawleyVary = []
1367	for i,refl in enumerate(PawleyRef):
1368	phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1369	if im:
1370	pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1371	else:
1372	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1373	if refl[5+im]:
1374	pawleyVary.append(pfx+'PWLref:'+str(i))
1375	GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary) #does G2mv.StoreEquivalence
1376	phaseVary += pawleyVary
1377
1378	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
1379
1380	def cellFill(pfx,SGData,parmDict,sigDict):
1381	'''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1382	from the parameter and sig dictionaries.
1383
1384	:param str pfx: parameter prefix ("n::", where n is a phase number)
1385	:param dict SGdata: a symmetry object
1386	:param dict parmDict: a dictionary of parameters
1387	:param dict sigDict: a dictionary of uncertainties on parameters
1388
1389	:returns: A,sigA where each is a list of six terms with the A terms
1390	'''
1391	if SGData['SGLaue'] in ['-1',]:
1392	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1393	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1394	elif SGData['SGLaue'] in ['2/m',]:
1395	if SGData['SGUniq'] == 'a':
1396	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1397	0,0,parmDict[pfx+'A5']]
1398	elif SGData['SGUniq'] == 'b':
1399	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1400	0,parmDict[pfx+'A4'],0]
1401	else:
1402	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1403	parmDict[pfx+'A3'],0,0]
1404	elif SGData['SGLaue'] in ['mmm',]:
1405	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1406	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1407	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1408	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1409	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1410	parmDict[pfx+'A0'],0,0]
1411	elif SGData['SGLaue'] in ['3R', '3mR']:
1412	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1413	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1414	elif SGData['SGLaue'] in ['m3m','m3']:
1415	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1416
1417	try:
1418	if SGData['SGLaue'] in ['-1',]:
1419	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1420	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1421	elif SGData['SGLaue'] in ['2/m',]:
1422	if SGData['SGUniq'] == 'a':
1423	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1424	0,0,sigDict[pfx+'A5']]
1425	elif SGData['SGUniq'] == 'b':
1426	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1427	0,sigDict[pfx+'A4'],0]
1428	else:
1429	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1430	sigDict[pfx+'A3'],0,0]
1431	elif SGData['SGLaue'] in ['mmm',]:
1432	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1433	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1434	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1435	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1436	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1437	elif SGData['SGLaue'] in ['3R', '3mR']:
1438	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1439	elif SGData['SGLaue'] in ['m3m','m3']:
1440	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1441	except KeyError:
1442	sigA = [0,0,0,0,0,0]
1443	return A,sigA
1444
1445	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1446	'needs a doc string'
1447	if phaseRest:
1448	Amat = G2lat.cell2AB(cell)[0]
1449	cx,ct,cs = AtPtrs[:3]
1450	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1451	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1452	['ChemComp','Sites'],['Texture','HKLs']]
1453	for name,rest in names:
1454	itemRest = phaseRest[name]
1455	if itemRest[rest] and itemRest['Use']:
1456	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1457	if name in ['Bond','Angle','Plane','Chiral']:
1458	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
1459	for indx,ops,obs,esd in itemRest[rest]:
1460	try:
1461	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1462	AtName = ''
1463	for i,Aname in enumerate(AtNames):
1464	AtName += Aname
1465	if ops[i] == '1':
1466	AtName += '-'
1467	else:
1468	AtName += '+('+ops[i]+')-'
1469	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1470	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1471	if name == 'Bond':
1472	calc = G2mth.getRestDist(XYZ,Amat)
1473	elif name == 'Angle':
1474	calc = G2mth.getRestAngle(XYZ,Amat)
1475	elif name == 'Plane':
1476	calc = G2mth.getRestPlane(XYZ,Amat)
1477	elif name == 'Chiral':
1478	calc = G2mth.getRestChiral(XYZ,Amat)
1479	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1480	except KeyError:
1481	del itemRest[rest]
1482	elif name in ['Torsion','Rama']:
1483	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
1484	coeffDict = itemRest['Coeff']
1485	for indx,ops,cofName,esd in itemRest[rest]:
1486	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1487	AtName = ''
1488	for i,Aname in enumerate(AtNames):
1489	AtName += Aname+'+('+ops[i]+')-'
1490	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1491	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1492	if name == 'Torsion':
1493	tor = G2mth.getRestTorsion(XYZ,Amat)
1494	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1495	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1496	else:
1497	phi,psi = G2mth.getRestRama(XYZ,Amat)
1498	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1499	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])
1500	elif name == 'ChemComp':
1501	print >>pFile,' atoms mul*frac factor prod'
1502	for indx,factors,obs,esd in itemRest[rest]:
1503	try:
1504	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1505	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1506	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1507	mulfrac = mul*frac
1508	calcs = mulfracfactors
1509	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1510	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1511	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1512	except KeyError:
1513	del itemRest[rest]
1514	elif name == 'Texture' and textureData['Order']:
1515	Start = False
1516	SHCoef = textureData['SH Coeff'][1]
1517	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1518	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1519	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
1520	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1521	PH = np.array(hkl)
1522	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1523	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1524	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1525	Z = ma.masked_greater(Z,0.0)
1526	num = ma.count(Z)
1527	sum = 0
1528	if num:
1529	sum = np.sum(Z)
1530	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)
1531
1532	def getCellEsd(pfx,SGData,A,covData):
1533	'needs a doc string'
1534	dpr = 180./np.pi
1535	rVsq = G2lat.calc_rVsq(A)
1536	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1537	cell = np.array(G2lat.Gmat2cell(g)) #real cell
1538	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1539	scos = cosd(cellst[3:6])
1540	ssin = sind(cellst[3:6])
1541	scot = scos/ssin
1542	rcos = cosd(cell[3:6])
1543	rsin = sind(cell[3:6])
1544	rcot = rcos/rsin
1545	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1546	varyList = covData['varyList']
1547	covMatrix = covData['covMatrix']
1548	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1549	Ax = np.array(A)
1550	Ax[3:] /= 2.
1551	drVdA = np.array([
1552	Ax[1]Ax[2]-Ax[5]*2,
1553	Ax[0]Ax[2]-Ax[4]*2,
1554	Ax[0]Ax[1]-Ax[3]*2,
1555	Ax[4]Ax[5]-Ax[2]Ax[3],
1556	Ax[3]Ax[5]-Ax[1]Ax[4],
1557	Ax[3]Ax[4]-Ax[0]Ax[5]])
1558	srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1559	Vol = 1/np.sqrt(rVsq)
1560	sigVol = Vol*3np.sqrt(srcvlsq)/2. #ok - checks with GSAS
1561	R123 = Ax[0]Ax[1]Ax[2]
1562	dsasdg = np.zeros((3,6))
1563	dadg = np.zeros((6,6))
1564	for i0 in range(3): #0 1 2
1565	i1 = (i0+1)%3 #1 2 0
1566	i2 = (i1+1)%3 #2 0 1
1567	i3 = 5-i2 #3 5 4
1568	i4 = 5-i1 #4 3 5
1569	i5 = 5-i0 #5 4 3
1570	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
1571	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
1572	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1573	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
1574	denom = np.sqrt(denmsq)
1575	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
1576	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
1577	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
1578	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
1579	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
1580	dadg[i5][i5] = -Ax[i0]/denom
1581	for i0 in range(3):
1582	i1 = (i0+1)%3
1583	i2 = (i1+1)%3
1584	i3 = 5-i2
1585	for ij in range(6):
1586	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1587	if ij == i0:
1588	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1589	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
1590	sigMat = np.inner(dadg,np.inner(dadg,vcov))
1591	var = np.diag(sigMat)
1592	CS = np.where(var>0.,np.sqrt(var),0.)
1593	return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1594
1595	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1596	'needs a doc string'
1597
1598	def PrintAtomsAndSig(General,Atoms,atomsSig):
1599	print >>pFile,'\n Atoms:'
1600	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1601	if General['Type'] == 'macromolecular':
1602	line = ' res no residue chain '+line
1603	cx,ct,cs,cia = General['AtomPtrs']
1604	print >>pFile,line
1605	print >>pFile,135*'-'
1606	fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1607	cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1608	noFXsig = {cx:[10' ','%10s'],cx+1:[10' ','%10s'],cx+2:[10' ','%10s'],cx+3:[8' ','%8s']}
1609	for atyp in General['AtomTypes']: #zero composition
1610	General['NoAtoms'][atyp] = 0.
1611	for i,at in enumerate(Atoms):
1612	General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5]) #new composition
1613	if General['Type'] == 'macromolecular':
1614	name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1615	valstr = ' values: '
1616	sigstr = ' sig : '
1617	else:
1618	name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1619	valstr = ' values:'
1620	sigstr = ' sig :'
1621	for ind in range(cx,cx+4):
1622	sigind = str(i)+':'+str(ind)
1623	valstr += fmt[ind]%(at[ind])
1624	if sigind in atomsSig:
1625	sigstr += fmt[ind]%(atomsSig[sigind])
1626	else:
1627	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1628	if at[cia] == 'I':
1629	valstr += fmt[cia+1]%(at[cia+1])
1630	if '%d:%d'%(i,cia+1) in atomsSig:
1631	sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1632	else:
1633	sigstr += 8*' '
1634	else:
1635	valstr += 8*' '
1636	sigstr += 8*' '
1637	for ind in range(cia+2,cia+8):
1638	sigind = str(i)+':'+str(ind)
1639	valstr += fmt[ind]%(at[ind])
1640	if sigind in atomsSig:
1641	sigstr += fmt[ind]%(atomsSig[sigind])
1642	else:
1643	sigstr += 8*' '
1644	print >>pFile,name
1645	print >>pFile,valstr
1646	print >>pFile,sigstr
1647
1648	def PrintMomentsAndSig(General,Atoms,atomsSig):
1649	print >>pFile,'\n Magnetic Moments:' #add magnitude & angle, etc.? TBD
1650	line = ' name Mx My Mz'
1651	cx,ct,cs,cia = General['AtomPtrs']
1652	cmx = cx+4
1653	AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1654	print >>pFile,line
1655	print >>pFile,135*'-'
1656	fmt = {0:'%7s',ct:'%7s',cmx:'%10.5f',cmx+1:'%10.5f',cmx+2:'%10.5f'}
1657	noFXsig = {cmx:[10' ','%10s'],cmx+1:[10' ','%10s'],cmx+2:[10*' ','%10s']}
1658	for i,at in enumerate(Atoms):
1659	if AtInfo[at[ct]]:
1660	name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1661	valstr = ' values:'
1662	sigstr = ' sig :'
1663	for ind in range(cmx,cmx+3):
1664	sigind = str(i)+':'+str(ind)
1665	valstr += fmt[ind]%(at[ind])
1666	if sigind in atomsSig:
1667	sigstr += fmt[ind]%(atomsSig[sigind])
1668	else:
1669	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1670	print >>pFile,name
1671	print >>pFile,valstr
1672	print >>pFile,sigstr
1673
1674	def PrintWavesAndSig(General,Atoms,wavesSig):
1675	cx,ct,cs,cia = General['AtomPtrs']
1676	print >>pFile,'\n Modulation waves'
1677	names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1678	'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1679	'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1680	print >>pFile,135*'-'
1681	for i,at in enumerate(Atoms):
1682	AtomSS = at[-1]['SS1']
1683	waveType = AtomSS['waveType']
1684	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1685	Waves = AtomSS[Stype]
1686	if len(Waves):
1687	print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:' \
1688	%(at[ct-1],at[cs],Stype,waveType)
1689	line = ''
1690	for iw,wave in enumerate(Waves):
1691	stiw = ':'+str(i)+':'+str(iw)
1692	namstr = ' names :'
1693	valstr = ' values:'
1694	sigstr = ' esds :'
1695	if Stype == 'Spos':
1696	nt = 6
1697	ot = 0
1698	if waveType in ['ZigZag','Block',] and not iw:
1699	nt = 5
1700	ot = 6
1701	for j in range(nt):
1702	name = names['Spos'][j+ot]
1703	namstr += '%12s'%(name)
1704	valstr += '%12.4f'%(wave[0][j])
1705	if name+stiw in wavesSig:
1706	sigstr += '%12.4f'%(wavesSig[name+stiw])
1707	else:
1708	sigstr += 12*' '
1709	elif Stype == 'Sfrac':
1710	ot = 0
1711	if 'Crenel' in waveType and not iw:
1712	ot = 2
1713	for j in range(2):
1714	name = names['Sfrac'][j+ot]
1715	namstr += '%12s'%(names['Sfrac'][j+ot])
1716	valstr += '%12.4f'%(wave[0][j])
1717	if name+stiw in wavesSig:
1718	sigstr += '%12.4f'%(wavesSig[name+stiw])
1719	else:
1720	sigstr += 12*' '
1721	elif Stype == 'Sadp':
1722	for j in range(12):
1723	name = names['Sadp'][j]
1724	namstr += '%10s'%(names['Sadp'][j])
1725	valstr += '%10.6f'%(wave[0][j])
1726	if name+stiw in wavesSig:
1727	sigstr += '%10.6f'%(wavesSig[name+stiw])
1728	else:
1729	sigstr += 10*' '
1730	elif Stype == 'Smag':
1731	for j in range(6):
1732	name = names['Smag'][j]
1733	namstr += '%12s'%(names['Smag'][j])
1734	valstr += '%12.4f'%(wave[0][j])
1735	if name+stiw in wavesSig:
1736	sigstr += '%12.4f'%(wavesSig[name+stiw])
1737	else:
1738	sigstr += 12*' '
1739
1740	print >>pFile,namstr
1741	print >>pFile,valstr
1742	print >>pFile,sigstr
1743
1744
1745	def PrintRBObjPOAndSig(rbfx,rbsx):
1746	namstr = ' names :'
1747	valstr = ' values:'
1748	sigstr = ' esds :'
1749	for i,px in enumerate(['Px:','Py:','Pz:']):
1750	name = pfx+rbfx+px+rbsx
1751	namstr += '%12s'%('Pos '+px[1])
1752	valstr += '%12.5f'%(parmDict[name])
1753	if name in sigDict:
1754	sigstr += '%12.5f'%(sigDict[name])
1755	else:
1756	sigstr += 12*' '
1757	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1758	name = pfx+rbfx+po+rbsx
1759	namstr += '%12s'%('Ori '+po[1])
1760	valstr += '%12.5f'%(parmDict[name])
1761	if name in sigDict:
1762	sigstr += '%12.5f'%(sigDict[name])
1763	else:
1764	sigstr += 12*' '
1765	print >>pFile,namstr
1766	print >>pFile,valstr
1767	print >>pFile,sigstr
1768
1769	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1770	namstr = ' names :'
1771	valstr = ' values:'
1772	sigstr = ' esds :'
1773	if 'N' not in TLS:
1774	print >>pFile,' Thermal motion:'
1775	if 'T' in TLS:
1776	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1777	name = pfx+rbfx+pt+rbsx
1778	namstr += '%12s'%(pt[:3])
1779	valstr += '%12.5f'%(parmDict[name])
1780	if name in sigDict:
1781	sigstr += '%12.5f'%(sigDict[name])
1782	else:
1783	sigstr += 12*' '
1784	print >>pFile,namstr
1785	print >>pFile,valstr
1786	print >>pFile,sigstr
1787	if 'L' in TLS:
1788	namstr = ' names :'
1789	valstr = ' values:'
1790	sigstr = ' esds :'
1791	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1792	name = pfx+rbfx+pt+rbsx
1793	namstr += '%12s'%(pt[:3])
1794	valstr += '%12.3f'%(parmDict[name])
1795	if name in sigDict:
1796	sigstr += '%12.3f'%(sigDict[name])
1797	else:
1798	sigstr += 12*' '
1799	print >>pFile,namstr
1800	print >>pFile,valstr
1801	print >>pFile,sigstr
1802	if 'S' in TLS:
1803	namstr = ' names :'
1804	valstr = ' values:'
1805	sigstr = ' esds :'
1806	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1807	name = pfx+rbfx+pt+rbsx
1808	namstr += '%12s'%(pt[:3])
1809	valstr += '%12.4f'%(parmDict[name])
1810	if name in sigDict:
1811	sigstr += '%12.4f'%(sigDict[name])
1812	else:
1813	sigstr += 12*' '
1814	print >>pFile,namstr
1815	print >>pFile,valstr
1816	print >>pFile,sigstr
1817	if 'U' in TLS:
1818	name = pfx+rbfx+'U:'+rbsx
1819	namstr = ' names :'+'%12s'%('Uiso')
1820	valstr = ' values:'+'%12.5f'%(parmDict[name])
1821	if name in sigDict:
1822	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1823	else:
1824	sigstr = ' esds :'+12*' '
1825	print >>pFile,namstr
1826	print >>pFile,valstr
1827	print >>pFile,sigstr
1828
1829	def PrintRBObjTorAndSig(rbsx):
1830	namstr = ' names :'
1831	valstr = ' values:'
1832	sigstr = ' esds :'
1833	nTors = len(RBObj['Torsions'])
1834	if nTors:
1835	print >>pFile,' Torsions:'
1836	for it in range(nTors):
1837	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1838	namstr += '%12s'%('Tor'+str(it))
1839	valstr += '%12.4f'%(parmDict[name])
1840	if name in sigDict:
1841	sigstr += '%12.4f'%(sigDict[name])
1842	print >>pFile,namstr
1843	print >>pFile,valstr
1844	print >>pFile,sigstr
1845
1846	def PrintSHtextureAndSig(textureData,SHtextureSig):
1847	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1848	names = ['omega','chi','phi']
1849	namstr = ' names :'
1850	ptstr = ' values:'
1851	sigstr = ' esds :'
1852	for name in names:
1853	namstr += '%12s'%(name)
1854	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1855	if 'Sample '+name in SHtextureSig:
1856	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1857	else:
1858	sigstr += 12*' '
1859	print >>pFile,namstr
1860	print >>pFile,ptstr
1861	print >>pFile,sigstr
1862	print >>pFile,'\n Texture coefficients:'
1863	SHcoeff = textureData['SH Coeff'][1]
1864	SHkeys = SHcoeff.keys()
1865	nCoeff = len(SHcoeff)
1866	nBlock = nCoeff/10+1
1867	iBeg = 0
1868	iFin = min(iBeg+10,nCoeff)
1869	for block in range(nBlock):
1870	namstr = ' names :'
1871	ptstr = ' values:'
1872	sigstr = ' esds :'
1873	for name in SHkeys[iBeg:iFin]:
1874	namstr += '%12s'%(name)
1875	ptstr += '%12.3f'%(SHcoeff[name])
1876	if name in SHtextureSig:
1877	sigstr += '%12.3f'%(SHtextureSig[name])
1878	else:
1879	sigstr += 12*' '
1880	print >>pFile,namstr
1881	print >>pFile,ptstr
1882	print >>pFile,sigstr
1883	iBeg += 10
1884	iFin = min(iBeg+10,nCoeff)
1885
1886	print >>pFile,'\n Phases:'
1887	for phase in Phases:
1888	print >>pFile,' Result for phase: ',phase
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	for bab in ['BabA','BabU']:
2280	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2281	if hapData['Babinet'][bab][1]:
2282	hapVary.append(pfx+bab)
2283
2284	if Print:
2285	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2286	print >>pFile,135*'-'
2287	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2288	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
2289	if hapData['Pref.Ori.'][0] == 'MD':
2290	Ax = hapData['Pref.Ori.'][3]
2291	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
2292	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
2293	else: #'SH' for spherical harmonics
2294	PrintSHPO(hapData['Pref.Ori.'])
2295	print >>pFile,' Penalty hkl list: '+str(controlDict[pfx+'SHhkl'])+' tolerance: %.2f'%(controlDict[pfx+'SHtoler'])
2296	PrintSize(hapData['Size'])
2297	PrintMuStrain(hapData['Mustrain'],SGData)
2298	PrintHStrain(hapData['HStrain'],SGData)
2299	if hapData['Babinet']['BabA'][0]:
2300	PrintBabinet(hapData['Babinet'])
2301	if resetRefList:
2302	refList = []
2303	Uniq = []
2304	Phi = []
2305	useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2306	if Phases[phase]['General'].get('Modulated',False):
2307	ifSuper = True
2308	HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2309	HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2310	for h,k,l,m,d in HKLd:
2311	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2312	mul *= 2 # for powder overlap of Friedel pairs
2313	if m or not ext or useExt:
2314	if 'C' in inst['Type'][0]:
2315	pos = G2lat.Dsp2pos(inst,d)
2316	if limits[0] < pos < limits[1]:
2317	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])
2318	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2319	Uniq.append(uniq)
2320	Phi.append(phi)
2321	elif 'T' in inst['Type'][0]:
2322	pos = G2lat.Dsp2pos(inst,d)
2323	if limits[0] < pos < limits[1]:
2324	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2325	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])
2326	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2327	Uniq.append(uniq)
2328	Phi.append(phi)
2329	else:
2330	ifSuper = False
2331	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2332	HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2333	for h,k,l,d in HKLd:
2334	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2335	mul *= 2 # for powder overlap of Friedel pairs
2336	if ext and not useExt:
2337	continue
2338	if 'C' in inst['Type'][0]:
2339	pos = G2lat.Dsp2pos(inst,d)
2340	if limits[0] < pos < limits[1]:
2341	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])
2342	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2343	Uniq.append(uniq)
2344	Phi.append(phi)
2345	elif 'T' in inst['Type'][0]:
2346	pos = G2lat.Dsp2pos(inst,d)
2347	if limits[0] < pos < limits[1]:
2348	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2349	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])
2350	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2351	Uniq.append(uniq)
2352	Phi.append(phi)
2353	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2354	elif 'HKLF' in histogram:
2355	inst = Histogram['Instrument Parameters'][0]
2356	hId = Histogram['hId']
2357	hfx = ':%d:'%(hId)
2358	for item in inst:
2359	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2360	hapDict[hfx+item] = inst[item][1]
2361	pfx = str(pId)+':'+str(hId)+':'
2362	hapDict[pfx+'Scale'] = hapData['Scale'][0]
2363	if hapData['Scale'][1]:
2364	hapVary.append(pfx+'Scale')
2365
2366	extApprox,extType,extParms = hapData['Extinction']
2367	controlDict[pfx+'EType'] = extType
2368	controlDict[pfx+'EApprox'] = extApprox
2369	if 'C' in inst['Type'][0]:
2370	controlDict[pfx+'Tbar'] = extParms['Tbar']
2371	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2372	if 'Primary' in extType:
2373	Ekey = ['Ep',]
2374	elif 'I & II' in extType:
2375	Ekey = ['Eg','Es']
2376	elif 'Secondary Type II' == extType:
2377	Ekey = ['Es',]
2378	elif 'Secondary Type I' == extType:
2379	Ekey = ['Eg',]
2380	else: #'None'
2381	Ekey = []
2382	for eKey in Ekey:
2383	hapDict[pfx+eKey] = extParms[eKey][0]
2384	if extParms[eKey][1]:
2385	hapVary.append(pfx+eKey)
2386	for bab in ['BabA','BabU']:
2387	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2388	if hapData['Babinet'][bab][1]:
2389	hapVary.append(pfx+bab)
2390	Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2391	if len(Twins) == 1:
2392	hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2393	if hapData.get('Flack',[0,False])[1]:
2394	hapVary.append(pfx+'Flack')
2395	sumTwFr = 0.
2396	controlDict[pfx+'TwinLaw'] = []
2397	controlDict[pfx+'TwinInv'] = []
2398	NTL = 0
2399	for it,twin in enumerate(Twins):
2400	if 'bool' in str(type(twin[0])):
2401	controlDict[pfx+'TwinInv'].append(twin[0])
2402	controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2403	else:
2404	NTL += 1
2405	controlDict[pfx+'TwinInv'].append(False)
2406	controlDict[pfx+'TwinLaw'].append(twin[0])
2407	if it:
2408	hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2409	sumTwFr += twin[1]
2410	else:
2411	hapDict[pfx+'TwinFr:0'] = twin[1][0]
2412	controlDict[pfx+'TwinNMN'] = twin[1][1]
2413	if Twins[0][1][1]:
2414	hapVary.append(pfx+'TwinFr:'+str(it))
2415	controlDict[pfx+'NTL'] = NTL
2416	#need to make constraint on TwinFr
2417	controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2418	if len(Twins) > 1: #force sum to unity
2419	hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2420	if Print:
2421	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2422	print >>pFile,135*'-'
2423	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2424	if extType != 'None':
2425	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox)
2426	text = ' Parameters :'
2427	for eKey in Ekey:
2428	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2429	print >>pFile,text
2430	if hapData['Babinet']['BabA'][0]:
2431	PrintBabinet(hapData['Babinet'])
2432	if not SGData['SGInv'] and len(Twins) == 1:
2433	print >>pFile,' Flack parameter: %10.3f'%(hapData['Flack'][0]),' Refine?',hapData['Flack'][1]
2434	if len(Twins) > 1:
2435	for it,twin in enumerate(Twins):
2436	if 'bool' in str(type(twin[0])):
2437	print >>pFile,' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? '%(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1]
2438	else:
2439	print >>pFile,' Twin law: %s'%(str(twin[0]).replace('\n',',')),' Twin fr.: %5.3f Refine? '%(hapDict[pfx+'TwinFr:'+str(it)]),Twins[0][1][1]
2440
2441	Histogram['Reflection Lists'] = phase
2442
2443	return hapVary,hapDict,controlDict
2444
2445	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2446	'needs a doc string'
2447
2448	def PrintSizeAndSig(hapData,sizeSig):
2449	line = '\n Size model: %9s'%(hapData[0])
2450	refine = False
2451	if hapData[0] in ['isotropic','uniaxial']:
2452	line += ' equatorial:%12.5f'%(hapData[1][0])
2453	if sizeSig[0][0]:
2454	line += ', sig:%8.4f'%(sizeSig[0][0])
2455	refine = True
2456	if hapData[0] == 'uniaxial':
2457	line += ' axial:%12.4f'%(hapData[1][1])
2458	if sizeSig[0][1]:
2459	refine = True
2460	line += ', sig:%8.4f'%(sizeSig[0][1])
2461	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2462	if sizeSig[0][2]:
2463	refine = True
2464	line += ', sig:%8.4f'%(sizeSig[0][2])
2465	if refine:
2466	print >>pFile,line
2467	else:
2468	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2469	if sizeSig[0][2]:
2470	refine = True
2471	line += ', sig:%8.4f'%(sizeSig[0][2])
2472	Snames = ['S11','S22','S33','S12','S13','S23']
2473	ptlbls = ' name :'
2474	ptstr = ' value :'
2475	sigstr = ' sig :'
2476	for i,name in enumerate(Snames):
2477	ptlbls += '%12s' % (name)
2478	ptstr += '%12.6f' % (hapData[4][i])
2479	if sizeSig[1][i]:
2480	refine = True
2481	sigstr += '%12.6f' % (sizeSig[1][i])
2482	else:
2483	sigstr += 12*' '
2484	if refine:
2485	print >>pFile,line
2486	print >>pFile,ptlbls
2487	print >>pFile,ptstr
2488	print >>pFile,sigstr
2489
2490	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2491	line = '\n Mustrain model: %9s'%(hapData[0])
2492	refine = False
2493	if hapData[0] in ['isotropic','uniaxial']:
2494	line += ' equatorial:%12.1f'%(hapData[1][0])
2495	if mustrainSig[0][0]:
2496	line += ', sig:%8.1f'%(mustrainSig[0][0])
2497	refine = True
2498	if hapData[0] == 'uniaxial':
2499	line += ' axial:%12.1f'%(hapData[1][1])
2500	if mustrainSig[0][1]:
2501	line += ', sig:%8.1f'%(mustrainSig[0][1])
2502	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2503	if mustrainSig[0][2]:
2504	refine = True
2505	line += ', sig:%8.3f'%(mustrainSig[0][2])
2506	if refine:
2507	print >>pFile,line
2508	else:
2509	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2510	if mustrainSig[0][2]:
2511	refine = True
2512	line += ', sig:%8.3f'%(mustrainSig[0][2])
2513	Snames = G2spc.MustrainNames(SGData)
2514	ptlbls = ' name :'
2515	ptstr = ' value :'
2516	sigstr = ' sig :'
2517	for i,name in enumerate(Snames):
2518	ptlbls += '%12s' % (name)
2519	ptstr += '%12.1f' % (hapData[4][i])
2520	if mustrainSig[1][i]:
2521	refine = True
2522	sigstr += '%12.1f' % (mustrainSig[1][i])
2523	else:
2524	sigstr += 12*' '
2525	if refine:
2526	print >>pFile,line
2527	print >>pFile,ptlbls
2528	print >>pFile,ptstr
2529	print >>pFile,sigstr
2530
2531	def PrintHStrainAndSig(hapData,strainSig,SGData):
2532	Hsnames = G2spc.HStrainNames(SGData)
2533	ptlbls = ' name :'
2534	ptstr = ' value :'
2535	sigstr = ' sig :'
2536	refine = False
2537	for i,name in enumerate(Hsnames):
2538	ptlbls += '%12s' % (name)
2539	ptstr += '%12.4g' % (hapData[0][i])
2540	if name in strainSig:
2541	refine = True
2542	sigstr += '%12.4g' % (strainSig[name])
2543	else:
2544	sigstr += 12*' '
2545	if refine:
2546	print >>pFile,'\n Hydrostatic/elastic strain: '
2547	print >>pFile,ptlbls
2548	print >>pFile,ptstr
2549	print >>pFile,sigstr
2550
2551	def PrintSHPOAndSig(pfx,hapData,POsig):
2552	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2553	ptlbls = ' names :'
2554	ptstr = ' values:'
2555	sigstr = ' sig :'
2556	for item in hapData[5]:
2557	ptlbls += '%12s'%(item)
2558	ptstr += '%12.3f'%(hapData[5][item])
2559	if pfx+item in POsig:
2560	sigstr += '%12.3f'%(POsig[pfx+item])
2561	else:
2562	sigstr += 12*' '
2563	print >>pFile,ptlbls
2564	print >>pFile,ptstr
2565	print >>pFile,sigstr
2566
2567	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2568	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2569	text = ''
2570	for item in hapData[2]:
2571	if pfx+item in ScalExtSig:
2572	text += ' %s: '%(item)
2573	text += '%12.2e'%(hapData[2][item][0])
2574	if pfx+item in ScalExtSig:
2575	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2576	print >>pFile,text
2577
2578	def PrintBabinetAndSig(pfx,hapData,BabSig):
2579	print >>pFile,'\n Babinet form factor modification: '
2580	ptlbls = ' names :'
2581	ptstr = ' values:'
2582	sigstr = ' sig :'
2583	for item in hapData:
2584	ptlbls += '%12s'%(item)
2585	ptstr += '%12.3f'%(hapData[item][0])
2586	if pfx+item in BabSig:
2587	sigstr += '%12.3f'%(BabSig[pfx+item])
2588	else:
2589	sigstr += 12*' '
2590	print >>pFile,ptlbls
2591	print >>pFile,ptstr
2592	print >>pFile,sigstr
2593
2594	def PrintTwinsAndSig(pfx,twinData,TwinSig):
2595	print >>pFile,'\n Twin Law fractions : '
2596	ptlbls = ' names :'
2597	ptstr = ' values:'
2598	sigstr = ' sig :'
2599	for it,item in enumerate(twinData):
2600	ptlbls += ' twin #%d'%(it)
2601	if it:
2602	ptstr += '%12.3f'%(item[1])
2603	else:
2604	ptstr += '%12.3f'%(item[1][0])
2605	if pfx+'TwinFr:'+str(it) in TwinSig:
2606	sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2607	else:
2608	sigstr += 12*' '
2609	print >>pFile,ptlbls
2610	print >>pFile,ptstr
2611	print >>pFile,sigstr
2612
2613
2614	PhFrExtPOSig = {}
2615	SizeMuStrSig = {}
2616	ScalExtSig = {}
2617	BabSig = {}
2618	TwinFrSig = {}
2619	wtFrSum = {}
2620	for phase in Phases:
2621	HistoPhase = Phases[phase]['Histograms']
2622	General = Phases[phase]['General']
2623	SGData = General['SGData']
2624	pId = Phases[phase]['pId']
2625	histoList = HistoPhase.keys()
2626	histoList.sort()
2627	for histogram in histoList:
2628	try:
2629	Histogram = Histograms[histogram]
2630	except KeyError:
2631	#skip if histogram not included e.g. in a sequential refinement
2632	continue
2633	hapData = HistoPhase[histogram]
2634	hId = Histogram['hId']
2635	pfx = str(pId)+':'+str(hId)+':'
2636	if hId not in wtFrSum:
2637	wtFrSum[hId] = 0.
2638	if 'PWDR' in histogram:
2639	for item in ['Scale','Extinction']:
2640	hapData[item][0] = parmDict[pfx+item]
2641	if pfx+item in sigDict:
2642	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2643	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2644	if hapData['Pref.Ori.'][0] == 'MD':
2645	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2646	if pfx+'MD' in sigDict:
2647	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2648	else: #'SH' spherical harmonics
2649	for item in hapData['Pref.Ori.'][5]:
2650	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2651	if pfx+item in sigDict:
2652	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2653	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2654	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2655	pfx+'HStrain':{}})
2656	for item in ['Mustrain','Size']:
2657	hapData[item][1][2] = parmDict[pfx+item+';mx']
2658	# hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2659	if pfx+item+';mx' in sigDict:
2660	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2661	if hapData[item][0] in ['isotropic','uniaxial']:
2662	hapData[item][1][0] = parmDict[pfx+item+';i']
2663	if item == 'Size':
2664	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2665	if pfx+item+';i' in sigDict:
2666	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2667	if hapData[item][0] == 'uniaxial':
2668	hapData[item][1][1] = parmDict[pfx+item+';a']
2669	if item == 'Size':
2670	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2671	if pfx+item+';a' in sigDict:
2672	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2673	else: #generalized for mustrain or ellipsoidal for size
2674	Nterms = len(hapData[item][4])
2675	for i in range(Nterms):
2676	sfx = ';'+str(i)
2677	hapData[item][4][i] = parmDict[pfx+item+sfx]
2678	if pfx+item+sfx in sigDict:
2679	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2680	names = G2spc.HStrainNames(SGData)
2681	for i,name in enumerate(names):
2682	hapData['HStrain'][0][i] = parmDict[pfx+name]
2683	if pfx+name in sigDict:
2684	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2685	for name in ['BabA','BabU']:
2686	hapData['Babinet'][name][0] = parmDict[pfx+name]
2687	if pfx+name in sigDict:
2688	BabSig[pfx+name] = sigDict[pfx+name]
2689
2690	elif 'HKLF' in histogram:
2691	for item in ['Scale','Flack']:
2692	if parmDict.get(pfx+item):
2693	hapData[item][0] = parmDict[pfx+item]
2694	if pfx+item in sigDict:
2695	ScalExtSig[pfx+item] = sigDict[pfx+item]
2696	for item in ['Ep','Eg','Es']:
2697	if parmDict.get(pfx+item):
2698	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2699	if pfx+item in sigDict:
2700	ScalExtSig[pfx+item] = sigDict[pfx+item]
2701	for item in ['BabA','BabU']:
2702	hapData['Babinet'][item][0] = parmDict[pfx+item]
2703	if pfx+item in sigDict:
2704	BabSig[pfx+item] = sigDict[pfx+item]
2705	if 'Twins' in hapData:
2706	it = 1
2707	sumTwFr = 0.
2708	while True:
2709	try:
2710	hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2711	if pfx+'TwinFr:'+str(it) in sigDict:
2712	TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2713	if it:
2714	sumTwFr += hapData['Twins'][it][1]
2715	it += 1
2716	except KeyError:
2717	break
2718	hapData['Twins'][0][1][0] = 1.-sumTwFr
2719
2720	if Print:
2721	for phase in Phases:
2722	HistoPhase = Phases[phase]['Histograms']
2723	General = Phases[phase]['General']
2724	SGData = General['SGData']
2725	pId = Phases[phase]['pId']
2726	histoList = HistoPhase.keys()
2727	histoList.sort()
2728	for histogram in histoList:
2729	try:
2730	Histogram = Histograms[histogram]
2731	except KeyError:
2732	#skip if histogram not included e.g. in a sequential refinement
2733	continue
2734	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2735	print >>pFile,130*'-'
2736	hapData = HistoPhase[histogram]
2737	hId = Histogram['hId']
2738	Histogram['Residuals'][str(pId)+'::Name'] = phase
2739	pfx = str(pId)+':'+str(hId)+':'
2740	hfx = ':%s:'%(hId)
2741	if 'PWDR' in histogram:
2742	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2743	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2744	print >>pFile,' Bragg intensity sum = %.3g'%(Histogram['Residuals'][pfx+'sumInt'])
2745
2746	if pfx+'Scale' in PhFrExtPOSig:
2747	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2748	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2749	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2750	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2751	if pfx+'Extinction' in PhFrExtPOSig:
2752	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2753	if hapData['Pref.Ori.'][0] == 'MD':
2754	if pfx+'MD' in PhFrExtPOSig:
2755	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2756	else:
2757	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2758	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2759	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2760	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2761	if len(BabSig):
2762	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2763
2764	elif 'HKLF' in histogram:
2765	Inst = Histogram['Instrument Parameters'][0]
2766	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)' \
2767	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2768	Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next'])
2769	if FFtables != None and 'N' not in Inst['Type'][0]:
2770	PrintFprime(FFtables,hfx,pFile)
2771	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2772	if pfx+'Scale' in ScalExtSig:
2773	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2774	if hapData['Extinction'][0] != 'None':
2775	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2776	if len(BabSig):
2777	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2778	if pfx+'Flack' in ScalExtSig:
2779	print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])
2780	if len(TwinFrSig):
2781	PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2782
2783	################################################################################
2784	##### Histogram data
2785	################################################################################
2786
2787	def GetHistogramData(Histograms,Print=True,pFile=None):
2788	'needs a doc string'
2789
2790	def GetBackgroundParms(hId,Background):
2791	Back = Background[0]
2792	DebyePeaks = Background[1]
2793	bakType,bakFlag = Back[:2]
2794	backVals = Back[3:]
2795	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2796	backDict = dict(zip(backNames,backVals))
2797	backVary = []
2798	if bakFlag:
2799	backVary = backNames
2800	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2801	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2802	debyeDict = {}
2803	debyeList = []
2804	for i in range(DebyePeaks['nDebye']):
2805	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2806	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2807	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2808	debyeVary = []
2809	for item in debyeList:
2810	if item[1]:
2811	debyeVary.append(item[0])
2812	backDict.update(debyeDict)
2813	backVary += debyeVary
2814	peakDict = {}
2815	peakList = []
2816	for i in range(DebyePeaks['nPeaks']):
2817	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2818	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2819	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2820	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2821	peakVary = []
2822	for item in peakList:
2823	if item[1]:
2824	peakVary.append(item[0])
2825	backDict.update(peakDict)
2826	backVary += peakVary
2827	return bakType,backDict,backVary
2828
2829	def GetInstParms(hId,Inst):
2830	dataType = Inst['Type'][0]
2831	instDict = {}
2832	insVary = []
2833	pfx = ':'+str(hId)+':'
2834	insKeys = Inst.keys()
2835	insKeys.sort()
2836	for item in insKeys:
2837	insName = pfx+item
2838	instDict[insName] = Inst[item][1]
2839	if len(Inst[item]) > 2 and Inst[item][2]:
2840	insVary.append(insName)
2841	if 'C' in dataType:
2842	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2843	return dataType,instDict,insVary
2844
2845	def GetSampleParms(hId,Sample):
2846	sampVary = []
2847	hfx = ':'+str(hId)+':'
2848	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2849	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2850	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2851	if key in Sample:
2852	sampDict[hfx+key] = Sample[key]
2853	Type = Sample['Type']
2854	if 'Bragg' in Type: #Bragg-Brentano
2855	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2856	sampDict[hfx+item] = Sample[item][0]
2857	if Sample[item][1]:
2858	sampVary.append(hfx+item)
2859	elif 'Debye' in Type: #Debye-Scherrer
2860	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2861	sampDict[hfx+item] = Sample[item][0]
2862	if Sample[item][1]:
2863	sampVary.append(hfx+item)
2864	return Type,sampDict,sampVary
2865
2866	def PrintBackground(Background):
2867	Back = Background[0]
2868	DebyePeaks = Background[1]
2869	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2870	line = ' Coefficients: '
2871	for i,back in enumerate(Back[3:]):
2872	line += '%10.3f'%(back)
2873	if i and not i%10:
2874	line += '\n'+15*' '
2875	print >>pFile,line
2876	if DebyePeaks['nDebye']:
2877	print >>pFile,'\n Debye diffuse scattering coefficients'
2878	parms = ['DebyeA','DebyeR','DebyeU']
2879	line = ' names : '
2880	for parm in parms:
2881	line += '%8s refine?'%(parm)
2882	print >>pFile,line
2883	for j,term in enumerate(DebyePeaks['debyeTerms']):
2884	line = ' term'+'%2d'%(j)+':'
2885	for i in range(3):
2886	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2887	print >>pFile,line
2888	if DebyePeaks['nPeaks']:
2889	print >>pFile,'\n Single peak coefficients'
2890	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2891	line = ' names : '
2892	for parm in parms:
2893	line += '%8s refine?'%(parm)
2894	print >>pFile,line
2895	for j,term in enumerate(DebyePeaks['peaksList']):
2896	line = ' peak'+'%2d'%(j)+':'
2897	for i in range(4):
2898	line += '%12.3f %5s'%(term[2i],bool(term[2i+1]))
2899	print >>pFile,line
2900
2901	def PrintInstParms(Inst):
2902	print >>pFile,'\n Instrument Parameters:'
2903	insKeys = Inst.keys()
2904	insKeys.sort()
2905	iBeg = 0
2906	Ok = True
2907	while Ok:
2908	ptlbls = ' name :'
2909	ptstr = ' value :'
2910	varstr = ' refine:'
2911	iFin = min(iBeg+9,len(insKeys))
2912	for item in insKeys[iBeg:iFin]:
2913	if item not in ['Type','Source']:
2914	ptlbls += '%12s' % (item)
2915	ptstr += '%12.6f' % (Inst[item][1])
2916	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2917	varstr += 12*' '
2918	else:
2919	varstr += '%12s' % (str(bool(Inst[item][2])))
2920	print >>pFile,ptlbls
2921	print >>pFile,ptstr
2922	print >>pFile,varstr
2923	iBeg = iFin
2924	if iBeg == len(insKeys):
2925	Ok = False
2926	else:
2927	print >>pFile,'\n'
2928
2929	def PrintSampleParms(Sample):
2930	print >>pFile,'\n Sample Parameters:'
2931	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2932	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2933	ptlbls = ' name :'
2934	ptstr = ' value :'
2935	varstr = ' refine:'
2936	if 'Bragg' in Sample['Type']:
2937	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2938	ptlbls += '%14s'%(item)
2939	ptstr += '%14.4f'%(Sample[item][0])
2940	varstr += '%14s'%(str(bool(Sample[item][1])))
2941
2942	elif 'Debye' in Type: #Debye-Scherrer
2943	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2944	ptlbls += '%14s'%(item)
2945	ptstr += '%14.4f'%(Sample[item][0])
2946	varstr += '%14s'%(str(bool(Sample[item][1])))
2947
2948	print >>pFile,ptlbls
2949	print >>pFile,ptstr
2950	print >>pFile,varstr
2951
2952	histDict = {}
2953	histVary = []
2954	controlDict = {}
2955	histoList = Histograms.keys()
2956	histoList.sort()
2957	for histogram in histoList:
2958	if 'PWDR' in histogram:
2959	Histogram = Histograms[histogram]
2960	hId = Histogram['hId']
2961	pfx = ':'+str(hId)+':'
2962	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2963	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2964	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2965	for excl in controlDict[pfx+'Exclude']:
2966	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2967	if controlDict[pfx+'Exclude']:
2968	ma.mask_rows(Histogram['Data'])
2969	Background = Histogram['Background']
2970	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2971	controlDict[pfx+'bakType'] = Type
2972	histDict.update(bakDict)
2973	histVary += bakVary
2974
2975	Inst = Histogram['Instrument Parameters'][0]
2976	Type,instDict,insVary = GetInstParms(hId,Inst)
2977	controlDict[pfx+'histType'] = Type
2978	if 'XC' in Type:
2979	if pfx+'Lam1' in instDict:
2980	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2981	else:
2982	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2983	histDict.update(instDict)
2984	histVary += insVary
2985
2986	Sample = Histogram['Sample Parameters']
2987	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2988	controlDict[pfx+'instType'] = Type
2989	histDict.update(sampDict)
2990	histVary += sampVary
2991
2992
2993	if Print:
2994	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2995	print >>pFile,135*'-'
2996	Units = {'C':' deg','T':' msec'}
2997	units = Units[controlDict[pfx+'histType'][2]]
2998	Limits = controlDict[pfx+'Limits']
2999	print >>pFile,' Instrument type: ',Sample['Type']
3000	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
3001	if len(controlDict[pfx+'Exclude']):
3002	excls = controlDict[pfx+'Exclude']
3003	for excl in excls:
3004	print >>pFile,' Excluded region: %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)
3005	PrintSampleParms(Sample)
3006	PrintInstParms(Inst)
3007	PrintBackground(Background)
3008	elif 'HKLF' in histogram:
3009	Histogram = Histograms[histogram]
3010	hId = Histogram['hId']
3011	pfx = ':'+str(hId)+':'
3012	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3013	Inst = Histogram['Instrument Parameters'][0]
3014	controlDict[pfx+'histType'] = Inst['Type'][0]
3015	if 'X' in Inst['Type'][0]:
3016	histDict[pfx+'Lam'] = Inst['Lam'][1]
3017	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3018	return histVary,histDict,controlDict
3019
3020	def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3021	'needs a doc string'
3022
3023	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3024	Back = Background[0]
3025	DebyePeaks = Background[1]
3026	lenBack = len(Back[3:])
3027	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
3028	for i in range(lenBack):
3029	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3030	if pfx+'Back;'+str(i) in sigDict:
3031	backSig[i] = sigDict[pfx+'Back;'+str(i)]
3032	if DebyePeaks['nDebye']:
3033	for i in range(DebyePeaks['nDebye']):
3034	names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3035	for j,name in enumerate(names):
3036	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3037	if name in sigDict:
3038	backSig[lenBack+3*i+j] = sigDict[name]
3039	if DebyePeaks['nPeaks']:
3040	for i in range(DebyePeaks['nPeaks']):
3041	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3042	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3043	for j,name in enumerate(names):
3044	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3045	if name in sigDict:
3046	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
3047	return backSig
3048
3049	def SetInstParms(pfx,Inst,parmDict,sigDict):
3050	instSig = {}
3051	insKeys = Inst.keys()
3052	insKeys.sort()
3053	for item in insKeys:
3054	insName = pfx+item
3055	Inst[item][1] = parmDict[insName]
3056	if insName in sigDict:
3057	instSig[item] = sigDict[insName]
3058	else:
3059	instSig[item] = 0
3060	return instSig
3061
3062	def SetSampleParms(pfx,Sample,parmDict,sigDict):
3063	if 'Bragg' in Sample['Type']: #Bragg-Brentano
3064	sampSig = [0 for i in range(5)]
3065	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3066	Sample[item][0] = parmDict[pfx+item]
3067	if pfx+item in sigDict:
3068	sampSig[i] = sigDict[pfx+item]
3069	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3070	sampSig = [0 for i in range(4)]
3071	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3072	Sample[item][0] = parmDict[pfx+item]
3073	if pfx+item in sigDict:
3074	sampSig[i] = sigDict[pfx+item]
3075	return sampSig
3076
3077	def PrintBackgroundSig(Background,backSig):
3078	Back = Background[0]
3079	DebyePeaks = Background[1]
3080	lenBack = len(Back[3:])
3081	valstr = ' value : '
3082	sigstr = ' sig : '
3083	refine = False
3084	for i,back in enumerate(Back[3:]):
3085	valstr += '%10.4g'%(back)
3086	if Back[1]:
3087	refine = True
3088	sigstr += '%10.4g'%(backSig[i])
3089	else:
3090	sigstr += 10*' '
3091	if refine:
3092	print >>pFile,'\n Background function: ',Back[0]
3093	print >>pFile,valstr
3094	print >>pFile,sigstr
3095	if DebyePeaks['nDebye']:
3096	ifAny = False
3097	ptfmt = "%12.3f"
3098	names = ' names :'
3099	ptstr = ' values:'
3100	sigstr = ' esds :'
3101	for item in sigDict:
3102	if 'Debye' in item:
3103	ifAny = True
3104	names += '%12s'%(item)
3105	ptstr += ptfmt%(parmDict[item])
3106	sigstr += ptfmt%(sigDict[item])
3107	if ifAny:
3108	print >>pFile,'\n Debye diffuse scattering coefficients'
3109	print >>pFile,names
3110	print >>pFile,ptstr
3111	print >>pFile,sigstr
3112	if DebyePeaks['nPeaks']:
3113	print >>pFile,'\n Single peak coefficients:'
3114	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3115	line = ' peak no. '
3116	for parm in parms:
3117	line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3118	print >>pFile,line
3119	for ip in range(DebyePeaks['nPeaks']):
3120	ptstr = ' %4d '%(ip)
3121	for parm in parms:
3122	fmt = '%14.3f'
3123	efmt = '%12.3f'
3124	if parm == 'BkPkpos':
3125	fmt = '%14.4f'
3126	efmt = '%12.4f'
3127	name = pfx+parm+';%d'%(ip)
3128	ptstr += fmt%(parmDict[name])
3129	if name in sigDict:
3130	ptstr += efmt%(sigDict[name])
3131	else:
3132	ptstr += 12*' '
3133	print >>pFile,ptstr
3134	sumBk = np.array(Histogram['sumBk'])
3135	print >>pFile,' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g' \
3136	%(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk))
3137
3138	def PrintInstParmsSig(Inst,instSig):
3139	refine = False
3140	insKeys = instSig.keys()
3141	insKeys.sort()
3142	iBeg = 0
3143	Ok = True
3144	while Ok:
3145	ptlbls = ' names :'
3146	ptstr = ' value :'
3147	sigstr = ' sig :'
3148	iFin = min(iBeg+9,len(insKeys))
3149	for name in insKeys[iBeg:iFin]:
3150	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3151	ptlbls += '%12s' % (name)
3152	ptstr += '%12.6f' % (Inst[name][1])
3153	if instSig[name]:
3154	refine = True
3155	sigstr += '%12.6f' % (instSig[name])
3156	else:
3157	sigstr += 12*' '
3158	if refine:
3159	print >>pFile,'\n Instrument Parameters:'
3160	print >>pFile,ptlbls
3161	print >>pFile,ptstr
3162	print >>pFile,sigstr
3163	iBeg = iFin
3164	if iBeg == len(insKeys):
3165	Ok = False
3166
3167	def PrintSampleParmsSig(Sample,sampleSig):
3168	ptlbls = ' names :'
3169	ptstr = ' values:'
3170	sigstr = ' sig :'
3171	refine = False
3172	if 'Bragg' in Sample['Type']:
3173	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3174	ptlbls += '%14s'%(item)
3175	ptstr += '%14.4f'%(Sample[item][0])
3176	if sampleSig[i]:
3177	refine = True
3178	sigstr += '%14.4f'%(sampleSig[i])
3179	else:
3180	sigstr += 14*' '
3181
3182	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3183	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3184	ptlbls += '%14s'%(item)
3185	ptstr += '%14.4f'%(Sample[item][0])
3186	if sampleSig[i]:
3187	refine = True
3188	sigstr += '%14.4f'%(sampleSig[i])
3189	else:
3190	sigstr += 14*' '
3191
3192	if refine:
3193	print >>pFile,'\n Sample Parameters:'
3194	print >>pFile,ptlbls
3195	print >>pFile,ptstr
3196	print >>pFile,sigstr
3197
3198	histoList = Histograms.keys()
3199	histoList.sort()
3200	for histogram in histoList:
3201	if 'PWDR' in histogram:
3202	Histogram = Histograms[histogram]
3203	hId = Histogram['hId']
3204	pfx = ':'+str(hId)+':'
3205	Background = Histogram['Background']
3206	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3207
3208	Inst = Histogram['Instrument Parameters'][0]
3209	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3210
3211	Sample = Histogram['Sample Parameters']
3212	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3213
3214	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
3215	print >>pFile,135*'-'
3216	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
3217	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
3218	print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
3219	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
3220	if Print:
3221	print >>pFile,' Instrument type: ',Sample['Type']
3222	if FFtables != None and 'N' not in Inst['Type'][0]:
3223	PrintFprime(FFtables,pfx,pFile)
3224	PrintSampleParmsSig(Sample,sampSig)
3225	PrintInstParmsSig(Inst,instSig)
3226	PrintBackgroundSig(Background,backSig)
3227

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

Download in other formats: