Context Navigation

source: trunk/GSASIIstrMath.py @ 1391

Last change on this file since 1391 was 1391, checked in by vondreele, 9 years ago
remove stray prints Breit-Wigner & Lynn & Seeger neutron anomalous form factors for CW data Add 'Ext' to HKLF reflection lists make sure wave & tbar are in TOF HKLF reflection lists add '_refln_f_sigma' & '_refln_f_squared_sigma' to possible cif reflection columns
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 105.6 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMath - structure math routines
4	-----------------------------------------
5	'''
6	########### SVN repository information ###################
7	# $Date: 2014-06-20 18:43:42 +0000 (Fri, 20 Jun 2014) $
8	# $Author: vondreele $
9	# $Revision: 1391 $
10	# $URL: trunk/GSASIIstrMath.py $
11	# $Id: GSASIIstrMath.py 1391 2014-06-20 18:43:42Z vondreele $
12	########### SVN repository information ###################
13	import time
14	import math
15	import copy
16	import numpy as np
17	import numpy.ma as ma
18	import numpy.linalg as nl
19	import scipy.optimize as so
20	import scipy.stats as st
21	import GSASIIpath
22	GSASIIpath.SetVersionNumber("$Revision: 1391 $")
23	import GSASIIElem as G2el
24	import GSASIIlattice as G2lat
25	import GSASIIspc as G2spc
26	import GSASIIpwd as G2pwd
27	import GSASIImapvars as G2mv
28	import GSASIImath as G2mth
29
30	sind = lambda x: np.sin(x*np.pi/180.)
31	cosd = lambda x: np.cos(x*np.pi/180.)
32	tand = lambda x: np.tan(x*np.pi/180.)
33	asind = lambda x: 180.*np.arcsin(x)/np.pi
34	acosd = lambda x: 180.*np.arccos(x)/np.pi
35	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
36
37	ateln2 = 8.0*math.log(2.0)
38
39	################################################################################
40	##### Rigid Body Models
41	################################################################################
42
43	def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
44	''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
45	current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
46	'''
47	atxIds = ['Ax:','Ay:','Az:']
48	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
49	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
50	if not RBIds['Vector'] and not RBIds['Residue']:
51	return
52	VRBIds = RBIds['Vector']
53	RRBIds = RBIds['Residue']
54	if Update:
55	RBData = rigidbodyDict
56	else:
57	RBData = copy.deepcopy(rigidbodyDict) # don't mess with original!
58	if RBIds['Vector']: # first update the vector magnitudes
59	VRBData = RBData['Vector']
60	for i,rbId in enumerate(VRBIds):
61	if VRBData[rbId]['useCount']:
62	for j in range(len(VRBData[rbId]['VectMag'])):
63	name = '::RBV;'+str(j)+':'+str(i)
64	VRBData[rbId]['VectMag'][j] = parmDict[name]
65	for phase in Phases:
66	Phase = Phases[phase]
67	General = Phase['General']
68	cell = General['Cell'][1:7]
69	Amat,Bmat = G2lat.cell2AB(cell)
70	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
71	pfx = str(Phase['pId'])+'::'
72	if Update:
73	RBModels = Phase['RBModels']
74	else:
75	RBModels = copy.deepcopy(Phase['RBModels']) # again don't mess with original!
76	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
77	jrb = VRBIds.index(RBObj['RBId'])
78	rbsx = str(irb)+':'+str(jrb)
79	for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
80	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
81	for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
82	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
83	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
84	TLS = RBObj['ThermalMotion']
85	if 'T' in TLS[0]:
86	for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
87	TLS[1][i] = parmDict[pfx+pt+rbsx]
88	if 'L' in TLS[0]:
89	for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
90	TLS[1][i+6] = parmDict[pfx+pt+rbsx]
91	if 'S' in TLS[0]:
92	for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
93	TLS[1][i+12] = parmDict[pfx+pt+rbsx]
94	if 'U' in TLS[0]:
95	TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
96	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
97	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
98	for i,x in enumerate(XYZ):
99	atId = RBObj['Ids'][i]
100	for j in [0,1,2]:
101	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
102	if UIJ[i][0] == 'A':
103	for j in range(6):
104	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
105	elif UIJ[i][0] == 'I':
106	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
107
108	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
109	jrb = RRBIds.index(RBObj['RBId'])
110	rbsx = str(irb)+':'+str(jrb)
111	for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
112	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
113	for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
114	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
115	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
116	TLS = RBObj['ThermalMotion']
117	if 'T' in TLS[0]:
118	for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
119	RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
120	if 'L' in TLS[0]:
121	for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
122	RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
123	if 'S' in TLS[0]:
124	for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
125	RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
126	if 'U' in TLS[0]:
127	RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
128	for itors,tors in enumerate(RBObj['Torsions']):
129	tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
130	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
131	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
132	for i,x in enumerate(XYZ):
133	atId = RBObj['Ids'][i]
134	for j in [0,1,2]:
135	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
136	if UIJ[i][0] == 'A':
137	for j in range(6):
138	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
139	elif UIJ[i][0] == 'I':
140	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
141
142	def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
143	'Needs a doc string'
144	atxIds = ['dAx:','dAy:','dAz:']
145	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
146	TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
147	LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
148	SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
149	PIds = ['Px:','Py:','Pz:']
150	OIds = ['Oa:','Oi:','Oj:','Ok:']
151	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
152	if not RBIds['Vector'] and not RBIds['Residue']:
153	return
154	VRBIds = RBIds['Vector']
155	RRBIds = RBIds['Residue']
156	RBData = rigidbodyDict
157	for item in parmDict:
158	if 'RB' in item:
159	dFdvDict[item] = 0. #NB: this is a vector which is no. refl. long & must be filled!
160	General = Phase['General']
161	cell = General['Cell'][1:7]
162	Amat,Bmat = G2lat.cell2AB(cell)
163	rpd = np.pi/180.
164	rpd2 = rpd**2
165	g = nl.inv(np.inner(Bmat,Bmat))
166	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
167	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
168	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
169	pfx = str(Phase['pId'])+'::'
170	RBModels = Phase['RBModels']
171	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
172	VModel = RBData['Vector'][RBObj['RBId']]
173	Q = RBObj['Orient'][0]
174	Pos = RBObj['Orig'][0]
175	jrb = VRBIds.index(RBObj['RBId'])
176	rbsx = str(irb)+':'+str(jrb)
177	dXdv = []
178	for iv in range(len(VModel['VectMag'])):
179	dCdv = []
180	for vec in VModel['rbVect'][iv]:
181	dCdv.append(G2mth.prodQVQ(Q,vec))
182	dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
183	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
184	for ia,atId in enumerate(RBObj['Ids']):
185	atNum = AtLookup[atId]
186	dx = 0.00001
187	for iv in range(len(VModel['VectMag'])):
188	for ix in [0,1,2]:
189	dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
190	for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
191	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
192	for iv in range(4):
193	Q[iv] -= dx
194	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
195	Q[iv] += 2.*dx
196	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
197	Q[iv] -= dx
198	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
199	for ix in [0,1,2]:
200	dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
201	X = G2mth.prodQVQ(Q,Cart[ia])
202	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
203	dFdu = G2lat.U6toUij(dFdu.T)
204	dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
205	dFdu = G2lat.UijtoU6(dFdu)
206	atNum = AtLookup[atId]
207	if 'T' in RBObj['ThermalMotion'][0]:
208	for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
209	dFdvDict[pfx+name+rbsx] += dFdu[i]
210	if 'L' in RBObj['ThermalMotion'][0]:
211	dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
212	dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
213	dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
214	dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
215	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
216	dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2(-dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
217	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
218	dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2(-dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
219	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
220	if 'S' in RBObj['ThermalMotion'][0]:
221	dFdvDict[pfx+'RBVS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
222	dFdvDict[pfx+'RBVS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
223	dFdvDict[pfx+'RBVS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
224	dFdvDict[pfx+'RBVS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
225	dFdvDict[pfx+'RBVS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
226	dFdvDict[pfx+'RBVS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
227	dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
228	dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
229	if 'U' in RBObj['ThermalMotion'][0]:
230	dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
231
232
233	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
234	Q = RBObj['Orient'][0]
235	Pos = RBObj['Orig'][0]
236	jrb = RRBIds.index(RBObj['RBId'])
237	torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
238	rbsx = str(irb)+':'+str(jrb)
239	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
240	for itors,tors in enumerate(RBObj['Torsions']): #derivative error?
241	tname = pfx+'RBRTr;'+str(itors)+':'+rbsx
242	orId,pvId = torData[itors][:2]
243	pivotVec = Cart[orId]-Cart[pvId]
244	QA = G2mth.AVdeg2Q(-0.001,pivotVec)
245	QB = G2mth.AVdeg2Q(0.001,pivotVec)
246	for ir in torData[itors][3]:
247	atNum = AtLookup[RBObj['Ids'][ir]]
248	rVec = Cart[ir]-Cart[pvId]
249	dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
250	dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
251	for ix in [0,1,2]:
252	dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
253	for ia,atId in enumerate(RBObj['Ids']):
254	atNum = AtLookup[atId]
255	dx = 0.00001
256	for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
257	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
258	for iv in range(4):
259	Q[iv] -= dx
260	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
261	Q[iv] += 2.*dx
262	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
263	Q[iv] -= dx
264	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
265	for ix in [0,1,2]:
266	dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
267	X = G2mth.prodQVQ(Q,Cart[ia])
268	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
269	dFdu = G2lat.U6toUij(dFdu.T)
270	dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
271	dFdu = G2lat.UijtoU6(dFdu)
272	atNum = AtLookup[atId]
273	if 'T' in RBObj['ThermalMotion'][0]:
274	for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
275	dFdvDict[pfx+name+rbsx] += dFdu[i]
276	if 'L' in RBObj['ThermalMotion'][0]:
277	dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
278	dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
279	dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
280	dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
281	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
282	dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2(dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
283	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
284	dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2(dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
285	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
286	if 'S' in RBObj['ThermalMotion'][0]:
287	dFdvDict[pfx+'RBRS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
288	dFdvDict[pfx+'RBRS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
289	dFdvDict[pfx+'RBRS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
290	dFdvDict[pfx+'RBRS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
291	dFdvDict[pfx+'RBRS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
292	dFdvDict[pfx+'RBRS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
293	dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
294	dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
295	if 'U' in RBObj['ThermalMotion'][0]:
296	dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
297
298	################################################################################
299	##### Penalty & restraint functions
300	################################################################################
301
302	def penaltyFxn(HistoPhases,parmDict,varyList):
303	'Needs a doc string'
304	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
305	pNames = []
306	pVals = []
307	pWt = []
308	negWt = {}
309	pWsum = {}
310	for phase in Phases:
311	pId = Phases[phase]['pId']
312	negWt[pId] = Phases[phase]['General']['Pawley neg wt']
313	General = Phases[phase]['General']
314	textureData = General['SH Texture']
315	SGData = General['SGData']
316	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
317	cell = General['Cell'][1:7]
318	Amat,Bmat = G2lat.cell2AB(cell)
319	if phase not in restraintDict:
320	continue
321	phaseRest = restraintDict[phase]
322	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
323	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
324	['ChemComp','Sites'],['Texture','HKLs']]
325	for name,rest in names:
326	pWsum[name] = 0.
327	itemRest = phaseRest[name]
328	if itemRest[rest] and itemRest['Use']:
329	wt = itemRest['wtFactor']
330	if name in ['Bond','Angle','Plane','Chiral']:
331	for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
332	pNames.append(str(pId)+':'+name+':'+str(i))
333	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
334	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
335	if name == 'Bond':
336	calc = G2mth.getRestDist(XYZ,Amat)
337	elif name == 'Angle':
338	calc = G2mth.getRestAngle(XYZ,Amat)
339	elif name == 'Plane':
340	calc = G2mth.getRestPlane(XYZ,Amat)
341	elif name == 'Chiral':
342	calc = G2mth.getRestChiral(XYZ,Amat)
343	pVals.append(obs-calc)
344	pWt.append(wt/esd**2)
345	pWsum[name] += wt((obs-calc)/esd)*2
346	elif name in ['Torsion','Rama']:
347	coeffDict = itemRest['Coeff']
348	for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
349	pNames.append(str(pId)+':'+name+':'+str(i))
350	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
351	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
352	if name == 'Torsion':
353	tor = G2mth.getRestTorsion(XYZ,Amat)
354	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
355	else:
356	phi,psi = G2mth.getRestRama(XYZ,Amat)
357	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
358	pVals.append(restr)
359	pWt.append(wt/esd**2)
360	pWsum[name] += wt(restr/esd)*2
361	elif name == 'ChemComp':
362	for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
363	pNames.append(str(pId)+':'+name+':'+str(i))
364	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
365	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
366	calc = np.sum(mulfracfactors)
367	pVals.append(obs-calc)
368	pWt.append(wt/esd**2)
369	pWsum[name] += wt((obs-calc)/esd)*2
370	elif name == 'Texture':
371	SHkeys = textureData['SH Coeff'][1].keys()
372	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
373	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
374	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
375	for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
376	PH = np.array(hkl)
377	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
378	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
379	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
380	Z1 = -ma.masked_greater(Z,0.0)
381	IndZ1 = np.array(ma.nonzero(Z1))
382	for ind in IndZ1.T:
383	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
384	pVals.append(Z1[ind[0]][ind[1]])
385	pWt.append(wt/esd1**2)
386	pWsum[name] += wt((obs-calc)/esd)*2
387	if ifesd2:
388	Z2 = 1.-Z
389	for ind in np.ndindex(grid,grid):
390	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
391	pVals.append(Z1[ind[0]][ind[1]])
392	pWt.append(wt/esd2**2)
393	pWsum[name] += wt((obs-calc)/esd)*2
394
395	pWsum['PWLref'] = 0.
396	for item in varyList:
397	if 'PWLref' in item and parmDict[item] < 0.:
398	pId = int(item.split(':')[0])
399	if negWt[pId]:
400	pNames.append(item)
401	pVals.append(-parmDict[item])
402	pWt.append(negWt[pId])
403	pWsum['PWLref'] += negWt[pId](-parmDict[item])*2
404	pVals = np.array(pVals)
405	pWt = np.array(pWt) #should this be np.sqrt?
406	return pNames,pVals,pWt,pWsum
407
408	def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
409	'Needs a doc string'
410	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
411	pDerv = np.zeros((len(varyList),len(pVal)))
412	for phase in Phases:
413	# if phase not in restraintDict:
414	# continue
415	pId = Phases[phase]['pId']
416	General = Phases[phase]['General']
417	SGData = General['SGData']
418	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
419	cell = General['Cell'][1:7]
420	Amat,Bmat = G2lat.cell2AB(cell)
421	textureData = General['SH Texture']
422
423	SHkeys = textureData['SH Coeff'][1].keys()
424	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
425	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
426	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
427	sam = SamSym[textureData['Model']]
428	phaseRest = restraintDict.get(phase,{})
429	names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
430	'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
431	'ChemComp':'Sites','Texture':'HKLs'}
432	lasthkl = np.array([0,0,0])
433	for ip,pName in enumerate(pNames):
434	pnames = pName.split(':')
435	if pId == int(pnames[0]):
436	name = pnames[1]
437	if 'PWL' in pName:
438	pDerv[varyList.index(pName)][ip] += 1.
439	continue
440	id = int(pnames[2])
441	itemRest = phaseRest[name]
442	if name in ['Bond','Angle','Plane','Chiral']:
443	indx,ops,obs,esd = itemRest[names[name]][id]
444	dNames = []
445	for ind in indx:
446	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
447	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
448	if name == 'Bond':
449	deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
450	elif name == 'Angle':
451	deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
452	elif name == 'Plane':
453	deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
454	elif name == 'Chiral':
455	deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
456	elif name in ['Torsion','Rama']:
457	coffDict = itemRest['Coeff']
458	indx,ops,cofName,esd = itemRest[names[name]][id]
459	dNames = []
460	for ind in indx:
461	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
462	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
463	if name == 'Torsion':
464	deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
465	else:
466	deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
467	elif name == 'ChemComp':
468	indx,factors,obs,esd = itemRest[names[name]][id]
469	dNames = []
470	for ind in indx:
471	dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
472	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
473	deriv = mul*factors
474	elif 'Texture' in name:
475	deriv = []
476	dNames = []
477	hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
478	hkl = np.array(hkl)
479	if np.any(lasthkl-hkl):
480	PH = np.array(hkl)
481	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
482	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
483	lasthkl = copy.copy(hkl)
484	if 'unit' in name:
485	pass
486	else:
487	gam = float(pnames[3])
488	psi = float(pnames[4])
489	for SHname in ODFln:
490	l,m,n = eval(SHname[1:])
491	Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
492	dNames += [str(pId)+'::'+SHname]
493	deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
494	for dName,drv in zip(dNames,deriv):
495	try:
496	ind = varyList.index(dName)
497	pDerv[ind][ip] += drv
498	except ValueError:
499	pass
500	return pDerv
501
502	################################################################################
503	##### Function & derivative calculations
504	################################################################################
505
506	def GetAtomFXU(pfx,calcControls,parmDict):
507	'Needs a doc string'
508	Natoms = calcControls['Natoms'][pfx]
509	Tdata = Natoms*[' ',]
510	Mdata = np.zeros(Natoms)
511	IAdata = Natoms*[' ',]
512	Fdata = np.zeros(Natoms)
513	FFdata = []
514	BLdata = []
515	Xdata = np.zeros((3,Natoms))
516	dXdata = np.zeros((3,Natoms))
517	Uisodata = np.zeros(Natoms)
518	Uijdata = np.zeros((6,Natoms))
519	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
520	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
521	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
522	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
523	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
524	for iatm in range(Natoms):
525	for key in keys:
526	parm = pfx+key+str(iatm)
527	if parm in parmDict:
528	keys[key][iatm] = parmDict[parm]
529	Fdata = np.where(Fdata,Fdata,1.e-8) #avoid divide by zero in derivative calc.?
530	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
531
532	def StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
533	''' Not Used: here only for comparison the StructureFactor2 - faster version
534	Compute structure factors for all h,k,l for phase
535	puts the result, F^2, in each ref[8] in refList
536	input:
537
538	:param dict refDict: where
539	'RefList' list where each ref = h,k,l,m,d,...
540	'FF' dict of form factors - filed in below
541	:param np.array G: reciprocal metric tensor
542	:param str pfx: phase id string
543	:param dict SGData: space group info. dictionary output from SpcGroup
544	:param dict calcControls:
545	:param dict ParmDict:
546
547	'''
548	twopi = 2.0*np.pi
549	twopisq = 2.0np.pi*2
550	phfx = pfx.split(':')[0]+hfx
551	ast = np.sqrt(np.diag(G))
552	Mast = twopisq*np.multiply.outer(ast,ast)
553	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
554	SGT = np.array([ops[1] for ops in SGData['SGOps']])
555	FFtables = calcControls['FFtables']
556	BLtables = calcControls['BLtables']
557	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
558	FF = np.zeros(len(Tdata))
559	if 'N' in calcControls[hfx+'histType']:
560	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
561	else:
562	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
563	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
564	Uij = np.array(G2lat.U6toUij(Uijdata))
565	bij = Mast*Uij.T
566	if not len(refDict['FF']):
567	if 'N' in calcControls[hfx+'histType']:
568	dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's
569	else:
570	dat = G2el.getFFvalues(FFtables,0.)
571	refDict['FF']['El'] = dat.keys()
572	refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
573	for iref,refl in enumerate(refDict['RefList']):
574	fbs = np.array([0,0])
575	H = refl[:3]
576	SQ = 1./(2.refl[4])*2
577	SQfactor = 4.0SQtwopisq
578	Bab = parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor)
579	if not np.any(refDict['FF']['FF'][iref]): #no form factors - 1st time thru StructureFactor
580	if 'N' in calcControls[hfx+'histType']:
581	dat = G2el.getBLvalues(BLtables)
582	refDict['FF']['FF'][iref] = dat.values()
583	else: #'X'
584	dat = G2el.getFFvalues(FFtables,SQ)
585	refDict['FF']['FF'][iref] = dat.values()
586	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
587	FF = refDict['FF']['FF'][iref][Tindx]
588	Uniq = np.inner(H,SGMT)
589	Phi = np.inner(H,SGT)
590	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
591	sinp = np.sin(phase)
592	cosp = np.cos(phase)
593	biso = -SQfactor*Uisodata
594	Tiso = np.where(biso<1.,np.exp(biso),1.0)
595	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
596	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
597	Tcorr = TisoTuijMdata*Fdata/len(Uniq)
598	fa = np.array([(FF+FP-Bab)cospTcorr,-FPPsinpTcorr])
599	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
600	if not SGData['SGInv']:
601	fb = np.array([(FF+FP-Bab)sinpTcorr,FPPcospTcorr])
602	fbs = np.sum(np.sum(fb,axis=1),axis=1)
603	fasq = fas**2
604	fbsq = fbs**2 #imaginary
605	refl[9] = np.sum(fasq)+np.sum(fbsq)
606	refl[10] = atan2d(fbs[0],fas[0])
607
608	def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
609	''' Compute structure factors for all h,k,l for phase
610	puts the result, F^2, in each ref[8] in refList
611	input:
612
613	:param dict refDict: where
614	'RefList' list where each ref = h,k,l,m,d,...
615	'FF' dict of form factors - filed in below
616	:param np.array G: reciprocal metric tensor
617	:param str pfx: phase id string
618	:param dict SGData: space group info. dictionary output from SpcGroup
619	:param dict calcControls:
620	:param dict ParmDict:
621
622	'''
623	twopi = 2.0*np.pi
624	twopisq = 2.0np.pi*2
625	phfx = pfx.split(':')[0]+hfx
626	ast = np.sqrt(np.diag(G))
627	Mast = twopisq*np.multiply.outer(ast,ast)
628	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
629	SGT = np.array([ops[1] for ops in SGData['SGOps']])
630	FFtables = calcControls['FFtables']
631	BLtables = calcControls['BLtables']
632	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
633	FF = np.zeros(len(Tdata))
634	if 'N' in calcControls[hfx+'histType']:
635	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
636	else:
637	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
638	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
639	Uij = np.array(G2lat.U6toUij(Uijdata))
640	bij = Mast*Uij.T
641	blkSize = 100 #no. of reflections in a block
642	nRef = refDict['RefList'].shape[0]
643	if not len(refDict['FF']): #no form factors - 1st time thru StructureFactor
644	if 'N' in calcControls[hfx+'histType']:
645	dat = G2el.getBLvalues(BLtables)
646	refDict['FF']['El'] = dat.keys()
647	refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
648	else: #'X'
649	dat = G2el.getFFvalues(FFtables,0.)
650	refDict['FF']['El'] = dat.keys()
651	refDict['FF']['FF'] = np.ones((nRef,len(dat)))
652	for iref,ref in enumerate(refDict['RefList']):
653	SQ = 1./(2.ref[4])*2
654	dat = G2el.getFFvalues(FFtables,SQ)
655	refDict['FF']['FF'][iref] *= dat.values()
656	#reflection processing begins here - big arrays!
657	iBeg = 0
658	while iBeg < nRef:
659	iFin = min(iBeg+blkSize,nRef)
660	refl = refDict['RefList'][iBeg:iFin]
661	H = refl.T[:3]
662	SQ = 1./(2.refl.T[4])*2
663	SQfactor = 4.0SQtwopisq
664	Bab = np.repeat(parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor),len(SGT))
665	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
666	FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
667	Uniq = np.reshape(np.inner(H.T,SGMT),(-1,3))
668	Phi = np.inner(H.T,SGT).flatten()
669	phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T)+Phi[:,np.newaxis])
670	sinp = np.sin(phase)
671	cosp = np.cos(phase)
672	biso = -SQfactor*Uisodata[:,np.newaxis]
673	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
674	HbH = -np.sum(Uniq.T*np.inner(bij,Uniq),axis=1)
675	Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
676	Tcorr = TisoTuijMdata*Fdata/len(SGMT)
677	fa = np.array([((FF+FP).T-Bab).TcospTcorr,-FPPsinpTcorr])
678	fa = np.reshape(fa,(2,len(refl),len(SGT),len(Mdata)))
679	fas = np.sum(np.sum(fa,axis=2),axis=2) #real
680	fbs = np.zeros_like(fas)
681	if not SGData['SGInv']:
682	fb = np.array([((FF+FP).T-Bab).TsinpTcorr,FPPcospTcorr])
683	fb = np.reshape(fb,(2,len(refl),len(SGT),len(Mdata)))
684	fbs = np.sum(np.sum(fb,axis=2),axis=2)
685	fasq = fas**2
686	fbsq = fbs**2 #imaginary
687	refl.T[9] = np.sum(fasq,axis=0)+np.sum(fbsq,axis=0)
688	refl.T[10] = atan2d(fbs[0],fas[0])
689	iBeg += blkSize
690
691	def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
692	'Needs a doc string'
693	twopi = 2.0*np.pi
694	twopisq = 2.0np.pi*2
695	phfx = pfx.split(':')[0]+hfx
696	ast = np.sqrt(np.diag(G))
697	Mast = twopisq*np.multiply.outer(ast,ast)
698	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
699	SGT = np.array([ops[1] for ops in SGData['SGOps']])
700	FFtables = calcControls['FFtables']
701	BLtables = calcControls['BLtables']
702	nRef = len(refDict['RefList'])
703	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
704	mSize = len(Mdata)
705	FF = np.zeros(len(Tdata))
706	if 'N' in calcControls[hfx+'histType']:
707	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
708	else:
709	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
710	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
711	Uij = np.array(G2lat.U6toUij(Uijdata))
712	bij = Mast*Uij.T
713	dFdvDict = {}
714	dFdfr = np.zeros((nRef,mSize))
715	dFdx = np.zeros((nRef,mSize,3))
716	dFdui = np.zeros((nRef,mSize))
717	dFdua = np.zeros((nRef,mSize,6))
718	dFdbab = np.zeros((nRef,2))
719	for iref,refl in enumerate(refDict['RefList']):
720	H = np.array(refl[:3])
721	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
722	SQfactor = 8.0SQnp.pi**2
723	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
724	Bab = parmDict[phfx+'BabA']*dBabdA
725	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
726	FF = refDict['FF']['FF'][iref].T[Tindx]
727	# FF = [refDict['FF'][iref][El] for El in Tdata]
728	Uniq = np.inner(H,SGMT)
729	Phi = np.inner(H,SGT)
730	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
731	sinp = np.sin(phase)
732	cosp = np.cos(phase)
733	occ = Mdata*Fdata/len(Uniq)
734	biso = -SQfactor*Uisodata
735	Tiso = np.where(biso<1.,np.exp(biso),1.0)
736	HbH = -np.inner(H,np.inner(bij,H))
737	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
738	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
739	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
740	Tcorr = Tiso*Tuij
741	fot = (FF+FP-Bab)occTcorr
742	fotp = FPPoccTcorr
743	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
744	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
745
746	fas = np.sum(np.sum(fa,axis=1),axis=1)
747	fbs = np.sum(np.sum(fb,axis=1),axis=1)
748	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
749	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
750	#sum below is over Uniq
751	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2) #Fdata != 0 ever avoids /0. problem
752	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
753	dfadui = np.sum(-SQfactor*fa,axis=2)
754	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
755	dfadba = np.sum(-cosp(occTcorr)[:,np.newaxis],axis=1)
756	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for
757	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
758	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
759	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
760	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
761	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
762	if not SGData['SGInv']:
763	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #Fdata != 0 ever avoids /0. problem
764	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
765	dfbdui = np.sum(-SQfactor*fb,axis=2)
766	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
767	dfbdba = np.sum(-sinp(occTcorr)[:,np.newaxis],axis=1)
768	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
769	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
770	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
771	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
772	dFdbab[iref] += 2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
773	#loop over atoms - each dict entry is list of derivatives for all the reflections
774	for i in range(len(Mdata)):
775	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
776	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
777	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
778	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
779	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
780	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
781	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
782	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
783	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
784	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
785	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
786	dFdvDict[pfx+'BabA'] = dFdbab.T[0]
787	dFdvDict[pfx+'BabU'] = dFdbab.T[1]
788	return dFdvDict
789
790	def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
791	''' Single crystal extinction function; puts correction in ref[11] and returns
792	corrections needed for derivatives
793	'''
794	ref[11] = 1.0
795	dervCor = 1.0
796	dervDict = {}
797	if calcControls[phfx+'EType'] != 'None':
798	cos2T = 1.0-0.5(parmDict[hfx+'Lam']/ref[4])*2 #cos(2theta)
799	if 'SXC' in parmDict[hfx+'Type']:
800	AV = 7.9406e5/parmDict[pfx+'Vol']**2
801	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
802	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
803	elif 'SNT' in parmDict[hfx+'Type']:
804	AV = 1.e7/parmDict[pfx+'Vol']**2
805	PL = 1./(4.refl[4]*2)
806	P12 = 1.0
807	elif 'SNC' in parmDict[hfx+'Type']:
808	AV = 1.e7/parmDict[pfx+'Vol']**2
809	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
810	P12 = 1.0
811
812	PLZ = AVP12parmDict[hfx+'Lam']*2ref[7]
813	if 'Primary' in calcControls[phfx+'EType']:
814	PLZ *= 1.5
815	else:
816	PLZ *= calcControls[phfx+'Tbar']
817
818	if 'Primary' in calcControls[phfx+'EType']:
819	PSIG = parmDict[phfx+'Ep']
820	elif 'I & II' in calcControls[phfx+'EType']:
821	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
822	elif 'Type II' in calcControls[phfx+'EType']:
823	PSIG = parmDict[phfx+'Es']
824	else: # 'Secondary Type I'
825	PSIG = parmDict[phfx+'Eg']/PL
826
827	AG = 0.58+0.48cos2T+0.24cos2T**2
828	AL = 0.025+0.285*cos2T
829	BG = 0.02-0.025*cos2T
830	BL = 0.15-0.2(0.75-cos2T)*2
831	if cos2T < 0.:
832	BL = -0.45*cos2T
833	CG = 2.
834	CL = 2.
835	PF = PLZ*PSIG
836
837	if 'Gaussian' in calcControls[phfx+'EApprox']:
838	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
839	extCor = np.sqrt(PF4)
840	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
841	else:
842	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
843	extCor = np.sqrt(PF4)
844	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
845
846	dervCor = (1.+PF)*PF3
847	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
848	dervDict[phfx+'Ep'] = -ref[7]PLZPF3
849	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
850	dervDict[phfx+'Es'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
851	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
852	dervDict[phfx+'Eg'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
853
854	ref[11] = 1./extCor
855	return dervCor,dervDict
856
857	def Dict2Values(parmdict, varylist):
858	'''Use before call to leastsq to setup list of values for the parameters
859	in parmdict, as selected by key in varylist'''
860	return [parmdict[key] for key in varylist]
861
862	def Values2Dict(parmdict, varylist, values):
863	''' Use after call to leastsq to update the parameter dictionary with
864	values corresponding to keys in varylist'''
865	parmdict.update(zip(varylist,values))
866
867	def GetNewCellParms(parmDict,varyList):
868	'Needs a doc string'
869	newCellDict = {}
870	Anames = ['A'+str(i) for i in range(6)]
871	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
872	for item in varyList:
873	keys = item.split(':')
874	if keys[2] in Ddict:
875	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
876	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
877	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
878	return newCellDict # is e.g. {'0::D11':A0+D11}
879
880	def ApplyXYZshifts(parmDict,varyList):
881	'''
882	takes atom x,y,z shift and applies it to corresponding atom x,y,z value
883
884	:param dict parmDict: parameter dictionary
885	:param list varyList: list of variables (not used!)
886	:returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
887
888	'''
889	newAtomDict = {}
890	for item in parmDict:
891	if 'dA' in item:
892	parm = ''.join(item.split('d'))
893	parmDict[parm] += parmDict[item]
894	newAtomDict[item] = [parm,parmDict[parm]]
895	return newAtomDict
896
897	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
898	'Spherical harmonics texture'
899	IFCoup = 'Bragg' in calcControls[hfx+'instType']
900	odfCor = 1.0
901	H = refl[:3]
902	cell = G2lat.Gmat2cell(g)
903	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
904	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
905	phi,beta = G2lat.CrsAng(H,cell,SGData)
906	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
907	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
908	for item in SHnames:
909	L,M,N = eval(item.strip('C'))
910	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
911	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
912	Lnorm = G2lat.Lnorm(L)
913	odfCor += parmDict[pfx+item]LnormKcl*Ksl
914	return odfCor
915
916	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
917	'Spherical harmonics texture derivatives'
918	FORPI = 4.0*np.pi
919	IFCoup = 'Bragg' in calcControls[hfx+'instType']
920	odfCor = 1.0
921	dFdODF = {}
922	dFdSA = [0,0,0]
923	H = refl[:3]
924	cell = G2lat.Gmat2cell(g)
925	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
926	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
927	phi,beta = G2lat.CrsAng(H,cell,SGData)
928	psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
929	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
930	for item in SHnames:
931	L,M,N = eval(item.strip('C'))
932	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
933	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
934	Lnorm = G2lat.Lnorm(L)
935	odfCor += parmDict[pfx+item]LnormKcl*Ksl
936	dFdODF[pfx+item] = LnormKclKsl
937	for i in range(3):
938	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
939	return odfCor,dFdODF,dFdSA
940
941	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
942	'spherical harmonics preferred orientation (cylindrical symmetry only)'
943	odfCor = 1.0
944	H = refl[:3]
945	cell = G2lat.Gmat2cell(g)
946	Sangl = [0.,0.,0.]
947	if 'Bragg' in calcControls[hfx+'instType']:
948	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
949	IFCoup = True
950	else:
951	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
952	IFCoup = False
953	phi,beta = G2lat.CrsAng(H,cell,SGData)
954	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
955	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
956	for item in SHnames:
957	L,N = eval(item.strip('C'))
958	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
959	odfCor += parmDict[phfx+item]LnormKcsl
960	return np.squeeze(odfCor)
961
962	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
963	'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
964	FORPI = 12.5663706143592
965	odfCor = 1.0
966	dFdODF = {}
967	H = refl[:3]
968	cell = G2lat.Gmat2cell(g)
969	Sangl = [0.,0.,0.]
970	if 'Bragg' in calcControls[hfx+'instType']:
971	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
972	IFCoup = True
973	else:
974	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
975	IFCoup = False
976	phi,beta = G2lat.CrsAng(H,cell,SGData)
977	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
978	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
979	for item in SHnames:
980	L,N = eval(item.strip('C'))
981	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
982	odfCor += parmDict[phfx+item]LnormKcsl
983	dFdODF[phfx+item] = Kcsl*Lnorm
984	return odfCor,dFdODF
985
986	def GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
987	'Needs a doc string'
988	POcorr = 1.0
989	if calcControls[phfx+'poType'] == 'MD':
990	MD = parmDict[phfx+'MD']
991	if MD != 1.0:
992	MDAxis = calcControls[phfx+'MDAxis']
993	sumMD = 0
994	for H in uniq:
995	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
996	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
997	sumMD += A**3
998	POcorr = sumMD/len(uniq)
999	else: #spherical harmonics
1000	if calcControls[phfx+'SHord']:
1001	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1002	return POcorr
1003
1004	def GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1005	'Needs a doc string'
1006	POcorr = 1.0
1007	POderv = {}
1008	if calcControls[phfx+'poType'] == 'MD':
1009	MD = parmDict[phfx+'MD']
1010	MDAxis = calcControls[phfx+'MDAxis']
1011	sumMD = 0
1012	sumdMD = 0
1013	for H in uniq:
1014	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1015	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1016	sumMD += A**3
1017	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
1018	POcorr = sumMD/len(uniq)
1019	POderv[phfx+'MD'] = sumdMD/len(uniq)
1020	else: #spherical harmonics
1021	if calcControls[phfx+'SHord']:
1022	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1023	return POcorr,POderv
1024
1025	def GetAbsorb(refl,hfx,calcControls,parmDict):
1026	'Needs a doc string'
1027	if 'Debye' in calcControls[hfx+'instType']:
1028	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1029	else:
1030	return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1031
1032	def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
1033	'Needs a doc string'
1034	if 'Debye' in calcControls[hfx+'instType']:
1035	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1036	else:
1037	return G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1038
1039	def GetIntensityCorr(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1040	'Needs a doc string' #need powder extinction!
1041	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
1042	if 'X' in parmDict[hfx+'Type']:
1043	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1044	Icorr *= GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1045	if pfx+'SHorder' in parmDict:
1046	Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1047	Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
1048	refl[11] = Icorr
1049
1050	def GetIntensityDerv(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1051	'Needs a doc string' #need powder extinction derivs!
1052	dIdsh = 1./parmDict[hfx+'Scale']
1053	dIdsp = 1./parmDict[phfx+'Scale']
1054	if 'X' in parmDict[hfx+'Type']:
1055	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1056	dIdPola /= pola
1057	else: #'N'
1058	dIdPola = 0.0
1059	POcorr,dIdPO = GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1060	for iPO in dIdPO:
1061	dIdPO[iPO] /= POcorr
1062	dFdODF = {}
1063	dFdSA = [0,0,0]
1064	if pfx+'SHorder' in parmDict:
1065	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1066	for iSH in dFdODF:
1067	dFdODF[iSH] /= odfCor
1068	for i in range(3):
1069	dFdSA[i] /= odfCor
1070	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
1071	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
1072
1073	def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
1074	'Needs a doc string'
1075	costh = cosd(refl[5]/2.)
1076	#crystallite size
1077	if calcControls[phfx+'SizeType'] == 'isotropic':
1078	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1079	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1080	H = np.array(refl[:3])
1081	P = np.array(calcControls[phfx+'SizeAxis'])
1082	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1083	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
1084	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1085	else: #ellipsoidal crystallites
1086	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1087	H = np.array(refl[:3])
1088	lenR = G2pwd.ellipseSize(H,Sij,GB)
1089	Sgam = 1.8wave/(np.picosth*lenR)
1090	#microstrain
1091	if calcControls[phfx+'MustrainType'] == 'isotropic':
1092	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
1093	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1094	H = np.array(refl[:3])
1095	P = np.array(calcControls[phfx+'MustrainAxis'])
1096	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1097	Si = parmDict[phfx+'Mustrain;i']
1098	Sa = parmDict[phfx+'Mustrain;a']
1099	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
1100	else: #generalized - P.W. Stephens model
1101	pwrs = calcControls[phfx+'MuPwrs']
1102	sum = 0
1103	for i,pwr in enumerate(pwrs):
1104	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1105	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
1106	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
1107	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
1108	sig /= ateln2
1109	return sig,gam
1110
1111	def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
1112	'Needs a doc string'
1113	gamDict = {}
1114	sigDict = {}
1115	costh = cosd(refl[5]/2.)
1116	tanth = tand(refl[5]/2.)
1117	#crystallite size derivatives
1118	if calcControls[phfx+'SizeType'] == 'isotropic':
1119	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1120	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
1121	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
1122	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1123	H = np.array(refl[:3])
1124	P = np.array(calcControls[phfx+'SizeAxis'])
1125	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1126	Si = parmDict[phfx+'Size;i']
1127	Sa = parmDict[phfx+'Size;a']
1128	gami = (1.8wave/np.pi)/(SiSa)
1129	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1130	Sgam = gami*sqtrm
1131	gam = Sgam/costh
1132	dsi = (gamiSicosP*2/(sqtrmcosth)-gam/Si)
1133	dsa = (gamiSasinP*2/(sqtrmcosth)-gam/Sa)
1134	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1135	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1136	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1137	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1138	else: #ellipsoidal crystallites
1139	const = 1.8wave/(np.picosth)
1140	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1141	H = np.array(refl[:3])
1142	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1143	Sgam = 1.8wave/(np.picosth*lenR)
1144	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1145	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1146	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1147	gamDict[phfx+'Size;mx'] = Sgam
1148	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
1149
1150	#microstrain derivatives
1151	if calcControls[phfx+'MustrainType'] == 'isotropic':
1152	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
1153	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
1154	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
1155	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1156	H = np.array(refl[:3])
1157	P = np.array(calcControls[phfx+'MustrainAxis'])
1158	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1159	Si = parmDict[phfx+'Mustrain;i']
1160	Sa = parmDict[phfx+'Mustrain;a']
1161	gami = 0.018SiSa*tanth/np.pi
1162	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1163	Mgam = gami/sqtrm
1164	dsi = -gamiSicosP2/sqtrm3
1165	dsa = -gamiSasinP2/sqtrm3
1166	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1167	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1168	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1169	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1170	else: #generalized - P.W. Stephens model
1171	pwrs = calcControls[phfx+'MuPwrs']
1172	const = 0.018refl[4]2tanth
1173	sum = 0
1174	for i,pwr in enumerate(pwrs):
1175	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1176	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1177	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain;mx']
1178	sigDict[phfx+'Mustrain:'+str(i)] = \
1179	2.constterm(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1180	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
1181	for i in range(len(pwrs)):
1182	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1183	gamDict[phfx+'Mustrain;mx'] = Mgam
1184	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1185	return sigDict,gamDict
1186
1187	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1188	'Needs a doc string'
1189	h,k,l = refl[:3]
1190	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1191	d = np.sqrt(dsq)
1192
1193	refl[4] = d
1194	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
1195	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1196	if 'Bragg' in calcControls[hfx+'instType']:
1197	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1198	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
1199	else: #Debye-Scherrer - simple but maybe not right
1200	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1201	return pos
1202
1203	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1204	'Needs a doc string'
1205	dpr = 180./np.pi
1206	h,k,l = refl[:3]
1207	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1208	dst = np.sqrt(dstsq)
1209	pos = refl[5]-parmDict[hfx+'Zero']
1210	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
1211	dpdw = const*dst
1212	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
1213	dpdA = constwave/(2.0*dst)
1214	dpdZ = 1.0
1215	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1216	if 'Bragg' in calcControls[hfx+'instType']:
1217	dpdSh = -4.constcosd(pos/2.0)
1218	dpdTr = -constsind(pos)100.0
1219	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1220	else: #Debye-Scherrer - simple but maybe not right
1221	dpdXd = -const*cosd(pos)
1222	dpdYd = -const*sind(pos)
1223	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1224
1225	def GetHStrainShift(refl,SGData,phfx,parmDict):
1226	'Needs a doc string'
1227	laue = SGData['SGLaue']
1228	uniq = SGData['SGUniq']
1229	h,k,l = refl[:3]
1230	if laue in ['m3','m3m']:
1231	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
1232	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
1233	elif laue in ['6/m','6/mmm','3m1','31m','3']:
1234	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
1235	elif laue in ['3R','3mR']:
1236	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
1237	elif laue in ['4/m','4/mmm']:
1238	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
1239	elif laue in ['mmm']:
1240	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
1241	elif laue in ['2/m']:
1242	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
1243	if uniq == 'a':
1244	Dij += parmDict[phfx+'D23']kl
1245	elif uniq == 'b':
1246	Dij += parmDict[phfx+'D13']hl
1247	elif uniq == 'c':
1248	Dij += parmDict[phfx+'D12']hk
1249	else:
1250	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
1251	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
1252	return -Dijrefl[4]2tand(refl[5]/2.0)
1253
1254	def GetHStrainShiftDerv(refl,SGData,phfx):
1255	'Needs a doc string'
1256	laue = SGData['SGLaue']
1257	uniq = SGData['SGUniq']
1258	h,k,l = refl[:3]
1259	if laue in ['m3','m3m']:
1260	dDijDict = {phfx+'D11':h2+k2+l**2,
1261	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
1262	elif laue in ['6/m','6/mmm','3m1','31m','3']:
1263	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
1264	elif laue in ['3R','3mR']:
1265	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
1266	elif laue in ['4/m','4/mmm']:
1267	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
1268	elif laue in ['mmm']:
1269	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
1270	elif laue in ['2/m']:
1271	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
1272	if uniq == 'a':
1273	dDijDict[phfx+'D23'] = k*l
1274	elif uniq == 'b':
1275	dDijDict[phfx+'D13'] = h*l
1276	elif uniq == 'c':
1277	dDijDict[phfx+'D12'] = h*k
1278	names.append()
1279	else:
1280	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
1281	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
1282	for item in dDijDict:
1283	dDijDict[item] = -refl[4]2tand(refl[5]/2.0)
1284	return dDijDict
1285
1286	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1287	'Needs a doc string'
1288	histoList = Histograms.keys()
1289	histoList.sort()
1290	for histogram in histoList:
1291	if 'PWDR' in histogram[:4]:
1292	Histogram = Histograms[histogram]
1293	hId = Histogram['hId']
1294	hfx = ':%d:'%(hId)
1295	Limits = calcControls[hfx+'Limits']
1296	shl = max(parmDict[hfx+'SH/L'],0.0005)
1297	Ka2 = False
1298	kRatio = 0.0
1299	if hfx+'Lam1' in parmDict.keys():
1300	Ka2 = True
1301	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1302	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1303	x,y,w,yc,yb,yd = Histogram['Data']
1304	xB = np.searchsorted(x,Limits[0])
1305	xF = np.searchsorted(x,Limits[1])
1306	ymb = np.array(y-yb)
1307	ymb = np.where(ymb,ymb,1.0)
1308	ycmb = np.array(yc-yb)
1309	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
1310	refLists = Histogram['Reflection Lists']
1311	for phase in refLists:
1312	Phase = Phases[phase]
1313	pId = Phase['pId']
1314	phfx = '%d:%d:'%(pId,hId)
1315	refDict = refLists[phase]
1316	sumFo = 0.0
1317	sumdF = 0.0
1318	sumFosq = 0.0
1319	sumdFsq = 0.0
1320	for refl in refDict['RefList']:
1321	if 'C' in calcControls[hfx+'histType']:
1322	yp = np.zeros_like(yb)
1323	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1324	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1325	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1326	iFin2 = iFin
1327	yp[iBeg:iFin] = refl[11]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
1328	if Ka2:
1329	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
1330	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1331	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1332	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1333	if not iBeg2+iFin2: #peak below low limit - skip peak
1334	continue
1335	elif not iBeg2-iFin2: #peak above high limit - done
1336	break
1337	yp[iBeg2:iFin2] += refl[11]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here
1338	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[11](1.+kRatio)),0.0))
1339	elif 'T' in calcControls[hfx+'histType']:
1340	print 'TOF Undefined at present'
1341	raise Exception #no TOF yet
1342	Fo = np.sqrt(np.abs(refl[8]))
1343	Fc = np.sqrt(np.abs(refl[9]))
1344	sumFo += Fo
1345	sumFosq += refl[8]**2
1346	sumdF += np.abs(Fo-Fc)
1347	sumdFsq += (refl[8]-refl[9])**2
1348	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1349	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1350	Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1351	Histogram['Residuals']['hId'] = hId
1352	elif 'HKLF' in histogram[:4]:
1353	Histogram = Histograms[histogram]
1354	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1355
1356	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1357	'Needs a doc string'
1358
1359	def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1360	U = parmDict[hfx+'U']
1361	V = parmDict[hfx+'V']
1362	W = parmDict[hfx+'W']
1363	X = parmDict[hfx+'X']
1364	Y = parmDict[hfx+'Y']
1365	tanPos = tand(refl[5]/2.0)
1366	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
1367	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
1368	sig = max(0.001,sig)
1369	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
1370	gam = max(0.001,gam)
1371	return sig,gam
1372
1373	hId = Histogram['hId']
1374	hfx = ':%d:'%(hId)
1375	bakType = calcControls[hfx+'bakType']
1376	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1377	yc = np.zeros_like(yb)
1378
1379	if 'C' in calcControls[hfx+'histType']:
1380	shl = max(parmDict[hfx+'SH/L'],0.002)
1381	Ka2 = False
1382	if hfx+'Lam1' in parmDict.keys():
1383	wave = parmDict[hfx+'Lam1']
1384	Ka2 = True
1385	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1386	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1387	else:
1388	wave = parmDict[hfx+'Lam']
1389	else:
1390	print 'TOF Undefined at present'
1391	raise ValueError
1392	for phase in Histogram['Reflection Lists']:
1393	refDict = Histogram['Reflection Lists'][phase]
1394	Phase = Phases[phase]
1395	pId = Phase['pId']
1396	pfx = '%d::'%(pId)
1397	phfx = '%d:%d:'%(pId,hId)
1398	hfx = ':%d:'%(hId)
1399	SGData = Phase['General']['SGData']
1400	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1401	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] #Do I want to modify by Dij?
1402	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1403	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
1404	Vst = np.sqrt(nl.det(G)) #V*
1405	if not Phase['General'].get('doPawley'):
1406	time0 = time.time()
1407	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1408	# print 'sf calc time: %.3fs'%(time.time()-time0)
1409	time0 = time.time()
1410	for iref,refl in enumerate(refDict['RefList']):
1411	if 'C' in calcControls[hfx+'histType']:
1412	h,k,l = refl[:3]
1413	Uniq = np.inner(refl[:3],SGMT)
1414	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
1415	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
1416	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
1417	refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam
1418	GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[11]
1419	refl[11] = VstLorenz
1420	if Phase['General'].get('doPawley'):
1421	try:
1422	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1423	refl[9] = parmDict[pInd]
1424	except KeyError:
1425	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
1426	continue
1427	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1428	iBeg = np.searchsorted(x,refl[5]-fmin)
1429	iFin = np.searchsorted(x,refl[5]+fmax)
1430	if not iBeg+iFin: #peak below low limit - skip peak
1431	continue
1432	elif not iBeg-iFin: #peak above high limit - done
1433	break
1434	yc[iBeg:iFin] += refl[11]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
1435	if Ka2:
1436	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
1437	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1438	iBeg = np.searchsorted(x,pos2-fmin)
1439	iFin = np.searchsorted(x,pos2+fmax)
1440	if not iBeg+iFin: #peak below low limit - skip peak
1441	continue
1442	elif not iBeg-iFin: #peak above high limit - done
1443	return yc,yb
1444	yc[iBeg:iFin] += refl[11]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #and here
1445	elif 'T' in calcControls[hfx+'histType']:
1446	print 'TOF Undefined at present'
1447	raise Exception #no TOF yet
1448	# print 'profile calc time: %.3fs'%(time.time()-time0)
1449	return yc,yb
1450
1451	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1452	'Needs a doc string'
1453
1454	def cellVaryDerv(pfx,SGData,dpdA):
1455	if SGData['SGLaue'] in ['-1',]:
1456	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1457	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1458	elif SGData['SGLaue'] in ['2/m',]:
1459	if SGData['SGUniq'] == 'a':
1460	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1461	elif SGData['SGUniq'] == 'b':
1462	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1463	else:
1464	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1465	elif SGData['SGLaue'] in ['mmm',]:
1466	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1467	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1468	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1469	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1470	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1471	elif SGData['SGLaue'] in ['3R', '3mR']:
1472	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
1473	elif SGData['SGLaue'] in ['m3m','m3']:
1474	return [[pfx+'A0',dpdA[0]]]
1475
1476	# create a list of dependent variables and set up a dictionary to hold their derivatives
1477	dependentVars = G2mv.GetDependentVars()
1478	depDerivDict = {}
1479	for j in dependentVars:
1480	depDerivDict[j] = np.zeros(shape=(len(x)))
1481	#print 'dependent vars',dependentVars
1482	lenX = len(x)
1483	hId = Histogram['hId']
1484	hfx = ':%d:'%(hId)
1485	bakType = calcControls[hfx+'bakType']
1486	dMdv = np.zeros(shape=(len(varylist),len(x)))
1487	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1488	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1489	bBpos =varylist.index(hfx+'Back:0')
1490	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1491	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1492	for name in varylist:
1493	if 'Debye' in name:
1494	id = int(name.split(':')[-1])
1495	parm = name[:int(name.rindex(':'))]
1496	ip = names.index(parm)
1497	dMdv[varylist.index(name)] = dMddb[3*id+ip]
1498	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1499	for name in varylist:
1500	if 'BkPk' in name:
1501	parm,id = name.split(';')
1502	id = int(id)
1503	if parm in names:
1504	ip = names.index(parm)
1505	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1506	cw = np.diff(x)
1507	cw = np.append(cw,cw[-1])
1508	if 'C' in calcControls[hfx+'histType']:
1509	shl = max(parmDict[hfx+'SH/L'],0.002)
1510	Ka2 = False
1511	if hfx+'Lam1' in parmDict.keys():
1512	wave = parmDict[hfx+'Lam1']
1513	Ka2 = True
1514	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1515	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1516	else:
1517	wave = parmDict[hfx+'Lam']
1518	else:
1519	print 'TOF Undefined at present'
1520	raise ValueError
1521	for phase in Histogram['Reflection Lists']:
1522	refDict = Histogram['Reflection Lists'][phase]
1523	Phase = Phases[phase]
1524	SGData = Phase['General']['SGData']
1525	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1526	pId = Phase['pId']
1527	pfx = '%d::'%(pId)
1528	phfx = '%d:%d:'%(pId,hId)
1529	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] #And modify here by Dij?
1530	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1531	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
1532	if not Phase['General'].get('doPawley'):
1533	time0 = time.time()
1534	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1535	# print 'sf-derv time %.3fs'%(time.time()-time0)
1536	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1537	time0 = time.time()
1538	for iref,refl in enumerate(refDict['RefList']):
1539	if 'C' in calcControls[hfx+'histType']: #CW powder
1540	h,k,l = refl[:3]
1541	Uniq = np.inner(refl[:3],SGMT)
1542	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1543	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1544	iBeg = np.searchsorted(x,refl[5]-fmin)
1545	iFin = np.searchsorted(x,refl[5]+fmax)
1546	if not iBeg+iFin: #peak below low limit - skip peak
1547	continue
1548	elif not iBeg-iFin: #peak above high limit - done
1549	break
1550	pos = refl[5]
1551	tanth = tand(pos/2.0)
1552	costh = cosd(pos/2.0)
1553	lenBF = iFin-iBeg
1554	dMdpk = np.zeros(shape=(6,lenBF))
1555	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1556	for i in range(5):
1557	dMdpk[i] += 100.cw[iBeg:iFin]refl[11]refl[9]dMdipk[i]
1558	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1559	if Ka2:
1560	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
1561	iBeg2 = np.searchsorted(x,pos2-fmin)
1562	iFin2 = np.searchsorted(x,pos2+fmax)
1563	if iBeg2-iFin2:
1564	lenBF2 = iFin2-iBeg2
1565	dMdpk2 = np.zeros(shape=(6,lenBF2))
1566	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1567	for i in range(5):
1568	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[11]refl[9]kRatio*dMdipk2[i]
1569	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[11]*dMdipk2[0]
1570	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1571	if Phase['General'].get('doPawley'):
1572	dMdpw = np.zeros(len(x))
1573	try:
1574	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1575	idx = varylist.index(pIdx)
1576	dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1577	if Ka2:
1578	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1579	dMdv[idx] = dMdpw
1580	except: # ValueError:
1581	pass
1582	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1583	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1584	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1585	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1586	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1587	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1588	hfx+'DisplaceY':[dpdY,'pos'],}
1589	if 'Bragg' in calcControls[hfx+'instType']:
1590	names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1591	hfx+'SurfRoughB':[dFdAb[1],'int'],})
1592	else:
1593	names.update({hfx+'Absorption':[dFdAb,'int'],})
1594	for name in names:
1595	item = names[name]
1596	if name in varylist:
1597	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1598	if Ka2:
1599	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1600	elif name in dependentVars:
1601	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1602	if Ka2:
1603	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1604	for iPO in dIdPO:
1605	if iPO in varylist:
1606	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1607	if Ka2:
1608	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1609	elif iPO in dependentVars:
1610	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1611	if Ka2:
1612	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1613	for i,name in enumerate(['omega','chi','phi']):
1614	aname = pfx+'SH '+name
1615	if aname in varylist:
1616	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1617	if Ka2:
1618	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1619	elif aname in dependentVars:
1620	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1621	if Ka2:
1622	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1623	for iSH in dFdODF:
1624	if iSH in varylist:
1625	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1626	if Ka2:
1627	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1628	elif iSH in dependentVars:
1629	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1630	if Ka2:
1631	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1632	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1633	for name,dpdA in cellDervNames:
1634	if name in varylist:
1635	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1636	if Ka2:
1637	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1638	elif name in dependentVars:
1639	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1640	if Ka2:
1641	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1642	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1643	for name in dDijDict:
1644	if name in varylist:
1645	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1646	if Ka2:
1647	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1648	elif name in dependentVars:
1649	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1650	if Ka2:
1651	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1652	sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
1653	for name in gamDict:
1654	if name in varylist:
1655	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1656	if Ka2:
1657	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1658	elif name in dependentVars:
1659	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1660	if Ka2:
1661	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1662	for name in sigDict:
1663	if name in varylist:
1664	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1665	if Ka2:
1666	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1667	elif name in dependentVars:
1668	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1669	if Ka2:
1670	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1671	for name in ['BabA','BabU']:
1672	if refl[9]:
1673	if phfx+name in varylist:
1674	dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1675	if Ka2:
1676	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1677	elif phfx+name in dependentVars:
1678	depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1679	if Ka2:
1680	depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1681	elif 'T' in calcControls[hfx+'histType']:
1682	print 'TOF Undefined at present'
1683	raise Exception #no TOF yet
1684	if not Phase['General'].get('doPawley'):
1685	#do atom derivatives - for RB,F,X & U so far
1686	corr = dervDict['int']/refl[9]
1687	if Ka2:
1688	corr2 = dervDict2['int']/refl[9]
1689	for name in varylist+dependentVars:
1690	if '::RBV;' in name:
1691	pass
1692	else:
1693	try:
1694	aname = name.split(pfx)[1][:2]
1695	if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1696	except IndexError:
1697	continue
1698	if name in varylist:
1699	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1700	if Ka2:
1701	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1702	elif name in dependentVars:
1703	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1704	if Ka2:
1705	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1706	# print 'profile derv time: %.3fs'%(time.time()-time0)
1707	# now process derivatives in constraints
1708	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1709	return dMdv
1710
1711	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1712	'''Loop over histograms and compute derivatives of the fitting
1713	model (M) with respect to all parameters. Results are returned in
1714	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
1715	n is the number of parameters and m is the number of data
1716	points. This can exceed memory when m gets large. This routine is
1717	used when refinement derivatives are selected as "analtytic
1718	Jacobian" in Controls.
1719
1720	:returns: Jacobian numpy.array dMdv for all histograms concatinated
1721	'''
1722	parmDict.update(zip(varylist,values))
1723	G2mv.Dict2Map(parmDict,varylist)
1724	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1725	nvar = len(varylist)
1726	dMdv = np.empty(0)
1727	histoList = Histograms.keys()
1728	histoList.sort()
1729	for histogram in histoList:
1730	if 'PWDR' in histogram[:4]:
1731	Histogram = Histograms[histogram]
1732	hId = Histogram['hId']
1733	hfx = ':%d:'%(hId)
1734	wtFactor = calcControls[hfx+'wtFactor']
1735	Limits = calcControls[hfx+'Limits']
1736	x,y,w,yc,yb,yd = Histogram['Data']
1737	W = wtFactor*w
1738	xB = np.searchsorted(x,Limits[0])
1739	xF = np.searchsorted(x,Limits[1])
1740	dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
1741	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1742	elif 'HKLF' in histogram[:4]:
1743	Histogram = Histograms[histogram]
1744	nobs = Histogram['Residuals']['Nobs']
1745	phase = Histogram['Reflection Lists']
1746	Phase = Phases[phase]
1747	hId = Histogram['hId']
1748	hfx = ':%d:'%(hId)
1749	wtFactor = calcControls[hfx+'wtFactor']
1750	pfx = '%d::'%(Phase['pId'])
1751	phfx = '%d:%d:'%(Phase['pId'],hId)
1752	SGData = Phase['General']['SGData']
1753	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1754	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1755	refDict = Histogram['Data']
1756	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1757	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1758	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1759	dependentVars = G2mv.GetDependentVars()
1760	depDerivDict = {}
1761	for j in dependentVars:
1762	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
1763	if calcControls['F**2']:
1764	for iref,ref in enumerate(refDict['RefList']):
1765	if ref[6] > 0:
1766	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
1767	w = 1.0/ref[6]
1768	if w*ref[5] >= calcControls['minF/sig']:
1769	for j,var in enumerate(varylist):
1770	if var in dFdvDict:
1771	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*dervCor
1772	for var in dependentVars:
1773	if var in dFdvDict:
1774	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*dervCor
1775	if phfx+'Scale' in varylist:
1776	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
1777	elif phfx+'Scale' in dependentVars:
1778	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
1779	for item in ['Ep','Es','Eg']:
1780	if phfx+item in varylist:
1781	dMdvh[varylist.index(phfx+item)][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
1782	elif phfx+item in dependentVars:
1783	depDerivDict[phfx+item][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
1784	for item in ['BabA','BabU']:
1785	if phfx+item in varylist:
1786	dMdvh[varylist.index(phfx+item)][iref] = wdervCordFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1787	elif phfx+item in dependentVars:
1788	depDerivDict[phfx+item][iref] = wdervCordFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1789	else:
1790	for iref,ref in enumerate(refDict['RefList']):
1791	if ref[5] > 0.:
1792	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
1793	Fo = np.sqrt(ref[5])
1794	Fc = np.sqrt(ref[7])
1795	w = 1.0/ref[6]
1796	if 2.0Fow*Fo >= calcControls['minF/sig']:
1797	for j,var in enumerate(varylist):
1798	if var in dFdvDict:
1799	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1800	for var in dependentVars:
1801	if var in dFdvDict:
1802	depDerivDict[var][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1803	if phfx+'Scale' in varylist:
1804	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
1805	elif phfx+'Scale' in dependentVars:
1806	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
1807	for item in ['Ep','Es','Eg']:
1808	if phfx+item in varylist:
1809	dMdvh[varylist.index(phfx+item)][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
1810	elif phfx+item in dependentVars:
1811	depDerivDict[phfx+item][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
1812	for item in ['BabA','BabU']:
1813	if phfx+item in varylist:
1814	dMdvh[varylist.index(phfx+item)][iref] = wdervCordFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1815	elif phfx+item in dependentVars:
1816	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1817	# now process derivatives in constraints
1818	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1819	else:
1820	continue #skip non-histogram entries
1821	if len(dMdv):
1822	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
1823	else:
1824	dMdv = np.sqrt(wtFactor)*dMdvh
1825
1826	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
1827	if np.any(pVals):
1828	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1829	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
1830
1831	return dMdv
1832
1833	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1834	'''Loop over histograms and compute derivatives of the fitting
1835	model (M) with respect to all parameters. For each histogram, the
1836	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
1837	number of parameters and m is the number of data points *in the
1838	histogram*. The (n by n) Hessian is computed from each Jacobian
1839	and it is returned. This routine is used when refinement
1840	derivatives are selected as "analtytic Hessian" in Controls.
1841
1842	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
1843	'''
1844	parmDict.update(zip(varylist,values))
1845	G2mv.Dict2Map(parmDict,varylist)
1846	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1847	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
1848	nvar = len(varylist)
1849	Hess = np.empty(0)
1850	histoList = Histograms.keys()
1851	histoList.sort()
1852	for histogram in histoList:
1853	if 'PWDR' in histogram[:4]:
1854	Histogram = Histograms[histogram]
1855	hId = Histogram['hId']
1856	hfx = ':%d:'%(hId)
1857	wtFactor = calcControls[hfx+'wtFactor']
1858	Limits = calcControls[hfx+'Limits']
1859	x,y,w,yc,yb,yd = Histogram['Data']
1860	W = wtFactor*w
1861	dy = y-yc
1862	xB = np.searchsorted(x,Limits[0])
1863	xF = np.searchsorted(x,Limits[1])
1864	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
1865	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1866	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
1867	Dy = dy[xB:xF][np.newaxis,:]
1868	dMdvh *= Wt
1869	if dlg:
1870	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1871	if len(Hess):
1872	Hess += np.inner(dMdvh,dMdvh)
1873	dMdvh = WtDy
1874	Vec += np.sum(dMdvh,axis=1)
1875	else:
1876	Hess = np.inner(dMdvh,dMdvh)
1877	dMdvh = WtDy
1878	Vec = np.sum(dMdvh,axis=1)
1879	elif 'HKLF' in histogram[:4]:
1880	Histogram = Histograms[histogram]
1881	nobs = Histogram['Residuals']['Nobs']
1882	phase = Histogram['Reflection Lists']
1883	Phase = Phases[phase]
1884	hId = Histogram['hId']
1885	hfx = ':%d:'%(hId)
1886	wtFactor = calcControls[hfx+'wtFactor']
1887	pfx = '%d::'%(Phase['pId'])
1888	phfx = '%d:%d:'%(Phase['pId'],hId)
1889	SGData = Phase['General']['SGData']
1890	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1891	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1892	refDict = Histogram['Data']
1893	time0 = time.time()
1894	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1895	print 'sf-deriv time: %.3f'%(time.time()-time0)
1896	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1897	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1898	dependentVars = G2mv.GetDependentVars()
1899	depDerivDict = {}
1900	for j in dependentVars:
1901	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
1902	wdf = np.zeros(len(refDict['RefList']))
1903	time0 = time.time()
1904	if calcControls['F**2']:
1905	for iref,ref in enumerate(refDict['RefList']):
1906	if ref[6] > 0:
1907	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
1908	w = 1.0/ref[6]
1909	if w*ref[5] >= calcControls['minF/sig']:
1910	wdf[iref] = w*(ref[5]-ref[7])
1911	for j,var in enumerate(varylist):
1912	if var in dFdvDict:
1913	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1914	for var in dependentVars:
1915	if var in dFdvDict:
1916	depDerivDict[var][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1917	if phfx+'Scale' in varylist:
1918	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
1919	elif phfx+'Scale' in dependentVars:
1920	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
1921	for item in ['Ep','Es','Eg']:
1922	if phfx+item in varylist:
1923	dMdvh[varylist.index(phfx+item)][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
1924	elif phfx+item in dependentVars:
1925	depDerivDict[phfx+item][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
1926	for item in ['BabA','BabU']:
1927	if phfx+item in varylist:
1928	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1929	elif phfx+item in dependentVars:
1930	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1931	else:
1932	for iref,ref in enumerate(refDict['RefList']):
1933	if ref[5] > 0.:
1934	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
1935	Fo = np.sqrt(ref[5])
1936	Fc = np.sqrt(ref[7])
1937	w = 1.0/ref[6]
1938	if 2.0Fow*Fo >= calcControls['minF/sig']:
1939	wdf[iref] = 2.0Fow*(Fo-Fc)
1940	for j,var in enumerate(varylist):
1941	if var in dFdvDict:
1942	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1943	for var in dependentVars:
1944	if var in dFdvDict:
1945	depDerivDict[var][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1946	if phfx+'Scale' in varylist:
1947	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
1948	elif phfx+'Scale' in dependentVars:
1949	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
1950	for item in ['Ep','Es','Eg']:
1951	if phfx+item in varylist:
1952	dMdvh[varylist.index(phfx+item)][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
1953	elif phfx+item in dependentVars:
1954	depDerivDict[phfx+item][iref] = wdervDict[phfx+item]parmDict[phfx+'Scale']
1955	for item in ['BabA','BabU']:
1956	if phfx+item in varylist:
1957	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1958	elif phfx+item in dependentVars:
1959	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1960	# now process derivatives in constraints
1961	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1962	print 'matrix build time: %.3f'%(time.time()-time0)
1963
1964	if dlg:
1965	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1966	if len(Hess):
1967	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
1968	Hess += wtFactor*np.inner(dMdvh,dMdvh)
1969	else:
1970	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
1971	Hess = wtFactor*np.inner(dMdvh,dMdvh)
1972	else:
1973	continue #skip non-histogram entries
1974	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
1975	if np.any(pVals):
1976	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1977	Vec += np.sum(dpdvpWtpVals,axis=1)
1978	Hess += np.inner(dpdv*pWt,dpdv)
1979	return Vec,Hess
1980
1981	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1982	'Needs a doc string'
1983	Values2Dict(parmDict, varylist, values)
1984	G2mv.Dict2Map(parmDict,varylist)
1985	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1986	M = np.empty(0)
1987	SumwYo = 0
1988	Nobs = 0
1989	ApplyRBModels(parmDict,Phases,rigidbodyDict)
1990	histoList = Histograms.keys()
1991	histoList.sort()
1992	for histogram in histoList:
1993	if 'PWDR' in histogram[:4]:
1994	Histogram = Histograms[histogram]
1995	hId = Histogram['hId']
1996	hfx = ':%d:'%(hId)
1997	wtFactor = calcControls[hfx+'wtFactor']
1998	Limits = calcControls[hfx+'Limits']
1999	x,y,w,yc,yb,yd = Histogram['Data']
2000	yc *= 0.0 #zero full calcd profiles
2001	yb *= 0.0
2002	yd *= 0.0
2003	xB = np.searchsorted(x,Limits[0])
2004	xF = np.searchsorted(x,Limits[1])
2005	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2006	varylist,Histogram,Phases,calcControls,pawleyLookup)
2007	yc[xB:xF] += yb[xB:xF]
2008	if not np.any(y): #fill dummy data
2009	rv = st.poisson(yc[xB:xF])
2010	y[xB:xF] = rv.rvs()
2011	Z = np.ones_like(yc[xB:xF])
2012	Z[1::2] *= -1
2013	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2014	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2015	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2016	W = wtFactor*w
2017	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2018	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2019	Nobs += Histogram['Residuals']['Nobs']
2020	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
2021	SumwYo += Histogram['Residuals']['sumwYo']
2022	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2023	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
2024	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2025	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
2026	sumYB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yc[xB:xF])*ma.abs(y[xB:xF]-yb[xB:xF])/y[xB:xF],0.))
2027	sumwYB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](ma.abs(y[xB:xF]-yc[xB:xF])ma.abs(y[xB:xF]-yb[xB:xF])/y[xB:xF])**2,0.))
2028	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2029	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2030	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2031	if dlg:
2032	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2033	M = np.concatenate((M,wdy))
2034	#end of PWDR processing
2035	elif 'HKLF' in histogram[:4]:
2036	Histogram = Histograms[histogram]
2037	Histogram['Residuals'] = {}
2038	phase = Histogram['Reflection Lists']
2039	Phase = Phases[phase]
2040	hId = Histogram['hId']
2041	hfx = ':%d:'%(hId)
2042	wtFactor = calcControls[hfx+'wtFactor']
2043	pfx = '%d::'%(Phase['pId'])
2044	phfx = '%d:%d:'%(Phase['pId'],hId)
2045	SGData = Phase['General']['SGData']
2046	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2047	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2048	refDict = Histogram['Data']
2049	time0 = time.time()
2050	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2051	print 'sf-calc time: %.3f'%(time.time()-time0)
2052	df = np.zeros(len(refDict['RefList']))
2053	sumwYo = 0
2054	sumFo = 0
2055	sumFo2 = 0
2056	sumdF = 0
2057	sumdF2 = 0
2058	nobs = 0
2059	if calcControls['F**2']:
2060	for i,ref in enumerate(refDict['RefList']):
2061	if ref[6] > 0:
2062	SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
2063	w = 1.0/ref[6]
2064	ref[7] = parmDict[phfx+'Scale']*ref[9]
2065	ref[7] *= ref[11] #correct Fc^2 for extinction
2066	ref[8] = ref[5]/parmDict[phfx+'Scale']
2067	if w*ref[5] >= calcControls['minF/sig']:
2068	sumFo2 += ref[5]
2069	Fo = np.sqrt(ref[5])
2070	sumFo += Fo
2071	sumFo2 += ref[5]
2072	sumdF += abs(Fo-np.sqrt(ref[7]))
2073	sumdF2 += abs(ref[5]-ref[7])
2074	nobs += 1
2075	df[i] = -w*(ref[5]-ref[7])
2076	sumwYo += (wref[5])*2
2077	else:
2078	for i,ref in enumerate(refDict['RefList']):
2079	if ref[5] > 0.:
2080	SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
2081	ref[7] = parmDict[phfx+'Scale']*ref[9]
2082	ref[7] *= ref[11] #correct Fc^2 for extinction
2083	Fo = np.sqrt(ref[5])
2084	Fc = np.sqrt(ref[7])
2085	w = 2.0*Fo/ref[6]
2086	if w*Fo >= calcControls['minF/sig']:
2087	sumFo += Fo
2088	sumFo2 += ref[5]
2089	sumdF += abs(Fo-Fc)
2090	sumdF2 += abs(ref[5]-ref[7])
2091	nobs += 1
2092	df[i] = -w*(Fo-Fc)
2093	sumwYo += (wFo)*2
2094	Histogram['Residuals']['Nobs'] = nobs
2095	Histogram['Residuals']['sumwYo'] = sumwYo
2096	SumwYo += sumwYo
2097	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/Histogram['Residuals']['sumwYo'])100.)
2098	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2099	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2100	Histogram['Residuals'][phfx+'Nref'] = nobs
2101	Nobs += nobs
2102	if dlg:
2103	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2104	M = np.concatenate((M,wtFactor*df))
2105	# end of HKLF processing
2106	Histograms['sumwYo'] = SumwYo
2107	Histograms['Nobs'] = Nobs
2108	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
2109	if dlg:
2110	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2111	if not GoOn:
2112	parmDict['saved values'] = values
2113	dlg.Destroy()
2114	raise Exception #Abort!!
2115	pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2116	if len(pVals):
2117	pSum = np.sum(pWtpVals*2)
2118	for name in pWsum:
2119	if pWsum:
2120	print ' Penalty function for %8s = %12.5g'%(name,pWsum[name])
2121	print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2122	Nobs += len(pVals)
2123	M = np.concatenate((M,np.sqrt(pWt)*pVals))
2124	return M
2125

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

Download in other formats: