Context Navigation

source: trunk/GSASIIstrIO.py @ 3163

Last change on this file since 3163 was 3163, checked in by vondreele, 6 years ago
fix typo in Tb neutron scattering factors; extra space fix lst output of neutron resonant form factors - computations were correct implement import of mcif files.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 148.4 KB

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

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

Download in other formats: