Context Navigation

GSASIIstrIO.py @ 2964

Last change on this file since 2964 was 2964, checked in by vondreele, 6 years ago
work on protein validation - not correct yet; gives plot tho. fix refinement of Dij values in a LeBail? refinement (esp. sequential LeBail?).
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 147.3 KB

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

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

Context Navigation

source: branch/2frame/GSASIIstrIO.py @ 2964

Download in other formats: