Context Navigation

source: trunk/GSASIIstrIO.py @ 1357

Last change on this file since 1357 was 1312, checked in by toby, 11 years ago
Parametric fit; work on PseudoVars? (more 2 come); minor code cleanups
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 111.2 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2014-05-01 02:42:51 +0000 (Thu, 01 May 2014) $
4	# $Author: toby $
5	# $Revision: 1312 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 1312 2014-05-01 02:42:51Z toby $
8	########### SVN repository information ###################
9	'''
10	GSASIIstrIO: structure I/O routines
11	-------------------------------------
12
13	'''
14	import sys
15	import os
16	import os.path as ospath
17	import time
18	import math
19	import copy
20	import random
21	import cPickle
22	import numpy as np
23	import numpy.ma as ma
24	import numpy.linalg as nl
25	import scipy.optimize as so
26	import GSASIIpath
27	GSASIIpath.SetVersionNumber("$Revision: 1312 $")
28	import GSASIIElem as G2el
29	import GSASIIgrid as G2gd
30	import GSASIIlattice as G2lat
31	import GSASIIspc as G2spc
32	import GSASIIobj as G2obj
33	import GSASIImapvars as G2mv
34	import GSASIImath as G2mth
35
36	sind = lambda x: np.sin(x*np.pi/180.)
37	cosd = lambda x: np.cos(x*np.pi/180.)
38	tand = lambda x: np.tan(x*np.pi/180.)
39	asind = lambda x: 180.*np.arcsin(x)/np.pi
40	acosd = lambda x: 180.*np.arccos(x)/np.pi
41	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
42
43	ateln2 = 8.0*math.log(2.0)
44
45	def GetControls(GPXfile):
46	''' Returns dictionary of control items found in GSASII gpx file
47
48	:param str GPXfile: full .gpx file name
49	:return: dictionary of control items
50	'''
51	Controls = copy.copy(G2obj.DefaultControls)
52	fl = open(GPXfile,'rb')
53	while True:
54	try:
55	data = cPickle.load(fl)
56	except EOFError:
57	break
58	datum = data[0]
59	if datum[0] == 'Controls':
60	Controls.update(datum[1])
61	fl.close()
62	return Controls
63
64	def GetConstraints(GPXfile):
65	'''Read the constraints from the GPX file and interpret them
66
67	called in :func:`CheckConstraints`, :func:`GSASIIstrMain.Refine`
68	and :func:`GSASIIstrMain.SeqRefine`.
69	'''
70	fl = open(GPXfile,'rb')
71	while True:
72	try:
73	data = cPickle.load(fl)
74	except EOFError:
75	break
76	datum = data[0]
77	if datum[0] == 'Constraints':
78	constList = []
79	for item in datum[1]:
80	if item.startswith('_'): continue
81	constList += datum[1][item]
82	fl.close()
83	constDict,fixedList,ignored = ProcessConstraints(constList)
84	if ignored:
85	print ignored,'Constraints were rejected. Was a constrained phase, histogram or atom deleted?'
86	return constDict,fixedList
87	fl.close()
88	raise Exception,"No constraints in GPX file"
89
90	def ProcessConstraints(constList):
91	"""Interpret the constraints in the constList input into a dictionary, etc.
92	All :class:`GSASIIobj.G2VarObj` objects are mapped to the appropriate
93	phase/hist/atoms based on the object internals (random Ids). If this can't be
94	done (if a phase has been deleted, etc.), the variable is ignored.
95	If the constraint cannot be used due to too many dropped variables,
96	it is counted as ignored.
97
98	:param list constList: a list of lists where each item in the outer list
99	specifies a constraint of some form, as described in the :mod:`GSASIIobj`
100	:ref:`Constraint definition <Constraint_definitions_table>`.
101
102	:returns: a tuple of (constDict,fixedList,ignored) where:
103
104	* constDict (list of dicts) contains the constraint relationships
105	* fixedList (list) contains the fixed values for each type
106	of constraint.
107	* ignored (int) counts the number of invalid constraint items
108	(should always be zero!)
109	"""
110	constDict = []
111	fixedList = []
112	ignored = 0
113	for item in constList:
114	if item[-1] == 'h':
115	# process a hold
116	fixedList.append('0')
117	var = str(item[0][1])
118	if '?' not in var:
119	constDict.append({var:0.0})
120	else:
121	ignored += 1
122	elif item[-1] == 'f':
123	# process a new variable
124	fixedList.append(None)
125	D = {}
126	varyFlag = item[-2]
127	varname = item[-3]
128	for term in item[:-3]:
129	var = str(term[1])
130	if '?' not in var:
131	D[var] = term[0]
132	if len(D) > 1:
133	# add extra dict terms for input variable name and vary flag
134	if varname is not None:
135	if varname.startswith('::'):
136	varname = varname[2:].replace(':',';')
137	else:
138	varname = varname.replace(':',';')
139	D['_name'] = '::' + varname
140	D['_vary'] = varyFlag == True # force to bool
141	constDict.append(D)
142	else:
143	ignored += 1
144	#constFlag[-1] = ['Vary']
145	elif item[-1] == 'c':
146	# process a contraint relationship
147	D = {}
148	for term in item[:-3]:
149	var = str(term[1])
150	if '?' not in var:
151	D[var] = term[0]
152	if len(D) >= 1:
153	fixedList.append(str(item[-3]))
154	constDict.append(D)
155	else:
156	ignored += 1
157	elif item[-1] == 'e':
158	# process an equivalence
159	firstmult = None
160	eqlist = []
161	for term in item[:-3]:
162	if term[0] == 0: term[0] = 1.0
163	var = str(term[1])
164	if '?' in var: continue
165	if firstmult is None:
166	firstmult = term[0]
167	firstvar = var
168	else:
169	eqlist.append([var,firstmult/term[0]])
170	if len(eqlist) > 0:
171	G2mv.StoreEquivalence(firstvar,eqlist)
172	else:
173	ignored += 1
174	else:
175	ignored += 1
176	return constDict,fixedList,ignored
177
178	def CheckConstraints(GPXfile):
179	'''Load constraints and related info and return any error or warning messages'''
180	# init constraints
181	G2mv.InitVars()
182	# get variables
183	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
184	if not Phases:
185	return 'Error: No Phases!',''
186	if not Histograms:
187	return 'Error: no diffraction data',''
188	rigidbodyDict = GetRigidBodies(GPXfile)
189	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
190	rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
191	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)
192	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
193	histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
194	varyList = rbVary+phaseVary+hapVary+histVary
195	constrDict,fixedList = GetConstraints(GPXfile)
196	return G2mv.CheckConstraints(varyList,constrDict,fixedList)
197
198	def GetRestraints(GPXfile):
199	'''Read the restraints from the GPX file.
200	Throws an exception if not found in the .GPX file
201	'''
202	fl = open(GPXfile,'rb')
203	while True:
204	try:
205	data = cPickle.load(fl)
206	except EOFError:
207	break
208	datum = data[0]
209	if datum[0] == 'Restraints':
210	restraintDict = datum[1]
211	fl.close()
212	return restraintDict
213
214	def GetRigidBodies(GPXfile):
215	'''Read the rigid body models from the GPX file
216	'''
217	fl = open(GPXfile,'rb')
218	while True:
219	try:
220	data = cPickle.load(fl)
221	except EOFError:
222	break
223	datum = data[0]
224	if datum[0] == 'Rigid bodies':
225	rigidbodyDict = datum[1]
226	fl.close()
227	return rigidbodyDict
228
229	def GetFprime(controlDict,Histograms):
230	'Needs a doc string'
231	FFtables = controlDict['FFtables']
232	if not FFtables:
233	return
234	histoList = Histograms.keys()
235	histoList.sort()
236	for histogram in histoList:
237	if histogram[:4] in ['PWDR','HKLF']:
238	Histogram = Histograms[histogram]
239	hId = Histogram['hId']
240	hfx = ':%d:'%(hId)
241	keV = controlDict[hfx+'keV']
242	for El in FFtables:
243	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
244	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
245	FFtables[El][hfx+'FP'] = FP
246	FFtables[El][hfx+'FPP'] = FPP
247
248	def GetPhaseNames(GPXfile):
249	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
250
251	:param str GPXfile: full .gpx file name
252	:return: list of phase names
253	'''
254	fl = open(GPXfile,'rb')
255	PhaseNames = []
256	while True:
257	try:
258	data = cPickle.load(fl)
259	except EOFError:
260	break
261	datum = data[0]
262	if 'Phases' == datum[0]:
263	for datus in data[1:]:
264	PhaseNames.append(datus[0])
265	fl.close()
266	return PhaseNames
267
268	def GetAllPhaseData(GPXfile,PhaseName):
269	''' Returns the entire dictionary for PhaseName from GSASII gpx file
270
271	:param str GPXfile: full .gpx file name
272	:param str PhaseName: phase name
273	:return: phase dictionary
274	'''
275	fl = open(GPXfile,'rb')
276	General = {}
277	Atoms = []
278	while True:
279	try:
280	data = cPickle.load(fl)
281	except EOFError:
282	break
283	datum = data[0]
284	if 'Phases' == datum[0]:
285	for datus in data[1:]:
286	if datus[0] == PhaseName:
287	break
288	fl.close()
289	return datus[1]
290
291	def GetHistograms(GPXfile,hNames):
292	""" Returns a dictionary of histograms found in GSASII gpx file
293
294	:param str GPXfile: full .gpx file name
295	:param str hNames: list of histogram names
296	:return: dictionary of histograms (types = PWDR & HKLF)
297
298	"""
299	fl = open(GPXfile,'rb')
300	Histograms = {}
301	while True:
302	try:
303	data = cPickle.load(fl)
304	except EOFError:
305	break
306	datum = data[0]
307	hist = datum[0]
308	if hist in hNames:
309	if 'PWDR' in hist[:4]:
310	PWDRdata = {}
311	PWDRdata.update(datum[1][0]) #weight factor
312	PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1])) #masked powder data arrays/clear previous masks
313	PWDRdata[data[2][0]] = data[2][1] #Limits & excluded regions (if any)
314	PWDRdata[data[3][0]] = data[3][1] #Background
315	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
316	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
317	try:
318	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
319	except IndexError:
320	PWDRdata['Reflection Lists'] = {}
321	PWDRdata['Residuals'] = {}
322
323	Histograms[hist] = PWDRdata
324	elif 'HKLF' in hist[:4]:
325	HKLFdata = {}
326	HKLFdata.update(datum[1][0]) #weight factor
327	#patch
328	if isinstance(datum[1][1],list):
329	RefData = {'RefList':[],'FF':{}}
330	for ref in datum[1][1]:
331	RefData['RefList'].append(ref[:11]+[ref[13],])
332	RefData['RefList'] = np.array(RefData['RefList'])
333	datum[1][1] = RefData
334	#end patch
335	datum[1][1]['FF'] = {}
336	HKLFdata['Data'] = datum[1][1]
337	HKLFdata[data[1][0]] = data[1][1] #Instrument parameters
338	HKLFdata['Reflection Lists'] = None
339	HKLFdata['Residuals'] = {}
340	Histograms[hist] = HKLFdata
341	fl.close()
342	return Histograms
343
344	def GetHistogramNames(GPXfile,hType):
345	""" Returns a list of histogram names found in GSASII gpx file
346
347	:param str GPXfile: full .gpx file name
348	:param str hType: list of histogram types
349	:return: list of histogram names (types = PWDR & HKLF)
350
351	"""
352	fl = open(GPXfile,'rb')
353	HistogramNames = []
354	while True:
355	try:
356	data = cPickle.load(fl)
357	except EOFError:
358	break
359	datum = data[0]
360	if datum[0][:4] in hType:
361	HistogramNames.append(datum[0])
362	fl.close()
363	return HistogramNames
364
365	def GetUsedHistogramsAndPhases(GPXfile):
366	''' Returns all histograms that are found in any phase
367	and any phase that uses a histogram. This also
368	assigns numbers to used phases and histograms by the
369	order they appear in the file.
370
371	:param str GPXfile: full .gpx file name
372	:returns: (Histograms,Phases)
373
374	* Histograms = dictionary of histograms as {name:data,...}
375	* Phases = dictionary of phases that use histograms
376
377	'''
378	phaseNames = GetPhaseNames(GPXfile)
379	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
380	allHistograms = GetHistograms(GPXfile,histoList)
381	phaseData = {}
382	for name in phaseNames:
383	phaseData[name] = GetAllPhaseData(GPXfile,name)
384	Histograms = {}
385	Phases = {}
386	for phase in phaseData:
387	Phase = phaseData[phase]
388	if Phase['Histograms']:
389	if phase not in Phases:
390	pId = phaseNames.index(phase)
391	Phase['pId'] = pId
392	Phases[phase] = Phase
393	for hist in Phase['Histograms']:
394	if 'Use' not in Phase['Histograms'][hist]: #patch
395	Phase['Histograms'][hist]['Use'] = True
396	if hist not in Histograms and Phase['Histograms'][hist]['Use']:
397	try:
398	Histograms[hist] = allHistograms[hist]
399	hId = histoList.index(hist)
400	Histograms[hist]['hId'] = hId
401	except KeyError: # would happen if a referenced histogram were
402	# renamed or deleted
403	print('For phase "'+str(phase)+
404	'" unresolved reference to histogram "'+str(hist)+'"')
405	G2obj.IndexAllIds(Histograms=Histograms,Phases=Phases)
406	return Histograms,Phases
407
408	def getBackupName(GPXfile,makeBack):
409	'''
410	Get the name for the backup .gpx file name
411
412	:param str GPXfile: full .gpx file name
413	:param bool makeBack: if True the name of a new file is returned, if
414	False the name of the last file that exists is returned
415	:returns: the name of a backup file
416
417	'''
418	GPXpath,GPXname = ospath.split(GPXfile)
419	if GPXpath == '': GPXpath = '.'
420	Name = ospath.splitext(GPXname)[0]
421	files = os.listdir(GPXpath)
422	last = 0
423	for name in files:
424	name = name.split('.')
425	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
426	if makeBack:
427	last = max(last,int(name[1].strip('bak'))+1)
428	else:
429	last = max(last,int(name[1].strip('bak')))
430	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
431	return GPXback
432
433	def GPXBackup(GPXfile,makeBack=True):
434	'''
435	makes a backup of the current .gpx file (?)
436
437	:param str GPXfile: full .gpx file name
438	:param bool makeBack: if True (default), the backup is written to
439	a new file; if False, the last backup is overwritten
440	:returns: the name of the backup file that was written
441	'''
442	import distutils.file_util as dfu
443	GPXback = getBackupName(GPXfile,makeBack)
444	dfu.copy_file(GPXfile,GPXback)
445	return GPXback
446
447	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
448	''' Updates gpxfile from all histograms that are found in any phase
449	and any phase that used a histogram. Also updates rigid body definitions.
450
451
452	:param str GPXfile: full .gpx file name
453	:param dict Histograms: dictionary of histograms as {name:data,...}
454	:param dict Phases: dictionary of phases that use histograms
455	:param dict RigidBodies: dictionary of rigid bodies
456	:param dict CovData: dictionary of refined variables, varyList, & covariance matrix
457	:param bool makeBack: True if new backup of .gpx file is to be made; else use the last one made
458
459	'''
460
461	GPXback = GPXBackup(GPXfile,makeBack)
462	print 'Read from file:',GPXback
463	print 'Save to file :',GPXfile
464	infile = open(GPXback,'rb')
465	outfile = open(GPXfile,'wb')
466	while True:
467	try:
468	data = cPickle.load(infile)
469	except EOFError:
470	break
471	datum = data[0]
472	# print 'read: ',datum[0]
473	if datum[0] == 'Phases':
474	for iphase in range(len(data)):
475	if data[iphase][0] in Phases:
476	phaseName = data[iphase][0]
477	data[iphase][1].update(Phases[phaseName])
478	elif datum[0] == 'Covariance':
479	data[0][1] = CovData
480	elif datum[0] == 'Rigid bodies':
481	data[0][1] = RigidBodies
482	try:
483	histogram = Histograms[datum[0]]
484	# print 'found ',datum[0]
485	data[0][1][1] = histogram['Data']
486	data[0][1][0].update(histogram['Residuals'])
487	for datus in data[1:]:
488	# print ' read: ',datus[0]
489	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
490	datus[1] = histogram[datus[0]]
491	except KeyError:
492	pass
493
494	cPickle.dump(data,outfile,1)
495	infile.close()
496	outfile.close()
497	print 'GPX file save successful'
498
499	def SetSeqResult(GPXfile,Histograms,SeqResult):
500	'''
501	Needs doc string
502
503	:param str GPXfile: full .gpx file name
504	'''
505	GPXback = GPXBackup(GPXfile)
506	print 'Read from file:',GPXback
507	print 'Save to file :',GPXfile
508	infile = open(GPXback,'rb')
509	outfile = open(GPXfile,'wb')
510	while True:
511	try:
512	data = cPickle.load(infile)
513	except EOFError:
514	break
515	datum = data[0]
516	if datum[0] == 'Sequential results':
517	data[0][1] = SeqResult
518	# reset the Copy Next flag, since it should not be needed twice in a row
519	if datum[0] == 'Controls':
520	data[0][1]['Copy2Next'] = False
521	try:
522	histogram = Histograms[datum[0]]
523	data[0][1][1] = list(histogram['Data'])
524	for datus in data[1:]:
525	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
526	datus[1] = histogram[datus[0]]
527	except KeyError:
528	pass
529
530	cPickle.dump(data,outfile,1)
531	infile.close()
532	outfile.close()
533	print 'GPX file save successful'
534
535	def ShowBanner(pFile=None):
536	'Print authorship, copyright and citation notice'
537	print >>pFile,80''
538	print >>pFile,' General Structure Analysis System-II Crystal Structure Refinement'
539	print >>pFile,' by Robert B. Von Dreele & Brian H. Toby'
540	print >>pFile,' Argonne National Laboratory(C), 2010'
541	print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,'
542	print >>pFile,' as Operator of Argonne National Laboratory.'
543	print >>pFile,' Please cite:'
544	print >>pFile,' B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)'
545
546	print >>pFile,80'','\n'
547
548	def ShowControls(Controls,pFile=None,SeqRef=False):
549	'Print controls information'
550	print >>pFile,' Least squares controls:'
551	print >>pFile,' Refinement type: ',Controls['deriv type']
552	if 'Hessian' in Controls['deriv type']:
553	print >>pFile,' Maximum number of cycles:',Controls['max cyc']
554	else:
555	print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
556	print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
557	if SeqRef:
558	print >>pFile,' Sequential refinement controls:'
559	print >>pFile,' Copy of histogram results to next: ',Controls['Copy2Next']
560	print >>pFile,' Process histograms in reverse order: ',Controls['Reverse Seq']
561
562	def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
563	'needs a doc string'
564	# if SGLaue in ['-1','2/m','mmm']:
565	# return #no Pawley symmetry required constraints
566	eqvDict = {}
567	for i,varyI in enumerate(pawleyVary):
568	eqvDict[varyI] = []
569	refI = int(varyI.split(':')[-1])
570	ih,ik,il = PawleyRef[refI][:3]
571	dspI = PawleyRef[refI][4]
572	for varyJ in pawleyVary[i+1:]:
573	refJ = int(varyJ.split(':')[-1])
574	jh,jk,jl = PawleyRef[refJ][:3]
575	dspJ = PawleyRef[refJ][4]
576	if SGLaue in ['4/m','4/mmm']:
577	isum = ih2+ik2
578	jsum = jh2+jk2
579	if abs(il) == abs(jl) and isum == jsum:
580	eqvDict[varyI].append(varyJ)
581	elif SGLaue in ['3R','3mR']:
582	isum = ih2+ik2+il**2
583	jsum = jh2+jk2+jl**2
584	isum2 = ihik+ihil+ik*il
585	jsum2 = jhjk+jhjl+jk*jl
586	if isum == jsum and isum2 == jsum2:
587	eqvDict[varyI].append(varyJ)
588	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
589	isum = ih2+ik2+ih*ik
590	jsum = jh2+jk2+jh*jk
591	if abs(il) == abs(jl) and isum == jsum:
592	eqvDict[varyI].append(varyJ)
593	elif SGLaue in ['m3','m3m']:
594	isum = ih2+ik2+il**2
595	jsum = jh2+jk2+jl**2
596	if isum == jsum:
597	eqvDict[varyI].append(varyJ)
598	elif abs(dspI-dspJ)/dspI < 1.e-4:
599	eqvDict[varyI].append(varyJ)
600	for item in pawleyVary:
601	if eqvDict[item]:
602	for item2 in pawleyVary:
603	if item2 in eqvDict[item]:
604	eqvDict[item2] = []
605	G2mv.StoreEquivalence(item,eqvDict[item])
606
607	def cellVary(pfx,SGData):
608	'needs a doc string'
609	if SGData['SGLaue'] in ['-1',]:
610	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
611	elif SGData['SGLaue'] in ['2/m',]:
612	if SGData['SGUniq'] == 'a':
613	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
614	elif SGData['SGUniq'] == 'b':
615	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
616	else:
617	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
618	elif SGData['SGLaue'] in ['mmm',]:
619	return [pfx+'A0',pfx+'A1',pfx+'A2']
620	elif SGData['SGLaue'] in ['4/m','4/mmm']:
621	return [pfx+'A0',pfx+'A2']
622	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
623	return [pfx+'A0',pfx+'A2']
624	elif SGData['SGLaue'] in ['3R', '3mR']:
625	return [pfx+'A0',pfx+'A3']
626	elif SGData['SGLaue'] in ['m3m','m3']:
627	return [pfx+'A0',]
628
629	################################################################################
630	##### Rigid Body Models and not General.get('doPawley')
631	################################################################################
632
633	def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
634	'needs a doc string'
635
636	def PrintResRBModel(RBModel):
637	atNames = RBModel['atNames']
638	rbRef = RBModel['rbRef']
639	rbSeq = RBModel['rbSeq']
640	print >>pFile,'Residue RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
641	'No. times used: ',RBModel['useCount']
642	print >>pFile,' At name x y z'
643	for name,xyz in zip(atNames,RBModel['rbXYZ']):
644	print >>pFile,' %8s %10.4f %10.4f %10.4f'%(name,xyz[0],xyz[1],xyz[2])
645	print >>pFile,'Orientation defined by:',atNames[rbRef[0]],' -> ',atNames[rbRef[1]], \
646	' & ',atNames[rbRef[0]],' -> ',atNames[rbRef[2]]
647	if rbSeq:
648	for i,rbseq in enumerate(rbSeq):
649	print >>pFile,'Torsion sequence ',i,' Bond: '+atNames[rbseq[0]],' - ', \
650	atNames[rbseq[1]],' riding: ',[atNames[i] for i in rbseq[3]]
651
652	def PrintVecRBModel(RBModel):
653	rbRef = RBModel['rbRef']
654	atTypes = RBModel['rbTypes']
655	print >>pFile,'Vector RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
656	'No. times used: ',RBModel['useCount']
657	for i in range(len(RBModel['VectMag'])):
658	print >>pFile,'Vector no.: ',i,' Magnitude: ', \
659	'%8.4f'%(RBModel['VectMag'][i]),' Refine? ',RBModel['VectRef'][i]
660	print >>pFile,' No. Type vx vy vz'
661	for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
662	print >>pFile,' %d %2s %10.4f %10.4f %10.4f'%(j,name,xyz[0],xyz[1],xyz[2])
663	print >>pFile,' No. Type x y z'
664	for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
665	print >>pFile,' %d %2s %10.4f %10.4f %10.4f'%(i,name,xyz[0],xyz[1],xyz[2])
666	print >>pFile,'Orientation defined by: atom ',rbRef[0],' -> atom ',rbRef[1], \
667	' & atom ',rbRef[0],' -> atom ',rbRef[2]
668	rbVary = []
669	rbDict = {}
670	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
671	if len(rbIds['Vector']):
672	for irb,item in enumerate(rbIds['Vector']):
673	if rigidbodyDict['Vector'][item]['useCount']:
674	RBmags = rigidbodyDict['Vector'][item]['VectMag']
675	RBrefs = rigidbodyDict['Vector'][item]['VectRef']
676	for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
677	pid = '::RBV;'+str(i)+':'+str(irb)
678	rbDict[pid] = mag
679	if ref:
680	rbVary.append(pid)
681	if Print:
682	print >>pFile,'\nVector rigid body model:'
683	PrintVecRBModel(rigidbodyDict['Vector'][item])
684	if len(rbIds['Residue']):
685	for item in rbIds['Residue']:
686	if rigidbodyDict['Residue'][item]['useCount']:
687	if Print:
688	print >>pFile,'\nResidue rigid body model:'
689	PrintResRBModel(rigidbodyDict['Residue'][item])
690	return rbVary,rbDict
691
692	def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
693	'needs a doc string'
694
695	def PrintRBVectandSig(VectRB,VectSig):
696	print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
697	namstr = ' names :'
698	valstr = ' values:'
699	sigstr = ' esds :'
700	for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
701	namstr += '%12s'%('Vect '+str(i))
702	valstr += '%12.4f'%(val)
703	if sig:
704	sigstr += '%12.4f'%(sig)
705	else:
706	sigstr += 12*' '
707	print >>pFile,namstr
708	print >>pFile,valstr
709	print >>pFile,sigstr
710
711	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
712	if not RBIds['Vector']:
713	return
714	for irb,item in enumerate(RBIds['Vector']):
715	if rigidbodyDict['Vector'][item]['useCount']:
716	VectSig = []
717	RBmags = rigidbodyDict['Vector'][item]['VectMag']
718	for i,mag in enumerate(RBmags):
719	name = '::RBV;'+str(i)+':'+str(irb)
720	mag = parmDict[name]
721	if name in sigDict:
722	VectSig.append(sigDict[name])
723	PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)
724
725	################################################################################
726	##### Phase data
727	################################################################################
728
729	def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None):
730	'needs a doc string'
731
732	def PrintFFtable(FFtable):
733	print >>pFile,'\n X-ray scattering factors:'
734	print >>pFile,' Symbol fa fb fc'
735	print >>pFile,99*'-'
736	for Ename in FFtable:
737	ffdata = FFtable[Ename]
738	fa = ffdata['fa']
739	fb = ffdata['fb']
740	print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \
741	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
742
743	def PrintBLtable(BLtable):
744	print >>pFile,'\n Neutron scattering factors:'
745	print >>pFile,' Symbol isotope mass b resonant terms'
746	print >>pFile,99*'-'
747	for Ename in BLtable:
748	bldata = BLtable[Ename]
749	isotope = bldata[0]
750	mass = bldata[1][0]
751	blen = bldata[1][1]
752	bres = []
753	if len(bldata[1]) > 2:
754	bres = bldata[1][2:]
755	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
756	for item in bres:
757	line += '%10.5g'%(item)
758	print >>pFile,line
759
760	def PrintRBObjects(resRBData,vecRBData):
761
762	def PrintRBThermals():
763	tlstr = ['11','22','33','12','13','23']
764	sstr = ['12','13','21','23','31','32','AA','BB']
765	TLS = RB['ThermalMotion'][1]
766	TLSvar = RB['ThermalMotion'][2]
767	if 'T' in RB['ThermalMotion'][0]:
768	print >>pFile,'TLS data'
769	text = ''
770	for i in range(6):
771	text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
772	print >>pFile,text
773	if 'L' in RB['ThermalMotion'][0]:
774	text = ''
775	for i in range(6,12):
776	text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
777	print >>pFile,text
778	if 'S' in RB['ThermalMotion'][0]:
779	text = ''
780	for i in range(12,20):
781	text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
782	print >>pFile,text
783	if 'U' in RB['ThermalMotion'][0]:
784	print >>pFile,'Uiso data'
785	text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])
786
787	if len(resRBData):
788	for RB in resRBData:
789	Oxyz = RB['Orig'][0]
790	Qrijk = RB['Orient'][0]
791	Angle = 2.0*acosd(Qrijk[0])
792	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
793	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
794	print >>pFile,'Orientation angle,vector:', \
795	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
796	Torsions = RB['Torsions']
797	if len(Torsions):
798	text = 'Torsions: '
799	for torsion in Torsions:
800	text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
801	print >>pFile,text
802	PrintRBThermals()
803	if len(vecRBData):
804	for RB in vecRBData:
805	Oxyz = RB['Orig'][0]
806	Qrijk = RB['Orient'][0]
807	Angle = 2.0*acosd(Qrijk[0])
808	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
809	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
810	print >>pFile,'Orientation angle,vector:', \
811	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
812	PrintRBThermals()
813
814	def PrintAtoms(General,Atoms):
815	cx,ct,cs,cia = General['AtomPtrs']
816	print >>pFile,'\n Atoms:'
817	line = ' name type refine? x y z '+ \
818	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
819	if General['Type'] == 'magnetic':
820	line += ' Mx My Mz'
821	elif General['Type'] == 'macromolecular':
822	line = ' res no residue chain'+line
823	print >>pFile,line
824	if General['Type'] == 'nuclear':
825	print >>pFile,135*'-'
826	for i,at in enumerate(Atoms):
827	line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
828	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
829	if at[cia] == 'I':
830	line += '%8.4f'%(at[cia+1])+48*' '
831	else:
832	line += 8*' '
833	for j in range(6):
834	line += '%8.4f'%(at[cia+1+j])
835	print >>pFile,line
836	elif General['Type'] == 'macromolecular':
837	print >>pFile,135*'-'
838	for i,at in enumerate(Atoms):
839	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])+ \
840	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
841	if at[cia] == 'I':
842	line += '%8.4f'%(at[cia+1])+48*' '
843	else:
844	line += 8*' '
845	for j in range(6):
846	line += '%8.4f'%(at[cia+1+j])
847	print >>pFile,line
848
849	def PrintTexture(textureData):
850	topstr = '\n Spherical harmonics texture: Order:' + \
851	str(textureData['Order'])
852	if textureData['Order']:
853	print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
854	else:
855	print >>pFile,topstr
856	return
857	names = ['omega','chi','phi']
858	line = '\n'
859	for name in names:
860	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
861	print >>pFile,line
862	print >>pFile,'\n Texture coefficients:'
863	ptlbls = ' names :'
864	ptstr = ' values:'
865	SHcoeff = textureData['SH Coeff'][1]
866	for item in SHcoeff:
867	ptlbls += '%12s'%(item)
868	ptstr += '%12.4f'%(SHcoeff[item])
869	print >>pFile,ptlbls
870	print >>pFile,ptstr
871
872	def MakeRBParms(rbKey,phaseVary,phaseDict):
873	rbid = str(rbids.index(RB['RBId']))
874	pfxRB = pfx+'RB'+rbKey+'P'
875	pstr = ['x','y','z']
876	ostr = ['a','i','j','k']
877	for i in range(3):
878	name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
879	phaseDict[name] = RB['Orig'][0][i]
880	if RB['Orig'][1]:
881	phaseVary += [name,]
882	pfxRB = pfx+'RB'+rbKey+'O'
883	for i in range(4):
884	name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
885	phaseDict[name] = RB['Orient'][0][i]
886	if RB['Orient'][1] == 'AV' and i:
887	phaseVary += [name,]
888	elif RB['Orient'][1] == 'A' and not i:
889	phaseVary += [name,]
890
891	def MakeRBThermals(rbKey,phaseVary,phaseDict):
892	rbid = str(rbids.index(RB['RBId']))
893	tlstr = ['11','22','33','12','13','23']
894	sstr = ['12','13','21','23','31','32','AA','BB']
895	if 'T' in RB['ThermalMotion'][0]:
896	pfxRB = pfx+'RB'+rbKey+'T'
897	for i in range(6):
898	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
899	phaseDict[name] = RB['ThermalMotion'][1][i]
900	if RB['ThermalMotion'][2][i]:
901	phaseVary += [name,]
902	if 'L' in RB['ThermalMotion'][0]:
903	pfxRB = pfx+'RB'+rbKey+'L'
904	for i in range(6):
905	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
906	phaseDict[name] = RB['ThermalMotion'][1][i+6]
907	if RB['ThermalMotion'][2][i+6]:
908	phaseVary += [name,]
909	if 'S' in RB['ThermalMotion'][0]:
910	pfxRB = pfx+'RB'+rbKey+'S'
911	for i in range(8):
912	name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
913	phaseDict[name] = RB['ThermalMotion'][1][i+12]
914	if RB['ThermalMotion'][2][i+12]:
915	phaseVary += [name,]
916	if 'U' in RB['ThermalMotion'][0]:
917	name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
918	phaseDict[name] = RB['ThermalMotion'][1][0]
919	if RB['ThermalMotion'][2][0]:
920	phaseVary += [name,]
921
922	def MakeRBTorsions(rbKey,phaseVary,phaseDict):
923	rbid = str(rbids.index(RB['RBId']))
924	pfxRB = pfx+'RB'+rbKey+'Tr;'
925	for i,torsion in enumerate(RB['Torsions']):
926	name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
927	phaseDict[name] = torsion[0]
928	if torsion[1]:
929	phaseVary += [name,]
930
931	if Print:
932	print >>pFile,'\n Phases:'
933	phaseVary = []
934	phaseDict = {}
935	phaseConstr = {}
936	pawleyLookup = {}
937	FFtables = {} #scattering factors - xrays
938	BLtables = {} # neutrons
939	Natoms = {}
940	AtMults = {}
941	AtIA = {}
942	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
943	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
944	atomIndx = {}
945	for name in PhaseData:
946	General = PhaseData[name]['General']
947	pId = PhaseData[name]['pId']
948	pfx = str(pId)+'::'
949	FFtable = G2el.GetFFtable(General['AtomTypes'])
950	BLtable = G2el.GetBLtable(General)
951	FFtables.update(FFtable)
952	BLtables.update(BLtable)
953	Atoms = PhaseData[name]['Atoms']
954	AtLookup = G2mth.FillAtomLookUp(Atoms)
955	PawleyRef = PhaseData[name].get('Pawley ref',[])
956	SGData = General['SGData']
957	SGtext = G2spc.SGPrint(SGData)
958	cell = General['Cell']
959	A = G2lat.cell2A(cell[1:7])
960	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
961	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
962	if cell[0]:
963	phaseVary += cellVary(pfx,SGData)
964	resRBData = PhaseData[name]['RBModels'].get('Residue',[])
965	if resRBData:
966	rbids = rbIds['Residue'] #NB: used in the MakeRB routines
967	for iRB,RB in enumerate(resRBData):
968	MakeRBParms('R',phaseVary,phaseDict)
969	MakeRBThermals('R',phaseVary,phaseDict)
970	MakeRBTorsions('R',phaseVary,phaseDict)
971
972	vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
973	if vecRBData:
974	rbids = rbIds['Vector'] #NB: used in the MakeRB routines
975	for iRB,RB in enumerate(vecRBData):
976	MakeRBParms('V',phaseVary,phaseDict)
977	MakeRBThermals('V',phaseVary,phaseDict)
978
979	Natoms[pfx] = 0
980	if Atoms and not General.get('doPawley'):
981	cx,ct,cs,cia = General['AtomPtrs']
982	if General['Type'] in ['nuclear','macromolecular']:
983	Natoms[pfx] = len(Atoms)
984	for i,at in enumerate(Atoms):
985	atomIndx[at[-1]] = [pfx,i] #lookup table for restraints
986	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
987	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
988	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
989	pfx+'AI/A:'+str(i):at[cia],})
990	if at[cia] == 'I':
991	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
992	else:
993	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
994	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
995	if 'F' in at[ct+1]:
996	phaseVary.append(pfx+'Afrac:'+str(i))
997	if 'X' in at[ct+1]:
998	xId,xCoef = G2spc.GetCSxinel(at[cs])
999	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1000	equivs = [[],[],[]]
1001	for j in range(3):
1002	if xId[j] > 0:
1003	phaseVary.append(names[j])
1004	equivs[xId[j]-1].append([names[j],xCoef[j]])
1005	for equiv in equivs:
1006	if len(equiv) > 1:
1007	name = equiv[0][0]
1008	for eqv in equiv[1:]:
1009	G2mv.StoreEquivalence(name,(eqv,))
1010	if 'U' in at[ct+1]:
1011	if at[cia] == 'I':
1012	phaseVary.append(pfx+'AUiso:'+str(i))
1013	else:
1014	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1015	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1016	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1017	equivs = [[],[],[],[],[],[]]
1018	for j in range(6):
1019	if uId[j] > 0:
1020	phaseVary.append(names[j])
1021	equivs[uId[j]-1].append([names[j],uCoef[j]])
1022	for equiv in equivs:
1023	if len(equiv) > 1:
1024	name = equiv[0][0]
1025	for eqv in equiv[1:]:
1026	G2mv.StoreEquivalence(name,(eqv,))
1027	# elif General['Type'] == 'magnetic':
1028	# elif General['Type'] == 'macromolecular':
1029	textureData = General['SH Texture']
1030	if textureData['Order']:
1031	phaseDict[pfx+'SHorder'] = textureData['Order']
1032	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1033	for item in ['omega','chi','phi']:
1034	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1035	if textureData['Sample '+item][0]:
1036	phaseVary.append(pfx+'SH '+item)
1037	for item in textureData['SH Coeff'][1]:
1038	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1039	if textureData['SH Coeff'][0]:
1040	phaseVary.append(pfx+item)
1041
1042	if Print:
1043	print >>pFile,'\n Phase name: ',General['Name']
1044	print >>pFile,135*'-'
1045	PrintFFtable(FFtable)
1046	PrintBLtable(BLtable)
1047	print >>pFile,''
1048	for line in SGtext: print >>pFile,line
1049	PrintRBObjects(resRBData,vecRBData)
1050	PrintAtoms(General,Atoms)
1051	print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
1052	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
1053	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
1054	PrintTexture(textureData)
1055	if name in RestraintDict:
1056	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1057	textureData,RestraintDict[name],pFile)
1058
1059	elif PawleyRef:
1060	pawleyVary = []
1061	for i,refl in enumerate(PawleyRef):
1062	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
1063	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1064	if refl[5]:
1065	pawleyVary.append(pfx+'PWLref:'+str(i))
1066	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
1067	phaseVary += pawleyVary
1068
1069	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1070
1071	def cellFill(pfx,SGData,parmDict,sigDict):
1072	'''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1073	from the parameter and sig dictionaries.
1074
1075	:param str pfx: parameter prefix ("n::", where n is a phase number)
1076	:param dict SGdata: a symmetry object
1077	:param dict parmDict: a dictionary of parameters
1078	:param dict sigDict: a dictionary of uncertainties on parameters
1079
1080	:returns: A,sigA where each is a list of six terms with the A terms
1081	'''
1082	if SGData['SGLaue'] in ['-1',]:
1083	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1084	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1085	elif SGData['SGLaue'] in ['2/m',]:
1086	if SGData['SGUniq'] == 'a':
1087	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1088	parmDict[pfx+'A3'],0,0]
1089	elif SGData['SGUniq'] == 'b':
1090	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1091	0,parmDict[pfx+'A4'],0]
1092	else:
1093	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1094	0,0,parmDict[pfx+'A5']]
1095	elif SGData['SGLaue'] in ['mmm',]:
1096	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1097	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1098	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1099	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1100	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1101	parmDict[pfx+'A0'],0,0]
1102	elif SGData['SGLaue'] in ['3R', '3mR']:
1103	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1104	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1105	elif SGData['SGLaue'] in ['m3m','m3']:
1106	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1107
1108	try:
1109	if SGData['SGLaue'] in ['-1',]:
1110	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1111	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1112	elif SGData['SGLaue'] in ['2/m',]:
1113	if SGData['SGUniq'] == 'a':
1114	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1115	sigDict[pfx+'A3'],0,0]
1116	elif SGData['SGUniq'] == 'b':
1117	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1118	0,sigDict[pfx+'A4'],0]
1119	else:
1120	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1121	0,0,sigDict[pfx+'A5']]
1122	elif SGData['SGLaue'] in ['mmm',]:
1123	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1124	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1125	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1126	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1127	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1128	elif SGData['SGLaue'] in ['3R', '3mR']:
1129	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1130	elif SGData['SGLaue'] in ['m3m','m3']:
1131	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1132	except KeyError:
1133	sigA = [0,0,0,0,0,0]
1134
1135	return A,sigA
1136
1137	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1138	'needs a doc string'
1139	if phaseRest:
1140	Amat = G2lat.cell2AB(cell)[0]
1141	cx,ct,cs = AtPtrs[:3]
1142	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1143	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1144	['ChemComp','Sites'],['Texture','HKLs']]
1145	for name,rest in names:
1146	itemRest = phaseRest[name]
1147	if itemRest[rest] and itemRest['Use']:
1148	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1149	if name in ['Bond','Angle','Plane','Chiral']:
1150	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
1151	for indx,ops,obs,esd in itemRest[rest]:
1152	try:
1153	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1154	AtName = ''
1155	for i,Aname in enumerate(AtNames):
1156	AtName += Aname
1157	if ops[i] == '1':
1158	AtName += '-'
1159	else:
1160	AtName += '+('+ops[i]+')-'
1161	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1162	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1163	if name == 'Bond':
1164	calc = G2mth.getRestDist(XYZ,Amat)
1165	elif name == 'Angle':
1166	calc = G2mth.getRestAngle(XYZ,Amat)
1167	elif name == 'Plane':
1168	calc = G2mth.getRestPlane(XYZ,Amat)
1169	elif name == 'Chiral':
1170	calc = G2mth.getRestChiral(XYZ,Amat)
1171	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1172	except KeyError:
1173	del itemRest[rest]
1174	elif name in ['Torsion','Rama']:
1175	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
1176	coeffDict = itemRest['Coeff']
1177	for indx,ops,cofName,esd in itemRest[rest]:
1178	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1179	AtName = ''
1180	for i,Aname in enumerate(AtNames):
1181	AtName += Aname+'+('+ops[i]+')-'
1182	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1183	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1184	if name == 'Torsion':
1185	tor = G2mth.getRestTorsion(XYZ,Amat)
1186	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1187	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1188	else:
1189	phi,psi = G2mth.getRestRama(XYZ,Amat)
1190	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1191	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])
1192	elif name == 'ChemComp':
1193	print >>pFile,' atoms mul*frac factor prod'
1194	for indx,factors,obs,esd in itemRest[rest]:
1195	try:
1196	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1197	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1198	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1199	mulfrac = mul*frac
1200	calcs = mulfracfactors
1201	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1202	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1203	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1204	except KeyError:
1205	del itemRest[rest]
1206	elif name == 'Texture' and textureData['Order']:
1207	Start = False
1208	SHCoef = textureData['SH Coeff'][1]
1209	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1210	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1211	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
1212	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1213	PH = np.array(hkl)
1214	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1215	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1216	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1217	Z = ma.masked_greater(Z,0.0)
1218	num = ma.count(Z)
1219	sum = 0
1220	if num:
1221	sum = np.sum(Z)
1222	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)
1223
1224	def getCellEsd(pfx,SGData,A,covData):
1225	'needs a doc string'
1226	dpr = 180./np.pi
1227	rVsq = G2lat.calc_rVsq(A)
1228	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1229	cell = np.array(G2lat.Gmat2cell(g)) #real cell
1230	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1231	scos = cosd(cellst[3:6])
1232	ssin = sind(cellst[3:6])
1233	scot = scos/ssin
1234	rcos = cosd(cell[3:6])
1235	rsin = sind(cell[3:6])
1236	rcot = rcos/rsin
1237	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1238	varyList = covData['varyList']
1239	covMatrix = covData['covMatrix']
1240	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1241	Ax = np.array(A)
1242	Ax[3:] /= 2.
1243	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
1244	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
1245	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1246	Vol = 1/np.sqrt(rVsq)
1247	sigVol = Vol*3np.sqrt(srcvlsq)/2.
1248	R123 = Ax[0]Ax[1]Ax[2]
1249	dsasdg = np.zeros((3,6))
1250	dadg = np.zeros((6,6))
1251	for i0 in range(3): #0 1 2
1252	i1 = (i0+1)%3 #1 2 0
1253	i2 = (i1+1)%3 #2 0 1
1254	i3 = 5-i2 #3 5 4
1255	i4 = 5-i1 #4 3 5
1256	i5 = 5-i0 #5 4 3
1257	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
1258	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
1259	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1260	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
1261	denom = np.sqrt(denmsq)
1262	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
1263	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
1264	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
1265	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
1266	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
1267	dadg[i5][i5] = -Ax[i0]/denom
1268	for i0 in range(3):
1269	i1 = (i0+1)%3
1270	i2 = (i1+1)%3
1271	i3 = 5-i2
1272	for ij in range(6):
1273	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1274	if ij == i0:
1275	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1276	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
1277	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1278	var = np.diag(sigMat)
1279	CS = np.where(var>0.,np.sqrt(var),0.)
1280	return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
1281
1282	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1283	'needs a doc string'
1284
1285	def PrintAtomsAndSig(General,Atoms,atomsSig):
1286	print >>pFile,'\n Atoms:'
1287	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1288	if General['Type'] == 'magnetic':
1289	line += ' Mx My Mz'
1290	elif General['Type'] == 'macromolecular':
1291	line = ' res no residue chain '+line
1292	print >>pFile,line
1293	if General['Type'] == 'nuclear':
1294	print >>pFile,135*'-'
1295	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1296	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1297	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
1298	for atyp in General['AtomTypes']: #zero composition
1299	General['NoAtoms'][atyp] = 0.
1300	for i,at in enumerate(Atoms):
1301	General['NoAtoms'][at[1]] += at[6]*float(at[8]) #new composition
1302	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1303	valstr = ' values:'
1304	sigstr = ' sig :'
1305	for ind in [3,4,5,6]:
1306	sigind = str(i)+':'+str(ind)
1307	valstr += fmt[ind]%(at[ind])
1308	if sigind in atomsSig:
1309	sigstr += fmt[ind]%(atomsSig[sigind])
1310	else:
1311	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1312	if at[9] == 'I':
1313	valstr += fmt[10]%(at[10])
1314	if str(i)+':10' in atomsSig:
1315	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1316	else:
1317	sigstr += 8*' '
1318	else:
1319	valstr += 8*' '
1320	sigstr += 8*' '
1321	for ind in [11,12,13,14,15,16]:
1322	sigind = str(i)+':'+str(ind)
1323	valstr += fmt[ind]%(at[ind])
1324	if sigind in atomsSig:
1325	sigstr += fmt[ind]%(atomsSig[sigind])
1326	else:
1327	sigstr += 8*' '
1328	print >>pFile,name
1329	print >>pFile,valstr
1330	print >>pFile,sigstr
1331
1332	def PrintRBObjPOAndSig(rbfx,rbsx):
1333	namstr = ' names :'
1334	valstr = ' values:'
1335	sigstr = ' esds :'
1336	for i,px in enumerate(['Px:','Py:','Pz:']):
1337	name = pfx+rbfx+px+rbsx
1338	namstr += '%12s'%('Pos '+px[1])
1339	valstr += '%12.5f'%(parmDict[name])
1340	if name in sigDict:
1341	sigstr += '%12.5f'%(sigDict[name])
1342	else:
1343	sigstr += 12*' '
1344	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1345	name = pfx+rbfx+po+rbsx
1346	namstr += '%12s'%('Ori '+po[1])
1347	valstr += '%12.5f'%(parmDict[name])
1348	if name in sigDict:
1349	sigstr += '%12.5f'%(sigDict[name])
1350	else:
1351	sigstr += 12*' '
1352	print >>pFile,namstr
1353	print >>pFile,valstr
1354	print >>pFile,sigstr
1355
1356	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1357	namstr = ' names :'
1358	valstr = ' values:'
1359	sigstr = ' esds :'
1360	if 'N' not in TLS:
1361	print >>pFile,' Thermal motion:'
1362	if 'T' in TLS:
1363	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1364	name = pfx+rbfx+pt+rbsx
1365	namstr += '%12s'%(pt[:3])
1366	valstr += '%12.5f'%(parmDict[name])
1367	if name in sigDict:
1368	sigstr += '%12.5f'%(sigDict[name])
1369	else:
1370	sigstr += 12*' '
1371	print >>pFile,namstr
1372	print >>pFile,valstr
1373	print >>pFile,sigstr
1374	if 'L' in TLS:
1375	namstr = ' names :'
1376	valstr = ' values:'
1377	sigstr = ' esds :'
1378	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1379	name = pfx+rbfx+pt+rbsx
1380	namstr += '%12s'%(pt[:3])
1381	valstr += '%12.3f'%(parmDict[name])
1382	if name in sigDict:
1383	sigstr += '%12.3f'%(sigDict[name])
1384	else:
1385	sigstr += 12*' '
1386	print >>pFile,namstr
1387	print >>pFile,valstr
1388	print >>pFile,sigstr
1389	if 'S' in TLS:
1390	namstr = ' names :'
1391	valstr = ' values:'
1392	sigstr = ' esds :'
1393	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1394	name = pfx+rbfx+pt+rbsx
1395	namstr += '%12s'%(pt[:3])
1396	valstr += '%12.4f'%(parmDict[name])
1397	if name in sigDict:
1398	sigstr += '%12.4f'%(sigDict[name])
1399	else:
1400	sigstr += 12*' '
1401	print >>pFile,namstr
1402	print >>pFile,valstr
1403	print >>pFile,sigstr
1404	if 'U' in TLS:
1405	name = pfx+rbfx+'U:'+rbsx
1406	namstr = ' names :'+'%12s'%('Uiso')
1407	valstr = ' values:'+'%12.5f'%(parmDict[name])
1408	if name in sigDict:
1409	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1410	else:
1411	sigstr = ' esds :'+12*' '
1412	print >>pFile,namstr
1413	print >>pFile,valstr
1414	print >>pFile,sigstr
1415
1416	def PrintRBObjTorAndSig(rbsx):
1417	namstr = ' names :'
1418	valstr = ' values:'
1419	sigstr = ' esds :'
1420	nTors = len(RBObj['Torsions'])
1421	if nTors:
1422	print >>pFile,' Torsions:'
1423	for it in range(nTors):
1424	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1425	namstr += '%12s'%('Tor'+str(it))
1426	valstr += '%12.4f'%(parmDict[name])
1427	if name in sigDict:
1428	sigstr += '%12.4f'%(sigDict[name])
1429	print >>pFile,namstr
1430	print >>pFile,valstr
1431	print >>pFile,sigstr
1432
1433	def PrintSHtextureAndSig(textureData,SHtextureSig):
1434	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1435	names = ['omega','chi','phi']
1436	namstr = ' names :'
1437	ptstr = ' values:'
1438	sigstr = ' esds :'
1439	for name in names:
1440	namstr += '%12s'%(name)
1441	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1442	if 'Sample '+name in SHtextureSig:
1443	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1444	else:
1445	sigstr += 12*' '
1446	print >>pFile,namstr
1447	print >>pFile,ptstr
1448	print >>pFile,sigstr
1449	print >>pFile,'\n Texture coefficients:'
1450	namstr = ' names :'
1451	ptstr = ' values:'
1452	sigstr = ' esds :'
1453	SHcoeff = textureData['SH Coeff'][1]
1454	for name in SHcoeff:
1455	namstr += '%12s'%(name)
1456	ptstr += '%12.3f'%(SHcoeff[name])
1457	if name in SHtextureSig:
1458	sigstr += '%12.3f'%(SHtextureSig[name])
1459	else:
1460	sigstr += 12*' '
1461	print >>pFile,namstr
1462	print >>pFile,ptstr
1463	print >>pFile,sigstr
1464
1465	print >>pFile,'\n Phases:'
1466	for phase in Phases:
1467	print >>pFile,' Result for phase: ',phase
1468	Phase = Phases[phase]
1469	General = Phase['General']
1470	SGData = General['SGData']
1471	Atoms = Phase['Atoms']
1472	AtLookup = G2mth.FillAtomLookUp(Atoms)
1473	cell = General['Cell']
1474	pId = Phase['pId']
1475	pfx = str(pId)+'::'
1476	if cell[0]:
1477	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1478	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1479	print >>pFile,' Reciprocal metric tensor: '
1480	ptfmt = "%15.9f"
1481	names = ['A11','A22','A33','A12','A13','A23']
1482	namstr = ' names :'
1483	ptstr = ' values:'
1484	sigstr = ' esds :'
1485	for name,a,siga in zip(names,A,sigA):
1486	namstr += '%15s'%(name)
1487	ptstr += ptfmt%(a)
1488	if siga:
1489	sigstr += ptfmt%(siga)
1490	else:
1491	sigstr += 15*' '
1492	print >>pFile,namstr
1493	print >>pFile,ptstr
1494	print >>pFile,sigstr
1495	cell[1:7] = G2lat.A2cell(A)
1496	cell[7] = G2lat.calc_V(A)
1497	print >>pFile,' New unit cell:'
1498	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1499	names = ['a','b','c','alpha','beta','gamma','Volume']
1500	namstr = ' names :'
1501	ptstr = ' values:'
1502	sigstr = ' esds :'
1503	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1504	namstr += '%12s'%(name)
1505	ptstr += fmt%(a)
1506	if siga:
1507	sigstr += fmt%(siga)
1508	else:
1509	sigstr += 12*' '
1510	print >>pFile,namstr
1511	print >>pFile,ptstr
1512	print >>pFile,sigstr
1513
1514	General['Mass'] = 0.
1515	if Phase['General'].get('doPawley'):
1516	pawleyRef = Phase['Pawley ref']
1517	for i,refl in enumerate(pawleyRef):
1518	key = pfx+'PWLref:'+str(i)
1519	refl[6] = parmDict[key]
1520	if key in sigDict:
1521	refl[7] = sigDict[key]
1522	else:
1523	refl[7] = 0
1524	else:
1525	VRBIds = RBIds['Vector']
1526	RRBIds = RBIds['Residue']
1527	RBModels = Phase['RBModels']
1528	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1529	jrb = VRBIds.index(RBObj['RBId'])
1530	rbsx = str(irb)+':'+str(jrb)
1531	print >>pFile,' Vector rigid body parameters:'
1532	PrintRBObjPOAndSig('RBV',rbsx)
1533	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1534	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1535	jrb = RRBIds.index(RBObj['RBId'])
1536	rbsx = str(irb)+':'+str(jrb)
1537	print >>pFile,' Residue rigid body parameters:'
1538	PrintRBObjPOAndSig('RBR',rbsx)
1539	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1540	PrintRBObjTorAndSig(rbsx)
1541	atomsSig = {}
1542	if General['Type'] == 'nuclear': #this needs macromolecular variant!
1543	for i,at in enumerate(Atoms):
1544	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1545	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1546	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1547	for ind in [3,4,5,6]:
1548	at[ind] = parmDict[names[ind]]
1549	if ind in [3,4,5]:
1550	name = names[ind].replace('A','dA')
1551	else:
1552	name = names[ind]
1553	if name in sigDict:
1554	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1555	if at[9] == 'I':
1556	at[10] = parmDict[names[10]]
1557	if names[10] in sigDict:
1558	atomsSig[str(i)+':10'] = sigDict[names[10]]
1559	else:
1560	for ind in [11,12,13,14,15,16]:
1561	at[ind] = parmDict[names[ind]]
1562	if names[ind] in sigDict:
1563	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1564	ind = General['AtomTypes'].index(at[1])
1565	General['Mass'] += General['AtomMass'][ind]at[6]at[8]
1566	PrintAtomsAndSig(General,Atoms,atomsSig)
1567
1568	textureData = General['SH Texture']
1569	if textureData['Order']:
1570	SHtextureSig = {}
1571	for name in ['omega','chi','phi']:
1572	aname = pfx+'SH '+name
1573	textureData['Sample '+name][1] = parmDict[aname]
1574	if aname in sigDict:
1575	SHtextureSig['Sample '+name] = sigDict[aname]
1576	for name in textureData['SH Coeff'][1]:
1577	aname = pfx+name
1578	textureData['SH Coeff'][1][name] = parmDict[aname]
1579	if aname in sigDict:
1580	SHtextureSig[name] = sigDict[aname]
1581	PrintSHtextureAndSig(textureData,SHtextureSig)
1582	if phase in RestraintDict:
1583	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1584	textureData,RestraintDict[phase],pFile)
1585
1586	################################################################################
1587	##### Histogram & Phase data
1588	################################################################################
1589
1590	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1591	'''Loads the HAP histogram/phase information into dicts
1592
1593	:param dict Phases: phase information
1594	:param dict Histograms: Histogram information
1595	:param bool Print: prints information as it is read
1596	:param file pFile: file object to print to (the default, None causes printing to the console)
1597	:param bool resetRefList: Should the contents of the Reflection List be initialized
1598	on loading. The default, True, initializes the Reflection List as it is loaded.
1599
1600	:returns: (hapVary,hapDict,controlDict)
1601	* hapVary: list of refined variables
1602	* hapDict: dict with refined variables and their values
1603	* controlDict: dict with computation controls (?)
1604	'''
1605
1606	def PrintSize(hapData):
1607	if hapData[0] in ['isotropic','uniaxial']:
1608	line = '\n Size model : %9s'%(hapData[0])
1609	line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1610	if hapData[0] == 'uniaxial':
1611	line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1612	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1613	print >>pFile,line
1614	else:
1615	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1616	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1617	Snames = ['S11','S22','S33','S12','S13','S23']
1618	ptlbls = ' names :'
1619	ptstr = ' values:'
1620	varstr = ' refine:'
1621	for i,name in enumerate(Snames):
1622	ptlbls += '%12s' % (name)
1623	ptstr += '%12.6f' % (hapData[4][i])
1624	varstr += '%12s' % (str(hapData[5][i]))
1625	print >>pFile,ptlbls
1626	print >>pFile,ptstr
1627	print >>pFile,varstr
1628
1629	def PrintMuStrain(hapData,SGData):
1630	if hapData[0] in ['isotropic','uniaxial']:
1631	line = '\n Mustrain model: %9s'%(hapData[0])
1632	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1633	if hapData[0] == 'uniaxial':
1634	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1635	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1636	print >>pFile,line
1637	else:
1638	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1639	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1640	Snames = G2spc.MustrainNames(SGData)
1641	ptlbls = ' names :'
1642	ptstr = ' values:'
1643	varstr = ' refine:'
1644	for i,name in enumerate(Snames):
1645	ptlbls += '%12s' % (name)
1646	ptstr += '%12.6f' % (hapData[4][i])
1647	varstr += '%12s' % (str(hapData[5][i]))
1648	print >>pFile,ptlbls
1649	print >>pFile,ptstr
1650	print >>pFile,varstr
1651
1652	def PrintHStrain(hapData,SGData):
1653	print >>pFile,'\n Hydrostatic/elastic strain: '
1654	Hsnames = G2spc.HStrainNames(SGData)
1655	ptlbls = ' names :'
1656	ptstr = ' values:'
1657	varstr = ' refine:'
1658	for i,name in enumerate(Hsnames):
1659	ptlbls += '%12s' % (name)
1660	ptstr += '%12.6f' % (hapData[0][i])
1661	varstr += '%12s' % (str(hapData[1][i]))
1662	print >>pFile,ptlbls
1663	print >>pFile,ptstr
1664	print >>pFile,varstr
1665
1666	def PrintSHPO(hapData):
1667	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1668	str(hapData[4])+' Refine? '+str(hapData[2])
1669	ptlbls = ' names :'
1670	ptstr = ' values:'
1671	for item in hapData[5]:
1672	ptlbls += '%12s'%(item)
1673	ptstr += '%12.3f'%(hapData[5][item])
1674	print >>pFile,ptlbls
1675	print >>pFile,ptstr
1676
1677	def PrintBabinet(hapData):
1678	print >>pFile,'\n Babinet form factor modification: '
1679	ptlbls = ' names :'
1680	ptstr = ' values:'
1681	varstr = ' refine:'
1682	for name in ['BabA','BabU']:
1683	ptlbls += '%12s' % (name)
1684	ptstr += '%12.6f' % (hapData[name][0])
1685	varstr += '%12s' % (str(hapData[name][1]))
1686	print >>pFile,ptlbls
1687	print >>pFile,ptstr
1688	print >>pFile,varstr
1689
1690
1691	hapDict = {}
1692	hapVary = []
1693	controlDict = {}
1694	poType = {}
1695	poAxes = {}
1696	spAxes = {}
1697	spType = {}
1698
1699	for phase in Phases:
1700	HistoPhase = Phases[phase]['Histograms']
1701	SGData = Phases[phase]['General']['SGData']
1702	cell = Phases[phase]['General']['Cell'][1:7]
1703	A = G2lat.cell2A(cell)
1704	pId = Phases[phase]['pId']
1705	histoList = HistoPhase.keys()
1706	histoList.sort()
1707	for histogram in histoList:
1708	try:
1709	Histogram = Histograms[histogram]
1710	except KeyError:
1711	#skip if histogram not included e.g. in a sequential refinement
1712	continue
1713	hapData = HistoPhase[histogram]
1714	hId = Histogram['hId']
1715	if 'PWDR' in histogram:
1716	limits = Histogram['Limits'][1]
1717	inst = Histogram['Instrument Parameters'][0]
1718	Zero = inst['Zero'][1]
1719	if 'C' in inst['Type'][1]:
1720	try:
1721	wave = inst['Lam'][1]
1722	except KeyError:
1723	wave = inst['Lam1'][1]
1724	dmin = wave/(2.0*sind(limits[1]/2.0))
1725	pfx = str(pId)+':'+str(hId)+':'
1726	for item in ['Scale','Extinction']:
1727	hapDict[pfx+item] = hapData[item][0]
1728	if hapData[item][1]:
1729	hapVary.append(pfx+item)
1730	names = G2spc.HStrainNames(SGData)
1731	for i,name in enumerate(names):
1732	hapDict[pfx+name] = hapData['HStrain'][0][i]
1733	if hapData['HStrain'][1][i]:
1734	hapVary.append(pfx+name)
1735	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1736	if hapData['Pref.Ori.'][0] == 'MD':
1737	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1738	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1739	if hapData['Pref.Ori.'][2]:
1740	hapVary.append(pfx+'MD')
1741	else: #'SH' spherical harmonics
1742	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1743	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1744	for item in hapData['Pref.Ori.'][5]:
1745	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1746	if hapData['Pref.Ori.'][2]:
1747	hapVary.append(pfx+item)
1748	for item in ['Mustrain','Size']:
1749	controlDict[pfx+item+'Type'] = hapData[item][0]
1750	hapDict[pfx+item+';mx'] = hapData[item][1][2]
1751	if hapData[item][2][2]:
1752	hapVary.append(pfx+item+';mx')
1753	if hapData[item][0] in ['isotropic','uniaxial']:
1754	hapDict[pfx+item+';i'] = hapData[item][1][0]
1755	if hapData[item][2][0]:
1756	hapVary.append(pfx+item+';i')
1757	if hapData[item][0] == 'uniaxial':
1758	controlDict[pfx+item+'Axis'] = hapData[item][3]
1759	hapDict[pfx+item+';a'] = hapData[item][1][1]
1760	if hapData[item][2][1]:
1761	hapVary.append(pfx+item+';a')
1762	else: #generalized for mustrain or ellipsoidal for size
1763	Nterms = len(hapData[item][4])
1764	if item == 'Mustrain':
1765	names = G2spc.MustrainNames(SGData)
1766	pwrs = []
1767	for name in names:
1768	h,k,l = name[1:]
1769	pwrs.append([int(h),int(k),int(l)])
1770	controlDict[pfx+'MuPwrs'] = pwrs
1771	for i in range(Nterms):
1772	sfx = ':'+str(i)
1773	hapDict[pfx+item+sfx] = hapData[item][4][i]
1774	if hapData[item][5][i]:
1775	hapVary.append(pfx+item+sfx)
1776	for bab in ['BabA','BabU']:
1777	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1778	if hapData['Babinet'][bab][1]:
1779	hapVary.append(pfx+bab)
1780
1781	if Print:
1782	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1783	print >>pFile,135*'-'
1784	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1785	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1786	if hapData['Pref.Ori.'][0] == 'MD':
1787	Ax = hapData['Pref.Ori.'][3]
1788	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1789	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1790	else: #'SH' for spherical harmonics
1791	PrintSHPO(hapData['Pref.Ori.'])
1792	PrintSize(hapData['Size'])
1793	PrintMuStrain(hapData['Mustrain'],SGData)
1794	PrintHStrain(hapData['HStrain'],SGData)
1795	if hapData['Babinet']['BabA'][0]:
1796	PrintBabinet(hapData['Babinet'])
1797	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1798	if resetRefList:
1799	refList = []
1800	Uniq = []
1801	Phi = []
1802	for h,k,l,d in HKLd:
1803	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1804	mul *= 2 # for powder overlap of Friedel pairs
1805	if ext:
1806	continue
1807	if 'C' in inst['Type'][0]:
1808	pos = 2.0asind(wave/(2.0d))+Zero
1809	if limits[0] < pos < limits[1]:
1810	refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,0.0])
1811	Uniq.append(uniq)
1812	Phi.append(phi)
1813	else:
1814	raise ValueError
1815	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{}}
1816	elif 'HKLF' in histogram:
1817	inst = Histogram['Instrument Parameters'][0]
1818	hId = Histogram['hId']
1819	hfx = ':%d:'%(hId)
1820	for item in inst:
1821	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
1822	hapDict[hfx+item] = inst[item][1]
1823	pfx = str(pId)+':'+str(hId)+':'
1824	hapDict[pfx+'Scale'] = hapData['Scale'][0]
1825	if hapData['Scale'][1]:
1826	hapVary.append(pfx+'Scale')
1827
1828	extApprox,extType,extParms = hapData['Extinction']
1829	controlDict[pfx+'EType'] = extType
1830	controlDict[pfx+'EApprox'] = extApprox
1831	controlDict[pfx+'Tbar'] = extParms['Tbar']
1832	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1833	if 'Primary' in extType:
1834	Ekey = ['Ep',]
1835	elif 'I & II' in extType:
1836	Ekey = ['Eg','Es']
1837	elif 'Secondary Type II' == extType:
1838	Ekey = ['Es',]
1839	elif 'Secondary Type I' == extType:
1840	Ekey = ['Eg',]
1841	else: #'None'
1842	Ekey = []
1843	for eKey in Ekey:
1844	hapDict[pfx+eKey] = extParms[eKey][0]
1845	if extParms[eKey][1]:
1846	hapVary.append(pfx+eKey)
1847	for bab in ['BabA','BabU']:
1848	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1849	if hapData['Babinet'][bab][1]:
1850	hapVary.append(pfx+bab)
1851	if Print:
1852	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1853	print >>pFile,135*'-'
1854	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1855	if extType != 'None':
1856	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1857	text = ' Parameters :'
1858	for eKey in Ekey:
1859	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1860	print >>pFile,text
1861	if hapData['Babinet']['BabA'][0]:
1862	PrintBabinet(hapData['Babinet'])
1863	Histogram['Reflection Lists'] = phase
1864
1865	return hapVary,hapDict,controlDict
1866
1867	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1868	'needs a doc string'
1869
1870	def PrintSizeAndSig(hapData,sizeSig):
1871	line = '\n Size model: %9s'%(hapData[0])
1872	refine = False
1873	if hapData[0] in ['isotropic','uniaxial']:
1874	line += ' equatorial:%12.5f'%(hapData[1][0])
1875	if sizeSig[0][0]:
1876	line += ', sig:%8.4f'%(sizeSig[0][0])
1877	refine = True
1878	if hapData[0] == 'uniaxial':
1879	line += ' axial:%12.4f'%(hapData[1][1])
1880	if sizeSig[0][1]:
1881	refine = True
1882	line += ', sig:%8.4f'%(sizeSig[0][1])
1883	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1884	if sizeSig[0][2]:
1885	refine = True
1886	line += ', sig:%8.4f'%(sizeSig[0][2])
1887	if refine:
1888	print >>pFile,line
1889	else:
1890	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1891	if sizeSig[0][2]:
1892	refine = True
1893	line += ', sig:%8.4f'%(sizeSig[0][2])
1894	Snames = ['S11','S22','S33','S12','S13','S23']
1895	ptlbls = ' name :'
1896	ptstr = ' value :'
1897	sigstr = ' sig :'
1898	for i,name in enumerate(Snames):
1899	ptlbls += '%12s' % (name)
1900	ptstr += '%12.6f' % (hapData[4][i])
1901	if sizeSig[1][i]:
1902	refine = True
1903	sigstr += '%12.6f' % (sizeSig[1][i])
1904	else:
1905	sigstr += 12*' '
1906	if refine:
1907	print >>pFile,line
1908	print >>pFile,ptlbls
1909	print >>pFile,ptstr
1910	print >>pFile,sigstr
1911
1912	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1913	line = '\n Mustrain model: %9s'%(hapData[0])
1914	refine = False
1915	if hapData[0] in ['isotropic','uniaxial']:
1916	line += ' equatorial:%12.1f'%(hapData[1][0])
1917	if mustrainSig[0][0]:
1918	line += ', sig:%8.1f'%(mustrainSig[0][0])
1919	refine = True
1920	if hapData[0] == 'uniaxial':
1921	line += ' axial:%12.1f'%(hapData[1][1])
1922	if mustrainSig[0][1]:
1923	line += ', sig:%8.1f'%(mustrainSig[0][1])
1924	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1925	if mustrainSig[0][2]:
1926	refine = True
1927	line += ', sig:%8.3f'%(mustrainSig[0][2])
1928	if refine:
1929	print >>pFile,line
1930	else:
1931	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1932	if mustrainSig[0][2]:
1933	refine = True
1934	line += ', sig:%8.3f'%(mustrainSig[0][2])
1935	Snames = G2spc.MustrainNames(SGData)
1936	ptlbls = ' name :'
1937	ptstr = ' value :'
1938	sigstr = ' sig :'
1939	for i,name in enumerate(Snames):
1940	ptlbls += '%12s' % (name)
1941	ptstr += '%12.6f' % (hapData[4][i])
1942	if mustrainSig[1][i]:
1943	refine = True
1944	sigstr += '%12.6f' % (mustrainSig[1][i])
1945	else:
1946	sigstr += 12*' '
1947	if refine:
1948	print >>pFile,line
1949	print >>pFile,ptlbls
1950	print >>pFile,ptstr
1951	print >>pFile,sigstr
1952
1953	def PrintHStrainAndSig(hapData,strainSig,SGData):
1954	Hsnames = G2spc.HStrainNames(SGData)
1955	ptlbls = ' name :'
1956	ptstr = ' value :'
1957	sigstr = ' sig :'
1958	refine = False
1959	for i,name in enumerate(Hsnames):
1960	ptlbls += '%12s' % (name)
1961	ptstr += '%12.6g' % (hapData[0][i])
1962	if name in strainSig:
1963	refine = True
1964	sigstr += '%12.6g' % (strainSig[name])
1965	else:
1966	sigstr += 12*' '
1967	if refine:
1968	print >>pFile,'\n Hydrostatic/elastic strain: '
1969	print >>pFile,ptlbls
1970	print >>pFile,ptstr
1971	print >>pFile,sigstr
1972
1973	def PrintSHPOAndSig(pfx,hapData,POsig):
1974	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1975	ptlbls = ' names :'
1976	ptstr = ' values:'
1977	sigstr = ' sig :'
1978	for item in hapData[5]:
1979	ptlbls += '%12s'%(item)
1980	ptstr += '%12.3f'%(hapData[5][item])
1981	if pfx+item in POsig:
1982	sigstr += '%12.3f'%(POsig[pfx+item])
1983	else:
1984	sigstr += 12*' '
1985	print >>pFile,ptlbls
1986	print >>pFile,ptstr
1987	print >>pFile,sigstr
1988
1989	def PrintExtAndSig(pfx,hapData,ScalExtSig):
1990	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
1991	text = ''
1992	for item in hapData[2]:
1993	if pfx+item in ScalExtSig:
1994	text += ' %s: '%(item)
1995	text += '%12.2e'%(hapData[2][item][0])
1996	if pfx+item in ScalExtSig:
1997	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
1998	print >>pFile,text
1999
2000	def PrintBabinetAndSig(pfx,hapData,BabSig):
2001	print >>pFile,'\n Babinet form factor modification: '
2002	ptlbls = ' names :'
2003	ptstr = ' values:'
2004	sigstr = ' sig :'
2005	for item in hapData:
2006	ptlbls += '%12s'%(item)
2007	ptstr += '%12.3f'%(hapData[item][0])
2008	if pfx+item in BabSig:
2009	sigstr += '%12.3f'%(BabSig[pfx+item])
2010	else:
2011	sigstr += 12*' '
2012	print >>pFile,ptlbls
2013	print >>pFile,ptstr
2014	print >>pFile,sigstr
2015
2016	PhFrExtPOSig = {}
2017	SizeMuStrSig = {}
2018	ScalExtSig = {}
2019	BabSig = {}
2020	wtFrSum = {}
2021	for phase in Phases:
2022	HistoPhase = Phases[phase]['Histograms']
2023	General = Phases[phase]['General']
2024	SGData = General['SGData']
2025	pId = Phases[phase]['pId']
2026	histoList = HistoPhase.keys()
2027	histoList.sort()
2028	for histogram in histoList:
2029	try:
2030	Histogram = Histograms[histogram]
2031	except KeyError:
2032	#skip if histogram not included e.g. in a sequential refinement
2033	continue
2034	hapData = HistoPhase[histogram]
2035	hId = Histogram['hId']
2036	pfx = str(pId)+':'+str(hId)+':'
2037	if hId not in wtFrSum:
2038	wtFrSum[hId] = 0.
2039	if 'PWDR' in histogram:
2040	for item in ['Scale','Extinction']:
2041	hapData[item][0] = parmDict[pfx+item]
2042	if pfx+item in sigDict:
2043	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2044	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2045	if hapData['Pref.Ori.'][0] == 'MD':
2046	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2047	if pfx+'MD' in sigDict:
2048	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2049	else: #'SH' spherical harmonics
2050	for item in hapData['Pref.Ori.'][5]:
2051	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2052	if pfx+item in sigDict:
2053	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2054	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2055	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2056	pfx+'HStrain':{}})
2057	for item in ['Mustrain','Size']:
2058	hapData[item][1][2] = parmDict[pfx+item+';mx']
2059	hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
2060	if pfx+item+';mx' in sigDict:
2061	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2062	if hapData[item][0] in ['isotropic','uniaxial']:
2063	hapData[item][1][0] = parmDict[pfx+item+';i']
2064	if item == 'Size':
2065	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2066	if pfx+item+';i' in sigDict:
2067	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2068	if hapData[item][0] == 'uniaxial':
2069	hapData[item][1][1] = parmDict[pfx+item+';a']
2070	if item == 'Size':
2071	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2072	if pfx+item+';a' in sigDict:
2073	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2074	else: #generalized for mustrain or ellipsoidal for size
2075	Nterms = len(hapData[item][4])
2076	for i in range(Nterms):
2077	sfx = ':'+str(i)
2078	hapData[item][4][i] = parmDict[pfx+item+sfx]
2079	if pfx+item+sfx in sigDict:
2080	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2081	names = G2spc.HStrainNames(SGData)
2082	for i,name in enumerate(names):
2083	hapData['HStrain'][0][i] = parmDict[pfx+name]
2084	if pfx+name in sigDict:
2085	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2086	for name in ['BabA','BabU']:
2087	hapData['Babinet'][name][0] = parmDict[pfx+name]
2088	if pfx+name in sigDict:
2089	BabSig[pfx+name] = sigDict[pfx+name]
2090
2091	elif 'HKLF' in histogram:
2092	for item in ['Scale',]:
2093	if parmDict.get(pfx+item):
2094	hapData[item][0] = parmDict[pfx+item]
2095	if pfx+item in sigDict:
2096	ScalExtSig[pfx+item] = sigDict[pfx+item]
2097	for item in ['Ep','Eg','Es']:
2098	if parmDict.get(pfx+item):
2099	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2100	if pfx+item in sigDict:
2101	ScalExtSig[pfx+item] = sigDict[pfx+item]
2102	for name in ['BabA','BabU']:
2103	hapData['Babinet'][name][0] = parmDict[pfx+name]
2104	if pfx+name in sigDict:
2105	BabSig[pfx+name] = sigDict[pfx+name]
2106
2107	if Print:
2108	for phase in Phases:
2109	HistoPhase = Phases[phase]['Histograms']
2110	General = Phases[phase]['General']
2111	SGData = General['SGData']
2112	pId = Phases[phase]['pId']
2113	histoList = HistoPhase.keys()
2114	histoList.sort()
2115	for histogram in histoList:
2116	try:
2117	Histogram = Histograms[histogram]
2118	except KeyError:
2119	#skip if histogram not included e.g. in a sequential refinement
2120	continue
2121	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2122	print >>pFile,130*'-'
2123	hapData = HistoPhase[histogram]
2124	hId = Histogram['hId']
2125	Histogram['Residuals'][str(pId)+'::Name'] = phase
2126	pfx = str(pId)+':'+str(hId)+':'
2127	if 'PWDR' in histogram:
2128	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2129	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2130
2131	if pfx+'Scale' in PhFrExtPOSig:
2132	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2133	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2134	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2135	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2136	if pfx+'Extinction' in PhFrExtPOSig:
2137	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2138	if hapData['Pref.Ori.'][0] == 'MD':
2139	if pfx+'MD' in PhFrExtPOSig:
2140	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2141	else:
2142	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2143	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2144	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2145	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2146	if len(BabSig):
2147	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2148
2149	elif 'HKLF' in histogram:
2150	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2151	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2152	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2153	if pfx+'Scale' in ScalExtSig:
2154	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2155	if hapData['Extinction'][0] != 'None':
2156	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2157	if len(BabSig):
2158	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2159
2160	################################################################################
2161	##### Histogram data
2162	################################################################################
2163
2164	def GetHistogramData(Histograms,Print=True,pFile=None):
2165	'needs a doc string'
2166
2167	def GetBackgroundParms(hId,Background):
2168	Back = Background[0]
2169	DebyePeaks = Background[1]
2170	bakType,bakFlag = Back[:2]
2171	backVals = Back[3:]
2172	backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2173	backDict = dict(zip(backNames,backVals))
2174	backVary = []
2175	if bakFlag:
2176	backVary = backNames
2177	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2178	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2179	debyeDict = {}
2180	debyeList = []
2181	for i in range(DebyePeaks['nDebye']):
2182	debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2183	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2184	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2185	debyeVary = []
2186	for item in debyeList:
2187	if item[1]:
2188	debyeVary.append(item[0])
2189	backDict.update(debyeDict)
2190	backVary += debyeVary
2191	peakDict = {}
2192	peakList = []
2193	for i in range(DebyePeaks['nPeaks']):
2194	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2195	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2196	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2197	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2198	peakVary = []
2199	for item in peakList:
2200	if item[1]:
2201	peakVary.append(item[0])
2202	backDict.update(peakDict)
2203	backVary += peakVary
2204	return bakType,backDict,backVary
2205
2206	def GetInstParms(hId,Inst):
2207	dataType = Inst['Type'][0]
2208	instDict = {}
2209	insVary = []
2210	pfx = ':'+str(hId)+':'
2211	insKeys = Inst.keys()
2212	insKeys.sort()
2213	for item in insKeys:
2214	insName = pfx+item
2215	instDict[insName] = Inst[item][1]
2216	if Inst[item][2]:
2217	insVary.append(insName)
2218	# instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2219	# instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2220	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2221	return dataType,instDict,insVary
2222
2223	def GetSampleParms(hId,Sample):
2224	sampVary = []
2225	hfx = ':'+str(hId)+':'
2226	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2227	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2228	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2229	if key in Sample:
2230	sampDict[hfx+key] = Sample[key]
2231	Type = Sample['Type']
2232	if 'Bragg' in Type: #Bragg-Brentano
2233	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2234	sampDict[hfx+item] = Sample[item][0]
2235	if Sample[item][1]:
2236	sampVary.append(hfx+item)
2237	elif 'Debye' in Type: #Debye-Scherrer
2238	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2239	sampDict[hfx+item] = Sample[item][0]
2240	if Sample[item][1]:
2241	sampVary.append(hfx+item)
2242	return Type,sampDict,sampVary
2243
2244	def PrintBackground(Background):
2245	Back = Background[0]
2246	DebyePeaks = Background[1]
2247	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2248	line = ' Coefficients: '
2249	for i,back in enumerate(Back[3:]):
2250	line += '%10.3f'%(back)
2251	if i and not i%10:
2252	line += '\n'+15*' '
2253	print >>pFile,line
2254	if DebyePeaks['nDebye']:
2255	print >>pFile,'\n Debye diffuse scattering coefficients'
2256	parms = ['DebyeA','DebyeR','DebyeU']
2257	line = ' names : '
2258	for parm in parms:
2259	line += '%8s refine?'%(parm)
2260	print >>pFile,line
2261	for j,term in enumerate(DebyePeaks['debyeTerms']):
2262	line = ' term'+'%2d'%(j)+':'
2263	for i in range(3):
2264	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2265	print >>pFile,line
2266	if DebyePeaks['nPeaks']:
2267	print >>pFile,'\n Single peak coefficients'
2268	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2269	line = ' names : '
2270	for parm in parms:
2271	line += '%8s refine?'%(parm)
2272	print >>pFile,line
2273	for j,term in enumerate(DebyePeaks['peaksList']):
2274	line = ' peak'+'%2d'%(j)+':'
2275	for i in range(4):
2276	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2277	print >>pFile,line
2278
2279	def PrintInstParms(Inst):
2280	print >>pFile,'\n Instrument Parameters:'
2281	ptlbls = ' name :'
2282	ptstr = ' value :'
2283	varstr = ' refine:'
2284	insKeys = Inst.keys()
2285	insKeys.sort()
2286	for item in insKeys:
2287	if item != 'Type':
2288	ptlbls += '%12s' % (item)
2289	ptstr += '%12.6f' % (Inst[item][1])
2290	if item in ['Lam1','Lam2','Azimuth']:
2291	varstr += 12*' '
2292	else:
2293	varstr += '%12s' % (str(bool(Inst[item][2])))
2294	print >>pFile,ptlbls
2295	print >>pFile,ptstr
2296	print >>pFile,varstr
2297
2298	def PrintSampleParms(Sample):
2299	print >>pFile,'\n Sample Parameters:'
2300	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2301	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2302	ptlbls = ' name :'
2303	ptstr = ' value :'
2304	varstr = ' refine:'
2305	if 'Bragg' in Sample['Type']:
2306	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2307	ptlbls += '%14s'%(item)
2308	ptstr += '%14.4f'%(Sample[item][0])
2309	varstr += '%14s'%(str(bool(Sample[item][1])))
2310
2311	elif 'Debye' in Type: #Debye-Scherrer
2312	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2313	ptlbls += '%14s'%(item)
2314	ptstr += '%14.4f'%(Sample[item][0])
2315	varstr += '%14s'%(str(bool(Sample[item][1])))
2316
2317	print >>pFile,ptlbls
2318	print >>pFile,ptstr
2319	print >>pFile,varstr
2320
2321	histDict = {}
2322	histVary = []
2323	controlDict = {}
2324	histoList = Histograms.keys()
2325	histoList.sort()
2326	for histogram in histoList:
2327	if 'PWDR' in histogram:
2328	Histogram = Histograms[histogram]
2329	hId = Histogram['hId']
2330	pfx = ':'+str(hId)+':'
2331	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2332	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2333	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2334	for excl in controlDict[pfx+'Exclude']:
2335	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2336	if controlDict[pfx+'Exclude']:
2337	ma.mask_rows(Histogram['Data'])
2338	Background = Histogram['Background']
2339	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2340	controlDict[pfx+'bakType'] = Type
2341	histDict.update(bakDict)
2342	histVary += bakVary
2343
2344	Inst = Histogram['Instrument Parameters'][0]
2345	Type,instDict,insVary = GetInstParms(hId,Inst)
2346	controlDict[pfx+'histType'] = Type
2347	if pfx+'Lam1' in instDict:
2348	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2349	else:
2350	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2351	histDict.update(instDict)
2352	histVary += insVary
2353
2354	Sample = Histogram['Sample Parameters']
2355	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2356	controlDict[pfx+'instType'] = Type
2357	histDict.update(sampDict)
2358	histVary += sampVary
2359
2360
2361	if Print:
2362	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2363	print >>pFile,135*'-'
2364	Units = {'C':' deg','T':' msec'}
2365	units = Units[controlDict[pfx+'histType'][2]]
2366	Limits = controlDict[pfx+'Limits']
2367	print >>pFile,' Instrument type: ',Sample['Type']
2368	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2369	if len(controlDict[pfx+'Exclude']):
2370	excls = controlDict[pfx+'Exclude']
2371	for excl in excls:
2372	print >>pFile,' Excluded region: %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)
2373	PrintSampleParms(Sample)
2374	PrintInstParms(Inst)
2375	PrintBackground(Background)
2376	elif 'HKLF' in histogram:
2377	Histogram = Histograms[histogram]
2378	hId = Histogram['hId']
2379	pfx = ':'+str(hId)+':'
2380	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2381	Inst = Histogram['Instrument Parameters'][0]
2382	controlDict[pfx+'histType'] = Inst['Type'][0]
2383	histDict[pfx+'Lam'] = Inst['Lam'][1]
2384	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
2385	return histVary,histDict,controlDict
2386
2387	def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2388	'needs a doc string'
2389
2390	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2391	Back = Background[0]
2392	DebyePeaks = Background[1]
2393	lenBack = len(Back[3:])
2394	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
2395	for i in range(lenBack):
2396	Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2397	if pfx+'Back:'+str(i) in sigDict:
2398	backSig[i] = sigDict[pfx+'Back:'+str(i)]
2399	if DebyePeaks['nDebye']:
2400	for i in range(DebyePeaks['nDebye']):
2401	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2402	for j,name in enumerate(names):
2403	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2404	if name in sigDict:
2405	backSig[lenBack+3*i+j] = sigDict[name]
2406	if DebyePeaks['nPeaks']:
2407	for i in range(DebyePeaks['nPeaks']):
2408	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2409	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2410	for j,name in enumerate(names):
2411	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2412	if name in sigDict:
2413	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
2414	return backSig
2415
2416	def SetInstParms(pfx,Inst,parmDict,sigDict):
2417	instSig = {}
2418	insKeys = Inst.keys()
2419	insKeys.sort()
2420	for item in insKeys:
2421	insName = pfx+item
2422	Inst[item][1] = parmDict[insName]
2423	if insName in sigDict:
2424	instSig[item] = sigDict[insName]
2425	else:
2426	instSig[item] = 0
2427	return instSig
2428
2429	def SetSampleParms(pfx,Sample,parmDict,sigDict):
2430	if 'Bragg' in Sample['Type']: #Bragg-Brentano
2431	sampSig = [0 for i in range(5)]
2432	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2433	Sample[item][0] = parmDict[pfx+item]
2434	if pfx+item in sigDict:
2435	sampSig[i] = sigDict[pfx+item]
2436	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2437	sampSig = [0 for i in range(4)]
2438	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2439	Sample[item][0] = parmDict[pfx+item]
2440	if pfx+item in sigDict:
2441	sampSig[i] = sigDict[pfx+item]
2442	return sampSig
2443
2444	def PrintBackgroundSig(Background,backSig):
2445	Back = Background[0]
2446	DebyePeaks = Background[1]
2447	lenBack = len(Back[3:])
2448	valstr = ' value : '
2449	sigstr = ' sig : '
2450	refine = False
2451	for i,back in enumerate(Back[3:]):
2452	valstr += '%10.4g'%(back)
2453	if Back[1]:
2454	refine = True
2455	sigstr += '%10.4g'%(backSig[i])
2456	else:
2457	sigstr += 10*' '
2458	if refine:
2459	print >>pFile,'\n Background function: ',Back[0]
2460	print >>pFile,valstr
2461	print >>pFile,sigstr
2462	if DebyePeaks['nDebye']:
2463	ifAny = False
2464	ptfmt = "%12.3f"
2465	names = ' names :'
2466	ptstr = ' values:'
2467	sigstr = ' esds :'
2468	for item in sigDict:
2469	if 'Debye' in item:
2470	ifAny = True
2471	names += '%12s'%(item)
2472	ptstr += ptfmt%(parmDict[item])
2473	sigstr += ptfmt%(sigDict[item])
2474	if ifAny:
2475	print >>pFile,'\n Debye diffuse scattering coefficients'
2476	print >>pFile,names
2477	print >>pFile,ptstr
2478	print >>pFile,sigstr
2479	if DebyePeaks['nPeaks']:
2480	ifAny = False
2481	ptfmt = "%14.3f"
2482	names = ' names :'
2483	ptstr = ' values:'
2484	sigstr = ' esds :'
2485	for item in sigDict:
2486	if 'BkPk' in item:
2487	ifAny = True
2488	names += '%14s'%(item)
2489	ptstr += ptfmt%(parmDict[item])
2490	sigstr += ptfmt%(sigDict[item])
2491	if ifAny:
2492	print >>pFile,'\n Single peak coefficients'
2493	print >>pFile,names
2494	print >>pFile,ptstr
2495	print >>pFile,sigstr
2496
2497	def PrintInstParmsSig(Inst,instSig):
2498	ptlbls = ' names :'
2499	ptstr = ' value :'
2500	sigstr = ' sig :'
2501	refine = False
2502	insKeys = instSig.keys()
2503	insKeys.sort()
2504	for name in insKeys:
2505	if name not in ['Type','Lam1','Lam2','Azimuth']:
2506	ptlbls += '%12s' % (name)
2507	ptstr += '%12.6f' % (Inst[name][1])
2508	if instSig[name]:
2509	refine = True
2510	sigstr += '%12.6f' % (instSig[name])
2511	else:
2512	sigstr += 12*' '
2513	if refine:
2514	print >>pFile,'\n Instrument Parameters:'
2515	print >>pFile,ptlbls
2516	print >>pFile,ptstr
2517	print >>pFile,sigstr
2518
2519	def PrintSampleParmsSig(Sample,sampleSig):
2520	ptlbls = ' names :'
2521	ptstr = ' values:'
2522	sigstr = ' sig :'
2523	refine = False
2524	if 'Bragg' in Sample['Type']:
2525	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2526	ptlbls += '%14s'%(item)
2527	ptstr += '%14.4f'%(Sample[item][0])
2528	if sampleSig[i]:
2529	refine = True
2530	sigstr += '%14.4f'%(sampleSig[i])
2531	else:
2532	sigstr += 14*' '
2533
2534	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2535	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2536	ptlbls += '%14s'%(item)
2537	ptstr += '%14.4f'%(Sample[item][0])
2538	if sampleSig[i]:
2539	refine = True
2540	sigstr += '%14.4f'%(sampleSig[i])
2541	else:
2542	sigstr += 14*' '
2543
2544	if refine:
2545	print >>pFile,'\n Sample Parameters:'
2546	print >>pFile,ptlbls
2547	print >>pFile,ptstr
2548	print >>pFile,sigstr
2549
2550	histoList = Histograms.keys()
2551	histoList.sort()
2552	for histogram in histoList:
2553	if 'PWDR' in histogram:
2554	Histogram = Histograms[histogram]
2555	hId = Histogram['hId']
2556	pfx = ':'+str(hId)+':'
2557	Background = Histogram['Background']
2558	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2559
2560	Inst = Histogram['Instrument Parameters'][0]
2561	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2562
2563	Sample = Histogram['Sample Parameters']
2564	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2565
2566	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2567	print >>pFile,135*'-'
2568	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2569	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
2570	print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
2571	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
2572	if Print:
2573	print >>pFile,' Instrument type: ',Sample['Type']
2574	PrintSampleParmsSig(Sample,sampSig)
2575	PrintInstParmsSig(Inst,instSig)
2576	PrintBackgroundSig(Background,backSig)
2577

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

Download in other formats: