Context Navigation

source: trunk/GSASIIstrIO.py @ 3984

Last change on this file since 3984 was 3925, checked in by toby, 6 years ago
bugfix for fixed background: conflict for single xtal
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 160.2 KB

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

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

Download in other formats: