Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: trunk/GSASIIstrIO.py

Last change on this file was 5494, checked in by vondreele, 3 days ago
fixes to spin RB stuff - now will refine (no derivatives yet)
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 196.4 KB

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