Context Navigation

source: trunk/GSASIIstrIO.py @ 3438

Last change on this file since 3438 was 3438, checked in by vondreele, 5 years ago
g2strIO fix magnetic hkl extinctions in reflection generation remove Uniq & Phi - never used change all config 'debug' out put to have 'DBG' in message - many places fix phoenix decode problem for Comments Display transformation math on Transform dialog box & labels on matrix & vectors use a general exception in getSelection (G2plot line 7141) to avoid crashes
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 150.1 KB

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

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

Download in other formats: