Context Navigation

source: trunk/GSASIIstrIO.py @ 4629

Last change on this file since 4629 was 4629, checked in by vondreele, 3 years ago
GetXsectionCoeff? - remove a 'U' from an open - deprecated close infile & outfile in SetUsedHistogramsAndPhases?
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 168.0 KB

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

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

Download in other formats: