Context Navigation

source: trunk/GSASIIstrIO.py @ 3789

Last change on this file since 3789 was 3789, checked in by toby, 4 years ago
fix transpose error for lattice constraints
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 150.6 KB

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

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

Download in other formats: