Context Navigation

source: trunk/GSASIIstrIO.py @ 3883

Last change on this file since 3883 was 3883, checked in by toby, 4 years ago
correct seq ref hist error
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 158.7 KB

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

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

Download in other formats: