Context Navigation

source: trunk/GSASIIstrIO.py @ 4369

Last change on this file since 4369 was 4197, checked in by vondreele, 6 years ago
put reference to Cromer & Liberman in Absorb & Fprime as well as refinement LST output & atmdata.py
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 161.0 KB

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

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

Download in other formats: