Context Navigation

source: trunk/GSASIIstrIO.py @ 3377

Last change on this file since 3377 was 3377, checked in by vondreele, 5 years ago
revise sequential refinement to retain previous runs & update those refined with new values. provide a "Clear seq. results" button
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 149.9 KB

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

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

Download in other formats: