Context Navigation

source: trunk/GSASIIstrIO.py @ 3924

Last change on this file since 3924 was 3906, checked in by toby, 6 years ago
implement fixed backgound file in Riteveld Fit; show Rietveld & peak fits with bkg file added to calc bkg, not subtracted
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 160.1 KB

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

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

Download in other formats: