Context Navigation

source: trunk/GSASIIstrIO.py @ 3873

Last change on this file since 3873 was 3873, checked in by toby, 4 years ago
trap error if unable to delete seq files
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 158.6 KB

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

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

Download in other formats: