Context Navigation

source: trunk/GSASIIstrIO.py @ 2748

Last change on this file since 2748 was 2748, checked in by vondreele, 7 years ago
Sigh, another LeBail? fix
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 148.0 KB

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

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

Download in other formats: