Context Navigation

source: trunk/GSASIIstrIO.py @ 3311

Last change on this file since 3311 was 3311, checked in by vondreele, 6 years ago
remove unused SytSym? argument from GetCSpqinel - not needed
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 149.0 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2018-03-06 15:43:03 +0000 (Tue, 06 Mar 2018) $
4	# $Author: vondreele $
5	# $Revision: 3311 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 3311 2018-03-06 15:43:03Z 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: 3311 $")
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,Waves[0]))
990	else:
991	continue
992	for iw,wave in enumerate(Waves[1:]):
993	line = ''
994	if Waves[0] in ['Block','ZigZag'] and Stype == 'Spos' and not iw:
995	for item in names[Stype][6:]:
996	line += '%8s '%(item)
997	else:
998	if Stype == 'Spos':
999	for item in names[Stype][:6]:
1000	line += '%8s '%(item)
1001	else:
1002	for item in names[Stype]:
1003	line += '%8s '%(item)
1004	pFile.write(line+'\n')
1005	line = ''
1006	for item in wave[0]:
1007	line += '%8.4f '%(item)
1008	line += ' Refine? '+str(wave[1])
1009	pFile.write(line+'\n')
1010
1011	def PrintTexture(textureData):
1012	topstr = '\n Spherical harmonics texture: Order:' + \
1013	str(textureData['Order'])
1014	if textureData['Order']:
1015	pFile.write('%s Refine? %s\n'%(topstr,textureData['SH Coeff'][0]))
1016	else:
1017	pFile.write(topstr+'\n')
1018	return
1019	names = ['omega','chi','phi']
1020	line = '\n'
1021	for name in names:
1022	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
1023	pFile.write(line+'\n')
1024	pFile.write('\n Texture coefficients:\n')
1025	SHcoeff = textureData['SH Coeff'][1]
1026	SHkeys = list(SHcoeff.keys())
1027	nCoeff = len(SHcoeff)
1028	nBlock = nCoeff//10+1
1029	iBeg = 0
1030	iFin = min(iBeg+10,nCoeff)
1031	for block in range(nBlock):
1032	ptlbls = ' names :'
1033	ptstr = ' values:'
1034	for item in SHkeys[iBeg:iFin]:
1035	ptlbls += '%12s'%(item)
1036	ptstr += '%12.4f'%(SHcoeff[item])
1037	pFile.write(ptlbls+'\n')
1038	pFile.write(ptstr+'\n')
1039	iBeg += 10
1040	iFin = min(iBeg+10,nCoeff)
1041
1042	def MakeRBParms(rbKey,phaseVary,phaseDict):
1043	rbid = str(rbids.index(RB['RBId']))
1044	pfxRB = pfx+'RB'+rbKey+'P'
1045	pstr = ['x','y','z']
1046	ostr = ['a','i','j','k']
1047	for i in range(3):
1048	name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1049	phaseDict[name] = RB['Orig'][0][i]
1050	if RB['Orig'][1]:
1051	phaseVary += [name,]
1052	pfxRB = pfx+'RB'+rbKey+'O'
1053	for i in range(4):
1054	name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1055	phaseDict[name] = RB['Orient'][0][i]
1056	if RB['Orient'][1] == 'AV' and i:
1057	phaseVary += [name,]
1058	elif RB['Orient'][1] == 'A' and not i:
1059	phaseVary += [name,]
1060
1061	def MakeRBThermals(rbKey,phaseVary,phaseDict):
1062	rbid = str(rbids.index(RB['RBId']))
1063	tlstr = ['11','22','33','12','13','23']
1064	sstr = ['12','13','21','23','31','32','AA','BB']
1065	if 'T' in RB['ThermalMotion'][0]:
1066	pfxRB = pfx+'RB'+rbKey+'T'
1067	for i in range(6):
1068	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1069	phaseDict[name] = RB['ThermalMotion'][1][i]
1070	if RB['ThermalMotion'][2][i]:
1071	phaseVary += [name,]
1072	if 'L' in RB['ThermalMotion'][0]:
1073	pfxRB = pfx+'RB'+rbKey+'L'
1074	for i in range(6):
1075	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1076	phaseDict[name] = RB['ThermalMotion'][1][i+6]
1077	if RB['ThermalMotion'][2][i+6]:
1078	phaseVary += [name,]
1079	if 'S' in RB['ThermalMotion'][0]:
1080	pfxRB = pfx+'RB'+rbKey+'S'
1081	for i in range(8):
1082	name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1083	phaseDict[name] = RB['ThermalMotion'][1][i+12]
1084	if RB['ThermalMotion'][2][i+12]:
1085	phaseVary += [name,]
1086	if 'U' in RB['ThermalMotion'][0]:
1087	name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1088	phaseDict[name] = RB['ThermalMotion'][1][0]
1089	if RB['ThermalMotion'][2][0]:
1090	phaseVary += [name,]
1091
1092	def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1093	rbid = str(rbids.index(RB['RBId']))
1094	pfxRB = pfx+'RB'+rbKey+'Tr;'
1095	for i,torsion in enumerate(RB['Torsions']):
1096	name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1097	phaseDict[name] = torsion[0]
1098	if torsion[1]:
1099	phaseVary += [name,]
1100
1101	if Print:
1102	pFile.write('\n Phases:\n')
1103	phaseVary = []
1104	phaseDict = {}
1105	pawleyLookup = {}
1106	FFtables = {} #scattering factors - xrays
1107	MFtables = {} #Mag. form factors
1108	BLtables = {} # neutrons
1109	Natoms = {}
1110	maxSSwave = {}
1111	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1112	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1113	atomIndx = {}
1114	for name in PhaseData:
1115	General = PhaseData[name]['General']
1116	pId = PhaseData[name]['pId']
1117	pfx = str(pId)+'::'
1118	FFtable = G2el.GetFFtable(General['AtomTypes'])
1119	BLtable = G2el.GetBLtable(General)
1120	FFtables.update(FFtable)
1121	BLtables.update(BLtable)
1122	phaseDict[pfx+'isMag'] = False
1123	SGData = General['SGData']
1124	SGtext,SGtable = G2spc.SGPrint(SGData)
1125	if General['Type'] == 'magnetic':
1126	MFtable = G2el.GetMFtable(General['AtomTypes'],General['Lande g'])
1127	MFtables.update(MFtable)
1128	phaseDict[pfx+'isMag'] = True
1129	SpnFlp = SGData['SpnFlp']
1130	Atoms = PhaseData[name]['Atoms']
1131	if Atoms and not General.get('doPawley'):
1132	cx,ct,cs,cia = General['AtomPtrs']
1133	AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1134	PawleyRef = PhaseData[name].get('Pawley ref',[])
1135	cell = General['Cell']
1136	A = G2lat.cell2A(cell[1:7])
1137	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1138	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1139	if cell[0]:
1140	phaseVary += cellVary(pfx,SGData) #also fills in symmetry required constraints
1141	SSGtext = [] #no superstructure
1142	im = 0
1143	if General.get('Modulated',False):
1144	im = 1 #refl offset
1145	Vec,vRef,maxH = General['SuperVec']
1146	phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
1147	SSGData = General['SSGData']
1148	SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
1149	if vRef:
1150	phaseVary += modVary(pfx,SSGData)
1151	resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1152	if resRBData:
1153	rbids = rbIds['Residue'] #NB: used in the MakeRB routines
1154	for iRB,RB in enumerate(resRBData):
1155	MakeRBParms('R',phaseVary,phaseDict)
1156	MakeRBThermals('R',phaseVary,phaseDict)
1157	MakeRBTorsions('R',phaseVary,phaseDict)
1158
1159	vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1160	if vecRBData:
1161	rbids = rbIds['Vector'] #NB: used in the MakeRB routines
1162	for iRB,RB in enumerate(vecRBData):
1163	MakeRBParms('V',phaseVary,phaseDict)
1164	MakeRBThermals('V',phaseVary,phaseDict)
1165
1166	Natoms[pfx] = 0
1167	maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0}
1168	if Atoms and not General.get('doPawley'):
1169	cx,ct,cs,cia = General['AtomPtrs']
1170	Natoms[pfx] = len(Atoms)
1171	for i,at in enumerate(Atoms):
1172	atomIndx[at[cia+8]] = [pfx,i] #lookup table for restraints
1173	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1174	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1175	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
1176	pfx+'AI/A:'+str(i):at[cia],})
1177	if at[cia] == 'I':
1178	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1179	else:
1180	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1181	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1182	if General['Type'] == 'magnetic':
1183	phaseDict.update({pfx+'AMx:'+str(i):at[cx+4],pfx+'AMy:'+str(i):at[cx+5],pfx+'AMz:'+str(i):at[cx+6]})
1184	if 'F' in at[ct+1]:
1185	phaseVary.append(pfx+'Afrac:'+str(i))
1186	if 'X' in at[ct+1]:
1187	try: #patch for sytsym name changes
1188	xId,xCoef = G2spc.GetCSxinel(at[cs])
1189	except KeyError:
1190	Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1191	at[cs] = Sytsym
1192	xId,xCoef = G2spc.GetCSxinel(at[cs])
1193	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1194	equivs = [[],[],[]]
1195	for j in range(3):
1196	if xId[j] > 0:
1197	phaseVary.append(names[j])
1198	equivs[xId[j]-1].append([names[j],xCoef[j]])
1199	for equiv in equivs:
1200	if len(equiv) > 1:
1201	name = equiv[0][0]
1202	coef = equiv[0][1]
1203	for eqv in equiv[1:]:
1204	eqv[1] /= coef
1205	G2mv.StoreEquivalence(name,(eqv,))
1206	if 'U' in at[ct+1]:
1207	if at[cia] == 'I':
1208	phaseVary.append(pfx+'AUiso:'+str(i))
1209	else:
1210	try: #patch for sytsym name changes
1211	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1212	except KeyError:
1213	Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1214	at[cs] = Sytsym
1215	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1216	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1217	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1218	equivs = [[],[],[],[],[],[]]
1219	for j in range(6):
1220	if uId[j] > 0:
1221	phaseVary.append(names[j])
1222	equivs[uId[j]-1].append([names[j],uCoef[j]])
1223	for equiv in equivs:
1224	if len(equiv) > 1:
1225	name = equiv[0][0]
1226	coef = equiv[0][1]
1227	for eqv in equiv[1:]:
1228	eqv[1] /= coef
1229	G2mv.StoreEquivalence(name,(eqv,))
1230	if 'M' in at[ct+1]:
1231	SytSym,Mul,Nop,dupDir = G2spc.SytSym(at[cx:cx+3],SGData)
1232	mId,mCoef = G2spc.GetCSpqinel(SpnFlp,dupDir)
1233	names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)]
1234	equivs = [[],[],[]]
1235	for j in range(3):
1236	if mId[j] > 0:
1237	phaseVary.append(names[j])
1238	equivs[mId[j]-1].append([names[j],mCoef[j]])
1239	for equiv in equivs:
1240	if len(equiv) > 1:
1241	name = equiv[0][0]
1242	coef = equiv[0][1]
1243	for eqv in equiv[1:]:
1244	eqv[1] /= coef
1245	G2mv.StoreEquivalence(name,(eqv,))
1246	if General.get('Modulated',False):
1247	AtomSS = at[-1]['SS1']
1248	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1249	Waves = AtomSS[Stype]
1250	if len(Waves):
1251	waveType = Waves[0]
1252	else:
1253	continue
1254	phaseDict[pfx+Stype[1].upper()+'waveType:'+str(i)] = waveType
1255	nx = 0
1256	for iw,wave in enumerate(Waves[1:]):
1257	if not iw:
1258	if waveType in ['ZigZag','Block']:
1259	nx = 1
1260	CSI = G2spc.GetSSfxuinel(waveType,Stype,1,at[cx:cx+3],SGData,SSGData)
1261	else:
1262	CSI = G2spc.GetSSfxuinel('Fourier',Stype,iw+1-nx,at[cx:cx+3],SGData,SSGData)
1263	uId,uCoef = CSI
1264	stiw = str(i)+':'+str(iw)
1265	if Stype == 'Spos':
1266	if waveType in ['ZigZag','Block',] and not iw:
1267	names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1268	equivs = [[],[], [],[],[]]
1269	else:
1270	names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1271	pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1272	equivs = [[],[],[], [],[],[]]
1273	elif Stype == 'Sadp':
1274	names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1275	pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1276	pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1277	pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1278	equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1279	elif Stype == 'Sfrac':
1280	equivs = [[],[]]
1281	if 'Crenel' in waveType and not iw:
1282	names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1283	else:
1284	names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1285	elif Stype == 'Smag':
1286	equivs = [[],[],[], [],[],[]]
1287	names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1288	pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1289	phaseDict.update(dict(zip(names,wave[0])))
1290	if wave[1]: #what do we do here for multiple terms in modulation constraints?
1291	for j in range(len(equivs)):
1292	if uId[j][0] > 0:
1293	phaseVary.append(names[j])
1294	equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1295	for equiv in equivs:
1296	if len(equiv) > 1:
1297	name = equiv[0][0]
1298	coef = equiv[0][1]
1299	for eqv in equiv[1:]:
1300	eqv[1] /= coef
1301	G2mv.StoreEquivalence(name,(eqv,))
1302	maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1303	textureData = General['SH Texture']
1304	if textureData['Order'] and not seqRef:
1305	phaseDict[pfx+'SHorder'] = textureData['Order']
1306	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1307	for item in ['omega','chi','phi']:
1308	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1309	if textureData['Sample '+item][0]:
1310	phaseVary.append(pfx+'SH '+item)
1311	for item in textureData['SH Coeff'][1]:
1312	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1313	if textureData['SH Coeff'][0]:
1314	phaseVary.append(pfx+item)
1315
1316	if Print:
1317	pFile.write('\n Phase name: %s\n'%General['Name'])
1318	pFile.write(135*'='+'\n')
1319	PrintFFtable(FFtable)
1320	PrintBLtable(BLtable)
1321	if General['Type'] == 'magnetic':
1322	PrintMFtable(MFtable)
1323	pFile.write('\n')
1324	#how do we print magnetic symmetry table? TBD
1325	if len(SSGtext): #if superstructure
1326	for line in SSGtext: pFile.write(line+'\n')
1327	if len(SSGtable):
1328	for item in SSGtable:
1329	line = ' %s '%(item)
1330	pFile.write(line+'\n')
1331	else:
1332	pFile.write(' ( 1) %s\n'%(SSGtable[0]))
1333	else:
1334	for line in SGtext: pFile.write(line+'\n')
1335	if len(SGtable):
1336	for item in SGtable:
1337	line = ' %s '%(item)
1338	pFile.write(line+'\n')
1339	else:
1340	pFile.write(' ( 1) %s\n'%(SGtable[0]))
1341	PrintRBObjects(resRBData,vecRBData)
1342	PrintAtoms(General,Atoms)
1343	if General['Type'] == 'magnetic':
1344	PrintMoments(General,Atoms)
1345	if General.get('Modulated',False):
1346	PrintWaves(General,Atoms)
1347	pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1348	(cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1349	if len(SSGtext): #if superstructure
1350	pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1351	(Vec[0],Vec[1],Vec[2],maxH,vRef))
1352	if not seqRef:
1353	PrintTexture(textureData)
1354	if name in RestraintDict:
1355	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1356	textureData,RestraintDict[name],pFile)
1357
1358	elif PawleyRef:
1359	if Print:
1360	pFile.write('\n Phase name: %s\n'%General['Name'])
1361	pFile.write(135*'='+'\n')
1362	pFile.write('\n')
1363	if len(SSGtext): #if superstructure
1364	for line in SSGtext: pFile.write(line+'\n')
1365	if len(SSGtable):
1366	for item in SSGtable:
1367	line = ' %s '%(item)
1368	pFile.write(line+'\n')
1369	else:
1370	pFile.write(' ( 1) %s\n'%SSGtable[0])
1371	else:
1372	for line in SGtext: pFile.write(line+'\n')
1373	if len(SGtable):
1374	for item in SGtable:
1375	line = ' %s '%(item)
1376	pFile.write(line+'\n')
1377	else:
1378	pFile.write(' ( 1) %s\n'%(SGtable[0]))
1379	pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1380	(cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1381	if len(SSGtext): #if superstructure
1382	pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1383	(Vec[0],Vec[1],Vec[2],maxH,vRef))
1384	pawleyVary = []
1385	for i,refl in enumerate(PawleyRef):
1386	phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1387	if im:
1388	pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1389	else:
1390	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1391	if refl[5+im]:
1392	pawleyVary.append(pfx+'PWLref:'+str(i))
1393	GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary) #does G2mv.StoreEquivalence
1394	phaseVary += pawleyVary
1395
1396	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
1397
1398	def cellFill(pfx,SGData,parmDict,sigDict):
1399	'''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1400	from the parameter and sig dictionaries.
1401
1402	:param str pfx: parameter prefix ("n::", where n is a phase number)
1403	:param dict SGdata: a symmetry object
1404	:param dict parmDict: a dictionary of parameters
1405	:param dict sigDict: a dictionary of uncertainties on parameters
1406
1407	:returns: A,sigA where each is a list of six terms with the A terms
1408	'''
1409	if SGData['SGLaue'] in ['-1',]:
1410	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1411	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1412	elif SGData['SGLaue'] in ['2/m',]:
1413	if SGData['SGUniq'] == 'a':
1414	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1415	0,0,parmDict[pfx+'A5']]
1416	elif SGData['SGUniq'] == 'b':
1417	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1418	0,parmDict[pfx+'A4'],0]
1419	else:
1420	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1421	parmDict[pfx+'A3'],0,0]
1422	elif SGData['SGLaue'] in ['mmm',]:
1423	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1424	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1425	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1426	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1427	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1428	parmDict[pfx+'A0'],0,0]
1429	elif SGData['SGLaue'] in ['3R', '3mR']:
1430	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1431	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1432	elif SGData['SGLaue'] in ['m3m','m3']:
1433	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1434
1435	try:
1436	if SGData['SGLaue'] in ['-1',]:
1437	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1438	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1439	elif SGData['SGLaue'] in ['2/m',]:
1440	if SGData['SGUniq'] == 'a':
1441	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1442	0,0,sigDict[pfx+'A5']]
1443	elif SGData['SGUniq'] == 'b':
1444	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1445	0,sigDict[pfx+'A4'],0]
1446	else:
1447	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1448	sigDict[pfx+'A3'],0,0]
1449	elif SGData['SGLaue'] in ['mmm',]:
1450	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1451	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1452	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1453	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1454	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1455	elif SGData['SGLaue'] in ['3R', '3mR']:
1456	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1457	elif SGData['SGLaue'] in ['m3m','m3']:
1458	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1459	except KeyError:
1460	sigA = [0,0,0,0,0,0]
1461	return A,sigA
1462
1463	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1464	'needs a doc string'
1465	if phaseRest:
1466	Amat = G2lat.cell2AB(cell)[0]
1467	cx,ct,cs = AtPtrs[:3]
1468	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1469	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1470	['ChemComp','Sites'],['Texture','HKLs']]
1471	for name,rest in names:
1472	itemRest = phaseRest[name]
1473	if rest in itemRest and itemRest[rest] and itemRest['Use']:
1474	pFile.write('\n %s restraint weight factor %10.3f Use: %s\n'%(name,itemRest['wtFactor'],str(itemRest['Use'])))
1475	if name in ['Bond','Angle','Plane','Chiral']:
1476	pFile.write(' calc obs sig delt/sig atoms(symOp): \n')
1477	for indx,ops,obs,esd in itemRest[rest]:
1478	try:
1479	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1480	AtName = ''
1481	for i,Aname in enumerate(AtNames):
1482	AtName += Aname
1483	if ops[i] == '1':
1484	AtName += '-'
1485	else:
1486	AtName += '+('+ops[i]+')-'
1487	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1488	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1489	if name == 'Bond':
1490	calc = G2mth.getRestDist(XYZ,Amat)
1491	elif name == 'Angle':
1492	calc = G2mth.getRestAngle(XYZ,Amat)
1493	elif name == 'Plane':
1494	calc = G2mth.getRestPlane(XYZ,Amat)
1495	elif name == 'Chiral':
1496	calc = G2mth.getRestChiral(XYZ,Amat)
1497	pFile.write(' %9.3f %9.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1]))
1498	except KeyError:
1499	del itemRest[rest]
1500	elif name in ['Torsion','Rama']:
1501	pFile.write(' atoms(symOp) calc obs sig delt/sig torsions: \n')
1502	coeffDict = itemRest['Coeff']
1503	for indx,ops,cofName,esd in itemRest[rest]:
1504	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1505	AtName = ''
1506	for i,Aname in enumerate(AtNames):
1507	AtName += Aname+'+('+ops[i]+')-'
1508	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1509	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1510	if name == 'Torsion':
1511	tor = G2mth.getRestTorsion(XYZ,Amat)
1512	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1513	pFile.write(' %8.3f %8.3f %.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1]))
1514	else:
1515	phi,psi = G2mth.getRestRama(XYZ,Amat)
1516	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1517	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]))
1518	elif name == 'ChemComp':
1519	pFile.write(' atoms mul*frac factor prod\n')
1520	for indx,factors,obs,esd in itemRest[rest]:
1521	try:
1522	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1523	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1524	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1525	mulfrac = mul*frac
1526	calcs = mulfracfactors
1527	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1528	pFile.write(' %10s %8.3f %8.3f %8.3f\n'%(atom,mf,fr,clc))
1529	pFile.write(' Sum: calc: %8.3f obs: %8.3f esd: %8.3f\n'%(np.sum(calcs),obs,esd))
1530	except KeyError:
1531	del itemRest[rest]
1532	elif name == 'Texture' and textureData['Order']:
1533	Start = False
1534	SHCoef = textureData['SH Coeff'][1]
1535	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1536	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1537	pFile.write (' HKL grid neg esd sum neg num neg use unit? unit esd \n')
1538	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1539	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1540	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1541	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1542	Z = ma.masked_greater(Z,0.0)
1543	num = ma.count(Z)
1544	sum = 0
1545	if num:
1546	sum = np.sum(Z)
1547	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))
1548
1549	def getCellEsd(pfx,SGData,A,covData):
1550	'needs a doc string'
1551	rVsq = G2lat.calc_rVsq(A)
1552	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1553	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1554	varyList = covData['varyList']
1555	covMatrix = covData['covMatrix']
1556	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1557	if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']:
1558	vcov[1,1] = vcov[3,3] = vcov[0,1] = vcov[1,0] = vcov[0,0]
1559	vcov[1,3] = vcov[3,1] = vcov[0,3] = vcov[3,0] = vcov[0,0]
1560	vcov[1,2] = vcov[2,1] = vcov[2,3] = vcov[3,2] = vcov[0,2]
1561	elif SGData['SGLaue'] in ['m3','m3m']:
1562	vcov[0:3,0:3] = vcov[0,0]
1563	elif SGData['SGLaue'] in ['4/m', '4/mmm']:
1564	vcov[0:2,0:2] = vcov[0,0]
1565	vcov[1,2] = vcov[2,1] = vcov[0,2]
1566	elif SGData['SGLaue'] in ['3R','3mR']:
1567	vcov[0:3,0:3] = vcov[0,0]
1568	# vcov[4,4] = vcov[5,5] = vcov[3,3]
1569	vcov[3:6,3:6] = vcov[3,3]
1570	vcov[0:3,3:6] = vcov[0,3]
1571	vcov[3:6,0:3] = vcov[3,0]
1572	delt = 1.e-9
1573	drVdA = np.zeros(6)
1574	for i in range(6):
1575	A[i] += delt
1576	drVdA[i] = G2lat.calc_rVsq(A)
1577	A[i] -= 2*delt
1578	drVdA[i] -= G2lat.calc_rVsq(A)
1579	A[i] += delt
1580	drVdA /= 2.*delt
1581	srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1582	Vol = 1/np.sqrt(rVsq)
1583	sigVol = Vol*3np.sqrt(srcvlsq)/2. #ok - checks with GSAS
1584
1585	dcdA = np.zeros((6,6))
1586	for i in range(6):
1587	pdcdA =np.zeros(6)
1588	A[i] += delt
1589	pdcdA += G2lat.A2cell(A)
1590	A[i] -= 2*delt
1591	pdcdA -= G2lat.A2cell(A)
1592	A[i] += delt
1593	dcdA[i] = pdcdA/(2.*delt)
1594
1595	sigMat = np.inner(dcdA,np.inner(dcdA,vcov))
1596	var = np.diag(sigMat)
1597	CS = np.where(var>0.,np.sqrt(var),0.)
1598	if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']:
1599	CS[3:6] = 0.0
1600	return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1601
1602	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1603	'''Called after a refinement to transfer parameters from the parameter dict to
1604	the phase(s) information read from a GPX file. Also prints values to the .lst file
1605	'''
1606
1607	def PrintAtomsAndSig(General,Atoms,atomsSig):
1608	pFile.write('\n Atoms:\n')
1609	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1610	if General['Type'] == 'macromolecular':
1611	line = ' res no residue chain '+line
1612	cx,ct,cs,cia = General['AtomPtrs']
1613	pFile.write(line+'\n')
1614	pFile.write(135*'-'+'\n')
1615	fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1616	cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1617	noFXsig = {cx:[10' ','%10s'],cx+1:[10' ','%10s'],cx+2:[10' ','%10s'],cx+3:[8' ','%8s']}
1618	for atyp in General['AtomTypes']: #zero composition
1619	General['NoAtoms'][atyp] = 0.
1620	for i,at in enumerate(Atoms):
1621	General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5]) #new composition
1622	if General['Type'] == 'macromolecular':
1623	name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1624	valstr = ' values: '
1625	sigstr = ' sig : '
1626	else:
1627	name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1628	valstr = ' values:'
1629	sigstr = ' sig :'
1630	for ind in range(cx,cx+4):
1631	sigind = str(i)+':'+str(ind)
1632	valstr += fmt[ind]%(at[ind])
1633	if sigind in atomsSig:
1634	sigstr += fmt[ind]%(atomsSig[sigind])
1635	else:
1636	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1637	if at[cia] == 'I':
1638	valstr += fmt[cia+1]%(at[cia+1])
1639	if '%d:%d'%(i,cia+1) in atomsSig:
1640	sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1641	else:
1642	sigstr += 8*' '
1643	else:
1644	valstr += 8*' '
1645	sigstr += 8*' '
1646	for ind in range(cia+2,cia+8):
1647	sigind = str(i)+':'+str(ind)
1648	valstr += fmt[ind]%(at[ind])
1649	if sigind in atomsSig:
1650	sigstr += fmt[ind]%(atomsSig[sigind])
1651	else:
1652	sigstr += 8*' '
1653	pFile.write(name+'\n')
1654	pFile.write(valstr+'\n')
1655	pFile.write(sigstr+'\n')
1656
1657	def PrintMomentsAndSig(General,Atoms,atomsSig):
1658	pFile.write('\n Magnetic Moments:\n') #add magnitude & angle, etc.? TBD
1659	line = ' name Mx My Mz'
1660	cx,ct,cs,cia = General['AtomPtrs']
1661	cmx = cx+4
1662	AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1663	pFile.write(line+'\n')
1664	pFile.write(135*'-'+'\n')
1665	fmt = {0:'%7s',ct:'%7s',cmx:'%10.5f',cmx+1:'%10.5f',cmx+2:'%10.5f'}
1666	noFXsig = {cmx:[10' ','%10s'],cmx+1:[10' ','%10s'],cmx+2:[10*' ','%10s']}
1667	for i,at in enumerate(Atoms):
1668	if AtInfo[at[ct]]:
1669	name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1670	valstr = ' values:'
1671	sigstr = ' sig :'
1672	for ind in range(cmx,cmx+3):
1673	sigind = str(i)+':'+str(ind)
1674	valstr += fmt[ind]%(at[ind])
1675	if sigind in atomsSig:
1676	sigstr += fmt[ind]%(atomsSig[sigind])
1677	else:
1678	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1679	pFile.write(name+'\n')
1680	pFile.write(valstr+'\n')
1681	pFile.write(sigstr+'\n')
1682
1683	def PrintWavesAndSig(General,Atoms,wavesSig):
1684	cx,ct,cs,cia = General['AtomPtrs']
1685	pFile.write('\n Modulation waves\n')
1686	names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1687	'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1688	'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1689	pFile.write(135*'-'+'\n')
1690	for i,at in enumerate(Atoms):
1691	AtomSS = at[-1]['SS1']
1692	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1693	Waves = AtomSS[Stype]
1694	if len(Waves):
1695	waveType = Waves[0]
1696	else:
1697	continue
1698	if len(Waves):
1699	pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
1700	(at[ct-1],at[cs],Stype,waveType))
1701	for iw,wave in enumerate(Waves[1:]):
1702	stiw = ':'+str(i)+':'+str(iw)
1703	namstr = ' names :'
1704	valstr = ' values:'
1705	sigstr = ' esds :'
1706	if Stype == 'Spos':
1707	nt = 6
1708	ot = 0
1709	if waveType in ['ZigZag','Block',] and not iw:
1710	nt = 5
1711	ot = 6
1712	for j in range(nt):
1713	name = names['Spos'][j+ot]
1714	namstr += '%12s'%(name)
1715	valstr += '%12.4f'%(wave[0][j])
1716	if name+stiw in wavesSig:
1717	sigstr += '%12.4f'%(wavesSig[name+stiw])
1718	else:
1719	sigstr += 12*' '
1720	elif Stype == 'Sfrac':
1721	ot = 0
1722	if 'Crenel' in waveType and not iw:
1723	ot = 2
1724	for j in range(2):
1725	name = names['Sfrac'][j+ot]
1726	namstr += '%12s'%(names['Sfrac'][j+ot])
1727	valstr += '%12.4f'%(wave[0][j])
1728	if name+stiw in wavesSig:
1729	sigstr += '%12.4f'%(wavesSig[name+stiw])
1730	else:
1731	sigstr += 12*' '
1732	elif Stype == 'Sadp':
1733	for j in range(12):
1734	name = names['Sadp'][j]
1735	namstr += '%10s'%(names['Sadp'][j])
1736	valstr += '%10.6f'%(wave[0][j])
1737	if name+stiw in wavesSig:
1738	sigstr += '%10.6f'%(wavesSig[name+stiw])
1739	else:
1740	sigstr += 10*' '
1741	elif Stype == 'Smag':
1742	for j in range(6):
1743	name = names['Smag'][j]
1744	namstr += '%12s'%(names['Smag'][j])
1745	valstr += '%12.4f'%(wave[0][j])
1746	if name+stiw in wavesSig:
1747	sigstr += '%12.4f'%(wavesSig[name+stiw])
1748	else:
1749	sigstr += 12*' '
1750
1751	pFile.write(namstr+'\n')
1752	pFile.write(valstr+'\n')
1753	pFile.write(sigstr+'\n')
1754
1755
1756	def PrintRBObjPOAndSig(rbfx,rbsx):
1757	namstr = ' names :'
1758	valstr = ' values:'
1759	sigstr = ' esds :'
1760	for i,px in enumerate(['Px:','Py:','Pz:']):
1761	name = pfx+rbfx+px+rbsx
1762	namstr += '%12s'%('Pos '+px[1])
1763	valstr += '%12.5f'%(parmDict[name])
1764	if name in sigDict:
1765	sigstr += '%12.5f'%(sigDict[name])
1766	else:
1767	sigstr += 12*' '
1768	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1769	name = pfx+rbfx+po+rbsx
1770	namstr += '%12s'%('Ori '+po[1])
1771	valstr += '%12.5f'%(parmDict[name])
1772	if name in sigDict:
1773	sigstr += '%12.5f'%(sigDict[name])
1774	else:
1775	sigstr += 12*' '
1776	pFile.write(namstr+'\n')
1777	pFile.write(valstr+'\n')
1778	pFile.write(sigstr+'\n')
1779
1780	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1781	namstr = ' names :'
1782	valstr = ' values:'
1783	sigstr = ' esds :'
1784	if 'N' not in TLS:
1785	pFile.write(' Thermal motion:\n')
1786	if 'T' in TLS:
1787	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1788	name = pfx+rbfx+pt+rbsx
1789	namstr += '%12s'%(pt[:3])
1790	valstr += '%12.5f'%(parmDict[name])
1791	if name in sigDict:
1792	sigstr += '%12.5f'%(sigDict[name])
1793	else:
1794	sigstr += 12*' '
1795	pFile.write(namstr+'\n')
1796	pFile.write(valstr+'\n')
1797	pFile.write(sigstr+'\n')
1798	if 'L' in TLS:
1799	namstr = ' names :'
1800	valstr = ' values:'
1801	sigstr = ' esds :'
1802	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1803	name = pfx+rbfx+pt+rbsx
1804	namstr += '%12s'%(pt[:3])
1805	valstr += '%12.3f'%(parmDict[name])
1806	if name in sigDict:
1807	sigstr += '%12.3f'%(sigDict[name])
1808	else:
1809	sigstr += 12*' '
1810	pFile.write(namstr+'\n')
1811	pFile.write(valstr+'\n')
1812	pFile.write(sigstr+'\n')
1813	if 'S' in TLS:
1814	namstr = ' names :'
1815	valstr = ' values:'
1816	sigstr = ' esds :'
1817	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1818	name = pfx+rbfx+pt+rbsx
1819	namstr += '%12s'%(pt[:3])
1820	valstr += '%12.4f'%(parmDict[name])
1821	if name in sigDict:
1822	sigstr += '%12.4f'%(sigDict[name])
1823	else:
1824	sigstr += 12*' '
1825	pFile.write(namstr+'\n')
1826	pFile.write(valstr+'\n')
1827	pFile.write(sigstr+'\n')
1828	if 'U' in TLS:
1829	name = pfx+rbfx+'U:'+rbsx
1830	namstr = ' names :'+'%12s'%('Uiso')
1831	valstr = ' values:'+'%12.5f'%(parmDict[name])
1832	if name in sigDict:
1833	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1834	else:
1835	sigstr = ' esds :'+12*' '
1836	pFile.write(namstr+'\n')
1837	pFile.write(valstr+'\n')
1838	pFile.write(sigstr+'\n')
1839
1840	def PrintRBObjTorAndSig(rbsx):
1841	namstr = ' names :'
1842	valstr = ' values:'
1843	sigstr = ' esds :'
1844	nTors = len(RBObj['Torsions'])
1845	if nTors:
1846	pFile.write(' Torsions:\n')
1847	for it in range(nTors):
1848	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1849	namstr += '%12s'%('Tor'+str(it))
1850	valstr += '%12.4f'%(parmDict[name])
1851	if name in sigDict:
1852	sigstr += '%12.4f'%(sigDict[name])
1853	pFile.write(namstr+'\n')
1854	pFile.write(valstr+'\n')
1855	pFile.write(sigstr+'\n')
1856
1857	def PrintSHtextureAndSig(textureData,SHtextureSig):
1858	pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
1859	names = ['omega','chi','phi']
1860	namstr = ' names :'
1861	ptstr = ' values:'
1862	sigstr = ' esds :'
1863	for name in names:
1864	namstr += '%12s'%(name)
1865	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1866	if 'Sample '+name in SHtextureSig:
1867	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1868	else:
1869	sigstr += 12*' '
1870	pFile.write(namstr+'\n')
1871	pFile.write(ptstr+'\n')
1872	pFile.write(sigstr+'\n')
1873	pFile.write('\n Texture coefficients:\n')
1874	SHcoeff = textureData['SH Coeff'][1]
1875	SHkeys = list(SHcoeff.keys())
1876	nCoeff = len(SHcoeff)
1877	nBlock = nCoeff//10+1
1878	iBeg = 0
1879	iFin = min(iBeg+10,nCoeff)
1880	for block in range(nBlock):
1881	namstr = ' names :'
1882	ptstr = ' values:'
1883	sigstr = ' esds :'
1884	for name in SHkeys[iBeg:iFin]:
1885	namstr += '%12s'%(name)
1886	ptstr += '%12.3f'%(SHcoeff[name])
1887	if name in SHtextureSig:
1888	sigstr += '%12.3f'%(SHtextureSig[name])
1889	else:
1890	sigstr += 12*' '
1891	pFile.write(namstr+'\n')
1892	pFile.write(ptstr+'\n')
1893	pFile.write(sigstr+'\n')
1894	iBeg += 10
1895	iFin = min(iBeg+10,nCoeff)
1896
1897	##########################################################################
1898	# SetPhaseData starts here
1899	pFile.write('\n Phases:\n')
1900	for phase in Phases:
1901	pFile.write(' Result for phase: %s\n'%phase)
1902	pFile.write(135*'='+'\n')
1903	Phase = Phases[phase]
1904	General = Phase['General']
1905	SGData = General['SGData']
1906	Atoms = Phase['Atoms']
1907	AtLookup = []
1908	if Atoms and not General.get('doPawley'):
1909	cx,ct,cs,cia = General['AtomPtrs']
1910	AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1911	cell = General['Cell']
1912	pId = Phase['pId']
1913	pfx = str(pId)+'::'
1914	if cell[0]:
1915	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1916	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1917	pFile.write(' Reciprocal metric tensor: \n')
1918	ptfmt = "%15.9f"
1919	names = ['A11','A22','A33','A12','A13','A23']
1920	namstr = ' names :'
1921	ptstr = ' values:'
1922	sigstr = ' esds :'
1923	for name,a,siga in zip(names,A,sigA):
1924	namstr += '%15s'%(name)
1925	ptstr += ptfmt%(a)
1926	if siga:
1927	sigstr += ptfmt%(siga)
1928	else:
1929	sigstr += 15*' '
1930	pFile.write(namstr+'\n')
1931	pFile.write(ptstr+'\n')
1932	pFile.write(sigstr+'\n')
1933	cell[1:7] = G2lat.A2cell(A)
1934	cell[7] = G2lat.calc_V(A)
1935	pFile.write(' New unit cell:\n')
1936	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1937	names = ['a','b','c','alpha','beta','gamma','Volume']
1938	namstr = ' names :'
1939	ptstr = ' values:'
1940	sigstr = ' esds :'
1941	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1942	namstr += '%12s'%(name)
1943	ptstr += fmt%(a)
1944	if siga:
1945	sigstr += fmt%(siga)
1946	else:
1947	sigstr += 12*' '
1948	pFile.write(namstr+'\n')
1949	pFile.write(ptstr+'\n')
1950	pFile.write(sigstr+'\n')
1951	ik = 6 #for Pawley stuff below
1952	if General.get('Modulated',False):
1953	ik = 7
1954	Vec,vRef,maxH = General['SuperVec']
1955	if vRef:
1956	pFile.write(' New modulation vector:\n')
1957	namstr = ' names :'
1958	ptstr = ' values:'
1959	sigstr = ' esds :'
1960	for var in ['mV0','mV1','mV2']:
1961	namstr += '%12s'%(pfx+var)
1962	ptstr += '%12.6f'%(parmDict[pfx+var])
1963	if pfx+var in sigDict:
1964	sigstr += '%12.6f'%(sigDict[pfx+var])
1965	else:
1966	sigstr += 12*' '
1967	pFile.write(namstr+'\n')
1968	pFile.write(ptstr+'\n')
1969	pFile.write(sigstr+'\n')
1970
1971	General['Mass'] = 0.
1972	if Phase['General'].get('doPawley'):
1973	pawleyRef = Phase['Pawley ref']
1974	for i,refl in enumerate(pawleyRef):
1975	key = pfx+'PWLref:'+str(i)
1976	refl[ik] = parmDict[key]
1977	if key in sigDict:
1978	refl[ik+1] = sigDict[key]
1979	else:
1980	refl[ik+1] = 0
1981	else:
1982	VRBIds = RBIds['Vector']
1983	RRBIds = RBIds['Residue']
1984	RBModels = Phase['RBModels']
1985	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1986	jrb = VRBIds.index(RBObj['RBId'])
1987	rbsx = str(irb)+':'+str(jrb)
1988	pFile.write(' Vector rigid body parameters:\n')
1989	PrintRBObjPOAndSig('RBV',rbsx)
1990	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1991	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1992	jrb = RRBIds.index(RBObj['RBId'])
1993	rbsx = str(irb)+':'+str(jrb)
1994	pFile.write(' Residue rigid body parameters:\n')
1995	PrintRBObjPOAndSig('RBR',rbsx)
1996	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1997	PrintRBObjTorAndSig(rbsx)
1998	atomsSig = {}
1999	wavesSig = {}
2000	cx,ct,cs,cia = General['AtomPtrs']
2001	for i,at in enumerate(Atoms):
2002	names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2003	cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2004	cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2005	cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2006	for ind in range(cx,cx+4):
2007	at[ind] = parmDict[names[ind]]
2008	if ind in range(cx,cx+3):
2009	name = names[ind].replace('A','dA')
2010	else:
2011	name = names[ind]
2012	if name in sigDict:
2013	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2014	if at[cia] == 'I':
2015	at[cia+1] = parmDict[names[cia+1]]
2016	if names[cia+1] in sigDict:
2017	atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2018	else:
2019	for ind in range(cia+2,cia+8):
2020	at[ind] = parmDict[names[ind]]
2021	if names[ind] in sigDict:
2022	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2023	if General['Type'] == 'magnetic':
2024	for ind in range(cx+4,cx+7):
2025	at[ind] = parmDict[names[ind]]
2026	if names[ind] in sigDict:
2027	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2028	ind = General['AtomTypes'].index(at[ct])
2029	General['Mass'] += General['AtomMass'][ind]at[cx+3]at[cx+5]
2030	if General.get('Modulated',False):
2031	AtomSS = at[-1]['SS1']
2032	for Stype in ['Sfrac','Spos','Sadp','Smag']:
2033	Waves = AtomSS[Stype]
2034	if len(Waves):
2035	waveType = Waves[0]
2036	else:
2037	continue
2038	for iw,wave in enumerate(Waves[1:]):
2039	stiw = str(i)+':'+str(iw)
2040	if Stype == 'Spos':
2041	if waveType in ['ZigZag','Block',] and not iw:
2042	names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2043	else:
2044	names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2045	'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2046	elif Stype == 'Sadp':
2047	names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2048	'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2049	'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2050	'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2051	elif Stype == 'Sfrac':
2052	if 'Crenel' in waveType and not iw:
2053	names = ['Fzero:'+stiw,'Fwid:'+stiw]
2054	else:
2055	names = ['Fsin:'+stiw,'Fcos:'+stiw]
2056	elif Stype == 'Smag':
2057	names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2058	'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2059	for iname,name in enumerate(names):
2060	AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2061	if pfx+name in sigDict:
2062	wavesSig[name] = sigDict[pfx+name]
2063
2064	PrintAtomsAndSig(General,Atoms,atomsSig)
2065	if General['Type'] == 'magnetic':
2066	PrintMomentsAndSig(General,Atoms,atomsSig)
2067	if General.get('Modulated',False):
2068	PrintWavesAndSig(General,Atoms,wavesSig)
2069
2070
2071	textureData = General['SH Texture']
2072	if textureData['Order']:
2073	SHtextureSig = {}
2074	for name in ['omega','chi','phi']:
2075	aname = pfx+'SH '+name
2076	textureData['Sample '+name][1] = parmDict[aname]
2077	if aname in sigDict:
2078	SHtextureSig['Sample '+name] = sigDict[aname]
2079	for name in textureData['SH Coeff'][1]:
2080	aname = pfx+name
2081	textureData['SH Coeff'][1][name] = parmDict[aname]
2082	if aname in sigDict:
2083	SHtextureSig[name] = sigDict[aname]
2084	PrintSHtextureAndSig(textureData,SHtextureSig)
2085	if phase in RestraintDict and not Phase['General'].get('doPawley'):
2086	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2087	textureData,RestraintDict[phase],pFile)
2088
2089	################################################################################
2090	##### Histogram & Phase data
2091	################################################################################
2092
2093	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2094	'''Loads the HAP histogram/phase information into dicts
2095
2096	:param dict Phases: phase information
2097	:param dict Histograms: Histogram information
2098	:param bool Print: prints information as it is read
2099	:param file pFile: file object to print to (the default, None causes printing to the console)
2100	:param bool resetRefList: Should the contents of the Reflection List be initialized
2101	on loading. The default, True, initializes the Reflection List as it is loaded.
2102
2103	:returns: (hapVary,hapDict,controlDict)
2104	* hapVary: list of refined variables
2105	* hapDict: dict with refined variables and their values
2106	* controlDict: dict with fixed parameters
2107	'''
2108
2109	def PrintSize(hapData):
2110	if hapData[0] in ['isotropic','uniaxial']:
2111	line = '\n Size model : %9s'%(hapData[0])
2112	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2113	if hapData[0] == 'uniaxial':
2114	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2115	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2116	pFile.write(line+'\n')
2117	else:
2118	pFile.write('\n Size model : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2119	(hapData[0],hapData[1][2],hapData[2][2]))
2120	Snames = ['S11','S22','S33','S12','S13','S23']
2121	ptlbls = ' names :'
2122	ptstr = ' values:'
2123	varstr = ' refine:'
2124	for i,name in enumerate(Snames):
2125	ptlbls += '%12s' % (name)
2126	ptstr += '%12.3f' % (hapData[4][i])
2127	varstr += '%12s' % (str(hapData[5][i]))
2128	pFile.write(ptlbls+'\n')
2129	pFile.write(ptstr+'\n')
2130	pFile.write(varstr+'\n')
2131
2132	def PrintMuStrain(hapData,SGData):
2133	if hapData[0] in ['isotropic','uniaxial']:
2134	line = '\n Mustrain model: %9s'%(hapData[0])
2135	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2136	if hapData[0] == 'uniaxial':
2137	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2138	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2139	pFile.write(line+'\n')
2140	else:
2141	pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2142	(hapData[0],hapData[1][2],hapData[2][2]))
2143	Snames = G2spc.MustrainNames(SGData)
2144	ptlbls = ' names :'
2145	ptstr = ' values:'
2146	varstr = ' refine:'
2147	for i,name in enumerate(Snames):
2148	ptlbls += '%12s' % (name)
2149	ptstr += '%12.1f' % (hapData[4][i])
2150	varstr += '%12s' % (str(hapData[5][i]))
2151	pFile.write(ptlbls+'\n')
2152	pFile.write(ptstr+'\n')
2153	pFile.write(varstr+'\n')
2154
2155	def PrintHStrain(hapData,SGData):
2156	pFile.write('\n Hydrostatic/elastic strain:\n')
2157	Hsnames = G2spc.HStrainNames(SGData)
2158	ptlbls = ' names :'
2159	ptstr = ' values:'
2160	varstr = ' refine:'
2161	for i,name in enumerate(Hsnames):
2162	ptlbls += '%12s' % (name)
2163	ptstr += '%12.4g' % (hapData[0][i])
2164	varstr += '%12s' % (str(hapData[1][i]))
2165	pFile.write(ptlbls+'\n')
2166	pFile.write(ptstr+'\n')
2167	pFile.write(varstr+'\n')
2168
2169	def PrintSHPO(hapData):
2170	pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2171	ptlbls = ' names :'
2172	ptstr = ' values:'
2173	for item in hapData[5]:
2174	ptlbls += '%12s'%(item)
2175	ptstr += '%12.3f'%(hapData[5][item])
2176	pFile.write(ptlbls+'\n')
2177	pFile.write(ptstr+'\n')
2178
2179	def PrintBabinet(hapData):
2180	pFile.write('\n Babinet form factor modification:\n')
2181	ptlbls = ' names :'
2182	ptstr = ' values:'
2183	varstr = ' refine:'
2184	for name in ['BabA','BabU']:
2185	ptlbls += '%12s' % (name)
2186	ptstr += '%12.6f' % (hapData[name][0])
2187	varstr += '%12s' % (str(hapData[name][1]))
2188	pFile.write(ptlbls+'\n')
2189	pFile.write(ptstr+'\n')
2190	pFile.write(varstr+'\n')
2191
2192	hapDict = {}
2193	hapVary = []
2194	controlDict = {}
2195
2196	for phase in Phases:
2197	HistoPhase = Phases[phase]['Histograms']
2198	SGData = Phases[phase]['General']['SGData']
2199	cell = Phases[phase]['General']['Cell'][1:7]
2200	A = G2lat.cell2A(cell)
2201	if Phases[phase]['General'].get('Modulated',False):
2202	SSGData = Phases[phase]['General']['SSGData']
2203	Vec,x,maxH = Phases[phase]['General']['SuperVec']
2204	pId = Phases[phase]['pId']
2205	histoList = list(HistoPhase.keys())
2206	histoList.sort()
2207	for histogram in histoList:
2208	try:
2209	Histogram = Histograms[histogram]
2210	except KeyError:
2211	#skip if histogram not included e.g. in a sequential refinement
2212	continue
2213	if not HistoPhase[histogram]['Use']: #remove previously created & now unused phase reflection list
2214	if phase in Histograms[histogram]['Reflection Lists']:
2215	del Histograms[histogram]['Reflection Lists'][phase]
2216	continue
2217	hapData = HistoPhase[histogram]
2218	hId = Histogram['hId']
2219	if 'PWDR' in histogram:
2220	limits = Histogram['Limits'][1]
2221	inst = Histogram['Instrument Parameters'][0] #TODO - grab table here if present
2222	if 'C' in inst['Type'][1]:
2223	try:
2224	wave = inst['Lam'][1]
2225	except KeyError:
2226	wave = inst['Lam1'][1]
2227	dmin = wave/(2.0*sind(limits[1]/2.0))
2228	elif 'T' in inst['Type'][0]:
2229	dmin = limits[0]/inst['difC'][1]
2230	pfx = str(pId)+':'+str(hId)+':'
2231	if Phases[phase]['General']['doPawley']:
2232	hapDict[pfx+'LeBail'] = False #Pawley supercedes LeBail
2233	hapDict[pfx+'newLeBail'] = True
2234	Tmin = G2lat.Dsp2pos(inst,dmin)
2235	if 'C' in inst['Type'][1]:
2236	limits[1] = min(limits[1],Tmin)
2237	else:
2238	limits[0] = max(limits[0],Tmin)
2239	else:
2240	hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2241	hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2242	if Phases[phase]['General']['Type'] == 'magnetic':
2243	dmin = max(dmin,Phases[phase]['General']['MagDmin'])
2244	for item in ['Scale','Extinction']:
2245	hapDict[pfx+item] = hapData[item][0]
2246	if hapData[item][1] and not hapDict[pfx+'LeBail']:
2247	hapVary.append(pfx+item)
2248	names = G2spc.HStrainNames(SGData)
2249	HSvals = []
2250	for i,name in enumerate(names):
2251	hapDict[pfx+name] = hapData['HStrain'][0][i]
2252	HSvals.append(hapDict[pfx+name])
2253	if hapData['HStrain'][1][i]:
2254	# if hapData['HStrain'][1][i] and not hapDict[pfx+'LeBail']:
2255	hapVary.append(pfx+name)
2256	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2257	if hapData['Pref.Ori.'][0] == 'MD':
2258	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2259	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2260	if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2261	hapVary.append(pfx+'MD')
2262	else: #'SH' spherical harmonics
2263	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2264	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2265	controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2266	controlDict[pfx+'SHhkl'] = []
2267	try: #patch for old Pref.Ori. items
2268	controlDict[pfx+'SHtoler'] = 0.1
2269	if hapData['Pref.Ori.'][6][0] != '':
2270	controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2271	controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2272	except IndexError:
2273	pass
2274	for item in hapData['Pref.Ori.'][5]:
2275	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2276	if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2277	hapVary.append(pfx+item)
2278	for item in ['Mustrain','Size']:
2279	controlDict[pfx+item+'Type'] = hapData[item][0]
2280	hapDict[pfx+item+';mx'] = hapData[item][1][2]
2281	if hapData[item][2][2]:
2282	hapVary.append(pfx+item+';mx')
2283	if hapData[item][0] in ['isotropic','uniaxial']:
2284	hapDict[pfx+item+';i'] = hapData[item][1][0]
2285	if hapData[item][2][0]:
2286	hapVary.append(pfx+item+';i')
2287	if hapData[item][0] == 'uniaxial':
2288	controlDict[pfx+item+'Axis'] = hapData[item][3]
2289	hapDict[pfx+item+';a'] = hapData[item][1][1]
2290	if hapData[item][2][1]:
2291	hapVary.append(pfx+item+';a')
2292	else: #generalized for mustrain or ellipsoidal for size
2293	Nterms = len(hapData[item][4])
2294	if item == 'Mustrain':
2295	names = G2spc.MustrainNames(SGData)
2296	pwrs = []
2297	for name in names:
2298	h,k,l = name[1:]
2299	pwrs.append([int(h),int(k),int(l)])
2300	controlDict[pfx+'MuPwrs'] = pwrs
2301	for i in range(Nterms):
2302	sfx = ';'+str(i)
2303	hapDict[pfx+item+sfx] = hapData[item][4][i]
2304	if hapData[item][5][i]:
2305	hapVary.append(pfx+item+sfx)
2306	if Phases[phase]['General']['Type'] != 'magnetic':
2307	for bab in ['BabA','BabU']:
2308	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2309	if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2310	hapVary.append(pfx+bab)
2311
2312	if Print:
2313	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2314	pFile.write(135*'='+'\n')
2315	if hapDict[pfx+'LeBail']:
2316	pFile.write(' Perform LeBail extraction\n')
2317	else:
2318	pFile.write(' Phase fraction : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2319	pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2320	if hapData['Pref.Ori.'][0] == 'MD':
2321	Ax = hapData['Pref.Ori.'][3]
2322	pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2323	(hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2324	else: #'SH' for spherical harmonics
2325	PrintSHPO(hapData['Pref.Ori.'])
2326	pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2327	PrintSize(hapData['Size'])
2328	PrintMuStrain(hapData['Mustrain'],SGData)
2329	PrintHStrain(hapData['HStrain'],SGData)
2330	if Phases[phase]['General']['Type'] != 'magnetic':
2331	if hapData['Babinet']['BabA'][0]:
2332	PrintBabinet(hapData['Babinet'])
2333	if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData['LeBail'] and hapData['newLeBail'])):
2334	if hapData.get('LeBail',True): #stop regeneating reflections for LeBail
2335	hapData['newLeBail'] = False
2336	refList = []
2337	Uniq = []
2338	Phi = []
2339	useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2340	if Phases[phase]['General'].get('Modulated',False):
2341	ifSuper = True
2342	HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2343	HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2344	for h,k,l,m,d in HKLd:
2345	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2346	mul *= 2 # for powder overlap of Friedel pairs
2347	if m or not ext or useExt:
2348	if 'C' in inst['Type'][0]:
2349	pos = G2lat.Dsp2pos(inst,d)
2350	if limits[0] < pos < limits[1]:
2351	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])
2352	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2353	Uniq.append(uniq)
2354	Phi.append(phi)
2355	elif 'T' in inst['Type'][0]:
2356	pos = G2lat.Dsp2pos(inst,d)
2357	if limits[0] < pos < limits[1]:
2358	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2359	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])
2360	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2361	#TODO - if tabulated put alp & bet in here
2362	Uniq.append(uniq)
2363	Phi.append(phi)
2364	else:
2365	ifSuper = False
2366	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2367	HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2368	for h,k,l,d in HKLd:
2369	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2370	mul *= 2 # for powder overlap of Friedel pairs
2371	if ext and not useExt:
2372	continue
2373	if 'C' in inst['Type'][0]:
2374	pos = G2lat.Dsp2pos(inst,d)
2375	if limits[0] < pos < limits[1]:
2376	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])
2377	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2378	Uniq.append(uniq)
2379	Phi.append(phi)
2380	elif 'T' in inst['Type'][0]:
2381	pos = G2lat.Dsp2pos(inst,d)
2382	if limits[0] < pos < limits[1]:
2383	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2384	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])
2385	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2386	Uniq.append(uniq)
2387	Phi.append(phi)
2388	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2389	elif 'HKLF' in histogram:
2390	inst = Histogram['Instrument Parameters'][0]
2391	hId = Histogram['hId']
2392	hfx = ':%d:'%(hId)
2393	for item in inst:
2394	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2395	hapDict[hfx+item] = inst[item][1]
2396	pfx = str(pId)+':'+str(hId)+':'
2397	hapDict[pfx+'Scale'] = hapData['Scale'][0]
2398	if hapData['Scale'][1]:
2399	hapVary.append(pfx+'Scale')
2400
2401	extApprox,extType,extParms = hapData['Extinction']
2402	controlDict[pfx+'EType'] = extType
2403	controlDict[pfx+'EApprox'] = extApprox
2404	if 'C' in inst['Type'][0]:
2405	controlDict[pfx+'Tbar'] = extParms['Tbar']
2406	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2407	if 'Primary' in extType:
2408	Ekey = ['Ep',]
2409	elif 'I & II' in extType:
2410	Ekey = ['Eg','Es']
2411	elif 'Secondary Type II' == extType:
2412	Ekey = ['Es',]
2413	elif 'Secondary Type I' == extType:
2414	Ekey = ['Eg',]
2415	else: #'None'
2416	Ekey = []
2417	for eKey in Ekey:
2418	hapDict[pfx+eKey] = extParms[eKey][0]
2419	if extParms[eKey][1]:
2420	hapVary.append(pfx+eKey)
2421	for bab in ['BabA','BabU']:
2422	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2423	if hapData['Babinet'][bab][1]:
2424	hapVary.append(pfx+bab)
2425	Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2426	if len(Twins) == 1:
2427	hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2428	if hapData.get('Flack',[0,False])[1]:
2429	hapVary.append(pfx+'Flack')
2430	sumTwFr = 0.
2431	controlDict[pfx+'TwinLaw'] = []
2432	controlDict[pfx+'TwinInv'] = []
2433	NTL = 0
2434	for it,twin in enumerate(Twins):
2435	if 'bool' in str(type(twin[0])):
2436	controlDict[pfx+'TwinInv'].append(twin[0])
2437	controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2438	else:
2439	NTL += 1
2440	controlDict[pfx+'TwinInv'].append(False)
2441	controlDict[pfx+'TwinLaw'].append(twin[0])
2442	if it:
2443	hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2444	sumTwFr += twin[1]
2445	else:
2446	hapDict[pfx+'TwinFr:0'] = twin[1][0]
2447	controlDict[pfx+'TwinNMN'] = twin[1][1]
2448	if Twins[0][1][1]:
2449	hapVary.append(pfx+'TwinFr:'+str(it))
2450	controlDict[pfx+'NTL'] = NTL
2451	#need to make constraint on TwinFr
2452	controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2453	if len(Twins) > 1: #force sum to unity
2454	hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2455	if Print:
2456	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2457	pFile.write(135*'='+'\n')
2458	pFile.write(' Scale factor : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2459	if extType != 'None':
2460	pFile.write(' Extinction Type: %15s approx: %10s\n'%(extType,extApprox))
2461	text = ' Parameters :'
2462	for eKey in Ekey:
2463	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2464	pFile.write(text+'\n')
2465	if hapData['Babinet']['BabA'][0]:
2466	PrintBabinet(hapData['Babinet'])
2467	if not SGData['SGInv'] and len(Twins) == 1:
2468	pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2469	if len(Twins) > 1:
2470	for it,twin in enumerate(Twins):
2471	if 'bool' in str(type(twin[0])):
2472	pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine?\n'%
2473	(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1])
2474	else:
2475	pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2476	(str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],Twins[0][1][1]))
2477
2478	Histogram['Reflection Lists'] = phase
2479
2480	return hapVary,hapDict,controlDict
2481
2482	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2483	'needs a doc string'
2484
2485	def PrintSizeAndSig(hapData,sizeSig):
2486	line = '\n Size model: %9s'%(hapData[0])
2487	refine = False
2488	if hapData[0] in ['isotropic','uniaxial']:
2489	line += ' equatorial:%12.5f'%(hapData[1][0])
2490	if sizeSig[0][0]:
2491	line += ', sig:%8.4f'%(sizeSig[0][0])
2492	refine = True
2493	if hapData[0] == 'uniaxial':
2494	line += ' axial:%12.4f'%(hapData[1][1])
2495	if sizeSig[0][1]:
2496	refine = True
2497	line += ', sig:%8.4f'%(sizeSig[0][1])
2498	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2499	if sizeSig[0][2]:
2500	refine = True
2501	line += ', sig:%8.4f'%(sizeSig[0][2])
2502	if refine:
2503	pFile.write(line+'\n')
2504	else:
2505	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2506	if sizeSig[0][2]:
2507	refine = True
2508	line += ', sig:%8.4f'%(sizeSig[0][2])
2509	Snames = ['S11','S22','S33','S12','S13','S23']
2510	ptlbls = ' name :'
2511	ptstr = ' value :'
2512	sigstr = ' sig :'
2513	for i,name in enumerate(Snames):
2514	ptlbls += '%12s' % (name)
2515	ptstr += '%12.6f' % (hapData[4][i])
2516	if sizeSig[1][i]:
2517	refine = True
2518	sigstr += '%12.6f' % (sizeSig[1][i])
2519	else:
2520	sigstr += 12*' '
2521	if refine:
2522	pFile.write(line+'\n')
2523	pFile.write(ptlbls+'\n')
2524	pFile.write(ptstr+'\n')
2525	pFile.write(sigstr+'\n')
2526
2527	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2528	line = '\n Mustrain model: %9s\n'%(hapData[0])
2529	refine = False
2530	if hapData[0] in ['isotropic','uniaxial']:
2531	line += ' equatorial:%12.1f'%(hapData[1][0])
2532	if mustrainSig[0][0]:
2533	line += ', sig:%8.1f'%(mustrainSig[0][0])
2534	refine = True
2535	if hapData[0] == 'uniaxial':
2536	line += ' axial:%12.1f'%(hapData[1][1])
2537	if mustrainSig[0][1]:
2538	line += ', sig:%8.1f'%(mustrainSig[0][1])
2539	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2540	if mustrainSig[0][2]:
2541	refine = True
2542	line += ', sig:%8.3f'%(mustrainSig[0][2])
2543	if refine:
2544	pFile.write(line+'\n')
2545	else:
2546	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2547	if mustrainSig[0][2]:
2548	refine = True
2549	line += ', sig:%8.3f'%(mustrainSig[0][2])
2550	Snames = G2spc.MustrainNames(SGData)
2551	ptlbls = ' name :'
2552	ptstr = ' value :'
2553	sigstr = ' sig :'
2554	for i,name in enumerate(Snames):
2555	ptlbls += '%12s' % (name)
2556	ptstr += '%12.1f' % (hapData[4][i])
2557	if mustrainSig[1][i]:
2558	refine = True
2559	sigstr += '%12.1f' % (mustrainSig[1][i])
2560	else:
2561	sigstr += 12*' '
2562	if refine:
2563	pFile.write(line+'\n')
2564	pFile.write(ptlbls+'\n')
2565	pFile.write(ptstr+'\n')
2566	pFile.write(sigstr+'\n')
2567
2568	def PrintHStrainAndSig(hapData,strainSig,SGData):
2569	Hsnames = G2spc.HStrainNames(SGData)
2570	ptlbls = ' name :'
2571	ptstr = ' value :'
2572	sigstr = ' sig :'
2573	refine = False
2574	for i,name in enumerate(Hsnames):
2575	ptlbls += '%12s' % (name)
2576	ptstr += '%12.4g' % (hapData[0][i])
2577	if name in strainSig:
2578	refine = True
2579	sigstr += '%12.4g' % (strainSig[name])
2580	else:
2581	sigstr += 12*' '
2582	if refine:
2583	pFile.write('\n Hydrostatic/elastic strain:\n')
2584	pFile.write(ptlbls+'\n')
2585	pFile.write(ptstr+'\n')
2586	pFile.write(sigstr+'\n')
2587
2588	def PrintSHPOAndSig(pfx,hapData,POsig):
2589	pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2590	ptlbls = ' names :'
2591	ptstr = ' values:'
2592	sigstr = ' sig :'
2593	for item in hapData[5]:
2594	ptlbls += '%12s'%(item)
2595	ptstr += '%12.3f'%(hapData[5][item])
2596	if pfx+item in POsig:
2597	sigstr += '%12.3f'%(POsig[pfx+item])
2598	else:
2599	sigstr += 12*' '
2600	pFile.write(ptlbls+'\n')
2601	pFile.write(ptstr+'\n')
2602	pFile.write(sigstr+'\n')
2603
2604	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2605	pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2606	text = ''
2607	for item in hapData[2]:
2608	if pfx+item in ScalExtSig:
2609	text += ' %s: '%(item)
2610	text += '%12.2e'%(hapData[2][item][0])
2611	if pfx+item in ScalExtSig:
2612	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2613	pFile.write(text+'\n')
2614
2615	def PrintBabinetAndSig(pfx,hapData,BabSig):
2616	pFile.write('\n Babinet form factor modification:\n')
2617	ptlbls = ' names :'
2618	ptstr = ' values:'
2619	sigstr = ' sig :'
2620	for item in hapData:
2621	ptlbls += '%12s'%(item)
2622	ptstr += '%12.3f'%(hapData[item][0])
2623	if pfx+item in BabSig:
2624	sigstr += '%12.3f'%(BabSig[pfx+item])
2625	else:
2626	sigstr += 12*' '
2627	pFile.write(ptlbls+'\n')
2628	pFile.write(ptstr+'\n')
2629	pFile.write(sigstr+'\n')
2630
2631	def PrintTwinsAndSig(pfx,twinData,TwinSig):
2632	pFile.write('\n Twin Law fractions :\n')
2633	ptlbls = ' names :'
2634	ptstr = ' values:'
2635	sigstr = ' sig :'
2636	for it,item in enumerate(twinData):
2637	ptlbls += ' twin #%d'%(it)
2638	if it:
2639	ptstr += '%12.3f'%(item[1])
2640	else:
2641	ptstr += '%12.3f'%(item[1][0])
2642	if pfx+'TwinFr:'+str(it) in TwinSig:
2643	sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2644	else:
2645	sigstr += 12*' '
2646	pFile.write(ptlbls+'\n')
2647	pFile.write(ptstr+'\n')
2648	pFile.write(sigstr+'\n')
2649
2650
2651	PhFrExtPOSig = {}
2652	SizeMuStrSig = {}
2653	ScalExtSig = {}
2654	BabSig = {}
2655	TwinFrSig = {}
2656	wtFrSum = {}
2657	for phase in Phases:
2658	HistoPhase = Phases[phase]['Histograms']
2659	General = Phases[phase]['General']
2660	SGData = General['SGData']
2661	pId = Phases[phase]['pId']
2662	histoList = list(HistoPhase.keys())
2663	histoList.sort()
2664	for histogram in histoList:
2665	try:
2666	Histogram = Histograms[histogram]
2667	except KeyError:
2668	#skip if histogram not included e.g. in a sequential refinement
2669	continue
2670	if not Phases[phase]['Histograms'][histogram]['Use']:
2671	#skip if phase absent from this histogram
2672	continue
2673	hapData = HistoPhase[histogram]
2674	hId = Histogram['hId']
2675	pfx = str(pId)+':'+str(hId)+':'
2676	if hId not in wtFrSum:
2677	wtFrSum[hId] = 0.
2678	if 'PWDR' in histogram:
2679	for item in ['Scale','Extinction']:
2680	hapData[item][0] = parmDict[pfx+item]
2681	if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2682	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2683	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2684	if hapData['Pref.Ori.'][0] == 'MD':
2685	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2686	if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2687	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2688	else: #'SH' spherical harmonics
2689	for item in hapData['Pref.Ori.'][5]:
2690	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2691	if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2692	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2693	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2694	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2695	pfx+'HStrain':{}})
2696	for item in ['Mustrain','Size']:
2697	hapData[item][1][2] = parmDict[pfx+item+';mx']
2698	# hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2699	if pfx+item+';mx' in sigDict:
2700	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2701	if hapData[item][0] in ['isotropic','uniaxial']:
2702	hapData[item][1][0] = parmDict[pfx+item+';i']
2703	if item == 'Size':
2704	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2705	if pfx+item+';i' in sigDict:
2706	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2707	if hapData[item][0] == 'uniaxial':
2708	hapData[item][1][1] = parmDict[pfx+item+';a']
2709	if item == 'Size':
2710	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2711	if pfx+item+';a' in sigDict:
2712	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2713	else: #generalized for mustrain or ellipsoidal for size
2714	Nterms = len(hapData[item][4])
2715	for i in range(Nterms):
2716	sfx = ';'+str(i)
2717	hapData[item][4][i] = parmDict[pfx+item+sfx]
2718	if pfx+item+sfx in sigDict:
2719	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2720	names = G2spc.HStrainNames(SGData)
2721	for i,name in enumerate(names):
2722	hapData['HStrain'][0][i] = parmDict[pfx+name]
2723	if pfx+name in sigDict:
2724	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2725	if Phases[phase]['General']['Type'] != 'magnetic':
2726	for name in ['BabA','BabU']:
2727	hapData['Babinet'][name][0] = parmDict[pfx+name]
2728	if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
2729	BabSig[pfx+name] = sigDict[pfx+name]
2730
2731	elif 'HKLF' in histogram:
2732	for item in ['Scale','Flack']:
2733	if parmDict.get(pfx+item):
2734	hapData[item][0] = parmDict[pfx+item]
2735	if pfx+item in sigDict:
2736	ScalExtSig[pfx+item] = sigDict[pfx+item]
2737	for item in ['Ep','Eg','Es']:
2738	if parmDict.get(pfx+item):
2739	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2740	if pfx+item in sigDict:
2741	ScalExtSig[pfx+item] = sigDict[pfx+item]
2742	for item in ['BabA','BabU']:
2743	hapData['Babinet'][item][0] = parmDict[pfx+item]
2744	if pfx+item in sigDict:
2745	BabSig[pfx+item] = sigDict[pfx+item]
2746	if 'Twins' in hapData:
2747	it = 1
2748	sumTwFr = 0.
2749	while True:
2750	try:
2751	hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2752	if pfx+'TwinFr:'+str(it) in sigDict:
2753	TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2754	if it:
2755	sumTwFr += hapData['Twins'][it][1]
2756	it += 1
2757	except KeyError:
2758	break
2759	hapData['Twins'][0][1][0] = 1.-sumTwFr
2760
2761	if Print:
2762	for phase in Phases:
2763	HistoPhase = Phases[phase]['Histograms']
2764	General = Phases[phase]['General']
2765	SGData = General['SGData']
2766	pId = Phases[phase]['pId']
2767	histoList = list(HistoPhase.keys())
2768	histoList.sort()
2769	for histogram in histoList:
2770	try:
2771	Histogram = Histograms[histogram]
2772	except KeyError:
2773	#skip if histogram not included e.g. in a sequential refinement
2774	continue
2775	hapData = HistoPhase[histogram]
2776	hId = Histogram['hId']
2777	Histogram['Residuals'][str(pId)+'::Name'] = phase
2778	pfx = str(pId)+':'+str(hId)+':'
2779	hfx = ':%s:'%(hId)
2780	if pfx+'Nref' not in Histogram['Residuals']: #skip not used phase in histogram
2781	continue
2782	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2783	pFile.write(135*'='+'\n')
2784	if 'PWDR' in histogram:
2785	pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
2786	(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
2787	pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
2788	pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
2789
2790	if parmDict[pfx+'LeBail']:
2791	pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
2792	else:
2793	if pfx+'Scale' in PhFrExtPOSig:
2794	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2795	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2796	pFile.write(' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f\n'%
2797	(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
2798	if pfx+'Extinction' in PhFrExtPOSig:
2799	pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
2800	if hapData['Pref.Ori.'][0] == 'MD':
2801	if pfx+'MD' in PhFrExtPOSig:
2802	pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
2803	else:
2804	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2805	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2806	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2807	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2808	if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
2809	if len(BabSig):
2810	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2811
2812	elif 'HKLF' in histogram:
2813	Inst = Histogram['Instrument Parameters'][0]
2814	pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
2815	(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2816	Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
2817	if FFtables != None and 'N' not in Inst['Type'][0]:
2818	PrintFprime(FFtables,hfx,pFile)
2819	pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
2820	if pfx+'Scale' in ScalExtSig:
2821	pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
2822	if hapData['Extinction'][0] != 'None':
2823	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2824	if len(BabSig):
2825	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2826	if pfx+'Flack' in ScalExtSig:
2827	pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
2828	if len(TwinFrSig):
2829	PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2830
2831	################################################################################
2832	##### Histogram data
2833	################################################################################
2834
2835	def GetHistogramData(Histograms,Print=True,pFile=None):
2836	'needs a doc string'
2837
2838	def GetBackgroundParms(hId,Background):
2839	Back = Background[0]
2840	DebyePeaks = Background[1]
2841	bakType,bakFlag = Back[:2]
2842	backVals = Back[3:]
2843	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2844	backDict = dict(zip(backNames,backVals))
2845	backVary = []
2846	if bakFlag:
2847	backVary = backNames
2848	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2849	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2850	debyeDict = {}
2851	debyeList = []
2852	for i in range(DebyePeaks['nDebye']):
2853	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2854	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2855	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2856	debyeVary = []
2857	for item in debyeList:
2858	if item[1]:
2859	debyeVary.append(item[0])
2860	backDict.update(debyeDict)
2861	backVary += debyeVary
2862	peakDict = {}
2863	peakList = []
2864	for i in range(DebyePeaks['nPeaks']):
2865	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2866	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2867	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2868	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2869	peakVary = []
2870	for item in peakList:
2871	if item[1]:
2872	peakVary.append(item[0])
2873	backDict.update(peakDict)
2874	backVary += peakVary
2875	return bakType,backDict,backVary
2876
2877	def GetInstParms(hId,Inst):
2878	#patch
2879	if 'Z' not in Inst:
2880	Inst['Z'] = [0.0,0.0,False]
2881	dataType = Inst['Type'][0]
2882	instDict = {}
2883	insVary = []
2884	pfx = ':'+str(hId)+':'
2885	insKeys = list(Inst.keys())
2886	insKeys.sort()
2887	for item in insKeys:
2888	insName = pfx+item
2889	instDict[insName] = Inst[item][1]
2890	if len(Inst[item]) > 2 and Inst[item][2]:
2891	insVary.append(insName)
2892	if 'C' in dataType:
2893	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2894	elif 'T' in dataType: #trap zero alp, bet coeff.
2895	if not instDict[pfx+'alpha']:
2896	instDict[pfx+'alpha'] = 1.0
2897	if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
2898	instDict[pfx+'beta-1'] = 1.0
2899	return dataType,instDict,insVary
2900
2901	def GetSampleParms(hId,Sample):
2902	sampVary = []
2903	hfx = ':'+str(hId)+':'
2904	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2905	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2906	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2907	if key in Sample:
2908	sampDict[hfx+key] = Sample[key]
2909	Type = Sample['Type']
2910	if 'Bragg' in Type: #Bragg-Brentano
2911	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2912	sampDict[hfx+item] = Sample[item][0]
2913	if Sample[item][1]:
2914	sampVary.append(hfx+item)
2915	elif 'Debye' in Type: #Debye-Scherrer
2916	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2917	sampDict[hfx+item] = Sample[item][0]
2918	if Sample[item][1]:
2919	sampVary.append(hfx+item)
2920	return Type,sampDict,sampVary
2921
2922	def PrintBackground(Background):
2923	Back = Background[0]
2924	DebyePeaks = Background[1]
2925	pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
2926	line = ' Coefficients: '
2927	for i,back in enumerate(Back[3:]):
2928	line += '%10.3f'%(back)
2929	if i and not i%10:
2930	line += '\n'+15*' '
2931	pFile.write(line+'\n')
2932	if DebyePeaks['nDebye']:
2933	pFile.write('\n Debye diffuse scattering coefficients\n')
2934	parms = ['DebyeA','DebyeR','DebyeU']
2935	line = ' names : '
2936	for parm in parms:
2937	line += '%8s refine?'%(parm)
2938	pFile.write(line+'\n')
2939	for j,term in enumerate(DebyePeaks['debyeTerms']):
2940	line = ' term'+'%2d'%(j)+':'
2941	for i in range(3):
2942	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2943	pFile.write(line+'\n')
2944	if DebyePeaks['nPeaks']:
2945	pFile.write('\n Single peak coefficients\n')
2946	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2947	line = ' names : '
2948	for parm in parms:
2949	line += '%8s refine?'%(parm)
2950	pFile.write(line+'\n')
2951	for j,term in enumerate(DebyePeaks['peaksList']):
2952	line = ' peak'+'%2d'%(j)+':'
2953	for i in range(4):
2954	line += '%12.3f %5s'%(term[2i],bool(term[2i+1]))
2955	pFile.write(line+'\n')
2956
2957	def PrintInstParms(Inst):
2958	pFile.write('\n Instrument Parameters:\n')
2959	insKeys = list(Inst.keys())
2960	insKeys.sort()
2961	iBeg = 0
2962	Ok = True
2963	while Ok:
2964	ptlbls = ' name :'
2965	ptstr = ' value :'
2966	varstr = ' refine:'
2967	iFin = min(iBeg+9,len(insKeys))
2968	for item in insKeys[iBeg:iFin]:
2969	if item not in ['Type','Source']:
2970	ptlbls += '%12s' % (item)
2971	ptstr += '%12.6f' % (Inst[item][1])
2972	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2973	varstr += 12*' '
2974	else:
2975	varstr += '%12s' % (str(bool(Inst[item][2])))
2976	pFile.write(ptlbls+'\n')
2977	pFile.write(ptstr+'\n')
2978	pFile.write(varstr+'\n')
2979	iBeg = iFin
2980	if iBeg == len(insKeys):
2981	Ok = False
2982	else:
2983	pFile.write('\n')
2984
2985	def PrintSampleParms(Sample):
2986	pFile.write('\n Sample Parameters:\n')
2987	pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
2988	(Sample['Omega'],Sample['Chi'],Sample['Phi']))
2989	ptlbls = ' name :'
2990	ptstr = ' value :'
2991	varstr = ' refine:'
2992	if 'Bragg' in Sample['Type']:
2993	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2994	ptlbls += '%14s'%(item)
2995	ptstr += '%14.4f'%(Sample[item][0])
2996	varstr += '%14s'%(str(bool(Sample[item][1])))
2997
2998	elif 'Debye' in Type: #Debye-Scherrer
2999	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3000	ptlbls += '%14s'%(item)
3001	ptstr += '%14.4f'%(Sample[item][0])
3002	varstr += '%14s'%(str(bool(Sample[item][1])))
3003
3004	pFile.write(ptlbls+'\n')
3005	pFile.write(ptstr+'\n')
3006	pFile.write(varstr+'\n')
3007
3008	histDict = {}
3009	histVary = []
3010	controlDict = {}
3011	histoList = list(Histograms.keys())
3012	histoList.sort()
3013	for histogram in histoList:
3014	if 'PWDR' in histogram:
3015	Histogram = Histograms[histogram]
3016	hId = Histogram['hId']
3017	pfx = ':'+str(hId)+':'
3018	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3019	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3020	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3021	for excl in controlDict[pfx+'Exclude']:
3022	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3023	if controlDict[pfx+'Exclude']:
3024	ma.mask_rows(Histogram['Data'])
3025	Background = Histogram['Background']
3026	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3027	controlDict[pfx+'bakType'] = Type
3028	histDict.update(bakDict)
3029	histVary += bakVary
3030
3031	Inst = Histogram['Instrument Parameters'] #TODO ? ignores tabulated alp,bet & delt for TOF
3032	if 'T' in Type and len(Inst[1]): #patch - back-to-back exponential contribution to TOF line shape is removed
3033	print ('Warning: tabulated profile coefficients are ignored')
3034	Type,instDict,insVary = GetInstParms(hId,Inst[0])
3035	controlDict[pfx+'histType'] = Type
3036	if 'XC' in Type:
3037	if pfx+'Lam1' in instDict:
3038	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3039	else:
3040	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
3041	histDict.update(instDict)
3042	histVary += insVary
3043
3044	Sample = Histogram['Sample Parameters']
3045	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3046	controlDict[pfx+'instType'] = Type
3047	histDict.update(sampDict)
3048	histVary += sampVary
3049
3050
3051	if Print:
3052	pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3053	pFile.write(135*'='+'\n')
3054	Units = {'C':' deg','T':' msec'}
3055	units = Units[controlDict[pfx+'histType'][2]]
3056	Limits = controlDict[pfx+'Limits']
3057	pFile.write(' Instrument type: %s\n'%Sample['Type'])
3058	pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3059	if len(controlDict[pfx+'Exclude']):
3060	excls = controlDict[pfx+'Exclude']
3061	for excl in excls:
3062	pFile.write(' Excluded region: %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units))
3063	PrintSampleParms(Sample)
3064	PrintInstParms(Inst[0])
3065	PrintBackground(Background)
3066	elif 'HKLF' in histogram:
3067	Histogram = Histograms[histogram]
3068	hId = Histogram['hId']
3069	pfx = ':'+str(hId)+':'
3070	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3071	Inst = Histogram['Instrument Parameters'][0]
3072	controlDict[pfx+'histType'] = Inst['Type'][0]
3073	if 'X' in Inst['Type'][0]:
3074	histDict[pfx+'Lam'] = Inst['Lam'][1]
3075	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3076	return histVary,histDict,controlDict
3077
3078	def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3079	'needs a doc string'
3080
3081	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3082	Back = Background[0]
3083	DebyePeaks = Background[1]
3084	lenBack = len(Back[3:])
3085	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
3086	for i in range(lenBack):
3087	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3088	if pfx+'Back;'+str(i) in sigDict:
3089	backSig[i] = sigDict[pfx+'Back;'+str(i)]
3090	if DebyePeaks['nDebye']:
3091	for i in range(DebyePeaks['nDebye']):
3092	names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3093	for j,name in enumerate(names):
3094	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3095	if name in sigDict:
3096	backSig[lenBack+3*i+j] = sigDict[name]
3097	if DebyePeaks['nPeaks']:
3098	for i in range(DebyePeaks['nPeaks']):
3099	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3100	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3101	for j,name in enumerate(names):
3102	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3103	if name in sigDict:
3104	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
3105	return backSig
3106
3107	def SetInstParms(pfx,Inst,parmDict,sigDict):
3108	instSig = {}
3109	insKeys = list(Inst.keys())
3110	insKeys.sort()
3111	for item in insKeys:
3112	insName = pfx+item
3113	Inst[item][1] = parmDict[insName]
3114	if insName in sigDict:
3115	instSig[item] = sigDict[insName]
3116	else:
3117	instSig[item] = 0
3118	return instSig
3119
3120	def SetSampleParms(pfx,Sample,parmDict,sigDict):
3121	if 'Bragg' in Sample['Type']: #Bragg-Brentano
3122	sampSig = [0 for i in range(5)]
3123	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3124	Sample[item][0] = parmDict[pfx+item]
3125	if pfx+item in sigDict:
3126	sampSig[i] = sigDict[pfx+item]
3127	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3128	sampSig = [0 for i in range(4)]
3129	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3130	Sample[item][0] = parmDict[pfx+item]
3131	if pfx+item in sigDict:
3132	sampSig[i] = sigDict[pfx+item]
3133	return sampSig
3134
3135	def PrintBackgroundSig(Background,backSig):
3136	Back = Background[0]
3137	DebyePeaks = Background[1]
3138	valstr = ' value : '
3139	sigstr = ' sig : '
3140	refine = False
3141	for i,back in enumerate(Back[3:]):
3142	valstr += '%10.4g'%(back)
3143	if Back[1]:
3144	refine = True
3145	sigstr += '%10.4g'%(backSig[i])
3146	else:
3147	sigstr += 10*' '
3148	if refine:
3149	pFile.write('\n Background function: %s\n'%Back[0])
3150	pFile.write(valstr+'\n')
3151	pFile.write(sigstr+'\n')
3152	if DebyePeaks['nDebye']:
3153	ifAny = False
3154	ptfmt = "%12.3f"
3155	names = ' names :'
3156	ptstr = ' values:'
3157	sigstr = ' esds :'
3158	for item in sigDict:
3159	if 'Debye' in item:
3160	ifAny = True
3161	names += '%12s'%(item)
3162	ptstr += ptfmt%(parmDict[item])
3163	sigstr += ptfmt%(sigDict[item])
3164	if ifAny:
3165	pFile.write('\n Debye diffuse scattering coefficients\n')
3166	pFile.write(names+'\n')
3167	pFile.write(ptstr+'\n')
3168	pFile.write(sigstr+'\n')
3169	if DebyePeaks['nPeaks']:
3170	pFile.write('\n Single peak coefficients:\n')
3171	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3172	line = ' peak no. '
3173	for parm in parms:
3174	line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3175	pFile.write(line+'\n')
3176	for ip in range(DebyePeaks['nPeaks']):
3177	ptstr = ' %4d '%(ip)
3178	for parm in parms:
3179	fmt = '%14.3f'
3180	efmt = '%12.3f'
3181	if parm == 'BkPkpos':
3182	fmt = '%14.4f'
3183	efmt = '%12.4f'
3184	name = pfx+parm+';%d'%(ip)
3185	ptstr += fmt%(parmDict[name])
3186	if name in sigDict:
3187	ptstr += efmt%(sigDict[name])
3188	else:
3189	ptstr += 12*' '
3190	pFile.write(ptstr+'\n')
3191	sumBk = np.array(Histogram['sumBk'])
3192	pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3193	(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3194
3195	def PrintInstParmsSig(Inst,instSig):
3196	refine = False
3197	insKeys = list(instSig.keys())
3198	insKeys.sort()
3199	iBeg = 0
3200	Ok = True
3201	while Ok:
3202	ptlbls = ' names :'
3203	ptstr = ' value :'
3204	sigstr = ' sig :'
3205	iFin = min(iBeg+9,len(insKeys))
3206	for name in insKeys[iBeg:iFin]:
3207	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3208	ptlbls += '%12s' % (name)
3209	ptstr += '%12.6f' % (Inst[name][1])
3210	if instSig[name]:
3211	refine = True
3212	sigstr += '%12.6f' % (instSig[name])
3213	else:
3214	sigstr += 12*' '
3215	if refine:
3216	pFile.write('\n Instrument Parameters:\n')
3217	pFile.write(ptlbls+'\n')
3218	pFile.write(ptstr+'\n')
3219	pFile.write(sigstr+'\n')
3220	iBeg = iFin
3221	if iBeg == len(insKeys):
3222	Ok = False
3223
3224	def PrintSampleParmsSig(Sample,sampleSig):
3225	ptlbls = ' names :'
3226	ptstr = ' values:'
3227	sigstr = ' sig :'
3228	refine = False
3229	if 'Bragg' in Sample['Type']:
3230	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3231	ptlbls += '%14s'%(item)
3232	ptstr += '%14.4f'%(Sample[item][0])
3233	if sampleSig[i]:
3234	refine = True
3235	sigstr += '%14.4f'%(sampleSig[i])
3236	else:
3237	sigstr += 14*' '
3238
3239	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3240	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3241	ptlbls += '%14s'%(item)
3242	ptstr += '%14.4f'%(Sample[item][0])
3243	if sampleSig[i]:
3244	refine = True
3245	sigstr += '%14.4f'%(sampleSig[i])
3246	else:
3247	sigstr += 14*' '
3248
3249	if refine:
3250	pFile.write('\n Sample Parameters:\n')
3251	pFile.write(ptlbls+'\n')
3252	pFile.write(ptstr+'\n')
3253	pFile.write(sigstr+'\n')
3254
3255	histoList = list(Histograms.keys())
3256	histoList.sort()
3257	for histogram in histoList:
3258	if 'PWDR' in histogram:
3259	Histogram = Histograms[histogram]
3260	hId = Histogram['hId']
3261	pfx = ':'+str(hId)+':'
3262	Background = Histogram['Background']
3263	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3264
3265	Inst = Histogram['Instrument Parameters'][0]
3266	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3267
3268	Sample = Histogram['Sample Parameters']
3269	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3270
3271	pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3272	pFile.write(135*'='+'\n')
3273	pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3274	pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3275	(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3276	pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3277	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3278	if Print:
3279	pFile.write(' Instrument type: %s\n'%Sample['Type'])
3280	if FFtables != None and 'N' not in Inst['Type'][0]:
3281	PrintFprime(FFtables,pfx,pFile)
3282	PrintSampleParmsSig(Sample,sampSig)
3283	PrintInstParmsSig(Inst,instSig)
3284	PrintBackgroundSig(Background,backSig)
3285

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

Download in other formats: