Context Navigation

source: trunk/GSASIIstrIO.py @ 1144

Last change on this file since 1144 was 1144, checked in by vondreele, 9 years ago
fix bug in GetControls?
Property svn:keywords set to `Date Author Revision URL Id`
File size: 112.8 KB

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

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

Download in other formats: