Context Navigation

source: trunk/GSASIIstrIO.py @ 1409

Last change on this file since 1409 was 1409, checked in by toby, 9 years ago
more problems with older GPX files
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 111.3 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2014-07-03 22:19:01 +0000 (Thu, 03 Jul 2014) $
4	# $Author: toby $
5	# $Revision: 1409 $
6	# $URL: trunk/GSASIIstrIO.py $
7	# $Id: GSASIIstrIO.py 1409 2014-07-03 22:19:01Z 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: 1409 $")
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]['Mass']
751	if 'SL' in bldata[1]:
752	blen = bldata[1]['SL'][0]
753	bres = []
754	else:
755	blen = 0
756	bres = bldata[1]['BW-LS']
757	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
758	for item in bres:
759	line += '%10.5g'%(item)
760	print >>pFile,line
761
762	def PrintRBObjects(resRBData,vecRBData):
763
764	def PrintRBThermals():
765	tlstr = ['11','22','33','12','13','23']
766	sstr = ['12','13','21','23','31','32','AA','BB']
767	TLS = RB['ThermalMotion'][1]
768	TLSvar = RB['ThermalMotion'][2]
769	if 'T' in RB['ThermalMotion'][0]:
770	print >>pFile,'TLS data'
771	text = ''
772	for i in range(6):
773	text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
774	print >>pFile,text
775	if 'L' in RB['ThermalMotion'][0]:
776	text = ''
777	for i in range(6,12):
778	text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
779	print >>pFile,text
780	if 'S' in RB['ThermalMotion'][0]:
781	text = ''
782	for i in range(12,20):
783	text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
784	print >>pFile,text
785	if 'U' in RB['ThermalMotion'][0]:
786	print >>pFile,'Uiso data'
787	text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])
788
789	if len(resRBData):
790	for RB in resRBData:
791	Oxyz = RB['Orig'][0]
792	Qrijk = RB['Orient'][0]
793	Angle = 2.0*acosd(Qrijk[0])
794	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
795	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
796	print >>pFile,'Orientation angle,vector:', \
797	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
798	Torsions = RB['Torsions']
799	if len(Torsions):
800	text = 'Torsions: '
801	for torsion in Torsions:
802	text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
803	print >>pFile,text
804	PrintRBThermals()
805	if len(vecRBData):
806	for RB in vecRBData:
807	Oxyz = RB['Orig'][0]
808	Qrijk = RB['Orient'][0]
809	Angle = 2.0*acosd(Qrijk[0])
810	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
811	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
812	print >>pFile,'Orientation angle,vector:', \
813	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
814	PrintRBThermals()
815
816	def PrintAtoms(General,Atoms):
817	cx,ct,cs,cia = General['AtomPtrs']
818	print >>pFile,'\n Atoms:'
819	line = ' name type refine? x y z '+ \
820	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
821	if General['Type'] == 'magnetic':
822	line += ' Mx My Mz'
823	elif General['Type'] == 'macromolecular':
824	line = ' res no residue chain'+line
825	print >>pFile,line
826	if General['Type'] == 'nuclear':
827	print >>pFile,135*'-'
828	for i,at in enumerate(Atoms):
829	line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
830	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
831	if at[cia] == 'I':
832	line += '%8.4f'%(at[cia+1])+48*' '
833	else:
834	line += 8*' '
835	for j in range(6):
836	line += '%8.4f'%(at[cia+1+j])
837	print >>pFile,line
838	elif General['Type'] == 'macromolecular':
839	print >>pFile,135*'-'
840	for i,at in enumerate(Atoms):
841	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])+ \
842	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
843	if at[cia] == 'I':
844	line += '%8.4f'%(at[cia+1])+48*' '
845	else:
846	line += 8*' '
847	for j in range(6):
848	line += '%8.4f'%(at[cia+1+j])
849	print >>pFile,line
850
851	def PrintTexture(textureData):
852	topstr = '\n Spherical harmonics texture: Order:' + \
853	str(textureData['Order'])
854	if textureData['Order']:
855	print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
856	else:
857	print >>pFile,topstr
858	return
859	names = ['omega','chi','phi']
860	line = '\n'
861	for name in names:
862	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
863	print >>pFile,line
864	print >>pFile,'\n Texture coefficients:'
865	ptlbls = ' names :'
866	ptstr = ' values:'
867	SHcoeff = textureData['SH Coeff'][1]
868	for item in SHcoeff:
869	ptlbls += '%12s'%(item)
870	ptstr += '%12.4f'%(SHcoeff[item])
871	print >>pFile,ptlbls
872	print >>pFile,ptstr
873
874	def MakeRBParms(rbKey,phaseVary,phaseDict):
875	rbid = str(rbids.index(RB['RBId']))
876	pfxRB = pfx+'RB'+rbKey+'P'
877	pstr = ['x','y','z']
878	ostr = ['a','i','j','k']
879	for i in range(3):
880	name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
881	phaseDict[name] = RB['Orig'][0][i]
882	if RB['Orig'][1]:
883	phaseVary += [name,]
884	pfxRB = pfx+'RB'+rbKey+'O'
885	for i in range(4):
886	name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
887	phaseDict[name] = RB['Orient'][0][i]
888	if RB['Orient'][1] == 'AV' and i:
889	phaseVary += [name,]
890	elif RB['Orient'][1] == 'A' and not i:
891	phaseVary += [name,]
892
893	def MakeRBThermals(rbKey,phaseVary,phaseDict):
894	rbid = str(rbids.index(RB['RBId']))
895	tlstr = ['11','22','33','12','13','23']
896	sstr = ['12','13','21','23','31','32','AA','BB']
897	if 'T' in RB['ThermalMotion'][0]:
898	pfxRB = pfx+'RB'+rbKey+'T'
899	for i in range(6):
900	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
901	phaseDict[name] = RB['ThermalMotion'][1][i]
902	if RB['ThermalMotion'][2][i]:
903	phaseVary += [name,]
904	if 'L' in RB['ThermalMotion'][0]:
905	pfxRB = pfx+'RB'+rbKey+'L'
906	for i in range(6):
907	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
908	phaseDict[name] = RB['ThermalMotion'][1][i+6]
909	if RB['ThermalMotion'][2][i+6]:
910	phaseVary += [name,]
911	if 'S' in RB['ThermalMotion'][0]:
912	pfxRB = pfx+'RB'+rbKey+'S'
913	for i in range(8):
914	name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
915	phaseDict[name] = RB['ThermalMotion'][1][i+12]
916	if RB['ThermalMotion'][2][i+12]:
917	phaseVary += [name,]
918	if 'U' in RB['ThermalMotion'][0]:
919	name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
920	phaseDict[name] = RB['ThermalMotion'][1][0]
921	if RB['ThermalMotion'][2][0]:
922	phaseVary += [name,]
923
924	def MakeRBTorsions(rbKey,phaseVary,phaseDict):
925	rbid = str(rbids.index(RB['RBId']))
926	pfxRB = pfx+'RB'+rbKey+'Tr;'
927	for i,torsion in enumerate(RB['Torsions']):
928	name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
929	phaseDict[name] = torsion[0]
930	if torsion[1]:
931	phaseVary += [name,]
932
933	if Print:
934	print >>pFile,'\n Phases:'
935	phaseVary = []
936	phaseDict = {}
937	phaseConstr = {}
938	pawleyLookup = {}
939	FFtables = {} #scattering factors - xrays
940	BLtables = {} # neutrons
941	Natoms = {}
942	AtMults = {}
943	AtIA = {}
944	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
945	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
946	atomIndx = {}
947	for name in PhaseData:
948	General = PhaseData[name]['General']
949	pId = PhaseData[name]['pId']
950	pfx = str(pId)+'::'
951	FFtable = G2el.GetFFtable(General['AtomTypes'])
952	BLtable = G2el.GetBLtable(General)
953	FFtables.update(FFtable)
954	BLtables.update(BLtable)
955	Atoms = PhaseData[name]['Atoms']
956	AtLookup = G2mth.FillAtomLookUp(Atoms)
957	PawleyRef = PhaseData[name].get('Pawley ref',[])
958	SGData = General['SGData']
959	SGtext = G2spc.SGPrint(SGData)
960	cell = General['Cell']
961	A = G2lat.cell2A(cell[1:7])
962	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
963	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
964	if cell[0]:
965	phaseVary += cellVary(pfx,SGData)
966	resRBData = PhaseData[name]['RBModels'].get('Residue',[])
967	if resRBData:
968	rbids = rbIds['Residue'] #NB: used in the MakeRB routines
969	for iRB,RB in enumerate(resRBData):
970	MakeRBParms('R',phaseVary,phaseDict)
971	MakeRBThermals('R',phaseVary,phaseDict)
972	MakeRBTorsions('R',phaseVary,phaseDict)
973
974	vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
975	if vecRBData:
976	rbids = rbIds['Vector'] #NB: used in the MakeRB routines
977	for iRB,RB in enumerate(vecRBData):
978	MakeRBParms('V',phaseVary,phaseDict)
979	MakeRBThermals('V',phaseVary,phaseDict)
980
981	Natoms[pfx] = 0
982	if Atoms and not General.get('doPawley'):
983	cx,ct,cs,cia = General['AtomPtrs']
984	if General['Type'] in ['nuclear','macromolecular']:
985	Natoms[pfx] = len(Atoms)
986	for i,at in enumerate(Atoms):
987	atomIndx[at[-1]] = [pfx,i] #lookup table for restraints
988	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
989	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
990	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
991	pfx+'AI/A:'+str(i):at[cia],})
992	if at[cia] == 'I':
993	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
994	else:
995	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
996	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
997	if 'F' in at[ct+1]:
998	phaseVary.append(pfx+'Afrac:'+str(i))
999	if 'X' in at[ct+1]:
1000	xId,xCoef = G2spc.GetCSxinel(at[cs])
1001	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1002	equivs = [[],[],[]]
1003	for j in range(3):
1004	if xId[j] > 0:
1005	phaseVary.append(names[j])
1006	equivs[xId[j]-1].append([names[j],xCoef[j]])
1007	for equiv in equivs:
1008	if len(equiv) > 1:
1009	name = equiv[0][0]
1010	for eqv in equiv[1:]:
1011	G2mv.StoreEquivalence(name,(eqv,))
1012	if 'U' in at[ct+1]:
1013	if at[cia] == 'I':
1014	phaseVary.append(pfx+'AUiso:'+str(i))
1015	else:
1016	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1017	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1018	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1019	equivs = [[],[],[],[],[],[]]
1020	for j in range(6):
1021	if uId[j] > 0:
1022	phaseVary.append(names[j])
1023	equivs[uId[j]-1].append([names[j],uCoef[j]])
1024	for equiv in equivs:
1025	if len(equiv) > 1:
1026	name = equiv[0][0]
1027	for eqv in equiv[1:]:
1028	G2mv.StoreEquivalence(name,(eqv,))
1029	# elif General['Type'] == 'magnetic':
1030	# elif General['Type'] == 'macromolecular':
1031	textureData = General['SH Texture']
1032	if textureData['Order']:
1033	phaseDict[pfx+'SHorder'] = textureData['Order']
1034	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1035	for item in ['omega','chi','phi']:
1036	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1037	if textureData['Sample '+item][0]:
1038	phaseVary.append(pfx+'SH '+item)
1039	for item in textureData['SH Coeff'][1]:
1040	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1041	if textureData['SH Coeff'][0]:
1042	phaseVary.append(pfx+item)
1043
1044	if Print:
1045	print >>pFile,'\n Phase name: ',General['Name']
1046	print >>pFile,135*'-'
1047	PrintFFtable(FFtable)
1048	PrintBLtable(BLtable)
1049	print >>pFile,''
1050	for line in SGtext: print >>pFile,line
1051	PrintRBObjects(resRBData,vecRBData)
1052	PrintAtoms(General,Atoms)
1053	print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
1054	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
1055	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
1056	PrintTexture(textureData)
1057	if name in RestraintDict:
1058	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1059	textureData,RestraintDict[name],pFile)
1060
1061	elif PawleyRef:
1062	pawleyVary = []
1063	for i,refl in enumerate(PawleyRef):
1064	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
1065	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1066	if refl[5]:
1067	pawleyVary.append(pfx+'PWLref:'+str(i))
1068	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
1069	phaseVary += pawleyVary
1070
1071	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1072
1073	def cellFill(pfx,SGData,parmDict,sigDict):
1074	'''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1075	from the parameter and sig dictionaries.
1076
1077	:param str pfx: parameter prefix ("n::", where n is a phase number)
1078	:param dict SGdata: a symmetry object
1079	:param dict parmDict: a dictionary of parameters
1080	:param dict sigDict: a dictionary of uncertainties on parameters
1081
1082	:returns: A,sigA where each is a list of six terms with the A terms
1083	'''
1084	if SGData['SGLaue'] in ['-1',]:
1085	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1086	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1087	elif SGData['SGLaue'] in ['2/m',]:
1088	if SGData['SGUniq'] == 'a':
1089	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1090	parmDict[pfx+'A3'],0,0]
1091	elif SGData['SGUniq'] == 'b':
1092	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1093	0,parmDict[pfx+'A4'],0]
1094	else:
1095	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1096	0,0,parmDict[pfx+'A5']]
1097	elif SGData['SGLaue'] in ['mmm',]:
1098	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1099	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1100	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1101	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1102	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1103	parmDict[pfx+'A0'],0,0]
1104	elif SGData['SGLaue'] in ['3R', '3mR']:
1105	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1106	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1107	elif SGData['SGLaue'] in ['m3m','m3']:
1108	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1109
1110	try:
1111	if SGData['SGLaue'] in ['-1',]:
1112	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1113	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1114	elif SGData['SGLaue'] in ['2/m',]:
1115	if SGData['SGUniq'] == 'a':
1116	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1117	sigDict[pfx+'A3'],0,0]
1118	elif SGData['SGUniq'] == 'b':
1119	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1120	0,sigDict[pfx+'A4'],0]
1121	else:
1122	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1123	0,0,sigDict[pfx+'A5']]
1124	elif SGData['SGLaue'] in ['mmm',]:
1125	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1126	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1127	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1128	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1129	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1130	elif SGData['SGLaue'] in ['3R', '3mR']:
1131	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1132	elif SGData['SGLaue'] in ['m3m','m3']:
1133	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1134	except KeyError:
1135	sigA = [0,0,0,0,0,0]
1136
1137	return A,sigA
1138
1139	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1140	'needs a doc string'
1141	if phaseRest:
1142	Amat = G2lat.cell2AB(cell)[0]
1143	cx,ct,cs = AtPtrs[:3]
1144	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1145	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1146	['ChemComp','Sites'],['Texture','HKLs']]
1147	for name,rest in names:
1148	itemRest = phaseRest[name]
1149	if itemRest[rest] and itemRest['Use']:
1150	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1151	if name in ['Bond','Angle','Plane','Chiral']:
1152	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
1153	for indx,ops,obs,esd in itemRest[rest]:
1154	try:
1155	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1156	AtName = ''
1157	for i,Aname in enumerate(AtNames):
1158	AtName += Aname
1159	if ops[i] == '1':
1160	AtName += '-'
1161	else:
1162	AtName += '+('+ops[i]+')-'
1163	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1164	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1165	if name == 'Bond':
1166	calc = G2mth.getRestDist(XYZ,Amat)
1167	elif name == 'Angle':
1168	calc = G2mth.getRestAngle(XYZ,Amat)
1169	elif name == 'Plane':
1170	calc = G2mth.getRestPlane(XYZ,Amat)
1171	elif name == 'Chiral':
1172	calc = G2mth.getRestChiral(XYZ,Amat)
1173	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1174	except KeyError:
1175	del itemRest[rest]
1176	elif name in ['Torsion','Rama']:
1177	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
1178	coeffDict = itemRest['Coeff']
1179	for indx,ops,cofName,esd in itemRest[rest]:
1180	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1181	AtName = ''
1182	for i,Aname in enumerate(AtNames):
1183	AtName += Aname+'+('+ops[i]+')-'
1184	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1185	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1186	if name == 'Torsion':
1187	tor = G2mth.getRestTorsion(XYZ,Amat)
1188	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1189	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1190	else:
1191	phi,psi = G2mth.getRestRama(XYZ,Amat)
1192	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1193	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])
1194	elif name == 'ChemComp':
1195	print >>pFile,' atoms mul*frac factor prod'
1196	for indx,factors,obs,esd in itemRest[rest]:
1197	try:
1198	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1199	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1200	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1201	mulfrac = mul*frac
1202	calcs = mulfracfactors
1203	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1204	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1205	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1206	except KeyError:
1207	del itemRest[rest]
1208	elif name == 'Texture' and textureData['Order']:
1209	Start = False
1210	SHCoef = textureData['SH Coeff'][1]
1211	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1212	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1213	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
1214	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1215	PH = np.array(hkl)
1216	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1217	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1218	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1219	Z = ma.masked_greater(Z,0.0)
1220	num = ma.count(Z)
1221	sum = 0
1222	if num:
1223	sum = np.sum(Z)
1224	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)
1225
1226	def getCellEsd(pfx,SGData,A,covData):
1227	'needs a doc string'
1228	dpr = 180./np.pi
1229	rVsq = G2lat.calc_rVsq(A)
1230	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1231	cell = np.array(G2lat.Gmat2cell(g)) #real cell
1232	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1233	scos = cosd(cellst[3:6])
1234	ssin = sind(cellst[3:6])
1235	scot = scos/ssin
1236	rcos = cosd(cell[3:6])
1237	rsin = sind(cell[3:6])
1238	rcot = rcos/rsin
1239	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1240	varyList = covData['varyList']
1241	covMatrix = covData['covMatrix']
1242	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1243	Ax = np.array(A)
1244	Ax[3:] /= 2.
1245	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
1246	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
1247	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1248	Vol = 1/np.sqrt(rVsq)
1249	sigVol = Vol*3np.sqrt(srcvlsq)/2.
1250	R123 = Ax[0]Ax[1]Ax[2]
1251	dsasdg = np.zeros((3,6))
1252	dadg = np.zeros((6,6))
1253	for i0 in range(3): #0 1 2
1254	i1 = (i0+1)%3 #1 2 0
1255	i2 = (i1+1)%3 #2 0 1
1256	i3 = 5-i2 #3 5 4
1257	i4 = 5-i1 #4 3 5
1258	i5 = 5-i0 #5 4 3
1259	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
1260	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
1261	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1262	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
1263	denom = np.sqrt(denmsq)
1264	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
1265	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
1266	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
1267	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
1268	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
1269	dadg[i5][i5] = -Ax[i0]/denom
1270	for i0 in range(3):
1271	i1 = (i0+1)%3
1272	i2 = (i1+1)%3
1273	i3 = 5-i2
1274	for ij in range(6):
1275	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1276	if ij == i0:
1277	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1278	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
1279	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1280	var = np.diag(sigMat)
1281	CS = np.where(var>0.,np.sqrt(var),0.)
1282	return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
1283
1284	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1285	'needs a doc string'
1286
1287	def PrintAtomsAndSig(General,Atoms,atomsSig):
1288	print >>pFile,'\n Atoms:'
1289	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1290	if General['Type'] == 'magnetic':
1291	line += ' Mx My Mz'
1292	elif General['Type'] == 'macromolecular':
1293	line = ' res no residue chain '+line
1294	print >>pFile,line
1295	if General['Type'] == 'nuclear':
1296	print >>pFile,135*'-'
1297	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1298	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1299	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
1300	for atyp in General['AtomTypes']: #zero composition
1301	General['NoAtoms'][atyp] = 0.
1302	for i,at in enumerate(Atoms):
1303	General['NoAtoms'][at[1]] += at[6]*float(at[8]) #new composition
1304	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1305	valstr = ' values:'
1306	sigstr = ' sig :'
1307	for ind in [3,4,5,6]:
1308	sigind = str(i)+':'+str(ind)
1309	valstr += fmt[ind]%(at[ind])
1310	if sigind in atomsSig:
1311	sigstr += fmt[ind]%(atomsSig[sigind])
1312	else:
1313	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1314	if at[9] == 'I':
1315	valstr += fmt[10]%(at[10])
1316	if str(i)+':10' in atomsSig:
1317	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1318	else:
1319	sigstr += 8*' '
1320	else:
1321	valstr += 8*' '
1322	sigstr += 8*' '
1323	for ind in [11,12,13,14,15,16]:
1324	sigind = str(i)+':'+str(ind)
1325	valstr += fmt[ind]%(at[ind])
1326	if sigind in atomsSig:
1327	sigstr += fmt[ind]%(atomsSig[sigind])
1328	else:
1329	sigstr += 8*' '
1330	print >>pFile,name
1331	print >>pFile,valstr
1332	print >>pFile,sigstr
1333
1334	def PrintRBObjPOAndSig(rbfx,rbsx):
1335	namstr = ' names :'
1336	valstr = ' values:'
1337	sigstr = ' esds :'
1338	for i,px in enumerate(['Px:','Py:','Pz:']):
1339	name = pfx+rbfx+px+rbsx
1340	namstr += '%12s'%('Pos '+px[1])
1341	valstr += '%12.5f'%(parmDict[name])
1342	if name in sigDict:
1343	sigstr += '%12.5f'%(sigDict[name])
1344	else:
1345	sigstr += 12*' '
1346	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1347	name = pfx+rbfx+po+rbsx
1348	namstr += '%12s'%('Ori '+po[1])
1349	valstr += '%12.5f'%(parmDict[name])
1350	if name in sigDict:
1351	sigstr += '%12.5f'%(sigDict[name])
1352	else:
1353	sigstr += 12*' '
1354	print >>pFile,namstr
1355	print >>pFile,valstr
1356	print >>pFile,sigstr
1357
1358	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1359	namstr = ' names :'
1360	valstr = ' values:'
1361	sigstr = ' esds :'
1362	if 'N' not in TLS:
1363	print >>pFile,' Thermal motion:'
1364	if 'T' in TLS:
1365	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1366	name = pfx+rbfx+pt+rbsx
1367	namstr += '%12s'%(pt[:3])
1368	valstr += '%12.5f'%(parmDict[name])
1369	if name in sigDict:
1370	sigstr += '%12.5f'%(sigDict[name])
1371	else:
1372	sigstr += 12*' '
1373	print >>pFile,namstr
1374	print >>pFile,valstr
1375	print >>pFile,sigstr
1376	if 'L' in TLS:
1377	namstr = ' names :'
1378	valstr = ' values:'
1379	sigstr = ' esds :'
1380	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1381	name = pfx+rbfx+pt+rbsx
1382	namstr += '%12s'%(pt[:3])
1383	valstr += '%12.3f'%(parmDict[name])
1384	if name in sigDict:
1385	sigstr += '%12.3f'%(sigDict[name])
1386	else:
1387	sigstr += 12*' '
1388	print >>pFile,namstr
1389	print >>pFile,valstr
1390	print >>pFile,sigstr
1391	if 'S' in TLS:
1392	namstr = ' names :'
1393	valstr = ' values:'
1394	sigstr = ' esds :'
1395	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1396	name = pfx+rbfx+pt+rbsx
1397	namstr += '%12s'%(pt[:3])
1398	valstr += '%12.4f'%(parmDict[name])
1399	if name in sigDict:
1400	sigstr += '%12.4f'%(sigDict[name])
1401	else:
1402	sigstr += 12*' '
1403	print >>pFile,namstr
1404	print >>pFile,valstr
1405	print >>pFile,sigstr
1406	if 'U' in TLS:
1407	name = pfx+rbfx+'U:'+rbsx
1408	namstr = ' names :'+'%12s'%('Uiso')
1409	valstr = ' values:'+'%12.5f'%(parmDict[name])
1410	if name in sigDict:
1411	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1412	else:
1413	sigstr = ' esds :'+12*' '
1414	print >>pFile,namstr
1415	print >>pFile,valstr
1416	print >>pFile,sigstr
1417
1418	def PrintRBObjTorAndSig(rbsx):
1419	namstr = ' names :'
1420	valstr = ' values:'
1421	sigstr = ' esds :'
1422	nTors = len(RBObj['Torsions'])
1423	if nTors:
1424	print >>pFile,' Torsions:'
1425	for it in range(nTors):
1426	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1427	namstr += '%12s'%('Tor'+str(it))
1428	valstr += '%12.4f'%(parmDict[name])
1429	if name in sigDict:
1430	sigstr += '%12.4f'%(sigDict[name])
1431	print >>pFile,namstr
1432	print >>pFile,valstr
1433	print >>pFile,sigstr
1434
1435	def PrintSHtextureAndSig(textureData,SHtextureSig):
1436	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1437	names = ['omega','chi','phi']
1438	namstr = ' names :'
1439	ptstr = ' values:'
1440	sigstr = ' esds :'
1441	for name in names:
1442	namstr += '%12s'%(name)
1443	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1444	if 'Sample '+name in SHtextureSig:
1445	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1446	else:
1447	sigstr += 12*' '
1448	print >>pFile,namstr
1449	print >>pFile,ptstr
1450	print >>pFile,sigstr
1451	print >>pFile,'\n Texture coefficients:'
1452	namstr = ' names :'
1453	ptstr = ' values:'
1454	sigstr = ' esds :'
1455	SHcoeff = textureData['SH Coeff'][1]
1456	for name in SHcoeff:
1457	namstr += '%12s'%(name)
1458	ptstr += '%12.3f'%(SHcoeff[name])
1459	if name in SHtextureSig:
1460	sigstr += '%12.3f'%(SHtextureSig[name])
1461	else:
1462	sigstr += 12*' '
1463	print >>pFile,namstr
1464	print >>pFile,ptstr
1465	print >>pFile,sigstr
1466
1467	print >>pFile,'\n Phases:'
1468	for phase in Phases:
1469	print >>pFile,' Result for phase: ',phase
1470	Phase = Phases[phase]
1471	General = Phase['General']
1472	SGData = General['SGData']
1473	Atoms = Phase['Atoms']
1474	AtLookup = G2mth.FillAtomLookUp(Atoms)
1475	cell = General['Cell']
1476	pId = Phase['pId']
1477	pfx = str(pId)+'::'
1478	if cell[0]:
1479	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1480	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1481	print >>pFile,' Reciprocal metric tensor: '
1482	ptfmt = "%15.9f"
1483	names = ['A11','A22','A33','A12','A13','A23']
1484	namstr = ' names :'
1485	ptstr = ' values:'
1486	sigstr = ' esds :'
1487	for name,a,siga in zip(names,A,sigA):
1488	namstr += '%15s'%(name)
1489	ptstr += ptfmt%(a)
1490	if siga:
1491	sigstr += ptfmt%(siga)
1492	else:
1493	sigstr += 15*' '
1494	print >>pFile,namstr
1495	print >>pFile,ptstr
1496	print >>pFile,sigstr
1497	cell[1:7] = G2lat.A2cell(A)
1498	cell[7] = G2lat.calc_V(A)
1499	print >>pFile,' New unit cell:'
1500	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1501	names = ['a','b','c','alpha','beta','gamma','Volume']
1502	namstr = ' names :'
1503	ptstr = ' values:'
1504	sigstr = ' esds :'
1505	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1506	namstr += '%12s'%(name)
1507	ptstr += fmt%(a)
1508	if siga:
1509	sigstr += fmt%(siga)
1510	else:
1511	sigstr += 12*' '
1512	print >>pFile,namstr
1513	print >>pFile,ptstr
1514	print >>pFile,sigstr
1515
1516	General['Mass'] = 0.
1517	if Phase['General'].get('doPawley'):
1518	pawleyRef = Phase['Pawley ref']
1519	for i,refl in enumerate(pawleyRef):
1520	key = pfx+'PWLref:'+str(i)
1521	refl[6] = parmDict[key]
1522	if key in sigDict:
1523	refl[7] = sigDict[key]
1524	else:
1525	refl[7] = 0
1526	else:
1527	VRBIds = RBIds['Vector']
1528	RRBIds = RBIds['Residue']
1529	RBModels = Phase['RBModels']
1530	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1531	jrb = VRBIds.index(RBObj['RBId'])
1532	rbsx = str(irb)+':'+str(jrb)
1533	print >>pFile,' Vector rigid body parameters:'
1534	PrintRBObjPOAndSig('RBV',rbsx)
1535	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1536	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1537	jrb = RRBIds.index(RBObj['RBId'])
1538	rbsx = str(irb)+':'+str(jrb)
1539	print >>pFile,' Residue rigid body parameters:'
1540	PrintRBObjPOAndSig('RBR',rbsx)
1541	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1542	PrintRBObjTorAndSig(rbsx)
1543	atomsSig = {}
1544	if General['Type'] == 'nuclear': #this needs macromolecular variant!
1545	for i,at in enumerate(Atoms):
1546	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1547	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1548	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1549	for ind in [3,4,5,6]:
1550	at[ind] = parmDict[names[ind]]
1551	if ind in [3,4,5]:
1552	name = names[ind].replace('A','dA')
1553	else:
1554	name = names[ind]
1555	if name in sigDict:
1556	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1557	if at[9] == 'I':
1558	at[10] = parmDict[names[10]]
1559	if names[10] in sigDict:
1560	atomsSig[str(i)+':10'] = sigDict[names[10]]
1561	else:
1562	for ind in [11,12,13,14,15,16]:
1563	at[ind] = parmDict[names[ind]]
1564	if names[ind] in sigDict:
1565	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1566	ind = General['AtomTypes'].index(at[1])
1567	General['Mass'] += General['AtomMass'][ind]at[6]at[8]
1568	PrintAtomsAndSig(General,Atoms,atomsSig)
1569
1570	textureData = General['SH Texture']
1571	if textureData['Order']:
1572	SHtextureSig = {}
1573	for name in ['omega','chi','phi']:
1574	aname = pfx+'SH '+name
1575	textureData['Sample '+name][1] = parmDict[aname]
1576	if aname in sigDict:
1577	SHtextureSig['Sample '+name] = sigDict[aname]
1578	for name in textureData['SH Coeff'][1]:
1579	aname = pfx+name
1580	textureData['SH Coeff'][1][name] = parmDict[aname]
1581	if aname in sigDict:
1582	SHtextureSig[name] = sigDict[aname]
1583	PrintSHtextureAndSig(textureData,SHtextureSig)
1584	if phase in RestraintDict:
1585	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1586	textureData,RestraintDict[phase],pFile)
1587
1588	################################################################################
1589	##### Histogram & Phase data
1590	################################################################################
1591
1592	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1593	'''Loads the HAP histogram/phase information into dicts
1594
1595	:param dict Phases: phase information
1596	:param dict Histograms: Histogram information
1597	:param bool Print: prints information as it is read
1598	:param file pFile: file object to print to (the default, None causes printing to the console)
1599	:param bool resetRefList: Should the contents of the Reflection List be initialized
1600	on loading. The default, True, initializes the Reflection List as it is loaded.
1601
1602	:returns: (hapVary,hapDict,controlDict)
1603	* hapVary: list of refined variables
1604	* hapDict: dict with refined variables and their values
1605	* controlDict: dict with computation controls (?)
1606	'''
1607
1608	def PrintSize(hapData):
1609	if hapData[0] in ['isotropic','uniaxial']:
1610	line = '\n Size model : %9s'%(hapData[0])
1611	line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1612	if hapData[0] == 'uniaxial':
1613	line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1614	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1615	print >>pFile,line
1616	else:
1617	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1618	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1619	Snames = ['S11','S22','S33','S12','S13','S23']
1620	ptlbls = ' names :'
1621	ptstr = ' values:'
1622	varstr = ' refine:'
1623	for i,name in enumerate(Snames):
1624	ptlbls += '%12s' % (name)
1625	ptstr += '%12.6f' % (hapData[4][i])
1626	varstr += '%12s' % (str(hapData[5][i]))
1627	print >>pFile,ptlbls
1628	print >>pFile,ptstr
1629	print >>pFile,varstr
1630
1631	def PrintMuStrain(hapData,SGData):
1632	if hapData[0] in ['isotropic','uniaxial']:
1633	line = '\n Mustrain model: %9s'%(hapData[0])
1634	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1635	if hapData[0] == 'uniaxial':
1636	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1637	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1638	print >>pFile,line
1639	else:
1640	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1641	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1642	Snames = G2spc.MustrainNames(SGData)
1643	ptlbls = ' names :'
1644	ptstr = ' values:'
1645	varstr = ' refine:'
1646	for i,name in enumerate(Snames):
1647	ptlbls += '%12s' % (name)
1648	ptstr += '%12.6f' % (hapData[4][i])
1649	varstr += '%12s' % (str(hapData[5][i]))
1650	print >>pFile,ptlbls
1651	print >>pFile,ptstr
1652	print >>pFile,varstr
1653
1654	def PrintHStrain(hapData,SGData):
1655	print >>pFile,'\n Hydrostatic/elastic strain: '
1656	Hsnames = G2spc.HStrainNames(SGData)
1657	ptlbls = ' names :'
1658	ptstr = ' values:'
1659	varstr = ' refine:'
1660	for i,name in enumerate(Hsnames):
1661	ptlbls += '%12s' % (name)
1662	ptstr += '%12.6f' % (hapData[0][i])
1663	varstr += '%12s' % (str(hapData[1][i]))
1664	print >>pFile,ptlbls
1665	print >>pFile,ptstr
1666	print >>pFile,varstr
1667
1668	def PrintSHPO(hapData):
1669	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1670	str(hapData[4])+' Refine? '+str(hapData[2])
1671	ptlbls = ' names :'
1672	ptstr = ' values:'
1673	for item in hapData[5]:
1674	ptlbls += '%12s'%(item)
1675	ptstr += '%12.3f'%(hapData[5][item])
1676	print >>pFile,ptlbls
1677	print >>pFile,ptstr
1678
1679	def PrintBabinet(hapData):
1680	print >>pFile,'\n Babinet form factor modification: '
1681	ptlbls = ' names :'
1682	ptstr = ' values:'
1683	varstr = ' refine:'
1684	for name in ['BabA','BabU']:
1685	ptlbls += '%12s' % (name)
1686	ptstr += '%12.6f' % (hapData[name][0])
1687	varstr += '%12s' % (str(hapData[name][1]))
1688	print >>pFile,ptlbls
1689	print >>pFile,ptstr
1690	print >>pFile,varstr
1691
1692
1693	hapDict = {}
1694	hapVary = []
1695	controlDict = {}
1696	poType = {}
1697	poAxes = {}
1698	spAxes = {}
1699	spType = {}
1700
1701	for phase in Phases:
1702	HistoPhase = Phases[phase]['Histograms']
1703	SGData = Phases[phase]['General']['SGData']
1704	cell = Phases[phase]['General']['Cell'][1:7]
1705	A = G2lat.cell2A(cell)
1706	pId = Phases[phase]['pId']
1707	histoList = HistoPhase.keys()
1708	histoList.sort()
1709	for histogram in histoList:
1710	try:
1711	Histogram = Histograms[histogram]
1712	except KeyError:
1713	#skip if histogram not included e.g. in a sequential refinement
1714	continue
1715	hapData = HistoPhase[histogram]
1716	hId = Histogram['hId']
1717	if 'PWDR' in histogram:
1718	limits = Histogram['Limits'][1]
1719	inst = Histogram['Instrument Parameters'][0]
1720	Zero = inst['Zero'][1]
1721	if 'C' in inst['Type'][1]:
1722	try:
1723	wave = inst['Lam'][1]
1724	except KeyError:
1725	wave = inst['Lam1'][1]
1726	dmin = wave/(2.0*sind(limits[1]/2.0))
1727	pfx = str(pId)+':'+str(hId)+':'
1728	for item in ['Scale','Extinction']:
1729	hapDict[pfx+item] = hapData[item][0]
1730	if hapData[item][1]:
1731	hapVary.append(pfx+item)
1732	names = G2spc.HStrainNames(SGData)
1733	for i,name in enumerate(names):
1734	hapDict[pfx+name] = hapData['HStrain'][0][i]
1735	if hapData['HStrain'][1][i]:
1736	hapVary.append(pfx+name)
1737	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1738	if hapData['Pref.Ori.'][0] == 'MD':
1739	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1740	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1741	if hapData['Pref.Ori.'][2]:
1742	hapVary.append(pfx+'MD')
1743	else: #'SH' spherical harmonics
1744	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1745	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1746	for item in hapData['Pref.Ori.'][5]:
1747	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1748	if hapData['Pref.Ori.'][2]:
1749	hapVary.append(pfx+item)
1750	for item in ['Mustrain','Size']:
1751	controlDict[pfx+item+'Type'] = hapData[item][0]
1752	hapDict[pfx+item+';mx'] = hapData[item][1][2]
1753	if hapData[item][2][2]:
1754	hapVary.append(pfx+item+';mx')
1755	if hapData[item][0] in ['isotropic','uniaxial']:
1756	hapDict[pfx+item+';i'] = hapData[item][1][0]
1757	if hapData[item][2][0]:
1758	hapVary.append(pfx+item+';i')
1759	if hapData[item][0] == 'uniaxial':
1760	controlDict[pfx+item+'Axis'] = hapData[item][3]
1761	hapDict[pfx+item+';a'] = hapData[item][1][1]
1762	if hapData[item][2][1]:
1763	hapVary.append(pfx+item+';a')
1764	else: #generalized for mustrain or ellipsoidal for size
1765	Nterms = len(hapData[item][4])
1766	if item == 'Mustrain':
1767	names = G2spc.MustrainNames(SGData)
1768	pwrs = []
1769	for name in names:
1770	h,k,l = name[1:]
1771	pwrs.append([int(h),int(k),int(l)])
1772	controlDict[pfx+'MuPwrs'] = pwrs
1773	for i in range(Nterms):
1774	sfx = ':'+str(i)
1775	hapDict[pfx+item+sfx] = hapData[item][4][i]
1776	if hapData[item][5][i]:
1777	hapVary.append(pfx+item+sfx)
1778	for bab in ['BabA','BabU']:
1779	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1780	if hapData['Babinet'][bab][1]:
1781	hapVary.append(pfx+bab)
1782
1783	if Print:
1784	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1785	print >>pFile,135*'-'
1786	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1787	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1788	if hapData['Pref.Ori.'][0] == 'MD':
1789	Ax = hapData['Pref.Ori.'][3]
1790	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1791	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1792	else: #'SH' for spherical harmonics
1793	PrintSHPO(hapData['Pref.Ori.'])
1794	PrintSize(hapData['Size'])
1795	PrintMuStrain(hapData['Mustrain'],SGData)
1796	PrintHStrain(hapData['HStrain'],SGData)
1797	if hapData['Babinet']['BabA'][0]:
1798	PrintBabinet(hapData['Babinet'])
1799	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1800	if resetRefList:
1801	refList = []
1802	Uniq = []
1803	Phi = []
1804	for h,k,l,d in HKLd:
1805	ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1806	mul *= 2 # for powder overlap of Friedel pairs
1807	if ext:
1808	continue
1809	if 'C' in inst['Type'][0]:
1810	pos = 2.0asind(wave/(2.0d))+Zero
1811	if limits[0] < pos < limits[1]:
1812	refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,0.0])
1813	Uniq.append(uniq)
1814	Phi.append(phi)
1815	else:
1816	raise ValueError
1817	Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{}}
1818	elif 'HKLF' in histogram:
1819	inst = Histogram['Instrument Parameters'][0]
1820	hId = Histogram['hId']
1821	hfx = ':%d:'%(hId)
1822	for item in inst:
1823	if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
1824	hapDict[hfx+item] = inst[item][1]
1825	pfx = str(pId)+':'+str(hId)+':'
1826	hapDict[pfx+'Scale'] = hapData['Scale'][0]
1827	if hapData['Scale'][1]:
1828	hapVary.append(pfx+'Scale')
1829
1830	extApprox,extType,extParms = hapData['Extinction']
1831	controlDict[pfx+'EType'] = extType
1832	controlDict[pfx+'EApprox'] = extApprox
1833	controlDict[pfx+'Tbar'] = extParms['Tbar']
1834	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1835	if 'Primary' in extType:
1836	Ekey = ['Ep',]
1837	elif 'I & II' in extType:
1838	Ekey = ['Eg','Es']
1839	elif 'Secondary Type II' == extType:
1840	Ekey = ['Es',]
1841	elif 'Secondary Type I' == extType:
1842	Ekey = ['Eg',]
1843	else: #'None'
1844	Ekey = []
1845	for eKey in Ekey:
1846	hapDict[pfx+eKey] = extParms[eKey][0]
1847	if extParms[eKey][1]:
1848	hapVary.append(pfx+eKey)
1849	for bab in ['BabA','BabU']:
1850	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1851	if hapData['Babinet'][bab][1]:
1852	hapVary.append(pfx+bab)
1853	if Print:
1854	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1855	print >>pFile,135*'-'
1856	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1857	if extType != 'None':
1858	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1859	text = ' Parameters :'
1860	for eKey in Ekey:
1861	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1862	print >>pFile,text
1863	if hapData['Babinet']['BabA'][0]:
1864	PrintBabinet(hapData['Babinet'])
1865	Histogram['Reflection Lists'] = phase
1866
1867	return hapVary,hapDict,controlDict
1868
1869	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1870	'needs a doc string'
1871
1872	def PrintSizeAndSig(hapData,sizeSig):
1873	line = '\n Size model: %9s'%(hapData[0])
1874	refine = False
1875	if hapData[0] in ['isotropic','uniaxial']:
1876	line += ' equatorial:%12.5f'%(hapData[1][0])
1877	if sizeSig[0][0]:
1878	line += ', sig:%8.4f'%(sizeSig[0][0])
1879	refine = True
1880	if hapData[0] == 'uniaxial':
1881	line += ' axial:%12.4f'%(hapData[1][1])
1882	if sizeSig[0][1]:
1883	refine = True
1884	line += ', sig:%8.4f'%(sizeSig[0][1])
1885	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1886	if sizeSig[0][2]:
1887	refine = True
1888	line += ', sig:%8.4f'%(sizeSig[0][2])
1889	if refine:
1890	print >>pFile,line
1891	else:
1892	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1893	if sizeSig[0][2]:
1894	refine = True
1895	line += ', sig:%8.4f'%(sizeSig[0][2])
1896	Snames = ['S11','S22','S33','S12','S13','S23']
1897	ptlbls = ' name :'
1898	ptstr = ' value :'
1899	sigstr = ' sig :'
1900	for i,name in enumerate(Snames):
1901	ptlbls += '%12s' % (name)
1902	ptstr += '%12.6f' % (hapData[4][i])
1903	if sizeSig[1][i]:
1904	refine = True
1905	sigstr += '%12.6f' % (sizeSig[1][i])
1906	else:
1907	sigstr += 12*' '
1908	if refine:
1909	print >>pFile,line
1910	print >>pFile,ptlbls
1911	print >>pFile,ptstr
1912	print >>pFile,sigstr
1913
1914	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1915	line = '\n Mustrain model: %9s'%(hapData[0])
1916	refine = False
1917	if hapData[0] in ['isotropic','uniaxial']:
1918	line += ' equatorial:%12.1f'%(hapData[1][0])
1919	if mustrainSig[0][0]:
1920	line += ', sig:%8.1f'%(mustrainSig[0][0])
1921	refine = True
1922	if hapData[0] == 'uniaxial':
1923	line += ' axial:%12.1f'%(hapData[1][1])
1924	if mustrainSig[0][1]:
1925	line += ', sig:%8.1f'%(mustrainSig[0][1])
1926	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1927	if mustrainSig[0][2]:
1928	refine = True
1929	line += ', sig:%8.3f'%(mustrainSig[0][2])
1930	if refine:
1931	print >>pFile,line
1932	else:
1933	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1934	if mustrainSig[0][2]:
1935	refine = True
1936	line += ', sig:%8.3f'%(mustrainSig[0][2])
1937	Snames = G2spc.MustrainNames(SGData)
1938	ptlbls = ' name :'
1939	ptstr = ' value :'
1940	sigstr = ' sig :'
1941	for i,name in enumerate(Snames):
1942	ptlbls += '%12s' % (name)
1943	ptstr += '%12.6f' % (hapData[4][i])
1944	if mustrainSig[1][i]:
1945	refine = True
1946	sigstr += '%12.6f' % (mustrainSig[1][i])
1947	else:
1948	sigstr += 12*' '
1949	if refine:
1950	print >>pFile,line
1951	print >>pFile,ptlbls
1952	print >>pFile,ptstr
1953	print >>pFile,sigstr
1954
1955	def PrintHStrainAndSig(hapData,strainSig,SGData):
1956	Hsnames = G2spc.HStrainNames(SGData)
1957	ptlbls = ' name :'
1958	ptstr = ' value :'
1959	sigstr = ' sig :'
1960	refine = False
1961	for i,name in enumerate(Hsnames):
1962	ptlbls += '%12s' % (name)
1963	ptstr += '%12.6g' % (hapData[0][i])
1964	if name in strainSig:
1965	refine = True
1966	sigstr += '%12.6g' % (strainSig[name])
1967	else:
1968	sigstr += 12*' '
1969	if refine:
1970	print >>pFile,'\n Hydrostatic/elastic strain: '
1971	print >>pFile,ptlbls
1972	print >>pFile,ptstr
1973	print >>pFile,sigstr
1974
1975	def PrintSHPOAndSig(pfx,hapData,POsig):
1976	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1977	ptlbls = ' names :'
1978	ptstr = ' values:'
1979	sigstr = ' sig :'
1980	for item in hapData[5]:
1981	ptlbls += '%12s'%(item)
1982	ptstr += '%12.3f'%(hapData[5][item])
1983	if pfx+item in POsig:
1984	sigstr += '%12.3f'%(POsig[pfx+item])
1985	else:
1986	sigstr += 12*' '
1987	print >>pFile,ptlbls
1988	print >>pFile,ptstr
1989	print >>pFile,sigstr
1990
1991	def PrintExtAndSig(pfx,hapData,ScalExtSig):
1992	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
1993	text = ''
1994	for item in hapData[2]:
1995	if pfx+item in ScalExtSig:
1996	text += ' %s: '%(item)
1997	text += '%12.2e'%(hapData[2][item][0])
1998	if pfx+item in ScalExtSig:
1999	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2000	print >>pFile,text
2001
2002	def PrintBabinetAndSig(pfx,hapData,BabSig):
2003	print >>pFile,'\n Babinet form factor modification: '
2004	ptlbls = ' names :'
2005	ptstr = ' values:'
2006	sigstr = ' sig :'
2007	for item in hapData:
2008	ptlbls += '%12s'%(item)
2009	ptstr += '%12.3f'%(hapData[item][0])
2010	if pfx+item in BabSig:
2011	sigstr += '%12.3f'%(BabSig[pfx+item])
2012	else:
2013	sigstr += 12*' '
2014	print >>pFile,ptlbls
2015	print >>pFile,ptstr
2016	print >>pFile,sigstr
2017
2018	PhFrExtPOSig = {}
2019	SizeMuStrSig = {}
2020	ScalExtSig = {}
2021	BabSig = {}
2022	wtFrSum = {}
2023	for phase in Phases:
2024	HistoPhase = Phases[phase]['Histograms']
2025	General = Phases[phase]['General']
2026	SGData = General['SGData']
2027	pId = Phases[phase]['pId']
2028	histoList = HistoPhase.keys()
2029	histoList.sort()
2030	for histogram in histoList:
2031	try:
2032	Histogram = Histograms[histogram]
2033	except KeyError:
2034	#skip if histogram not included e.g. in a sequential refinement
2035	continue
2036	hapData = HistoPhase[histogram]
2037	hId = Histogram['hId']
2038	pfx = str(pId)+':'+str(hId)+':'
2039	if hId not in wtFrSum:
2040	wtFrSum[hId] = 0.
2041	if 'PWDR' in histogram:
2042	for item in ['Scale','Extinction']:
2043	hapData[item][0] = parmDict[pfx+item]
2044	if pfx+item in sigDict:
2045	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2046	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2047	if hapData['Pref.Ori.'][0] == 'MD':
2048	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2049	if pfx+'MD' in sigDict:
2050	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2051	else: #'SH' spherical harmonics
2052	for item in hapData['Pref.Ori.'][5]:
2053	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2054	if pfx+item in sigDict:
2055	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2056	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2057	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2058	pfx+'HStrain':{}})
2059	for item in ['Mustrain','Size']:
2060	hapData[item][1][2] = parmDict[pfx+item+';mx']
2061	hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
2062	if pfx+item+';mx' in sigDict:
2063	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2064	if hapData[item][0] in ['isotropic','uniaxial']:
2065	hapData[item][1][0] = parmDict[pfx+item+';i']
2066	if item == 'Size':
2067	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2068	if pfx+item+';i' in sigDict:
2069	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2070	if hapData[item][0] == 'uniaxial':
2071	hapData[item][1][1] = parmDict[pfx+item+';a']
2072	if item == 'Size':
2073	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2074	if pfx+item+';a' in sigDict:
2075	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2076	else: #generalized for mustrain or ellipsoidal for size
2077	Nterms = len(hapData[item][4])
2078	for i in range(Nterms):
2079	sfx = ':'+str(i)
2080	hapData[item][4][i] = parmDict[pfx+item+sfx]
2081	if pfx+item+sfx in sigDict:
2082	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2083	names = G2spc.HStrainNames(SGData)
2084	for i,name in enumerate(names):
2085	hapData['HStrain'][0][i] = parmDict[pfx+name]
2086	if pfx+name in sigDict:
2087	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2088	for name in ['BabA','BabU']:
2089	hapData['Babinet'][name][0] = parmDict[pfx+name]
2090	if pfx+name in sigDict:
2091	BabSig[pfx+name] = sigDict[pfx+name]
2092
2093	elif 'HKLF' in histogram:
2094	for item in ['Scale',]:
2095	if parmDict.get(pfx+item):
2096	hapData[item][0] = parmDict[pfx+item]
2097	if pfx+item in sigDict:
2098	ScalExtSig[pfx+item] = sigDict[pfx+item]
2099	for item in ['Ep','Eg','Es']:
2100	if parmDict.get(pfx+item):
2101	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2102	if pfx+item in sigDict:
2103	ScalExtSig[pfx+item] = sigDict[pfx+item]
2104	for name in ['BabA','BabU']:
2105	hapData['Babinet'][name][0] = parmDict[pfx+name]
2106	if pfx+name in sigDict:
2107	BabSig[pfx+name] = sigDict[pfx+name]
2108
2109	if Print:
2110	for phase in Phases:
2111	HistoPhase = Phases[phase]['Histograms']
2112	General = Phases[phase]['General']
2113	SGData = General['SGData']
2114	pId = Phases[phase]['pId']
2115	histoList = HistoPhase.keys()
2116	histoList.sort()
2117	for histogram in histoList:
2118	try:
2119	Histogram = Histograms[histogram]
2120	except KeyError:
2121	#skip if histogram not included e.g. in a sequential refinement
2122	continue
2123	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2124	print >>pFile,130*'-'
2125	hapData = HistoPhase[histogram]
2126	hId = Histogram['hId']
2127	Histogram['Residuals'][str(pId)+'::Name'] = phase
2128	pfx = str(pId)+':'+str(hId)+':'
2129	if 'PWDR' in histogram:
2130	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2131	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2132
2133	if pfx+'Scale' in PhFrExtPOSig:
2134	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2135	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2136	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2137	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2138	if pfx+'Extinction' in PhFrExtPOSig:
2139	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2140	if hapData['Pref.Ori.'][0] == 'MD':
2141	if pfx+'MD' in PhFrExtPOSig:
2142	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2143	else:
2144	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2145	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2146	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2147	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2148	if len(BabSig):
2149	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2150
2151	elif 'HKLF' in histogram:
2152	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2153	%(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2154	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2155	if pfx+'Scale' in ScalExtSig:
2156	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2157	if hapData['Extinction'][0] != 'None':
2158	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2159	if len(BabSig):
2160	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2161
2162	################################################################################
2163	##### Histogram data
2164	################################################################################
2165
2166	def GetHistogramData(Histograms,Print=True,pFile=None):
2167	'needs a doc string'
2168
2169	def GetBackgroundParms(hId,Background):
2170	Back = Background[0]
2171	DebyePeaks = Background[1]
2172	bakType,bakFlag = Back[:2]
2173	backVals = Back[3:]
2174	backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2175	backDict = dict(zip(backNames,backVals))
2176	backVary = []
2177	if bakFlag:
2178	backVary = backNames
2179	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2180	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2181	debyeDict = {}
2182	debyeList = []
2183	for i in range(DebyePeaks['nDebye']):
2184	debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2185	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2186	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2187	debyeVary = []
2188	for item in debyeList:
2189	if item[1]:
2190	debyeVary.append(item[0])
2191	backDict.update(debyeDict)
2192	backVary += debyeVary
2193	peakDict = {}
2194	peakList = []
2195	for i in range(DebyePeaks['nPeaks']):
2196	peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2197	':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2198	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2199	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2200	peakVary = []
2201	for item in peakList:
2202	if item[1]:
2203	peakVary.append(item[0])
2204	backDict.update(peakDict)
2205	backVary += peakVary
2206	return bakType,backDict,backVary
2207
2208	def GetInstParms(hId,Inst):
2209	dataType = Inst['Type'][0]
2210	instDict = {}
2211	insVary = []
2212	pfx = ':'+str(hId)+':'
2213	insKeys = Inst.keys()
2214	insKeys.sort()
2215	for item in insKeys:
2216	insName = pfx+item
2217	instDict[insName] = Inst[item][1]
2218	if len(Inst[item]) > 2 and Inst[item][2]:
2219	insVary.append(insName)
2220	# instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2221	# instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2222	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2223	return dataType,instDict,insVary
2224
2225	def GetSampleParms(hId,Sample):
2226	sampVary = []
2227	hfx = ':'+str(hId)+':'
2228	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2229	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2230	for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2231	if key in Sample:
2232	sampDict[hfx+key] = Sample[key]
2233	Type = Sample['Type']
2234	if 'Bragg' in Type: #Bragg-Brentano
2235	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2236	sampDict[hfx+item] = Sample[item][0]
2237	if Sample[item][1]:
2238	sampVary.append(hfx+item)
2239	elif 'Debye' in Type: #Debye-Scherrer
2240	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2241	sampDict[hfx+item] = Sample[item][0]
2242	if Sample[item][1]:
2243	sampVary.append(hfx+item)
2244	return Type,sampDict,sampVary
2245
2246	def PrintBackground(Background):
2247	Back = Background[0]
2248	DebyePeaks = Background[1]
2249	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2250	line = ' Coefficients: '
2251	for i,back in enumerate(Back[3:]):
2252	line += '%10.3f'%(back)
2253	if i and not i%10:
2254	line += '\n'+15*' '
2255	print >>pFile,line
2256	if DebyePeaks['nDebye']:
2257	print >>pFile,'\n Debye diffuse scattering coefficients'
2258	parms = ['DebyeA','DebyeR','DebyeU']
2259	line = ' names : '
2260	for parm in parms:
2261	line += '%8s refine?'%(parm)
2262	print >>pFile,line
2263	for j,term in enumerate(DebyePeaks['debyeTerms']):
2264	line = ' term'+'%2d'%(j)+':'
2265	for i in range(3):
2266	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2267	print >>pFile,line
2268	if DebyePeaks['nPeaks']:
2269	print >>pFile,'\n Single peak coefficients'
2270	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2271	line = ' names : '
2272	for parm in parms:
2273	line += '%8s refine?'%(parm)
2274	print >>pFile,line
2275	for j,term in enumerate(DebyePeaks['peaksList']):
2276	line = ' peak'+'%2d'%(j)+':'
2277	for i in range(4):
2278	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2279	print >>pFile,line
2280
2281	def PrintInstParms(Inst):
2282	print >>pFile,'\n Instrument Parameters:'
2283	ptlbls = ' name :'
2284	ptstr = ' value :'
2285	varstr = ' refine:'
2286	insKeys = Inst.keys()
2287	insKeys.sort()
2288	for item in insKeys:
2289	if item != 'Type':
2290	ptlbls += '%12s' % (item)
2291	ptstr += '%12.6f' % (Inst[item][1])
2292	if item in ['Lam1','Lam2','Azimuth']:
2293	varstr += 12*' '
2294	else:
2295	varstr += '%12s' % (str(bool(Inst[item][2])))
2296	print >>pFile,ptlbls
2297	print >>pFile,ptstr
2298	print >>pFile,varstr
2299
2300	def PrintSampleParms(Sample):
2301	print >>pFile,'\n Sample Parameters:'
2302	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2303	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2304	ptlbls = ' name :'
2305	ptstr = ' value :'
2306	varstr = ' refine:'
2307	if 'Bragg' in Sample['Type']:
2308	for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2309	ptlbls += '%14s'%(item)
2310	ptstr += '%14.4f'%(Sample[item][0])
2311	varstr += '%14s'%(str(bool(Sample[item][1])))
2312
2313	elif 'Debye' in Type: #Debye-Scherrer
2314	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2315	ptlbls += '%14s'%(item)
2316	ptstr += '%14.4f'%(Sample[item][0])
2317	varstr += '%14s'%(str(bool(Sample[item][1])))
2318
2319	print >>pFile,ptlbls
2320	print >>pFile,ptstr
2321	print >>pFile,varstr
2322
2323	histDict = {}
2324	histVary = []
2325	controlDict = {}
2326	histoList = Histograms.keys()
2327	histoList.sort()
2328	for histogram in histoList:
2329	if 'PWDR' in histogram:
2330	Histogram = Histograms[histogram]
2331	hId = Histogram['hId']
2332	pfx = ':'+str(hId)+':'
2333	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2334	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2335	controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2336	for excl in controlDict[pfx+'Exclude']:
2337	Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2338	if controlDict[pfx+'Exclude']:
2339	ma.mask_rows(Histogram['Data'])
2340	Background = Histogram['Background']
2341	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2342	controlDict[pfx+'bakType'] = Type
2343	histDict.update(bakDict)
2344	histVary += bakVary
2345
2346	Inst = Histogram['Instrument Parameters'][0]
2347	Type,instDict,insVary = GetInstParms(hId,Inst)
2348	controlDict[pfx+'histType'] = Type
2349	if pfx+'Lam1' in instDict:
2350	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2351	else:
2352	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2353	histDict.update(instDict)
2354	histVary += insVary
2355
2356	Sample = Histogram['Sample Parameters']
2357	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2358	controlDict[pfx+'instType'] = Type
2359	histDict.update(sampDict)
2360	histVary += sampVary
2361
2362
2363	if Print:
2364	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2365	print >>pFile,135*'-'
2366	Units = {'C':' deg','T':' msec'}
2367	units = Units[controlDict[pfx+'histType'][2]]
2368	Limits = controlDict[pfx+'Limits']
2369	print >>pFile,' Instrument type: ',Sample['Type']
2370	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2371	if len(controlDict[pfx+'Exclude']):
2372	excls = controlDict[pfx+'Exclude']
2373	for excl in excls:
2374	print >>pFile,' Excluded region: %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)
2375	PrintSampleParms(Sample)
2376	PrintInstParms(Inst)
2377	PrintBackground(Background)
2378	elif 'HKLF' in histogram:
2379	Histogram = Histograms[histogram]
2380	hId = Histogram['hId']
2381	pfx = ':'+str(hId)+':'
2382	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2383	Inst = Histogram['Instrument Parameters'][0]
2384	controlDict[pfx+'histType'] = Inst['Type'][0]
2385	histDict[pfx+'Lam'] = Inst['Lam'][1]
2386	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
2387	return histVary,histDict,controlDict
2388
2389	def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2390	'needs a doc string'
2391
2392	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2393	Back = Background[0]
2394	DebyePeaks = Background[1]
2395	lenBack = len(Back[3:])
2396	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
2397	for i in range(lenBack):
2398	Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2399	if pfx+'Back:'+str(i) in sigDict:
2400	backSig[i] = sigDict[pfx+'Back:'+str(i)]
2401	if DebyePeaks['nDebye']:
2402	for i in range(DebyePeaks['nDebye']):
2403	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2404	for j,name in enumerate(names):
2405	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2406	if name in sigDict:
2407	backSig[lenBack+3*i+j] = sigDict[name]
2408	if DebyePeaks['nPeaks']:
2409	for i in range(DebyePeaks['nPeaks']):
2410	names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2411	pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2412	for j,name in enumerate(names):
2413	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2414	if name in sigDict:
2415	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
2416	return backSig
2417
2418	def SetInstParms(pfx,Inst,parmDict,sigDict):
2419	instSig = {}
2420	insKeys = Inst.keys()
2421	insKeys.sort()
2422	for item in insKeys:
2423	insName = pfx+item
2424	Inst[item][1] = parmDict[insName]
2425	if insName in sigDict:
2426	instSig[item] = sigDict[insName]
2427	else:
2428	instSig[item] = 0
2429	return instSig
2430
2431	def SetSampleParms(pfx,Sample,parmDict,sigDict):
2432	if 'Bragg' in Sample['Type']: #Bragg-Brentano
2433	sampSig = [0 for i in range(5)]
2434	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2435	Sample[item][0] = parmDict[pfx+item]
2436	if pfx+item in sigDict:
2437	sampSig[i] = sigDict[pfx+item]
2438	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2439	sampSig = [0 for i in range(4)]
2440	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2441	Sample[item][0] = parmDict[pfx+item]
2442	if pfx+item in sigDict:
2443	sampSig[i] = sigDict[pfx+item]
2444	return sampSig
2445
2446	def PrintBackgroundSig(Background,backSig):
2447	Back = Background[0]
2448	DebyePeaks = Background[1]
2449	lenBack = len(Back[3:])
2450	valstr = ' value : '
2451	sigstr = ' sig : '
2452	refine = False
2453	for i,back in enumerate(Back[3:]):
2454	valstr += '%10.4g'%(back)
2455	if Back[1]:
2456	refine = True
2457	sigstr += '%10.4g'%(backSig[i])
2458	else:
2459	sigstr += 10*' '
2460	if refine:
2461	print >>pFile,'\n Background function: ',Back[0]
2462	print >>pFile,valstr
2463	print >>pFile,sigstr
2464	if DebyePeaks['nDebye']:
2465	ifAny = False
2466	ptfmt = "%12.3f"
2467	names = ' names :'
2468	ptstr = ' values:'
2469	sigstr = ' esds :'
2470	for item in sigDict:
2471	if 'Debye' in item:
2472	ifAny = True
2473	names += '%12s'%(item)
2474	ptstr += ptfmt%(parmDict[item])
2475	sigstr += ptfmt%(sigDict[item])
2476	if ifAny:
2477	print >>pFile,'\n Debye diffuse scattering coefficients'
2478	print >>pFile,names
2479	print >>pFile,ptstr
2480	print >>pFile,sigstr
2481	if DebyePeaks['nPeaks']:
2482	ifAny = False
2483	ptfmt = "%14.3f"
2484	names = ' names :'
2485	ptstr = ' values:'
2486	sigstr = ' esds :'
2487	for item in sigDict:
2488	if 'BkPk' in item:
2489	ifAny = True
2490	names += '%14s'%(item)
2491	ptstr += ptfmt%(parmDict[item])
2492	sigstr += ptfmt%(sigDict[item])
2493	if ifAny:
2494	print >>pFile,'\n Single peak coefficients'
2495	print >>pFile,names
2496	print >>pFile,ptstr
2497	print >>pFile,sigstr
2498
2499	def PrintInstParmsSig(Inst,instSig):
2500	ptlbls = ' names :'
2501	ptstr = ' value :'
2502	sigstr = ' sig :'
2503	refine = False
2504	insKeys = instSig.keys()
2505	insKeys.sort()
2506	for name in insKeys:
2507	if name not in ['Type','Lam1','Lam2','Azimuth']:
2508	ptlbls += '%12s' % (name)
2509	ptstr += '%12.6f' % (Inst[name][1])
2510	if instSig[name]:
2511	refine = True
2512	sigstr += '%12.6f' % (instSig[name])
2513	else:
2514	sigstr += 12*' '
2515	if refine:
2516	print >>pFile,'\n Instrument Parameters:'
2517	print >>pFile,ptlbls
2518	print >>pFile,ptstr
2519	print >>pFile,sigstr
2520
2521	def PrintSampleParmsSig(Sample,sampleSig):
2522	ptlbls = ' names :'
2523	ptstr = ' values:'
2524	sigstr = ' sig :'
2525	refine = False
2526	if 'Bragg' in Sample['Type']:
2527	for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2528	ptlbls += '%14s'%(item)
2529	ptstr += '%14.4f'%(Sample[item][0])
2530	if sampleSig[i]:
2531	refine = True
2532	sigstr += '%14.4f'%(sampleSig[i])
2533	else:
2534	sigstr += 14*' '
2535
2536	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2537	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2538	ptlbls += '%14s'%(item)
2539	ptstr += '%14.4f'%(Sample[item][0])
2540	if sampleSig[i]:
2541	refine = True
2542	sigstr += '%14.4f'%(sampleSig[i])
2543	else:
2544	sigstr += 14*' '
2545
2546	if refine:
2547	print >>pFile,'\n Sample Parameters:'
2548	print >>pFile,ptlbls
2549	print >>pFile,ptstr
2550	print >>pFile,sigstr
2551
2552	histoList = Histograms.keys()
2553	histoList.sort()
2554	for histogram in histoList:
2555	if 'PWDR' in histogram:
2556	Histogram = Histograms[histogram]
2557	hId = Histogram['hId']
2558	pfx = ':'+str(hId)+':'
2559	Background = Histogram['Background']
2560	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2561
2562	Inst = Histogram['Instrument Parameters'][0]
2563	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2564
2565	Sample = Histogram['Sample Parameters']
2566	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2567
2568	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2569	print >>pFile,135*'-'
2570	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2571	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
2572	print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
2573	(Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
2574	if Print:
2575	print >>pFile,' Instrument type: ',Sample['Type']
2576	PrintSampleParmsSig(Sample,sampSig)
2577	PrintInstParmsSig(Inst,instSig)
2578	PrintBackgroundSig(Background,backSig)
2579

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

Download in other formats: