Context Navigation

source: trunk/GSASIIstrIO.py @ 5037

Last change on this file since 5037 was 5025, checked in by vondreele, 4 years ago
Add shell plots to 3D reflection plots. Fix escape character issues in CopyRietveld2Origin Revise LeBail? operation - now works in sequential refinements. Remove "Fit LeBail?" from Calculate menu. OnLeBail? does 10 cycles instead of 3 in the "zero cycle" step. Initial FOSQ set to Scale*Phase fr. instead of 1.0 - improves initial LeBaill? fit. Put newLeBail flag in Controls - as more global. Allow LeBail? refinements in sequential refinements. Remove a wx.Yield() - obsolete Disable deletion of selected parameters in a LeBail? refinement & let the LS handle any singularities that may ensue Derivatives no longer controlled by LeBail? flag.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 176.8 KB

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

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

Download in other formats: