Context Navigation

source: trunk/GSASIIstrIO.py @ 2560

Last change on this file since 2560 was 2560, checked in by vondreele, 7 years ago
make A2Gmat return np.arrays put in spinners for Background & Container multipliers fix a couple of Status problems fix(?) lattice parameter errors for hexagonal & trigonal
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 145.0 KB

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

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

Download in other formats: