Context Navigation

source: trunk/GSASIIstruct.py @ 918

Last change on this file since 918 was 918, checked in by vondreele, 10 years ago
add delt/sig column in some restraint tables change testDeriv file name
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 221.0 KB

Line
1	# -- coding: utf-8 --
2	#GSASIIstructure - structure computation routines
3	########### SVN repository information ###################
4	# $Date: 2013-05-17 16:27:25 +0000 (Fri, 17 May 2013) $
5	# $Author: vondreele $
6	# $Revision: 918 $
7	# $URL: trunk/GSASIIstruct.py $
8	# $Id: GSASIIstruct.py 918 2013-05-17 16:27:25Z vondreele $
9	########### SVN repository information ###################
10	import sys
11	import os
12	import os.path as ospath
13	import time
14	import math
15	import copy
16	import random
17	import cPickle
18	import numpy as np
19	import numpy.ma as ma
20	import numpy.linalg as nl
21	import scipy.optimize as so
22	import GSASIIpath
23	GSASIIpath.SetVersionNumber("$Revision: 918 $")
24	import GSASIIElem as G2el
25	import GSASIIlattice as G2lat
26	import GSASIIspc as G2spc
27	import GSASIIpwd as G2pwd
28	import GSASIImapvars as G2mv
29	import GSASIImath as G2mth
30
31	sind = lambda x: np.sin(x*np.pi/180.)
32	cosd = lambda x: np.cos(x*np.pi/180.)
33	tand = lambda x: np.tan(x*np.pi/180.)
34	asind = lambda x: 180.*np.arcsin(x)/np.pi
35	acosd = lambda x: 180.*np.arccos(x)/np.pi
36	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
37
38	ateln2 = 8.0*math.log(2.0)
39	DEBUG = True
40
41	def GetControls(GPXfile):
42	''' Returns dictionary of control items found in GSASII gpx file
43	input:
44	GPXfile = .gpx full file name
45	return:
46	Controls = dictionary of control items
47	'''
48	Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
49	'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
50	fl = open(GPXfile,'rb')
51	while True:
52	try:
53	data = cPickle.load(fl)
54	except EOFError:
55	break
56	datum = data[0]
57	if datum[0] == 'Controls':
58	Controls.update(datum[1])
59	fl.close()
60	return Controls
61
62	def GetConstraints(GPXfile):
63	'''Read the constraints from the GPX file and interpret them
64	'''
65	constList = []
66	fl = open(GPXfile,'rb')
67	while True:
68	try:
69	data = cPickle.load(fl)
70	except EOFError:
71	break
72	datum = data[0]
73	if datum[0] == 'Constraints':
74	constDict = datum[1]
75	for item in constDict:
76	constList += constDict[item]
77	fl.close()
78	constDict,fixedList,ignored = ProcessConstraints(constList)
79	if ignored:
80	print ignored,'old-style Constraints were rejected'
81	return constDict,fixedList
82
83	def ProcessConstraints(constList):
84	"interpret constraints"
85	constDict = []
86	fixedList = []
87	ignored = 0
88	for item in constList:
89	if item[-1] == 'h':
90	# process a hold
91	fixedList.append('0')
92	constDict.append({item[0][1]:0.0})
93	elif item[-1] == 'f':
94	# process a new variable
95	fixedList.append(None)
96	constDict.append({})
97	for term in item[:-3]:
98	constDict[-1][term[1]] = term[0]
99	#constFlag[-1] = ['Vary']
100	elif item[-1] == 'c':
101	# process a contraint relationship
102	fixedList.append(str(item[-3]))
103	constDict.append({})
104	for term in item[:-3]:
105	constDict[-1][term[1]] = term[0]
106	#constFlag[-1] = ['VaryFree']
107	elif item[-1] == 'e':
108	# process an equivalence
109	firstmult = None
110	eqlist = []
111	for term in item[:-3]:
112	if term[0] == 0: term[0] = 1.0
113	if firstmult is None:
114	firstmult,firstvar = term
115	else:
116	eqlist.append([term[1],firstmult/term[0]])
117	G2mv.StoreEquivalence(firstvar,eqlist)
118	else:
119	ignored += 1
120	return constDict,fixedList,ignored
121
122	def CheckConstraints(GPXfile):
123	'''Load constraints and related info and return any error or warning messages'''
124	# init constraints
125	G2mv.InitVars()
126	# get variables
127	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
128	if not Phases:
129	return 'Error: No Phases!',''
130	if not Histograms:
131	return 'Error: no diffraction data',''
132	rigidbodyDict = GetRigidBodies(GPXfile)
133	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
134	rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
135	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)
136	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
137	histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
138	varyList = rbVary+phaseVary+hapVary+histVary
139	constrDict,fixedList = GetConstraints(GPXfile)
140	return G2mv.CheckConstraints(varyList,constrDict,fixedList)
141
142	def GetRestraints(GPXfile):
143	'''Read the restraints from the GPX file
144	'''
145	fl = open(GPXfile,'rb')
146	while True:
147	try:
148	data = cPickle.load(fl)
149	except EOFError:
150	break
151	datum = data[0]
152	if datum[0] == 'Restraints':
153	restraintDict = datum[1]
154	fl.close()
155	return restraintDict
156
157	def GetRigidBodies(GPXfile):
158	'''Read the rigid body models from the GPX file
159	'''
160	fl = open(GPXfile,'rb')
161	while True:
162	try:
163	data = cPickle.load(fl)
164	except EOFError:
165	break
166	datum = data[0]
167	if datum[0] == 'Rigid bodies':
168	rigidbodyDict = datum[1]
169	fl.close()
170	return rigidbodyDict
171
172	def GetPhaseNames(GPXfile):
173	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
174	input:
175	GPXfile = gpx full file name
176	return:
177	PhaseNames = list of phase names
178	'''
179	fl = open(GPXfile,'rb')
180	PhaseNames = []
181	while True:
182	try:
183	data = cPickle.load(fl)
184	except EOFError:
185	break
186	datum = data[0]
187	if 'Phases' == datum[0]:
188	for datus in data[1:]:
189	PhaseNames.append(datus[0])
190	fl.close()
191	return PhaseNames
192
193	def GetAllPhaseData(GPXfile,PhaseName):
194	''' Returns the entire dictionary for PhaseName from GSASII gpx file
195	input:
196	GPXfile = gpx full file name
197	PhaseName = phase name
198	return:
199	phase dictionary
200	'''
201	fl = open(GPXfile,'rb')
202	General = {}
203	Atoms = []
204	while True:
205	try:
206	data = cPickle.load(fl)
207	except EOFError:
208	break
209	datum = data[0]
210	if 'Phases' == datum[0]:
211	for datus in data[1:]:
212	if datus[0] == PhaseName:
213	break
214	fl.close()
215	return datus[1]
216
217	def GetHistograms(GPXfile,hNames):
218	""" Returns a dictionary of histograms found in GSASII gpx file
219	input:
220	GPXfile = .gpx full file name
221	hNames = list of histogram names
222	return:
223	Histograms = dictionary of histograms (types = PWDR & HKLF)
224	"""
225	fl = open(GPXfile,'rb')
226	Histograms = {}
227	while True:
228	try:
229	data = cPickle.load(fl)
230	except EOFError:
231	break
232	datum = data[0]
233	hist = datum[0]
234	if hist in hNames:
235	if 'PWDR' in hist[:4]:
236	PWDRdata = {}
237	try:
238	PWDRdata.update(datum[1][0]) #weight factor
239	except ValueError:
240	PWDRdata['wtFactor'] = 1.0 #patch
241	PWDRdata['Data'] = datum[1][1] #powder data arrays
242	PWDRdata[data[2][0]] = data[2][1] #Limits
243	PWDRdata[data[3][0]] = data[3][1] #Background
244	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
245	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
246	try:
247	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
248	except IndexError:
249	PWDRdata['Reflection Lists'] = {}
250
251	Histograms[hist] = PWDRdata
252	elif 'HKLF' in hist[:4]:
253	HKLFdata = {}
254	try:
255	HKLFdata.update(datum[1][0]) #weight factor
256	except ValueError:
257	HKLFdata['wtFactor'] = 1.0 #patch
258	HKLFdata['Data'] = datum[1][1]
259	HKLFdata[data[1][0]] = data[1][1] #Instrument parameters
260	HKLFdata['Reflection Lists'] = None
261	Histograms[hist] = HKLFdata
262	fl.close()
263	return Histograms
264
265	def GetHistogramNames(GPXfile,hType):
266	""" Returns a list of histogram names found in GSASII gpx file
267	input:
268	GPXfile = .gpx full file name
269	hType = list ['PWDR','HKLF']
270	return:
271	HistogramNames = list of histogram names (types = PWDR & HKLF)
272	"""
273	fl = open(GPXfile,'rb')
274	HistogramNames = []
275	while True:
276	try:
277	data = cPickle.load(fl)
278	except EOFError:
279	break
280	datum = data[0]
281	if datum[0][:4] in hType:
282	HistogramNames.append(datum[0])
283	fl.close()
284	return HistogramNames
285
286	def GetUsedHistogramsAndPhases(GPXfile):
287	''' Returns all histograms that are found in any phase
288	and any phase that uses a histogram
289	input:
290	GPXfile = .gpx full file name
291	return:
292	Histograms = dictionary of histograms as {name:data,...}
293	Phases = dictionary of phases that use histograms
294	'''
295	phaseNames = GetPhaseNames(GPXfile)
296	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
297	allHistograms = GetHistograms(GPXfile,histoList)
298	phaseData = {}
299	for name in phaseNames:
300	phaseData[name] = GetAllPhaseData(GPXfile,name)
301	Histograms = {}
302	Phases = {}
303	for phase in phaseData:
304	Phase = phaseData[phase]
305	if Phase['Histograms']:
306	if phase not in Phases:
307	pId = phaseNames.index(phase)
308	Phase['pId'] = pId
309	Phases[phase] = Phase
310	for hist in Phase['Histograms']:
311	if 'Use' not in Phase['Histograms'][hist]: #patch
312	Phase['Histograms'][hist]['Use'] = True
313	if hist not in Histograms and Phase['Histograms'][hist]['Use']:
314	Histograms[hist] = allHistograms[hist]
315	hId = histoList.index(hist)
316	Histograms[hist]['hId'] = hId
317	return Histograms,Phases
318
319	def getBackupName(GPXfile,makeBack):
320	GPXpath,GPXname = ospath.split(GPXfile)
321	if GPXpath == '': GPXpath = '.'
322	Name = ospath.splitext(GPXname)[0]
323	files = os.listdir(GPXpath)
324	last = 0
325	for name in files:
326	name = name.split('.')
327	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
328	if makeBack:
329	last = max(last,int(name[1].strip('bak'))+1)
330	else:
331	last = max(last,int(name[1].strip('bak')))
332	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
333	return GPXback
334
335	def GPXBackup(GPXfile,makeBack=True):
336	import distutils.file_util as dfu
337	GPXback = getBackupName(GPXfile,makeBack)
338	dfu.copy_file(GPXfile,GPXback)
339	return GPXback
340
341	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
342	''' Updates gpxfile from all histograms that are found in any phase
343	and any phase that used a histogram. Also updates rigid body definitions.
344	input:
345	GPXfile = .gpx full file name
346	Histograms = dictionary of histograms as {name:data,...}
347	Phases = dictionary of phases that use histograms
348	RigidBodies = dictionary of rigid bodies
349	CovData = dictionary of refined variables, varyList, & covariance matrix
350	makeBack = True if new backup of .gpx file is to be made; else use the last one made
351	'''
352
353	GPXback = GPXBackup(GPXfile,makeBack)
354	print 'Read from file:',GPXback
355	print 'Save to file :',GPXfile
356	infile = open(GPXback,'rb')
357	outfile = open(GPXfile,'wb')
358	while True:
359	try:
360	data = cPickle.load(infile)
361	except EOFError:
362	break
363	datum = data[0]
364	# print 'read: ',datum[0]
365	if datum[0] == 'Phases':
366	for iphase in range(len(data)):
367	if data[iphase][0] in Phases:
368	phaseName = data[iphase][0]
369	data[iphase][1].update(Phases[phaseName])
370	elif datum[0] == 'Covariance':
371	data[0][1] = CovData
372	elif datum[0] == 'Rigid bodies':
373	data[0][1] = RigidBodies
374	try:
375	histogram = Histograms[datum[0]]
376	# print 'found ',datum[0]
377	data[0][1][1] = histogram['Data']
378	for datus in data[1:]:
379	# print ' read: ',datus[0]
380	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
381	datus[1] = histogram[datus[0]]
382	except KeyError:
383	pass
384
385	cPickle.dump(data,outfile,1)
386	infile.close()
387	outfile.close()
388	print 'GPX file save successful'
389
390	def SetSeqResult(GPXfile,Histograms,SeqResult):
391	GPXback = GPXBackup(GPXfile)
392	print 'Read from file:',GPXback
393	print 'Save to file :',GPXfile
394	infile = open(GPXback,'rb')
395	outfile = open(GPXfile,'wb')
396	while True:
397	try:
398	data = cPickle.load(infile)
399	except EOFError:
400	break
401	datum = data[0]
402	if datum[0] == 'Sequential results':
403	data[0][1] = SeqResult
404	try:
405	histogram = Histograms[datum[0]]
406	data[0][1][1] = histogram['Data']
407	for datus in data[1:]:
408	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
409	datus[1] = histogram[datus[0]]
410	except KeyError:
411	pass
412
413	cPickle.dump(data,outfile,1)
414	infile.close()
415	outfile.close()
416	print 'GPX file save successful'
417
418	def ShowBanner(pFile=None):
419	print >>pFile,80''
420	print >>pFile,' General Structure Analysis System-II Crystal Structure Refinement'
421	print >>pFile,' by Robert B. Von Dreele & Brian H. Toby'
422	print >>pFile,' Argonne National Laboratory(C), 2010'
423	print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,'
424	print >>pFile,' as Operator of Argonne National Laboratory.'
425	print >>pFile,' Please cite:'
426	print >>pFile,' B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)'
427
428	print >>pFile,80'','\n'
429
430	def ShowControls(Controls,pFile=None):
431	print >>pFile,' Least squares controls:'
432	print >>pFile,' Refinement type: ',Controls['deriv type']
433	if 'Hessian' in Controls['deriv type']:
434	print >>pFile,' Maximum number of cycles:',Controls['max cyc']
435	else:
436	print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
437	print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
438
439	def GetFFtable(General):
440	''' returns a dictionary of form factor data for atom types found in General
441	input:
442	General = dictionary of phase info.; includes AtomTypes
443	return:
444	FFtable = dictionary of form factor data; key is atom type
445	'''
446	atomTypes = General['AtomTypes']
447	FFtable = {}
448	for El in atomTypes:
449	FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
450	for item in FFs:
451	if item['Symbol'] == El.upper():
452	FFtable[El] = item
453	return FFtable
454
455	def GetBLtable(General):
456	''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
457	input:
458	General = dictionary of phase info.; includes AtomTypes & Isotopes
459	return:
460	BLtable = dictionary of scattering length data; key is atom type
461	'''
462	atomTypes = General['AtomTypes']
463	BLtable = {}
464	isotopes = General['Isotopes']
465	isotope = General['Isotope']
466	for El in atomTypes:
467	BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
468	return BLtable
469
470	def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
471	# if SGLaue in ['-1','2/m','mmm']:
472	# return #no Pawley symmetry required constraints
473	eqvDict = {}
474	for i,varyI in enumerate(pawleyVary):
475	eqvDict[varyI] = []
476	refI = int(varyI.split(':')[-1])
477	ih,ik,il = PawleyRef[refI][:3]
478	dspI = PawleyRef[refI][4]
479	for varyJ in pawleyVary[i+1:]:
480	refJ = int(varyJ.split(':')[-1])
481	jh,jk,jl = PawleyRef[refJ][:3]
482	dspJ = PawleyRef[refJ][4]
483	if SGLaue in ['4/m','4/mmm']:
484	isum = ih2+ik2
485	jsum = jh2+jk2
486	if abs(il) == abs(jl) and isum == jsum:
487	eqvDict[varyI].append(varyJ)
488	elif SGLaue in ['3R','3mR']:
489	isum = ih2+ik2+il**2
490	jsum = jh2+jk2jl*2
491	isum2 = ihik+ihil+ik*il
492	jsum2 = jhjk+jhjl+jk*jl
493	if isum == jsum and isum2 == jsum2:
494	eqvDict[varyI].append(varyJ)
495	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
496	isum = ih2+ik2+ih*ik
497	jsum = jh2+jk2+jh*jk
498	if abs(il) == abs(jl) and isum == jsum:
499	eqvDict[varyI].append(varyJ)
500	elif SGLaue in ['m3','m3m']:
501	isum = ih2+ik2+il**2
502	jsum = jh2+jk2+jl**2
503	if isum == jsum:
504	eqvDict[varyI].append(varyJ)
505	elif abs(dspI-dspJ)/dspI < 1.e-4:
506	eqvDict[varyI].append(varyJ)
507	for item in pawleyVary:
508	if eqvDict[item]:
509	for item2 in pawleyVary:
510	if item2 in eqvDict[item]:
511	eqvDict[item2] = []
512	G2mv.StoreEquivalence(item,eqvDict[item])
513
514	def cellVary(pfx,SGData):
515	if SGData['SGLaue'] in ['-1',]:
516	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
517	elif SGData['SGLaue'] in ['2/m',]:
518	if SGData['SGUniq'] == 'a':
519	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
520	elif SGData['SGUniq'] == 'b':
521	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
522	else:
523	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
524	elif SGData['SGLaue'] in ['mmm',]:
525	return [pfx+'A0',pfx+'A1',pfx+'A2']
526	elif SGData['SGLaue'] in ['4/m','4/mmm']:
527	return [pfx+'A0',pfx+'A2']
528	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
529	return [pfx+'A0',pfx+'A2']
530	elif SGData['SGLaue'] in ['3R', '3mR']:
531	return [pfx+'A0',pfx+'A3']
532	elif SGData['SGLaue'] in ['m3m','m3']:
533	return [pfx+'A0',]
534
535	################################################################################
536	##### Rigid Body Models
537	################################################################################
538
539	def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
540
541	def PrintResRBModel(RBModel):
542	atNames = RBModel['atNames']
543	rbRef = RBModel['rbRef']
544	rbSeq = RBModel['rbSeq']
545	print >>pFile,'Residue RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
546	'No. times used: ',RBModel['useCount']
547	print >>pFile,' At name x y z'
548	for name,xyz in zip(atNames,RBModel['rbXYZ']):
549	print >>pFile,' %8s %10.4f %10.4f %10.4f'%(name,xyz[0],xyz[1],xyz[2])
550	print >>pFile,'Orientation defined by:',atNames[rbRef[0]],' -> ',atNames[rbRef[1]], \
551	' & ',atNames[rbRef[0]],' -> ',atNames[rbRef[2]]
552	if rbSeq:
553	for i,rbseq in enumerate(rbSeq):
554	print >>pFile,'Torsion sequence ',i,' Bond: '+atNames[rbseq[0]],' - ', \
555	atNames[rbseq[1]],' riding: ',[atNames[i] for i in rbseq[3]]
556
557	def PrintVecRBModel(RBModel):
558	rbRef = RBModel['rbRef']
559	atTypes = RBModel['rbTypes']
560	print >>pFile,'Vector RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
561	'No. times used: ',RBModel['useCount']
562	for i in range(len(RBModel['VectMag'])):
563	print >>pFile,'Vector no.: ',i,' Magnitude: ', \
564	'%8.4f'%(RBModel['VectMag'][i]),' Refine? ',RBModel['VectRef'][i]
565	print >>pFile,' No. Type vx vy vz'
566	for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
567	print >>pFile,' %d %2s %10.4f %10.4f %10.4f'%(j,name,xyz[0],xyz[1],xyz[2])
568	print >>pFile,' No. Type x y z'
569	for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
570	print >>pFile,' %d %2s %10.4f %10.4f %10.4f'%(i,name,xyz[0],xyz[1],xyz[2])
571	print >>pFile,'Orientation defined by: atom ',rbRef[0],' -> atom ',rbRef[1], \
572	' & atom ',rbRef[0],' -> atom ',rbRef[2]
573	rbVary = []
574	rbDict = {}
575	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
576	if len(rbIds['Vector']):
577	for irb,item in enumerate(rbIds['Vector']):
578	if rigidbodyDict['Vector'][item]['useCount']:
579	RBmags = rigidbodyDict['Vector'][item]['VectMag']
580	RBrefs = rigidbodyDict['Vector'][item]['VectRef']
581	for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
582	pid = '::RBV;'+str(i)+':'+str(irb)
583	rbDict[pid] = mag
584	if ref:
585	rbVary.append(pid)
586	if Print:
587	print >>pFile,'\nVector rigid body model:'
588	PrintVecRBModel(rigidbodyDict['Vector'][item])
589	if len(rbIds['Residue']):
590	for item in rbIds['Residue']:
591	if rigidbodyDict['Residue'][item]['useCount']:
592	if Print:
593	print >>pFile,'\nResidue rigid body model:'
594	PrintResRBModel(rigidbodyDict['Residue'][item])
595	return rbVary,rbDict
596
597	def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
598
599	def PrintRBVectandSig(VectRB,VectSig):
600	print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
601	namstr = ' names :'
602	valstr = ' values:'
603	sigstr = ' esds :'
604	for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
605	namstr += '%12s'%('Vect '+str(i))
606	valstr += '%12.4f'%(val)
607	if sig:
608	sigstr += '%12.4f'%(sig)
609	else:
610	sigstr += 12*' '
611	print >>pFile,namstr
612	print >>pFile,valstr
613	print >>pFile,sigstr
614
615	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
616	if not RBIds['Vector']:
617	return
618	for irb,item in enumerate(RBIds['Vector']):
619	if rigidbodyDict['Vector'][item]['useCount']:
620	VectSig = []
621	RBmags = rigidbodyDict['Vector'][item]['VectMag']
622	for i,mag in enumerate(RBmags):
623	name = '::RBV;'+str(i)+':'+str(irb)
624	mag = parmDict[name]
625	if name in sigDict:
626	VectSig.append(sigDict[name])
627	PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)
628
629	def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
630	''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
631	current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
632	'''
633	atxIds = ['Ax:','Ay:','Az:']
634	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
635	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
636	if not RBIds['Vector'] and not RBIds['Residue']:
637	return
638	VRBIds = RBIds['Vector']
639	RRBIds = RBIds['Residue']
640	if Update:
641	RBData = rigidbodyDict
642	else:
643	RBData = copy.deepcopy(rigidbodyDict) # don't mess with original!
644	if RBIds['Vector']: # first update the vector magnitudes
645	VRBData = RBData['Vector']
646	for i,rbId in enumerate(VRBIds):
647	for j in range(len(VRBData[rbId]['VectMag'])):
648	name = '::RBV;'+str(j)+':'+str(i)
649	VRBData[rbId]['VectMag'][j] = parmDict[name]
650	for phase in Phases:
651	Phase = Phases[phase]
652	General = Phase['General']
653	cell = General['Cell'][1:7]
654	Amat,Bmat = G2lat.cell2AB(cell)
655	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
656	pfx = str(Phase['pId'])+'::'
657	if Update:
658	RBModels = Phase['RBModels']
659	else:
660	RBModels = copy.deepcopy(Phase['RBModels']) # again don't mess with original!
661	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
662	jrb = VRBIds.index(RBObj['RBId'])
663	rbsx = str(irb)+':'+str(jrb)
664	for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
665	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
666	for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
667	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
668	TLS = RBObj['ThermalMotion']
669	if 'T' in TLS[0]:
670	for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
671	TLS[1][i] = parmDict[pfx+pt+rbsx]
672	if 'L' in TLS[0]:
673	for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
674	TLS[1][i+6] = parmDict[pfx+pt+rbsx]
675	if 'S' in TLS[0]:
676	for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
677	TLS[1][i+12] = parmDict[pfx+pt+rbsx]
678	if 'U' in TLS[0]:
679	TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
680	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
681	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
682	for i,x in enumerate(XYZ):
683	atId = RBObj['Ids'][i]
684	for j in [0,1,2]:
685	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
686	if UIJ[i][0] == 'A':
687	for j in range(6):
688	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
689	elif UIJ[i][0] == 'I':
690	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
691
692	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
693	jrb = RRBIds.index(RBObj['RBId'])
694	rbsx = str(irb)+':'+str(jrb)
695	for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
696	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
697	for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
698	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
699	TLS = RBObj['ThermalMotion']
700	if 'T' in TLS[0]:
701	for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
702	RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
703	if 'L' in TLS[0]:
704	for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
705	RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
706	if 'S' in TLS[0]:
707	for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
708	RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
709	if 'U' in TLS[0]:
710	RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
711	for itors,tors in enumerate(RBObj['Torsions']):
712	tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
713	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
714	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
715	for i,x in enumerate(XYZ):
716	atId = RBObj['Ids'][i]
717	for j in [0,1,2]:
718	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
719	if UIJ[i][0] == 'A':
720	for j in range(6):
721	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
722	elif UIJ[i][0] == 'I':
723	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
724
725	def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
726	atxIds = ['dAx:','dAy:','dAz:']
727	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
728	TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
729	LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
730	SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
731	PIds = ['Px:','Py:','Pz:']
732	OIds = ['Oa:','Oi:','Oj:','Ok:']
733	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
734	if not RBIds['Vector'] and not RBIds['Residue']:
735	return
736	VRBIds = RBIds['Vector']
737	RRBIds = RBIds['Residue']
738	RBData = rigidbodyDict
739	for item in parmDict:
740	if 'RB' in item:
741	dFdvDict[item] = 0. #NB: this is a vector which is no. refl. long & must be filled!
742	General = Phase['General']
743	cell = General['Cell'][1:7]
744	Amat,Bmat = G2lat.cell2AB(cell)
745	rpd = np.pi/180.
746	rpd2 = rpd**2
747	g = nl.inv(np.inner(Bmat,Bmat))
748	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
749	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
750	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
751	pfx = str(Phase['pId'])+'::'
752	RBModels = Phase['RBModels']
753	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
754	VModel = RBData['Vector'][RBObj['RBId']]
755	Q = RBObj['Orient'][0]
756	Pos = RBObj['Orig'][0]
757	jrb = VRBIds.index(RBObj['RBId'])
758	rbsx = str(irb)+':'+str(jrb)
759	dXdv = []
760	for iv in range(len(VModel['VectMag'])):
761	dCdv = []
762	for vec in VModel['rbVect'][iv]:
763	dCdv.append(G2mth.prodQVQ(Q,vec))
764	dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
765	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
766	for ia,atId in enumerate(RBObj['Ids']):
767	atNum = AtLookup[atId]
768	dx = 0.0001
769	for iv in range(len(VModel['VectMag'])):
770	for ix in [0,1,2]:
771	dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
772	for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
773	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
774	for iv in range(4):
775	Q[iv] -= dx
776	XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
777	Q[iv] += 2.*dx
778	XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
779	Q[iv] -= dx
780	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
781	for ix in [0,1,2]:
782	dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
783	X = G2mth.prodQVQ(Q,Cart[ia])
784	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
785	dFdu = G2lat.U6toUij(dFdu.T)
786	dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
787	dFdu = G2lat.UijtoU6(dFdu)
788	atNum = AtLookup[atId]
789	if 'T' in RBObj['ThermalMotion'][0]:
790	for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
791	dFdvDict[pfx+name+rbsx] += dFdu[i]
792	if 'L' in RBObj['ThermalMotion'][0]:
793	dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
794	dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
795	dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
796	dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
797	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
798	dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2(-dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
799	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
800	dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2(-dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
801	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
802	if 'S' in RBObj['ThermalMotion'][0]:
803	dFdvDict[pfx+'RBVS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
804	dFdvDict[pfx+'RBVS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
805	dFdvDict[pfx+'RBVS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
806	dFdvDict[pfx+'RBVS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
807	dFdvDict[pfx+'RBVS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
808	dFdvDict[pfx+'RBVS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
809	dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
810	dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
811	if 'U' in RBObj['ThermalMotion'][0]:
812	dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
813
814
815	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
816	Q = RBObj['Orient'][0]
817	Pos = RBObj['Orig'][0]
818	jrb = RRBIds.index(RBObj['RBId'])
819	torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
820	rbsx = str(irb)+':'+str(jrb)
821	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
822	for itors,tors in enumerate(RBObj['Torsions']): #derivative error?
823	tname = pfx+'RBRTr;'+str(itors)+':'+rbsx
824	orId,pvId = torData[itors][:2]
825	pivotVec = Cart[orId]-Cart[pvId]
826	QA = G2mth.AVdeg2Q(-0.001,pivotVec)
827	QB = G2mth.AVdeg2Q(0.001,pivotVec)
828	for ir in torData[itors][3]:
829	atNum = AtLookup[RBObj['Ids'][ir]]
830	rVec = Cart[ir]-Cart[pvId]
831	dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
832	dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
833	for ix in [0,1,2]:
834	dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
835	for ia,atId in enumerate(RBObj['Ids']):
836	atNum = AtLookup[atId]
837	dx = 0.0001
838	for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
839	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
840	for iv in range(4):
841	Q[iv] -= dx
842	XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
843	Q[iv] += 2.*dx
844	XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
845	Q[iv] -= dx
846	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
847	for ix in [0,1,2]:
848	dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
849	X = G2mth.prodQVQ(Q,Cart[ia])
850	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
851	dFdu = G2lat.U6toUij(dFdu.T)
852	dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
853	dFdu = G2lat.UijtoU6(dFdu)
854	atNum = AtLookup[atId]
855	if 'T' in RBObj['ThermalMotion'][0]:
856	for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
857	dFdvDict[pfx+name+rbsx] += dFdu[i]
858	if 'L' in RBObj['ThermalMotion'][0]:
859	dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
860	dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
861	dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
862	dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
863	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
864	dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2(dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
865	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
866	dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2(dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
867	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
868	if 'S' in RBObj['ThermalMotion'][0]:
869	dFdvDict[pfx+'RBRS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
870	dFdvDict[pfx+'RBRS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
871	dFdvDict[pfx+'RBRS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
872	dFdvDict[pfx+'RBRS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
873	dFdvDict[pfx+'RBRS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
874	dFdvDict[pfx+'RBRS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
875	dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
876	dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
877	if 'U' in RBObj['ThermalMotion'][0]:
878	dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
879
880
881	################################################################################
882	##### Phase data
883	################################################################################
884
885	def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None):
886
887	def PrintFFtable(FFtable):
888	print >>pFile,'\n X-ray scattering factors:'
889	print >>pFile,' Symbol fa fb fc'
890	print >>pFile,99*'-'
891	for Ename in FFtable:
892	ffdata = FFtable[Ename]
893	fa = ffdata['fa']
894	fb = ffdata['fb']
895	print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \
896	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
897
898	def PrintBLtable(BLtable):
899	print >>pFile,'\n Neutron scattering factors:'
900	print >>pFile,' Symbol isotope mass b resonant terms'
901	print >>pFile,99*'-'
902	for Ename in BLtable:
903	bldata = BLtable[Ename]
904	isotope = bldata[0]
905	mass = bldata[1][0]
906	blen = bldata[1][1]
907	bres = []
908	if len(bldata[1]) > 2:
909	bres = bldata[1][2:]
910	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
911	for item in bres:
912	line += '%10.5g'%(item)
913	print >>pFile,line
914
915	def PrintRBObjects(resRBData,vecRBData):
916
917	def PrintRBThermals():
918	tlstr = ['11','22','33','12','13','23']
919	sstr = ['12','13','21','23','31','32','AA','BB']
920	TLS = RB['ThermalMotion'][1]
921	TLSvar = RB['ThermalMotion'][2]
922	if 'T' in RB['ThermalMotion'][0]:
923	print >>pFile,'TLS data'
924	text = ''
925	for i in range(6):
926	text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
927	print >>pFile,text
928	if 'L' in RB['ThermalMotion'][0]:
929	text = ''
930	for i in range(6,12):
931	text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
932	print >>pFile,text
933	if 'S' in RB['ThermalMotion'][0]:
934	text = ''
935	for i in range(12,20):
936	text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
937	print >>pFile,text
938	if 'U' in RB['ThermalMotion'][0]:
939	print >>pFile,'Uiso data'
940	text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])
941
942	if len(resRBData):
943	for RB in resRBData:
944	Oxyz = RB['Orig'][0]
945	Qrijk = RB['Orient'][0]
946	Angle = 2.0*acosd(Qrijk[0])
947	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
948	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
949	print >>pFile,'Orientation angle,vector:', \
950	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
951	Torsions = RB['Torsions']
952	if len(Torsions):
953	text = 'Torsions: '
954	for torsion in Torsions:
955	text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
956	print >>pFile,text
957	PrintRBThermals()
958	if len(vecRBData):
959	for RB in vecRBData:
960	Oxyz = RB['Orig'][0]
961	Qrijk = RB['Orient'][0]
962	Angle = 2.0*acosd(Qrijk[0])
963	print >>pFile,'\nRBObject ',RB['RBname'],' at ', \
964	'%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
965	print >>pFile,'Orientation angle,vector:', \
966	'%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
967	PrintRBThermals()
968
969	def PrintAtoms(General,Atoms):
970	cx,ct,cs,cia = General['AtomPtrs']
971	print >>pFile,'\n Atoms:'
972	line = ' name type refine? x y z '+ \
973	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
974	if General['Type'] == 'magnetic':
975	line += ' Mx My Mz'
976	elif General['Type'] == 'macromolecular':
977	line = ' res no residue chain'+line
978	print >>pFile,line
979	if General['Type'] == 'nuclear':
980	print >>pFile,135*'-'
981	for i,at in enumerate(Atoms):
982	line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
983	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
984	if at[cia] == 'I':
985	line += '%8.4f'%(at[cia+1])+48*' '
986	else:
987	line += 8*' '
988	for j in range(6):
989	line += '%8.4f'%(at[cia+1+j])
990	print >>pFile,line
991	elif General['Type'] == 'macromolecular':
992	print >>pFile,135*'-'
993	for i,at in enumerate(Atoms):
994	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])+ \
995	'%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
996	if at[cia] == 'I':
997	line += '%8.4f'%(at[cia+1])+48*' '
998	else:
999	line += 8*' '
1000	for j in range(6):
1001	line += '%8.4f'%(at[cia+1+j])
1002	print >>pFile,line
1003
1004	def PrintTexture(textureData):
1005	topstr = '\n Spherical harmonics texture: Order:' + \
1006	str(textureData['Order'])
1007	if textureData['Order']:
1008	print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
1009	else:
1010	print >>pFile,topstr
1011	return
1012	names = ['omega','chi','phi']
1013	line = '\n'
1014	for name in names:
1015	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
1016	print >>pFile,line
1017	print >>pFile,'\n Texture coefficients:'
1018	ptlbls = ' names :'
1019	ptstr = ' values:'
1020	SHcoeff = textureData['SH Coeff'][1]
1021	for item in SHcoeff:
1022	ptlbls += '%12s'%(item)
1023	ptstr += '%12.4f'%(SHcoeff[item])
1024	print >>pFile,ptlbls
1025	print >>pFile,ptstr
1026
1027	def MakeRBParms(rbKey,phaseVary,phaseDict):
1028	rbid = str(rbids.index(RB['RBId']))
1029	pfxRB = pfx+'RB'+rbKey+'P'
1030	pstr = ['x','y','z']
1031	ostr = ['a','i','j','k']
1032	for i in range(3):
1033	name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1034	phaseDict[name] = RB['Orig'][0][i]
1035	if RB['Orig'][1]:
1036	phaseVary += [name,]
1037	pfxRB = pfx+'RB'+rbKey+'O'
1038	for i in range(4):
1039	name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1040	phaseDict[name] = RB['Orient'][0][i]
1041	if RB['Orient'][1] == 'AV' and i:
1042	phaseVary += [name,]
1043	elif RB['Orient'][1] == 'A' and not i:
1044	phaseVary += [name,]
1045
1046	def MakeRBThermals(rbKey,phaseVary,phaseDict):
1047	rbid = str(rbids.index(RB['RBId']))
1048	tlstr = ['11','22','33','12','13','23']
1049	sstr = ['12','13','21','23','31','32','AA','BB']
1050	if 'T' in RB['ThermalMotion'][0]:
1051	pfxRB = pfx+'RB'+rbKey+'T'
1052	for i in range(6):
1053	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1054	phaseDict[name] = RB['ThermalMotion'][1][i]
1055	if RB['ThermalMotion'][2][i]:
1056	phaseVary += [name,]
1057	if 'L' in RB['ThermalMotion'][0]:
1058	pfxRB = pfx+'RB'+rbKey+'L'
1059	for i in range(6):
1060	name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1061	phaseDict[name] = RB['ThermalMotion'][1][i+6]
1062	if RB['ThermalMotion'][2][i+6]:
1063	phaseVary += [name,]
1064	if 'S' in RB['ThermalMotion'][0]:
1065	pfxRB = pfx+'RB'+rbKey+'S'
1066	for i in range(8):
1067	name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1068	phaseDict[name] = RB['ThermalMotion'][1][i+12]
1069	if RB['ThermalMotion'][2][i+12]:
1070	phaseVary += [name,]
1071	if 'U' in RB['ThermalMotion'][0]:
1072	name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1073	phaseDict[name] = RB['ThermalMotion'][1][0]
1074	if RB['ThermalMotion'][2][0]:
1075	phaseVary += [name,]
1076
1077	def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1078	rbid = str(rbids.index(RB['RBId']))
1079	pfxRB = pfx+'RB'+rbKey+'Tr;'
1080	for i,torsion in enumerate(RB['Torsions']):
1081	name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1082	phaseDict[name] = torsion[0]
1083	if torsion[1]:
1084	phaseVary += [name,]
1085
1086	if Print:
1087	print >>pFile,'\n Phases:'
1088	phaseVary = []
1089	phaseDict = {}
1090	phaseConstr = {}
1091	pawleyLookup = {}
1092	FFtables = {} #scattering factors - xrays
1093	BLtables = {} # neutrons
1094	Natoms = {}
1095	AtMults = {}
1096	AtIA = {}
1097	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1098	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1099	atomIndx = {}
1100	for name in PhaseData:
1101	General = PhaseData[name]['General']
1102	pId = PhaseData[name]['pId']
1103	pfx = str(pId)+'::'
1104	FFtable = GetFFtable(General)
1105	BLtable = GetBLtable(General)
1106	FFtables.update(FFtable)
1107	BLtables.update(BLtable)
1108	Atoms = PhaseData[name]['Atoms']
1109	AtLookup = G2mth.FillAtomLookUp(Atoms)
1110	PawleyRef = PhaseData[name].get('Pawley ref',[])
1111	SGData = General['SGData']
1112	SGtext = G2spc.SGPrint(SGData)
1113	cell = General['Cell']
1114	A = G2lat.cell2A(cell[1:7])
1115	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1116	pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1117	if cell[0]:
1118	phaseVary += cellVary(pfx,SGData)
1119	resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1120	if resRBData:
1121	rbids = rbIds['Residue'] #NB: used in the MakeRB routines
1122	for iRB,RB in enumerate(resRBData):
1123	MakeRBParms('R',phaseVary,phaseDict)
1124	MakeRBThermals('R',phaseVary,phaseDict)
1125	MakeRBTorsions('R',phaseVary,phaseDict)
1126
1127	vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1128	if vecRBData:
1129	rbids = rbIds['Vector'] #NB: used in the MakeRB routines
1130	for iRB,RB in enumerate(vecRBData):
1131	MakeRBParms('V',phaseVary,phaseDict)
1132	MakeRBThermals('V',phaseVary,phaseDict)
1133
1134	Natoms[pfx] = 0
1135	if Atoms and not General.get('doPawley'):
1136	cx,ct,cs,cia = General['AtomPtrs']
1137	if General['Type'] in ['nuclear','macromolecular']:
1138	Natoms[pfx] = len(Atoms)
1139	for i,at in enumerate(Atoms):
1140	atomIndx[at[-1]] = [pfx,i] #lookup table for restraints
1141	phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1142	pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1143	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
1144	pfx+'AI/A:'+str(i):at[cia],})
1145	if at[cia] == 'I':
1146	phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1147	else:
1148	phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1149	pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1150	if 'F' in at[ct+1]:
1151	phaseVary.append(pfx+'Afrac:'+str(i))
1152	if 'X' in at[ct+1]:
1153	xId,xCoef = G2spc.GetCSxinel(at[cs])
1154	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1155	equivs = [[],[],[]]
1156	for j in range(3):
1157	if xId[j] > 0:
1158	phaseVary.append(names[j])
1159	equivs[xId[j]-1].append([names[j],xCoef[j]])
1160	for equiv in equivs:
1161	if len(equiv) > 1:
1162	name = equiv[0][0]
1163	for eqv in equiv[1:]:
1164	G2mv.StoreEquivalence(name,(eqv,))
1165	if 'U' in at[ct+1]:
1166	if at[9] == 'I':
1167	phaseVary.append(pfx+'AUiso:'+str(i))
1168	else:
1169	uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1170	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1171	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1172	equivs = [[],[],[],[],[],[]]
1173	for j in range(6):
1174	if uId[j] > 0:
1175	phaseVary.append(names[j])
1176	equivs[uId[j]-1].append([names[j],uCoef[j]])
1177	for equiv in equivs:
1178	if len(equiv) > 1:
1179	name = equiv[0][0]
1180	for eqv in equiv[1:]:
1181	G2mv.StoreEquivalence(name,(eqv,))
1182	# elif General['Type'] == 'magnetic':
1183	# elif General['Type'] == 'macromolecular':
1184	textureData = General['SH Texture']
1185	if textureData['Order']:
1186	phaseDict[pfx+'SHorder'] = textureData['Order']
1187	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1188	for item in ['omega','chi','phi']:
1189	phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1190	if textureData['Sample '+item][0]:
1191	phaseVary.append(pfx+'SH '+item)
1192	for item in textureData['SH Coeff'][1]:
1193	phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1194	if textureData['SH Coeff'][0]:
1195	phaseVary.append(pfx+item)
1196
1197	if Print:
1198	print >>pFile,'\n Phase name: ',General['Name']
1199	print >>pFile,135*'-'
1200	PrintFFtable(FFtable)
1201	PrintBLtable(BLtable)
1202	print >>pFile,''
1203	for line in SGtext: print >>pFile,line
1204	PrintRBObjects(resRBData,vecRBData)
1205	PrintAtoms(General,Atoms)
1206	print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
1207	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
1208	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
1209	PrintTexture(textureData)
1210	if name in RestraintDict:
1211	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1212	textureData,RestraintDict[name],pFile)
1213
1214	elif PawleyRef:
1215	pawleyVary = []
1216	for i,refl in enumerate(PawleyRef):
1217	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
1218	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1219	if refl[5]:
1220	pawleyVary.append(pfx+'PWLref:'+str(i))
1221	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
1222	phaseVary += pawleyVary
1223
1224	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1225
1226	def cellFill(pfx,SGData,parmDict,sigDict):
1227	if SGData['SGLaue'] in ['-1',]:
1228	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1229	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1230	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1231	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1232	elif SGData['SGLaue'] in ['2/m',]:
1233	if SGData['SGUniq'] == 'a':
1234	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1235	parmDict[pfx+'A3'],0,0]
1236	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1237	sigDict[pfx+'A3'],0,0]
1238	elif SGData['SGUniq'] == 'b':
1239	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1240	0,parmDict[pfx+'A4'],0]
1241	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1242	0,sigDict[pfx+'A4'],0]
1243	else:
1244	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1245	0,0,parmDict[pfx+'A5']]
1246	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1247	0,0,sigDict[pfx+'A5']]
1248	elif SGData['SGLaue'] in ['mmm',]:
1249	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1250	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1251	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1252	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1253	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1254	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1255	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1256	parmDict[pfx+'A0'],0,0]
1257	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1258	elif SGData['SGLaue'] in ['3R', '3mR']:
1259	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1260	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1261	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1262	elif SGData['SGLaue'] in ['m3m','m3']:
1263	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1264	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1265	return A,sigA
1266
1267	def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1268	if phaseRest:
1269	Amat = G2lat.cell2AB(cell)[0]
1270	cx,ct,cs = AtPtrs[:3]
1271	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1272	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1273	['ChemComp','Sites'],['Texture','HKLs']]
1274	for name,rest in names:
1275	itemRest = phaseRest[name]
1276	if itemRest[rest] and itemRest['Use']:
1277	print >>pFile,'\n %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1278	if name in ['Bond','Angle','Plane','Chiral']:
1279	print >>pFile,' calc obs sig delt/sig atoms(symOp): '
1280	for indx,ops,obs,esd in itemRest[rest]:
1281	try:
1282	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1283	AtName = ''
1284	for i,Aname in enumerate(AtNames):
1285	AtName += Aname
1286	if ops[i] == '1':
1287	AtName += '-'
1288	else:
1289	AtName += '+('+ops[i]+')-'
1290	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1291	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1292	if name == 'Bond':
1293	calc = G2mth.getRestDist(XYZ,Amat)
1294	elif name == 'Angle':
1295	calc = G2mth.getRestAngle(XYZ,Amat)
1296	elif name == 'Plane':
1297	calc = G2mth.getRestPlane(XYZ,Amat)
1298	elif name == 'Chiral':
1299	calc = G2mth.getRestChiral(XYZ,Amat)
1300	print >>pFile,' %9.3f %9.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1301	except KeyError:
1302	del itemRest[rest]
1303	elif name in ['Torsion','Rama']:
1304	print >>pFile,' atoms(symOp) calc obs sig delt/sig torsions: '
1305	coeffDict = itemRest['Coeff']
1306	for indx,ops,cofName,esd in itemRest[rest]:
1307	AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1308	AtName = ''
1309	for i,Aname in enumerate(AtNames):
1310	AtName += Aname+'+('+ops[i]+')-'
1311	XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1312	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1313	if name == 'Torsion':
1314	tor = G2mth.getRestTorsion(XYZ,Amat)
1315	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1316	print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1317	else:
1318	phi,psi = G2mth.getRestRama(XYZ,Amat)
1319	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
1320	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])
1321	elif name == 'ChemComp':
1322	print >>pFile,' atoms mul*frac factor prod'
1323	for indx,factors,obs,esd in itemRest[rest]:
1324	try:
1325	atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1326	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1327	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1328	mulfrac = mul*frac
1329	calcs = mulfracfactors
1330	for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1331	print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1332	print >>pFile,' Sum: calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1333	except KeyError:
1334	del itemRest[rest]
1335	elif name == 'Texture' and textureData['Order']:
1336	Start = False
1337	SHCoef = textureData['SH Coeff'][1]
1338	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1339	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1340	print ' HKL grid neg esd sum neg num neg use unit? unit esd '
1341	for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1342	PH = np.array(hkl)
1343	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1344	ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1345	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1346	Z = ma.masked_greater(Z,0.0)
1347	num = ma.count(Z)
1348	sum = 0
1349	if num:
1350	sum = np.sum(Z)
1351	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)
1352
1353	def getCellEsd(pfx,SGData,A,covData):
1354	dpr = 180./np.pi
1355	rVsq = G2lat.calc_rVsq(A)
1356	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
1357	cell = np.array(G2lat.Gmat2cell(g)) #real cell
1358	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1359	scos = cosd(cellst[3:6])
1360	ssin = sind(cellst[3:6])
1361	scot = scos/ssin
1362	rcos = cosd(cell[3:6])
1363	rsin = sind(cell[3:6])
1364	rcot = rcos/rsin
1365	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1366	varyList = covData['varyList']
1367	covMatrix = covData['covMatrix']
1368	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1369	Ax = np.array(A)
1370	Ax[3:] /= 2.
1371	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
1372	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
1373	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1374	Vol = 1/np.sqrt(rVsq)
1375	sigVol = Vol*3np.sqrt(srcvlsq)/2.
1376	R123 = Ax[0]Ax[1]Ax[2]
1377	dsasdg = np.zeros((3,6))
1378	dadg = np.zeros((6,6))
1379	for i0 in range(3): #0 1 2
1380	i1 = (i0+1)%3 #1 2 0
1381	i2 = (i1+1)%3 #2 0 1
1382	i3 = 5-i2 #3 5 4
1383	i4 = 5-i1 #4 3 5
1384	i5 = 5-i0 #5 4 3
1385	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
1386	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
1387	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1388	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
1389	denom = np.sqrt(denmsq)
1390	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
1391	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
1392	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
1393	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
1394	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
1395	dadg[i5][i5] = -Ax[i0]/denom
1396	for i0 in range(3):
1397	i1 = (i0+1)%3
1398	i2 = (i1+1)%3
1399	i3 = 5-i2
1400	for ij in range(6):
1401	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1402	if ij == i0:
1403	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1404	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
1405	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1406	var = np.diag(sigMat)
1407	CS = np.where(var>0.,np.sqrt(var),0.)
1408	cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
1409	return cellSig
1410
1411	def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1412
1413	def PrintAtomsAndSig(General,Atoms,atomsSig):
1414	print >>pFile,'\n Atoms:'
1415	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
1416	if General['Type'] == 'magnetic':
1417	line += ' Mx My Mz'
1418	elif General['Type'] == 'macromolecular':
1419	line = ' res no residue chain '+line
1420	print >>pFile,line
1421	if General['Type'] == 'nuclear':
1422	print >>pFile,135*'-'
1423	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1424	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1425	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
1426	for atyp in General['AtomTypes']: #zero composition
1427	General['NoAtoms'][atyp] = 0.
1428	for i,at in enumerate(Atoms):
1429	General['NoAtoms'][at[1]] += at[6]*float(at[8]) #new composition
1430	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1431	valstr = ' values:'
1432	sigstr = ' sig :'
1433	for ind in [3,4,5,6]:
1434	sigind = str(i)+':'+str(ind)
1435	valstr += fmt[ind]%(at[ind])
1436	if sigind in atomsSig:
1437	sigstr += fmt[ind]%(atomsSig[sigind])
1438	else:
1439	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1440	if at[9] == 'I':
1441	valstr += fmt[10]%(at[10])
1442	if str(i)+':10' in atomsSig:
1443	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1444	else:
1445	sigstr += 8*' '
1446	else:
1447	valstr += 8*' '
1448	sigstr += 8*' '
1449	for ind in [11,12,13,14,15,16]:
1450	sigind = str(i)+':'+str(ind)
1451	valstr += fmt[ind]%(at[ind])
1452	if sigind in atomsSig:
1453	sigstr += fmt[ind]%(atomsSig[sigind])
1454	else:
1455	sigstr += 8*' '
1456	print >>pFile,name
1457	print >>pFile,valstr
1458	print >>pFile,sigstr
1459
1460	def PrintRBObjPOAndSig(rbfx,rbsx):
1461	namstr = ' names :'
1462	valstr = ' values:'
1463	sigstr = ' esds :'
1464	for i,px in enumerate(['Px:','Py:','Pz:']):
1465	name = pfx+rbfx+px+rbsx
1466	namstr += '%12s'%('Pos '+px[1])
1467	valstr += '%12.5f'%(parmDict[name])
1468	if name in sigDict:
1469	sigstr += '%12.5f'%(sigDict[name])
1470	else:
1471	sigstr += 12*' '
1472	for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1473	name = pfx+rbfx+po+rbsx
1474	namstr += '%12s'%('Ori '+po[1])
1475	valstr += '%12.5f'%(parmDict[name])
1476	if name in sigDict:
1477	sigstr += '%12.5f'%(sigDict[name])
1478	else:
1479	sigstr += 12*' '
1480	print >>pFile,namstr
1481	print >>pFile,valstr
1482	print >>pFile,sigstr
1483
1484	def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1485	namstr = ' names :'
1486	valstr = ' values:'
1487	sigstr = ' esds :'
1488	if 'N' not in TLS:
1489	print >>pFile,' Thermal motion:'
1490	if 'T' in TLS:
1491	for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1492	name = pfx+rbfx+pt+rbsx
1493	namstr += '%12s'%(pt[:3])
1494	valstr += '%12.5f'%(parmDict[name])
1495	if name in sigDict:
1496	sigstr += '%12.5f'%(sigDict[name])
1497	else:
1498	sigstr += 12*' '
1499	print >>pFile,namstr
1500	print >>pFile,valstr
1501	print >>pFile,sigstr
1502	if 'L' in TLS:
1503	namstr = ' names :'
1504	valstr = ' values:'
1505	sigstr = ' esds :'
1506	for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1507	name = pfx+rbfx+pt+rbsx
1508	namstr += '%12s'%(pt[:3])
1509	valstr += '%12.3f'%(parmDict[name])
1510	if name in sigDict:
1511	sigstr += '%12.3f'%(sigDict[name])
1512	else:
1513	sigstr += 12*' '
1514	print >>pFile,namstr
1515	print >>pFile,valstr
1516	print >>pFile,sigstr
1517	if 'S' in TLS:
1518	namstr = ' names :'
1519	valstr = ' values:'
1520	sigstr = ' esds :'
1521	for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1522	name = pfx+rbfx+pt+rbsx
1523	namstr += '%12s'%(pt[:3])
1524	valstr += '%12.4f'%(parmDict[name])
1525	if name in sigDict:
1526	sigstr += '%12.4f'%(sigDict[name])
1527	else:
1528	sigstr += 12*' '
1529	print >>pFile,namstr
1530	print >>pFile,valstr
1531	print >>pFile,sigstr
1532	if 'U' in TLS:
1533	name = pfx+rbfx+'U:'+rbsx
1534	namstr = ' names :'+'%12s'%('Uiso')
1535	valstr = ' values:'+'%12.5f'%(parmDict[name])
1536	if name in sigDict:
1537	sigstr = ' esds :'+'%12.5f'%(sigDict[name])
1538	else:
1539	sigstr = ' esds :'+12*' '
1540	print >>pFile,namstr
1541	print >>pFile,valstr
1542	print >>pFile,sigstr
1543
1544	def PrintRBObjTorAndSig(rbsx):
1545	namstr = ' names :'
1546	valstr = ' values:'
1547	sigstr = ' esds :'
1548	nTors = len(RBObj['Torsions'])
1549	if nTors:
1550	print >>pFile,' Torsions:'
1551	for it in range(nTors):
1552	name = pfx+'RBRTr;'+str(it)+':'+rbsx
1553	namstr += '%12s'%('Tor'+str(it))
1554	valstr += '%12.4f'%(parmDict[name])
1555	if name in sigDict:
1556	sigstr += '%12.4f'%(sigDict[name])
1557	print >>pFile,namstr
1558	print >>pFile,valstr
1559	print >>pFile,sigstr
1560
1561	def PrintSHtextureAndSig(textureData,SHtextureSig):
1562	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1563	names = ['omega','chi','phi']
1564	namstr = ' names :'
1565	ptstr = ' values:'
1566	sigstr = ' esds :'
1567	for name in names:
1568	namstr += '%12s'%(name)
1569	ptstr += '%12.3f'%(textureData['Sample '+name][1])
1570	if 'Sample '+name in SHtextureSig:
1571	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1572	else:
1573	sigstr += 12*' '
1574	print >>pFile,namstr
1575	print >>pFile,ptstr
1576	print >>pFile,sigstr
1577	print >>pFile,'\n Texture coefficients:'
1578	namstr = ' names :'
1579	ptstr = ' values:'
1580	sigstr = ' esds :'
1581	SHcoeff = textureData['SH Coeff'][1]
1582	for name in SHcoeff:
1583	namstr += '%12s'%(name)
1584	ptstr += '%12.3f'%(SHcoeff[name])
1585	if name in SHtextureSig:
1586	sigstr += '%12.3f'%(SHtextureSig[name])
1587	else:
1588	sigstr += 12*' '
1589	print >>pFile,namstr
1590	print >>pFile,ptstr
1591	print >>pFile,sigstr
1592
1593	print >>pFile,'\n Phases:'
1594	for phase in Phases:
1595	print >>pFile,' Result for phase: ',phase
1596	Phase = Phases[phase]
1597	General = Phase['General']
1598	SGData = General['SGData']
1599	Atoms = Phase['Atoms']
1600	AtLookup = G2mth.FillAtomLookUp(Atoms)
1601	cell = General['Cell']
1602	pId = Phase['pId']
1603	pfx = str(pId)+'::'
1604	if cell[0]:
1605	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1606	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
1607	print >>pFile,' Reciprocal metric tensor: '
1608	ptfmt = "%15.9f"
1609	names = ['A11','A22','A33','A12','A13','A23']
1610	namstr = ' names :'
1611	ptstr = ' values:'
1612	sigstr = ' esds :'
1613	for name,a,siga in zip(names,A,sigA):
1614	namstr += '%15s'%(name)
1615	ptstr += ptfmt%(a)
1616	if siga:
1617	sigstr += ptfmt%(siga)
1618	else:
1619	sigstr += 15*' '
1620	print >>pFile,namstr
1621	print >>pFile,ptstr
1622	print >>pFile,sigstr
1623	cell[1:7] = G2lat.A2cell(A)
1624	cell[7] = G2lat.calc_V(A)
1625	print >>pFile,' New unit cell:'
1626	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1627	names = ['a','b','c','alpha','beta','gamma','Volume']
1628	namstr = ' names :'
1629	ptstr = ' values:'
1630	sigstr = ' esds :'
1631	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1632	namstr += '%12s'%(name)
1633	ptstr += fmt%(a)
1634	if siga:
1635	sigstr += fmt%(siga)
1636	else:
1637	sigstr += 12*' '
1638	print >>pFile,namstr
1639	print >>pFile,ptstr
1640	print >>pFile,sigstr
1641
1642	General['Mass'] = 0.
1643	if Phase['General'].get('doPawley'):
1644	pawleyRef = Phase['Pawley ref']
1645	for i,refl in enumerate(pawleyRef):
1646	key = pfx+'PWLref:'+str(i)
1647	refl[6] = parmDict[key]
1648	if key in sigDict:
1649	refl[7] = sigDict[key]
1650	else:
1651	refl[7] = 0
1652	else:
1653	VRBIds = RBIds['Vector']
1654	RRBIds = RBIds['Residue']
1655	RBModels = Phase['RBModels']
1656	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1657	jrb = VRBIds.index(RBObj['RBId'])
1658	rbsx = str(irb)+':'+str(jrb)
1659	print >>pFile,' Vector rigid body parameters:'
1660	PrintRBObjPOAndSig('RBV',rbsx)
1661	PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1662	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1663	jrb = RRBIds.index(RBObj['RBId'])
1664	rbsx = str(irb)+':'+str(jrb)
1665	print >>pFile,' Residue rigid body parameters:'
1666	PrintRBObjPOAndSig('RBR',rbsx)
1667	PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1668	PrintRBObjTorAndSig(rbsx)
1669	atomsSig = {}
1670	if General['Type'] == 'nuclear': #this needs macromolecular variant!
1671	for i,at in enumerate(Atoms):
1672	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1673	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1674	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1675	for ind in [3,4,5,6]:
1676	at[ind] = parmDict[names[ind]]
1677	if ind in [3,4,5]:
1678	name = names[ind].replace('A','dA')
1679	else:
1680	name = names[ind]
1681	if name in sigDict:
1682	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1683	if at[9] == 'I':
1684	at[10] = parmDict[names[10]]
1685	if names[10] in sigDict:
1686	atomsSig[str(i)+':10'] = sigDict[names[10]]
1687	else:
1688	for ind in [11,12,13,14,15,16]:
1689	at[ind] = parmDict[names[ind]]
1690	if names[ind] in sigDict:
1691	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1692	ind = General['AtomTypes'].index(at[1])
1693	General['Mass'] += General['AtomMass'][ind]at[6]at[8]
1694	PrintAtomsAndSig(General,Atoms,atomsSig)
1695
1696	textureData = General['SH Texture']
1697	if textureData['Order']:
1698	SHtextureSig = {}
1699	for name in ['omega','chi','phi']:
1700	aname = pfx+'SH '+name
1701	textureData['Sample '+name][1] = parmDict[aname]
1702	if aname in sigDict:
1703	SHtextureSig['Sample '+name] = sigDict[aname]
1704	for name in textureData['SH Coeff'][1]:
1705	aname = pfx+name
1706	textureData['SH Coeff'][1][name] = parmDict[aname]
1707	if aname in sigDict:
1708	SHtextureSig[name] = sigDict[aname]
1709	PrintSHtextureAndSig(textureData,SHtextureSig)
1710	if phase in RestraintDict:
1711	PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1712	textureData,RestraintDict[phase],pFile)
1713
1714	################################################################################
1715	##### Histogram & Phase data
1716	################################################################################
1717
1718	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None):
1719
1720	def PrintSize(hapData):
1721	if hapData[0] in ['isotropic','uniaxial']:
1722	line = '\n Size model : %9s'%(hapData[0])
1723	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1724	if hapData[0] == 'uniaxial':
1725	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1726	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1727	print >>pFile,line
1728	else:
1729	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1730	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1731	Snames = ['S11','S22','S33','S12','S13','S23']
1732	ptlbls = ' names :'
1733	ptstr = ' values:'
1734	varstr = ' refine:'
1735	for i,name in enumerate(Snames):
1736	ptlbls += '%12s' % (name)
1737	ptstr += '%12.6f' % (hapData[4][i])
1738	varstr += '%12s' % (str(hapData[5][i]))
1739	print >>pFile,ptlbls
1740	print >>pFile,ptstr
1741	print >>pFile,varstr
1742
1743	def PrintMuStrain(hapData,SGData):
1744	if hapData[0] in ['isotropic','uniaxial']:
1745	line = '\n Mustrain model: %9s'%(hapData[0])
1746	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1747	if hapData[0] == 'uniaxial':
1748	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1749	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1750	print >>pFile,line
1751	else:
1752	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1753	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1754	Snames = G2spc.MustrainNames(SGData)
1755	ptlbls = ' names :'
1756	ptstr = ' values:'
1757	varstr = ' refine:'
1758	for i,name in enumerate(Snames):
1759	ptlbls += '%12s' % (name)
1760	ptstr += '%12.6f' % (hapData[4][i])
1761	varstr += '%12s' % (str(hapData[5][i]))
1762	print >>pFile,ptlbls
1763	print >>pFile,ptstr
1764	print >>pFile,varstr
1765
1766	def PrintHStrain(hapData,SGData):
1767	print >>pFile,'\n Hydrostatic/elastic strain: '
1768	Hsnames = G2spc.HStrainNames(SGData)
1769	ptlbls = ' names :'
1770	ptstr = ' values:'
1771	varstr = ' refine:'
1772	for i,name in enumerate(Hsnames):
1773	ptlbls += '%12s' % (name)
1774	ptstr += '%12.6f' % (hapData[0][i])
1775	varstr += '%12s' % (str(hapData[1][i]))
1776	print >>pFile,ptlbls
1777	print >>pFile,ptstr
1778	print >>pFile,varstr
1779
1780	def PrintSHPO(hapData):
1781	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1782	str(hapData[4])+' Refine? '+str(hapData[2])
1783	ptlbls = ' names :'
1784	ptstr = ' values:'
1785	for item in hapData[5]:
1786	ptlbls += '%12s'%(item)
1787	ptstr += '%12.3f'%(hapData[5][item])
1788	print >>pFile,ptlbls
1789	print >>pFile,ptstr
1790
1791	def PrintBabinet(hapData):
1792	print >>pFile,'\n Babinet form factor modification: '
1793	ptlbls = ' names :'
1794	ptstr = ' values:'
1795	varstr = ' refine:'
1796	for name in ['BabA','BabU']:
1797	ptlbls += '%12s' % (name)
1798	ptstr += '%12.6f' % (hapData[name][0])
1799	varstr += '%12s' % (str(hapData[name][1]))
1800	print >>pFile,ptlbls
1801	print >>pFile,ptstr
1802	print >>pFile,varstr
1803
1804
1805	hapDict = {}
1806	hapVary = []
1807	controlDict = {}
1808	poType = {}
1809	poAxes = {}
1810	spAxes = {}
1811	spType = {}
1812
1813	for phase in Phases:
1814	HistoPhase = Phases[phase]['Histograms']
1815	SGData = Phases[phase]['General']['SGData']
1816	cell = Phases[phase]['General']['Cell'][1:7]
1817	A = G2lat.cell2A(cell)
1818	pId = Phases[phase]['pId']
1819	histoList = HistoPhase.keys()
1820	histoList.sort()
1821	for histogram in histoList:
1822	try:
1823	Histogram = Histograms[histogram]
1824	except KeyError:
1825	#skip if histogram not included e.g. in a sequential refinement
1826	continue
1827	hapData = HistoPhase[histogram]
1828	hId = Histogram['hId']
1829	if 'PWDR' in histogram:
1830	limits = Histogram['Limits'][1]
1831	inst = Histogram['Instrument Parameters'][0]
1832	Zero = inst['Zero'][1]
1833	if 'C' in inst['Type'][1]:
1834	try:
1835	wave = inst['Lam'][1]
1836	except KeyError:
1837	wave = inst['Lam1'][1]
1838	dmin = wave/(2.0*sind(limits[1]/2.0))
1839	pfx = str(pId)+':'+str(hId)+':'
1840	for item in ['Scale','Extinction']:
1841	hapDict[pfx+item] = hapData[item][0]
1842	if hapData[item][1]:
1843	hapVary.append(pfx+item)
1844	names = G2spc.HStrainNames(SGData)
1845	for i,name in enumerate(names):
1846	hapDict[pfx+name] = hapData['HStrain'][0][i]
1847	if hapData['HStrain'][1][i]:
1848	hapVary.append(pfx+name)
1849	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1850	if hapData['Pref.Ori.'][0] == 'MD':
1851	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1852	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1853	if hapData['Pref.Ori.'][2]:
1854	hapVary.append(pfx+'MD')
1855	else: #'SH' spherical harmonics
1856	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1857	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1858	for item in hapData['Pref.Ori.'][5]:
1859	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1860	if hapData['Pref.Ori.'][2]:
1861	hapVary.append(pfx+item)
1862	for item in ['Mustrain','Size']:
1863	controlDict[pfx+item+'Type'] = hapData[item][0]
1864	hapDict[pfx+item+';mx'] = hapData[item][1][2]
1865	if hapData[item][2][2]:
1866	hapVary.append(pfx+item+';mx')
1867	if hapData[item][0] in ['isotropic','uniaxial']:
1868	hapDict[pfx+item+';i'] = hapData[item][1][0]
1869	if hapData[item][2][0]:
1870	hapVary.append(pfx+item+';i')
1871	if hapData[item][0] == 'uniaxial':
1872	controlDict[pfx+item+'Axis'] = hapData[item][3]
1873	hapDict[pfx+item+';a'] = hapData[item][1][1]
1874	if hapData[item][2][1]:
1875	hapVary.append(pfx+item+';a')
1876	else: #generalized for mustrain or ellipsoidal for size
1877	Nterms = len(hapData[item][4])
1878	if item == 'Mustrain':
1879	names = G2spc.MustrainNames(SGData)
1880	pwrs = []
1881	for name in names:
1882	h,k,l = name[1:]
1883	pwrs.append([int(h),int(k),int(l)])
1884	controlDict[pfx+'MuPwrs'] = pwrs
1885	for i in range(Nterms):
1886	sfx = ':'+str(i)
1887	hapDict[pfx+item+sfx] = hapData[item][4][i]
1888	if hapData[item][5][i]:
1889	hapVary.append(pfx+item+sfx)
1890	for bab in ['BabA','BabU']:
1891	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1892	if hapData['Babinet'][bab][1]:
1893	hapVary.append(pfx+bab)
1894
1895	if Print:
1896	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1897	print >>pFile,135*'-'
1898	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1899	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1900	if hapData['Pref.Ori.'][0] == 'MD':
1901	Ax = hapData['Pref.Ori.'][3]
1902	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1903	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1904	else: #'SH' for spherical harmonics
1905	PrintSHPO(hapData['Pref.Ori.'])
1906	PrintSize(hapData['Size'])
1907	PrintMuStrain(hapData['Mustrain'],SGData)
1908	PrintHStrain(hapData['HStrain'],SGData)
1909	if hapData['Babinet']['BabA'][0]:
1910	PrintBabinet(hapData['Babinet'])
1911	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1912	refList = []
1913	for h,k,l,d in HKLd:
1914	ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1915	mul *= 2 # for powder overlap of Friedel pairs
1916	if ext:
1917	continue
1918	if 'C' in inst['Type'][0]:
1919	pos = 2.0asind(wave/(2.0d))+Zero
1920	if limits[0] < pos < limits[1]:
1921	refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1922	#last item should contain form factor values by atom type
1923	else:
1924	raise ValueError
1925	Histogram['Reflection Lists'][phase] = refList
1926	elif 'HKLF' in histogram:
1927	inst = Histogram['Instrument Parameters'][0]
1928	hId = Histogram['hId']
1929	hfx = ':%d:'%(hId)
1930	for item in inst:
1931	hapDict[hfx+item] = inst[item][1]
1932	pfx = str(pId)+':'+str(hId)+':'
1933	hapDict[pfx+'Scale'] = hapData['Scale'][0]
1934	if hapData['Scale'][1]:
1935	hapVary.append(pfx+'Scale')
1936
1937	extApprox,extType,extParms = hapData['Extinction']
1938	controlDict[pfx+'EType'] = extType
1939	controlDict[pfx+'EApprox'] = extApprox
1940	controlDict[pfx+'Tbar'] = extParms['Tbar']
1941	controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1942	if 'Primary' in extType:
1943	Ekey = ['Ep',]
1944	elif 'I & II' in extType:
1945	Ekey = ['Eg','Es']
1946	elif 'Secondary Type II' == extType:
1947	Ekey = ['Es',]
1948	elif 'Secondary Type I' == extType:
1949	Ekey = ['Eg',]
1950	else: #'None'
1951	Ekey = []
1952	for eKey in Ekey:
1953	hapDict[pfx+eKey] = extParms[eKey][0]
1954	if extParms[eKey][1]:
1955	hapVary.append(pfx+eKey)
1956	for bab in ['BabA','BabU']:
1957	hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1958	if hapData['Babinet'][bab][1]:
1959	hapVary.append(pfx+bab)
1960	if Print:
1961	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1962	print >>pFile,135*'-'
1963	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1964	if extType != 'None':
1965	print >>pFile,' Extinction Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1966	text = ' Parameters :'
1967	for eKey in Ekey:
1968	text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1969	print >>pFile,text
1970	if hapData['Babinet']['BabA'][0]:
1971	PrintBabinet(hapData['Babinet'])
1972	Histogram['Reflection Lists'] = phase
1973
1974	return hapVary,hapDict,controlDict
1975
1976	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1977
1978	def PrintSizeAndSig(hapData,sizeSig):
1979	line = '\n Size model: %9s'%(hapData[0])
1980	refine = False
1981	if hapData[0] in ['isotropic','uniaxial']:
1982	line += ' equatorial:%12.4f'%(hapData[1][0])
1983	if sizeSig[0][0]:
1984	line += ', sig:%8.4f'%(sizeSig[0][0])
1985	refine = True
1986	if hapData[0] == 'uniaxial':
1987	line += ' axial:%12.4f'%(hapData[1][1])
1988	if sizeSig[0][1]:
1989	refine = True
1990	line += ', sig:%8.4f'%(sizeSig[0][1])
1991	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1992	if sizeSig[0][2]:
1993	refine = True
1994	line += ', sig:%8.4f'%(sizeSig[0][2])
1995	if refine:
1996	print >>pFile,line
1997	else:
1998	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1999	if sizeSig[0][2]:
2000	refine = True
2001	line += ', sig:%8.4f'%(sizeSig[0][2])
2002	Snames = ['S11','S22','S33','S12','S13','S23']
2003	ptlbls = ' name :'
2004	ptstr = ' value :'
2005	sigstr = ' sig :'
2006	for i,name in enumerate(Snames):
2007	ptlbls += '%12s' % (name)
2008	ptstr += '%12.6f' % (hapData[4][i])
2009	if sizeSig[1][i]:
2010	refine = True
2011	sigstr += '%12.6f' % (sizeSig[1][i])
2012	else:
2013	sigstr += 12*' '
2014	if refine:
2015	print >>pFile,line
2016	print >>pFile,ptlbls
2017	print >>pFile,ptstr
2018	print >>pFile,sigstr
2019
2020	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2021	line = '\n Mustrain model: %9s'%(hapData[0])
2022	refine = False
2023	if hapData[0] in ['isotropic','uniaxial']:
2024	line += ' equatorial:%12.1f'%(hapData[1][0])
2025	if mustrainSig[0][0]:
2026	line += ', sig:%8.1f'%(mustrainSig[0][0])
2027	refine = True
2028	if hapData[0] == 'uniaxial':
2029	line += ' axial:%12.1f'%(hapData[1][1])
2030	if mustrainSig[0][1]:
2031	line += ', sig:%8.1f'%(mustrainSig[0][1])
2032	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2033	if mustrainSig[0][2]:
2034	refine = True
2035	line += ', sig:%8.3f'%(mustrainSig[0][2])
2036	if refine:
2037	print >>pFile,line
2038	else:
2039	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2040	if mustrainSig[0][2]:
2041	refine = True
2042	line += ', sig:%8.3f'%(mustrainSig[0][2])
2043	Snames = G2spc.MustrainNames(SGData)
2044	ptlbls = ' name :'
2045	ptstr = ' value :'
2046	sigstr = ' sig :'
2047	for i,name in enumerate(Snames):
2048	ptlbls += '%12s' % (name)
2049	ptstr += '%12.6f' % (hapData[4][i])
2050	if mustrainSig[1][i]:
2051	refine = True
2052	sigstr += '%12.6f' % (mustrainSig[1][i])
2053	else:
2054	sigstr += 12*' '
2055	if refine:
2056	print >>pFile,line
2057	print >>pFile,ptlbls
2058	print >>pFile,ptstr
2059	print >>pFile,sigstr
2060
2061	def PrintHStrainAndSig(hapData,strainSig,SGData):
2062	Hsnames = G2spc.HStrainNames(SGData)
2063	ptlbls = ' name :'
2064	ptstr = ' value :'
2065	sigstr = ' sig :'
2066	refine = False
2067	for i,name in enumerate(Hsnames):
2068	ptlbls += '%12s' % (name)
2069	ptstr += '%12.6g' % (hapData[0][i])
2070	if name in strainSig:
2071	refine = True
2072	sigstr += '%12.6g' % (strainSig[name])
2073	else:
2074	sigstr += 12*' '
2075	if refine:
2076	print >>pFile,'\n Hydrostatic/elastic strain: '
2077	print >>pFile,ptlbls
2078	print >>pFile,ptstr
2079	print >>pFile,sigstr
2080
2081	def PrintSHPOAndSig(pfx,hapData,POsig):
2082	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2083	ptlbls = ' names :'
2084	ptstr = ' values:'
2085	sigstr = ' sig :'
2086	for item in hapData[5]:
2087	ptlbls += '%12s'%(item)
2088	ptstr += '%12.3f'%(hapData[5][item])
2089	if pfx+item in POsig:
2090	sigstr += '%12.3f'%(POsig[pfx+item])
2091	else:
2092	sigstr += 12*' '
2093	print >>pFile,ptlbls
2094	print >>pFile,ptstr
2095	print >>pFile,sigstr
2096
2097	def PrintExtAndSig(pfx,hapData,ScalExtSig):
2098	print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2099	text = ''
2100	for item in hapData[2]:
2101	if pfx+item in ScalExtSig:
2102	text += ' %s: '%(item)
2103	text += '%12.2e'%(hapData[2][item][0])
2104	if pfx+item in ScalExtSig:
2105	text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2106	print >>pFile,text
2107
2108	def PrintBabinetAndSig(pfx,hapData,BabSig):
2109	print >>pFile,'\n Babinet form factor modification: '
2110	ptlbls = ' names :'
2111	ptstr = ' values:'
2112	sigstr = ' sig :'
2113	for item in hapData:
2114	ptlbls += '%12s'%(item)
2115	ptstr += '%12.3f'%(hapData[item][0])
2116	if pfx+item in BabSig:
2117	sigstr += '%12.3f'%(BabSig[pfx+item])
2118	else:
2119	sigstr += 12*' '
2120	print >>pFile,ptlbls
2121	print >>pFile,ptstr
2122	print >>pFile,sigstr
2123
2124	PhFrExtPOSig = {}
2125	SizeMuStrSig = {}
2126	ScalExtSig = {}
2127	BabSig = {}
2128	wtFrSum = {}
2129	for phase in Phases:
2130	HistoPhase = Phases[phase]['Histograms']
2131	General = Phases[phase]['General']
2132	SGData = General['SGData']
2133	pId = Phases[phase]['pId']
2134	histoList = HistoPhase.keys()
2135	histoList.sort()
2136	for histogram in histoList:
2137	try:
2138	Histogram = Histograms[histogram]
2139	except KeyError:
2140	#skip if histogram not included e.g. in a sequential refinement
2141	continue
2142	hapData = HistoPhase[histogram]
2143	hId = Histogram['hId']
2144	pfx = str(pId)+':'+str(hId)+':'
2145	if hId not in wtFrSum:
2146	wtFrSum[hId] = 0.
2147	if 'PWDR' in histogram:
2148	for item in ['Scale','Extinction']:
2149	hapData[item][0] = parmDict[pfx+item]
2150	if pfx+item in sigDict:
2151	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2152	wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2153	if hapData['Pref.Ori.'][0] == 'MD':
2154	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2155	if pfx+'MD' in sigDict:
2156	PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2157	else: #'SH' spherical harmonics
2158	for item in hapData['Pref.Ori.'][5]:
2159	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2160	if pfx+item in sigDict:
2161	PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2162	SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2163	pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2164	pfx+'HStrain':{}})
2165	for item in ['Mustrain','Size']:
2166	hapData[item][1][2] = parmDict[pfx+item+';mx']
2167	hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
2168	if pfx+item+';mx' in sigDict:
2169	SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2170	if hapData[item][0] in ['isotropic','uniaxial']:
2171	hapData[item][1][0] = parmDict[pfx+item+';i']
2172	if item == 'Size':
2173	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2174	if pfx+item+';i' in sigDict:
2175	SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2176	if hapData[item][0] == 'uniaxial':
2177	hapData[item][1][1] = parmDict[pfx+item+';a']
2178	if item == 'Size':
2179	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
2180	if pfx+item+';a' in sigDict:
2181	SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2182	else: #generalized for mustrain or ellipsoidal for size
2183	Nterms = len(hapData[item][4])
2184	for i in range(Nterms):
2185	sfx = ':'+str(i)
2186	hapData[item][4][i] = parmDict[pfx+item+sfx]
2187	if pfx+item+sfx in sigDict:
2188	SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2189	names = G2spc.HStrainNames(SGData)
2190	for i,name in enumerate(names):
2191	hapData['HStrain'][0][i] = parmDict[pfx+name]
2192	if pfx+name in sigDict:
2193	SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2194	for name in ['BabA','BabU']:
2195	hapData['Babinet'][name][0] = parmDict[pfx+name]
2196	if pfx+name in sigDict:
2197	BabSig[pfx+name] = sigDict[pfx+name]
2198
2199	elif 'HKLF' in histogram:
2200	for item in ['Scale',]:
2201	if parmDict.get(pfx+item):
2202	hapData[item][0] = parmDict[pfx+item]
2203	if pfx+item in sigDict:
2204	ScalExtSig[pfx+item] = sigDict[pfx+item]
2205	for item in ['Ep','Eg','Es']:
2206	if parmDict.get(pfx+item):
2207	hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2208	if pfx+item in sigDict:
2209	ScalExtSig[pfx+item] = sigDict[pfx+item]
2210	for name in ['BabA','BabU']:
2211	hapData['Babinet'][name][0] = parmDict[pfx+name]
2212	if pfx+name in sigDict:
2213	BabSig[pfx+name] = sigDict[pfx+name]
2214
2215	if Print:
2216	for phase in Phases:
2217	HistoPhase = Phases[phase]['Histograms']
2218	General = Phases[phase]['General']
2219	SGData = General['SGData']
2220	pId = Phases[phase]['pId']
2221	histoList = HistoPhase.keys()
2222	histoList.sort()
2223	for histogram in histoList:
2224	try:
2225	Histogram = Histograms[histogram]
2226	except KeyError:
2227	#skip if histogram not included e.g. in a sequential refinement
2228	continue
2229	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2230	print >>pFile,130*'-'
2231	hapData = HistoPhase[histogram]
2232	hId = Histogram['hId']
2233	pfx = str(pId)+':'+str(hId)+':'
2234	if 'PWDR' in histogram:
2235	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2236	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2237
2238	if pfx+'Scale' in PhFrExtPOSig:
2239	wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2240	sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2241	print >>pFile,' Phase fraction : %10.5f, sig %10.5f Weight fraction : %8.5f, sig %10.5f' \
2242	%(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2243	if pfx+'Extinction' in PhFrExtPOSig:
2244	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2245	if hapData['Pref.Ori.'][0] == 'MD':
2246	if pfx+'MD' in PhFrExtPOSig:
2247	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2248	else:
2249	PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2250	PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2251	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2252	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2253	if len(BabSig):
2254	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2255
2256	elif 'HKLF' in histogram:
2257	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
2258	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2259	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2260	if pfx+'Scale' in ScalExtSig:
2261	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2262	if hapData['Extinction'][0] != 'None':
2263	PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2264	if len(BabSig):
2265	PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2266
2267	################################################################################
2268	##### Histogram data
2269	################################################################################
2270
2271	def GetHistogramData(Histograms,Print=True,pFile=None):
2272
2273	def GetBackgroundParms(hId,Background):
2274	Back = Background[0]
2275	DebyePeaks = Background[1]
2276	bakType,bakFlag = Back[:2]
2277	backVals = Back[3:]
2278	backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2279	backDict = dict(zip(backNames,backVals))
2280	backVary = []
2281	if bakFlag:
2282	backVary = backNames
2283	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2284	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2285	debyeDict = {}
2286	debyeList = []
2287	for i in range(DebyePeaks['nDebye']):
2288	debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2289	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2290	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2291	debyeVary = []
2292	for item in debyeList:
2293	if item[1]:
2294	debyeVary.append(item[0])
2295	backDict.update(debyeDict)
2296	backVary += debyeVary
2297	peakDict = {}
2298	peakList = []
2299	for i in range(DebyePeaks['nPeaks']):
2300	peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
2301	':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
2302	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2303	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2304	peakVary = []
2305	for item in peakList:
2306	if item[1]:
2307	peakVary.append(item[0])
2308	backDict.update(peakDict)
2309	backVary += peakVary
2310	return bakType,backDict,backVary
2311
2312	def GetInstParms(hId,Inst):
2313	dataType = Inst['Type'][0]
2314	instDict = {}
2315	insVary = []
2316	pfx = ':'+str(hId)+':'
2317	insKeys = Inst.keys()
2318	insKeys.sort()
2319	for item in insKeys:
2320	insName = pfx+item
2321	instDict[insName] = Inst[item][1]
2322	if Inst[item][2]:
2323	insVary.append(insName)
2324	# instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2325	# instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2326	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2327	return dataType,instDict,insVary
2328
2329	def GetSampleParms(hId,Sample):
2330	sampVary = []
2331	hfx = ':'+str(hId)+':'
2332	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2333	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2334	Type = Sample['Type']
2335	if 'Bragg' in Type: #Bragg-Brentano
2336	for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering?
2337	sampDict[hfx+item] = Sample[item][0]
2338	if Sample[item][1]:
2339	sampVary.append(hfx+item)
2340	elif 'Debye' in Type: #Debye-Scherrer
2341	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2342	sampDict[hfx+item] = Sample[item][0]
2343	if Sample[item][1]:
2344	sampVary.append(hfx+item)
2345	return Type,sampDict,sampVary
2346
2347	def PrintBackground(Background):
2348	Back = Background[0]
2349	DebyePeaks = Background[1]
2350	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2351	line = ' Coefficients: '
2352	for i,back in enumerate(Back[3:]):
2353	line += '%10.3f'%(back)
2354	if i and not i%10:
2355	line += '\n'+15*' '
2356	print >>pFile,line
2357	if DebyePeaks['nDebye']:
2358	print >>pFile,'\n Debye diffuse scattering coefficients'
2359	parms = ['DebyeA','DebyeR','DebyeU']
2360	line = ' names : '
2361	for parm in parms:
2362	line += '%8s refine?'%(parm)
2363	print >>pFile,line
2364	for j,term in enumerate(DebyePeaks['debyeTerms']):
2365	line = ' term'+'%2d'%(j)+':'
2366	for i in range(3):
2367	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2368	print >>pFile,line
2369	if DebyePeaks['nPeaks']:
2370	print >>pFile,'\n Single peak coefficients'
2371	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2372	line = ' names : '
2373	for parm in parms:
2374	line += '%8s refine?'%(parm)
2375	print >>pFile,line
2376	for j,term in enumerate(DebyePeaks['peaksList']):
2377	line = ' peak'+'%2d'%(j)+':'
2378	for i in range(4):
2379	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
2380	print >>pFile,line
2381
2382	def PrintInstParms(Inst):
2383	print >>pFile,'\n Instrument Parameters:'
2384	ptlbls = ' name :'
2385	ptstr = ' value :'
2386	varstr = ' refine:'
2387	insKeys = Inst.keys()
2388	insKeys.sort()
2389	for item in insKeys:
2390	if item != 'Type':
2391	ptlbls += '%12s' % (item)
2392	ptstr += '%12.6f' % (Inst[item][1])
2393	if item in ['Lam1','Lam2','Azimuth']:
2394	varstr += 12*' '
2395	else:
2396	varstr += '%12s' % (str(bool(Inst[item][2])))
2397	print >>pFile,ptlbls
2398	print >>pFile,ptstr
2399	print >>pFile,varstr
2400
2401	def PrintSampleParms(Sample):
2402	print >>pFile,'\n Sample Parameters:'
2403	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2404	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
2405	ptlbls = ' name :'
2406	ptstr = ' value :'
2407	varstr = ' refine:'
2408	if 'Bragg' in Sample['Type']:
2409	for item in ['Scale','Shift','Transparency']:
2410	ptlbls += '%14s'%(item)
2411	ptstr += '%14.4f'%(Sample[item][0])
2412	varstr += '%14s'%(str(bool(Sample[item][1])))
2413
2414	elif 'Debye' in Type: #Debye-Scherrer
2415	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2416	ptlbls += '%14s'%(item)
2417	ptstr += '%14.4f'%(Sample[item][0])
2418	varstr += '%14s'%(str(bool(Sample[item][1])))
2419
2420	print >>pFile,ptlbls
2421	print >>pFile,ptstr
2422	print >>pFile,varstr
2423
2424	histDict = {}
2425	histVary = []
2426	controlDict = {}
2427	histoList = Histograms.keys()
2428	histoList.sort()
2429	for histogram in histoList:
2430	if 'PWDR' in histogram:
2431	Histogram = Histograms[histogram]
2432	hId = Histogram['hId']
2433	pfx = ':'+str(hId)+':'
2434	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2435	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2436
2437	Background = Histogram['Background']
2438	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2439	controlDict[pfx+'bakType'] = Type
2440	histDict.update(bakDict)
2441	histVary += bakVary
2442
2443	Inst = Histogram['Instrument Parameters'][0]
2444	Type,instDict,insVary = GetInstParms(hId,Inst)
2445	controlDict[pfx+'histType'] = Type
2446	if pfx+'Lam1' in instDict:
2447	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2448	else:
2449	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
2450	histDict.update(instDict)
2451	histVary += insVary
2452
2453	Sample = Histogram['Sample Parameters']
2454	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2455	controlDict[pfx+'instType'] = Type
2456	histDict.update(sampDict)
2457	histVary += sampVary
2458
2459
2460	if Print:
2461	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2462	print >>pFile,135*'-'
2463	Units = {'C':' deg','T':' msec'}
2464	units = Units[controlDict[pfx+'histType'][2]]
2465	Limits = controlDict[pfx+'Limits']
2466	print >>pFile,' Instrument type: ',Sample['Type']
2467	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2468	PrintSampleParms(Sample)
2469	PrintInstParms(Inst)
2470	PrintBackground(Background)
2471	elif 'HKLF' in histogram:
2472	Histogram = Histograms[histogram]
2473	hId = Histogram['hId']
2474	pfx = ':'+str(hId)+':'
2475	controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
2476	Inst = Histogram['Instrument Parameters'][0]
2477	controlDict[pfx+'histType'] = Inst['Type'][0]
2478	histDict[pfx+'Lam'] = Inst['Lam'][1]
2479	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
2480	return histVary,histDict,controlDict
2481
2482	def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2483
2484	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2485	Back = Background[0]
2486	DebyePeaks = Background[1]
2487	lenBack = len(Back[3:])
2488	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
2489	for i in range(lenBack):
2490	Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2491	if pfx+'Back:'+str(i) in sigDict:
2492	backSig[i] = sigDict[pfx+'Back:'+str(i)]
2493	if DebyePeaks['nDebye']:
2494	for i in range(DebyePeaks['nDebye']):
2495	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2496	for j,name in enumerate(names):
2497	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2498	if name in sigDict:
2499	backSig[lenBack+3*i+j] = sigDict[name]
2500	if DebyePeaks['nPeaks']:
2501	for i in range(DebyePeaks['nPeaks']):
2502	names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
2503	pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
2504	for j,name in enumerate(names):
2505	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2506	if name in sigDict:
2507	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
2508	return backSig
2509
2510	def SetInstParms(pfx,Inst,parmDict,sigDict):
2511	instSig = {}
2512	insKeys = Inst.keys()
2513	insKeys.sort()
2514	for item in insKeys:
2515	insName = pfx+item
2516	Inst[item][1] = parmDict[insName]
2517	if insName in sigDict:
2518	instSig[item] = sigDict[insName]
2519	else:
2520	instSig[item] = 0
2521	return instSig
2522
2523	def SetSampleParms(pfx,Sample,parmDict,sigDict):
2524	if 'Bragg' in Sample['Type']: #Bragg-Brentano
2525	sampSig = [0 for i in range(3)]
2526	for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering?
2527	Sample[item][0] = parmDict[pfx+item]
2528	if pfx+item in sigDict:
2529	sampSig[i] = sigDict[pfx+item]
2530	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2531	sampSig = [0 for i in range(4)]
2532	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2533	Sample[item][0] = parmDict[pfx+item]
2534	if pfx+item in sigDict:
2535	sampSig[i] = sigDict[pfx+item]
2536	return sampSig
2537
2538	def PrintBackgroundSig(Background,backSig):
2539	Back = Background[0]
2540	DebyePeaks = Background[1]
2541	lenBack = len(Back[3:])
2542	valstr = ' value : '
2543	sigstr = ' sig : '
2544	refine = False
2545	for i,back in enumerate(Back[3:]):
2546	valstr += '%10.4g'%(back)
2547	if Back[1]:
2548	refine = True
2549	sigstr += '%10.4g'%(backSig[i])
2550	else:
2551	sigstr += 10*' '
2552	if refine:
2553	print >>pFile,'\n Background function: ',Back[0]
2554	print >>pFile,valstr
2555	print >>pFile,sigstr
2556	if DebyePeaks['nDebye']:
2557	ifAny = False
2558	ptfmt = "%12.3f"
2559	names = ' names :'
2560	ptstr = ' values:'
2561	sigstr = ' esds :'
2562	for item in sigDict:
2563	if 'Debye' in item:
2564	ifAny = True
2565	names += '%12s'%(item)
2566	ptstr += ptfmt%(parmDict[item])
2567	sigstr += ptfmt%(sigDict[item])
2568	if ifAny:
2569	print >>pFile,'\n Debye diffuse scattering coefficients'
2570	print >>pFile,names
2571	print >>pFile,ptstr
2572	print >>pFile,sigstr
2573	if DebyePeaks['nPeaks']:
2574	ifAny = False
2575	ptfmt = "%14.3f"
2576	names = ' names :'
2577	ptstr = ' values:'
2578	sigstr = ' esds :'
2579	for item in sigDict:
2580	if 'BkPk' in item:
2581	ifAny = True
2582	names += '%14s'%(item)
2583	ptstr += ptfmt%(parmDict[item])
2584	sigstr += ptfmt%(sigDict[item])
2585	if ifAny:
2586	print >>pFile,'\n Single peak coefficients'
2587	print >>pFile,names
2588	print >>pFile,ptstr
2589	print >>pFile,sigstr
2590
2591	def PrintInstParmsSig(Inst,instSig):
2592	ptlbls = ' names :'
2593	ptstr = ' value :'
2594	sigstr = ' sig :'
2595	refine = False
2596	insKeys = instSig.keys()
2597	insKeys.sort()
2598	for name in insKeys:
2599	if name not in ['Type','Lam1','Lam2','Azimuth']:
2600	ptlbls += '%12s' % (name)
2601	ptstr += '%12.6f' % (Inst[name][1])
2602	if instSig[name]:
2603	refine = True
2604	sigstr += '%12.6f' % (instSig[name])
2605	else:
2606	sigstr += 12*' '
2607	if refine:
2608	print >>pFile,'\n Instrument Parameters:'
2609	print >>pFile,ptlbls
2610	print >>pFile,ptstr
2611	print >>pFile,sigstr
2612
2613	def PrintSampleParmsSig(Sample,sampleSig):
2614	ptlbls = ' names :'
2615	ptstr = ' values:'
2616	sigstr = ' sig :'
2617	refine = False
2618	if 'Bragg' in Sample['Type']:
2619	for i,item in enumerate(['Scale','Shift','Transparency']):
2620	ptlbls += '%14s'%(item)
2621	ptstr += '%14.4f'%(Sample[item][0])
2622	if sampleSig[i]:
2623	refine = True
2624	sigstr += '%14.4f'%(sampleSig[i])
2625	else:
2626	sigstr += 14*' '
2627
2628	elif 'Debye' in Sample['Type']: #Debye-Scherrer
2629	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2630	ptlbls += '%14s'%(item)
2631	ptstr += '%14.4f'%(Sample[item][0])
2632	if sampleSig[i]:
2633	refine = True
2634	sigstr += '%14.4f'%(sampleSig[i])
2635	else:
2636	sigstr += 14*' '
2637
2638	if refine:
2639	print >>pFile,'\n Sample Parameters:'
2640	print >>pFile,ptlbls
2641	print >>pFile,ptstr
2642	print >>pFile,sigstr
2643
2644	histoList = Histograms.keys()
2645	histoList.sort()
2646	for histogram in histoList:
2647	if 'PWDR' in histogram:
2648	Histogram = Histograms[histogram]
2649	hId = Histogram['hId']
2650	pfx = ':'+str(hId)+':'
2651	Background = Histogram['Background']
2652	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2653
2654	Inst = Histogram['Instrument Parameters'][0]
2655	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2656
2657	Sample = Histogram['Sample Parameters']
2658	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2659
2660	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2661	print >>pFile,135*'-'
2662	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
2663	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2664	if Print:
2665	print >>pFile,' Instrument type: ',Sample['Type']
2666	PrintSampleParmsSig(Sample,sampSig)
2667	PrintInstParmsSig(Inst,instSig)
2668	PrintBackgroundSig(Background,backSig)
2669
2670	################################################################################
2671	##### Penalty & restraint functions
2672	################################################################################
2673
2674	def penaltyFxn(HistoPhases,parmDict,varyList):
2675	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2676	pNames = []
2677	pVals = []
2678	pWt = []
2679	negWt = {}
2680	for phase in Phases:
2681	pId = Phases[phase]['pId']
2682	negWt[pId] = Phases[phase]['General']['Pawley neg wt']
2683	General = Phases[phase]['General']
2684	textureData = General['SH Texture']
2685	SGData = General['SGData']
2686	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
2687	cell = General['Cell'][1:7]
2688	Amat,Bmat = G2lat.cell2AB(cell)
2689	if phase not in restraintDict:
2690	continue
2691	phaseRest = restraintDict[phase]
2692	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
2693	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
2694	['ChemComp','Sites'],['Texture','HKLs']]
2695	for name,rest in names:
2696	itemRest = phaseRest[name]
2697	if itemRest[rest] and itemRest['Use']:
2698	wt = itemRest['wtFactor']
2699	if name in ['Bond','Angle','Plane','Chiral']:
2700	for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
2701	pNames.append(str(pId)+':'+name+':'+str(i))
2702	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2703	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2704	if name == 'Bond':
2705	calc = G2mth.getRestDist(XYZ,Amat)
2706	elif name == 'Angle':
2707	calc = G2mth.getRestAngle(XYZ,Amat)
2708	elif name == 'Plane':
2709	calc = G2mth.getRestPlane(XYZ,Amat)
2710	elif name == 'Chiral':
2711	calc = G2mth.getRestChiral(XYZ,Amat)
2712	pVals.append(obs-calc)
2713	pWt.append(wt/esd**2)
2714	elif name in ['Torsion','Rama']:
2715	coeffDict = itemRest['Coeff']
2716	for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
2717	pNames.append(str(pId)+':'+name+':'+str(i))
2718	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2719	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2720	if name == 'Torsion':
2721	tor = G2mth.getRestTorsion(XYZ,Amat)
2722	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
2723	else:
2724	phi,psi = G2mth.getRestRama(XYZ,Amat)
2725	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
2726	pVals.append(obs-calc)
2727	pWt.append(wt/esd**2)
2728	elif name == 'ChemComp':
2729	for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
2730	pNames.append(str(pId)+':'+name+':'+str(i))
2731	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
2732	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
2733	calc = np.sum(mulfracfactors)
2734	pVals.append(obs-calc)
2735	pWt.append(wt/esd**2)
2736	elif name == 'Texture':
2737	SHkeys = textureData['SH Coeff'][1].keys()
2738	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
2739	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2740	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2741	for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
2742	PH = np.array(hkl)
2743	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2744	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
2745	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
2746	Z1 = -ma.masked_greater(Z,0.0)
2747	IndZ1 = np.array(ma.nonzero(Z1))
2748	for ind in IndZ1.T:
2749	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
2750	pVals.append(Z1[ind[0]][ind[1]])
2751	pWt.append(wt/esd1**2)
2752	if ifesd2:
2753	Z2 = 1.-Z
2754	for ind in np.ndindex(grid,grid):
2755	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
2756	pVals.append(Z1[ind[0]][ind[1]])
2757	pWt.append(wt/esd2**2)
2758
2759	for item in varyList:
2760	if 'PWLref' in item and parmDict[item] < 0.:
2761	pId = int(item.split(':')[0])
2762	if negWt[pId]:
2763	pNames.append(item)
2764	pVals.append(-parmDict[item])
2765	pWt.append(negWt[pId])
2766	pVals = np.array(pVals)
2767	pWt = np.array(pWt) #should this be np.sqrt?
2768	return pNames,pVals,pWt
2769
2770	def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
2771	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2772	pDerv = np.zeros((len(varyList),len(pVal)))
2773	for phase in Phases:
2774	# if phase not in restraintDict:
2775	# continue
2776	pId = Phases[phase]['pId']
2777	General = Phases[phase]['General']
2778	SGData = General['SGData']
2779	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
2780	cell = General['Cell'][1:7]
2781	Amat,Bmat = G2lat.cell2AB(cell)
2782	textureData = General['SH Texture']
2783
2784	SHkeys = textureData['SH Coeff'][1].keys()
2785	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
2786	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2787	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2788	sam = SamSym[textureData['Model']]
2789	phaseRest = restraintDict.get(phase,{})
2790	names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
2791	'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
2792	'ChemComp':'Sites','Texture':'HKLs'}
2793	lasthkl = np.array([0,0,0])
2794	for ip,pName in enumerate(pNames):
2795	pnames = pName.split(':')
2796	if pId == int(pnames[0]):
2797	name = pnames[1]
2798	if 'PWL' in pName:
2799	pDerv[varyList.index(pName)][ip] += 1.
2800	continue
2801	id = int(pnames[2])
2802	itemRest = phaseRest[name]
2803	if name in ['Bond','Angle','Plane','Chiral']:
2804	indx,ops,obs,esd = itemRest[names[name]][id]
2805	dNames = []
2806	for ind in indx:
2807	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
2808	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2809	if name == 'Bond':
2810	deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
2811	elif name == 'Angle':
2812	deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
2813	elif name == 'Plane':
2814	deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
2815	elif name == 'Chiral':
2816	deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
2817	elif name in ['Torsion','Rama']:
2818	coffDict = itemRest['Coeff']
2819	indx,ops,cofName,esd = itemRest[names[name]][id]
2820	dNames = []
2821	for ind in indx:
2822	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
2823	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2824	if name == 'Torsion':
2825	deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
2826	else:
2827	deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
2828	elif name == 'ChemComp':
2829	indx,factors,obs,esd = itemRest[names[name]][id]
2830	dNames = []
2831	for ind in indx:
2832	dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
2833	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
2834	deriv = mul*factors
2835	elif 'Texture' in name:
2836	deriv = []
2837	dNames = []
2838	hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
2839	hkl = np.array(hkl)
2840	if np.any(lasthkl-hkl):
2841	PH = np.array(hkl)
2842	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2843	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
2844	lasthkl = copy.copy(hkl)
2845	if 'unit' in name:
2846	pass
2847	else:
2848	gam = float(pnames[3])
2849	psi = float(pnames[4])
2850	for SHname in ODFln:
2851	l,m,n = eval(SHname[1:])
2852	Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
2853	dNames += [str(pId)+'::'+SHname]
2854	deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
2855	for dName,drv in zip(dNames,deriv):
2856	try:
2857	ind = varyList.index(dName)
2858	pDerv[ind][ip] += drv
2859	except ValueError:
2860	pass
2861	return pDerv
2862
2863	################################################################################
2864	##### Function & derivative calculations
2865	################################################################################
2866
2867	def GetAtomFXU(pfx,calcControls,parmDict):
2868	Natoms = calcControls['Natoms'][pfx]
2869	Tdata = Natoms*[' ',]
2870	Mdata = np.zeros(Natoms)
2871	IAdata = Natoms*[' ',]
2872	Fdata = np.zeros(Natoms)
2873	FFdata = []
2874	BLdata = []
2875	Xdata = np.zeros((3,Natoms))
2876	dXdata = np.zeros((3,Natoms))
2877	Uisodata = np.zeros(Natoms)
2878	Uijdata = np.zeros((6,Natoms))
2879	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
2880	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
2881	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
2882	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
2883	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
2884	for iatm in range(Natoms):
2885	for key in keys:
2886	parm = pfx+key+str(iatm)
2887	if parm in parmDict:
2888	keys[key][iatm] = parmDict[parm]
2889	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
2890
2891	def getFFvalues(FFtables,SQ):
2892	FFvals = {}
2893	for El in FFtables:
2894	FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
2895	return FFvals
2896
2897	def getBLvalues(BLtables):
2898	BLvals = {}
2899	for El in BLtables:
2900	BLvals[El] = BLtables[El][1][1]
2901	return BLvals
2902
2903	def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2904	''' Compute structure factors for all h,k,l for phase
2905	input:
2906	refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2907	G: reciprocal metric tensor
2908	pfx: phase id string
2909	SGData: space group info. dictionary output from SpcGroup
2910	calcControls:
2911	ParmDict:
2912	puts result F^2 in each ref[8] in refList
2913	'''
2914	twopi = 2.0*np.pi
2915	twopisq = 2.0np.pi*2
2916	phfx = pfx.split(':')[0]+hfx
2917	ast = np.sqrt(np.diag(G))
2918	Mast = twopisq*np.multiply.outer(ast,ast)
2919	FFtables = calcControls['FFtables']
2920	BLtables = calcControls['BLtables']
2921	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2922	FF = np.zeros(len(Tdata))
2923	if 'N' in parmDict[hfx+'Type']:
2924	FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
2925	else:
2926	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2927	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2928	maxPos = len(SGData['SGOps'])
2929	Uij = np.array(G2lat.U6toUij(Uijdata))
2930	bij = Mast*Uij.T
2931	for refl in refList:
2932	fbs = np.array([0,0])
2933	H = refl[:3]
2934	SQ = 1./(2.refl[4])*2
2935	SQfactor = 4.0SQtwopisq
2936	Bab = parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor)
2937	if not len(refl[-1]): #no form factors
2938	if 'N' in parmDict[hfx+'Type']:
2939	refl[-1] = getBLvalues(BLtables)
2940	else: #'X'
2941	refl[-1] = getFFvalues(FFtables,SQ)
2942	for i,El in enumerate(Tdata):
2943	FF[i] = refl[-1][El]
2944	Uniq = refl[11]
2945	phi = refl[12]
2946	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
2947	sinp = np.sin(phase)
2948	cosp = np.cos(phase)
2949	occ = Mdata*Fdata/len(Uniq)
2950	biso = -SQfactor*Uisodata
2951	Tiso = np.where(biso<1.,np.exp(biso),1.0)
2952	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
2953	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
2954	Tcorr = Tiso*Tuij
2955	fa = np.array([(FF+FP-Bab)occcospTcorr,-FPPoccsinpTcorr])
2956	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
2957	if not SGData['SGInv']:
2958	fb = np.array([(FF+FP-Bab)occsinpTcorr,FPPocccospTcorr])
2959	fbs = np.sum(np.sum(fb,axis=1),axis=1)
2960	fasq = fas**2
2961	fbsq = fbs**2 #imaginary
2962	refl[9] = np.sum(fasq)+np.sum(fbsq)
2963	refl[10] = atan2d(fbs[0],fas[0])
2964
2965	def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2966	twopi = 2.0*np.pi
2967	twopisq = 2.0np.pi*2
2968	phfx = pfx.split(':')[0]+hfx
2969	ast = np.sqrt(np.diag(G))
2970	Mast = twopisq*np.multiply.outer(ast,ast)
2971	FFtables = calcControls['FFtables']
2972	BLtables = calcControls['BLtables']
2973	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2974	FF = np.zeros(len(Tdata))
2975	if 'N' in parmDict[hfx+'Type']:
2976	FP = 0.
2977	FPP = 0.
2978	else:
2979	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2980	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2981	maxPos = len(SGData['SGOps'])
2982	Uij = np.array(G2lat.U6toUij(Uijdata))
2983	bij = Mast*Uij.T
2984	dFdvDict = {}
2985	dFdfr = np.zeros((len(refList),len(Mdata)))
2986	dFdx = np.zeros((len(refList),len(Mdata),3))
2987	dFdui = np.zeros((len(refList),len(Mdata)))
2988	dFdua = np.zeros((len(refList),len(Mdata),6))
2989	dFdbab = np.zeros((len(refList),2))
2990	for iref,refl in enumerate(refList):
2991	H = np.array(refl[:3])
2992	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
2993	SQfactor = 8.0SQnp.pi**2
2994	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2995	Bab = parmDict[phfx+'BabA']*dBabdA
2996	for i,El in enumerate(Tdata):
2997	FF[i] = refl[-1][El]
2998	Uniq = refl[11]
2999	phi = refl[12]
3000	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
3001	sinp = np.sin(phase)
3002	cosp = np.cos(phase)
3003	occ = Mdata*Fdata/len(Uniq)
3004	biso = -SQfactor*Uisodata
3005	Tiso = np.where(biso<1.,np.exp(biso),1.0)
3006	HbH = -np.inner(H,np.inner(bij,H))
3007	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
3008	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
3009	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
3010	Tcorr = Tiso*Tuij
3011	fot = (FF+FP-Bab)occTcorr
3012	fotp = FPPoccTcorr
3013	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
3014	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
3015
3016	fas = np.sum(np.sum(fa,axis=1),axis=1)
3017	fbs = np.sum(np.sum(fb,axis=1),axis=1)
3018	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
3019	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
3020	#sum below is over Uniq
3021	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
3022	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
3023	dfadui = np.sum(-SQfactor*fa,axis=2)
3024	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
3025	dfadba = np.sum(-cosp(occTcorr)[:,np.newaxis],axis=1)
3026	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for
3027	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
3028	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
3029	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
3030	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
3031	dFdbab[iref] = np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
3032	if not SGData['SGInv']:
3033	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom
3034	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
3035	dfbdui = np.sum(-SQfactor*fb,axis=2)
3036	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
3037	dfbdba = np.sum(-sinp(occTcorr)[:,np.newaxis],axis=1)
3038	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
3039	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
3040	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
3041	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
3042	dFdbab[iref] += np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
3043	#loop over atoms - each dict entry is list of derivatives for all the reflections
3044	for i in range(len(Mdata)):
3045	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
3046	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
3047	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
3048	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
3049	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
3050	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
3051	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
3052	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
3053	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
3054	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
3055	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
3056	dFdvDict[pfx+'BabA'] = dFdbab.T[0]
3057	dFdvDict[pfx+'BabU'] = dFdbab.T[1]
3058	return dFdvDict
3059
3060	def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
3061	''' Single crystal extinction function; puts correction in ref[13] and returns
3062	corrections needed for derivatives
3063	'''
3064	ref[13] = 1.0
3065	dervCor = 1.0
3066	dervDict = {}
3067	if calcControls[phfx+'EType'] != 'None':
3068	cos2T = 1.0-0.5(parmDict[hfx+'Lam']/ref[4])*2 #cos(2theta)
3069	if 'SXC' in parmDict[hfx+'Type']:
3070	AV = 7.9406e5/parmDict[pfx+'Vol']**2
3071	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
3072	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
3073	elif 'SNT' in parmDict[hfx+'Type']:
3074	AV = 1.e7/parmDict[pfx+'Vol']**2
3075	PL = 1./(4.refl[4]*2)
3076	P12 = 1.0
3077	elif 'SNC' in parmDict[hfx+'Type']:
3078	AV = 1.e7/parmDict[pfx+'Vol']**2
3079	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
3080	P12 = 1.0
3081
3082	PLZ = AVP12parmDict[hfx+'Lam']*2ref[7]
3083	if 'Primary' in calcControls[phfx+'EType']:
3084	PLZ *= 1.5
3085	else:
3086	PLZ *= calcControls[phfx+'Tbar']
3087
3088	if 'Primary' in calcControls[phfx+'EType']:
3089	PSIG = parmDict[phfx+'Ep']
3090	elif 'I & II' in calcControls[phfx+'EType']:
3091	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
3092	elif 'Type II' in calcControls[phfx+'EType']:
3093	PSIG = parmDict[phfx+'Es']
3094	else: # 'Secondary Type I'
3095	PSIG = parmDict[phfx+'Eg']/PL
3096
3097	AG = 0.58+0.48cos2T+0.24cos2T**2
3098	AL = 0.025+0.285*cos2T
3099	BG = 0.02-0.025*cos2T
3100	BL = 0.15-0.2(0.75-cos2T)*2
3101	if cos2T < 0.:
3102	BL = -0.45*cos2T
3103	CG = 2.
3104	CL = 2.
3105	PF = PLZ*PSIG
3106
3107	if 'Gaussian' in calcControls[phfx+'EApprox']:
3108	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
3109	extCor = np.sqrt(PF4)
3110	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
3111	else:
3112	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
3113	extCor = np.sqrt(PF4)
3114	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
3115
3116	dervCor = (1.+PF)PF3
3117	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
3118	dervDict[phfx+'Ep'] = -ref[7]PLZPF3
3119	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
3120	dervDict[phfx+'Es'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
3121	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
3122	dervDict[phfx+'Eg'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
3123
3124	ref[13] = 1./extCor
3125	return dervCor,dervDict
3126
3127
3128	def Dict2Values(parmdict, varylist):
3129	'''Use before call to leastsq to setup list of values for the parameters
3130	in parmdict, as selected by key in varylist'''
3131	return [parmdict[key] for key in varylist]
3132
3133	def Values2Dict(parmdict, varylist, values):
3134	''' Use after call to leastsq to update the parameter dictionary with
3135	values corresponding to keys in varylist'''
3136	parmdict.update(zip(varylist,values))
3137
3138	def GetNewCellParms(parmDict,varyList):
3139	newCellDict = {}
3140	Anames = ['A'+str(i) for i in range(6)]
3141	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
3142	for item in varyList:
3143	keys = item.split(':')
3144	if keys[2] in Ddict:
3145	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
3146	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
3147	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
3148	return newCellDict # is e.g. {'0::D11':A0+D11}
3149
3150	def ApplyXYZshifts(parmDict,varyList):
3151	''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
3152	input:
3153	parmDict - parameter dictionary
3154	varyList - list of variables
3155	returns:
3156	newAtomDict - dictionary of new atomic coordinate names & values;
3157	key is parameter shift name
3158	'''
3159	newAtomDict = {}
3160	for item in parmDict:
3161	if 'dA' in item:
3162	parm = ''.join(item.split('d'))
3163	parmDict[parm] += parmDict[item]
3164	newAtomDict[item] = [parm,parmDict[parm]]
3165	return newAtomDict
3166
3167	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
3168	#Spherical harmonics texture
3169	IFCoup = 'Bragg' in calcControls[hfx+'instType']
3170	odfCor = 1.0
3171	H = refl[:3]
3172	cell = G2lat.Gmat2cell(g)
3173	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
3174	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
3175	phi,beta = G2lat.CrsAng(H,cell,SGData)
3176	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
3177	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
3178	for item in SHnames:
3179	L,M,N = eval(item.strip('C'))
3180	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
3181	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
3182	Lnorm = G2lat.Lnorm(L)
3183	odfCor += parmDict[pfx+item]LnormKcl*Ksl
3184	return odfCor
3185
3186	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
3187	#Spherical harmonics texture derivatives
3188	FORPI = 4.0*np.pi
3189	IFCoup = 'Bragg' in calcControls[hfx+'instType']
3190	odfCor = 1.0
3191	dFdODF = {}
3192	dFdSA = [0,0,0]
3193	H = refl[:3]
3194	cell = G2lat.Gmat2cell(g)
3195	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
3196	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
3197	phi,beta = G2lat.CrsAng(H,cell,SGData)
3198	psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
3199	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
3200	for item in SHnames:
3201	L,M,N = eval(item.strip('C'))
3202	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
3203	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
3204	Lnorm = G2lat.Lnorm(L)
3205	odfCor += parmDict[pfx+item]LnormKcl*Ksl
3206	dFdODF[pfx+item] = LnormKclKsl
3207	for i in range(3):
3208	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
3209	return odfCor,dFdODF,dFdSA
3210
3211	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
3212	#sphericaql harmonics preferred orientation (cylindrical symmetry only)
3213	odfCor = 1.0
3214	H = refl[:3]
3215	cell = G2lat.Gmat2cell(g)
3216	Sangl = [0.,0.,0.]
3217	if 'Bragg' in calcControls[hfx+'instType']:
3218	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
3219	IFCoup = True
3220	else:
3221	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
3222	IFCoup = False
3223	phi,beta = G2lat.CrsAng(H,cell,SGData)
3224	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
3225	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
3226	for item in SHnames:
3227	L,N = eval(item.strip('C'))
3228	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
3229	odfCor += parmDict[phfx+item]LnormKcsl
3230	return odfCor
3231
3232	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
3233	#spherical harmonics preferred orientation derivatives (cylindrical symmetry only)
3234	FORPI = 12.5663706143592
3235	odfCor = 1.0
3236	dFdODF = {}
3237	H = refl[:3]
3238	cell = G2lat.Gmat2cell(g)
3239	Sangl = [0.,0.,0.]
3240	if 'Bragg' in calcControls[hfx+'instType']:
3241	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
3242	IFCoup = True
3243	else:
3244	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
3245	IFCoup = False
3246	phi,beta = G2lat.CrsAng(H,cell,SGData)
3247	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
3248	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
3249	for item in SHnames:
3250	L,N = eval(item.strip('C'))
3251	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
3252	odfCor += parmDict[phfx+item]LnormKcsl
3253	dFdODF[phfx+item] = Kcsl*Lnorm
3254	return odfCor,dFdODF
3255
3256	def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
3257	POcorr = 1.0
3258	if calcControls[phfx+'poType'] == 'MD':
3259	MD = parmDict[phfx+'MD']
3260	if MD != 1.0:
3261	MDAxis = calcControls[phfx+'MDAxis']
3262	sumMD = 0
3263	for H in refl[11]:
3264	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
3265	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
3266	sumMD += A**3
3267	POcorr = sumMD/len(refl[11])
3268	else: #spherical harmonics
3269	if calcControls[phfx+'SHord']:
3270	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
3271	return POcorr
3272
3273	def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
3274	POcorr = 1.0
3275	POderv = {}
3276	if calcControls[phfx+'poType'] == 'MD':
3277	MD = parmDict[phfx+'MD']
3278	MDAxis = calcControls[phfx+'MDAxis']
3279	sumMD = 0
3280	sumdMD = 0
3281	for H in refl[11]:
3282	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
3283	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
3284	sumMD += A**3
3285	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
3286	POcorr = sumMD/len(refl[11])
3287	POderv[phfx+'MD'] = sumdMD/len(refl[11])
3288	else: #spherical harmonics
3289	if calcControls[phfx+'SHord']:
3290	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
3291	return POcorr,POderv
3292
3293	def GetAbsorb(refl,hfx,calcControls,parmDict):
3294	if 'Debye' in calcControls[hfx+'instType']:
3295	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
3296	else:
3297	return 1.0
3298
3299	def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
3300	if 'Debye' in calcControls[hfx+'instType']:
3301	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
3302	else:
3303	return 0.0
3304
3305	def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
3306	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
3307	if 'X' in parmDict[hfx+'Type']:
3308	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
3309	Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
3310	if pfx+'SHorder' in parmDict:
3311	Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
3312	Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
3313	refl[13] = Icorr
3314
3315	def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
3316	dIdsh = 1./parmDict[hfx+'Scale']
3317	dIdsp = 1./parmDict[phfx+'Scale']
3318	if 'X' in parmDict[hfx+'Type']:
3319	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
3320	dIdPola /= pola
3321	else: #'N'
3322	dIdPola = 0.0
3323	POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
3324	for iPO in dIdPO:
3325	dIdPO[iPO] /= POcorr
3326	dFdODF = {}
3327	dFdSA = [0,0,0]
3328	if pfx+'SHorder' in parmDict:
3329	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
3330	for iSH in dFdODF:
3331	dFdODF[iSH] /= odfCor
3332	for i in range(3):
3333	dFdSA[i] /= odfCor
3334	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
3335	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
3336
3337	def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
3338	costh = cosd(refl[5]/2.)
3339	#crystallite size
3340	if calcControls[phfx+'SizeType'] == 'isotropic':
3341	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
3342	elif calcControls[phfx+'SizeType'] == 'uniaxial':
3343	H = np.array(refl[:3])
3344	P = np.array(calcControls[phfx+'SizeAxis'])
3345	cosP,sinP = G2lat.CosSinAngle(H,P,G)
3346	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
3347	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
3348	else: #ellipsoidal crystallites
3349	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
3350	H = np.array(refl[:3])
3351	lenR = G2pwd.ellipseSize(H,Sij,GB)
3352	Sgam = 1.8wave/(np.picosth*lenR)
3353	#microstrain
3354	if calcControls[phfx+'MustrainType'] == 'isotropic':
3355	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
3356	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
3357	H = np.array(refl[:3])
3358	P = np.array(calcControls[phfx+'MustrainAxis'])
3359	cosP,sinP = G2lat.CosSinAngle(H,P,G)
3360	Si = parmDict[phfx+'Mustrain;i']
3361	Sa = parmDict[phfx+'Mustrain;a']
3362	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
3363	else: #generalized - P.W. Stephens model
3364	pwrs = calcControls[phfx+'MuPwrs']
3365	sum = 0
3366	for i,pwr in enumerate(pwrs):
3367	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
3368	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
3369	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
3370	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
3371	sig /= ateln2
3372	return sig,gam
3373
3374	def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
3375	gamDict = {}
3376	sigDict = {}
3377	costh = cosd(refl[5]/2.)
3378	tanth = tand(refl[5]/2.)
3379	#crystallite size derivatives
3380	if calcControls[phfx+'SizeType'] == 'isotropic':
3381	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
3382	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
3383	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
3384	elif calcControls[phfx+'SizeType'] == 'uniaxial':
3385	H = np.array(refl[:3])
3386	P = np.array(calcControls[phfx+'SizeAxis'])
3387	cosP,sinP = G2lat.CosSinAngle(H,P,G)
3388	Si = parmDict[phfx+'Size;i']
3389	Sa = parmDict[phfx+'Size;a']
3390	gami = (1.8wave/np.pi)/(SiSa)
3391	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
3392	Sgam = gami*sqtrm
3393	gam = Sgam/costh
3394	dsi = (gamiSicosP*2/(sqtrmcosth)-gam/Si)
3395	dsa = (gamiSasinP*2/(sqtrmcosth)-gam/Sa)
3396	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
3397	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
3398	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
3399	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
3400	else: #ellipsoidal crystallites
3401	const = 1.8wave/(np.picosth)
3402	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
3403	H = np.array(refl[:3])
3404	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
3405	Sgam = 1.8wave/(np.picosth*lenR)
3406	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
3407	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
3408	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
3409	gamDict[phfx+'Size;mx'] = Sgam
3410	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
3411
3412	#microstrain derivatives
3413	if calcControls[phfx+'MustrainType'] == 'isotropic':
3414	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
3415	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
3416	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
3417	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
3418	H = np.array(refl[:3])
3419	P = np.array(calcControls[phfx+'MustrainAxis'])
3420	cosP,sinP = G2lat.CosSinAngle(H,P,G)
3421	Si = parmDict[phfx+'Mustrain;i']
3422	Sa = parmDict[phfx+'Mustrain;a']
3423	gami = 0.018SiSa*tanth/np.pi
3424	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
3425	Mgam = gami/sqtrm
3426	dsi = -gamiSicosP2/sqtrm3
3427	dsa = -gamiSasinP2/sqtrm3
3428	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
3429	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
3430	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
3431	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
3432	else: #generalized - P.W. Stephens model
3433	pwrs = calcControls[phfx+'MuPwrs']
3434	const = 0.018refl[4]2tanth
3435	sum = 0
3436	for i,pwr in enumerate(pwrs):
3437	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
3438	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
3439	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain;mx']
3440	sigDict[phfx+'Mustrain:'+str(i)] = \
3441	2.constterm(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
3442	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
3443	for i in range(len(pwrs)):
3444	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
3445	gamDict[phfx+'Mustrain;mx'] = Mgam
3446	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
3447	return sigDict,gamDict
3448
3449	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
3450	h,k,l = refl[:3]
3451	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
3452	d = np.sqrt(dsq)
3453
3454	refl[4] = d
3455	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
3456	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
3457	if 'Bragg' in calcControls[hfx+'instType']:
3458	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
3459	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
3460	else: #Debye-Scherrer - simple but maybe not right
3461	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
3462	return pos
3463
3464	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
3465	dpr = 180./np.pi
3466	h,k,l = refl[:3]
3467	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
3468	dst = np.sqrt(dstsq)
3469	pos = refl[5]-parmDict[hfx+'Zero']
3470	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
3471	dpdw = const*dst
3472	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
3473	dpdA = constwave/(2.0*dst)
3474	dpdZ = 1.0
3475	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
3476	if 'Bragg' in calcControls[hfx+'instType']:
3477	dpdSh = -4.constcosd(pos/2.0)
3478	dpdTr = -constsind(pos)100.0
3479	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
3480	else: #Debye-Scherrer - simple but maybe not right
3481	dpdXd = -const*cosd(pos)
3482	dpdYd = -const*sind(pos)
3483	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
3484
3485	def GetHStrainShift(refl,SGData,phfx,parmDict):
3486	laue = SGData['SGLaue']
3487	uniq = SGData['SGUniq']
3488	h,k,l = refl[:3]
3489	if laue in ['m3','m3m']:
3490	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
3491	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
3492	elif laue in ['6/m','6/mmm','3m1','31m','3']:
3493	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
3494	elif laue in ['3R','3mR']:
3495	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
3496	elif laue in ['4/m','4/mmm']:
3497	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
3498	elif laue in ['mmm']:
3499	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
3500	elif laue in ['2/m']:
3501	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
3502	if uniq == 'a':
3503	Dij += parmDict[phfx+'D23']kl
3504	elif uniq == 'b':
3505	Dij += parmDict[phfx+'D13']hl
3506	elif uniq == 'c':
3507	Dij += parmDict[phfx+'D12']hk
3508	else:
3509	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
3510	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
3511	return -Dijrefl[4]2tand(refl[5]/2.0)
3512
3513	def GetHStrainShiftDerv(refl,SGData,phfx):
3514	laue = SGData['SGLaue']
3515	uniq = SGData['SGUniq']
3516	h,k,l = refl[:3]
3517	if laue in ['m3','m3m']:
3518	dDijDict = {phfx+'D11':h2+k2+l**2,
3519	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
3520	elif laue in ['6/m','6/mmm','3m1','31m','3']:
3521	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
3522	elif laue in ['3R','3mR']:
3523	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
3524	elif laue in ['4/m','4/mmm']:
3525	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
3526	elif laue in ['mmm']:
3527	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
3528	elif laue in ['2/m']:
3529	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
3530	if uniq == 'a':
3531	dDijDict[phfx+'D23'] = k*l
3532	elif uniq == 'b':
3533	dDijDict[phfx+'D13'] = h*l
3534	elif uniq == 'c':
3535	dDijDict[phfx+'D12'] = h*k
3536	names.append()
3537	else:
3538	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
3539	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
3540	for item in dDijDict:
3541	dDijDict[item] = -refl[4]2tand(refl[5]/2.0)
3542	return dDijDict
3543
3544	def GetFprime(controlDict,Histograms):
3545	FFtables = controlDict['FFtables']
3546	if not FFtables:
3547	return
3548	histoList = Histograms.keys()
3549	histoList.sort()
3550	for histogram in histoList:
3551	if histogram[:4] in ['PWDR','HKLF']:
3552	Histogram = Histograms[histogram]
3553	hId = Histogram['hId']
3554	hfx = ':%d:'%(hId)
3555	keV = controlDict[hfx+'keV']
3556	for El in FFtables:
3557	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
3558	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
3559	FFtables[El][hfx+'FP'] = FP
3560	FFtables[El][hfx+'FPP'] = FPP
3561
3562	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
3563	histoList = Histograms.keys()
3564	histoList.sort()
3565	for histogram in histoList:
3566	if 'PWDR' in histogram[:4]:
3567	Histogram = Histograms[histogram]
3568	hId = Histogram['hId']
3569	hfx = ':%d:'%(hId)
3570	Limits = calcControls[hfx+'Limits']
3571	shl = max(parmDict[hfx+'SH/L'],0.0005)
3572	Ka2 = False
3573	kRatio = 0.0
3574	if hfx+'Lam1' in parmDict.keys():
3575	Ka2 = True
3576	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
3577	kRatio = parmDict[hfx+'I(L2)/I(L1)']
3578	x,y,w,yc,yb,yd = Histogram['Data']
3579	xB = np.searchsorted(x,Limits[0])
3580	xF = np.searchsorted(x,Limits[1])
3581	ymb = np.array(y-yb)
3582	ymb = np.where(ymb,ymb,1.0)
3583	ycmb = np.array(yc-yb)
3584	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
3585	refLists = Histogram['Reflection Lists']
3586	for phase in refLists:
3587	Phase = Phases[phase]
3588	pId = Phase['pId']
3589	phfx = '%d:%d:'%(pId,hId)
3590	refList = refLists[phase]
3591	sumFo = 0.0
3592	sumdF = 0.0
3593	sumFosq = 0.0
3594	sumdFsq = 0.0
3595	for refl in refList:
3596	if 'C' in calcControls[hfx+'histType']:
3597	yp = np.zeros_like(yb)
3598	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3599	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
3600	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
3601	iFin2 = iFin
3602	yp[iBeg:iFin] = refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
3603	if Ka2:
3604	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
3605	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
3606	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
3607	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
3608	if not iBeg2+iFin2: #peak below low limit - skip peak
3609	continue
3610	elif not iBeg2-iFin2: #peak above high limit - done
3611	break
3612	yp[iBeg2:iFin2] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here
3613	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[13](1.+kRatio)),0.0))
3614	elif 'T' in calcControls[hfx+'histType']:
3615	print 'TOF Undefined at present'
3616	raise Exception #no TOF yet
3617	Fo = np.sqrt(np.abs(refl[8]))
3618	Fc = np.sqrt(np.abs(refl[9]))
3619	sumFo += Fo
3620	sumFosq += refl[8]**2
3621	sumdF += np.abs(Fo-Fc)
3622	sumdFsq += (refl[8]-refl[9])**2
3623	Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3624	Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3625	Histogram[phfx+'Nref'] = len(refList)
3626
3627	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3628
3629	def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
3630	U = parmDict[hfx+'U']
3631	V = parmDict[hfx+'V']
3632	W = parmDict[hfx+'W']
3633	X = parmDict[hfx+'X']
3634	Y = parmDict[hfx+'Y']
3635	tanPos = tand(refl[5]/2.0)
3636	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
3637	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
3638	sig = max(0.001,sig)
3639	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
3640	gam = max(0.001,gam)
3641	return sig,gam
3642
3643	hId = Histogram['hId']
3644	hfx = ':%d:'%(hId)
3645	bakType = calcControls[hfx+'bakType']
3646	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
3647	yc = np.zeros_like(yb)
3648
3649	if 'C' in calcControls[hfx+'histType']:
3650	shl = max(parmDict[hfx+'SH/L'],0.002)
3651	Ka2 = False
3652	if hfx+'Lam1' in parmDict.keys():
3653	wave = parmDict[hfx+'Lam1']
3654	Ka2 = True
3655	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
3656	kRatio = parmDict[hfx+'I(L2)/I(L1)']
3657	else:
3658	wave = parmDict[hfx+'Lam']
3659	else:
3660	print 'TOF Undefined at present'
3661	raise ValueError
3662	for phase in Histogram['Reflection Lists']:
3663	refList = Histogram['Reflection Lists'][phase]
3664	Phase = Phases[phase]
3665	pId = Phase['pId']
3666	pfx = '%d::'%(pId)
3667	phfx = '%d:%d:'%(pId,hId)
3668	hfx = ':%d:'%(hId)
3669	SGData = Phase['General']['SGData']
3670	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
3671	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3672	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
3673	Vst = np.sqrt(nl.det(G)) #V*
3674	if not Phase['General'].get('doPawley'):
3675	time0 = time.time()
3676	StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
3677	# print 'sf calc time: %.3fs'%(time.time()-time0)
3678	time0 = time.time()
3679	for refl in refList:
3680	if 'C' in calcControls[hfx+'histType']:
3681	h,k,l = refl[:3]
3682	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
3683	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
3684	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
3685	refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam
3686	GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13]
3687	refl[13] = VstLorenz
3688	if Phase['General'].get('doPawley'):
3689	try:
3690	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3691	refl[9] = parmDict[pInd]
3692	except KeyError:
3693	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
3694	continue
3695	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3696	iBeg = np.searchsorted(x,refl[5]-fmin)
3697	iFin = np.searchsorted(x,refl[5]+fmax)
3698	if not iBeg+iFin: #peak below low limit - skip peak
3699	continue
3700	elif not iBeg-iFin: #peak above high limit - done
3701	break
3702	yc[iBeg:iFin] += refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
3703	if Ka2:
3704	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
3705	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
3706	iBeg = np.searchsorted(x,pos2-fmin)
3707	iFin = np.searchsorted(x,pos2+fmax)
3708	if not iBeg+iFin: #peak below low limit - skip peak
3709	continue
3710	elif not iBeg-iFin: #peak above high limit - done
3711	return yc,yb
3712	yc[iBeg:iFin] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here
3713	elif 'T' in calcControls[hfx+'histType']:
3714	print 'TOF Undefined at present'
3715	raise Exception #no TOF yet
3716	# print 'profile calc time: %.3fs'%(time.time()-time0)
3717	return yc,yb
3718
3719	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
3720
3721	def cellVaryDerv(pfx,SGData,dpdA):
3722	if SGData['SGLaue'] in ['-1',]:
3723	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
3724	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
3725	elif SGData['SGLaue'] in ['2/m',]:
3726	if SGData['SGUniq'] == 'a':
3727	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
3728	elif SGData['SGUniq'] == 'b':
3729	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
3730	else:
3731	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
3732	elif SGData['SGLaue'] in ['mmm',]:
3733	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
3734	elif SGData['SGLaue'] in ['4/m','4/mmm']:
3735	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
3736	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
3737	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
3738	elif SGData['SGLaue'] in ['3R', '3mR']:
3739	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
3740	elif SGData['SGLaue'] in ['m3m','m3']:
3741	return [[pfx+'A0',dpdA[0]]]
3742
3743	# create a list of dependent variables and set up a dictionary to hold their derivatives
3744	dependentVars = G2mv.GetDependentVars()
3745	depDerivDict = {}
3746	for j in dependentVars:
3747	depDerivDict[j] = np.zeros(shape=(len(x)))
3748	#print 'dependent vars',dependentVars
3749	lenX = len(x)
3750	hId = Histogram['hId']
3751	hfx = ':%d:'%(hId)
3752	bakType = calcControls[hfx+'bakType']
3753	dMdv = np.zeros(shape=(len(varylist),len(x)))
3754	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
3755	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
3756	bBpos =varylist.index(hfx+'Back:0')
3757	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
3758	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
3759	for name in varylist:
3760	if 'Debye' in name:
3761	id = int(name.split(':')[-1])
3762	parm = name[:int(name.rindex(':'))]
3763	ip = names.index(parm)
3764	dMdv[varylist.index(name)] = dMddb[3*id+ip]
3765	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
3766	for name in varylist:
3767	if 'BkPk' in name:
3768	id = int(name.split(':')[-1])
3769	parm = name[:int(name.rindex(':'))]
3770	ip = names.index(parm)
3771	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
3772	cw = np.diff(x)
3773	cw = np.append(cw,cw[-1])
3774	if 'C' in calcControls[hfx+'histType']:
3775	shl = max(parmDict[hfx+'SH/L'],0.002)
3776	Ka2 = False
3777	if hfx+'Lam1' in parmDict.keys():
3778	wave = parmDict[hfx+'Lam1']
3779	Ka2 = True
3780	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
3781	kRatio = parmDict[hfx+'I(L2)/I(L1)']
3782	else:
3783	wave = parmDict[hfx+'Lam']
3784	else:
3785	print 'TOF Undefined at present'
3786	raise ValueError
3787	for phase in Histogram['Reflection Lists']:
3788	refList = Histogram['Reflection Lists'][phase]
3789	Phase = Phases[phase]
3790	SGData = Phase['General']['SGData']
3791	pId = Phase['pId']
3792	pfx = '%d::'%(pId)
3793	phfx = '%d:%d:'%(pId,hId)
3794	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
3795	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3796	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
3797	if not Phase['General'].get('doPawley'):
3798	time0 = time.time()
3799	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
3800	# print 'sf-derv time %.3fs'%(time.time()-time0)
3801	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
3802	time0 = time.time()
3803	for iref,refl in enumerate(refList):
3804	if 'C' in calcControls[hfx+'histType']: #CW powder
3805	h,k,l = refl[:3]
3806	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3807	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3808	iBeg = np.searchsorted(x,refl[5]-fmin)
3809	iFin = np.searchsorted(x,refl[5]+fmax)
3810	if not iBeg+iFin: #peak below low limit - skip peak
3811	continue
3812	elif not iBeg-iFin: #peak above high limit - done
3813	break
3814	pos = refl[5]
3815	tanth = tand(pos/2.0)
3816	costh = cosd(pos/2.0)
3817	lenBF = iFin-iBeg
3818	dMdpk = np.zeros(shape=(6,lenBF))
3819	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
3820	for i in range(5):
3821	dMdpk[i] += 100.cw[iBeg:iFin]refl[13]refl[9]dMdipk[i]
3822	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
3823	if Ka2:
3824	po