Context Navigation

source: branch/logging/GSASIIstrMath.py @ 2927

Last change on this file since 2927 was 1456, checked in by vondreele, 11 years ago
consolidate HKLF parts of dervRefine & HessRefine? into new routine dervHKLF
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 101.5 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMath - structure math routines
4	-----------------------------------------
5	'''
6	########### SVN repository information ###################
7	# $Date: 2014-08-04 19:29:15 +0000 (Mon, 04 Aug 2014) $
8	# $Author: vondreele $
9	# $Revision: 1456 $
10	# $URL: branch/logging/GSASIIstrMath.py $
11	# $Id: GSASIIstrMath.py 1456 2014-08-04 19:29:15Z 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: 1456 $")
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 'NC' 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	if 'NT' in calcControls[hfx+'histType']:
575	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[12])
576	fbs = np.array([0,0])
577	H = refl[:3]
578	SQ = 1./(2.refl[4])*2
579	SQfactor = 4.0SQtwopisq
580	Bab = parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor)
581	if not np.any(refDict['FF']['FF'][iref]): #no form factors - 1st time thru StructureFactor
582	if 'N' in calcControls[hfx+'histType']:
583	dat = G2el.getBLvalues(BLtables)
584	refDict['FF']['FF'][iref] = dat.values()
585	else: #'X'
586	dat = G2el.getFFvalues(FFtables,SQ)
587	refDict['FF']['FF'][iref] = dat.values()
588	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
589	FF = refDict['FF']['FF'][iref][Tindx]
590	Uniq = np.inner(H,SGMT)
591	Phi = np.inner(H,SGT)
592	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
593	sinp = np.sin(phase)
594	cosp = np.cos(phase)
595	biso = -SQfactor*Uisodata
596	Tiso = np.where(biso<1.,np.exp(biso),1.0)
597	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
598	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
599	Tcorr = TisoTuijMdata*Fdata/len(Uniq)
600	fa = np.array([(FF+FP-Bab)cospTcorr,-FPPsinpTcorr])
601	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
602	if not SGData['SGInv']:
603	fb = np.array([(FF+FP-Bab)sinpTcorr,FPPcospTcorr])
604	fbs = np.sum(np.sum(fb,axis=1),axis=1)
605	fasq = fas**2
606	fbsq = fbs**2 #imaginary
607	refl[9] = np.sum(fasq)+np.sum(fbsq)
608	refl[10] = atan2d(fbs[0],fas[0])
609
610	def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
611	''' Compute structure factors for all h,k,l for phase
612	puts the result, F^2, in each ref[8] in refList
613	input:
614
615	:param dict refDict: where
616	'RefList' list where each ref = h,k,l,m,d,...
617	'FF' dict of form factors - filed in below
618	:param np.array G: reciprocal metric tensor
619	:param str pfx: phase id string
620	:param dict SGData: space group info. dictionary output from SpcGroup
621	:param dict calcControls:
622	:param dict ParmDict:
623
624	'''
625	twopi = 2.0*np.pi
626	twopisq = 2.0np.pi*2
627	phfx = pfx.split(':')[0]+hfx
628	ast = np.sqrt(np.diag(G))
629	Mast = twopisq*np.multiply.outer(ast,ast)
630	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
631	SGT = np.array([ops[1] for ops in SGData['SGOps']])
632	FFtables = calcControls['FFtables']
633	BLtables = calcControls['BLtables']
634	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
635	FF = np.zeros(len(Tdata))
636	if 'NC' in calcControls[hfx+'histType']:
637	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
638	elif 'X' in calcControls[hfx+'histType']:
639	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
640	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
641	Uij = np.array(G2lat.U6toUij(Uijdata))
642	bij = Mast*Uij.T
643	blkSize = 100 #no. of reflections in a block
644	nRef = refDict['RefList'].shape[0]
645	if not len(refDict['FF']): #no form factors - 1st time thru StructureFactor
646	if 'N' in calcControls[hfx+'histType']:
647	dat = G2el.getBLvalues(BLtables)
648	refDict['FF']['El'] = dat.keys()
649	refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
650	else: #'X'
651	dat = G2el.getFFvalues(FFtables,0.)
652	refDict['FF']['El'] = dat.keys()
653	refDict['FF']['FF'] = np.ones((nRef,len(dat)))
654	for iref,ref in enumerate(refDict['RefList']):
655	SQ = 1./(2.ref[4])*2
656	dat = G2el.getFFvalues(FFtables,SQ)
657	refDict['FF']['FF'][iref] *= dat.values()
658	#reflection processing begins here - big arrays!
659	iBeg = 0
660	while iBeg < nRef:
661	iFin = min(iBeg+blkSize,nRef)
662	refl = refDict['RefList'][iBeg:iFin]
663	H = refl.T[:3]
664	SQ = 1./(2.refl.T[4])*2
665	SQfactor = 4.0SQtwopisq
666	if 'T' in calcControls[hfx+'histType']:
667	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
668	FP = np.repeat(FP.T,len(SGT),axis=0)
669	FPP = np.repeat(FPP.T,len(SGT),axis=0)
670	Bab = np.repeat(parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor),len(SGT))
671	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
672	FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
673	Uniq = np.reshape(np.inner(H.T,SGMT),(-1,3))
674	Phi = np.inner(H.T,SGT).flatten()
675	phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T)+Phi[:,np.newaxis])
676	sinp = np.sin(phase)
677	cosp = np.cos(phase)
678	biso = -SQfactor*Uisodata[:,np.newaxis]
679	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
680	HbH = -np.sum(Uniq.T*np.inner(bij,Uniq),axis=1)
681	Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
682	Tcorr = TisoTuijMdata*Fdata/len(SGMT)
683	fa = np.array([((FF+FP).T-Bab).TcospTcorr,-FPPsinpTcorr])
684	fa = np.reshape(fa,(2,len(refl),len(SGT),len(Mdata)))
685	fas = np.sum(np.sum(fa,axis=2),axis=2) #real
686	fbs = np.zeros_like(fas)
687	if not SGData['SGInv']:
688	fb = np.array([((FF+FP).T-Bab).TsinpTcorr,FPPcospTcorr])
689	fb = np.reshape(fb,(2,len(refl),len(SGT),len(Mdata)))
690	fbs = np.sum(np.sum(fb,axis=2),axis=2)
691	fasq = fas**2
692	fbsq = fbs**2 #imaginary
693	refl.T[9] = np.sum(fasq,axis=0)+np.sum(fbsq,axis=0)
694	refl.T[10] = atan2d(fbs[0],fas[0])
695	iBeg += blkSize
696
697	def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
698	'Needs a doc string'
699	twopi = 2.0*np.pi
700	twopisq = 2.0np.pi*2
701	phfx = pfx.split(':')[0]+hfx
702	ast = np.sqrt(np.diag(G))
703	Mast = twopisq*np.multiply.outer(ast,ast)
704	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
705	SGT = np.array([ops[1] for ops in SGData['SGOps']])
706	FFtables = calcControls['FFtables']
707	BLtables = calcControls['BLtables']
708	nRef = len(refDict['RefList'])
709	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
710	mSize = len(Mdata)
711	FF = np.zeros(len(Tdata))
712	if 'N' in calcControls[hfx+'histType']:
713	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
714	else:
715	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
716	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
717	Uij = np.array(G2lat.U6toUij(Uijdata))
718	bij = Mast*Uij.T
719	dFdvDict = {}
720	dFdfr = np.zeros((nRef,mSize))
721	dFdx = np.zeros((nRef,mSize,3))
722	dFdui = np.zeros((nRef,mSize))
723	dFdua = np.zeros((nRef,mSize,6))
724	dFdbab = np.zeros((nRef,2))
725	for iref,refl in enumerate(refDict['RefList']):
726	H = np.array(refl[:3])
727	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
728	SQfactor = 8.0SQnp.pi**2
729	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
730	Bab = parmDict[phfx+'BabA']*dBabdA
731	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
732	FF = refDict['FF']['FF'][iref].T[Tindx]
733	# FF = [refDict['FF'][iref][El] for El in Tdata]
734	Uniq = np.inner(H,SGMT)
735	Phi = np.inner(H,SGT)
736	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
737	sinp = np.sin(phase)
738	cosp = np.cos(phase)
739	occ = Mdata*Fdata/len(Uniq)
740	biso = -SQfactor*Uisodata
741	Tiso = np.where(biso<1.,np.exp(biso),1.0)
742	HbH = -np.inner(H,np.inner(bij,H))
743	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
744	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
745	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
746	Tcorr = Tiso*Tuij
747	fot = (FF+FP-Bab)occTcorr
748	fotp = FPPoccTcorr
749	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
750	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
751
752	fas = np.sum(np.sum(fa,axis=1),axis=1)
753	fbs = np.sum(np.sum(fb,axis=1),axis=1)
754	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
755	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
756	#sum below is over Uniq
757	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2) #Fdata != 0 ever avoids /0. problem
758	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
759	dfadui = np.sum(-SQfactor*fa,axis=2)
760	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
761	dfadba = np.sum(-cosp(occTcorr)[:,np.newaxis],axis=1)
762	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for
763	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
764	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
765	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
766	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
767	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
768	if not SGData['SGInv']:
769	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #Fdata != 0 ever avoids /0. problem
770	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
771	dfbdui = np.sum(-SQfactor*fb,axis=2)
772	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
773	dfbdba = np.sum(-sinp(occTcorr)[:,np.newaxis],axis=1)
774	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
775	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
776	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
777	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
778	dFdbab[iref] += 2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
779	#loop over atoms - each dict entry is list of derivatives for all the reflections
780	for i in range(len(Mdata)):
781	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
782	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
783	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
784	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
785	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
786	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
787	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
788	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
789	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
790	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
791	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
792	dFdvDict[pfx+'BabA'] = dFdbab.T[0]
793	dFdvDict[pfx+'BabU'] = dFdbab.T[1]
794	return dFdvDict
795
796	def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
797	''' Single crystal extinction function; puts correction in ref[11] and returns
798	corrections needed for derivatives
799	'''
800	ref[11] = 1.0
801	dervCor = 1.0
802	dervDict = {}
803	if calcControls[phfx+'EType'] != 'None':
804	SQ = 1/(4.ref[4]*2)
805	if 'C' in parmDict[hfx+'Type']:
806	cos2T = 1.0-2.SQparmDict[hfx+'Lam']**2 #cos(2theta)
807	else: #'T'
808	cos2T = 1.0-2.SQref[12]**2 #cos(2theta)
809	if 'SXC' in parmDict[hfx+'Type']:
810	AV = 7.9406e5/parmDict[pfx+'Vol']**2
811	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
812	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
813	PLZ = AVP12ref[7]parmDict[hfx+'Lam']*2
814	elif 'SNT' in parmDict[hfx+'Type']:
815	AV = 1.e7/parmDict[pfx+'Vol']**2
816	PL = SQ
817	PLZ = AVref[7]ref[12]**2
818	elif 'SNC' in parmDict[hfx+'Type']:
819	AV = 1.e7/parmDict[pfx+'Vol']**2
820	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
821	PLZ = AVref[7]parmDict[hfx+'Lam']**2
822
823	if 'Primary' in calcControls[phfx+'EType']:
824	PLZ *= 1.5
825	else:
826	if 'C' in parmDict[hfx+'Type']:
827	PLZ *= calcControls[phfx+'Tbar']
828	else: #'T'
829	PLZ *= ref[13]
830	if 'Primary' in calcControls[phfx+'EType']:
831	PLZ *= 1.5
832	PSIG = parmDict[phfx+'Ep']
833	elif 'I & II' in calcControls[phfx+'EType']:
834	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
835	elif 'Type II' in calcControls[phfx+'EType']:
836	PSIG = parmDict[phfx+'Es']
837	else: # 'Secondary Type I'
838	PSIG = parmDict[phfx+'Eg']/PL
839
840	AG = 0.58+0.48cos2T+0.24cos2T**2
841	AL = 0.025+0.285*cos2T
842	BG = 0.02-0.025*cos2T
843	BL = 0.15-0.2(0.75-cos2T)*2
844	if cos2T < 0.:
845	BL = -0.45*cos2T
846	CG = 2.
847	CL = 2.
848	PF = PLZ*PSIG
849
850	if 'Gaussian' in calcControls[phfx+'EApprox']:
851	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
852	extCor = np.sqrt(PF4)
853	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
854	else:
855	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
856	extCor = np.sqrt(PF4)
857	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
858
859	dervCor = (1.+PF)*PF3 #extinction corr for other derivatives
860	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
861	dervDict[phfx+'Ep'] = -ref[7]PLZPF3
862	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
863	dervDict[phfx+'Es'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
864	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
865	dervDict[phfx+'Eg'] = -ref[7]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
866
867	ref[11] = 1./extCor
868	return dervCor,dervDict
869
870	def Dict2Values(parmdict, varylist):
871	'''Use before call to leastsq to setup list of values for the parameters
872	in parmdict, as selected by key in varylist'''
873	return [parmdict[key] for key in varylist]
874
875	def Values2Dict(parmdict, varylist, values):
876	''' Use after call to leastsq to update the parameter dictionary with
877	values corresponding to keys in varylist'''
878	parmdict.update(zip(varylist,values))
879
880	def GetNewCellParms(parmDict,varyList):
881	'Needs a doc string'
882	newCellDict = {}
883	Anames = ['A'+str(i) for i in range(6)]
884	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
885	for item in varyList:
886	keys = item.split(':')
887	if keys[2] in Ddict:
888	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
889	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
890	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
891	return newCellDict # is e.g. {'0::D11':A0+D11}
892
893	def ApplyXYZshifts(parmDict,varyList):
894	'''
895	takes atom x,y,z shift and applies it to corresponding atom x,y,z value
896
897	:param dict parmDict: parameter dictionary
898	:param list varyList: list of variables (not used!)
899	:returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
900
901	'''
902	newAtomDict = {}
903	for item in parmDict:
904	if 'dA' in item:
905	parm = ''.join(item.split('d'))
906	parmDict[parm] += parmDict[item]
907	newAtomDict[item] = [parm,parmDict[parm]]
908	return newAtomDict
909
910	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
911	'Spherical harmonics texture'
912	IFCoup = 'Bragg' in calcControls[hfx+'instType']
913	odfCor = 1.0
914	H = refl[:3]
915	cell = G2lat.Gmat2cell(g)
916	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
917	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
918	phi,beta = G2lat.CrsAng(H,cell,SGData)
919	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
920	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
921	for item in SHnames:
922	L,M,N = eval(item.strip('C'))
923	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
924	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
925	Lnorm = G2lat.Lnorm(L)
926	odfCor += parmDict[pfx+item]LnormKcl*Ksl
927	return odfCor
928
929	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
930	'Spherical harmonics texture derivatives'
931	FORPI = 4.0*np.pi
932	IFCoup = 'Bragg' in calcControls[hfx+'instType']
933	odfCor = 1.0
934	dFdODF = {}
935	dFdSA = [0,0,0]
936	H = refl[:3]
937	cell = G2lat.Gmat2cell(g)
938	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
939	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
940	phi,beta = G2lat.CrsAng(H,cell,SGData)
941	psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
942	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
943	for item in SHnames:
944	L,M,N = eval(item.strip('C'))
945	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
946	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
947	Lnorm = G2lat.Lnorm(L)
948	odfCor += parmDict[pfx+item]LnormKcl*Ksl
949	dFdODF[pfx+item] = LnormKclKsl
950	for i in range(3):
951	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
952	return odfCor,dFdODF,dFdSA
953
954	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
955	'spherical harmonics preferred orientation (cylindrical symmetry only)'
956	odfCor = 1.0
957	H = refl[:3]
958	cell = G2lat.Gmat2cell(g)
959	Sangl = [0.,0.,0.]
960	if 'Bragg' in calcControls[hfx+'instType']:
961	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
962	IFCoup = True
963	else:
964	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
965	IFCoup = False
966	phi,beta = G2lat.CrsAng(H,cell,SGData)
967	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
968	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
969	for item in SHnames:
970	L,N = eval(item.strip('C'))
971	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
972	odfCor += parmDict[phfx+item]LnormKcsl
973	return np.squeeze(odfCor)
974
975	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
976	'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
977	FORPI = 12.5663706143592
978	odfCor = 1.0
979	dFdODF = {}
980	H = refl[:3]
981	cell = G2lat.Gmat2cell(g)
982	Sangl = [0.,0.,0.]
983	if 'Bragg' in calcControls[hfx+'instType']:
984	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
985	IFCoup = True
986	else:
987	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
988	IFCoup = False
989	phi,beta = G2lat.CrsAng(H,cell,SGData)
990	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
991	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
992	for item in SHnames:
993	L,N = eval(item.strip('C'))
994	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
995	odfCor += parmDict[phfx+item]LnormKcsl
996	dFdODF[phfx+item] = Kcsl*Lnorm
997	return odfCor,dFdODF
998
999	def GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1000	'Needs a doc string'
1001	POcorr = 1.0
1002	if calcControls[phfx+'poType'] == 'MD':
1003	MD = parmDict[phfx+'MD']
1004	if MD != 1.0:
1005	MDAxis = calcControls[phfx+'MDAxis']
1006	sumMD = 0
1007	for H in uniq:
1008	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1009	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1010	sumMD += A**3
1011	POcorr = sumMD/len(uniq)
1012	else: #spherical harmonics
1013	if calcControls[phfx+'SHord']:
1014	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1015	return POcorr
1016
1017	def GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1018	'Needs a doc string'
1019	POcorr = 1.0
1020	POderv = {}
1021	if calcControls[phfx+'poType'] == 'MD':
1022	MD = parmDict[phfx+'MD']
1023	MDAxis = calcControls[phfx+'MDAxis']
1024	sumMD = 0
1025	sumdMD = 0
1026	for H in uniq:
1027	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1028	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1029	sumMD += A**3
1030	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
1031	POcorr = sumMD/len(uniq)
1032	POderv[phfx+'MD'] = sumdMD/len(uniq)
1033	else: #spherical harmonics
1034	if calcControls[phfx+'SHord']:
1035	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1036	return POcorr,POderv
1037
1038	def GetAbsorb(refl,hfx,calcControls,parmDict):
1039	'Needs a doc string'
1040	if 'Debye' in calcControls[hfx+'instType']:
1041	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1042	else:
1043	return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1044
1045	def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
1046	'Needs a doc string'
1047	if 'Debye' in calcControls[hfx+'instType']:
1048	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1049	else:
1050	return G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1051
1052	def GetIntensityCorr(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1053	'Needs a doc string' #need powder extinction!
1054	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
1055	if 'X' in parmDict[hfx+'Type']:
1056	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1057	Icorr *= GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1058	if pfx+'SHorder' in parmDict:
1059	Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1060	Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
1061	refl[11] = Icorr
1062
1063	def GetIntensityDerv(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1064	'Needs a doc string' #need powder extinction derivs!
1065	dIdsh = 1./parmDict[hfx+'Scale']
1066	dIdsp = 1./parmDict[phfx+'Scale']
1067	if 'X' in parmDict[hfx+'Type']:
1068	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1069	dIdPola /= pola
1070	else: #'N'
1071	dIdPola = 0.0
1072	POcorr,dIdPO = GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1073	for iPO in dIdPO:
1074	dIdPO[iPO] /= POcorr
1075	dFdODF = {}
1076	dFdSA = [0,0,0]
1077	if pfx+'SHorder' in parmDict:
1078	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1079	for iSH in dFdODF:
1080	dFdODF[iSH] /= odfCor
1081	for i in range(3):
1082	dFdSA[i] /= odfCor
1083	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
1084	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
1085
1086	def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
1087	'Needs a doc string'
1088	costh = cosd(refl[5]/2.)
1089	#crystallite size
1090	if calcControls[phfx+'SizeType'] == 'isotropic':
1091	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1092	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1093	H = np.array(refl[:3])
1094	P = np.array(calcControls[phfx+'SizeAxis'])
1095	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1096	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
1097	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1098	else: #ellipsoidal crystallites
1099	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1100	H = np.array(refl[:3])
1101	lenR = G2pwd.ellipseSize(H,Sij,GB)
1102	Sgam = 1.8wave/(np.picosth*lenR)
1103	#microstrain
1104	if calcControls[phfx+'MustrainType'] == 'isotropic':
1105	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
1106	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1107	H = np.array(refl[:3])
1108	P = np.array(calcControls[phfx+'MustrainAxis'])
1109	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1110	Si = parmDict[phfx+'Mustrain;i']
1111	Sa = parmDict[phfx+'Mustrain;a']
1112	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
1113	else: #generalized - P.W. Stephens model
1114	pwrs = calcControls[phfx+'MuPwrs']
1115	sum = 0
1116	for i,pwr in enumerate(pwrs):
1117	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1118	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
1119	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
1120	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
1121	sig /= ateln2
1122	return sig,gam
1123
1124	def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
1125	'Needs a doc string'
1126	gamDict = {}
1127	sigDict = {}
1128	costh = cosd(refl[5]/2.)
1129	tanth = tand(refl[5]/2.)
1130	#crystallite size derivatives
1131	if calcControls[phfx+'SizeType'] == 'isotropic':
1132	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1133	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
1134	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
1135	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1136	H = np.array(refl[:3])
1137	P = np.array(calcControls[phfx+'SizeAxis'])
1138	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1139	Si = parmDict[phfx+'Size;i']
1140	Sa = parmDict[phfx+'Size;a']
1141	gami = (1.8wave/np.pi)/(SiSa)
1142	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1143	Sgam = gami*sqtrm
1144	gam = Sgam/costh
1145	dsi = (gamiSicosP*2/(sqtrmcosth)-gam/Si)
1146	dsa = (gamiSasinP*2/(sqtrmcosth)-gam/Sa)
1147	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1148	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1149	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1150	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1151	else: #ellipsoidal crystallites
1152	const = 1.8wave/(np.picosth)
1153	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1154	H = np.array(refl[:3])
1155	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1156	Sgam = 1.8wave/(np.picosth*lenR)
1157	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1158	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1159	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1160	gamDict[phfx+'Size;mx'] = Sgam
1161	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
1162
1163	#microstrain derivatives
1164	if calcControls[phfx+'MustrainType'] == 'isotropic':
1165	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5]/2.)/np.pi
1166	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
1167	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
1168	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1169	H = np.array(refl[:3])
1170	P = np.array(calcControls[phfx+'MustrainAxis'])
1171	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1172	Si = parmDict[phfx+'Mustrain;i']
1173	Sa = parmDict[phfx+'Mustrain;a']
1174	gami = 0.018SiSa*tanth/np.pi
1175	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1176	Mgam = gami/sqtrm
1177	dsi = -gamiSicosP2/sqtrm3
1178	dsa = -gamiSasinP2/sqtrm3
1179	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1180	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1181	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1182	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1183	else: #generalized - P.W. Stephens model
1184	pwrs = calcControls[phfx+'MuPwrs']
1185	const = 0.018refl[4]2tanth
1186	sum = 0
1187	for i,pwr in enumerate(pwrs):
1188	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
1189	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1190	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain;mx']
1191	sigDict[phfx+'Mustrain:'+str(i)] = \
1192	2.constterm(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1193	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
1194	for i in range(len(pwrs)):
1195	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1196	gamDict[phfx+'Mustrain;mx'] = Mgam
1197	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1198	return sigDict,gamDict
1199
1200	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1201	'Needs a doc string'
1202	h,k,l = refl[:3]
1203	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1204	d = np.sqrt(dsq)
1205
1206	refl[4] = d
1207	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
1208	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1209	if 'Bragg' in calcControls[hfx+'instType']:
1210	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1211	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
1212	else: #Debye-Scherrer - simple but maybe not right
1213	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1214	return pos
1215
1216	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1217	'Needs a doc string'
1218	dpr = 180./np.pi
1219	h,k,l = refl[:3]
1220	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1221	dst = np.sqrt(dstsq)
1222	pos = refl[5]-parmDict[hfx+'Zero']
1223	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
1224	dpdw = const*dst
1225	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
1226	dpdA = constwave/(2.0*dst)
1227	dpdZ = 1.0
1228	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1229	if 'Bragg' in calcControls[hfx+'instType']:
1230	dpdSh = -4.constcosd(pos/2.0)
1231	dpdTr = -constsind(pos)100.0
1232	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1233	else: #Debye-Scherrer - simple but maybe not right
1234	dpdXd = -const*cosd(pos)
1235	dpdYd = -const*sind(pos)
1236	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1237
1238	def GetHStrainShift(refl,SGData,phfx,parmDict):
1239	'Needs a doc string'
1240	laue = SGData['SGLaue']
1241	uniq = SGData['SGUniq']
1242	h,k,l = refl[:3]
1243	if laue in ['m3','m3m']:
1244	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
1245	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
1246	elif laue in ['6/m','6/mmm','3m1','31m','3']:
1247	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
1248	elif laue in ['3R','3mR']:
1249	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
1250	elif laue in ['4/m','4/mmm']:
1251	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
1252	elif laue in ['mmm']:
1253	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
1254	elif laue in ['2/m']:
1255	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
1256	if uniq == 'a':
1257	Dij += parmDict[phfx+'D23']kl
1258	elif uniq == 'b':
1259	Dij += parmDict[phfx+'D13']hl
1260	elif uniq == 'c':
1261	Dij += parmDict[phfx+'D12']hk
1262	else:
1263	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
1264	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
1265	return -Dijrefl[4]2tand(refl[5]/2.0)
1266
1267	def GetHStrainShiftDerv(refl,SGData,phfx):
1268	'Needs a doc string'
1269	laue = SGData['SGLaue']
1270	uniq = SGData['SGUniq']
1271	h,k,l = refl[:3]
1272	if laue in ['m3','m3m']:
1273	dDijDict = {phfx+'D11':h2+k2+l**2,
1274	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
1275	elif laue in ['6/m','6/mmm','3m1','31m','3']:
1276	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
1277	elif laue in ['3R','3mR']:
1278	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
1279	elif laue in ['4/m','4/mmm']:
1280	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
1281	elif laue in ['mmm']:
1282	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
1283	elif laue in ['2/m']:
1284	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
1285	if uniq == 'a':
1286	dDijDict[phfx+'D23'] = k*l
1287	elif uniq == 'b':
1288	dDijDict[phfx+'D13'] = h*l
1289	elif uniq == 'c':
1290	dDijDict[phfx+'D12'] = h*k
1291	names.append()
1292	else:
1293	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
1294	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
1295	for item in dDijDict:
1296	dDijDict[item] = -refl[4]2tand(refl[5]/2.0)
1297	return dDijDict
1298
1299	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1300	'Needs a doc string'
1301	histoList = Histograms.keys()
1302	histoList.sort()
1303	for histogram in histoList:
1304	if 'PWDR' in histogram[:4]:
1305	Histogram = Histograms[histogram]
1306	hId = Histogram['hId']
1307	hfx = ':%d:'%(hId)
1308	Limits = calcControls[hfx+'Limits']
1309	shl = max(parmDict[hfx+'SH/L'],0.0005)
1310	Ka2 = False
1311	kRatio = 0.0
1312	if hfx+'Lam1' in parmDict.keys():
1313	Ka2 = True
1314	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1315	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1316	x,y,w,yc,yb,yd = Histogram['Data']
1317	xB = np.searchsorted(x,Limits[0])
1318	xF = np.searchsorted(x,Limits[1])
1319	ymb = np.array(y-yb)
1320	ymb = np.where(ymb,ymb,1.0)
1321	ycmb = np.array(yc-yb)
1322	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
1323	refLists = Histogram['Reflection Lists']
1324	for phase in refLists:
1325	Phase = Phases[phase]
1326	pId = Phase['pId']
1327	phfx = '%d:%d:'%(pId,hId)
1328	refDict = refLists[phase]
1329	sumFo = 0.0
1330	sumdF = 0.0
1331	sumFosq = 0.0
1332	sumdFsq = 0.0
1333	for refl in refDict['RefList']:
1334	if 'C' in calcControls[hfx+'histType']:
1335	yp = np.zeros_like(yb)
1336	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1337	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1338	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1339	iFin2 = iFin
1340	yp[iBeg:iFin] = refl[11]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
1341	if Ka2:
1342	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
1343	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1344	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1345	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1346	if not iBeg2+iFin2: #peak below low limit - skip peak
1347	continue
1348	elif not iBeg2-iFin2: #peak above high limit - done
1349	break
1350	yp[iBeg2:iFin2] += refl[11]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here
1351	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[11](1.+kRatio)),0.0))
1352	elif 'T' in calcControls[hfx+'histType']:
1353	print 'TOF Undefined at present'
1354	raise Exception #no TOF yet
1355	Fo = np.sqrt(np.abs(refl[8]))
1356	Fc = np.sqrt(np.abs(refl[9]))
1357	sumFo += Fo
1358	sumFosq += refl[8]**2
1359	sumdF += np.abs(Fo-Fc)
1360	sumdFsq += (refl[8]-refl[9])**2
1361	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1362	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1363	Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1364	Histogram['Residuals']['hId'] = hId
1365	elif 'HKLF' in histogram[:4]:
1366	Histogram = Histograms[histogram]
1367	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1368
1369	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1370	'Needs a doc string'
1371
1372	def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1373	U = parmDict[hfx+'U']
1374	V = parmDict[hfx+'V']
1375	W = parmDict[hfx+'W']
1376	X = parmDict[hfx+'X']
1377	Y = parmDict[hfx+'Y']
1378	tanPos = tand(refl[5]/2.0)
1379	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
1380	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
1381	sig = max(0.001,sig)
1382	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
1383	gam = max(0.001,gam)
1384	return sig,gam
1385
1386	hId = Histogram['hId']
1387	hfx = ':%d:'%(hId)
1388	bakType = calcControls[hfx+'bakType']
1389	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1390	yc = np.zeros_like(yb)
1391
1392	if 'C' in calcControls[hfx+'histType']:
1393	shl = max(parmDict[hfx+'SH/L'],0.002)
1394	Ka2 = False
1395	if hfx+'Lam1' in parmDict.keys():
1396	wave = parmDict[hfx+'Lam1']
1397	Ka2 = True
1398	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1399	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1400	else:
1401	wave = parmDict[hfx+'Lam']
1402	else:
1403	print 'TOF Undefined at present - might be nothing need be done here'
1404	for phase in Histogram['Reflection Lists']:
1405	refDict = Histogram['Reflection Lists'][phase]
1406	Phase = Phases[phase]
1407	pId = Phase['pId']
1408	pfx = '%d::'%(pId)
1409	phfx = '%d:%d:'%(pId,hId)
1410	hfx = ':%d:'%(hId)
1411	SGData = Phase['General']['SGData']
1412	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1413	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] #Do I want to modify by Dij?
1414	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1415	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
1416	Vst = np.sqrt(nl.det(G)) #V*
1417	if not Phase['General'].get('doPawley'):
1418	time0 = time.time()
1419	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1420	# print 'sf calc time: %.3fs'%(time.time()-time0)
1421	time0 = time.time()
1422	for iref,refl in enumerate(refDict['RefList']):
1423	if 'C' in calcControls[hfx+'histType']:
1424	h,k,l = refl[:3]
1425	Uniq = np.inner(refl[:3],SGMT)
1426	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
1427	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
1428	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
1429	refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam
1430	GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[11]
1431	refl[11] = VstLorenz
1432	if Phase['General'].get('doPawley'):
1433	try:
1434	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1435	refl[9] = parmDict[pInd]
1436	except KeyError:
1437	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
1438	continue
1439	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1440	iBeg = np.searchsorted(x,refl[5]-fmin)
1441	iFin = np.searchsorted(x,refl[5]+fmax)
1442	if not iBeg+iFin: #peak below low limit - skip peak
1443	continue
1444	elif not iBeg-iFin: #peak above high limit - done
1445	break
1446	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
1447	if Ka2:
1448	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
1449	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1450	iBeg = np.searchsorted(x,pos2-fmin)
1451	iFin = np.searchsorted(x,pos2+fmax)
1452	if not iBeg+iFin: #peak below low limit - skip peak
1453	continue
1454	elif not iBeg-iFin: #peak above high limit - done
1455	return yc,yb
1456	yc[iBeg:iFin] += refl[11]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #and here
1457	elif 'T' in calcControls[hfx+'histType']:
1458	print 'TOF Undefined at present'
1459	raise Exception #no TOF yet
1460	# print 'profile calc time: %.3fs'%(time.time()-time0)
1461	return yc,yb
1462
1463	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1464	'Needs a doc string'
1465
1466	def cellVaryDerv(pfx,SGData,dpdA):
1467	if SGData['SGLaue'] in ['-1',]:
1468	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1469	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1470	elif SGData['SGLaue'] in ['2/m',]:
1471	if SGData['SGUniq'] == 'a':
1472	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1473	elif SGData['SGUniq'] == 'b':
1474	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1475	else:
1476	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1477	elif SGData['SGLaue'] in ['mmm',]:
1478	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1479	elif SGData['SGLaue'] in ['4/m','4/mmm']:
1480	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1481	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1482	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1483	elif SGData['SGLaue'] in ['3R', '3mR']:
1484	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
1485	elif SGData['SGLaue'] in ['m3m','m3']:
1486	return [[pfx+'A0',dpdA[0]]]
1487
1488	# create a list of dependent variables and set up a dictionary to hold their derivatives
1489	dependentVars = G2mv.GetDependentVars()
1490	depDerivDict = {}
1491	for j in dependentVars:
1492	depDerivDict[j] = np.zeros(shape=(len(x)))
1493	#print 'dependent vars',dependentVars
1494	lenX = len(x)
1495	hId = Histogram['hId']
1496	hfx = ':%d:'%(hId)
1497	bakType = calcControls[hfx+'bakType']
1498	dMdv = np.zeros(shape=(len(varylist),len(x)))
1499	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1500	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1501	bBpos =varylist.index(hfx+'Back:0')
1502	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1503	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1504	for name in varylist:
1505	if 'Debye' in name:
1506	id = int(name.split(':')[-1])
1507	parm = name[:int(name.rindex(':'))]
1508	ip = names.index(parm)
1509	dMdv[varylist.index(name)] = dMddb[3*id+ip]
1510	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1511	for name in varylist:
1512	if 'BkPk' in name:
1513	parm,id = name.split(';')
1514	id = int(id)
1515	if parm in names:
1516	ip = names.index(parm)
1517	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1518	cw = np.diff(x)
1519	cw = np.append(cw,cw[-1])
1520	if 'C' in calcControls[hfx+'histType']:
1521	shl = max(parmDict[hfx+'SH/L'],0.002)
1522	Ka2 = False
1523	if hfx+'Lam1' in parmDict.keys():
1524	wave = parmDict[hfx+'Lam1']
1525	Ka2 = True
1526	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
1527	kRatio = parmDict[hfx+'I(L2)/I(L1)']
1528	else:
1529	wave = parmDict[hfx+'Lam']
1530	else:
1531	print 'TOF Undefined at present'
1532	raise ValueError
1533	for phase in Histogram['Reflection Lists']:
1534	refDict = Histogram['Reflection Lists'][phase]
1535	Phase = Phases[phase]
1536	SGData = Phase['General']['SGData']
1537	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1538	pId = Phase['pId']
1539	pfx = '%d::'%(pId)
1540	phfx = '%d:%d:'%(pId,hId)
1541	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] #And modify here by Dij?
1542	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1543	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
1544	if not Phase['General'].get('doPawley'):
1545	time0 = time.time()
1546	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1547	# print 'sf-derv time %.3fs'%(time.time()-time0)
1548	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1549	time0 = time.time()
1550	for iref,refl in enumerate(refDict['RefList']):
1551	if 'C' in calcControls[hfx+'histType']: #CW powder
1552	h,k,l = refl[:3]
1553	Uniq = np.inner(refl[:3],SGMT)
1554	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1555	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1556	iBeg = np.searchsorted(x,refl[5]-fmin)
1557	iFin = np.searchsorted(x,refl[5]+fmax)
1558	if not iBeg+iFin: #peak below low limit - skip peak
1559	continue
1560	elif not iBeg-iFin: #peak above high limit - done
1561	break
1562	pos = refl[5]
1563	tanth = tand(pos/2.0)
1564	costh = cosd(pos/2.0)
1565	lenBF = iFin-iBeg
1566	dMdpk = np.zeros(shape=(6,lenBF))
1567	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1568	for i in range(5):
1569	dMdpk[i] += 100.cw[iBeg:iFin]refl[11]refl[9]dMdipk[i]
1570	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1571	if Ka2:
1572	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
1573	iBeg2 = np.searchsorted(x,pos2-fmin)
1574	iFin2 = np.searchsorted(x,pos2+fmax)
1575	if iBeg2-iFin2:
1576	lenBF2 = iFin2-iBeg2
1577	dMdpk2 = np.zeros(shape=(6,lenBF2))
1578	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1579	for i in range(5):
1580	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[11]refl[9]kRatio*dMdipk2[i]
1581	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[11]*dMdipk2[0]
1582	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1583	if Phase['General'].get('doPawley'):
1584	dMdpw = np.zeros(len(x))
1585	try:
1586	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1587	idx = varylist.index(pIdx)
1588	dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1589	if Ka2:
1590	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1591	dMdv[idx] = dMdpw
1592	except: # ValueError:
1593	pass
1594	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1595	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1596	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1597	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1598	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1599	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1600	hfx+'DisplaceY':[dpdY,'pos'],}
1601	if 'Bragg' in calcControls[hfx+'instType']:
1602	names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1603	hfx+'SurfRoughB':[dFdAb[1],'int'],})
1604	else:
1605	names.update({hfx+'Absorption':[dFdAb,'int'],})
1606	for name in names:
1607	item = names[name]
1608	if name in varylist:
1609	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1610	if Ka2:
1611	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1612	elif name in dependentVars:
1613	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1614	if Ka2:
1615	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1616	for iPO in dIdPO:
1617	if iPO in varylist:
1618	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1619	if Ka2:
1620	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1621	elif iPO in dependentVars:
1622	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1623	if Ka2:
1624	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1625	for i,name in enumerate(['omega','chi','phi']):
1626	aname = pfx+'SH '+name
1627	if aname in varylist:
1628	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1629	if Ka2:
1630	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1631	elif aname in dependentVars:
1632	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1633	if Ka2:
1634	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1635	for iSH in dFdODF:
1636	if iSH in varylist:
1637	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1638	if Ka2:
1639	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1640	elif iSH in dependentVars:
1641	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1642	if Ka2:
1643	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1644	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1645	for name,dpdA in cellDervNames:
1646	if name in varylist:
1647	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1648	if Ka2:
1649	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1650	elif name in dependentVars:
1651	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1652	if Ka2:
1653	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1654	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1655	for name in dDijDict:
1656	if name in varylist:
1657	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1658	if Ka2:
1659	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1660	elif name in dependentVars:
1661	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1662	if Ka2:
1663	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1664	sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
1665	for name in gamDict:
1666	if name in varylist:
1667	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1668	if Ka2:
1669	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1670	elif name in dependentVars:
1671	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1672	if Ka2:
1673	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1674	for name in sigDict:
1675	if name in varylist:
1676	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1677	if Ka2:
1678	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1679	elif name in dependentVars:
1680	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1681	if Ka2:
1682	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1683	for name in ['BabA','BabU']:
1684	if refl[9]:
1685	if phfx+name in varylist:
1686	dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1687	if Ka2:
1688	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1689	elif phfx+name in dependentVars:
1690	depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1691	if Ka2:
1692	depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1693	elif 'T' in calcControls[hfx+'histType']:
1694	print 'TOF Undefined at present'
1695	raise Exception #no TOF yet
1696	if not Phase['General'].get('doPawley'):
1697	#do atom derivatives - for RB,F,X & U so far
1698	corr = dervDict['int']/refl[9]
1699	if Ka2:
1700	corr2 = dervDict2['int']/refl[9]
1701	for name in varylist+dependentVars:
1702	if '::RBV;' in name:
1703	pass
1704	else:
1705	try:
1706	aname = name.split(pfx)[1][:2]
1707	if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1708	except IndexError:
1709	continue
1710	if name in varylist:
1711	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1712	if Ka2:
1713	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1714	elif name in dependentVars:
1715	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1716	if Ka2:
1717	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1718	# print 'profile derv time: %.3fs'%(time.time()-time0)
1719	# now process derivatives in constraints
1720	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1721	return dMdv
1722
1723	def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
1724	'''Loop over reflections ina HKLF histogram and compute derivatives of the fitting
1725	model (M) with respect to all parameters. Independent and dependant dM/dp arrays
1726	are returned to either dervRefine or HessRefine.
1727
1728	:returns:
1729	'''
1730	nobs = Histogram['Residuals']['Nobs']
1731	hId = Histogram['hId']
1732	hfx = ':%d:'%(hId)
1733	pfx = '%d::'%(Phase['pId'])
1734	phfx = '%d:%d:'%(Phase['pId'],hId)
1735	SGData = Phase['General']['SGData']
1736	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1737	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
1738	refDict = Histogram['Data']
1739	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1740	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1741	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1742	dependentVars = G2mv.GetDependentVars()
1743	depDerivDict = {}
1744	for j in dependentVars:
1745	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
1746	wdf = np.zeros(len(refDict['RefList']))
1747	if calcControls['F**2']:
1748	for iref,ref in enumerate(refDict['RefList']):
1749	if ref[6] > 0:
1750	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
1751	w = 1.0/ref[6]
1752	if w*ref[5] >= calcControls['minF/sig']:
1753	wdf[iref] = w*(ref[5]-ref[7])
1754	for j,var in enumerate(varylist):
1755	if var in dFdvDict:
1756	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1757	for var in dependentVars:
1758	if var in dFdvDict:
1759	depDerivDict[var][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1760	if phfx+'Scale' in varylist:
1761	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
1762	elif phfx+'Scale' in dependentVars:
1763	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
1764	for item in ['Ep','Es','Eg']:
1765	if phfx+item in varylist and dervDict:
1766	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/dervCor
1767	elif phfx+item in dependentVars and dervDict:
1768	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/dervCor
1769	for item in ['BabA','BabU']:
1770	if phfx+item in varylist:
1771	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1772	elif phfx+item in dependentVars:
1773	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1774	else:
1775	for iref,ref in enumerate(refDict['RefList']):
1776	if ref[5] > 0.:
1777	dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
1778	Fo = np.sqrt(ref[5])
1779	Fc = np.sqrt(ref[7])
1780	w = 1.0/ref[6]
1781	if 2.0Fow*Fo >= calcControls['minF/sig']:
1782	wdf[iref] = 2.0Fow*(Fo-Fc)
1783	for j,var in enumerate(varylist):
1784	if var in dFdvDict:
1785	dMdvh[j][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1786	for var in dependentVars:
1787	if var in dFdvDict:
1788	depDerivDict[var][iref] = wdFdvDict[var][iref]dervCor*parmDict[phfx+'Scale']
1789	if phfx+'Scale' in varylist:
1790	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[9]dervCor
1791	elif phfx+'Scale' in dependentVars:
1792	depDerivDict[phfx+'Scale'][iref] = wref[9]dervCor
1793	for item in ['Ep','Es','Eg']:
1794	if phfx+item in varylist and dervDict:
1795	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/dervCor #correct
1796	elif phfx+item in dependentVars and dervDict:
1797	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/dervCor
1798	for item in ['BabA','BabU']:
1799	if phfx+item in varylist:
1800	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1801	elif phfx+item in dependentVars:
1802	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*dervCor
1803	return dMdvh,depDerivDict,wdf
1804
1805
1806	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1807	'''Loop over histograms and compute derivatives of the fitting
1808	model (M) with respect to all parameters. Results are returned in
1809	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
1810	n is the number of parameters and m is the number of data
1811	points. This can exceed memory when m gets large. This routine is
1812	used when refinement derivatives are selected as "analtytic
1813	Jacobian" in Controls.
1814
1815	:returns: Jacobian numpy.array dMdv for all histograms concatinated
1816	'''
1817	parmDict.update(zip(varylist,values))
1818	G2mv.Dict2Map(parmDict,varylist)
1819	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1820	nvar = len(varylist)
1821	dMdv = np.empty(0)
1822	histoList = Histograms.keys()
1823	histoList.sort()
1824	for histogram in histoList:
1825	if 'PWDR' in histogram[:4]:
1826	Histogram = Histograms[histogram]
1827	hId = Histogram['hId']
1828	hfx = ':%d:'%(hId)
1829	wtFactor = calcControls[hfx+'wtFactor']
1830	Limits = calcControls[hfx+'Limits']
1831	x,y,w,yc,yb,yd = Histogram['Data']
1832	W = wtFactor*w
1833	xB = np.searchsorted(x,Limits[0])
1834	xF = np.searchsorted(x,Limits[1])
1835	dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
1836	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1837	elif 'HKLF' in histogram[:4]:
1838	Histogram = Histograms[histogram]
1839	phase = Histogram['Reflection Lists']
1840	Phase = Phases[phase]
1841	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
1842	hfx = ':%d:'%(Histogram['hId'])
1843	wtFactor = calcControls[hfx+'wtFactor']
1844	# now process derivatives in constraints
1845	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1846	else:
1847	continue #skip non-histogram entries
1848	if len(dMdv):
1849	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
1850	else:
1851	dMdv = np.sqrt(wtFactor)*dMdvh
1852
1853	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
1854	if np.any(pVals):
1855	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1856	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
1857
1858	return dMdv
1859
1860	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1861	'''Loop over histograms and compute derivatives of the fitting
1862	model (M) with respect to all parameters. For each histogram, the
1863	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
1864	number of parameters and m is the number of data points *in the
1865	histogram*. The (n by n) Hessian is computed from each Jacobian
1866	and it is returned. This routine is used when refinement
1867	derivatives are selected as "analtytic Hessian" in Controls.
1868
1869	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
1870	'''
1871	parmDict.update(zip(varylist,values))
1872	G2mv.Dict2Map(parmDict,varylist)
1873	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1874	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
1875	nvar = len(varylist)
1876	Hess = np.empty(0)
1877	histoList = Histograms.keys()
1878	histoList.sort()
1879	for histogram in histoList:
1880	if 'PWDR' in histogram[:4]:
1881	Histogram = Histograms[histogram]
1882	hId = Histogram['hId']
1883	hfx = ':%d:'%(hId)
1884	wtFactor = calcControls[hfx+'wtFactor']
1885	Limits = calcControls[hfx+'Limits']
1886	x,y,w,yc,yb,yd = Histogram['Data']
1887	W = wtFactor*w
1888	dy = y-yc
1889	xB = np.searchsorted(x,Limits[0])
1890	xF = np.searchsorted(x,Limits[1])
1891	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
1892	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1893	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
1894	Dy = dy[xB:xF][np.newaxis,:]
1895	dMdvh *= Wt
1896	if dlg:
1897	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1898	if len(Hess):
1899	Hess += np.inner(dMdvh,dMdvh)
1900	dMdvh = WtDy
1901	Vec += np.sum(dMdvh,axis=1)
1902	else:
1903	Hess = np.inner(dMdvh,dMdvh)
1904	dMdvh = WtDy
1905	Vec = np.sum(dMdvh,axis=1)
1906	elif 'HKLF' in histogram[:4]:
1907	Histogram = Histograms[histogram]
1908	phase = Histogram['Reflection Lists']
1909	Phase = Phases[phase]
1910	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
1911	hId = Histogram['hId']
1912	hfx = ':%d:'%(Histogram['hId'])
1913	wtFactor = calcControls[hfx+'wtFactor']
1914	# now process derivatives in constraints
1915	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1916	# print 'matrix build time: %.3f'%(time.time()-time0)
1917
1918	if dlg:
1919	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1920	if len(Hess):
1921	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
1922	Hess += wtFactor*np.inner(dMdvh,dMdvh)
1923	else:
1924	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
1925	Hess = wtFactor*np.inner(dMdvh,dMdvh)
1926	else:
1927	continue #skip non-histogram entries
1928	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
1929	if np.any(pVals):
1930	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1931	Vec += np.sum(dpdvpWtpVals,axis=1)
1932	Hess += np.inner(dpdv*pWt,dpdv)
1933	return Vec,Hess
1934
1935	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1936	'Needs a doc string'
1937	Values2Dict(parmDict, varylist, values)
1938	G2mv.Dict2Map(parmDict,varylist)
1939	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1940	M = np.empty(0)
1941	SumwYo = 0
1942	Nobs = 0
1943	ApplyRBModels(parmDict,Phases,rigidbodyDict)
1944	histoList = Histograms.keys()
1945	histoList.sort()
1946	for histogram in histoList:
1947	if 'PWDR' in histogram[:4]:
1948	Histogram = Histograms[histogram]
1949	hId = Histogram['hId']
1950	hfx = ':%d:'%(hId)
1951	wtFactor = calcControls[hfx+'wtFactor']
1952	Limits = calcControls[hfx+'Limits']
1953	x,y,w,yc,yb,yd = Histogram['Data']
1954	yc *= 0.0 #zero full calcd profiles
1955	yb *= 0.0
1956	yd *= 0.0
1957	xB = np.searchsorted(x,Limits[0])
1958	xF = np.searchsorted(x,Limits[1])
1959	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
1960	varylist,Histogram,Phases,calcControls,pawleyLookup)
1961	yc[xB:xF] += yb[xB:xF]
1962	if not np.any(y): #fill dummy data
1963	rv = st.poisson(yc[xB:xF])
1964	y[xB:xF] = rv.rvs()
1965	Z = np.ones_like(yc[xB:xF])
1966	Z[1::2] *= -1
1967	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
1968	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
1969	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
1970	W = wtFactor*w
1971	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
1972	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
1973	Nobs += Histogram['Residuals']['Nobs']
1974	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
1975	SumwYo += Histogram['Residuals']['sumwYo']
1976	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
1977	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
1978	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
1979	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
1980	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.))
1981	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.))
1982	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
1983	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
1984	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
1985	if dlg:
1986	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1987	M = np.concatenate((M,wdy))
1988	#end of PWDR processing
1989	elif 'HKLF' in histogram[:4]:
1990	Histogram = Histograms[histogram]
1991	Histogram['Residuals'] = {}
1992	phase = Histogram['Reflection Lists']
1993	Phase = Phases[phase]
1994	hId = Histogram['hId']
1995	hfx = ':%d:'%(hId)
1996	wtFactor = calcControls[hfx+'wtFactor']
1997	pfx = '%d::'%(Phase['pId'])
1998	phfx = '%d:%d:'%(Phase['pId'],hId)
1999	SGData = Phase['General']['SGData']
2000	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2001	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2002	refDict = Histogram['Data']
2003	time0 = time.time()
2004	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2005	# print 'sf-calc time: %.3f'%(time.time()-time0)
2006	df = np.zeros(len(refDict['RefList']))
2007	sumwYo = 0
2008	sumFo = 0
2009	sumFo2 = 0
2010	sumdF = 0
2011	sumdF2 = 0
2012	nobs = 0
2013	if calcControls['F**2']:
2014	for i,ref in enumerate(refDict['RefList']):
2015	if ref[6] > 0:
2016	SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
2017	w = 1.0/ref[6]
2018	ref[7] = parmDict[phfx+'Scale']ref[9]ref[11] #correct Fc^2 for extinction
2019	ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2020	if w*ref[5] >= calcControls['minF/sig']:
2021	Fo = np.sqrt(ref[5])
2022	sumFo += Fo
2023	sumFo2 += ref[5]
2024	sumdF += abs(Fo-np.sqrt(ref[7]))
2025	sumdF2 += abs(ref[5]-ref[7])
2026	nobs += 1
2027	df[i] = -w*(ref[5]-ref[7])
2028	sumwYo += (wref[5])*2
2029	else:
2030	for i,ref in enumerate(refDict['RefList']):
2031	if ref[5] > 0.:
2032	SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[11]
2033	ref[7] = parmDict[phfx+'Scale']ref[9]ref[11] #correct Fc^2 for extinctio
2034	ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
2035	Fo = np.sqrt(ref[5])
2036	Fc = np.sqrt(ref[7])
2037	w = 2.0*Fo/ref[6]
2038	if w*Fo >= calcControls['minF/sig']:
2039	sumFo += Fo
2040	sumFo2 += ref[5]
2041	sumdF += abs(Fo-Fc)
2042	sumdF2 += abs(ref[5]-ref[7])
2043	nobs += 1
2044	df[i] = -w*(Fo-Fc)
2045	sumwYo += (wFo)*2
2046	Histogram['Residuals']['Nobs'] = nobs
2047	Histogram['Residuals']['sumwYo'] = sumwYo
2048	SumwYo += sumwYo
2049	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/Histogram['Residuals']['sumwYo'])100.)
2050	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2051	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2052	Histogram['Residuals'][phfx+'Nref'] = nobs
2053	Nobs += nobs
2054	if dlg:
2055	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2056	M = np.concatenate((M,wtFactor*df))
2057	# end of HKLF processing
2058	Histograms['sumwYo'] = SumwYo
2059	Histograms['Nobs'] = Nobs
2060	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
2061	if dlg:
2062	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
2063	if not GoOn:
2064	parmDict['saved values'] = values
2065	dlg.Destroy()
2066	raise Exception #Abort!!
2067	pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
2068	if len(pVals):
2069	pSum = np.sum(pWtpVals*2)
2070	for name in pWsum:
2071	if pWsum:
2072	print ' Penalty function for %8s = %12.5g'%(name,pWsum[name])
2073	print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
2074	Nobs += len(pVals)
2075	M = np.concatenate((M,np.sqrt(pWt)*pVals))
2076	return M
2077

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

Download in other formats: