Context Navigation

source: trunk/GSASIIstrIO.py @ 4823

Last change on this file since 4823 was 4823, checked in by toby, 2 years ago
fix naming of new variable constraints; label them & improve constraint help text
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 171.8 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2021-02-20 19:25:23 +0000 (Sat, 20 Feb 2021) $
4	# $Author: toby $
5	# $Revision: 4823 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 4823 2021-02-20 19:25:23Z toby $
8	########### SVN repository information ###################
9	'''
10	GSASIIstrIO: structure I/O routines
11	-------------------------------------
12
13	Contains routines for reading from GPX files and printing to the .LST file.
14	Used for refinements and in G2scriptable.
15
16	Should not contain any wxpython references as this should be able to be used
17	in non-GUI settings.
18
19	'''
20	from __future__ import division, print_function
21	import platform
22	import os
23	import os.path as ospath
24	import time
25	import math
26	import copy
27	if '2' in platform.python_version_tuple()[0]:
28	import cPickle
29	else:
30	import pickle as cPickle
31	import numpy as np
32	import numpy.ma as ma
33	import GSASIIpath
34	GSASIIpath.SetVersionNumber("$Revision: 4823 $")
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	namedVarList = []
194	for item in constList:
195	if item[-1] == 'h':
196	# process a hold
197	fixedList.append('0')
198	var = str(item[0][1])
199	if '?' not in var:
200	constDict.append({var:0.0})
201	else:
202	ignored += 1
203	elif item[-1] == 'f':
204	# process a new variable
205	fixedList.append(None)
206	D = {}
207	varyFlag = item[-2]
208	varname = item[-3]
209	for term in item[:-3]:
210	var = str(term[1])
211	if '?' not in var:
212	D[var] = term[0]
213	if len(D) > 1:
214	# add extra dict terms for input variable name and vary flag
215	if varname is not None:
216	varname = str(varname) # in case this is a G2VarObj
217	if varname.startswith(':'):
218	D['_name'] = varname
219	else:
220	D['_name'] = '::nv-' + varname
221	D['_name'] = G2obj.MakeUniqueLabel(D['_name'],namedVarList)
222	D['_vary'] = varyFlag == True # force to bool
223	constDict.append(D)
224	else:
225	ignored += 1
226	#constFlag[-1] = ['Vary']
227	elif item[-1] == 'c':
228	# process a contraint relationship
229	D = {}
230	for term in item[:-3]:
231	var = str(term[1])
232	if '?' not in var:
233	D[var] = term[0]
234	if len(D) >= 1:
235	fixedList.append(str(item[-3]))
236	constDict.append(D)
237	else:
238	ignored += 1
239	elif item[-1] == 'e':
240	# process an equivalence
241	firstmult = None
242	eqlist = []
243	for term in item[:-3]:
244	if term[0] == 0: term[0] = 1.0
245	var = str(term[1])
246	if '?' in var: continue
247	if firstmult is None:
248	firstmult = term[0]
249	firstvar = var
250	else:
251	eqlist.append([var,firstmult/term[0]])
252	if len(eqlist) > 0:
253	G2mv.StoreEquivalence(firstvar,eqlist,False)
254	else:
255	ignored += 1
256	else:
257	ignored += 1
258	return constDict,fixedList,ignored
259
260	def ReadCheckConstraints(GPXfile, seqHist=None):
261	'''Load constraints and related info and return any error or warning messages
262	This is done from the GPX file rather than the tree.
263
264	:param dict seqHist: defines a specific histogram to be loaded for a sequential
265	refinement, if None (default) all are loaded.
266	'''
267	# init constraints
268	G2mv.InitVars()
269	# get variables
270	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
271	if not Phases:
272	return 'Error: No phases or no histograms in phases!',''
273	if not Histograms:
274	return 'Error: no diffraction data',''
275	constrDict,fixedList = GetConstraints(GPXfile) # load user constraints before internally generated ones
276	if seqHist: Histograms = {seqHist:Histograms[seqHist]} # sequential fit: only need one histogram
277	rigidbodyDict = GetRigidBodies(GPXfile)
278	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
279	rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
280	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
281	GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
282	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False,resetRefList=False)
283	histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
284	varyList = rbVary+phaseVary+hapVary+histVary
285	msg = G2mv.EvaluateMultipliers(constrDict,phaseDict,hapDict,histDict)
286	if msg:
287	return 'Unable to interpret multiplier(s): '+msg,''
288	errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
289	if errmsg:
290	# print some diagnostic info on the constraints
291	G2fil.G2Print('Error in constraints:\n'+errmsg+
292	'\nRefinement not possible due to conflict in constraints, see below:\n')
293	G2fil.G2Print(G2mv.VarRemapShow(varyList,True),mode='error')
294	return errmsg, warnmsg
295
296	def makeTwinFrConstr(Phases,Histograms,hapVary):
297	TwConstr = []
298	TwFixed = []
299	for Phase in Phases:
300	pId = Phases[Phase]['pId']
301	for Histogram in Phases[Phase]['Histograms']:
302	try:
303	hId = Histograms[Histogram]['hId']
304	phfx = '%d:%d:'%(pId,hId)
305	if phfx+'TwinFr:0' in hapVary:
306	TwFixed.append('1.0') #constraint value
307	nTwin = len(Phases[Phase]['Histograms'][Histogram]['Twins'])
308	TwConstr.append({phfx+'TwinFr:'+str(i):'1.0' for i in range(nTwin)})
309	except KeyError: #unused histograms?
310	pass
311	return TwConstr,TwFixed
312
313	def GetRestraints(GPXfile):
314	'''Read the restraints from the GPX file.
315	Throws an exception if not found in the .GPX file
316	'''
317	IndexGPX(GPXfile)
318	fl = open(GPXfile,'rb')
319	pos = gpxIndex.get('Restraints')
320	if pos is None:
321	raise Exception("No Restraints in GPX file")
322	fl.seek(pos)
323	datum = cPickleLoad(fl)[0]
324	fl.close()
325	return datum[1]
326
327	def GetRigidBodies(GPXfile):
328	'''Read the rigid body models from the GPX file
329	'''
330	IndexGPX(GPXfile)
331	fl = open(GPXfile,'rb')
332	pos = gpxIndex.get('Rigid bodies')
333	if pos is None:
334	raise Exception("No Rigid bodies in GPX file")
335	fl.seek(pos)
336	datum = cPickleLoad(fl)[0]
337	fl.close()
338	return datum[1]
339
340	def GetFprime(controlDict,Histograms):
341	'Needs a doc string'
342	FFtables = controlDict['FFtables']
343	if not FFtables:
344	return
345	histoList = list(Histograms.keys())
346	histoList.sort()
347	for histogram in histoList:
348	if histogram[:4] in ['PWDR','HKLF']:
349	Histogram = Histograms[histogram]
350	hId = Histogram['hId']
351	hfx = ':%d:'%(hId)
352	if 'X' in controlDict[hfx+'histType']:
353	keV = controlDict[hfx+'keV']
354	for El in FFtables:
355	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
356	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
357	FFtables[El][hfx+'FP'] = FP
358	FFtables[El][hfx+'FPP'] = FPP
359
360	def PrintFprime(FFtables,pfx,pFile):
361	pFile.write('\n Resonant form factors:(ref: D.T. Cromer & D.A. Liberman (1981), Acta Cryst. A37, 267-268.)\n')
362	Elstr = ' Element:'
363	FPstr = " f' :"
364	FPPstr = ' f" :'
365	for El in FFtables:
366	Elstr += ' %8s'%(El)
367	FPstr += ' %8.3f'%(FFtables[El][pfx+'FP'])
368	FPPstr += ' %8.3f'%(FFtables[El][pfx+'FPP'])
369	pFile.write(Elstr+'\n')
370	pFile.write(FPstr+'\n')
371	pFile.write(FPPstr+'\n')
372
373	def PrintBlength(BLtables,wave,pFile):
374	pFile.write('\n Resonant neutron scattering lengths:\n')
375	Elstr = ' Element:'
376	FPstr = " b' :"
377	FPPstr = ' b" :'
378	for El in BLtables:
379	BP,BPP = G2el.BlenResCW([El,],BLtables,wave)
380	Elstr += ' %8s'%(El)
381	FPstr += ' %8.3f'%(BP)
382	FPPstr += ' %8.3f'%(BPP)
383	pFile.write(Elstr+'\n')
384	pFile.write(FPstr+'\n')
385	pFile.write(FPPstr+'\n')
386
387	def PrintISOmodes(pFile,Phases,parmDict,sigDict):
388	'''Prints the values for the ISODISTORT modes into the project's
389	.lst file after a refinement.
390	'''
391	for phase in Phases:
392	if 'ISODISTORT' not in Phases[phase]: continue
393	data = Phases[phase]
394	ISO = data['ISODISTORT']
395	if 'G2VarList' in ISO:
396	deltaList = []
397	for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
398	dvar = gv.varname()
399	var = dvar.replace('::dA','::A')
400	albl = Ilbl[:Ilbl.rfind('_')]
401	v = Ilbl[Ilbl.rfind('_')+1:]
402	pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
403	if var in parmDict:
404	cval = parmDict[var]
405	else:
406	print('PrintISOmodes Error: Atom not found',"No value found for parameter "+str(var))
407	return
408	deltaList.append(cval-pval)
409	modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
410
411	if 'G2OccVarList' in ISO:
412	deltaOccList = []
413	for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
414	var = gv.varname()
415	albl = Ilbl[:Ilbl.rfind('_')]
416	pval = ISO['BaseOcc'][albl]
417	if var in parmDict:
418	cval = parmDict[var]
419	else:
420	print('PrintISOmodes Error: Atom not found',"No value found for parameter "+str(var))
421	return
422	deltaOccList.append(cval-pval)
423	modeOccVals = np.inner(ISO['Var2OccMatrix'],deltaOccList)
424
425	if 'G2VarList' in ISO:
426	pFile.write('\n ISODISTORT Displacive Modes for phase {}\n'.format(data['General'].get('Name','')))
427	l = str(max([len(i) for i in ISO['IsoModeList']])+3)
428	fmt = ' {:'+l+'}{}'
429	for var,val,norm,G2mode in zip(
430	ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] ):
431	try:
432	value = G2py3.FormatSigFigs(val/norm)
433	if str(G2mode) in sigDict:
434	value = G2mth.ValEsd(val/norm,sigDict[str(G2mode)]/norm)
435	except TypeError:
436	value = '?'
437	pFile.write(fmt.format(var,value)+'\n')
438
439	if 'G2OccVarList' in ISO:
440	pFile.write('\n ISODISTORT Occupancy Modes for phase {}\n'.format(data['General'].get('Name','')))
441	l = str(max([len(i) for i in ISO['OccModeList']])+3)
442	fmt = ' {:'+l+'}{}'
443	for var,val,norm,G2mode in zip(
444	ISO['OccModeList'],modeOccVals,ISO['OccNormList'],ISO['G2OccModeList'] ):
445	try:
446	value = G2py3.FormatSigFigs(val/norm)
447	if str(G2mode) in sigDict:
448	value = G2mth.ValEsd(val/norm,sigDict[str(G2mode)]/norm)
449	except TypeError:
450	value = '?'
451	pFile.write(fmt.format(var,value)+'\n')
452
453	def GetPhaseNames(GPXfile):
454	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
455
456	:param str GPXfile: full .gpx file name
457	:return: list of phase names
458	'''
459	IndexGPX(GPXfile)
460	fl = open(GPXfile,'rb')
461	pos = gpxIndex.get('Phases')
462	if pos is None:
463	raise Exception("No Phases in GPX file")
464	fl.seek(pos)
465	data = cPickleLoad(fl)
466	fl.close()
467	return [datus[0] for datus in data[1:]]
468
469	def GetAllPhaseData(GPXfile,PhaseName):
470	''' Returns the entire dictionary for PhaseName from GSASII gpx file
471
472	:param str GPXfile: full .gpx file name
473	:param str PhaseName: phase name
474	:return: phase dictionary or None if PhaseName is not present
475	'''
476	IndexGPX(GPXfile)
477	fl = open(GPXfile,'rb')
478	pos = gpxIndex.get('Phases')
479	if pos is None:
480	raise Exception("No Phases in GPX file")
481	fl.seek(pos)
482	data = cPickleLoad(fl)
483	fl.close()
484
485	for datus in data[1:]:
486	if datus[0] == PhaseName:
487	return datus[1]
488
489	def GetHistograms(GPXfile,hNames):
490	""" Returns a dictionary of histograms found in GSASII gpx file
491
492	:param str GPXfile: full .gpx file name
493	:param str hNames: list of histogram names
494	:return: dictionary of histograms (types = PWDR & HKLF)
495
496	"""
497	IndexGPX(GPXfile)
498	fl = open(GPXfile,'rb')
499	Histograms = {}
500	for hist in hNames:
501	pos = gpxIndex.get(hist)
502	if pos is None:
503	raise Exception("Histogram {} not found in GPX file".format(hist))
504	fl.seek(pos)
505	data = cPickleLoad(fl)
506	datum = data[0]
507	if 'PWDR' in hist[:4]:
508	PWDRdata = {}
509	PWDRdata.update(datum[1][0]) #weight factor
510	PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1])) #masked powder data arrays/clear previous masks
511	PWDRdata[data[2][0]] = data[2][1] #Limits & excluded regions (if any)
512	PWDRdata[data[3][0]] = data[3][1] #Background
513	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
514	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
515	try:
516	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
517	except IndexError:
518	PWDRdata['Reflection Lists'] = {}
519	PWDRdata['Residuals'] = {}
520	Histograms[hist] = PWDRdata
521	elif 'HKLF' in hist[:4]:
522	HKLFdata = {}
523	HKLFdata.update(datum[1][0]) #weight factor
524	#patch
525	if 'list' in str(type(datum[1][1])):
526	#if isinstance(datum[1][1],list):
527	RefData = {'RefList':[],'FF':{}}
528	for ref in datum[1][1]:
529	RefData['RefList'].append(ref[:11]+[ref[13],])
530	RefData['RefList'] = np.array(RefData['RefList'])
531	datum[1][1] = RefData
532	#end patch
533	datum[1][1]['FF'] = {}
534	HKLFdata['Data'] = datum[1][1]
535	HKLFdata[data[1][0]] = data[1][1] #Instrument parameters
536	HKLFdata['Reflection Lists'] = None
537	HKLFdata['Residuals'] = {}
538	Histograms[hist] = HKLFdata
539	fl.close()
540	return Histograms
541
542	def GetHistogramNames(GPXfile,hTypes):
543	""" Returns a list of histogram names found in a GSAS-II .gpx file that
544	match specifed histogram types. Names are returned in the order they
545	appear in the file.
546
547	:param str GPXfile: full .gpx file name
548	:param str hTypes: list of histogram types
549	:return: list of histogram names (types = PWDR & HKLF)
550
551	"""
552	IndexGPX(GPXfile)
553	return [n[0] for n in gpxNamelist if n[0][:4] in hTypes]
554
555	def GetUsedHistogramsAndPhases(GPXfile):
556	''' Returns all histograms that are found in any phase
557	and any phase that uses a histogram. This also
558	assigns numbers to used phases and histograms by the
559	order they appear in the file.
560
561	:param str GPXfile: full .gpx file name
562	:returns: (Histograms,Phases)
563
564	* Histograms = dictionary of histograms as {name:data,...}
565	* Phases = dictionary of phases that use histograms
566
567	'''
568	phaseNames = GetPhaseNames(GPXfile)
569	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
570	allHistograms = GetHistograms(GPXfile,histoList)
571	phaseData = {}
572	for name in phaseNames:
573	phaseData[name] = GetAllPhaseData(GPXfile,name)
574	Histograms = {}
575	Phases = {}
576	for phase in phaseData:
577	Phase = phaseData[phase]
578	if Phase['General']['Type'] == 'faulted': continue #don't use faulted phases!
579	if Phase['Histograms']:
580	for hist in Phase['Histograms']:
581	if 'Use' not in Phase['Histograms'][hist]: #patch
582	Phase['Histograms'][hist]['Use'] = True
583	if Phase['Histograms'][hist]['Use'] and phase not in Phases:
584	pId = phaseNames.index(phase)
585	Phase['pId'] = pId
586	Phases[phase] = Phase
587	if hist not in Histograms and Phase['Histograms'][hist]['Use']:
588	try:
589	Histograms[hist] = allHistograms[hist]
590	hId = histoList.index(hist)
591	Histograms[hist]['hId'] = hId
592	except KeyError: # would happen if a referenced histogram were
593	# renamed or deleted
594	G2fil.G2Print('Warning: For phase "'+phase+
595	'" unresolved reference to histogram "'+hist+'"')
596	# load the fix background info into the histograms
597	for hist in Histograms:
598	if 'Background' not in Histograms[hist]: continue
599	fixedBkg = Histograms[hist]['Background'][1].get('background PWDR')
600	if fixedBkg:
601	if not fixedBkg[0]: continue
602	# patch: add refinement flag, if needed
603	if len(fixedBkg) == 2: fixedBkg += [False]
604	h = Histograms[hist]['Background'][1]
605	try:
606	Limits = Histograms[hist]['Limits'][1]
607	x = Histograms[hist]['Data'][0]
608	xB = np.searchsorted(x,Limits[0])
609	xF = np.searchsorted(x,Limits[1])+1
610	h['fixback'] = allHistograms[fixedBkg[0]]['Data'][1][xB:xF]
611	except KeyError: # would happen if a referenced histogram were renamed or deleted
612	G2fil.G2Print('Warning: For hist "{}" unresolved background reference to hist "{}"'
613	.format(hist,fixedBkg[0]))
614	G2obj.IndexAllIds(Histograms=Histograms,Phases=Phases)
615	return Histograms,Phases
616
617	def getBackupName(GPXfile,makeBack):
618	'''
619	Get the name for the backup .gpx file name
620
621	:param str GPXfile: full .gpx file name
622	:param bool makeBack: if True the name of a new file is returned, if
623	False the name of the last file that exists is returned
624	:returns: the name of a backup file
625
626	'''
627	GPXpath,GPXname = ospath.split(GPXfile)
628	if GPXpath == '': GPXpath = '.'
629	Name = ospath.splitext(GPXname)[0]
630	files = os.listdir(GPXpath)
631	last = 0
632	for name in files:
633	name = name.split('.')
634	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
635	if makeBack:
636	last = max(last,int(name[1].strip('bak'))+1)
637	else:
638	last = max(last,int(name[1].strip('bak')))
639	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
640	return GPXback
641
642	def GPXBackup(GPXfile,makeBack=True):
643	'''
644	makes a backup of the specified .gpx file
645
646	:param str GPXfile: full .gpx file name
647	:param bool makeBack: if True (default), the backup is written to
648	a new file; if False, the last backup is overwritten
649	:returns: the name of the backup file that was written
650	'''
651	import distutils.file_util as dfu
652	GPXback = getBackupName(GPXfile,makeBack)
653	tries = 0
654	while True:
655	try:
656	dfu.copy_file(GPXfile,GPXback)
657	break
658	except:
659	tries += 1
660	if tries > 10:
661	return GPXfile #failed!
662	time.sleep(1) #just wait a second!
663	return GPXback
664
665	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,parmFrozenList,makeBack=True):
666	''' Updates gpxfile from all histograms that are found in any phase
667	and any phase that used a histogram. Also updates rigid body definitions.
668	This is used for non-sequential fits, but not for sequential fitting.
669
670	:param str GPXfile: full .gpx file name
671	:param dict Histograms: dictionary of histograms as {name:data,...}
672	:param dict Phases: dictionary of phases that use histograms
673	:param dict RigidBodies: dictionary of rigid bodies
674	:param dict CovData: dictionary of refined variables, varyList, & covariance matrix
675	:param list parmFrozenList: list of parameters (as str) that are frozen
676	due to limits; converted to :class:`GSASIIobj.G2VarObj` objects.
677	:param bool makeBack: True if new backup of .gpx file is to be made; else
678	use the last one made
679	'''
680
681	import distutils.file_util as dfu
682	GPXback = GPXBackup(GPXfile,makeBack)
683	G2fil.G2Print ('Read from file:'+GPXback)
684	G2fil.G2Print ('Save to file :'+GPXfile)
685	infile = open(GPXback,'rb')
686	outfile = open(GPXfile,'wb')
687	while True:
688	try:
689	data = cPickleLoad(infile)
690	except EOFError:
691	break
692	datum = data[0]
693	# print 'read: ',datum[0]
694	if datum[0] == 'Phases':
695	for iphase in range(len(data)):
696	if data[iphase][0] in Phases:
697	phaseName = data[iphase][0]
698	data[iphase][1].update(Phases[phaseName])
699	elif datum[0] == 'Covariance':
700	data[0][1] = CovData
701	elif datum[0] == 'Rigid bodies':
702	data[0][1] = RigidBodies
703	elif datum[0] == 'Controls':
704	Controls = data[0][1]
705	if 'parmFrozen' not in Controls:
706	Controls['parmFrozen'] = {}
707	Controls['parmFrozen']['FrozenList'] = [
708	i if type(i) is G2obj.G2VarObj
709	else G2obj.G2VarObj(i)
710	for i in parmFrozenList]
711	try:
712	histogram = Histograms[datum[0]]
713	# print 'found ',datum[0]
714	data[0][1][1] = histogram['Data']
715	data[0][1][0].update(histogram['Residuals'])
716	for datus in data[1:]:
717	# print ' read: ',datus[0]
718	if datus[0] in ['Instrument Parameters','Sample Parameters','Reflection Lists']:
719	datus[1] = histogram[datus[0]]
720	if datus[0] == 'Background': # remove fixed background from file
721	d1 = {key:histogram['Background'][1][key]
722	for key in histogram['Background'][1]
723	if not key.startswith('_fixed')}
724	datus[1] = copy.deepcopy(histogram['Background'])
725	datus[1][1] = d1
726	except KeyError:
727	pass
728	try:
729	cPickle.dump(data,outfile,1)
730	except AttributeError:
731	G2fil.G2Print ('ERROR - bad data in least squares result')
732	infile.close()
733	outfile.close()
734	dfu.copy_file(GPXback,GPXfile)
735	G2fil.G2Print ('GPX file save failed - old version retained',mode='error')
736	return
737
738	infile.close()
739	outfile.close()
740
741	G2fil.G2Print ('GPX file save successful')
742
743	def GetSeqResult(GPXfile):
744	'''
745	Returns the sequential results table information from a GPX file.
746	Called at the beginning of :meth:`GSASIIstrMain.SeqRefine`
747
748	:param str GPXfile: full .gpx file name
749	:returns: a dict containing the sequential results table
750	'''
751	IndexGPX(GPXfile)
752	pos = gpxIndex.get('Sequential results')
753	if pos is None:
754	return {}
755	fl = open(GPXfile,'rb')
756	fl.seek(pos)
757	datum = cPickleLoad(fl)[0]
758	fl.close()
759	return datum[1]
760
761	def SetupSeqSavePhases(GPXfile):
762	'''Initialize the files used to save intermediate results from
763	sequential fits.
764	'''
765	IndexGPX(GPXfile)
766	# load initial Phase results from GPX
767	fl = open(GPXfile,'rb')
768	pos = gpxIndex.get('Phases')
769	if pos is None:
770	raise Exception("No Phases in GPX file")
771	fl.seek(pos)
772	data = cPickleLoad(fl)
773	fl.close()
774	# create GPX-like file to store latest Phase info; init with start vals
775	GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
776	fp = open(GPXphase,'wb')
777	cPickle.dump(data,fp,1)
778	fp.close()
779	# create empty file for histogram info
780	GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
781	fp = open(GPXhist,'wb')
782	fp.close()
783
784	def SaveUpdatedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,parmFrozen):
785	'''
786	Save phase and histogram information into "pseudo-gpx" files. The phase
787	information is overwritten each time this is called, but histogram information is
788	appended after each sequential step.
789
790	:param str GPXfile: full .gpx file name
791	:param dict Histograms: dictionary of histograms as {name:data,...}
792	:param dict Phases: dictionary of phases that use histograms
793	:param dict RigidBodies: dictionary of rigid bodies
794	:param dict CovData: dictionary of refined variables, varyList, & covariance matrix
795	:param dict parmFrozen: dict with frozen parameters for all phases
796	and histograms (specified as str values)
797	'''
798
799	import distutils.file_util as dfu
800
801	GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
802	fp = open(GPXphase,'rb')
803	data = cPickleLoad(fp) # first block in file should be Phases
804	if data[0][0] != 'Phases':
805	raise Exception('Unexpected block in {} file. How did this happen?'
806	.format(GPXphase))
807	fp.close()
808	# update previous phase info
809	for datum in data[1:]:
810	if datum[0] in Phases:
811	datum[1].update(Phases[datum[0]])
812	# save latest Phase/refinement info
813	fp = open(GPXphase,'wb')
814	cPickle.dump(data,fp,1)
815	cPickle.dump([['Covariance',CovData]],fp,1)
816	cPickle.dump([['Rigid bodies',RigidBodies]],fp,1)
817	cPickle.dump([['parmFrozen',parmFrozen]],fp,1)
818	fp.close()
819	# create an entry that looks like a PWDR tree item
820	for key in Histograms:
821	if key.startswith('PWDR '):
822	break
823	else:
824	raise Exception('No PWDR entry in Histogram dict!')
825	histname = key
826	hist = copy.deepcopy(Histograms[key])
827	xfer_dict = {'Index Peak List': [[], []],
828	'Comments': [],
829	'Unit Cells List': [],
830	'Peak List': {'peaks': [], 'sigDict': {}},
831	}
832	histData = hist['Data']
833	del hist['Data']
834	for key in ('Limits','Background','Instrument Parameters',
835	'Sample Parameters','Reflection Lists'):
836	xfer_dict[key] = hist[key]
837	if key == 'Background': # remove fixed background from file
838	xfer_dict['Background'][1] = {k:hist['Background'][1][k]
839	for k in hist['Background'][1]
840	if not k.startswith('_fixed')}
841	del hist[key]
842	# xform into a gpx-type entry
843	data = []
844	data.append([histname,[hist,histData,histname]])
845	for key in ['Comments','Limits','Background','Instrument Parameters',
846	'Sample Parameters','Peak List','Index Peak List',
847	'Unit Cells List','Reflection Lists']:
848	data.append([key,xfer_dict[key]])
849	# append histogram to histogram info
850	GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
851	fp = open(GPXhist,'ab')
852	cPickle.dump(data,fp,1)
853	fp.close()
854	return
855
856	def SetSeqResult(GPXfile,Histograms,SeqResult):
857	'''
858	Places the sequential results information into a GPX file
859	after a refinement has been completed.
860	Called at the end of :meth:`GSASIIstrMain.SeqRefine`
861
862	:param str GPXfile: full .gpx file name
863	'''
864	GPXback = GPXBackup(GPXfile)
865	G2fil.G2Print ('Read from file:'+GPXback)
866	G2fil.G2Print ('Save to file :'+GPXfile)
867	GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
868	fp = open(GPXphase,'rb')
869	data = cPickleLoad(fp) # first block in file should be Phases
870	if data[0][0] != 'Phases':
871	raise Exception('Unexpected block in {} file. How did this happen?'.format(GPXphase))
872	Phases = {}
873	for name,vals in data[1:]:
874	Phases[name] = vals
875	name,CovData = cPickleLoad(fp)[0] # 2nd block in file should be Covariance
876	name,RigidBodies = cPickleLoad(fp)[0] # 3rd block in file should be Rigid Bodies
877	name,parmFrozenDict = cPickleLoad(fp)[0] # 4th block in file should be frozen parameters
878	fp.close()
879	GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
880	hist = open(GPXhist,'rb')
881	# build an index to the GPXhist file
882	histIndex = {}
883	while True:
884	loc = hist.tell()
885	try:
886	datum = cPickleLoad(hist)[0]
887	except EOFError:
888	break
889	histIndex[datum[0]] = loc
890
891	infile = open(GPXback,'rb')
892	outfile = open(GPXfile,'wb')
893	while True:
894	try:
895	data = cPickleLoad(infile)
896	except EOFError:
897	break
898	datum = data[0]
899	if datum[0] == 'Sequential results':
900	data[0][1] = SeqResult
901	elif datum[0] == 'Phases':
902	for pdata in data[1:]:
903	if pdata[0] in Phases:
904	pdata[1].update(Phases[pdata[0]])
905	elif datum[0] == 'Covariance':
906	data[0][1] = CovData
907	elif datum[0] == 'Rigid bodies':
908	data[0][1] = RigidBodies
909	elif datum[0] == 'Controls': # reset the Copy Next flag after a sequential fit
910	Controls = data[0][1]
911	Controls['Copy2Next'] = False
912	for key in parmFrozenDict:
913	Controls['parmFrozen'][key] = [
914	i if type(i) is G2obj.G2VarObj
915	else G2obj.G2VarObj(i)
916	for i in parmFrozenDict[key]]
917	elif datum[0] in histIndex:
918	hist.seek(histIndex[datum[0]])
919	hdata = cPickleLoad(hist)
920	if data[0][0] != hdata[0][0]:
921	G2fil.G2Print('Error! Updating {} with {}'.format(data[0][0],hdata[0][0]))
922	data[0] = hdata[0]
923	xferItems = ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']
924	hItems = {name:j+1 for j,(name,val) in enumerate(hdata[1:]) if name in xferItems}
925	for j,(name,val) in enumerate(data[1:]):
926	if name not in xferItems: continue
927	data[j+1][1] = hdata[hItems[name]][1]
928	cPickle.dump(data,outfile,1)
929	hist.close()
930	infile.close()
931	outfile.close()
932	# clean up tmp files
933	try:
934	os.remove(GPXphase)
935	except:
936	G2fil.G2Print('Warning: unable to delete {}'.format(GPXphase))
937	try:
938	os.remove(GPXhist)
939	except:
940	G2fil.G2Print('Warning: unable to delete {}'.format(GPXhist))
941	G2fil.G2Print ('GPX file merge completed')
942
943	#==============================================================================
944	# Refinement routines
945	#==============================================================================
946	def ShowBanner(pFile=None):
947	'Print authorship, copyright and citation notice'
948	pFile.write(80''+'\n')
949	pFile.write(' General Structure Analysis System-II Crystal Structure Refinement\n')
950	pFile.write(' by Robert B. Von Dreele & Brian H. Toby\n')
951	pFile.write(' Argonne National Laboratory(C), 2010\n')
952	pFile.write(' This product includes software developed by the UChicago Argonne, LLC,\n')
953	pFile.write(' as Operator of Argonne National Laboratory.\n')
954	pFile.write(' Please cite:\n')
955	pFile.write(' B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)\n')
956
957	pFile.write(80''+'\n')
958
959	def ShowControls(Controls,pFile=None,SeqRef=False,preFrozenCount=0):
960	'Print controls information'
961	pFile.write(' Least squares controls:\n')
962	pFile.write(' Refinement type: %s\n'%Controls['deriv type'])
963	if 'Hessian' in Controls['deriv type']:
964	pFile.write(' Maximum number of cycles: %d\n'%Controls['max cyc'])
965	else:
966	pFile.write(' Minimum delta-M/M for convergence: %.2g\n'%(Controls['min dM/M']))
967	pFile.write(' Regularize hydrogens (if any): %s\n'%Controls.get('HatomFix',False))
968	pFile.write(' Initial shift factor: %.3f\n'%(Controls['shift factor']))
969	if SeqRef:
970	pFile.write(' Sequential refinement controls:\n')
971	pFile.write(' Copy of histogram results to next: %s\n'%(Controls['Copy2Next']))
972	pFile.write(' Process histograms in reverse order: %s\n'%(Controls['Reverse Seq']))
973	if preFrozenCount:
974	pFile.write('\n Starting refinement with {} Frozen variables\n\n'.format(preFrozenCount))
975
976	def GetPawleyConstr(SGLaue,PawleyRef,im,pawleyVary):
977	'needs a doc string'
978	# if SGLaue in ['-1','2/m','mmm']:
979	# return #no Pawley symmetry required constraints
980	eqvDict = {}
981	for i,varyI in enumerate(pawleyVary):
982	eqvDict[varyI] = []
983	refI = int(varyI.split(':')[-1])
984	ih,ik,il = PawleyRef[refI][:3]
985	dspI = PawleyRef[refI][4+im]
986	for varyJ in pawleyVary[i+1:]:
987	refJ = int(varyJ.split(':')[-1])
988	jh,jk,jl = PawleyRef[refJ][:3]
989	dspJ = PawleyRef[refJ][4+im]
990	if SGLaue in ['4/m','4/mmm']:
991	isum = ih2+ik2
992	jsum = jh2+jk2
993	if abs(il) == abs(jl) and isum == jsum:
994	eqvDict[varyI].append(varyJ)
995	elif SGLaue in ['3R','3mR']:
996	isum = ih2+ik2+il**2
997	jsum = jh2+jk2+jl**2
998	isum2 = ihik+ihil+ik*il
999	jsum2 = jhjk+jhjl+jk*jl
1000	if isum == jsum and isum2 == jsum2:
1001	eqvDict[varyI].append(varyJ)
1002	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1003	isum = ih2+ik2+ih*ik
1004	jsum = jh2+jk2+jh*jk
1005	if abs(il) == abs(jl) and isum == jsum:
1006	eqvDict[varyI].append(varyJ)
1007	elif SGLaue in ['m3','m3m']:
1008	isum = ih2+ik2+il**2
1009	jsum = jh2+jk2+jl**2
1010	if isum == jsum:
1011	eqvDict[varyI].append(varyJ)
1012	elif abs(dspI-dspJ)/dspI < 1.e-4:
1013	eqvDict[varyI].append(varyJ)
1014	for item in pawleyVary:
1015	if eqvDict[item]:
1016	for item2 in pawleyVary:
1017	if item2 in eqvDict[item]:
1018	eqvDict[item2] = []
1019	G2mv.StoreEquivalence(item,eqvDict[item])
1020
1021	def cellVary(pfx,SGData):
1022	'''Creates equivalences for a phase based on the Laue class.
1023	Returns a list of A tensor terms that are non-zero.
1024	'''
1025	if SGData['SGLaue'] in ['-1',]:
1026	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1027	elif SGData['SGLaue'] in ['2/m',]:
1028	if SGData['SGUniq'] == 'a':
1029	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
1030	elif SGData['SGUniq'] == 'b':
1031	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
1032	else:
1033	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
1034	elif SGData['SGLaue'] in ['mmm',]:
1035	return [pfx+'A0',pfx+'A1',pfx+'A2']
1036	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1037	G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',))
1038	return [pfx+'A0',pfx+'A1',pfx+'A2']
1039	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1040	G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A3',))
1041	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
1042	elif SGData['SGLaue'] in ['3R', '3mR']:
1043	G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A2',))
1044	G2mv.StoreEquivalence(pfx+'A3',(pfx+'A4',pfx+'A5',))
1045	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1046	elif SGData['SGLaue'] in ['m3m','m3']:
1047	G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A2',))
1048	return [pfx+'A0',pfx+'A1',pfx+'A2']
1049
1050	def modVary(pfx,SSGData):
1051	vary = []
1052	for i,item in enumerate(SSGData['modSymb']):
1053	if item in ['a','b','g']:
1054	vary.append(pfx+'mV%d'%(i))
1055	return vary
1056
1057	################################################################################
1058	##### Rigid Body Models and not General.get('doPawley')
1059	################################################################################
1060
1061	def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
1062	'Get Rigid body info from tree entry and print it to .LST file'
1063
1064	def PrintResRBModel(RBModel):
1065	pFile.write('Residue RB name: %s no.atoms: %d, No. times used: %d\n'%
1066	(RBModel['RBname'],len(RBModel['rbTypes']),RBModel['useCount']))
1067	for i in WriteResRBModel(RBModel):
1068	pFile.write(i)
1069
1070	def PrintVecRBModel(RBModel):
1071	pFile.write('Vector RB name: %s no.atoms: %d No. times used: %d\n'%
1072	(RBModel['RBname'],len(RBModel['rbTypes']),RBModel['useCount']))
1073	for i in WriteVecRBModel(RBModel):
1074	pFile.write(i)
1075	pFile.write('Orientation defined by: atom %s -> atom %s & atom %s -> atom %s\n'%
1076	(RBModel['rbRef'][0],RBModel['rbRef'][1],RBModel['rbRef'][0],RBModel['rbRef'][2]))
1077
1078	rbVary = []
1079	rbDict = {}
1080	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
1081	if len(rbIds['Vector']):
1082	for irb,item in enumerate(rbIds['Vector']):
1083	if rigidbodyDict['Vector'][item]['useCount']:
1084	RBmags = rigidbodyDict['Vector'][item]['VectMag']
1085	RBrefs = rigidbodyDict['Vector'][item]['VectRef']
1086	for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
1087	pid = '::RBV;'+str(i)+':'+str(irb)
1088	rbDict[pid] = mag
1089	if ref:
1090	rbVary.append(pid)
1091	if Print:
1092	pFile.write('\nVector rigid body model:\n')
1093	PrintVecRBModel(rigidbodyDict['Vector'][item])
1094	if len(rbIds['Residue']):
1095	for item in rbIds['Residue']:
1096	if rigidbodyDict['Residue'][item]['useCount']:
1097	if Print:
1098	pFile.write('\nResidue rigid body model:\n')
1099	PrintResRBModel(rigidbodyDict['Residue'][item])
1100	return rbVary,rbDict
1101
1102	def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
1103	'needs a doc string'
1104
1105	def PrintRBVectandSig(VectRB,VectSig):
1106	pFile.write('\n Rigid body vector magnitudes for %s:\n'%VectRB['RBname'])
1107	namstr = ' names :'
1108	valstr = ' values:'
1109	sigstr = ' esds :'
1110	for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
1111	namstr += '%12s'%('Vect '+str(i))
1112	valstr += '%12.4f'%(val)
1113	if sig:
1114	sigstr += '%12.4f'%(sig)
1115	else:
1116	sigstr += 12*' '
1117	pFile.write(namstr+'\n')
1118	pFile.write(valstr+'\n')
1119	pFile.write(sigstr+'\n')
1120
1121	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
1122	if not RBIds['Vector']:
1123	return
1124	for irb,item in enumerate(RBIds['Vector']):
1125	if rigidbodyDict['Vector'][item]['useCount']:
1126	VectSig = []
1127	RBmags = rigidbodyDict['Vector'][item]['VectMag']
1128	for i,mag in enumerate(RBmags):
1129	name = '::RBV;'+str(i)+':'+str(irb)
1130	if name in sigDict:
1131	VectSig.append(sigDict[name])
1132	PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)
1133
1134	################################################################################
1135	##### Phase data
1136	################################################################################
1137	def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None,seqRef=False):
1138	'''Setup the phase information for a structural refinement, used for
1139	regular and sequential refinements, optionally printing information
1140	to the .lst file (if Print is True)
1141	'''
1142
1143	def PrintFFtable(FFtable):
1144	pFile.write('\n X-ray scattering factors:\n')
1145	pFile.write(' Symbol fa fb fc\n')
1146	pFile.write(99*'-'+'\n')
1147	for Ename in FFtable:
1148	ffdata = FFtable[Ename]
1149	fa = ffdata['fa']
1150	fb = ffdata['fb']
1151	pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1152	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc']))
1153
1154	def PrintMFtable(MFtable):
1155	pFile.write('\n <j0> Magnetic scattering factors:\n')
1156	pFile.write(' Symbol mfa mfb mfc\n')
1157	pFile.write(99*'-'+'\n')
1158	for Ename in MFtable:
1159	mfdata = MFtable[Ename]
1160	fa = mfdata['mfa']
1161	fb = mfdata['mfb']
1162	pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1163	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['mfc']))
1164	pFile.write('\n <j2> Magnetic scattering factors:\n')
1165	pFile.write(' Symbol nfa nfb nfc\n')
1166	pFile.write(99*'-'+'\n')
1167	for Ename in MFtable:
1168	mfdata = MFtable[Ename]
1169	fa = mfdata['nfa']
1170	fb = mfdata['nfb']
1171	pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1172	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['nfc']))
1173
1174	def PrintBLtable(BLtable):
1175	pFile.write('\n Neutron scattering factors:\n')
1176	pFile.write(' Symbol isotope mass b resonant terms\n')
1177	pFile.write(99*'-'+'\n')
1178	for Ename in BLtable:
1179	bldata = BLtable[Ename]
1180	isotope = bldata[0]
1181	mass = bldata[1]['Mass']
1182	if 'BW-LS' in bldata[1]:
1183	bres = bldata[1]['BW-LS']
1184	blen = 0
1185	else:
1186	blen = bldata[1]['SL'][0]
1187	bres = []
1188	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
1189	for item in bres:
1190	line += '%10.5g'%(item)
1191	pFile.write(line+'\n')
1192
1193	def PrintRBObjects(resRBData,vecRBData):
1194
1195	def PrintRBThermals():
1196	tlstr = ['11','22','33','12','13','23']
1197	sstr = ['12','13','21','23','31','32','AA','BB']
1198	TLS = RB['ThermalMotion'][1]
1199	TLSvar = RB['ThermalMotion'][2]
1200	if 'T' in RB['ThermalMotion'][0]:
1201	pFile.write('TLS data\n')
1202	text = ''
1203	for i in range(6):
1204	text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
1205	pFile.write(text+'\n')
1206	if 'L' in RB['ThermalMotion'][0]:
1207	text = ''
1208	for i in range(6,12):
1209	text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
1210	pFile.write(text+'\n')
1211	if 'S' in RB['ThermalMotion'][0]:
1212	text = ''
1213	for i in range(12,20):
1214	text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
1215	pFile.write(text+'\n')
1216	if 'U' in RB['ThermalMotion'][0]:
1217	pFile.write('Uiso data\n')
1218	text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])
1219	pFile.write(text+'\n')
1220
1221	if len(resRBData):
1222	for RB in resRBData:
1223	Oxyz = RB['Orig'][0]
1224	Qrijk = RB['Orient'][0]
1225	Angle = 2.0*acosd(Qrijk[0])
1226	pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'%
1227	(RB['RBname'],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1]))
1228	pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'%
1229	(Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1]))
1230	Torsions = RB['Torsions']
1231	if len(Torsions):
1232	text = 'Torsions: '
1233	for torsion in Torsions:
1234	text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
1235	pFile.write(text+'\n')
1236	PrintRBThermals()
1237	if len(vecRBData):
1238	for RB in vecRBData:
1239	Oxyz = RB['Orig'][0]
1240	Qrijk = RB['Orient'][0]
1241	Angle = 2.0*acosd(Qrijk[0])
1242	pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'%
1243	(RB['RBname'],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1]))
1244	pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'%
1245	(Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1]))
1246	PrintRBThermals()
1247
1248	def PrintAtoms(General,Atoms):
1249	cx,ct,cs,cia = General['AtomPtrs']
1250	pFile.write('\n Atoms:\n')
1251	line = ' name type refine? x y z '+ \
1252	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
1253	if General['Type'] == 'macromolecular':
1254	line = ' res no residue chain'+line
1255	pFile.write(line+'\n')
1256	if General['Type'] in ['nuclear','magnetic','faulted',]:
1257	pFile.write(135*'-'+'\n')
1258	for i,at in enumerate(Atoms):
1259	line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
1260	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
1261	if at[cia] == 'I':
1262	line += '%8.5f'%(at[cia+1])+48*' '
1263	else:
1264	line += 8*' '
1265	for j in range(6):
1266	line += '%8.5f'%(at[cia+2+j])
1267	pFile.write(line+'\n')
1268	elif General['Type'] == 'macromolecular':
1269	pFile.write(135*'-'+'\n')
1270	for i,at in enumerate(Atoms):
1271	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])+ \
1272	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
1273	if at[cia] == 'I':
1274	line += '%8.4f'%(at[cia+1])+48*' '
1275	else:
1276	line += 8*' '
1277	for j in range(6):
1278	line += '%8.4f'%(at[cia+2+j])
1279	pFile.write(line+'\n')
1280
1281	def PrintMoments(General,Atoms):
1282	cx,ct,cs,cia = General['AtomPtrs']
1283	cmx = cx+4
1284	AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1285	pFile.write('\n Magnetic moments:\n')
1286	line = ' name type refine? Mx My Mz '
1287	pFile.write(line+'\n')
1288	pFile.write(135*'-'+'\n')
1289	for i,at in enumerate(Atoms):
1290	if AtInfo[at[ct]]:
1291	line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cmx])+'%10.5f'%(at[cmx+1])+ \
1292	'%10.5f'%(at[cmx+2])
1293	pFile.write(line+'\n')
1294
1295
1296	def PrintWaves(General,Atoms):
1297	cx,ct,cs,cia = General['AtomPtrs']
1298	pFile.write('\n Modulation waves\n')
1299	names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1300	'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1301	'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1302	pFile.write(135*'-'+'\n')
1303	for i,at in enumerate(Atoms):
1304	AtomSS = at[-1]['SS1']
1305	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1306	Waves = AtomSS[Stype]
1307	if len(Waves):
1308	pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
1309	(at[ct-1],at[cs],Stype,Waves[0]))
1310	else:
1311	continue
1312	for iw,wave in enumerate(Waves[1:]):
1313	line = ''
1314	if Waves[0] in ['Block','ZigZag'] and Stype == 'Spos' and not iw:
1315	for item in names[Stype][6:]:
1316	line += '%8s '%(item)
1317	else:
1318	if Stype == 'Spos':
1319	for item in names[Stype][:6]:
1320	line += '%8s '%(item)
1321	else:
1322	for item in names[Stype]:
1323	line += '%8s '%(item)
1324	pFile.write(line+'\n')
1325	line = ''
1326	for item in wave[0]:
1327	line += '%8.4f '%(item)
1328	line += ' Refine? '+str(wave[1])
1329	pFile.write(line+'\n')
1330
1331	def PrintTexture(textureData):
1332	topstr = '\n Spherical harmonics texture: Order:' + \
1333	str(textureData['Order'])
1334	if textureData['Order']:
1335	pFile.write('%s Refine? %s\n'%(topstr,textureData['SH Coeff'][0]))
1336	else:
1337	pFile.write(topstr+'\n')
1338	return
1339	names = ['omega','chi','phi']
1340	line = '\n'
1341	for name in names:
1342	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
1343	pFile.write(line+'\n')
1344	pFile.write('\n Texture coefficients:\n')
1345	SHcoeff = textureData['SH Coeff'][1]
1346	SHkeys = list(SHcoeff.keys())
1347	nCoeff = len(SHcoeff)
1348	nBlock = nCoeff//10+1
1349	iBeg = 0
1350	iFin = min(iBeg+10,nCoeff)
1351	for block in range(nBlock):
1352	ptlbls = ' names :'
1353	ptstr = ' values:'
1354	for item in SHkeys[iBeg:iFin]:
1355	ptlbls += '%12s'%(item)
1356	ptstr += '%12.4f'%(SHcoeff[item])
1357	pFile.write(ptlbls+'\n')
1358	pFile.write(ptstr+'\n')
1359	iBeg += 10
1360	iFin = min(iBeg+10,nCoeff)
1361
1362	def MakeRBParms(rbKey,phaseVary,phaseDict):
1363	rbid = str(rbids.index(RB['RBId']))
1364	pfxRB = pfx+'RB'+rbKey+'P'
1365	pstr = ['x','y','z']
1366	ostr = ['a','i','j','k']
1367	for i in range(3):
1368	name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1369	phaseDict[name] = RB['Orig'][0][i]
1370	if RB['Orig'][1]:
1371	phaseVary += [name,]
1372	pfxRB = pfx+'RB'+rbKey+'O'
1373	A,V = G2mth.Q2AV(RB['Orig'][0])
1374	fixAxis = [0, np.abs(V).argmax()+1]
1375	for i in range(4):
1376	name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1377	phaseDict[name] = RB['Orient'][0][i]
1378	if RB['Orient'][1] == 'AV' and i:
1379	phaseVary += [name,]
1380	elif RB['Orient'][1] == 'A' and not i:
1381	phaseVary += [name,]
1382	elif RB['Orient'][1] == 'V' and i not in fixAxis:
1383	phaseVary += [name,]
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	for i in WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
2092	pFile.write(i+'\n')
2093
2094	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
2095	for i in WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
2096	pFile.write(i)
2097
2098	def PrintRBObjTorAndSig(rbsx):
2099	nTors = len(RBObj['Torsions'])
2100	if nTors:
2101	for i in WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
2102	pFile.write(i)
2103
2104	def PrintSHtextureAndSig(textureData,SHtextureSig):
2105	Tindx = 1.0
2106	Tvar = 0.0
2107	pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
2108	names = ['omega','chi','phi']
2109	namstr = ' names :'
2110	ptstr = ' values:'
2111	sigstr = ' esds :'
2112	for name in names:
2113	namstr += '%12s'%(name)
2114	ptstr += '%12.3f'%(textureData['Sample '+name][1])
2115	if 'Sample '+name in SHtextureSig:
2116	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2117	else:
2118	sigstr += 12*' '
2119	pFile.write(namstr+'\n')
2120	pFile.write(ptstr+'\n')
2121	pFile.write(sigstr+'\n')
2122	pFile.write('\n Texture coefficients:\n')
2123	SHcoeff = textureData['SH Coeff'][1]
2124	SHkeys = list(SHcoeff.keys())
2125	nCoeff = len(SHcoeff)
2126	nBlock = nCoeff//10+1
2127	iBeg = 0
2128	iFin = min(iBeg+10,nCoeff)
2129	for block in range(nBlock):
2130	namstr = ' names :'
2131	ptstr = ' values:'
2132	sigstr = ' esds :'
2133	for name in SHkeys[iBeg:iFin]:
2134	namstr += '%12s'%(name)
2135	ptstr += '%12.3f'%(SHcoeff[name])
2136	l = 2.0*eval(name.strip('C'))[0]+1
2137	Tindx += SHcoeff[name]**2/l
2138	if name in SHtextureSig:
2139	Tvar += (2.SHcoeff[name]SHtextureSig[name]/l)**2
2140	sigstr += '%12.3f'%(SHtextureSig[name])
2141	else:
2142	sigstr += 12*' '
2143	pFile.write(namstr+'\n')
2144	pFile.write(ptstr+'\n')
2145	pFile.write(sigstr+'\n')
2146	iBeg += 10
2147	iFin = min(iBeg+10,nCoeff)
2148	pFile.write(' Texture index J = %.3f(%d)'%(Tindx,int(1000*np.sqrt(Tvar))))
2149
2150	##########################################################################
2151	# SetPhaseData starts here
2152	pFile.write('\n Phases:\n')
2153	for phase in Phases:
2154	pFile.write(' Result for phase: %s\n'%phase)
2155	pFile.write(135*'='+'\n')
2156	Phase = Phases[phase]
2157	General = Phase['General']
2158	SGData = General['SGData']
2159	Atoms = Phase['Atoms']
2160	AtLookup = []
2161	if Atoms and not General.get('doPawley'):
2162	cx,ct,cs,cia = General['AtomPtrs']
2163	AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
2164	cell = General['Cell']
2165	pId = Phase['pId']
2166	pfx = str(pId)+'::'
2167	if cell[0]:
2168	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
2169	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
2170	pFile.write(' Reciprocal metric tensor: \n')
2171	ptfmt = "%15.9f"
2172	names = ['A11','A22','A33','A12','A13','A23']
2173	namstr = ' names :'
2174	ptstr = ' values:'
2175	sigstr = ' esds :'
2176	for name,a,siga in zip(names,A,sigA):
2177	namstr += '%15s'%(name)
2178	ptstr += ptfmt%(a)
2179	if siga:
2180	sigstr += ptfmt%(siga)
2181	else:
2182	sigstr += 15*' '
2183	pFile.write(namstr+'\n')
2184	pFile.write(ptstr+'\n')
2185	pFile.write(sigstr+'\n')
2186	cell[1:7] = G2lat.A2cell(A)
2187	cell[7] = G2lat.calc_V(A)
2188	pFile.write(' New unit cell:\n')
2189	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
2190	names = ['a','b','c','alpha','beta','gamma','Volume']
2191	namstr = ' names :'
2192	ptstr = ' values:'
2193	sigstr = ' esds :'
2194	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
2195	namstr += '%12s'%(name)
2196	ptstr += fmt%(a)
2197	if siga:
2198	sigstr += fmt%(siga)
2199	else:
2200	sigstr += 12*' '
2201	pFile.write(namstr+'\n')
2202	pFile.write(ptstr+'\n')
2203	pFile.write(sigstr+'\n')
2204	ik = 6 #for Pawley stuff below
2205	if General.get('Modulated',False):
2206	ik = 7
2207	Vec,vRef,maxH = General['SuperVec']
2208	if vRef:
2209	pFile.write(' New modulation vector:\n')
2210	namstr = ' names :'
2211	ptstr = ' values:'
2212	sigstr = ' esds :'
2213	for iv,var in enumerate(['mV0','mV1','mV2']):
2214	namstr += '%12s'%(pfx+var)
2215	ptstr += '%12.6f'%(parmDict[pfx+var])
2216	if pfx+var in sigDict:
2217	Vec[iv] = parmDict[pfx+var]
2218	sigstr += '%12.6f'%(sigDict[pfx+var])
2219	else:
2220	sigstr += 12*' '
2221	pFile.write(namstr+'\n')
2222	pFile.write(ptstr+'\n')
2223	pFile.write(sigstr+'\n')
2224
2225	General['Mass'] = 0.
2226	if Phase['General'].get('doPawley'):
2227	pawleyRef = Phase['Pawley ref']
2228	for i,refl in enumerate(pawleyRef):
2229	key = pfx+'PWLref:'+str(i)
2230	refl[ik] = parmDict[key]
2231	if key in sigDict:
2232	refl[ik+1] = sigDict[key]
2233	else:
2234	refl[ik+1] = 0
2235	else:
2236	VRBIds = RBIds['Vector']
2237	RRBIds = RBIds['Residue']
2238	RBModels = Phase['RBModels']
2239	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
2240	jrb = VRBIds.index(RBObj['RBId'])
2241	rbsx = str(irb)+':'+str(jrb)
2242	pFile.write(' Vector rigid body parameters:\n')
2243	PrintRBObjPOAndSig('RBV',rbsx)
2244	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
2245	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2246	jrb = RRBIds.index(RBObj['RBId'])
2247	rbsx = str(irb)+':'+str(jrb)
2248	pFile.write(' Residue rigid body parameters:\n')
2249	PrintRBObjPOAndSig('RBR',rbsx)
2250	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2251	PrintRBObjTorAndSig(rbsx)
2252	atomsSig = {}
2253	wavesSig = {}
2254	cx,ct,cs,cia = General['AtomPtrs']
2255	for i,at in enumerate(Atoms):
2256	names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2257	cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2258	cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2259	cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2260	for ind in range(cx,cx+4):
2261	at[ind] = parmDict[names[ind]]
2262	if ind in range(cx,cx+3):
2263	name = names[ind].replace('A','dA')
2264	else:
2265	name = names[ind]
2266	if name in sigDict:
2267	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2268	if at[cia] == 'I':
2269	at[cia+1] = parmDict[names[cia+1]]
2270	if names[cia+1] in sigDict:
2271	atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2272	else:
2273	for ind in range(cia+2,cia+8):
2274	at[ind] = parmDict[names[ind]]
2275	if names[ind] in sigDict:
2276	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2277	if General['Type'] == 'magnetic':
2278	for ind in range(cx+4,cx+7):
2279	at[ind] = parmDict[names[ind]]
2280	if names[ind] in sigDict:
2281	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2282	ind = General['AtomTypes'].index(at[ct])
2283	General['Mass'] += General['AtomMass'][ind]at[cx+3]at[cx+5]
2284	if General.get('Modulated',False):
2285	AtomSS = at[-1]['SS1']
2286	for Stype in ['Sfrac','Spos','Sadp','Smag']:
2287	Waves = AtomSS[Stype]
2288	if len(Waves):
2289	waveType = Waves[0]
2290	else:
2291	continue
2292	for iw,wave in enumerate(Waves[1:]):
2293	stiw = str(i)+':'+str(iw)
2294	if Stype == 'Spos':
2295	if waveType in ['ZigZag','Block',] and not iw:
2296	names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2297	else:
2298	names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2299	'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2300	elif Stype == 'Sadp':
2301	names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2302	'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2303	'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2304	'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2305	elif Stype == 'Sfrac':
2306	if 'Crenel' in waveType and not iw:
2307	names = ['Fzero:'+stiw,'Fwid:'+stiw]
2308	else:
2309	names = ['Fsin:'+stiw,'Fcos:'+stiw]
2310	elif Stype == 'Smag':
2311	names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2312	'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2313	for iname,name in enumerate(names):
2314	AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2315	if pfx+name in sigDict:
2316	wavesSig[name] = sigDict[pfx+name]
2317
2318	PrintAtomsAndSig(General,Atoms,atomsSig)
2319	if General['Type'] == 'magnetic':
2320	PrintMomentsAndSig(General,Atoms,atomsSig)
2321	if General.get('Modulated',False):
2322	PrintWavesAndSig(General,Atoms,wavesSig)
2323
2324	density = G2mth.getDensity(General)[0]
2325	pFile.write('\n Density: {:.4f} g/cm**3\n'.format(density))
2326
2327
2328	textureData = General['SH Texture']
2329	if textureData['Order']:
2330	SHtextureSig = {}
2331	for name in ['omega','chi','phi']:
2332	aname = pfx+'SH '+name
2333	textureData['Sample '+name][1] = parmDict[aname]
2334	if aname in sigDict:
2335	SHtextureSig['Sample '+name] = sigDict[aname]
2336	for name in textureData['SH Coeff'][1]:
2337	aname = pfx+name
2338	textureData['SH Coeff'][1][name] = parmDict[aname]
2339	if aname in sigDict:
2340	SHtextureSig[name] = sigDict[aname]
2341	PrintSHtextureAndSig(textureData,SHtextureSig)
2342	if phase in RestraintDict and not Phase['General'].get('doPawley'):
2343	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2344	textureData,RestraintDict[phase],pFile)
2345
2346	################################################################################
2347	##### Histogram & Phase data
2348	################################################################################
2349
2350	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2351	'''Loads the HAP histogram/phase information into dicts
2352
2353	:param dict Phases: phase information
2354	:param dict Histograms: Histogram information
2355	:param bool Print: prints information as it is read
2356	:param file pFile: file object to print to (the default, None causes printing to the console)
2357	:param bool resetRefList: Should the contents of the Reflection List be initialized
2358	on loading. The default, True, initializes the Reflection List as it is loaded.
2359
2360	:returns: (hapVary,hapDict,controlDict)
2361	* hapVary: list of refined variables
2362	* hapDict: dict with refined variables and their values
2363	* controlDict: dict with fixed parameters
2364	'''
2365
2366	def PrintSize(hapData):
2367	if hapData[0] in ['isotropic','uniaxial']:
2368	line = '\n Size model : %9s'%(hapData[0])
2369	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2370	if hapData[0] == 'uniaxial':
2371	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2372	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2373	pFile.write(line+'\n')
2374	else:
2375	pFile.write('\n Size model : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2376	(hapData[0],hapData[1][2],hapData[2][2]))
2377	Snames = ['S11','S22','S33','S12','S13','S23']
2378	ptlbls = ' names :'
2379	ptstr = ' values:'
2380	varstr = ' refine:'
2381	for i,name in enumerate(Snames):
2382	ptlbls += '%12s' % (name)
2383	ptstr += '%12.3f' % (hapData[4][i])
2384	varstr += '%12s' % (str(hapData[5][i]))
2385	pFile.write(ptlbls+'\n')
2386	pFile.write(ptstr+'\n')
2387	pFile.write(varstr+'\n')
2388
2389	def PrintMuStrain(hapData,SGData):
2390	if hapData[0] in ['isotropic','uniaxial']:
2391	line = '\n Mustrain model: %9s'%(hapData[0])
2392	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2393	if hapData[0] == 'uniaxial':
2394	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2395	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2396	pFile.write(line+'\n')
2397	else:
2398	pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2399	(hapData[0],hapData[1][2],hapData[2][2]))
2400	Snames = G2spc.MustrainNames(SGData)
2401	ptlbls = ' names :'
2402	ptstr = ' values:'
2403	varstr = ' refine:'
2404	for i,name in enumerate(Snames):
2405	ptlbls += '%12s' % (name)
2406	ptstr += '%12.1f' % (hapData[4][i])
2407	varstr += '%12s' % (str(hapData[5][i]))
2408	pFile.write(ptlbls+'\n')
2409	pFile.write(ptstr+'\n')
2410	pFile.write(varstr+'\n')
2411
2412	def PrintHStrain(hapData,SGData):
2413	pFile.write('\n Hydrostatic/elastic strain:\n')
2414	Hsnames = G2spc.HStrainNames(SGData)
2415	ptlbls = ' names :'
2416	ptstr = ' values:'
2417	varstr = ' refine:'
2418	for i,name in enumerate(Hsnames):
2419	ptlbls += '%12s' % (name)
2420	ptstr += '%12.4g' % (hapData[0][i])
2421	varstr += '%12s' % (str(hapData[1][i]))
2422	pFile.write(ptlbls+'\n')
2423	pFile.write(ptstr+'\n')
2424	pFile.write(varstr+'\n')
2425
2426	def PrintSHPO(hapData):
2427	pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2428	ptlbls = ' names :'
2429	ptstr = ' values:'
2430	for item in hapData[5]:
2431	ptlbls += '%12s'%(item)
2432	ptstr += '%12.3f'%(hapData[5][item])
2433	pFile.write(ptlbls+'\n')
2434	pFile.write(ptstr+'\n')
2435
2436	def PrintBabinet(hapData):
2437	pFile.write('\n Babinet form factor modification:\n')
2438	ptlbls = ' names :'
2439	ptstr = ' values:'
2440	varstr = ' refine:'
2441	for name in ['BabA','BabU']:
2442	ptlbls += '%12s' % (name)
2443	ptstr += '%12.6f' % (hapData[name][0])
2444	varstr += '%12s' % (str(hapData[name][1]))
2445	pFile.write(ptlbls+'\n')
2446	pFile.write(ptstr+'\n')
2447	pFile.write(varstr+'\n')
2448
2449	hapDict = {}
2450	hapVary = []
2451	controlDict = {}
2452
2453	for phase in Phases:
2454	HistoPhase = Phases[phase]['Histograms']
2455	SGData = Phases[phase]['General']['SGData']
2456	cell = Phases[phase]['General']['Cell'][1:7]
2457	A = G2lat.cell2A(cell)
2458	if Phases[phase]['General'].get('Modulated',False):
2459	SSGData = Phases[phase]['General']['SSGData']
2460	Vec,x,maxH = Phases[phase]['General']['SuperVec']
2461	pId = Phases[phase]['pId']
2462	for histogram in Histograms:
2463	if histogram not in HistoPhase and phase in Histograms[histogram]['Reflection Lists']:
2464	#remove previously created reflection list if histogram is removed from phase
2465	#print("removing ",phase,"from",histogram)
2466	del Histograms[histogram]['Reflection Lists'][phase]
2467	histoList = list(HistoPhase.keys())
2468	histoList.sort()
2469	for histogram in histoList:
2470	try:
2471	Histogram = Histograms[histogram]
2472	except KeyError:
2473	#skip if histogram not included e.g. in a sequential refinement
2474	continue
2475	if not HistoPhase[histogram]['Use']: #remove previously created & now unused phase reflection list
2476	if phase in Histograms[histogram]['Reflection Lists']:
2477	del Histograms[histogram]['Reflection Lists'][phase]
2478	continue
2479	hapData = HistoPhase[histogram]
2480	hId = Histogram['hId']
2481	if 'PWDR' in histogram:
2482	limits = Histogram['Limits'][1]
2483	inst = Histogram['Instrument Parameters'][0] #TODO - grab table here if present
2484	if 'C' in inst['Type'][1]:
2485	try:
2486	wave = inst['Lam'][1]
2487	except KeyError:
2488	wave = inst['Lam1'][1]
2489	dmin = wave/(2.0*sind(limits[1]/2.0))
2490	elif 'T' in inst['Type'][0]:
2491	dmin = limits[0]/inst['difC'][1]
2492	else:
2493	wave = inst['Lam'][1]
2494	dmin = wave/(2.0*sind(limits[1]/2.0))
2495	pfx = str(pId)+':'+str(hId)+':'
2496	if Phases[phase]['General']['doPawley']:
2497	hapDict[pfx+'LeBail'] = False #Pawley supercedes LeBail
2498	hapDict[pfx+'newLeBail'] = True
2499	Tmin = G2lat.Dsp2pos(inst,dmin)
2500	if 'T' in inst['Type'][1]:
2501	limits[0] = max(limits[0],Tmin)
2502	else:
2503	limits[1] = min(limits[1],Tmin)
2504	else:
2505	hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2506	hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2507	if Phases[phase]['General']['Type'] == 'magnetic':
2508	dmin = max(dmin,Phases[phase]['General'].get('MagDmin',0.))
2509	for item in ['Scale','Extinction']:
2510	hapDict[pfx+item] = hapData[item][0]
2511	if hapData[item][1] and not hapDict[pfx+'LeBail']:
2512	hapVary.append(pfx+item)
2513	names = G2spc.HStrainNames(SGData)
2514	HSvals = []
2515	for i,name in enumerate(names):
2516	hapDict[pfx+name] = hapData['HStrain'][0][i]
2517	HSvals.append(hapDict[pfx+name])
2518	if hapData['HStrain'][1][i]:
2519	hapVary.append(pfx+name)
2520	if 'Layer Disp' in hapData:
2521	hapDict[pfx+'LayerDisp'] = hapData['Layer Disp'][0]
2522	if hapData['Layer Disp'][1]:
2523	hapVary.append(pfx+'LayerDisp')
2524	else:
2525	hapDict[pfx+'LayerDisp'] = 0.0
2526	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2527	if hapData['Pref.Ori.'][0] == 'MD':
2528	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2529	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2530	if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2531	hapVary.append(pfx+'MD')
2532	else: #'SH' spherical harmonics
2533	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2534	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2535	controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2536	controlDict[pfx+'SHhkl'] = []
2537	try: #patch for old Pref.Ori. items
2538	controlDict[pfx+'SHtoler'] = 0.1
2539	if hapData['Pref.Ori.'][6][0] != '':
2540	controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2541	controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2542	except IndexError:
2543	pass
2544	for item in hapData['Pref.Ori.'][5]:
2545	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2546	if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2547	hapVary.append(pfx+item)
2548	for item in ['Mustrain','Size']:
2549	controlDict[pfx+item+'Type'] = hapData[item][0]
2550	hapDict[pfx+item+';mx'] = hapData[item][1][2]
2551	if hapData[item][2][2]:
2552	hapVary.append(pfx+item+';mx')
2553	if hapData[item][0] in ['isotropic','uniaxial']:
2554	hapDict[pfx+item+';i'] = hapData[item][1][0]
2555	if hapData[item][2][0]:
2556	hapVary.append(pfx+item+';i')
2557	if hapData[item][0] == 'uniaxial':
2558	controlDict[pfx+item+'Axis'] = hapData[item][3]
2559	hapDict[pfx+item+';a'] = hapData[item][1][1]
2560	if hapData[item][2][1]:
2561	hapVary.append(pfx+item+';a')
2562	else: #generalized for mustrain or ellipsoidal for size
2563	Nterms = len(hapData[item][4])
2564	if item == 'Mustrain':
2565	names = G2spc.MustrainNames(SGData)
2566	pwrs = []
2567	for name in names:
2568	h,k,l = name[1:]
2569	pwrs.append([int(h),int(k),int(l)])
2570	controlDict[pfx+'MuPwrs'] = pwrs
2571	for i in range(Nterms):
2572	sfx = ';'+str(i)
2573	hapDict[pfx+item+sfx] = hapData[item][4][i]
2574	if hapData[item][5][i]:
2575	hapVary.append(pfx+item+sfx)
2576	if Phases[phase]['General']['Type'] != 'magnetic':
2577	for bab in ['BabA','BabU']:
2578	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2579	if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2580	hapVary.append(pfx+bab)
2581
2582	if Print:
2583	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2584	pFile.write(135*'='+'\n')
2585	if hapDict[pfx+'LeBail']:
2586	pFile.write(' Perform LeBail extraction\n')
2587	else:
2588	pFile.write(' Phase fraction : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2589	pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2590	if hapData['Pref.Ori.'][0] == 'MD':
2591	Ax = hapData['Pref.Ori.'][3]
2592	pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2593	(hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2594	else: #'SH' for spherical harmonics
2595	PrintSHPO(hapData['Pref.Ori.'])
2596	pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2597	PrintSize(hapData['Size'])
2598	PrintMuStrain(hapData['Mustrain'],SGData)
2599	PrintHStrain(hapData['HStrain'],SGData)
2600	if 'Layer Disp' in hapData:
2601	pFile.write(' Layer Displacement: %10.3f Refine? %s\n'%(hapData['Layer Disp'][0],hapData['Layer Disp'][1]))
2602	if Phases[phase]['General']['Type'] != 'magnetic':
2603	if hapData['Babinet']['BabA'][0]:
2604	PrintBabinet(hapData['Babinet'])
2605	if phase in Histogram['Reflection Lists'] and 'RefList' not in Histogram['Reflection Lists'][phase] and hapData.get('LeBail',False):
2606	hapData['newLeBail'] = True
2607	if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData.get('LeBail',False) and hapData['newLeBail'])):
2608	if hapData.get('LeBail',True): #stop regeneating reflections for LeBail
2609	hapData['newLeBail'] = False
2610	refList = []
2611	# Uniq = []
2612	# Phi = []
2613	useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2614	if Phases[phase]['General'].get('Modulated',False):
2615	ifSuper = True
2616	HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2617	HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2618	for h,k,l,m,d in HKLd:
2619	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2620	mul *= 2 # for powder overlap of Friedel pairs
2621	if m or not ext or useExt:
2622	if 'C' in inst['Type'][0]:
2623	pos = G2lat.Dsp2pos(inst,d)
2624	if limits[0] < pos < limits[1]:
2625	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])
2626	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2627	# Uniq.append(uniq)
2628	# Phi.append(phi)
2629	elif 'T' in inst['Type'][0]:
2630	pos = G2lat.Dsp2pos(inst,d)
2631	if limits[0] < pos < limits[1]:
2632	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2633	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])
2634	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2635	#TODO - if tabulated put alp & bet in here
2636	# Uniq.append(uniq)
2637	# Phi.append(phi)
2638	elif 'B' in inst['Type'][0]:
2639	pos = G2lat.Dsp2pos(inst,d)
2640	if limits[0] < pos < limits[1]:
2641	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])
2642	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2643	else:
2644	ifSuper = False
2645	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2646	HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2647	for h,k,l,d in HKLd:
2648	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2649	if ext and 'N' in inst['Type'][0] and 'MagSpGrp' in SGData:
2650	ext = G2spc.checkMagextc([h,k,l],SGData)
2651	mul *= 2 # for powder overlap of Friedel pairs
2652	if ext and not useExt:
2653	continue
2654	if 'C' in inst['Type'][0]:
2655	pos = G2lat.Dsp2pos(inst,d)
2656	if limits[0] < pos < limits[1]:
2657	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])
2658	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2659	# Uniq.append(uniq)
2660	# Phi.append(phi)
2661	elif 'T' in inst['Type'][0]:
2662	pos = G2lat.Dsp2pos(inst,d)
2663	if limits[0] < pos < limits[1]:
2664	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2665	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])
2666	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2667	# Uniq.append(uniq)
2668	# Phi.append(phi)
2669	elif 'B' in inst['Type'][0]:
2670	pos = G2lat.Dsp2pos(inst,d)
2671	if limits[0] < pos < limits[1]:
2672	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])
2673	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2674	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2675	elif 'HKLF' in histogram:
2676	inst = Histogram['Instrument Parameters'][0]
2677	hId = Histogram['hId']
2678	hfx = ':%d:'%(hId)
2679	for item in inst:
2680	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2681	hapDict[hfx+item] = inst[item][1]
2682	pfx = str(pId)+':'+str(hId)+':'
2683	hapDict[pfx+'Scale'] = hapData['Scale'][0]
2684	if hapData['Scale'][1]:
2685	hapVary.append(pfx+'Scale')
2686
2687	extApprox,extType,extParms = hapData['Extinction']
2688	controlDict[pfx+'EType'] = extType
2689	controlDict[pfx+'EApprox'] = extApprox
2690	if 'C' in inst['Type'][0]:
2691	controlDict[pfx+'Tbar'] = extParms['Tbar']
2692	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2693	if 'Primary' in extType:
2694	Ekey = ['Ep',]
2695	elif 'I & II' in extType:
2696	Ekey = ['Eg','Es']
2697	elif 'Secondary Type II' == extType:
2698	Ekey = ['Es',]
2699	elif 'Secondary Type I' == extType:
2700	Ekey = ['Eg',]
2701	else: #'None'
2702	Ekey = []
2703	for eKey in Ekey:
2704	hapDict[pfx+eKey] = extParms[eKey][0]
2705	if extParms[eKey][1]:
2706	hapVary.append(pfx+eKey)
2707	for bab in ['BabA','BabU']:
2708	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2709	if hapData['Babinet'][bab][1]:
2710	hapVary.append(pfx+bab)
2711	Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2712	if len(Twins) == 1:
2713	hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2714	if hapData.get('Flack',[0,False])[1]:
2715	hapVary.append(pfx+'Flack')
2716	sumTwFr = 0.
2717	controlDict[pfx+'TwinLaw'] = []
2718	controlDict[pfx+'TwinInv'] = []
2719	NTL = 0
2720	for it,twin in enumerate(Twins):
2721	if 'bool' in str(type(twin[0])):
2722	controlDict[pfx+'TwinInv'].append(twin[0])
2723	controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2724	else:
2725	NTL += 1
2726	controlDict[pfx+'TwinInv'].append(False)
2727	controlDict[pfx+'TwinLaw'].append(twin[0])
2728	if it:
2729	hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2730	sumTwFr += twin[1]
2731	else:
2732	hapDict[pfx+'TwinFr:0'] = twin[1][0]
2733	controlDict[pfx+'TwinNMN'] = twin[1][1]
2734	if Twins[0][1][1]:
2735	hapVary.append(pfx+'TwinFr:'+str(it))
2736	controlDict[pfx+'NTL'] = NTL
2737	#need to make constraint on TwinFr
2738	controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2739	if len(Twins) > 1: #force sum to unity
2740	hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2741	if Print:
2742	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2743	pFile.write(135*'='+'\n')
2744	pFile.write(' Scale factor : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2745	if extType != 'None':
2746	pFile.write(' Extinction Type: %15s approx: %10s\n'%(extType,extApprox))
2747	text = ' Parameters :'
2748	for eKey in Ekey:
2749	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2750	pFile.write(text+'\n')
2751	if hapData['Babinet']['BabA'][0]:
2752	PrintBabinet(hapData['Babinet'])
2753	if not SGData['SGInv'] and len(Twins) == 1:
2754	pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2755	if len(Twins) > 1:
2756	for it,twin in enumerate(Twins):
2757	if 'bool' in str(type(twin[0])):
2758	pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'%
2759	(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1])))
2760	else:
2761	pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2762	(str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1])))
2763
2764	Histogram['Reflection Lists'] = phase
2765
2766	return hapVary,hapDict,controlDict
2767
2768	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,Print=True,pFile=None):
2769	'needs a doc string'
2770
2771	def PrintSizeAndSig(hapData,sizeSig):
2772	line = '\n Size model: %9s'%(hapData[0])
2773	refine = False
2774	if hapData[0] in ['isotropic','uniaxial']:
2775	line += ' equatorial:%12.5f'%(hapData[1][0])
2776	if sizeSig[0][0]:
2777	line += ', sig:%8.4f'%(sizeSig[0][0])
2778	refine = True
2779	if hapData[0] == 'uniaxial':
2780	line += ' axial:%12.4f'%(hapData[1][1])
2781	if sizeSig[0][1]:
2782	refine = True
2783	line += ', sig:%8.4f'%(sizeSig[0][1])
2784	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2785	if sizeSig[0][2]:
2786	refine = True
2787	line += ', sig:%8.4f'%(sizeSig[0][2])
2788	if refine:
2789	pFile.write(line+'\n')
2790	else:
2791	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2792	if sizeSig[0][2]:
2793	refine = True
2794	line += ', sig:%8.4f'%(sizeSig[0][2])
2795	Snames = ['S11','S22','S33','S12','S13','S23']
2796	ptlbls = ' name :'
2797	ptstr = ' value :'
2798	sigstr = ' sig :'
2799	for i,name in enumerate(Snames):
2800	ptlbls += '%12s' % (name)
2801	ptstr += '%12.6f' % (hapData[4][i])
2802	if sizeSig[1][i]:
2803	refine = True
2804	sigstr += '%12.6f' % (sizeSig[1][i])
2805	else:
2806	sigstr += 12*' '
2807	if refine:
2808	pFile.write(line+'\n')
2809	pFile.write(ptlbls+'\n')
2810	pFile.write(ptstr+'\n')
2811	pFile.write(sigstr+'\n')
2812
2813	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2814	line = '\n Mustrain model: %9s\n'%(hapData[0])
2815	refine = False
2816	if hapData[0] in ['isotropic','uniaxial']:
2817	line += ' equatorial:%12.1f'%(hapData[1][0])
2818	if mustrainSig[0][0]:
2819	line += ', sig:%8.1f'%(mustrainSig[0][0])
2820	refine = True
2821	if hapData[0] == 'uniaxial':
2822	line += ' axial:%12.1f'%(hapData[1][1])
2823	if mustrainSig[0][1]:
2824	line += ', sig:%8.1f'%(mustrainSig[0][1])
2825	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2826	if mustrainSig[0][2]:
2827	refine = True
2828	line += ', sig:%8.3f'%(mustrainSig[0][2])
2829	if refine:
2830	pFile.write(line+'\n')
2831	else:
2832	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2833	if mustrainSig[0][2]:
2834	refine = True
2835	line += ', sig:%8.3f'%(mustrainSig[0][2])
2836	Snames = G2spc.MustrainNames(SGData)
2837	ptlbls = ' name :'
2838	ptstr = ' value :'
2839	sigstr = ' sig :'
2840	for i,name in enumerate(Snames):
2841	ptlbls += '%12s' % (name)
2842	ptstr += '%12.1f' % (hapData[4][i])
2843	if mustrainSig[1][i]:
2844	refine = True
2845	sigstr += '%12.1f' % (mustrainSig[1][i])
2846	else:
2847	sigstr += 12*' '
2848	if refine:
2849	pFile.write(line+'\n')
2850	pFile.write(ptlbls+'\n')
2851	pFile.write(ptstr+'\n')
2852	pFile.write(sigstr+'\n')
2853
2854	def PrintHStrainAndSig(hapData,strainSig,SGData):
2855	Hsnames = G2spc.HStrainNames(SGData)
2856	ptlbls = ' name :'
2857	ptstr = ' value :'
2858	sigstr = ' sig :'
2859	refine = False
2860	for i,name in enumerate(Hsnames):
2861	ptlbls += '%12s' % (name)
2862	ptstr += '%12.4g' % (hapData[0][i])
2863	if name in strainSig:
2864	refine = True
2865	sigstr += '%12.4g' % (strainSig[name])
2866	else:
2867	sigstr += 12*' '
2868	if refine:
2869	pFile.write('\n Hydrostatic/elastic strain:\n')
2870	pFile.write(ptlbls+'\n')
2871	pFile.write(ptstr+'\n')
2872	pFile.write(sigstr+'\n')
2873
2874	def PrintSHPOAndSig(pfx,hapData,POsig):
2875	Tindx = 1.0
2876	Tvar = 0.0
2877	pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2878	ptlbls = ' names :'
2879	ptstr = ' values:'
2880	sigstr = ' sig :'
2881	for item in hapData[5]:
2882	ptlbls += '%12s'%(item)
2883	ptstr += '%12.3f'%(hapData[5][item])
2884	l = 2.0*eval(item.strip('C'))[0]+1
2885	Tindx += hapData[5][item]**2/l
2886	if pfx+item in POsig:
2887	Tvar += (2.hapData[5][item]POsig[pfx+item]/l)**2
2888	sigstr += '%12.3f'%(POsig[pfx+item])
2889	else:
2890	sigstr += 12*' '
2891	pFile.write(ptlbls+'\n')
2892	pFile.write(ptstr+'\n')
2893	pFile.write(sigstr+'\n')
2894	pFile.write('\n Texture index J = %.3f(%d)\n'%(Tindx,int(1000*np.sqrt(Tvar))))
2895
2896	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2897	pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2898	text = ''
2899	for item in hapData[2]:
2900	if pfx+item in ScalExtSig:
2901	text += ' %s: '%(item)
2902	text += '%12.2e'%(hapData[2][item][0])
2903	if pfx+item in ScalExtSig:
2904	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2905	pFile.write(text+'\n')
2906
2907	def PrintBabinetAndSig(pfx,hapData,BabSig):
2908	pFile.write('\n Babinet form factor modification:\n')
2909	ptlbls = ' names :'
2910	ptstr = ' values:'
2911	sigstr = ' sig :'
2912	for item in hapData:
2913	ptlbls += '%12s'%(item)
2914	ptstr += '%12.3f'%(hapData[item][0])
2915	if pfx+item in BabSig:
2916	sigstr += '%12.3f'%(BabSig[pfx+item])
2917	else:
2918	sigstr += 12*' '
2919	pFile.write(ptlbls+'\n')
2920	pFile.write(ptstr+'\n')
2921	pFile.write(sigstr+'\n')
2922
2923	def PrintTwinsAndSig(pfx,twinData,TwinSig):
2924	pFile.write('\n Twin Law fractions :\n')
2925	ptlbls = ' names :'
2926	ptstr = ' values:'
2927	sigstr = ' sig :'
2928	for it,item in enumerate(twinData):
2929	ptlbls += ' twin #%d'%(it)
2930	if it:
2931	ptstr += '%12.3f'%(item[1])
2932	else:
2933	ptstr += '%12.3f'%(item[1][0])
2934	if pfx+'TwinFr:'+str(it) in TwinSig:
2935	sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2936	else:
2937	sigstr += 12*' '
2938	pFile.write(ptlbls+'\n')
2939	pFile.write(ptstr+'\n')
2940	pFile.write(sigstr+'\n')
2941
2942
2943	PhFrExtPOSig = {}
2944	SizeMuStrSig = {}
2945	ScalExtSig = {}
2946	BabSig = {}
2947	TwinFrSig = {}
2948	wtFrSum = {}
2949	for phase in Phases:
2950	HistoPhase = Phases[phase]['Histograms']
2951	General = Phases[phase]['General']
2952	SGData = General['SGData']
2953	pId = Phases[phase]['pId']
2954	histoList = list(HistoPhase.keys())
2955	histoList.sort()
2956	for histogram in histoList:
2957	try:
2958	Histogram = Histograms[histogram]
2959	except KeyError:
2960	#skip if histogram not included e.g. in a sequential refinement
2961	continue
2962	if not Phases[phase]['Histograms'][histogram]['Use']:
2963	#skip if phase absent from this histogram
2964	continue
2965	hapData = HistoPhase[histogram]
2966	hId = Histogram['hId']
2967	pfx = str(pId)+':'+str(hId)+':'
2968	if hId not in wtFrSum:
2969	wtFrSum[hId] = 0.
2970	if 'PWDR' in histogram:
2971	for item in ['Scale','Extinction']:
2972	hapData[item][0] = parmDict[pfx+item]
2973	if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2974	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2975	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2976	if hapData['Pref.Ori.'][0] == 'MD':
2977	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2978	if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2979	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2980	else: #'SH' spherical harmonics
2981	for item in hapData['Pref.Ori.'][5]:
2982	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2983	if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2984	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2985	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2986	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2987	pfx+'HStrain':{}})
2988	for item in ['Mustrain','Size']:
2989	hapData[item][1][2] = parmDict[pfx+item+';mx']
2990	# hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2991	if pfx+item+';mx' in sigDict:
2992	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2993	if hapData[item][0] in ['isotropic','uniaxial']:
2994	hapData[item][1][0] = parmDict[pfx+item+';i']
2995	if item == 'Size':
2996	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2997	if pfx+item+';i' in sigDict:
2998	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2999	if hapData[item][0] == 'uniaxial':
3000	hapData[item][1][1] = parmDict[pfx+item+';a']
3001	if item == 'Size':
3002	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
3003	if pfx+item+';a' in sigDict:
3004	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
3005	else: #generalized for mustrain or ellipsoidal for size
3006	Nterms = len(hapData[item][4])
3007	for i in range(Nterms):
3008	sfx = ';'+str(i)
3009	hapData[item][4][i] = parmDict[pfx+item+sfx]
3010	if pfx+item+sfx in sigDict:
3011	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
3012	names = G2spc.HStrainNames(SGData)
3013	for i,name in enumerate(names):
3014	hapData['HStrain'][0][i] = parmDict[pfx+name]
3015	if pfx+name in sigDict:
3016	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
3017	if 'Layer Disp' in hapData:
3018	hapData['Layer Disp'][0] = parmDict[pfx+'LayerDisp']
3019	if pfx+'LayerDisp' in sigDict:
3020	SizeMuStrSig[pfx+'LayerDisp'] = sigDict[pfx+'LayerDisp']
3021	if Phases[phase]['General']['Type'] != 'magnetic':
3022	for name in ['BabA','BabU']:
3023	hapData['Babinet'][name][0] = parmDict[pfx+name]
3024	if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
3025	BabSig[pfx+name] = sigDict[pfx+name]
3026
3027	elif 'HKLF' in histogram:
3028	for item in ['Scale','Flack']:
3029	if parmDict.get(pfx+item):
3030	hapData[item][0] = parmDict[pfx+item]
3031	if pfx+item in sigDict:
3032	ScalExtSig[pfx+item] = sigDict[pfx+item]
3033	for item in ['Ep','Eg','Es']:
3034	if parmDict.get(pfx+item):
3035	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
3036	if pfx+item in sigDict:
3037	ScalExtSig[pfx+item] = sigDict[pfx+item]
3038	for item in ['BabA','BabU']:
3039	hapData['Babinet'][item][0] = parmDict[pfx+item]
3040	if pfx+item in sigDict:
3041	BabSig[pfx+item] = sigDict[pfx+item]
3042	if 'Twins' in hapData:
3043	it = 1
3044	sumTwFr = 0.
3045	while True:
3046	try:
3047	hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
3048	if pfx+'TwinFr:'+str(it) in sigDict:
3049	TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
3050	if it:
3051	sumTwFr += hapData['Twins'][it][1]
3052	it += 1
3053	except KeyError:
3054	break
3055	hapData['Twins'][0][1][0] = 1.-sumTwFr
3056
3057	if Print:
3058	for phase in Phases:
3059	HistoPhase = Phases[phase]['Histograms']
3060	General = Phases[phase]['General']
3061	SGData = General['SGData']
3062	pId = Phases[phase]['pId']
3063	histoList = list(HistoPhase.keys())
3064	histoList.sort()
3065	for histogram in histoList:
3066	try:
3067	Histogram = Histograms[histogram]
3068	except KeyError:
3069	#skip if histogram not included e.g. in a sequential refinement
3070	continue
3071	hapData = HistoPhase[histogram]
3072	hId = Histogram['hId']
3073	Histogram['Residuals'][str(pId)+'::Name'] = phase
3074	pfx = str(pId)+':'+str(hId)+':'
3075	hfx = ':%s:'%(hId)
3076	if pfx+'Nref' not in Histogram['Residuals']: #skip not used phase in histogram
3077	continue
3078	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
3079	pFile.write(135*'='+'\n')
3080	Inst = Histogram['Instrument Parameters'][0]
3081	if 'PWDR' in histogram:
3082	pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
3083	(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
3084	pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
3085	pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
3086
3087	if parmDict[pfx+'LeBail']:
3088	pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
3089	else:
3090	# if calcControls != None: #skipped in seqRefine
3091	# if 'X'in Inst['Type'][0]:
3092	# PrintFprime(calcControls['FFtables'],hfx,pFile)
3093	# elif 'NC' in Inst['Type'][0]:
3094	# PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3095	if pfx+'Scale' in PhFrExtPOSig:
3096	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
3097	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
3098	pFile.write(' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f\n'%
3099	(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
3100	if pfx+'Extinction' in PhFrExtPOSig:
3101	pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
3102	if hapData['Pref.Ori.'][0] == 'MD':
3103	if pfx+'MD' in PhFrExtPOSig:
3104	pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
3105	else:
3106	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
3107	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
3108	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
3109	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
3110	if pfx+'LayerDisp' in SizeMuStrSig:
3111	pFile.write(' Layer displacement : %10.3f, sig %10.3f\n'%(hapData['Layer Disp'][0],SizeMuStrSig[pfx+'LayerDisp']))
3112	if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
3113	if len(BabSig):
3114	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3115
3116	elif 'HKLF' in histogram:
3117	pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
3118	(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
3119	Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
3120	if calcControls != None: #skipped in seqRefine
3121	if 'X'in Inst['Type'][0]:
3122	PrintFprime(calcControls['FFtables'],hfx,pFile)
3123	elif 'NC' in Inst['Type'][0]:
3124	PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3125	pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
3126	if pfx+'Scale' in ScalExtSig:
3127	pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
3128	if hapData['Extinction'][0] != 'None':
3129	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
3130	if len(BabSig):
3131	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3132	if pfx+'Flack' in ScalExtSig:
3133	pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
3134	if len(TwinFrSig):
3135	PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
3136
3137	################################################################################
3138	##### Histogram data
3139	################################################################################
3140
3141	def GetHistogramData(Histograms,Print=True,pFile=None):
3142	'needs a doc string'
3143
3144	def GetBackgroundParms(hId,Background):
3145	Back = Background[0]
3146	DebyePeaks = Background[1]
3147	bakType,bakFlag = Back[:2]
3148	backVals = Back[3:]
3149	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
3150	backDict = dict(zip(backNames,backVals))
3151	backVary = []
3152	if bakFlag:
3153	backVary = backNames
3154	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
3155	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
3156	debyeDict = {}
3157	debyeList = []
3158	for i in range(DebyePeaks['nDebye']):
3159	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
3160	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
3161	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
3162	debyeVary = []
3163	for item in debyeList:
3164	if item[1]:
3165	debyeVary.append(item[0])
3166	backDict.update(debyeDict)
3167	backVary += debyeVary
3168	peakDict = {}
3169	peakList = []
3170	for i in range(DebyePeaks['nPeaks']):
3171	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
3172	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
3173	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
3174	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
3175	peakVary = []
3176	for item in peakList:
3177	if item[1]:
3178	peakVary.append(item[0])
3179	backDict.update(peakDict)
3180	backVary += peakVary
3181	if 'background PWDR' in Background[1]:
3182	backDict[':'+str(hId)+':Back File'] = Background[1]['background PWDR'][0]
3183	backDict[':'+str(hId)+':BF mult'] = Background[1]['background PWDR'][1]
3184	try:
3185	if Background[1]['background PWDR'][0] and Background[1]['background PWDR'][2]:
3186	backVary.append(':'+str(hId)+':BF mult')
3187	except IndexError: # old version without refine flag
3188	pass
3189	return bakType,backDict,backVary
3190
3191	def GetInstParms(hId,Inst):
3192	#patch
3193	if 'Z' not in Inst:
3194	Inst['Z'] = [0.0,0.0,False]
3195	dataType = Inst['Type'][0]
3196	instDict = {}
3197	insVary = []
3198	pfx = ':'+str(hId)+':'
3199	insKeys = list(Inst.keys())
3200	insKeys.sort()
3201	for item in insKeys:
3202	if item in 'Azimuth':
3203	continue
3204	insName = pfx+item
3205	instDict[insName] = Inst[item][1]
3206	if len(Inst[item]) > 2 and Inst[item][2]:
3207	insVary.append(insName)
3208	if 'C' in dataType:
3209	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
3210	elif 'T' in dataType: #trap zero alp, bet coeff.
3211	if not instDict[pfx+'alpha']:
3212	instDict[pfx+'alpha'] = 1.0
3213	if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3214	instDict[pfx+'beta-1'] = 1.0
3215	elif 'B' in dataType: #trap zero alp, bet coeff.
3216	if not instDict[pfx+'alpha-0'] and not instDict[pfx+'alpha-1']:
3217	instDict[pfx+'alpha-1'] = 1.0
3218	if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3219	instDict[pfx+'beta-1'] = 1.0
3220	return dataType,instDict,insVary
3221
3222	def GetSampleParms(hId,Sample):
3223	sampVary = []
3224	hfx = ':'+str(hId)+':'
3225	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
3226	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi'],hfx+'Azimuth':Sample['Azimuth']}
3227	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
3228	if key in Sample:
3229	sampDict[hfx+key] = Sample[key]
3230	Type = Sample['Type']
3231	if 'Bragg' in Type: #Bragg-Brentano
3232	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3233	sampDict[hfx+item] = Sample[item][0]
3234	if Sample[item][1]:
3235	sampVary.append(hfx+item)
3236	elif 'Debye' in Type: #Debye-Scherrer
3237	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3238	sampDict[hfx+item] = Sample[item][0]
3239	if Sample[item][1]:
3240	sampVary.append(hfx+item)
3241	return Type,sampDict,sampVary
3242
3243	def PrintBackground(Background):
3244	Back = Background[0]
3245	DebyePeaks = Background[1]
3246	pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
3247	line = ' Coefficients: '
3248	for i,back in enumerate(Back[3:]):
3249	line += '%10.3f'%(back)
3250	if i and not i%10:
3251	line += '\n'+15*' '
3252	pFile.write(line+'\n')
3253	if DebyePeaks['nDebye']:
3254	pFile.write('\n Debye diffuse scattering coefficients\n')
3255	parms = ['DebyeA','DebyeR','DebyeU']
3256	line = ' names : '
3257	for parm in parms:
3258	line += '%8s refine?'%(parm)
3259	pFile.write(line+'\n')
3260	for j,term in enumerate(DebyePeaks['debyeTerms']):
3261	line = ' term'+'%2d'%(j)+':'
3262	for i in range(3):
3263	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
3264	pFile.write(line+'\n')
3265	if DebyePeaks['nPeaks']:
3266	pFile.write('\n Single peak coefficients\n')
3267	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3268	line = ' names : '
3269	for parm in parms:
3270	line += '%8s refine?'%(parm)
3271	pFile.write(line+'\n')
3272	for j,term in enumerate(DebyePeaks['peaksList']):
3273	line = ' peak'+'%2d'%(j)+':'
3274	for i in range(4):
3275	line += '%12.3f %5s'%(term[2i],bool(term[2i+1]))
3276	pFile.write(line+'\n')
3277	if 'background PWDR' in DebyePeaks:
3278	try:
3279	pFile.write(' Fixed background file: %s; mult: %.3f Refine? %s\n'%(DebyePeaks['background PWDR'][0],
3280	DebyePeaks['background PWDR'][1],DebyePeaks['background PWDR'][2]))
3281	except IndexError: #old version without refine flag
3282	pass
3283
3284	def PrintInstParms(Inst):
3285	pFile.write('\n Instrument Parameters:\n')
3286	insKeys = list(Inst.keys())
3287	insKeys.sort()
3288	iBeg = 0
3289	Ok = True
3290	while Ok:
3291	ptlbls = ' name :'
3292	ptstr = ' value :'
3293	varstr = ' refine:'
3294	iFin = min(iBeg+9,len(insKeys))
3295	for item in insKeys[iBeg:iFin]:
3296	if item not in ['Type','Source']:
3297	ptlbls += '%12s' % (item)
3298	ptstr += '%12.6f' % (Inst[item][1])
3299	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
3300	varstr += 12*' '
3301	else:
3302	varstr += '%12s' % (str(bool(Inst[item][2])))
3303	pFile.write(ptlbls+'\n')
3304	pFile.write(ptstr+'\n')
3305	pFile.write(varstr+'\n')
3306	iBeg = iFin
3307	if iBeg == len(insKeys):
3308	Ok = False
3309	else:
3310	pFile.write('\n')
3311
3312	def PrintSampleParms(Sample):
3313	pFile.write('\n Sample Parameters:\n')
3314	pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3315	(Sample['Omega'],Sample['Chi'],Sample['Phi']))
3316	ptlbls = ' name :'
3317	ptstr = ' value :'
3318	varstr = ' refine:'
3319	if 'Bragg' in Sample['Type']:
3320	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3321	ptlbls += '%14s'%(item)
3322	ptstr += '%14.4f'%(Sample[item][0])
3323	varstr += '%14s'%(str(bool(Sample[item][1])))
3324
3325	elif 'Debye' in Type: #Debye-Scherrer
3326	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3327	ptlbls += '%14s'%(item)
3328	ptstr += '%14.4f'%(Sample[item][0])
3329	varstr += '%14s'%(str(bool(Sample[item][1])))
3330
3331	pFile.write(ptlbls+'\n')
3332	pFile.write(ptstr+'\n')
3333	pFile.write(varstr+'\n')
3334
3335	histDict = {}
3336	histVary = []
3337	controlDict = {}
3338	histoList = list(Histograms.keys())
3339	histoList.sort()
3340	for histogram in histoList:
3341	if 'PWDR' in histogram:
3342	Histogram = Histograms[histogram]
3343	hId = Histogram['hId']
3344	pfx = ':'+str(hId)+':'
3345	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3346	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3347	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3348	for excl in controlDict[pfx+'Exclude']:
3349	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3350	if controlDict[pfx+'Exclude']:
3351	ma.mask_rows(Histogram['Data'])
3352	Background = Histogram['Background']
3353	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3354	controlDict[pfx+'bakType'] = Type
3355	histDict.update(bakDict)
3356	histVary += bakVary
3357
3358	Inst = Histogram['Instrument Parameters'] #TODO ? ignores tabulated alp,bet & delt for TOF
3359	if 'T' in Type and len(Inst[1]): #patch - back-to-back exponential contribution to TOF line shape is removed
3360	G2fil.G2Print ('Warning: tabulated profile coefficients are ignored')
3361	Type,instDict,insVary = GetInstParms(hId,Inst[0])
3362	controlDict[pfx+'histType'] = Type
3363	if 'XC' in Type or 'XB' in Type:
3364	if pfx+'Lam1' in instDict:
3365	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3366	else:
3367	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
3368	histDict.update(instDict)
3369	histVary += insVary
3370
3371	Sample = Histogram['Sample Parameters']
3372	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3373	controlDict[pfx+'instType'] = Type
3374	histDict.update(sampDict)
3375	histVary += sampVary
3376
3377
3378	if Print:
3379	pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3380	pFile.write(135*'='+'\n')
3381	Units = {'C':' deg','T':' msec','B':' deg'}
3382	units = Units[controlDict[pfx+'histType'][2]]
3383	Limits = controlDict[pfx+'Limits']
3384	pFile.write(' Instrument type: %s\n'%Sample['Type'])
3385	pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3386	if len(controlDict[pfx+'Exclude']):
3387	excls = controlDict[pfx+'Exclude']
3388	for excl in excls:
3389	pFile.write(' Excluded region: %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units))
3390	PrintSampleParms(Sample)
3391	PrintInstParms(Inst[0])
3392	PrintBackground(Background)
3393	elif 'HKLF' in histogram:
3394	Histogram = Histograms[histogram]
3395	hId = Histogram['hId']
3396	pfx = ':'+str(hId)+':'
3397	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3398	Inst = Histogram['Instrument Parameters'][0]
3399	controlDict[pfx+'histType'] = Inst['Type'][0]
3400	if 'X' in Inst['Type'][0]:
3401	histDict[pfx+'Lam'] = Inst['Lam'][1]
3402	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3403	elif 'NC' in Inst['Type'][0] or 'NB' in Inst['Type'][0]:
3404	histDict[pfx+'Lam'] = Inst['Lam'][1]
3405	return histVary,histDict,controlDict
3406
3407	def SetHistogramData(parmDict,sigDict,Histograms,calcControls,Print=True,pFile=None,seq=False):
3408	'Shows histogram data after a refinement'
3409
3410	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3411	Back = Background[0]
3412	DebyePeaks = Background[1]
3413	lenBack = len(Back[3:])
3414	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
3415	for i in range(lenBack):
3416	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3417	if pfx+'Back;'+str(i) in sigDict:
3418	backSig[i] = sigDict[pfx+'Back;'+str(i)]
3419	if DebyePeaks['nDebye']:
3420	for i in range(DebyePeaks['nDebye']):
3421	names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3422	for j,name in enumerate(names):
3423	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3424	if name in sigDict:
3425	backSig[lenBack+3*i+j] = sigDict[name]
3426	if DebyePeaks['nPeaks']:
3427	for i in range(DebyePeaks['nPeaks']):
3428	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3429	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3430	for j,name in enumerate(names):
3431	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3432	if name in sigDict:
3433	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
3434	if pfx+'BF mult' in sigDict:
3435	DebyePeaks['background PWDR'][1] = parmDict[pfx+'BF mult']
3436	backSig.append(sigDict[pfx+'BF mult'])
3437
3438	return backSig
3439
3440	def SetInstParms(pfx,Inst,parmDict,sigDict):
3441	instSig = {}
3442	insKeys = list(Inst.keys())
3443	insKeys.sort()
3444	for item in insKeys:
3445	insName = pfx+item
3446	Inst[item][1] = parmDict[insName]
3447	if insName in sigDict:
3448	instSig[item] = sigDict[insName]
3449	else:
3450	instSig[item] = 0
3451	return instSig
3452
3453	def SetSampleParms(pfx,Sample,parmDict,sigDict):
3454	if 'Bragg' in Sample['Type']: #Bragg-Brentano
3455	sampSig = [0 for i in range(5)]
3456	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3457	Sample[item][0] = parmDict[pfx+item]
3458	if pfx+item in sigDict:
3459	sampSig[i] = sigDict[pfx+item]
3460	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3461	sampSig = [0 for i in range(4)]
3462	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3463	Sample[item][0] = parmDict[pfx+item]
3464	if pfx+item in sigDict:
3465	sampSig[i] = sigDict[pfx+item]
3466	return sampSig
3467
3468	def PrintBackgroundSig(Background,backSig):
3469	Back = Background[0]
3470	DebyePeaks = Background[1]
3471	valstr = ' value : '
3472	sigstr = ' sig : '
3473	refine = False
3474	for i,back in enumerate(Back[3:]):
3475	valstr += '%10.4g'%(back)
3476	if Back[1]:
3477	refine = True
3478	sigstr += '%10.4g'%(backSig[i])
3479	else:
3480	sigstr += 10*' '
3481	if refine:
3482	pFile.write('\n Background function: %s\n'%Back[0])
3483	pFile.write(valstr+'\n')
3484	pFile.write(sigstr+'\n')
3485	if DebyePeaks['nDebye']:
3486	ifAny = False
3487	ptfmt = "%12.3f"
3488	names = ' names :'
3489	ptstr = ' values:'
3490	sigstr = ' esds :'
3491	for item in sigDict:
3492	if 'Debye' in item:
3493	ifAny = True
3494	names += '%12s'%(item)
3495	ptstr += ptfmt%(parmDict[item])
3496	sigstr += ptfmt%(sigDict[item])
3497	if ifAny:
3498	pFile.write('\n Debye diffuse scattering coefficients\n')
3499	pFile.write(names+'\n')
3500	pFile.write(ptstr+'\n')
3501	pFile.write(sigstr+'\n')
3502	if DebyePeaks['nPeaks']:
3503	pFile.write('\n Single peak coefficients:\n')
3504	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3505	line = ' peak no. '
3506	for parm in parms:
3507	line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3508	pFile.write(line+'\n')
3509	for ip in range(DebyePeaks['nPeaks']):
3510	ptstr = ' %4d '%(ip)
3511	for parm in parms:
3512	fmt = '%14.3f'
3513	efmt = '%12.3f'
3514	if parm == 'BkPkpos':
3515	fmt = '%14.4f'
3516	efmt = '%12.4f'
3517	name = pfx+parm+';%d'%(ip)
3518	ptstr += fmt%(parmDict[name])
3519	if name in sigDict:
3520	ptstr += efmt%(sigDict[name])
3521	else:
3522	ptstr += 12*' '
3523	pFile.write(ptstr+'\n')
3524	if 'background PWDR' in DebyePeaks and DebyePeaks['background PWDR'][2]:
3525	pFile.write(' Fixed background scale: %.3f(%d)\n'%(DebyePeaks['background PWDR'][1],int(1000*backSig[-1])))
3526	sumBk = np.array(Histogram['sumBk'])
3527	pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3528	(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3529
3530	def PrintInstParmsSig(Inst,instSig):
3531	refine = False
3532	insKeys = list(instSig.keys())
3533	insKeys.sort()
3534	iBeg = 0
3535	Ok = True
3536	while Ok:
3537	ptlbls = ' names :'
3538	ptstr = ' value :'
3539	sigstr = ' sig :'
3540	iFin = min(iBeg+9,len(insKeys))
3541	for name in insKeys[iBeg:iFin]:
3542	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3543	ptlbls += '%12s' % (name)
3544	ptstr += '%12.6f' % (Inst[name][1])
3545	if instSig[name]:
3546	refine = True
3547	sigstr += '%12.6f' % (instSig[name])
3548	else:
3549	sigstr += 12*' '
3550	if refine:
3551	pFile.write('\n Instrument Parameters:\n')
3552	pFile.write(ptlbls+'\n')
3553	pFile.write(ptstr+'\n')
3554	pFile.write(sigstr+'\n')
3555	iBeg = iFin
3556	if iBeg == len(insKeys):
3557	Ok = False
3558
3559	def PrintSampleParmsSig(Sample,sampleSig):
3560	ptlbls = ' names :'
3561	ptstr = ' values:'
3562	sigstr = ' sig :'
3563	refine = False
3564	if 'Bragg' in Sample['Type']:
3565	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3566	ptlbls += '%14s'%(item)
3567	ptstr += '%14.4f'%(Sample[item][0])
3568	if sampleSig[i]:
3569	refine = True
3570	sigstr += '%14.4f'%(sampleSig[i])
3571	else:
3572	sigstr += 14*' '
3573
3574	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3575	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3576	ptlbls += '%14s'%(item)
3577	ptstr += '%14.4f'%(Sample[item][0])
3578	if sampleSig[i]:
3579	refine = True
3580	sigstr += '%14.4f'%(sampleSig[i])
3581	else:
3582	sigstr += 14*' '
3583
3584	if refine:
3585	pFile.write('\n Sample Parameters:\n')
3586	pFile.write(ptlbls+'\n')
3587	pFile.write(ptstr+'\n')
3588	pFile.write(sigstr+'\n')
3589
3590	histoList = list(Histograms.keys())
3591	histoList.sort()
3592	for histogram in histoList:
3593	if 'PWDR' in histogram:
3594	Histogram = Histograms[histogram]
3595	hId = Histogram['hId']
3596	pfx = ':'+str(hId)+':'
3597	Background = Histogram['Background']
3598	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3599
3600	Inst = Histogram['Instrument Parameters'][0]
3601	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3602
3603	Sample = Histogram['Sample Parameters']
3604	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3605
3606	if not seq:
3607	pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3608	pFile.write(135*'='+'\n')
3609	pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3610	pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3611	(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3612	pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3613	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3614	if Print:
3615	pFile.write(' Instrument type: %s\n'%Sample['Type'])
3616	if calcControls != None: #skipped in seqRefine
3617	if 'X' in Inst['Type'][0]:
3618	PrintFprime(calcControls['FFtables'],pfx,pFile)
3619	elif 'NC' in Inst['Type'][0]:
3620	PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3621	PrintSampleParmsSig(Sample,sampSig)
3622	PrintInstParmsSig(Inst,instSig)
3623	PrintBackgroundSig(Background,backSig)
3624
3625	def WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
3626	'''Cribbed version of PrintRBObjPOAndSig but returns lists of strings.
3627	Moved so it can be used in ExportCIF
3628	'''
3629	namstr = ' names :'
3630	valstr = ' values:'
3631	sigstr = ' esds :'
3632	for i,px in enumerate(['Px:','Py:','Pz:']):
3633	name = pfx+rbfx+px+rbsx
3634	namstr += '%12s'%('Pos '+px[1])
3635	valstr += '%12.5f'%(parmDict[name])
3636	if name in sigDict:
3637	sigstr += '%12.5f'%(sigDict[name])
3638	else:
3639	sigstr += 12*' '
3640	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
3641	name = pfx+rbfx+po+rbsx
3642	namstr += '%12s'%('Ori '+po[1])
3643	valstr += '%12.5f'%(parmDict[name])
3644	if name in sigDict:
3645	sigstr += '%12.5f'%(sigDict[name])
3646	else:
3647	sigstr += 12*' '
3648	return (namstr,valstr,sigstr)
3649
3650	def WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
3651	'''Cribbed version of PrintRBObjTLSAndSig but returns lists of strings.
3652	Moved so it can be used in ExportCIF
3653	'''
3654	out = []
3655	namstr = ' names :'
3656	valstr = ' values:'
3657	sigstr = ' esds :'
3658	if 'N' not in TLS:
3659	out.append(' Thermal motion:\n')
3660	if 'T' in TLS:
3661	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
3662	name = pfx+rbfx+pt+rbsx
3663	namstr += '%12s'%(pt[:3])
3664	valstr += '%12.5f'%(parmDict[name])
3665	if name in sigDict:
3666	sigstr += '%12.5f'%(sigDict[name])
3667	else:
3668	sigstr += 12*' '
3669	out.append(namstr+'\n')
3670	out.append(valstr+'\n')
3671	out.append(sigstr+'\n')
3672	if 'L' in TLS:
3673	namstr = ' names :'
3674	valstr = ' values:'
3675	sigstr = ' esds :'
3676	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
3677	name = pfx+rbfx+pt+rbsx
3678	namstr += '%12s'%(pt[:3])
3679	valstr += '%12.3f'%(parmDict[name])
3680	if name in sigDict:
3681	sigstr += '%12.3f'%(sigDict[name])
3682	else:
3683	sigstr += 12*' '
3684	out.append(namstr+'\n')
3685	out.append(valstr+'\n')
3686	out.append(sigstr+'\n')
3687	if 'S' in TLS:
3688	namstr = ' names :'
3689	valstr = ' values:'
3690	sigstr = ' esds :'
3691	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
3692	name = pfx+rbfx+pt+rbsx
3693	namstr += '%12s'%(pt[:3])
3694	valstr += '%12.4f'%(parmDict[name])
3695	if name in sigDict:
3696	sigstr += '%12.4f'%(sigDict[name])
3697	else:
3698	sigstr += 12*' '
3699	out.append(namstr+'\n')
3700	out.append(valstr+'\n')
3701	out.append(sigstr+'\n')
3702	if 'U' in TLS:
3703	name = pfx+rbfx+'U:'+rbsx
3704	namstr = ' names :'+'%12s'%('Uiso')
3705	valstr = ' values:'+'%12.5f'%(parmDict[name])
3706	if name in sigDict:
3707	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
3708	else:
3709	sigstr = ' esds :'+12*' '
3710	out.append(namstr+'\n')
3711	out.append(valstr+'\n')
3712	out.append(sigstr+'\n')
3713	return out
3714
3715	def WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
3716	'''Cribbed version of PrintRBObjTorAndSig but returns lists of strings.
3717	Moved so it can be used in ExportCIF
3718	'''
3719	out = []
3720	namstr = ' names :'
3721	valstr = ' values:'
3722	sigstr = ' esds :'
3723	out.append(' Torsions:\n')
3724	for it in range(nTors):
3725	name = pfx+'RBRTr;'+str(it)+':'+rbsx
3726	namstr += '%12s'%('Tor'+str(it))
3727	valstr += '%12.4f'%(parmDict[name])
3728	if name in sigDict:
3729	sigstr += '%12.4f'%(sigDict[name])
3730	out.append(namstr+'\n')
3731	out.append(valstr+'\n')
3732	out.append(sigstr+'\n')
3733	return out
3734
3735	def WriteResRBModel(RBModel):
3736	'''Write description of a residue rigid body. Code shifted from
3737	PrintResRBModel to make usable from G2export_CIF
3738	'''
3739	out = []
3740	atNames = RBModel['atNames']
3741	rbRef = RBModel['rbRef']
3742	rbSeq = RBModel['rbSeq']
3743	out.append(' At name x y z\n')
3744	for name,xyz in zip(atNames,RBModel['rbXYZ']):
3745	out.append(' %8s %10.4f %10.4f %10.4f\n'%(name,xyz[0],xyz[1],xyz[2]))
3746	out.append(' Orientation defined by: %s -> %s & %s -> %s\n'%
3747	(atNames[rbRef[0]],atNames[rbRef[1]],atNames[rbRef[0]],atNames[rbRef[2]]))
3748	if rbSeq:
3749	for i,rbseq in enumerate(rbSeq):
3750	# nameLst = [atNames[i] for i in rbseq[3]]
3751	out.append(' Torsion sequence %d Bond: %s - %s riding: %s\n'%
3752	(i,atNames[rbseq[0]],atNames[rbseq[1]],str([atNames[i] for i in rbseq[3]])))
3753	return out
3754
3755	def WriteVecRBModel(RBModel,sigDict={},irb=None):
3756	'''Write description of a vector rigid body. Code shifted from
3757	PrintVecRBModel to make usable from G2export_CIF
3758	'''
3759	out = []
3760	rbRef = RBModel['rbRef']
3761	atTypes = RBModel['rbTypes']
3762	for i in range(len(RBModel['VectMag'])):
3763	if sigDict and irb is not None:
3764	name = '::RBV;'+str(i)+':'+str(irb)
3765	s = G2mth.ValEsd(RBModel['VectMag'][i],sigDict.get(name,-.0001))
3766	out.append('Vector no.: %d Magnitude: %s\n'%(i,s))
3767	else:
3768	out.append('Vector no.: %d Magnitude: %8.4f Refine? %s\n'%
3769	(i,RBModel['VectMag'][i],RBModel['VectRef'][i]))
3770	out.append(' No. Type vx vy vz\n')
3771	for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
3772	out.append(' %d %2s %10.4f %10.4f %10.4f\n'%(j,name,xyz[0],xyz[1],xyz[2]))
3773	out.append(' No. Type x y z\n')
3774	for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
3775	out.append(' %d %2s %10.4f %10.4f %10.4f\n'%(i,name,xyz[0],xyz[1],xyz[2]))
3776	return out

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

Download in other formats: