Context Navigation

source: trunk/GSASIIstrIO.py @ 3760

Last change on this file since 3760 was 3760, checked in by vondreele, 4 years ago
remove AwaveType?,PwaveType?,FwaveType? & MwaveType? from constraints lists. fix for nonFourier refine results listing
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 150.5 KB

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

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

Download in other formats: