Context Navigation

source: trunk/GSASIIstrIO.py @ 3711

Last change on this file since 3711 was 3711, checked in by toby, 5 years ago
major constraint update: move conflicting equivs to constraints; allow a formula for multiplier; update docs extensively. New routine EvaluateMultipliers? needed to evaluate formulae
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 150.2 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2018-10-29 21:10:08 +0000 (Mon, 29 Oct 2018) $
4	# $Author: toby $
5	# $Revision: 3711 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 3711 2018-10-29 21:10:08Z toby $
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: 3711 $")
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(constDict,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
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):
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 resetRefList and (not hapDict[pfx+'LeBail'] or (hapData['LeBail'] and hapData['newLeBail'])):
2356	if hapData.get('LeBail',True): #stop regeneating reflections for LeBail
2357	hapData['newLeBail'] = False
2358	refList = []
2359	# Uniq = []
2360	# Phi = []
2361	useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2362	if Phases[phase]['General'].get('Modulated',False):
2363	ifSuper = True
2364	HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2365	HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2366	for h,k,l,m,d in HKLd:
2367	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2368	mul *= 2 # for powder overlap of Friedel pairs
2369	if m or not ext or useExt:
2370	if 'C' in inst['Type'][0]:
2371	pos = G2lat.Dsp2pos(inst,d)
2372	if limits[0] < pos < limits[1]:
2373	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])
2374	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2375	# Uniq.append(uniq)
2376	# Phi.append(phi)
2377	elif 'T' in inst['Type'][0]:
2378	pos = G2lat.Dsp2pos(inst,d)
2379	if limits[0] < pos < limits[1]:
2380	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2381	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])
2382	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2383	#TODO - if tabulated put alp & bet in here
2384	# Uniq.append(uniq)
2385	# Phi.append(phi)
2386	else:
2387	ifSuper = False
2388	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2389	HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2390	for h,k,l,d in HKLd:
2391	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2392	if ext and 'N' in inst['Type'][0] and 'MagSpGrp' in SGData:
2393	ext = G2spc.checkMagextc([h,k,l],SGData)
2394	mul *= 2 # for powder overlap of Friedel pairs
2395	if ext:# and not useExt:
2396	continue
2397	if 'C' in inst['Type'][0]:
2398	pos = G2lat.Dsp2pos(inst,d)
2399	if limits[0] < pos < limits[1]:
2400	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])
2401	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2402	# Uniq.append(uniq)
2403	# Phi.append(phi)
2404	elif 'T' in inst['Type'][0]:
2405	pos = G2lat.Dsp2pos(inst,d)
2406	if limits[0] < pos < limits[1]:
2407	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2408	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])
2409	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2410	# Uniq.append(uniq)
2411	# Phi.append(phi)
2412	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2413	elif 'HKLF' in histogram:
2414	inst = Histogram['Instrument Parameters'][0]
2415	hId = Histogram['hId']
2416	hfx = ':%d:'%(hId)
2417	for item in inst:
2418	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2419	hapDict[hfx+item] = inst[item][1]
2420	pfx = str(pId)+':'+str(hId)+':'
2421	hapDict[pfx+'Scale'] = hapData['Scale'][0]
2422	if hapData['Scale'][1]:
2423	hapVary.append(pfx+'Scale')
2424
2425	extApprox,extType,extParms = hapData['Extinction']
2426	controlDict[pfx+'EType'] = extType
2427	controlDict[pfx+'EApprox'] = extApprox
2428	if 'C' in inst['Type'][0]:
2429	controlDict[pfx+'Tbar'] = extParms['Tbar']
2430	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2431	if 'Primary' in extType:
2432	Ekey = ['Ep',]
2433	elif 'I & II' in extType:
2434	Ekey = ['Eg','Es']
2435	elif 'Secondary Type II' == extType:
2436	Ekey = ['Es',]
2437	elif 'Secondary Type I' == extType:
2438	Ekey = ['Eg',]
2439	else: #'None'
2440	Ekey = []
2441	for eKey in Ekey:
2442	hapDict[pfx+eKey] = extParms[eKey][0]
2443	if extParms[eKey][1]:
2444	hapVary.append(pfx+eKey)
2445	for bab in ['BabA','BabU']:
2446	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2447	if hapData['Babinet'][bab][1]:
2448	hapVary.append(pfx+bab)
2449	Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2450	if len(Twins) == 1:
2451	hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2452	if hapData.get('Flack',[0,False])[1]:
2453	hapVary.append(pfx+'Flack')
2454	sumTwFr = 0.
2455	controlDict[pfx+'TwinLaw'] = []
2456	controlDict[pfx+'TwinInv'] = []
2457	NTL = 0
2458	for it,twin in enumerate(Twins):
2459	if 'bool' in str(type(twin[0])):
2460	controlDict[pfx+'TwinInv'].append(twin[0])
2461	controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2462	else:
2463	NTL += 1
2464	controlDict[pfx+'TwinInv'].append(False)
2465	controlDict[pfx+'TwinLaw'].append(twin[0])
2466	if it:
2467	hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2468	sumTwFr += twin[1]
2469	else:
2470	hapDict[pfx+'TwinFr:0'] = twin[1][0]
2471	controlDict[pfx+'TwinNMN'] = twin[1][1]
2472	if Twins[0][1][1]:
2473	hapVary.append(pfx+'TwinFr:'+str(it))
2474	controlDict[pfx+'NTL'] = NTL
2475	#need to make constraint on TwinFr
2476	controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2477	if len(Twins) > 1: #force sum to unity
2478	hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2479	if Print:
2480	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2481	pFile.write(135*'='+'\n')
2482	pFile.write(' Scale factor : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2483	if extType != 'None':
2484	pFile.write(' Extinction Type: %15s approx: %10s\n'%(extType,extApprox))
2485	text = ' Parameters :'
2486	for eKey in Ekey:
2487	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2488	pFile.write(text+'\n')
2489	if hapData['Babinet']['BabA'][0]:
2490	PrintBabinet(hapData['Babinet'])
2491	if not SGData['SGInv'] and len(Twins) == 1:
2492	pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2493	if len(Twins) > 1:
2494	for it,twin in enumerate(Twins):
2495	if 'bool' in str(type(twin[0])):
2496	pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine?\n'%
2497	(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1])
2498	else:
2499	pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2500	(str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],Twins[0][1][1]))
2501
2502	Histogram['Reflection Lists'] = phase
2503
2504	return hapVary,hapDict,controlDict
2505
2506	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2507	'needs a doc string'
2508
2509	def PrintSizeAndSig(hapData,sizeSig):
2510	line = '\n Size model: %9s'%(hapData[0])
2511	refine = False
2512	if hapData[0] in ['isotropic','uniaxial']:
2513	line += ' equatorial:%12.5f'%(hapData[1][0])
2514	if sizeSig[0][0]:
2515	line += ', sig:%8.4f'%(sizeSig[0][0])
2516	refine = True
2517	if hapData[0] == 'uniaxial':
2518	line += ' axial:%12.4f'%(hapData[1][1])
2519	if sizeSig[0][1]:
2520	refine = True
2521	line += ', sig:%8.4f'%(sizeSig[0][1])
2522	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2523	if sizeSig[0][2]:
2524	refine = True
2525	line += ', sig:%8.4f'%(sizeSig[0][2])
2526	if refine:
2527	pFile.write(line+'\n')
2528	else:
2529	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2530	if sizeSig[0][2]:
2531	refine = True
2532	line += ', sig:%8.4f'%(sizeSig[0][2])
2533	Snames = ['S11','S22','S33','S12','S13','S23']
2534	ptlbls = ' name :'
2535	ptstr = ' value :'
2536	sigstr = ' sig :'
2537	for i,name in enumerate(Snames):
2538	ptlbls += '%12s' % (name)
2539	ptstr += '%12.6f' % (hapData[4][i])
2540	if sizeSig[1][i]:
2541	refine = True
2542	sigstr += '%12.6f' % (sizeSig[1][i])
2543	else:
2544	sigstr += 12*' '
2545	if refine:
2546	pFile.write(line+'\n')
2547	pFile.write(ptlbls+'\n')
2548	pFile.write(ptstr+'\n')
2549	pFile.write(sigstr+'\n')
2550
2551	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2552	line = '\n Mustrain model: %9s\n'%(hapData[0])
2553	refine = False
2554	if hapData[0] in ['isotropic','uniaxial']:
2555	line += ' equatorial:%12.1f'%(hapData[1][0])
2556	if mustrainSig[0][0]:
2557	line += ', sig:%8.1f'%(mustrainSig[0][0])
2558	refine = True
2559	if hapData[0] == 'uniaxial':
2560	line += ' axial:%12.1f'%(hapData[1][1])
2561	if mustrainSig[0][1]:
2562	line += ', sig:%8.1f'%(mustrainSig[0][1])
2563	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2564	if mustrainSig[0][2]:
2565	refine = True
2566	line += ', sig:%8.3f'%(mustrainSig[0][2])
2567	if refine:
2568	pFile.write(line+'\n')
2569	else:
2570	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2571	if mustrainSig[0][2]:
2572	refine = True
2573	line += ', sig:%8.3f'%(mustrainSig[0][2])
2574	Snames = G2spc.MustrainNames(SGData)
2575	ptlbls = ' name :'
2576	ptstr = ' value :'
2577	sigstr = ' sig :'
2578	for i,name in enumerate(Snames):
2579	ptlbls += '%12s' % (name)
2580	ptstr += '%12.1f' % (hapData[4][i])
2581	if mustrainSig[1][i]:
2582	refine = True
2583	sigstr += '%12.1f' % (mustrainSig[1][i])
2584	else:
2585	sigstr += 12*' '
2586	if refine:
2587	pFile.write(line+'\n')
2588	pFile.write(ptlbls+'\n')
2589	pFile.write(ptstr+'\n')
2590	pFile.write(sigstr+'\n')
2591
2592	def PrintHStrainAndSig(hapData,strainSig,SGData):
2593	Hsnames = G2spc.HStrainNames(SGData)
2594	ptlbls = ' name :'
2595	ptstr = ' value :'
2596	sigstr = ' sig :'
2597	refine = False
2598	for i,name in enumerate(Hsnames):
2599	ptlbls += '%12s' % (name)
2600	ptstr += '%12.4g' % (hapData[0][i])
2601	if name in strainSig:
2602	refine = True
2603	sigstr += '%12.4g' % (strainSig[name])
2604	else:
2605	sigstr += 12*' '
2606	if refine:
2607	pFile.write('\n Hydrostatic/elastic strain:\n')
2608	pFile.write(ptlbls+'\n')
2609	pFile.write(ptstr+'\n')
2610	pFile.write(sigstr+'\n')
2611
2612	def PrintSHPOAndSig(pfx,hapData,POsig):
2613	pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2614	ptlbls = ' names :'
2615	ptstr = ' values:'
2616	sigstr = ' sig :'
2617	for item in hapData[5]:
2618	ptlbls += '%12s'%(item)
2619	ptstr += '%12.3f'%(hapData[5][item])
2620	if pfx+item in POsig:
2621	sigstr += '%12.3f'%(POsig[pfx+item])
2622	else:
2623	sigstr += 12*' '
2624	pFile.write(ptlbls+'\n')
2625	pFile.write(ptstr+'\n')
2626	pFile.write(sigstr+'\n')
2627
2628	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2629	pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2630	text = ''
2631	for item in hapData[2]:
2632	if pfx+item in ScalExtSig:
2633	text += ' %s: '%(item)
2634	text += '%12.2e'%(hapData[2][item][0])
2635	if pfx+item in ScalExtSig:
2636	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2637	pFile.write(text+'\n')
2638
2639	def PrintBabinetAndSig(pfx,hapData,BabSig):
2640	pFile.write('\n Babinet form factor modification:\n')
2641	ptlbls = ' names :'
2642	ptstr = ' values:'
2643	sigstr = ' sig :'
2644	for item in hapData:
2645	ptlbls += '%12s'%(item)
2646	ptstr += '%12.3f'%(hapData[item][0])
2647	if pfx+item in BabSig:
2648	sigstr += '%12.3f'%(BabSig[pfx+item])
2649	else:
2650	sigstr += 12*' '
2651	pFile.write(ptlbls+'\n')
2652	pFile.write(ptstr+'\n')
2653	pFile.write(sigstr+'\n')
2654
2655	def PrintTwinsAndSig(pfx,twinData,TwinSig):
2656	pFile.write('\n Twin Law fractions :\n')
2657	ptlbls = ' names :'
2658	ptstr = ' values:'
2659	sigstr = ' sig :'
2660	for it,item in enumerate(twinData):
2661	ptlbls += ' twin #%d'%(it)
2662	if it:
2663	ptstr += '%12.3f'%(item[1])
2664	else:
2665	ptstr += '%12.3f'%(item[1][0])
2666	if pfx+'TwinFr:'+str(it) in TwinSig:
2667	sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2668	else:
2669	sigstr += 12*' '
2670	pFile.write(ptlbls+'\n')
2671	pFile.write(ptstr+'\n')
2672	pFile.write(sigstr+'\n')
2673
2674
2675	PhFrExtPOSig = {}
2676	SizeMuStrSig = {}
2677	ScalExtSig = {}
2678	BabSig = {}
2679	TwinFrSig = {}
2680	wtFrSum = {}
2681	for phase in Phases:
2682	HistoPhase = Phases[phase]['Histograms']
2683	General = Phases[phase]['General']
2684	SGData = General['SGData']
2685	pId = Phases[phase]['pId']
2686	histoList = list(HistoPhase.keys())
2687	histoList.sort()
2688	for histogram in histoList:
2689	try:
2690	Histogram = Histograms[histogram]
2691	except KeyError:
2692	#skip if histogram not included e.g. in a sequential refinement
2693	continue
2694	if not Phases[phase]['Histograms'][histogram]['Use']:
2695	#skip if phase absent from this histogram
2696	continue
2697	hapData = HistoPhase[histogram]
2698	hId = Histogram['hId']
2699	pfx = str(pId)+':'+str(hId)+':'
2700	if hId not in wtFrSum:
2701	wtFrSum[hId] = 0.
2702	if 'PWDR' in histogram:
2703	for item in ['Scale','Extinction']:
2704	hapData[item][0] = parmDict[pfx+item]
2705	if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2706	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2707	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2708	if hapData['Pref.Ori.'][0] == 'MD':
2709	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2710	if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2711	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2712	else: #'SH' spherical harmonics
2713	for item in hapData['Pref.Ori.'][5]:
2714	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2715	if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2716	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2717	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2718	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2719	pfx+'HStrain':{}})
2720	for item in ['Mustrain','Size']:
2721	hapData[item][1][2] = parmDict[pfx+item+';mx']
2722	# hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2723	if pfx+item+';mx' in sigDict:
2724	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2725	if hapData[item][0] in ['isotropic','uniaxial']:
2726	hapData[item][1][0] = parmDict[pfx+item+';i']
2727	if item == 'Size':
2728	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2729	if pfx+item+';i' in sigDict:
2730	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2731	if hapData[item][0] == 'uniaxial':
2732	hapData[item][1][1] = parmDict[pfx+item+';a']
2733	if item == 'Size':
2734	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2735	if pfx+item+';a' in sigDict:
2736	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2737	else: #generalized for mustrain or ellipsoidal for size
2738	Nterms = len(hapData[item][4])
2739	for i in range(Nterms):
2740	sfx = ';'+str(i)
2741	hapData[item][4][i] = parmDict[pfx+item+sfx]
2742	if pfx+item+sfx in sigDict:
2743	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2744	names = G2spc.HStrainNames(SGData)
2745	for i,name in enumerate(names):
2746	hapData['HStrain'][0][i] = parmDict[pfx+name]
2747	if pfx+name in sigDict:
2748	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2749	if Phases[phase]['General']['Type'] != 'magnetic':
2750	for name in ['BabA','BabU']:
2751	hapData['Babinet'][name][0] = parmDict[pfx+name]
2752	if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
2753	BabSig[pfx+name] = sigDict[pfx+name]
2754
2755	elif 'HKLF' in histogram:
2756	for item in ['Scale','Flack']:
2757	if parmDict.get(pfx+item):
2758	hapData[item][0] = parmDict[pfx+item]
2759	if pfx+item in sigDict:
2760	ScalExtSig[pfx+item] = sigDict[pfx+item]
2761	for item in ['Ep','Eg','Es']:
2762	if parmDict.get(pfx+item):
2763	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2764	if pfx+item in sigDict:
2765	ScalExtSig[pfx+item] = sigDict[pfx+item]
2766	for item in ['BabA','BabU']:
2767	hapData['Babinet'][item][0] = parmDict[pfx+item]
2768	if pfx+item in sigDict:
2769	BabSig[pfx+item] = sigDict[pfx+item]
2770	if 'Twins' in hapData:
2771	it = 1
2772	sumTwFr = 0.
2773	while True:
2774	try:
2775	hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2776	if pfx+'TwinFr:'+str(it) in sigDict:
2777	TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2778	if it:
2779	sumTwFr += hapData['Twins'][it][1]
2780	it += 1
2781	except KeyError:
2782	break
2783	hapData['Twins'][0][1][0] = 1.-sumTwFr
2784
2785	if Print:
2786	for phase in Phases:
2787	HistoPhase = Phases[phase]['Histograms']
2788	General = Phases[phase]['General']
2789	SGData = General['SGData']
2790	pId = Phases[phase]['pId']
2791	histoList = list(HistoPhase.keys())
2792	histoList.sort()
2793	for histogram in histoList:
2794	try:
2795	Histogram = Histograms[histogram]
2796	except KeyError:
2797	#skip if histogram not included e.g. in a sequential refinement
2798	continue
2799	hapData = HistoPhase[histogram]
2800	hId = Histogram['hId']
2801	Histogram['Residuals'][str(pId)+'::Name'] = phase
2802	pfx = str(pId)+':'+str(hId)+':'
2803	hfx = ':%s:'%(hId)
2804	if pfx+'Nref' not in Histogram['Residuals']: #skip not used phase in histogram
2805	continue
2806	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2807	pFile.write(135*'='+'\n')
2808	if 'PWDR' in histogram:
2809	pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
2810	(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
2811	pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
2812	pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
2813
2814	if parmDict[pfx+'LeBail']:
2815	pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
2816	else:
2817	if pfx+'Scale' in PhFrExtPOSig:
2818	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2819	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2820	pFile.write(' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f\n'%
2821	(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
2822	if pfx+'Extinction' in PhFrExtPOSig:
2823	pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
2824	if hapData['Pref.Ori.'][0] == 'MD':
2825	if pfx+'MD' in PhFrExtPOSig:
2826	pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
2827	else:
2828	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2829	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2830	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2831	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2832	if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
2833	if len(BabSig):
2834	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2835
2836	elif 'HKLF' in histogram:
2837	Inst = Histogram['Instrument Parameters'][0]
2838	pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
2839	(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2840	Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
2841	if FFtables != None and 'N' not in Inst['Type'][0]:
2842	PrintFprime(FFtables,hfx,pFile)
2843	pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
2844	if pfx+'Scale' in ScalExtSig:
2845	pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
2846	if hapData['Extinction'][0] != 'None':
2847	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2848	if len(BabSig):
2849	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2850	if pfx+'Flack' in ScalExtSig:
2851	pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
2852	if len(TwinFrSig):
2853	PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2854
2855	################################################################################
2856	##### Histogram data
2857	################################################################################
2858
2859	def GetHistogramData(Histograms,Print=True,pFile=None):
2860	'needs a doc string'
2861
2862	def GetBackgroundParms(hId,Background):
2863	Back = Background[0]
2864	DebyePeaks = Background[1]
2865	bakType,bakFlag = Back[:2]
2866	backVals = Back[3:]
2867	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2868	backDict = dict(zip(backNames,backVals))
2869	backVary = []
2870	if bakFlag:
2871	backVary = backNames
2872	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2873	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2874	debyeDict = {}
2875	debyeList = []
2876	for i in range(DebyePeaks['nDebye']):
2877	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2878	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2879	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2880	debyeVary = []
2881	for item in debyeList:
2882	if item[1]:
2883	debyeVary.append(item[0])
2884	backDict.update(debyeDict)
2885	backVary += debyeVary
2886	peakDict = {}
2887	peakList = []
2888	for i in range(DebyePeaks['nPeaks']):
2889	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2890	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2891	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2892	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2893	peakVary = []
2894	for item in peakList:
2895	if item[1]:
2896	peakVary.append(item[0])
2897	backDict.update(peakDict)
2898	backVary += peakVary
2899	return bakType,backDict,backVary
2900
2901	def GetInstParms(hId,Inst):
2902	#patch
2903	if 'Z' not in Inst:
2904	Inst['Z'] = [0.0,0.0,False]
2905	dataType = Inst['Type'][0]
2906	instDict = {}
2907	insVary = []
2908	pfx = ':'+str(hId)+':'
2909	insKeys = list(Inst.keys())
2910	insKeys.sort()
2911	for item in insKeys:
2912	insName = pfx+item
2913	instDict[insName] = Inst[item][1]
2914	if len(Inst[item]) > 2 and Inst[item][2]:
2915	insVary.append(insName)
2916	if 'C' in dataType:
2917	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2918	elif 'T' in dataType: #trap zero alp, bet coeff.
2919	if not instDict[pfx+'alpha']:
2920	instDict[pfx+'alpha'] = 1.0
2921	if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
2922	instDict[pfx+'beta-1'] = 1.0
2923	return dataType,instDict,insVary
2924
2925	def GetSampleParms(hId,Sample):
2926	sampVary = []
2927	hfx = ':'+str(hId)+':'
2928	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2929	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2930	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2931	if key in Sample:
2932	sampDict[hfx+key] = Sample[key]
2933	Type = Sample['Type']
2934	if 'Bragg' in Type: #Bragg-Brentano
2935	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2936	sampDict[hfx+item] = Sample[item][0]
2937	if Sample[item][1]:
2938	sampVary.append(hfx+item)
2939	elif 'Debye' in Type: #Debye-Scherrer
2940	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2941	sampDict[hfx+item] = Sample[item][0]
2942	if Sample[item][1]:
2943	sampVary.append(hfx+item)
2944	return Type,sampDict,sampVary
2945
2946	def PrintBackground(Background):
2947	Back = Background[0]
2948	DebyePeaks = Background[1]
2949	pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
2950	line = ' Coefficients: '
2951	for i,back in enumerate(Back[3:]):
2952	line += '%10.3f'%(back)
2953	if i and not i%10:
2954	line += '\n'+15*' '
2955	pFile.write(line+'\n')
2956	if DebyePeaks['nDebye']:
2957	pFile.write('\n Debye diffuse scattering coefficients\n')
2958	parms = ['DebyeA','DebyeR','DebyeU']
2959	line = ' names : '
2960	for parm in parms:
2961	line += '%8s refine?'%(parm)
2962	pFile.write(line+'\n')
2963	for j,term in enumerate(DebyePeaks['debyeTerms']):
2964	line = ' term'+'%2d'%(j)+':'
2965	for i in range(3):
2966	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2967	pFile.write(line+'\n')
2968	if DebyePeaks['nPeaks']:
2969	pFile.write('\n Single peak coefficients\n')
2970	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2971	line = ' names : '
2972	for parm in parms:
2973	line += '%8s refine?'%(parm)
2974	pFile.write(line+'\n')
2975	for j,term in enumerate(DebyePeaks['peaksList']):
2976	line = ' peak'+'%2d'%(j)+':'
2977	for i in range(4):
2978	line += '%12.3f %5s'%(term[2i],bool(term[2i+1]))
2979	pFile.write(line+'\n')
2980
2981	def PrintInstParms(Inst):
2982	pFile.write('\n Instrument Parameters:\n')
2983	insKeys = list(Inst.keys())
2984	insKeys.sort()
2985	iBeg = 0
2986	Ok = True
2987	while Ok:
2988	ptlbls = ' name :'
2989	ptstr = ' value :'
2990	varstr = ' refine:'
2991	iFin = min(iBeg+9,len(insKeys))
2992	for item in insKeys[iBeg:iFin]:
2993	if item not in ['Type','Source']:
2994	ptlbls += '%12s' % (item)
2995	ptstr += '%12.6f' % (Inst[item][1])
2996	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2997	varstr += 12*' '
2998	else:
2999	varstr += '%12s' % (str(bool(Inst[item][2])))
3000	pFile.write(ptlbls+'\n')
3001	pFile.write(ptstr+'\n')
3002	pFile.write(varstr+'\n')
3003	iBeg = iFin
3004	if iBeg == len(insKeys):
3005	Ok = False
3006	else:
3007	pFile.write('\n')
3008
3009	def PrintSampleParms(Sample):
3010	pFile.write('\n Sample Parameters:\n')
3011	pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3012	(Sample['Omega'],Sample['Chi'],Sample['Phi']))
3013	ptlbls = ' name :'
3014	ptstr = ' value :'
3015	varstr = ' refine:'
3016	if 'Bragg' in Sample['Type']:
3017	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3018	ptlbls += '%14s'%(item)
3019	ptstr += '%14.4f'%(Sample[item][0])
3020	varstr += '%14s'%(str(bool(Sample[item][1])))
3021
3022	elif 'Debye' in Type: #Debye-Scherrer
3023	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3024	ptlbls += '%14s'%(item)
3025	ptstr += '%14.4f'%(Sample[item][0])
3026	varstr += '%14s'%(str(bool(Sample[item][1])))
3027
3028	pFile.write(ptlbls+'\n')
3029	pFile.write(ptstr+'\n')
3030	pFile.write(varstr+'\n')
3031
3032	histDict = {}
3033	histVary = []
3034	controlDict = {}
3035	histoList = list(Histograms.keys())
3036	histoList.sort()
3037	for histogram in histoList:
3038	if 'PWDR' in histogram:
3039	Histogram = Histograms[histogram]
3040	hId = Histogram['hId']
3041	pfx = ':'+str(hId)+':'
3042	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3043	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3044	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3045	for excl in controlDict[pfx+'Exclude']:
3046	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3047	if controlDict[pfx+'Exclude']:
3048	ma.mask_rows(Histogram['Data'])
3049	Background = Histogram['Background']
3050	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3051	controlDict[pfx+'bakType'] = Type
3052	histDict.update(bakDict)
3053	histVary += bakVary
3054
3055	Inst = Histogram['Instrument Parameters'] #TODO ? ignores tabulated alp,bet & delt for TOF
3056	if 'T' in Type and len(Inst[1]): #patch - back-to-back exponential contribution to TOF line shape is removed
3057	print ('Warning: tabulated profile coefficients are ignored')
3058	Type,instDict,insVary = GetInstParms(hId,Inst[0])
3059	controlDict[pfx+'histType'] = Type
3060	if 'XC' in Type:
3061	if pfx+'Lam1' in instDict:
3062	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3063	else:
3064	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
3065	histDict.update(instDict)
3066	histVary += insVary
3067
3068	Sample = Histogram['Sample Parameters']
3069	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3070	controlDict[pfx+'instType'] = Type
3071	histDict.update(sampDict)
3072	histVary += sampVary
3073
3074
3075	if Print:
3076	pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3077	pFile.write(135*'='+'\n')
3078	Units = {'C':' deg','T':' msec'}
3079	units = Units[controlDict[pfx+'histType'][2]]
3080	Limits = controlDict[pfx+'Limits']
3081	pFile.write(' Instrument type: %s\n'%Sample['Type'])
3082	pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3083	if len(controlDict[pfx+'Exclude']):
3084	excls = controlDict[pfx+'Exclude']
3085	for excl in excls:
3086	pFile.write(' Excluded region: %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units))
3087	PrintSampleParms(Sample)
3088	PrintInstParms(Inst[0])
3089	PrintBackground(Background)
3090	elif 'HKLF' in histogram:
3091	Histogram = Histograms[histogram]
3092	hId = Histogram['hId']
3093	pfx = ':'+str(hId)+':'
3094	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3095	Inst = Histogram['Instrument Parameters'][0]
3096	controlDict[pfx+'histType'] = Inst['Type'][0]
3097	if 'X' in Inst['Type'][0]:
3098	histDict[pfx+'Lam'] = Inst['Lam'][1]
3099	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3100	return histVary,histDict,controlDict
3101
3102	def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3103	'needs a doc string'
3104
3105	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3106	Back = Background[0]
3107	DebyePeaks = Background[1]
3108	lenBack = len(Back[3:])
3109	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
3110	for i in range(lenBack):
3111	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3112	if pfx+'Back;'+str(i) in sigDict:
3113	backSig[i] = sigDict[pfx+'Back;'+str(i)]
3114	if DebyePeaks['nDebye']:
3115	for i in range(DebyePeaks['nDebye']):
3116	names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3117	for j,name in enumerate(names):
3118	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3119	if name in sigDict:
3120	backSig[lenBack+3*i+j] = sigDict[name]
3121	if DebyePeaks['nPeaks']:
3122	for i in range(DebyePeaks['nPeaks']):
3123	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3124	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3125	for j,name in enumerate(names):
3126	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3127	if name in sigDict:
3128	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
3129	return backSig
3130
3131	def SetInstParms(pfx,Inst,parmDict,sigDict):
3132	instSig = {}
3133	insKeys = list(Inst.keys())
3134	insKeys.sort()
3135	for item in insKeys:
3136	insName = pfx+item
3137	Inst[item][1] = parmDict[insName]
3138	if insName in sigDict:
3139	instSig[item] = sigDict[insName]
3140	else:
3141	instSig[item] = 0
3142	return instSig
3143
3144	def SetSampleParms(pfx,Sample,parmDict,sigDict):
3145	if 'Bragg' in Sample['Type']: #Bragg-Brentano
3146	sampSig = [0 for i in range(5)]
3147	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3148	Sample[item][0] = parmDict[pfx+item]
3149	if pfx+item in sigDict:
3150	sampSig[i] = sigDict[pfx+item]
3151	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3152	sampSig = [0 for i in range(4)]
3153	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3154	Sample[item][0] = parmDict[pfx+item]
3155	if pfx+item in sigDict:
3156	sampSig[i] = sigDict[pfx+item]
3157	return sampSig
3158
3159	def PrintBackgroundSig(Background,backSig):
3160	Back = Background[0]
3161	DebyePeaks = Background[1]
3162	valstr = ' value : '
3163	sigstr = ' sig : '
3164	refine = False
3165	for i,back in enumerate(Back[3:]):
3166	valstr += '%10.4g'%(back)
3167	if Back[1]:
3168	refine = True
3169	sigstr += '%10.4g'%(backSig[i])
3170	else:
3171	sigstr += 10*' '
3172	if refine:
3173	pFile.write('\n Background function: %s\n'%Back[0])
3174	pFile.write(valstr+'\n')
3175	pFile.write(sigstr+'\n')
3176	if DebyePeaks['nDebye']:
3177	ifAny = False
3178	ptfmt = "%12.3f"
3179	names = ' names :'
3180	ptstr = ' values:'
3181	sigstr = ' esds :'
3182	for item in sigDict:
3183	if 'Debye' in item:
3184	ifAny = True
3185	names += '%12s'%(item)
3186	ptstr += ptfmt%(parmDict[item])
3187	sigstr += ptfmt%(sigDict[item])
3188	if ifAny:
3189	pFile.write('\n Debye diffuse scattering coefficients\n')
3190	pFile.write(names+'\n')
3191	pFile.write(ptstr+'\n')
3192	pFile.write(sigstr+'\n')
3193	if DebyePeaks['nPeaks']:
3194	pFile.write('\n Single peak coefficients:\n')
3195	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3196	line = ' peak no. '
3197	for parm in parms:
3198	line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3199	pFile.write(line+'\n')
3200	for ip in range(DebyePeaks['nPeaks']):
3201	ptstr = ' %4d '%(ip)
3202	for parm in parms:
3203	fmt = '%14.3f'
3204	efmt = '%12.3f'
3205	if parm == 'BkPkpos':
3206	fmt = '%14.4f'
3207	efmt = '%12.4f'
3208	name = pfx+parm+';%d'%(ip)
3209	ptstr += fmt%(parmDict[name])
3210	if name in sigDict:
3211	ptstr += efmt%(sigDict[name])
3212	else:
3213	ptstr += 12*' '
3214	pFile.write(ptstr+'\n')
3215	sumBk = np.array(Histogram['sumBk'])
3216	pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3217	(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3218
3219	def PrintInstParmsSig(Inst,instSig):
3220	refine = False
3221	insKeys = list(instSig.keys())
3222	insKeys.sort()
3223	iBeg = 0
3224	Ok = True
3225	while Ok:
3226	ptlbls = ' names :'
3227	ptstr = ' value :'
3228	sigstr = ' sig :'
3229	iFin = min(iBeg+9,len(insKeys))
3230	for name in insKeys[iBeg:iFin]:
3231	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3232	ptlbls += '%12s' % (name)
3233	ptstr += '%12.6f' % (Inst[name][1])
3234	if instSig[name]:
3235	refine = True
3236	sigstr += '%12.6f' % (instSig[name])
3237	else:
3238	sigstr += 12*' '
3239	if refine:
3240	pFile.write('\n Instrument Parameters:\n')
3241	pFile.write(ptlbls+'\n')
3242	pFile.write(ptstr+'\n')
3243	pFile.write(sigstr+'\n')
3244	iBeg = iFin
3245	if iBeg == len(insKeys):
3246	Ok = False
3247
3248	def PrintSampleParmsSig(Sample,sampleSig):
3249	ptlbls = ' names :'
3250	ptstr = ' values:'
3251	sigstr = ' sig :'
3252	refine = False
3253	if 'Bragg' in Sample['Type']:
3254	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3255	ptlbls += '%14s'%(item)
3256	ptstr += '%14.4f'%(Sample[item][0])
3257	if sampleSig[i]:
3258	refine = True
3259	sigstr += '%14.4f'%(sampleSig[i])
3260	else:
3261	sigstr += 14*' '
3262
3263	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3264	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3265	ptlbls += '%14s'%(item)
3266	ptstr += '%14.4f'%(Sample[item][0])
3267	if sampleSig[i]:
3268	refine = True
3269	sigstr += '%14.4f'%(sampleSig[i])
3270	else:
3271	sigstr += 14*' '
3272
3273	if refine:
3274	pFile.write('\n Sample Parameters:\n')
3275	pFile.write(ptlbls+'\n')
3276	pFile.write(ptstr+'\n')
3277	pFile.write(sigstr+'\n')
3278
3279	histoList = list(Histograms.keys())
3280	histoList.sort()
3281	for histogram in histoList:
3282	if 'PWDR' in histogram:
3283	Histogram = Histograms[histogram]
3284	hId = Histogram['hId']
3285	pfx = ':'+str(hId)+':'
3286	Background = Histogram['Background']
3287	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3288
3289	Inst = Histogram['Instrument Parameters'][0]
3290	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3291
3292	Sample = Histogram['Sample Parameters']
3293	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3294
3295	pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3296	pFile.write(135*'='+'\n')
3297	pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3298	pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3299	(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3300	pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3301	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3302	if Print:
3303	pFile.write(' Instrument type: %s\n'%Sample['Type'])
3304	if FFtables != None and 'N' not in Inst['Type'][0]:
3305	PrintFprime(FFtables,pfx,pFile)
3306	PrintSampleParmsSig(Sample,sampSig)
3307	PrintInstParmsSig(Inst,instSig)
3308	PrintBackgroundSig(Background,backSig)
3309

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

Download in other formats: