Context Navigation

source: trunk/GSASIIstrIO.py @ 3340

Last change on this file since 3340 was 3340, checked in by vondreele, 6 years ago
make testDeriv py3 compatible fix problem with Reload draw atoms - mag moments add \|Mag\| to lst output of magnetic moments correct errors in mag moment structure factor & derivatives
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 149.3 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2018-04-11 18:32:49 +0000 (Wed, 11 Apr 2018) $
4	# $Author: vondreele $
5	# $Revision: 3340 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 3340 2018-04-11 18:32:49Z 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: 3340 $")
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	cell = General['Cell'][1:7]
1659	G = G2lat.fillgmat(cell)
1660	ast = np.sqrt(np.diag(G))
1661	GS = G/np.outer(ast,ast)
1662	pFile.write('\n Magnetic Moments:\n') #add magnitude & angle, etc.? TBD
1663	line = ' name Mx My Mz \|Mag\|'
1664	cx,ct,cs,cia = General['AtomPtrs']
1665	cmx = cx+4
1666	AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1667	pFile.write(line+'\n')
1668	pFile.write(135*'-'+'\n')
1669	fmt = {0:'%7s',ct:'%7s',cmx:'%10.3f',cmx+1:'%10.3f',cmx+2:'%10.3f'}
1670	noFXsig = {cmx:[10' ','%10s'],cmx+1:[10' ','%10s'],cmx+2:[10*' ','%10s']}
1671	for i,at in enumerate(Atoms):
1672	if AtInfo[at[ct]]:
1673	name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1674	valstr = ' values:'
1675	sigstr = ' sig :'
1676	for ind in range(cmx,cmx+3):
1677	sigind = str(i)+':'+str(ind)
1678	valstr += fmt[ind]%(at[ind])
1679	if sigind in atomsSig:
1680	sigstr += fmt[ind]%(atomsSig[sigind])
1681	else:
1682	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1683	mag = np.array(at[cmx:cmx+3])
1684	Mag = np.sqrt(np.inner(mag,np.inner(mag,GS)))
1685	valstr += '%10.3f'%Mag
1686	sigstr += 10*' '
1687	pFile.write(name+'\n')
1688	pFile.write(valstr+'\n')
1689	pFile.write(sigstr+'\n')
1690
1691	def PrintWavesAndSig(General,Atoms,wavesSig):
1692	cx,ct,cs,cia = General['AtomPtrs']
1693	pFile.write('\n Modulation waves\n')
1694	names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1695	'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1696	'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1697	pFile.write(135*'-'+'\n')
1698	for i,at in enumerate(Atoms):
1699	AtomSS = at[-1]['SS1']
1700	for Stype in ['Sfrac','Spos','Sadp','Smag']:
1701	Waves = AtomSS[Stype]
1702	if len(Waves):
1703	waveType = Waves[0]
1704	else:
1705	continue
1706	if len(Waves):
1707	pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
1708	(at[ct-1],at[cs],Stype,waveType))
1709	for iw,wave in enumerate(Waves[1:]):
1710	stiw = ':'+str(i)+':'+str(iw)
1711	namstr = ' names :'
1712	valstr = ' values:'
1713	sigstr = ' esds :'
1714	if Stype == 'Spos':
1715	nt = 6
1716	ot = 0
1717	if waveType in ['ZigZag','Block',] and not iw:
1718	nt = 5
1719	ot = 6
1720	for j in range(nt):
1721	name = names['Spos'][j+ot]
1722	namstr += '%12s'%(name)
1723	valstr += '%12.4f'%(wave[0][j])
1724	if name+stiw in wavesSig:
1725	sigstr += '%12.4f'%(wavesSig[name+stiw])
1726	else:
1727	sigstr += 12*' '
1728	elif Stype == 'Sfrac':
1729	ot = 0
1730	if 'Crenel' in waveType and not iw:
1731	ot = 2
1732	for j in range(2):
1733	name = names['Sfrac'][j+ot]
1734	namstr += '%12s'%(names['Sfrac'][j+ot])
1735	valstr += '%12.4f'%(wave[0][j])
1736	if name+stiw in wavesSig:
1737	sigstr += '%12.4f'%(wavesSig[name+stiw])
1738	else:
1739	sigstr += 12*' '
1740	elif Stype == 'Sadp':
1741	for j in range(12):
1742	name = names['Sadp'][j]
1743	namstr += '%10s'%(names['Sadp'][j])
1744	valstr += '%10.6f'%(wave[0][j])
1745	if name+stiw in wavesSig:
1746	sigstr += '%10.6f'%(wavesSig[name+stiw])
1747	else:
1748	sigstr += 10*' '
1749	elif Stype == 'Smag':
1750	for j in range(6):
1751	name = names['Smag'][j]
1752	namstr += '%12s'%(names['Smag'][j])
1753	valstr += '%12.4f'%(wave[0][j])
1754	if name+stiw in wavesSig:
1755	sigstr += '%12.4f'%(wavesSig[name+stiw])
1756	else:
1757	sigstr += 12*' '
1758
1759	pFile.write(namstr+'\n')
1760	pFile.write(valstr+'\n')
1761	pFile.write(sigstr+'\n')
1762
1763
1764	def PrintRBObjPOAndSig(rbfx,rbsx):
1765	namstr = ' names :'
1766	valstr = ' values:'
1767	sigstr = ' esds :'
1768	for i,px in enumerate(['Px:','Py:','Pz:']):
1769	name = pfx+rbfx+px+rbsx
1770	namstr += '%12s'%('Pos '+px[1])
1771	valstr += '%12.5f'%(parmDict[name])
1772	if name in sigDict:
1773	sigstr += '%12.5f'%(sigDict[name])
1774	else:
1775	sigstr += 12*' '
1776	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1777	name = pfx+rbfx+po+rbsx
1778	namstr += '%12s'%('Ori '+po[1])
1779	valstr += '%12.5f'%(parmDict[name])
1780	if name in sigDict:
1781	sigstr += '%12.5f'%(sigDict[name])
1782	else:
1783	sigstr += 12*' '
1784	pFile.write(namstr+'\n')
1785	pFile.write(valstr+'\n')
1786	pFile.write(sigstr+'\n')
1787
1788	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1789	namstr = ' names :'
1790	valstr = ' values:'
1791	sigstr = ' esds :'
1792	if 'N' not in TLS:
1793	pFile.write(' Thermal motion:\n')
1794	if 'T' in TLS:
1795	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1796	name = pfx+rbfx+pt+rbsx
1797	namstr += '%12s'%(pt[:3])
1798	valstr += '%12.5f'%(parmDict[name])
1799	if name in sigDict:
1800	sigstr += '%12.5f'%(sigDict[name])
1801	else:
1802	sigstr += 12*' '
1803	pFile.write(namstr+'\n')
1804	pFile.write(valstr+'\n')
1805	pFile.write(sigstr+'\n')
1806	if 'L' in TLS:
1807	namstr = ' names :'
1808	valstr = ' values:'
1809	sigstr = ' esds :'
1810	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1811	name = pfx+rbfx+pt+rbsx
1812	namstr += '%12s'%(pt[:3])
1813	valstr += '%12.3f'%(parmDict[name])
1814	if name in sigDict:
1815	sigstr += '%12.3f'%(sigDict[name])
1816	else:
1817	sigstr += 12*' '
1818	pFile.write(namstr+'\n')
1819	pFile.write(valstr+'\n')
1820	pFile.write(sigstr+'\n')
1821	if 'S' in TLS:
1822	namstr = ' names :'
1823	valstr = ' values:'
1824	sigstr = ' esds :'
1825	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1826	name = pfx+rbfx+pt+rbsx
1827	namstr += '%12s'%(pt[:3])
1828	valstr += '%12.4f'%(parmDict[name])
1829	if name in sigDict:
1830	sigstr += '%12.4f'%(sigDict[name])
1831	else:
1832	sigstr += 12*' '
1833	pFile.write(namstr+'\n')
1834	pFile.write(valstr+'\n')
1835	pFile.write(sigstr+'\n')
1836	if 'U' in TLS:
1837	name = pfx+rbfx+'U:'+rbsx
1838	namstr = ' names :'+'%12s'%('Uiso')
1839	valstr = ' values:'+'%12.5f'%(parmDict[name])
1840	if name in sigDict:
1841	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1842	else:
1843	sigstr = ' esds :'+12*' '
1844	pFile.write(namstr+'\n')
1845	pFile.write(valstr+'\n')
1846	pFile.write(sigstr+'\n')
1847
1848	def PrintRBObjTorAndSig(rbsx):
1849	namstr = ' names :'
1850	valstr = ' values:'
1851	sigstr = ' esds :'
1852	nTors = len(RBObj['Torsions'])
1853	if nTors:
1854	pFile.write(' Torsions:\n')
1855	for it in range(nTors):
1856	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1857	namstr += '%12s'%('Tor'+str(it))
1858	valstr += '%12.4f'%(parmDict[name])
1859	if name in sigDict:
1860	sigstr += '%12.4f'%(sigDict[name])
1861	pFile.write(namstr+'\n')
1862	pFile.write(valstr+'\n')
1863	pFile.write(sigstr+'\n')
1864
1865	def PrintSHtextureAndSig(textureData,SHtextureSig):
1866	pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
1867	names = ['omega','chi','phi']
1868	namstr = ' names :'
1869	ptstr = ' values:'
1870	sigstr = ' esds :'
1871	for name in names:
1872	namstr += '%12s'%(name)
1873	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1874	if 'Sample '+name in SHtextureSig:
1875	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1876	else:
1877	sigstr += 12*' '
1878	pFile.write(namstr+'\n')
1879	pFile.write(ptstr+'\n')
1880	pFile.write(sigstr+'\n')
1881	pFile.write('\n Texture coefficients:\n')
1882	SHcoeff = textureData['SH Coeff'][1]
1883	SHkeys = list(SHcoeff.keys())
1884	nCoeff = len(SHcoeff)
1885	nBlock = nCoeff//10+1
1886	iBeg = 0
1887	iFin = min(iBeg+10,nCoeff)
1888	for block in range(nBlock):
1889	namstr = ' names :'
1890	ptstr = ' values:'
1891	sigstr = ' esds :'
1892	for name in SHkeys[iBeg:iFin]:
1893	namstr += '%12s'%(name)
1894	ptstr += '%12.3f'%(SHcoeff[name])
1895	if name in SHtextureSig:
1896	sigstr += '%12.3f'%(SHtextureSig[name])
1897	else:
1898	sigstr += 12*' '
1899	pFile.write(namstr+'\n')
1900	pFile.write(ptstr+'\n')
1901	pFile.write(sigstr+'\n')
1902	iBeg += 10
1903	iFin = min(iBeg+10,nCoeff)
1904
1905	##########################################################################
1906	# SetPhaseData starts here
1907	pFile.write('\n Phases:\n')
1908	for phase in Phases:
1909	pFile.write(' Result for phase: %s\n'%phase)
1910	pFile.write(135*'='+'\n')
1911	Phase = Phases[phase]
1912	General = Phase['General']
1913	SGData = General['SGData']
1914	Atoms = Phase['Atoms']
1915	AtLookup = []
1916	if Atoms and not General.get('doPawley'):
1917	cx,ct,cs,cia = General['AtomPtrs']
1918	AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1919	cell = General['Cell']
1920	pId = Phase['pId']
1921	pfx = str(pId)+'::'
1922	if cell[0]:
1923	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1924	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1925	pFile.write(' Reciprocal metric tensor: \n')
1926	ptfmt = "%15.9f"
1927	names = ['A11','A22','A33','A12','A13','A23']
1928	namstr = ' names :'
1929	ptstr = ' values:'
1930	sigstr = ' esds :'
1931	for name,a,siga in zip(names,A,sigA):
1932	namstr += '%15s'%(name)
1933	ptstr += ptfmt%(a)
1934	if siga:
1935	sigstr += ptfmt%(siga)
1936	else:
1937	sigstr += 15*' '
1938	pFile.write(namstr+'\n')
1939	pFile.write(ptstr+'\n')
1940	pFile.write(sigstr+'\n')
1941	cell[1:7] = G2lat.A2cell(A)
1942	cell[7] = G2lat.calc_V(A)
1943	pFile.write(' New unit cell:\n')
1944	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1945	names = ['a','b','c','alpha','beta','gamma','Volume']
1946	namstr = ' names :'
1947	ptstr = ' values:'
1948	sigstr = ' esds :'
1949	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1950	namstr += '%12s'%(name)
1951	ptstr += fmt%(a)
1952	if siga:
1953	sigstr += fmt%(siga)
1954	else:
1955	sigstr += 12*' '
1956	pFile.write(namstr+'\n')
1957	pFile.write(ptstr+'\n')
1958	pFile.write(sigstr+'\n')
1959	ik = 6 #for Pawley stuff below
1960	if General.get('Modulated',False):
1961	ik = 7
1962	Vec,vRef,maxH = General['SuperVec']
1963	if vRef:
1964	pFile.write(' New modulation vector:\n')
1965	namstr = ' names :'
1966	ptstr = ' values:'
1967	sigstr = ' esds :'
1968	for var in ['mV0','mV1','mV2']:
1969	namstr += '%12s'%(pfx+var)
1970	ptstr += '%12.6f'%(parmDict[pfx+var])
1971	if pfx+var in sigDict:
1972	sigstr += '%12.6f'%(sigDict[pfx+var])
1973	else:
1974	sigstr += 12*' '
1975	pFile.write(namstr+'\n')
1976	pFile.write(ptstr+'\n')
1977	pFile.write(sigstr+'\n')
1978
1979	General['Mass'] = 0.
1980	if Phase['General'].get('doPawley'):
1981	pawleyRef = Phase['Pawley ref']
1982	for i,refl in enumerate(pawleyRef):
1983	key = pfx+'PWLref:'+str(i)
1984	refl[ik] = parmDict[key]
1985	if key in sigDict:
1986	refl[ik+1] = sigDict[key]
1987	else:
1988	refl[ik+1] = 0
1989	else:
1990	VRBIds = RBIds['Vector']
1991	RRBIds = RBIds['Residue']
1992	RBModels = Phase['RBModels']
1993	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1994	jrb = VRBIds.index(RBObj['RBId'])
1995	rbsx = str(irb)+':'+str(jrb)
1996	pFile.write(' Vector rigid body parameters:\n')
1997	PrintRBObjPOAndSig('RBV',rbsx)
1998	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1999	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2000	jrb = RRBIds.index(RBObj['RBId'])
2001	rbsx = str(irb)+':'+str(jrb)
2002	pFile.write(' Residue rigid body parameters:\n')
2003	PrintRBObjPOAndSig('RBR',rbsx)
2004	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2005	PrintRBObjTorAndSig(rbsx)
2006	atomsSig = {}
2007	wavesSig = {}
2008	cx,ct,cs,cia = General['AtomPtrs']
2009	for i,at in enumerate(Atoms):
2010	names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2011	cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2012	cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2013	cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2014	for ind in range(cx,cx+4):
2015	at[ind] = parmDict[names[ind]]
2016	if ind in range(cx,cx+3):
2017	name = names[ind].replace('A','dA')
2018	else:
2019	name = names[ind]
2020	if name in sigDict:
2021	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2022	if at[cia] == 'I':
2023	at[cia+1] = parmDict[names[cia+1]]
2024	if names[cia+1] in sigDict:
2025	atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2026	else:
2027	for ind in range(cia+2,cia+8):
2028	at[ind] = parmDict[names[ind]]
2029	if names[ind] in sigDict:
2030	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2031	if General['Type'] == 'magnetic':
2032	for ind in range(cx+4,cx+7):
2033	at[ind] = parmDict[names[ind]]
2034	if names[ind] in sigDict:
2035	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2036	ind = General['AtomTypes'].index(at[ct])
2037	General['Mass'] += General['AtomMass'][ind]at[cx+3]at[cx+5]
2038	if General.get('Modulated',False):
2039	AtomSS = at[-1]['SS1']
2040	for Stype in ['Sfrac','Spos','Sadp','Smag']:
2041	Waves = AtomSS[Stype]
2042	if len(Waves):
2043	waveType = Waves[0]
2044	else:
2045	continue
2046	for iw,wave in enumerate(Waves[1:]):
2047	stiw = str(i)+':'+str(iw)
2048	if Stype == 'Spos':
2049	if waveType in ['ZigZag','Block',] and not iw:
2050	names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2051	else:
2052	names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2053	'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2054	elif Stype == 'Sadp':
2055	names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2056	'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2057	'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2058	'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2059	elif Stype == 'Sfrac':
2060	if 'Crenel' in waveType and not iw:
2061	names = ['Fzero:'+stiw,'Fwid:'+stiw]
2062	else:
2063	names = ['Fsin:'+stiw,'Fcos:'+stiw]
2064	elif Stype == 'Smag':
2065	names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2066	'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2067	for iname,name in enumerate(names):
2068	AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2069	if pfx+name in sigDict:
2070	wavesSig[name] = sigDict[pfx+name]
2071
2072	PrintAtomsAndSig(General,Atoms,atomsSig)
2073	if General['Type'] == 'magnetic':
2074	PrintMomentsAndSig(General,Atoms,atomsSig)
2075	if General.get('Modulated',False):
2076	PrintWavesAndSig(General,Atoms,wavesSig)
2077
2078
2079	textureData = General['SH Texture']
2080	if textureData['Order']:
2081	SHtextureSig = {}
2082	for name in ['omega','chi','phi']:
2083	aname = pfx+'SH '+name
2084	textureData['Sample '+name][1] = parmDict[aname]
2085	if aname in sigDict:
2086	SHtextureSig['Sample '+name] = sigDict[aname]
2087	for name in textureData['SH Coeff'][1]:
2088	aname = pfx+name
2089	textureData['SH Coeff'][1][name] = parmDict[aname]
2090	if aname in sigDict:
2091	SHtextureSig[name] = sigDict[aname]
2092	PrintSHtextureAndSig(textureData,SHtextureSig)
2093	if phase in RestraintDict and not Phase['General'].get('doPawley'):
2094	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2095	textureData,RestraintDict[phase],pFile)
2096
2097	################################################################################
2098	##### Histogram & Phase data
2099	################################################################################
2100
2101	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2102	'''Loads the HAP histogram/phase information into dicts
2103
2104	:param dict Phases: phase information
2105	:param dict Histograms: Histogram information
2106	:param bool Print: prints information as it is read
2107	:param file pFile: file object to print to (the default, None causes printing to the console)
2108	:param bool resetRefList: Should the contents of the Reflection List be initialized
2109	on loading. The default, True, initializes the Reflection List as it is loaded.
2110
2111	:returns: (hapVary,hapDict,controlDict)
2112	* hapVary: list of refined variables
2113	* hapDict: dict with refined variables and their values
2114	* controlDict: dict with fixed parameters
2115	'''
2116
2117	def PrintSize(hapData):
2118	if hapData[0] in ['isotropic','uniaxial']:
2119	line = '\n Size model : %9s'%(hapData[0])
2120	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2121	if hapData[0] == 'uniaxial':
2122	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2123	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2124	pFile.write(line+'\n')
2125	else:
2126	pFile.write('\n Size model : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2127	(hapData[0],hapData[1][2],hapData[2][2]))
2128	Snames = ['S11','S22','S33','S12','S13','S23']
2129	ptlbls = ' names :'
2130	ptstr = ' values:'
2131	varstr = ' refine:'
2132	for i,name in enumerate(Snames):
2133	ptlbls += '%12s' % (name)
2134	ptstr += '%12.3f' % (hapData[4][i])
2135	varstr += '%12s' % (str(hapData[5][i]))
2136	pFile.write(ptlbls+'\n')
2137	pFile.write(ptstr+'\n')
2138	pFile.write(varstr+'\n')
2139
2140	def PrintMuStrain(hapData,SGData):
2141	if hapData[0] in ['isotropic','uniaxial']:
2142	line = '\n Mustrain model: %9s'%(hapData[0])
2143	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2144	if hapData[0] == 'uniaxial':
2145	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2146	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2147	pFile.write(line+'\n')
2148	else:
2149	pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2150	(hapData[0],hapData[1][2],hapData[2][2]))
2151	Snames = G2spc.MustrainNames(SGData)
2152	ptlbls = ' names :'
2153	ptstr = ' values:'
2154	varstr = ' refine:'
2155	for i,name in enumerate(Snames):
2156	ptlbls += '%12s' % (name)
2157	ptstr += '%12.1f' % (hapData[4][i])
2158	varstr += '%12s' % (str(hapData[5][i]))
2159	pFile.write(ptlbls+'\n')
2160	pFile.write(ptstr+'\n')
2161	pFile.write(varstr+'\n')
2162
2163	def PrintHStrain(hapData,SGData):
2164	pFile.write('\n Hydrostatic/elastic strain:\n')
2165	Hsnames = G2spc.HStrainNames(SGData)
2166	ptlbls = ' names :'
2167	ptstr = ' values:'
2168	varstr = ' refine:'
2169	for i,name in enumerate(Hsnames):
2170	ptlbls += '%12s' % (name)
2171	ptstr += '%12.4g' % (hapData[0][i])
2172	varstr += '%12s' % (str(hapData[1][i]))
2173	pFile.write(ptlbls+'\n')
2174	pFile.write(ptstr+'\n')
2175	pFile.write(varstr+'\n')
2176
2177	def PrintSHPO(hapData):
2178	pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2179	ptlbls = ' names :'
2180	ptstr = ' values:'
2181	for item in hapData[5]:
2182	ptlbls += '%12s'%(item)
2183	ptstr += '%12.3f'%(hapData[5][item])
2184	pFile.write(ptlbls+'\n')
2185	pFile.write(ptstr+'\n')
2186
2187	def PrintBabinet(hapData):
2188	pFile.write('\n Babinet form factor modification:\n')
2189	ptlbls = ' names :'
2190	ptstr = ' values:'
2191	varstr = ' refine:'
2192	for name in ['BabA','BabU']:
2193	ptlbls += '%12s' % (name)
2194	ptstr += '%12.6f' % (hapData[name][0])
2195	varstr += '%12s' % (str(hapData[name][1]))
2196	pFile.write(ptlbls+'\n')
2197	pFile.write(ptstr+'\n')
2198	pFile.write(varstr+'\n')
2199
2200	hapDict = {}
2201	hapVary = []
2202	controlDict = {}
2203
2204	for phase in Phases:
2205	HistoPhase = Phases[phase]['Histograms']
2206	SGData = Phases[phase]['General']['SGData']
2207	cell = Phases[phase]['General']['Cell'][1:7]
2208	A = G2lat.cell2A(cell)
2209	if Phases[phase]['General'].get('Modulated',False):
2210	SSGData = Phases[phase]['General']['SSGData']
2211	Vec,x,maxH = Phases[phase]['General']['SuperVec']
2212	pId = Phases[phase]['pId']
2213	histoList = list(HistoPhase.keys())
2214	histoList.sort()
2215	for histogram in histoList:
2216	try:
2217	Histogram = Histograms[histogram]
2218	except KeyError:
2219	#skip if histogram not included e.g. in a sequential refinement
2220	continue
2221	if not HistoPhase[histogram]['Use']: #remove previously created & now unused phase reflection list
2222	if phase in Histograms[histogram]['Reflection Lists']:
2223	del Histograms[histogram]['Reflection Lists'][phase]
2224	continue
2225	hapData = HistoPhase[histogram]
2226	hId = Histogram['hId']
2227	if 'PWDR' in histogram:
2228	limits = Histogram['Limits'][1]
2229	inst = Histogram['Instrument Parameters'][0] #TODO - grab table here if present
2230	if 'C' in inst['Type'][1]:
2231	try:
2232	wave = inst['Lam'][1]
2233	except KeyError:
2234	wave = inst['Lam1'][1]
2235	dmin = wave/(2.0*sind(limits[1]/2.0))
2236	elif 'T' in inst['Type'][0]:
2237	dmin = limits[0]/inst['difC'][1]
2238	pfx = str(pId)+':'+str(hId)+':'
2239	if Phases[phase]['General']['doPawley']:
2240	hapDict[pfx+'LeBail'] = False #Pawley supercedes LeBail
2241	hapDict[pfx+'newLeBail'] = True
2242	Tmin = G2lat.Dsp2pos(inst,dmin)
2243	if 'C' in inst['Type'][1]:
2244	limits[1] = min(limits[1],Tmin)
2245	else:
2246	limits[0] = max(limits[0],Tmin)
2247	else:
2248	hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2249	hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2250	if Phases[phase]['General']['Type'] == 'magnetic':
2251	dmin = max(dmin,Phases[phase]['General']['MagDmin'])
2252	for item in ['Scale','Extinction']:
2253	hapDict[pfx+item] = hapData[item][0]
2254	if hapData[item][1] and not hapDict[pfx+'LeBail']:
2255	hapVary.append(pfx+item)
2256	names = G2spc.HStrainNames(SGData)
2257	HSvals = []
2258	for i,name in enumerate(names):
2259	hapDict[pfx+name] = hapData['HStrain'][0][i]
2260	HSvals.append(hapDict[pfx+name])
2261	if hapData['HStrain'][1][i]:
2262	# if hapData['HStrain'][1][i] and not hapDict[pfx+'LeBail']:
2263	hapVary.append(pfx+name)
2264	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2265	if hapData['Pref.Ori.'][0] == 'MD':
2266	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2267	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2268	if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2269	hapVary.append(pfx+'MD')
2270	else: #'SH' spherical harmonics
2271	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2272	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2273	controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2274	controlDict[pfx+'SHhkl'] = []
2275	try: #patch for old Pref.Ori. items
2276	controlDict[pfx+'SHtoler'] = 0.1
2277	if hapData['Pref.Ori.'][6][0] != '':
2278	controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2279	controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2280	except IndexError:
2281	pass
2282	for item in hapData['Pref.Ori.'][5]:
2283	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2284	if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2285	hapVary.append(pfx+item)
2286	for item in ['Mustrain','Size']:
2287	controlDict[pfx+item+'Type'] = hapData[item][0]
2288	hapDict[pfx+item+';mx'] = hapData[item][1][2]
2289	if hapData[item][2][2]:
2290	hapVary.append(pfx+item+';mx')
2291	if hapData[item][0] in ['isotropic','uniaxial']:
2292	hapDict[pfx+item+';i'] = hapData[item][1][0]
2293	if hapData[item][2][0]:
2294	hapVary.append(pfx+item+';i')
2295	if hapData[item][0] == 'uniaxial':
2296	controlDict[pfx+item+'Axis'] = hapData[item][3]
2297	hapDict[pfx+item+';a'] = hapData[item][1][1]
2298	if hapData[item][2][1]:
2299	hapVary.append(pfx+item+';a')
2300	else: #generalized for mustrain or ellipsoidal for size
2301	Nterms = len(hapData[item][4])
2302	if item == 'Mustrain':
2303	names = G2spc.MustrainNames(SGData)
2304	pwrs = []
2305	for name in names:
2306	h,k,l = name[1:]
2307	pwrs.append([int(h),int(k),int(l)])
2308	controlDict[pfx+'MuPwrs'] = pwrs
2309	for i in range(Nterms):
2310	sfx = ';'+str(i)
2311	hapDict[pfx+item+sfx] = hapData[item][4][i]
2312	if hapData[item][5][i]:
2313	hapVary.append(pfx+item+sfx)
2314	if Phases[phase]['General']['Type'] != 'magnetic':
2315	for bab in ['BabA','BabU']:
2316	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2317	if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2318	hapVary.append(pfx+bab)
2319
2320	if Print:
2321	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2322	pFile.write(135*'='+'\n')
2323	if hapDict[pfx+'LeBail']:
2324	pFile.write(' Perform LeBail extraction\n')
2325	else:
2326	pFile.write(' Phase fraction : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2327	pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2328	if hapData['Pref.Ori.'][0] == 'MD':
2329	Ax = hapData['Pref.Ori.'][3]
2330	pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2331	(hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2332	else: #'SH' for spherical harmonics
2333	PrintSHPO(hapData['Pref.Ori.'])
2334	pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2335	PrintSize(hapData['Size'])
2336	PrintMuStrain(hapData['Mustrain'],SGData)
2337	PrintHStrain(hapData['HStrain'],SGData)
2338	if Phases[phase]['General']['Type'] != 'magnetic':
2339	if hapData['Babinet']['BabA'][0]:
2340	PrintBabinet(hapData['Babinet'])
2341	if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData['LeBail'] and hapData['newLeBail'])):
2342	if hapData.get('LeBail',True): #stop regeneating reflections for LeBail
2343	hapData['newLeBail'] = False
2344	refList = []
2345	Uniq = []
2346	Phi = []
2347	useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2348	if Phases[phase]['General'].get('Modulated',False):
2349	ifSuper = True
2350	HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2351	HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2352	for h,k,l,m,d in HKLd:
2353	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2354	mul *= 2 # for powder overlap of Friedel pairs
2355	if m or not ext or useExt:
2356	if 'C' in inst['Type'][0]:
2357	pos = G2lat.Dsp2pos(inst,d)
2358	if limits[0] < pos < limits[1]:
2359	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])
2360	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2361	Uniq.append(uniq)
2362	Phi.append(phi)
2363	elif 'T' in inst['Type'][0]:
2364	pos = G2lat.Dsp2pos(inst,d)
2365	if limits[0] < pos < limits[1]:
2366	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2367	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])
2368	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2369	#TODO - if tabulated put alp & bet in here
2370	Uniq.append(uniq)
2371	Phi.append(phi)
2372	else:
2373	ifSuper = False
2374	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2375	HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2376	for h,k,l,d in HKLd:
2377	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2378	mul *= 2 # for powder overlap of Friedel pairs
2379	if ext and not useExt:
2380	continue
2381	if 'C' in inst['Type'][0]:
2382	pos = G2lat.Dsp2pos(inst,d)
2383	if limits[0] < pos < limits[1]:
2384	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])
2385	#... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2386	Uniq.append(uniq)
2387	Phi.append(phi)
2388	elif 'T' in inst['Type'][0]:
2389	pos = G2lat.Dsp2pos(inst,d)
2390	if limits[0] < pos < limits[1]:
2391	wave = inst['difC'][1]d/(252.816inst['fltPath'][0])
2392	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])
2393	# ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2394	Uniq.append(uniq)
2395	Phi.append(phi)
2396	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2397	elif 'HKLF' in histogram:
2398	inst = Histogram['Instrument Parameters'][0]
2399	hId = Histogram['hId']
2400	hfx = ':%d:'%(hId)
2401	for item in inst:
2402	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2403	hapDict[hfx+item] = inst[item][1]
2404	pfx = str(pId)+':'+str(hId)+':'
2405	hapDict[pfx+'Scale'] = hapData['Scale'][0]
2406	if hapData['Scale'][1]:
2407	hapVary.append(pfx+'Scale')
2408
2409	extApprox,extType,extParms = hapData['Extinction']
2410	controlDict[pfx+'EType'] = extType
2411	controlDict[pfx+'EApprox'] = extApprox
2412	if 'C' in inst['Type'][0]:
2413	controlDict[pfx+'Tbar'] = extParms['Tbar']
2414	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2415	if 'Primary' in extType:
2416	Ekey = ['Ep',]
2417	elif 'I & II' in extType:
2418	Ekey = ['Eg','Es']
2419	elif 'Secondary Type II' == extType:
2420	Ekey = ['Es',]
2421	elif 'Secondary Type I' == extType:
2422	Ekey = ['Eg',]
2423	else: #'None'
2424	Ekey = []
2425	for eKey in Ekey:
2426	hapDict[pfx+eKey] = extParms[eKey][0]
2427	if extParms[eKey][1]:
2428	hapVary.append(pfx+eKey)
2429	for bab in ['BabA','BabU']:
2430	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2431	if hapData['Babinet'][bab][1]:
2432	hapVary.append(pfx+bab)
2433	Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2434	if len(Twins) == 1:
2435	hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2436	if hapData.get('Flack',[0,False])[1]:
2437	hapVary.append(pfx+'Flack')
2438	sumTwFr = 0.
2439	controlDict[pfx+'TwinLaw'] = []
2440	controlDict[pfx+'TwinInv'] = []
2441	NTL = 0
2442	for it,twin in enumerate(Twins):
2443	if 'bool' in str(type(twin[0])):
2444	controlDict[pfx+'TwinInv'].append(twin[0])
2445	controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2446	else:
2447	NTL += 1
2448	controlDict[pfx+'TwinInv'].append(False)
2449	controlDict[pfx+'TwinLaw'].append(twin[0])
2450	if it:
2451	hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2452	sumTwFr += twin[1]
2453	else:
2454	hapDict[pfx+'TwinFr:0'] = twin[1][0]
2455	controlDict[pfx+'TwinNMN'] = twin[1][1]
2456	if Twins[0][1][1]:
2457	hapVary.append(pfx+'TwinFr:'+str(it))
2458	controlDict[pfx+'NTL'] = NTL
2459	#need to make constraint on TwinFr
2460	controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2461	if len(Twins) > 1: #force sum to unity
2462	hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2463	if Print:
2464	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2465	pFile.write(135*'='+'\n')
2466	pFile.write(' Scale factor : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2467	if extType != 'None':
2468	pFile.write(' Extinction Type: %15s approx: %10s\n'%(extType,extApprox))
2469	text = ' Parameters :'
2470	for eKey in Ekey:
2471	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2472	pFile.write(text+'\n')
2473	if hapData['Babinet']['BabA'][0]:
2474	PrintBabinet(hapData['Babinet'])
2475	if not SGData['SGInv'] and len(Twins) == 1:
2476	pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2477	if len(Twins) > 1:
2478	for it,twin in enumerate(Twins):
2479	if 'bool' in str(type(twin[0])):
2480	pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine?\n'%
2481	(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1])
2482	else:
2483	pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2484	(str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],Twins[0][1][1]))
2485
2486	Histogram['Reflection Lists'] = phase
2487
2488	return hapVary,hapDict,controlDict
2489
2490	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2491	'needs a doc string'
2492
2493	def PrintSizeAndSig(hapData,sizeSig):
2494	line = '\n Size model: %9s'%(hapData[0])
2495	refine = False
2496	if hapData[0] in ['isotropic','uniaxial']:
2497	line += ' equatorial:%12.5f'%(hapData[1][0])
2498	if sizeSig[0][0]:
2499	line += ', sig:%8.4f'%(sizeSig[0][0])
2500	refine = True
2501	if hapData[0] == 'uniaxial':
2502	line += ' axial:%12.4f'%(hapData[1][1])
2503	if sizeSig[0][1]:
2504	refine = True
2505	line += ', sig:%8.4f'%(sizeSig[0][1])
2506	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2507	if sizeSig[0][2]:
2508	refine = True
2509	line += ', sig:%8.4f'%(sizeSig[0][2])
2510	if refine:
2511	pFile.write(line+'\n')
2512	else:
2513	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2514	if sizeSig[0][2]:
2515	refine = True
2516	line += ', sig:%8.4f'%(sizeSig[0][2])
2517	Snames = ['S11','S22','S33','S12','S13','S23']
2518	ptlbls = ' name :'
2519	ptstr = ' value :'
2520	sigstr = ' sig :'
2521	for i,name in enumerate(Snames):
2522	ptlbls += '%12s' % (name)
2523	ptstr += '%12.6f' % (hapData[4][i])
2524	if sizeSig[1][i]:
2525	refine = True
2526	sigstr += '%12.6f' % (sizeSig[1][i])
2527	else:
2528	sigstr += 12*' '
2529	if refine:
2530	pFile.write(line+'\n')
2531	pFile.write(ptlbls+'\n')
2532	pFile.write(ptstr+'\n')
2533	pFile.write(sigstr+'\n')
2534
2535	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2536	line = '\n Mustrain model: %9s\n'%(hapData[0])
2537	refine = False
2538	if hapData[0] in ['isotropic','uniaxial']:
2539	line += ' equatorial:%12.1f'%(hapData[1][0])
2540	if mustrainSig[0][0]:
2541	line += ', sig:%8.1f'%(mustrainSig[0][0])
2542	refine = True
2543	if hapData[0] == 'uniaxial':
2544	line += ' axial:%12.1f'%(hapData[1][1])
2545	if mustrainSig[0][1]:
2546	line += ', sig:%8.1f'%(mustrainSig[0][1])
2547	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2548	if mustrainSig[0][2]:
2549	refine = True
2550	line += ', sig:%8.3f'%(mustrainSig[0][2])
2551	if refine:
2552	pFile.write(line+'\n')
2553	else:
2554	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2555	if mustrainSig[0][2]:
2556	refine = True
2557	line += ', sig:%8.3f'%(mustrainSig[0][2])
2558	Snames = G2spc.MustrainNames(SGData)
2559	ptlbls = ' name :'
2560	ptstr = ' value :'
2561	sigstr = ' sig :'
2562	for i,name in enumerate(Snames):
2563	ptlbls += '%12s' % (name)
2564	ptstr += '%12.1f' % (hapData[4][i])
2565	if mustrainSig[1][i]:
2566	refine = True
2567	sigstr += '%12.1f' % (mustrainSig[1][i])
2568	else:
2569	sigstr += 12*' '
2570	if refine:
2571	pFile.write(line+'\n')
2572	pFile.write(ptlbls+'\n')
2573	pFile.write(ptstr+'\n')
2574	pFile.write(sigstr+'\n')
2575
2576	def PrintHStrainAndSig(hapData,strainSig,SGData):
2577	Hsnames = G2spc.HStrainNames(SGData)
2578	ptlbls = ' name :'
2579	ptstr = ' value :'
2580	sigstr = ' sig :'
2581	refine = False
2582	for i,name in enumerate(Hsnames):
2583	ptlbls += '%12s' % (name)
2584	ptstr += '%12.4g' % (hapData[0][i])
2585	if name in strainSig:
2586	refine = True
2587	sigstr += '%12.4g' % (strainSig[name])
2588	else:
2589	sigstr += 12*' '
2590	if refine:
2591	pFile.write('\n Hydrostatic/elastic strain:\n')
2592	pFile.write(ptlbls+'\n')
2593	pFile.write(ptstr+'\n')
2594	pFile.write(sigstr+'\n')
2595
2596	def PrintSHPOAndSig(pfx,hapData,POsig):
2597	pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2598	ptlbls = ' names :'
2599	ptstr = ' values:'
2600	sigstr = ' sig :'
2601	for item in hapData[5]:
2602	ptlbls += '%12s'%(item)
2603	ptstr += '%12.3f'%(hapData[5][item])
2604	if pfx+item in POsig:
2605	sigstr += '%12.3f'%(POsig[pfx+item])
2606	else:
2607	sigstr += 12*' '
2608	pFile.write(ptlbls+'\n')
2609	pFile.write(ptstr+'\n')
2610	pFile.write(sigstr+'\n')
2611
2612	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2613	pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2614	text = ''
2615	for item in hapData[2]:
2616	if pfx+item in ScalExtSig:
2617	text += ' %s: '%(item)
2618	text += '%12.2e'%(hapData[2][item][0])
2619	if pfx+item in ScalExtSig:
2620	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2621	pFile.write(text+'\n')
2622
2623	def PrintBabinetAndSig(pfx,hapData,BabSig):
2624	pFile.write('\n Babinet form factor modification:\n')
2625	ptlbls = ' names :'
2626	ptstr = ' values:'
2627	sigstr = ' sig :'
2628	for item in hapData:
2629	ptlbls += '%12s'%(item)
2630	ptstr += '%12.3f'%(hapData[item][0])
2631	if pfx+item in BabSig:
2632	sigstr += '%12.3f'%(BabSig[pfx+item])
2633	else:
2634	sigstr += 12*' '
2635	pFile.write(ptlbls+'\n')
2636	pFile.write(ptstr+'\n')
2637	pFile.write(sigstr+'\n')
2638
2639	def PrintTwinsAndSig(pfx,twinData,TwinSig):
2640	pFile.write('\n Twin Law fractions :\n')
2641	ptlbls = ' names :'
2642	ptstr = ' values:'
2643	sigstr = ' sig :'
2644	for it,item in enumerate(twinData):
2645	ptlbls += ' twin #%d'%(it)
2646	if it:
2647	ptstr += '%12.3f'%(item[1])
2648	else:
2649	ptstr += '%12.3f'%(item[1][0])
2650	if pfx+'TwinFr:'+str(it) in TwinSig:
2651	sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2652	else:
2653	sigstr += 12*' '
2654	pFile.write(ptlbls+'\n')
2655	pFile.write(ptstr+'\n')
2656	pFile.write(sigstr+'\n')
2657
2658
2659	PhFrExtPOSig = {}
2660	SizeMuStrSig = {}
2661	ScalExtSig = {}
2662	BabSig = {}
2663	TwinFrSig = {}
2664	wtFrSum = {}
2665	for phase in Phases:
2666	HistoPhase = Phases[phase]['Histograms']
2667	General = Phases[phase]['General']
2668	SGData = General['SGData']
2669	pId = Phases[phase]['pId']
2670	histoList = list(HistoPhase.keys())
2671	histoList.sort()
2672	for histogram in histoList:
2673	try:
2674	Histogram = Histograms[histogram]
2675	except KeyError:
2676	#skip if histogram not included e.g. in a sequential refinement
2677	continue
2678	if not Phases[phase]['Histograms'][histogram]['Use']:
2679	#skip if phase absent from this histogram
2680	continue
2681	hapData = HistoPhase[histogram]
2682	hId = Histogram['hId']
2683	pfx = str(pId)+':'+str(hId)+':'
2684	if hId not in wtFrSum:
2685	wtFrSum[hId] = 0.
2686	if 'PWDR' in histogram:
2687	for item in ['Scale','Extinction']:
2688	hapData[item][0] = parmDict[pfx+item]
2689	if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2690	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2691	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2692	if hapData['Pref.Ori.'][0] == 'MD':
2693	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2694	if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2695	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2696	else: #'SH' spherical harmonics
2697	for item in hapData['Pref.Ori.'][5]:
2698	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2699	if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2700	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2701	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2702	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2703	pfx+'HStrain':{}})
2704	for item in ['Mustrain','Size']:
2705	hapData[item][1][2] = parmDict[pfx+item+';mx']
2706	# hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2707	if pfx+item+';mx' in sigDict:
2708	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2709	if hapData[item][0] in ['isotropic','uniaxial']:
2710	hapData[item][1][0] = parmDict[pfx+item+';i']
2711	if item == 'Size':
2712	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2713	if pfx+item+';i' in sigDict:
2714	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2715	if hapData[item][0] == 'uniaxial':
2716	hapData[item][1][1] = parmDict[pfx+item+';a']
2717	if item == 'Size':
2718	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2719	if pfx+item+';a' in sigDict:
2720	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2721	else: #generalized for mustrain or ellipsoidal for size
2722	Nterms = len(hapData[item][4])
2723	for i in range(Nterms):
2724	sfx = ';'+str(i)
2725	hapData[item][4][i] = parmDict[pfx+item+sfx]
2726	if pfx+item+sfx in sigDict:
2727	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2728	names = G2spc.HStrainNames(SGData)
2729	for i,name in enumerate(names):
2730	hapData['HStrain'][0][i] = parmDict[pfx+name]
2731	if pfx+name in sigDict:
2732	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2733	if Phases[phase]['General']['Type'] != 'magnetic':
2734	for name in ['BabA','BabU']:
2735	hapData['Babinet'][name][0] = parmDict[pfx+name]
2736	if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
2737	BabSig[pfx+name] = sigDict[pfx+name]
2738
2739	elif 'HKLF' in histogram:
2740	for item in ['Scale','Flack']:
2741	if parmDict.get(pfx+item):
2742	hapData[item][0] = parmDict[pfx+item]
2743	if pfx+item in sigDict:
2744	ScalExtSig[pfx+item] = sigDict[pfx+item]
2745	for item in ['Ep','Eg','Es']:
2746	if parmDict.get(pfx+item):
2747	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2748	if pfx+item in sigDict:
2749	ScalExtSig[pfx+item] = sigDict[pfx+item]
2750	for item in ['BabA','BabU']:
2751	hapData['Babinet'][item][0] = parmDict[pfx+item]
2752	if pfx+item in sigDict:
2753	BabSig[pfx+item] = sigDict[pfx+item]
2754	if 'Twins' in hapData:
2755	it = 1
2756	sumTwFr = 0.
2757	while True:
2758	try:
2759	hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2760	if pfx+'TwinFr:'+str(it) in sigDict:
2761	TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2762	if it:
2763	sumTwFr += hapData['Twins'][it][1]
2764	it += 1
2765	except KeyError:
2766	break
2767	hapData['Twins'][0][1][0] = 1.-sumTwFr
2768
2769	if Print:
2770	for phase in Phases:
2771	HistoPhase = Phases[phase]['Histograms']
2772	General = Phases[phase]['General']
2773	SGData = General['SGData']
2774	pId = Phases[phase]['pId']
2775	histoList = list(HistoPhase.keys())
2776	histoList.sort()
2777	for histogram in histoList:
2778	try:
2779	Histogram = Histograms[histogram]
2780	except KeyError:
2781	#skip if histogram not included e.g. in a sequential refinement
2782	continue
2783	hapData = HistoPhase[histogram]
2784	hId = Histogram['hId']
2785	Histogram['Residuals'][str(pId)+'::Name'] = phase
2786	pfx = str(pId)+':'+str(hId)+':'
2787	hfx = ':%s:'%(hId)
2788	if pfx+'Nref' not in Histogram['Residuals']: #skip not used phase in histogram
2789	continue
2790	pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2791	pFile.write(135*'='+'\n')
2792	if 'PWDR' in histogram:
2793	pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
2794	(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
2795	pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
2796	pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
2797
2798	if parmDict[pfx+'LeBail']:
2799	pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
2800	else:
2801	if pfx+'Scale' in PhFrExtPOSig:
2802	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2803	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2804	pFile.write(' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f\n'%
2805	(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
2806	if pfx+'Extinction' in PhFrExtPOSig:
2807	pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
2808	if hapData['Pref.Ori.'][0] == 'MD':
2809	if pfx+'MD' in PhFrExtPOSig:
2810	pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
2811	else:
2812	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2813	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2814	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2815	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2816	if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
2817	if len(BabSig):
2818	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2819
2820	elif 'HKLF' in histogram:
2821	Inst = Histogram['Instrument Parameters'][0]
2822	pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
2823	(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2824	Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
2825	if FFtables != None and 'N' not in Inst['Type'][0]:
2826	PrintFprime(FFtables,hfx,pFile)
2827	pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
2828	if pfx+'Scale' in ScalExtSig:
2829	pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
2830	if hapData['Extinction'][0] != 'None':
2831	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2832	if len(BabSig):
2833	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2834	if pfx+'Flack' in ScalExtSig:
2835	pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
2836	if len(TwinFrSig):
2837	PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2838
2839	################################################################################
2840	##### Histogram data
2841	################################################################################
2842
2843	def GetHistogramData(Histograms,Print=True,pFile=None):
2844	'needs a doc string'
2845
2846	def GetBackgroundParms(hId,Background):
2847	Back = Background[0]
2848	DebyePeaks = Background[1]
2849	bakType,bakFlag = Back[:2]
2850	backVals = Back[3:]
2851	backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2852	backDict = dict(zip(backNames,backVals))
2853	backVary = []
2854	if bakFlag:
2855	backVary = backNames
2856	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2857	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2858	debyeDict = {}
2859	debyeList = []
2860	for i in range(DebyePeaks['nDebye']):
2861	debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2862	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2863	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2864	debyeVary = []
2865	for item in debyeList:
2866	if item[1]:
2867	debyeVary.append(item[0])
2868	backDict.update(debyeDict)
2869	backVary += debyeVary
2870	peakDict = {}
2871	peakList = []
2872	for i in range(DebyePeaks['nPeaks']):
2873	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2874	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2875	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2876	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2877	peakVary = []
2878	for item in peakList:
2879	if item[1]:
2880	peakVary.append(item[0])
2881	backDict.update(peakDict)
2882	backVary += peakVary
2883	return bakType,backDict,backVary
2884
2885	def GetInstParms(hId,Inst):
2886	#patch
2887	if 'Z' not in Inst:
2888	Inst['Z'] = [0.0,0.0,False]
2889	dataType = Inst['Type'][0]
2890	instDict = {}
2891	insVary = []
2892	pfx = ':'+str(hId)+':'
2893	insKeys = list(Inst.keys())
2894	insKeys.sort()
2895	for item in insKeys:
2896	insName = pfx+item
2897	instDict[insName] = Inst[item][1]
2898	if len(Inst[item]) > 2 and Inst[item][2]:
2899	insVary.append(insName)
2900	if 'C' in dataType:
2901	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2902	elif 'T' in dataType: #trap zero alp, bet coeff.
2903	if not instDict[pfx+'alpha']:
2904	instDict[pfx+'alpha'] = 1.0
2905	if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
2906	instDict[pfx+'beta-1'] = 1.0
2907	return dataType,instDict,insVary
2908
2909	def GetSampleParms(hId,Sample):
2910	sampVary = []
2911	hfx = ':'+str(hId)+':'
2912	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2913	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2914	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2915	if key in Sample:
2916	sampDict[hfx+key] = Sample[key]
2917	Type = Sample['Type']
2918	if 'Bragg' in Type: #Bragg-Brentano
2919	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2920	sampDict[hfx+item] = Sample[item][0]
2921	if Sample[item][1]:
2922	sampVary.append(hfx+item)
2923	elif 'Debye' in Type: #Debye-Scherrer
2924	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2925	sampDict[hfx+item] = Sample[item][0]
2926	if Sample[item][1]:
2927	sampVary.append(hfx+item)
2928	return Type,sampDict,sampVary
2929
2930	def PrintBackground(Background):
2931	Back = Background[0]
2932	DebyePeaks = Background[1]
2933	pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
2934	line = ' Coefficients: '
2935	for i,back in enumerate(Back[3:]):
2936	line += '%10.3f'%(back)
2937	if i and not i%10:
2938	line += '\n'+15*' '
2939	pFile.write(line+'\n')
2940	if DebyePeaks['nDebye']:
2941	pFile.write('\n Debye diffuse scattering coefficients\n')
2942	parms = ['DebyeA','DebyeR','DebyeU']
2943	line = ' names : '
2944	for parm in parms:
2945	line += '%8s refine?'%(parm)
2946	pFile.write(line+'\n')
2947	for j,term in enumerate(DebyePeaks['debyeTerms']):
2948	line = ' term'+'%2d'%(j)+':'
2949	for i in range(3):
2950	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2951	pFile.write(line+'\n')
2952	if DebyePeaks['nPeaks']:
2953	pFile.write('\n Single peak coefficients\n')
2954	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2955	line = ' names : '
2956	for parm in parms:
2957	line += '%8s refine?'%(parm)
2958	pFile.write(line+'\n')
2959	for j,term in enumerate(DebyePeaks['peaksList']):
2960	line = ' peak'+'%2d'%(j)+':'
2961	for i in range(4):
2962	line += '%12.3f %5s'%(term[2i],bool(term[2i+1]))
2963	pFile.write(line+'\n')
2964
2965	def PrintInstParms(Inst):
2966	pFile.write('\n Instrument Parameters:\n')
2967	insKeys = list(Inst.keys())
2968	insKeys.sort()
2969	iBeg = 0
2970	Ok = True
2971	while Ok:
2972	ptlbls = ' name :'
2973	ptstr = ' value :'
2974	varstr = ' refine:'
2975	iFin = min(iBeg+9,len(insKeys))
2976	for item in insKeys[iBeg:iFin]:
2977	if item not in ['Type','Source']:
2978	ptlbls += '%12s' % (item)
2979	ptstr += '%12.6f' % (Inst[item][1])
2980	if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2981	varstr += 12*' '
2982	else:
2983	varstr += '%12s' % (str(bool(Inst[item][2])))
2984	pFile.write(ptlbls+'\n')
2985	pFile.write(ptstr+'\n')
2986	pFile.write(varstr+'\n')
2987	iBeg = iFin
2988	if iBeg == len(insKeys):
2989	Ok = False
2990	else:
2991	pFile.write('\n')
2992
2993	def PrintSampleParms(Sample):
2994	pFile.write('\n Sample Parameters:\n')
2995	pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
2996	(Sample['Omega'],Sample['Chi'],Sample['Phi']))
2997	ptlbls = ' name :'
2998	ptstr = ' value :'
2999	varstr = ' refine:'
3000	if 'Bragg' in Sample['Type']:
3001	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3002	ptlbls += '%14s'%(item)
3003	ptstr += '%14.4f'%(Sample[item][0])
3004	varstr += '%14s'%(str(bool(Sample[item][1])))
3005
3006	elif 'Debye' in Type: #Debye-Scherrer
3007	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3008	ptlbls += '%14s'%(item)
3009	ptstr += '%14.4f'%(Sample[item][0])
3010	varstr += '%14s'%(str(bool(Sample[item][1])))
3011
3012	pFile.write(ptlbls+'\n')
3013	pFile.write(ptstr+'\n')
3014	pFile.write(varstr+'\n')
3015
3016	histDict = {}
3017	histVary = []
3018	controlDict = {}
3019	histoList = list(Histograms.keys())
3020	histoList.sort()
3021	for histogram in histoList:
3022	if 'PWDR' in histogram:
3023	Histogram = Histograms[histogram]
3024	hId = Histogram['hId']
3025	pfx = ':'+str(hId)+':'
3026	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3027	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3028	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3029	for excl in controlDict[pfx+'Exclude']:
3030	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3031	if controlDict[pfx+'Exclude']:
3032	ma.mask_rows(Histogram['Data'])
3033	Background = Histogram['Background']
3034	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3035	controlDict[pfx+'bakType'] = Type
3036	histDict.update(bakDict)
3037	histVary += bakVary
3038
3039	Inst = Histogram['Instrument Parameters'] #TODO ? ignores tabulated alp,bet & delt for TOF
3040	if 'T' in Type and len(Inst[1]): #patch - back-to-back exponential contribution to TOF line shape is removed
3041	print ('Warning: tabulated profile coefficients are ignored')
3042	Type,instDict,insVary = GetInstParms(hId,Inst[0])
3043	controlDict[pfx+'histType'] = Type
3044	if 'XC' in Type:
3045	if pfx+'Lam1' in instDict:
3046	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3047	else:
3048	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
3049	histDict.update(instDict)
3050	histVary += insVary
3051
3052	Sample = Histogram['Sample Parameters']
3053	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3054	controlDict[pfx+'instType'] = Type
3055	histDict.update(sampDict)
3056	histVary += sampVary
3057
3058
3059	if Print:
3060	pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3061	pFile.write(135*'='+'\n')
3062	Units = {'C':' deg','T':' msec'}
3063	units = Units[controlDict[pfx+'histType'][2]]
3064	Limits = controlDict[pfx+'Limits']
3065	pFile.write(' Instrument type: %s\n'%Sample['Type'])
3066	pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3067	if len(controlDict[pfx+'Exclude']):
3068	excls = controlDict[pfx+'Exclude']
3069	for excl in excls:
3070	pFile.write(' Excluded region: %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units))
3071	PrintSampleParms(Sample)
3072	PrintInstParms(Inst[0])
3073	PrintBackground(Background)
3074	elif 'HKLF' in histogram:
3075	Histogram = Histograms[histogram]
3076	hId = Histogram['hId']
3077	pfx = ':'+str(hId)+':'
3078	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3079	Inst = Histogram['Instrument Parameters'][0]
3080	controlDict[pfx+'histType'] = Inst['Type'][0]
3081	if 'X' in Inst['Type'][0]:
3082	histDict[pfx+'Lam'] = Inst['Lam'][1]
3083	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3084	return histVary,histDict,controlDict
3085
3086	def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3087	'needs a doc string'
3088
3089	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3090	Back = Background[0]
3091	DebyePeaks = Background[1]
3092	lenBack = len(Back[3:])
3093	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
3094	for i in range(lenBack):
3095	Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3096	if pfx+'Back;'+str(i) in sigDict:
3097	backSig[i] = sigDict[pfx+'Back;'+str(i)]
3098	if DebyePeaks['nDebye']:
3099	for i in range(DebyePeaks['nDebye']):
3100	names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3101	for j,name in enumerate(names):
3102	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3103	if name in sigDict:
3104	backSig[lenBack+3*i+j] = sigDict[name]
3105	if DebyePeaks['nPeaks']:
3106	for i in range(DebyePeaks['nPeaks']):
3107	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3108	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3109	for j,name in enumerate(names):
3110	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3111	if name in sigDict:
3112	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
3113	return backSig
3114
3115	def SetInstParms(pfx,Inst,parmDict,sigDict):
3116	instSig = {}
3117	insKeys = list(Inst.keys())
3118	insKeys.sort()
3119	for item in insKeys:
3120	insName = pfx+item
3121	Inst[item][1] = parmDict[insName]
3122	if insName in sigDict:
3123	instSig[item] = sigDict[insName]
3124	else:
3125	instSig[item] = 0
3126	return instSig
3127
3128	def SetSampleParms(pfx,Sample,parmDict,sigDict):
3129	if 'Bragg' in Sample['Type']: #Bragg-Brentano
3130	sampSig = [0 for i in range(5)]
3131	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3132	Sample[item][0] = parmDict[pfx+item]
3133	if pfx+item in sigDict:
3134	sampSig[i] = sigDict[pfx+item]
3135	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3136	sampSig = [0 for i in range(4)]
3137	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3138	Sample[item][0] = parmDict[pfx+item]
3139	if pfx+item in sigDict:
3140	sampSig[i] = sigDict[pfx+item]
3141	return sampSig
3142
3143	def PrintBackgroundSig(Background,backSig):
3144	Back = Background[0]
3145	DebyePeaks = Background[1]
3146	valstr = ' value : '
3147	sigstr = ' sig : '
3148	refine = False
3149	for i,back in enumerate(Back[3:]):
3150	valstr += '%10.4g'%(back)
3151	if Back[1]:
3152	refine = True
3153	sigstr += '%10.4g'%(backSig[i])
3154	else:
3155	sigstr += 10*' '
3156	if refine:
3157	pFile.write('\n Background function: %s\n'%Back[0])
3158	pFile.write(valstr+'\n')
3159	pFile.write(sigstr+'\n')
3160	if DebyePeaks['nDebye']:
3161	ifAny = False
3162	ptfmt = "%12.3f"
3163	names = ' names :'
3164	ptstr = ' values:'
3165	sigstr = ' esds :'
3166	for item in sigDict:
3167	if 'Debye' in item:
3168	ifAny = True
3169	names += '%12s'%(item)
3170	ptstr += ptfmt%(parmDict[item])
3171	sigstr += ptfmt%(sigDict[item])
3172	if ifAny:
3173	pFile.write('\n Debye diffuse scattering coefficients\n')
3174	pFile.write(names+'\n')
3175	pFile.write(ptstr+'\n')
3176	pFile.write(sigstr+'\n')
3177	if DebyePeaks['nPeaks']:
3178	pFile.write('\n Single peak coefficients:\n')
3179	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3180	line = ' peak no. '
3181	for parm in parms:
3182	line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3183	pFile.write(line+'\n')
3184	for ip in range(DebyePeaks['nPeaks']):
3185	ptstr = ' %4d '%(ip)
3186	for parm in parms:
3187	fmt = '%14.3f'
3188	efmt = '%12.3f'
3189	if parm == 'BkPkpos':
3190	fmt = '%14.4f'
3191	efmt = '%12.4f'
3192	name = pfx+parm+';%d'%(ip)
3193	ptstr += fmt%(parmDict[name])
3194	if name in sigDict:
3195	ptstr += efmt%(sigDict[name])
3196	else:
3197	ptstr += 12*' '
3198	pFile.write(ptstr+'\n')
3199	sumBk = np.array(Histogram['sumBk'])
3200	pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3201	(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3202
3203	def PrintInstParmsSig(Inst,instSig):
3204	refine = False
3205	insKeys = list(instSig.keys())
3206	insKeys.sort()
3207	iBeg = 0
3208	Ok = True
3209	while Ok:
3210	ptlbls = ' names :'
3211	ptstr = ' value :'
3212	sigstr = ' sig :'
3213	iFin = min(iBeg+9,len(insKeys))
3214	for name in insKeys[iBeg:iFin]:
3215	if name not in ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3216	ptlbls += '%12s' % (name)
3217	ptstr += '%12.6f' % (Inst[name][1])
3218	if instSig[name]:
3219	refine = True
3220	sigstr += '%12.6f' % (instSig[name])
3221	else:
3222	sigstr += 12*' '
3223	if refine:
3224	pFile.write('\n Instrument Parameters:\n')
3225	pFile.write(ptlbls+'\n')
3226	pFile.write(ptstr+'\n')
3227	pFile.write(sigstr+'\n')
3228	iBeg = iFin
3229	if iBeg == len(insKeys):
3230	Ok = False
3231
3232	def PrintSampleParmsSig(Sample,sampleSig):
3233	ptlbls = ' names :'
3234	ptstr = ' values:'
3235	sigstr = ' sig :'
3236	refine = False
3237	if 'Bragg' in Sample['Type']:
3238	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3239	ptlbls += '%14s'%(item)
3240	ptstr += '%14.4f'%(Sample[item][0])
3241	if sampleSig[i]:
3242	refine = True
3243	sigstr += '%14.4f'%(sampleSig[i])
3244	else:
3245	sigstr += 14*' '
3246
3247	elif 'Debye' in Sample['Type']: #Debye-Scherrer
3248	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3249	ptlbls += '%14s'%(item)
3250	ptstr += '%14.4f'%(Sample[item][0])
3251	if sampleSig[i]:
3252	refine = True
3253	sigstr += '%14.4f'%(sampleSig[i])
3254	else:
3255	sigstr += 14*' '
3256
3257	if refine:
3258	pFile.write('\n Sample Parameters:\n')
3259	pFile.write(ptlbls+'\n')
3260	pFile.write(ptstr+'\n')
3261	pFile.write(sigstr+'\n')
3262
3263	histoList = list(Histograms.keys())
3264	histoList.sort()
3265	for histogram in histoList:
3266	if 'PWDR' in histogram:
3267	Histogram = Histograms[histogram]
3268	hId = Histogram['hId']
3269	pfx = ':'+str(hId)+':'
3270	Background = Histogram['Background']
3271	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3272
3273	Inst = Histogram['Instrument Parameters'][0]
3274	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3275
3276	Sample = Histogram['Sample Parameters']
3277	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3278
3279	pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3280	pFile.write(135*'='+'\n')
3281	pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3282	pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3283	(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3284	pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3285	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3286	if Print:
3287	pFile.write(' Instrument type: %s\n'%Sample['Type'])
3288	if FFtables != None and 'N' not in Inst['Type'][0]:
3289	PrintFprime(FFtables,pfx,pFile)
3290	PrintSampleParmsSig(Sample,sampSig)
3291	PrintInstParmsSig(Inst,instSig)
3292	PrintBackgroundSig(Background,backSig)
3293

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

Download in other formats: