Context Navigation

source: trunk/GSASIIstrIO.py @ 4917

Last change on this file since 4917 was 4917, checked in by toby, 2 years ago
split G2 & pillow readers for TIF; document
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 177.3 KB

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

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

Download in other formats: