Context Navigation

source: trunk/GSASIIstrMath.py @ 2047

Last change on this file since 2047 was 2046, checked in by vondreele, 10 years ago
some work on SS derivs.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 160.5 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMath - structure math routines
4	-----------------------------------------
5	'''
6	########### SVN repository information ###################
7	# $Date: 2015-11-11 21:48:12 +0000 (Wed, 11 Nov 2015) $
8	# $Author: toby $
9	# $Revision: 2046 $
10	# $URL: trunk/GSASIIstrMath.py $
11	# $Id: GSASIIstrMath.py 2046 2015-11-11 21:48:12Z toby $
12	########### SVN repository information ###################
13	import time
14	import copy
15	import numpy as np
16	import numpy.ma as ma
17	import numpy.linalg as nl
18	import scipy.optimize as so
19	import scipy.stats as st
20	import GSASIIpath
21	GSASIIpath.SetVersionNumber("$Revision: 2046 $")
22	import GSASIIElem as G2el
23	import GSASIIlattice as G2lat
24	import GSASIIspc as G2spc
25	import GSASIIpwd as G2pwd
26	import GSASIImapvars as G2mv
27	import GSASIImath as G2mth
28	import GSASIIobj as G2obj
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*np.log(2.0)
38	twopi = 2.0*np.pi
39	twopisq = 2.0np.pi*2
40
41	################################################################################
42	##### Rigid Body Models
43	################################################################################
44
45	def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
46	''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
47	current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
48	'''
49	atxIds = ['Ax:','Ay:','Az:']
50	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
51	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
52	if not RBIds['Vector'] and not RBIds['Residue']:
53	return
54	VRBIds = RBIds['Vector']
55	RRBIds = RBIds['Residue']
56	if Update:
57	RBData = rigidbodyDict
58	else:
59	RBData = copy.deepcopy(rigidbodyDict) # don't mess with original!
60	if RBIds['Vector']: # first update the vector magnitudes
61	VRBData = RBData['Vector']
62	for i,rbId in enumerate(VRBIds):
63	if VRBData[rbId]['useCount']:
64	for j in range(len(VRBData[rbId]['VectMag'])):
65	name = '::RBV;'+str(j)+':'+str(i)
66	VRBData[rbId]['VectMag'][j] = parmDict[name]
67	for phase in Phases:
68	Phase = Phases[phase]
69	General = Phase['General']
70	cx,ct,cs,cia = General['AtomPtrs']
71	cell = General['Cell'][1:7]
72	Amat,Bmat = G2lat.cell2AB(cell)
73	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
74	pfx = str(Phase['pId'])+'::'
75	if Update:
76	RBModels = Phase['RBModels']
77	else:
78	RBModels = copy.deepcopy(Phase['RBModels']) # again don't mess with original!
79	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
80	jrb = VRBIds.index(RBObj['RBId'])
81	rbsx = str(irb)+':'+str(jrb)
82	for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
83	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
84	for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
85	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
86	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
87	TLS = RBObj['ThermalMotion']
88	if 'T' in TLS[0]:
89	for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
90	TLS[1][i] = parmDict[pfx+pt+rbsx]
91	if 'L' in TLS[0]:
92	for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
93	TLS[1][i+6] = parmDict[pfx+pt+rbsx]
94	if 'S' in TLS[0]:
95	for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
96	TLS[1][i+12] = parmDict[pfx+pt+rbsx]
97	if 'U' in TLS[0]:
98	TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
99	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
100	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
101	for i,x in enumerate(XYZ):
102	atId = RBObj['Ids'][i]
103	for j in [0,1,2]:
104	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
105	if UIJ[i][0] == 'A':
106	for j in range(6):
107	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
108	elif UIJ[i][0] == 'I':
109	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
110
111	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
112	jrb = RRBIds.index(RBObj['RBId'])
113	rbsx = str(irb)+':'+str(jrb)
114	for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
115	RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
116	for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
117	RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
118	RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
119	TLS = RBObj['ThermalMotion']
120	if 'T' in TLS[0]:
121	for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
122	RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
123	if 'L' in TLS[0]:
124	for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
125	RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
126	if 'S' in TLS[0]:
127	for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
128	RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
129	if 'U' in TLS[0]:
130	RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
131	for itors,tors in enumerate(RBObj['Torsions']):
132	tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
133	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
134	UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
135	for i,x in enumerate(XYZ):
136	atId = RBObj['Ids'][i]
137	for j in [0,1,2]:
138	parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
139	if UIJ[i][0] == 'A':
140	for j in range(6):
141	parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
142	elif UIJ[i][0] == 'I':
143	parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
144
145	def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
146	'Needs a doc string'
147	atxIds = ['dAx:','dAy:','dAz:']
148	atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
149	TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
150	LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
151	SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
152	PIds = ['Px:','Py:','Pz:']
153	OIds = ['Oa:','Oi:','Oj:','Ok:']
154	RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) #these are lists of rbIds
155	if not RBIds['Vector'] and not RBIds['Residue']:
156	return
157	VRBIds = RBIds['Vector']
158	RRBIds = RBIds['Residue']
159	RBData = rigidbodyDict
160	for item in parmDict:
161	if 'RB' in item:
162	dFdvDict[item] = 0. #NB: this is a vector which is no. refl. long & must be filled!
163	General = Phase['General']
164	cx,ct,cs,cia = General['AtomPtrs']
165	cell = General['Cell'][1:7]
166	Amat,Bmat = G2lat.cell2AB(cell)
167	rpd = np.pi/180.
168	rpd2 = rpd**2
169	g = nl.inv(np.inner(Bmat,Bmat))
170	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
171	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
172	AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
173	pfx = str(Phase['pId'])+'::'
174	RBModels = Phase['RBModels']
175	for irb,RBObj in enumerate(RBModels.get('Vector',[])):
176	VModel = RBData['Vector'][RBObj['RBId']]
177	Q = RBObj['Orient'][0]
178	Pos = RBObj['Orig'][0]
179	jrb = VRBIds.index(RBObj['RBId'])
180	rbsx = str(irb)+':'+str(jrb)
181	dXdv = []
182	for iv in range(len(VModel['VectMag'])):
183	dCdv = []
184	for vec in VModel['rbVect'][iv]:
185	dCdv.append(G2mth.prodQVQ(Q,vec))
186	dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
187	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
188	for ia,atId in enumerate(RBObj['Ids']):
189	atNum = AtLookup[atId]
190	dx = 0.00001
191	for iv in range(len(VModel['VectMag'])):
192	for ix in [0,1,2]:
193	dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
194	for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
195	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
196	for iv in range(4):
197	Q[iv] -= dx
198	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
199	Q[iv] += 2.*dx
200	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
201	Q[iv] -= dx
202	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
203	for ix in [0,1,2]:
204	dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
205	X = G2mth.prodQVQ(Q,Cart[ia])
206	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
207	dFdu = G2lat.U6toUij(dFdu.T)
208	dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
209	dFdu = G2lat.UijtoU6(dFdu)
210	atNum = AtLookup[atId]
211	if 'T' in RBObj['ThermalMotion'][0]:
212	for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
213	dFdvDict[pfx+name+rbsx] += dFdu[i]
214	if 'L' in RBObj['ThermalMotion'][0]:
215	dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
216	dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
217	dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
218	dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
219	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
220	dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2(-dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
221	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
222	dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2(-dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
223	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
224	if 'S' in RBObj['ThermalMotion'][0]:
225	dFdvDict[pfx+'RBVS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
226	dFdvDict[pfx+'RBVS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
227	dFdvDict[pfx+'RBVS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
228	dFdvDict[pfx+'RBVS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
229	dFdvDict[pfx+'RBVS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
230	dFdvDict[pfx+'RBVS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
231	dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
232	dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
233	if 'U' in RBObj['ThermalMotion'][0]:
234	dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
235
236
237	for irb,RBObj in enumerate(RBModels.get('Residue',[])):
238	Q = RBObj['Orient'][0]
239	Pos = RBObj['Orig'][0]
240	jrb = RRBIds.index(RBObj['RBId'])
241	torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
242	rbsx = str(irb)+':'+str(jrb)
243	XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
244	for itors,tors in enumerate(RBObj['Torsions']): #derivative error?
245	tname = pfx+'RBRTr;'+str(itors)+':'+rbsx
246	orId,pvId = torData[itors][:2]
247	pivotVec = Cart[orId]-Cart[pvId]
248	QA = G2mth.AVdeg2Q(-0.001,pivotVec)
249	QB = G2mth.AVdeg2Q(0.001,pivotVec)
250	for ir in torData[itors][3]:
251	atNum = AtLookup[RBObj['Ids'][ir]]
252	rVec = Cart[ir]-Cart[pvId]
253	dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
254	dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
255	for ix in [0,1,2]:
256	dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
257	for ia,atId in enumerate(RBObj['Ids']):
258	atNum = AtLookup[atId]
259	dx = 0.00001
260	for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
261	dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
262	for iv in range(4):
263	Q[iv] -= dx
264	XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
265	Q[iv] += 2.*dx
266	XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
267	Q[iv] -= dx
268	dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
269	for ix in [0,1,2]:
270	dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
271	X = G2mth.prodQVQ(Q,Cart[ia])
272	dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
273	dFdu = G2lat.U6toUij(dFdu.T)
274	dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
275	dFdu = G2lat.UijtoU6(dFdu)
276	atNum = AtLookup[atId]
277	if 'T' in RBObj['ThermalMotion'][0]:
278	for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
279	dFdvDict[pfx+name+rbsx] += dFdu[i]
280	if 'L' in RBObj['ThermalMotion'][0]:
281	dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2(dFdu[1]X[2]*2+dFdu[2]X[1]*2-dFdu[5]X[1]*X[2])
282	dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2(dFdu[0]X[2]*2+dFdu[2]X[0]*2-dFdu[4]X[0]*X[2])
283	dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2(dFdu[0]X[1]*2+dFdu[1]X[0]*2-dFdu[3]X[0]*X[1])
284	dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2(-dFdu[3]X[2]*2-2.dFdu[2]X[0]X[1]+
285	dFdu[4]X[1]X[2]+dFdu[5]X[0]X[2])
286	dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2(dFdu[4]X[1]*2-2.dFdu[1]X[0]X[2]+
287	dFdu[3]X[1]X[2]+dFdu[5]X[0]X[1])
288	dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2(dFdu[5]X[0]*2-2.dFdu[0]X[1]X[2]+
289	dFdu[3]X[0]X[2]+dFdu[4]X[0]X[1])
290	if 'S' in RBObj['ThermalMotion'][0]:
291	dFdvDict[pfx+'RBRS12:'+rbsx] += rpd(dFdu[5]X[1]-2.dFdu[1]X[2])
292	dFdvDict[pfx+'RBRS13:'+rbsx] += rpd(-dFdu[5]X[2]+2.dFdu[2]X[1])
293	dFdvDict[pfx+'RBRS21:'+rbsx] += rpd(-dFdu[4]X[0]+2.dFdu[0]X[2])
294	dFdvDict[pfx+'RBRS23:'+rbsx] += rpd(dFdu[4]X[2]-2.dFdu[2]X[0])
295	dFdvDict[pfx+'RBRS31:'+rbsx] += rpd(dFdu[3]X[0]-2.dFdu[0]X[1])
296	dFdvDict[pfx+'RBRS32:'+rbsx] += rpd(-dFdu[3]X[1]+2.dFdu[1]X[0])
297	dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd(dFdu[4]X[1]-dFdu[3]*X[2])
298	dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd(dFdu[5]X[0]-dFdu[3]*X[2])
299	if 'U' in RBObj['ThermalMotion'][0]:
300	dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
301
302	################################################################################
303	##### Penalty & restraint functions
304	################################################################################
305
306	def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
307	'Needs a doc string'
308	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
309	pNames = []
310	pVals = []
311	pWt = []
312	negWt = {}
313	pWsum = {}
314	for phase in Phases:
315	pId = Phases[phase]['pId']
316	negWt[pId] = Phases[phase]['General']['Pawley neg wt']
317	General = Phases[phase]['General']
318	cx,ct,cs,cia = General['AtomPtrs']
319	textureData = General['SH Texture']
320	SGData = General['SGData']
321	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
322	cell = General['Cell'][1:7]
323	Amat,Bmat = G2lat.cell2AB(cell)
324	if phase not in restraintDict:
325	continue
326	phaseRest = restraintDict[phase]
327	names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
328	['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
329	['ChemComp','Sites'],['Texture','HKLs'],]
330	for name,rest in names:
331	pWsum[name] = 0.
332	itemRest = phaseRest[name]
333	if itemRest[rest] and itemRest['Use']:
334	wt = itemRest['wtFactor']
335	if name in ['Bond','Angle','Plane','Chiral']:
336	for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
337	pNames.append(str(pId)+':'+name+':'+str(i))
338	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
339	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
340	if name == 'Bond':
341	calc = G2mth.getRestDist(XYZ,Amat)
342	elif name == 'Angle':
343	calc = G2mth.getRestAngle(XYZ,Amat)
344	elif name == 'Plane':
345	calc = G2mth.getRestPlane(XYZ,Amat)
346	elif name == 'Chiral':
347	calc = G2mth.getRestChiral(XYZ,Amat)
348	pVals.append(obs-calc)
349	pWt.append(wt/esd**2)
350	pWsum[name] += wt((obs-calc)/esd)*2
351	elif name in ['Torsion','Rama']:
352	coeffDict = itemRest['Coeff']
353	for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
354	pNames.append(str(pId)+':'+name+':'+str(i))
355	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
356	XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
357	if name == 'Torsion':
358	tor = G2mth.getRestTorsion(XYZ,Amat)
359	restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
360	else:
361	phi,psi = G2mth.getRestRama(XYZ,Amat)
362	restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
363	pVals.append(restr)
364	pWt.append(wt/esd**2)
365	pWsum[name] += wt(restr/esd)*2
366	elif name == 'ChemComp':
367	for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
368	pNames.append(str(pId)+':'+name+':'+str(i))
369	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
370	frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
371	calc = np.sum(mulfracfactors)
372	pVals.append(obs-calc)
373	pWt.append(wt/esd**2)
374	pWsum[name] += wt((obs-calc)/esd)*2
375	elif name == 'Texture':
376	SHkeys = textureData['SH Coeff'][1].keys()
377	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
378	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
379	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
380	for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
381	PH = np.array(hkl)
382	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
383	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
384	R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
385	Z1 = -ma.masked_greater(Z,0.0)
386	IndZ1 = np.array(ma.nonzero(Z1))
387	for ind in IndZ1.T:
388	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
389	pVals.append(Z1[ind[0]][ind[1]])
390	pWt.append(wt/esd1**2)
391	pWsum[name] += wt(-Z1[ind[0]][ind[1]]/esd1)*2
392	if ifesd2:
393	Z2 = 1.-Z
394	for ind in np.ndindex(grid,grid):
395	pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
396	pVals.append(Z1[ind[0]][ind[1]])
397	pWt.append(wt/esd2**2)
398	pWsum[name] += wt(Z2/esd2)*2
399
400	for phase in Phases:
401	name = 'SH-Pref.Ori.'
402	pId = Phases[phase]['pId']
403	General = Phases[phase]['General']
404	SGData = General['SGData']
405	cell = General['Cell'][1:7]
406	pWsum[name] = 0.0
407	for hist in Phases[phase]['Histograms']:
408	if hist in Histograms and 'PWDR' in hist:
409	hId = Histograms[hist]['hId']
410	phfx = '%d:%d:'%(pId,hId)
411	if calcControls[phfx+'poType'] == 'SH':
412	toler = calcControls[phfx+'SHtoler']
413	wt = 1./toler**2
414	HKLs = np.array(calcControls[phfx+'SHhkl'])
415	SHnames = calcControls[phfx+'SHnames']
416	SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
417	for i,PH in enumerate(HKLs):
418	phi,beta = G2lat.CrsAng(PH,cell,SGData)
419	SH3Coef = {}
420	for item in SHcof:
421	L,N = eval(item.strip('C'))
422	SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]
423	ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
424	X = np.linspace(0,90.0,26)
425	Y = -ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)
426	IndY = ma.nonzero(Y)
427	for ind in IndY[0]:
428	pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
429	pVals.append(Y[ind])
430	pWt.append(wt)
431	pWsum[name] += wt(Y[ind])*2
432	pWsum['PWLref'] = 0.
433	for item in varyList:
434	if 'PWLref' in item and parmDict[item] < 0.:
435	pId = int(item.split(':')[0])
436	if negWt[pId]:
437	pNames.append(item)
438	pVals.append(-parmDict[item])
439	pWt.append(negWt[pId])
440	pWsum['PWLref'] += negWt[pId](-parmDict[item])*2
441	pVals = np.array(pVals)
442	pWt = np.array(pWt) #should this be np.sqrt?
443	return pNames,pVals,pWt,pWsum
444
445	def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
446	'Needs a doc string'
447	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
448	pDerv = np.zeros((len(varyList),len(pVal)))
449	for phase in Phases:
450	# if phase not in restraintDict:
451	# continue
452	pId = Phases[phase]['pId']
453	General = Phases[phase]['General']
454	cx,ct,cs,cia = General['AtomPtrs']
455	SGData = General['SGData']
456	AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
457	cell = General['Cell'][1:7]
458	Amat,Bmat = G2lat.cell2AB(cell)
459	textureData = General['SH Texture']
460
461	SHkeys = textureData['SH Coeff'][1].keys()
462	SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
463	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
464	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
465	sam = SamSym[textureData['Model']]
466	phaseRest = restraintDict.get(phase,{})
467	names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
468	'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
469	'ChemComp':'Sites','Texture':'HKLs'}
470	lasthkl = np.array([0,0,0])
471	for ip,pName in enumerate(pNames):
472	pnames = pName.split(':')
473	if pId == int(pnames[0]):
474	name = pnames[1]
475	if 'PWL' in pName:
476	pDerv[varyList.index(pName)][ip] += 1.
477	continue
478	elif 'SH-' in pName:
479	continue
480	id = int(pnames[2])
481	itemRest = phaseRest[name]
482	if name in ['Bond','Angle','Plane','Chiral']:
483	indx,ops,obs,esd = itemRest[names[name]][id]
484	dNames = []
485	for ind in indx:
486	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
487	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
488	if name == 'Bond':
489	deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
490	elif name == 'Angle':
491	deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
492	elif name == 'Plane':
493	deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
494	elif name == 'Chiral':
495	deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
496	elif name in ['Torsion','Rama']:
497	coffDict = itemRest['Coeff']
498	indx,ops,cofName,esd = itemRest[names[name]][id]
499	dNames = []
500	for ind in indx:
501	dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
502	XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
503	if name == 'Torsion':
504	deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
505	else:
506	deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
507	elif name == 'ChemComp':
508	indx,factors,obs,esd = itemRest[names[name]][id]
509	dNames = []
510	for ind in indx:
511	dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
512	mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
513	deriv = mul*factors
514	elif 'Texture' in name:
515	deriv = []
516	dNames = []
517	hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
518	hkl = np.array(hkl)
519	if np.any(lasthkl-hkl):
520	PH = np.array(hkl)
521	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
522	ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
523	lasthkl = copy.copy(hkl)
524	if 'unit' in name:
525	pass
526	else:
527	gam = float(pnames[3])
528	psi = float(pnames[4])
529	for SHname in ODFln:
530	l,m,n = eval(SHname[1:])
531	Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
532	dNames += [str(pId)+'::'+SHname]
533	deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
534	for dName,drv in zip(dNames,deriv):
535	try:
536	ind = varyList.index(dName)
537	pDerv[ind][ip] += drv
538	except ValueError:
539	pass
540
541	lasthkl = np.array([0,0,0])
542	for ip,pName in enumerate(pNames):
543	deriv = []
544	dNames = []
545	pnames = pName.split(':')
546	if 'SH-' in pName and pId == int(pnames[0]):
547	hId = int(pnames[1])
548	phfx = '%d:%d:'%(pId,hId)
549	psi = float(pnames[4])
550	HKLs = calcControls[phfx+'SHhkl']
551	SHnames = calcControls[phfx+'SHnames']
552	SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
553	hkl = np.array(HKLs[int(pnames[3])])
554	if np.any(lasthkl-hkl):
555	PH = np.array(hkl)
556	phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
557	SH3Coef = {}
558	for item in SHcof:
559	L,N = eval(item.strip('C'))
560	SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]
561	ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
562	lasthkl = copy.copy(hkl)
563	for SHname in SHnames:
564	l,n = eval(SHname[1:])
565	SH3name = 'C%d,0,%d'%(l,n)
566	Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
567	dNames += [phfx+SHname]
568	deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
569	for dName,drv in zip(dNames,deriv):
570	try:
571	ind = varyList.index(dName)
572	pDerv[ind][ip] += drv
573	except ValueError:
574	pass
575	return pDerv
576
577	################################################################################
578	##### Function & derivative calculations
579	################################################################################
580
581	def GetAtomFXU(pfx,calcControls,parmDict):
582	'Needs a doc string'
583	Natoms = calcControls['Natoms'][pfx]
584	Tdata = Natoms*[' ',]
585	Mdata = np.zeros(Natoms)
586	IAdata = Natoms*[' ',]
587	Fdata = np.zeros(Natoms)
588	FFdata = []
589	BLdata = []
590	Xdata = np.zeros((3,Natoms))
591	dXdata = np.zeros((3,Natoms))
592	Uisodata = np.zeros(Natoms)
593	Uijdata = np.zeros((6,Natoms))
594	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
595	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
596	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
597	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
598	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
599	for iatm in range(Natoms):
600	for key in keys:
601	parm = pfx+key+str(iatm)
602	if parm in parmDict:
603	keys[key][iatm] = parmDict[parm]
604	Fdata = np.where(Fdata,Fdata,1.e-8) #avoid divide by zero in derivative calc.
605	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
606
607	def GetAtomSSFXU(pfx,calcControls,parmDict):
608	'Needs a doc string'
609	Natoms = calcControls['Natoms'][pfx]
610	maxSSwave = calcControls['maxSSwave'][pfx]
611	Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
612	'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
613	XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
614	FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
615	USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
616	MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
617	waveTypes = []
618	keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
619	'Tzero:':XSSdata[0],'Xslope:':XSSdata[1],'Yslope:':XSSdata[2],'Zslope:':XSSdata[3],
620	'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
621	'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
622	'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
623	'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
624	for iatm in range(Natoms):
625	waveTypes.append(parmDict[pfx+'waveType:'+str(iatm)])
626	for key in keys:
627	for m in range(Nwave[key[0]]):
628	parm = pfx+key+str(iatm)+':%d'%(m)
629	if parm in parmDict:
630	keys[key][m][iatm] = parmDict[parm]
631	return np.array(waveTypes),FSSdata,XSSdata,USSdata,MSSdata
632
633	def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):
634
635	Natoms = calcControls['Natoms'][pfx]
636	maxSSwave = calcControls['maxSSwave'][pfx]
637	Smult = np.zeros((Natoms,len(SGOps)))
638	TauT = np.zeros((Natoms,len(SGOps)))
639	for ix,xyz in enumerate(XData.T):
640	for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):
641	sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)
642	Smult[ix][isym] = sdet
643	TauT[ix][isym] = tauT
644	return Smult,TauT
645
646	def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
647	''' Compute structure factors for all h,k,l for phase
648	puts the result, F^2, in each ref[8] in refList
649	input:
650
651	:param dict refDict: where
652	'RefList' list where each ref = h,k,l,m,d,...
653	'FF' dict of form factors - filed in below
654	:param np.array G: reciprocal metric tensor
655	:param str pfx: phase id string
656	:param dict SGData: space group info. dictionary output from SpcGroup
657	:param dict calcControls:
658	:param dict ParmDict:
659
660	'''
661	phfx = pfx.split(':')[0]+hfx
662	ast = np.sqrt(np.diag(G))
663	Mast = twopisq*np.multiply.outer(ast,ast)
664	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
665	SGT = np.array([ops[1] for ops in SGData['SGOps']])
666	FFtables = calcControls['FFtables']
667	BLtables = calcControls['BLtables']
668	Flack = 1.0
669	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
670	Flack = 1.-2.*parmDict[phfx+'Flack']
671	TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
672	TwDict = refDict.get('TwDict',{})
673	if 'S' in calcControls[hfx+'histType']:
674	NTL = calcControls[phfx+'NTL']
675	NM = calcControls[phfx+'TwinNMN']+1
676	TwinLaw = calcControls[phfx+'TwinLaw']
677	TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
678	TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
679	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
680	FF = np.zeros(len(Tdata))
681	if 'NC' in calcControls[hfx+'histType']:
682	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
683	elif 'X' in calcControls[hfx+'histType']:
684	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
685	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
686	Uij = np.array(G2lat.U6toUij(Uijdata))
687	bij = Mast*Uij.T
688	blkSize = 100 #no. of reflections in a block
689	nRef = refDict['RefList'].shape[0]
690	if not len(refDict['FF']): #no form factors - 1st time thru StructureFactor
691	if 'N' in calcControls[hfx+'histType']:
692	dat = G2el.getBLvalues(BLtables)
693	refDict['FF']['El'] = dat.keys()
694	refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
695	else: #'X'
696	dat = G2el.getFFvalues(FFtables,0.)
697	refDict['FF']['El'] = dat.keys()
698	refDict['FF']['FF'] = np.ones((nRef,len(dat)))
699	for iref,ref in enumerate(refDict['RefList']):
700	SQ = 1./(2.ref[4])*2
701	dat = G2el.getFFvalues(FFtables,SQ)
702	refDict['FF']['FF'][iref] *= dat.values()
703	#reflection processing begins here - big arrays!
704	iBeg = 0
705	while iBeg < nRef:
706	iFin = min(iBeg+blkSize,nRef)
707	refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems)
708	H = refl.T[:3] #array(blkSize,3)
709	H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,nTwins,3) or (blkSize,3)
710	TwMask = np.any(H,axis=-1)
711	if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
712	for ir in range(blkSize):
713	iref = ir+iBeg
714	if iref in TwDict:
715	for i in TwDict[iref]:
716	for n in range(NTL):
717	H[ir][i+nNM] = np.inner(TwinLaw[nNM],np.array(TwDict[iref][i])TwinInv[i+nNM])
718	TwMask = np.any(H,axis=-1)
719	SQ = 1./(2.refl.T[4])*2 #array(blkSize)
720	SQfactor = 4.0SQtwopisq #ditto prev.
721	if 'T' in calcControls[hfx+'histType']:
722	if 'P' in calcControls[hfx+'histType']:
723	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
724	else:
725	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
726	FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
727	FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
728	Bab = np.repeat(parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor),len(SGT)*len(TwinLaw))
729	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
730	FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
731	Uniq = np.inner(H,SGMT)
732	Phi = np.inner(H,SGT)
733	phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
734	sinp = np.sin(phase)
735	cosp = np.cos(phase)
736	biso = -SQfactor*Uisodata[:,np.newaxis]
737	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
738	HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
739	Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
740	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/len(SGMT)
741	if 'T' in calcControls[hfx+'histType']:
742	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-np.reshape(FlackFPP,sinp.shape)sinp*Tcorr])
743	fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinpTcorr,np.reshape(FlackFPP,cosp.shape)cosp*Tcorr])
744	else:
745	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-FlackFPPsinp*Tcorr])
746	fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinpTcorr,FlackFPPcosp*Tcorr])
747	fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl
748	fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl
749	if SGData['SGInv']: #centrosymmetric; B=0
750	fbs[0] *= 0.
751	if 'P' in calcControls[hfx+'histType']:
752	refl.T[9] = np.sum(fas2,axis=0)+np.sum(fbs2,axis=0)
753	refl.T[10] = atan2d(fbs[0],fas[0]) #ignore f' & f"
754	else:
755	if len(TwinLaw) > 1:
756	refl.T[9] = np.sum(fas[:,:,0],axis=0)2+np.sum(fbs[:,:,0],axis=0)2 #FcT from primary twin element
757	refl.T[7] = np.sum(TwinFrnp.sum(TwMask[np.newaxis,:,:]fas,axis=0)**2,axis=-1)+ \
758	np.sum(TwinFrnp.sum(TwMask[np.newaxis,:,:]fbs,axis=0)**2,axis=-1) #Fc sum over twins
759	refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f"
760	else:
761	refl.T[9] = np.sum(fas,axis=0)2+np.sum(fbs,axis=0)2
762	refl.T[7] = np.copy(refl.T[9])
763	refl.T[10] = atan2d(fbs[0],fas[0]) #ignore f' & f"
764	# refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
765	iBeg += blkSize
766
767	def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
768	'Needs a doc string'
769	phfx = pfx.split(':')[0]+hfx
770	ast = np.sqrt(np.diag(G))
771	Mast = twopisq*np.multiply.outer(ast,ast)
772	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
773	SGT = np.array([ops[1] for ops in SGData['SGOps']])
774	FFtables = calcControls['FFtables']
775	BLtables = calcControls['BLtables']
776	TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
777	TwDict = refDict.get('TwDict',{})
778	if 'S' in calcControls[hfx+'histType']:
779	NTL = calcControls[phfx+'NTL']
780	NM = calcControls[phfx+'TwinNMN']+1
781	TwinLaw = calcControls[phfx+'TwinLaw']
782	TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
783	TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
784	nTwin = len(TwinLaw)
785	nRef = len(refDict['RefList'])
786	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
787	mSize = len(Mdata)
788	FF = np.zeros(len(Tdata))
789	if 'NC' in calcControls[hfx+'histType']:
790	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
791	elif 'X' in calcControls[hfx+'histType']:
792	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
793	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
794	Uij = np.array(G2lat.U6toUij(Uijdata))
795	bij = Mast*Uij.T
796	dFdvDict = {}
797	if nTwin > 1:
798	dFdfr = np.zeros((nRef,nTwin,mSize))
799	dFdx = np.zeros((nRef,nTwin,mSize,3))
800	dFdui = np.zeros((nRef,nTwin,mSize))
801	dFdua = np.zeros((nRef,nTwin,mSize,6))
802	dFdbab = np.zeros((nRef,nTwin,2))
803	dFdfl = np.zeros((nRef,nTwin))
804	dFdtw = np.zeros((nRef,nTwin))
805	else:
806	dFdfr = np.zeros((nRef,mSize))
807	dFdx = np.zeros((nRef,mSize,3))
808	dFdui = np.zeros((nRef,mSize))
809	dFdua = np.zeros((nRef,mSize,6))
810	dFdbab = np.zeros((nRef,2))
811	dFdfl = np.zeros((nRef))
812	dFdtw = np.zeros((nRef))
813	Flack = 1.0
814	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
815	Flack = 1.-2.*parmDict[phfx+'Flack']
816	time0 = time.time()
817	nref = len(refDict['RefList'])/100
818	for iref,refl in enumerate(refDict['RefList']):
819	if 'T' in calcControls[hfx+'histType']:
820	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
821	H = np.array(refl[:3])
822	H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3)
823	TwMask = np.any(H,axis=-1)
824	if TwinLaw.shape[0] > 1 and TwDict:
825	if iref in TwDict:
826	for i in TwDict[iref]:
827	for n in range(NTL):
828	H[i+nNM] = np.inner(TwinLaw[nNM],np.array(TwDict[iref][i])TwinInv[i+nNM])
829	TwMask = np.any(H,axis=-1)
830	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
831	SQfactor = 8.0SQnp.pi**2
832	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
833	Bab = parmDict[phfx+'BabA']*dBabdA
834	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
835	FF = refDict['FF']['FF'][iref].T[Tindx].T
836	Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3)
837	Phi = np.inner(H,SGT)
838	phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
839	sinp = np.sin(phase)
840	cosp = np.cos(phase)
841	occ = Mdata*Fdata/len(SGT)
842	biso = -SQfactor*Uisodata[:,np.newaxis]
843	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
844	HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
845	Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
846	Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
847	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
848	Tcorr = (np.reshape(Tiso,Tuij.shape)Tuij).Tocc
849	fot = (FF+FP-Bab)*Tcorr
850	fotp = FPP*Tcorr
851	fa = np.array([((FF+FP).T-Bab).TcospTcorr,-FlackFPPsinp*Tcorr])
852	fb = np.array([((FF+FP).T-Bab).TsinpTcorr,FlackFPPcosp*Tcorr])
853	fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins)
854	fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl
855	fax = np.array([-fotsinp,-fotpcosp]) #positions array(2,ntwi,nEqv,nAtoms)
856	fbx = np.array([fotcosp,-fotpsinp])
857	#sum below is over Uniq
858	dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
859	dfadba = np.sum(-cosp*Tcorr[:,np.newaxis],axis=1)
860	dfadui = np.sum(-SQfactor*fa,axis=-2)
861	if nTwin > 1:
862	dfadx = np.array([np.sum(twopiUniq[it]np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
863	dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
864	# array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
865	else:
866	dfadx = np.sum(twopiUniqnp.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
867	dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
868	# array(2,nAtom,3) & array(2,nAtom,6)
869	if not SGData['SGInv']:
870	dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem
871	dfbdba = np.sum(-sinp*Tcorr[:,np.newaxis],axis=1)
872	dfbdui = np.sum(-SQfactor*fb,axis=-2)
873	if len(TwinLaw) > 1:
874	dfbdx = np.array([np.sum(twopiUniq[it]np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
875	dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
876	else:
877	dfadfl = np.sum(-FPPTcorrsinp)
878	dfbdfl = np.sum(FPPTcorrcosp)
879	dfbdx = np.sum(twopiUniqnp.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)
880	dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)
881	else:
882	dfbdfr = np.zeros_like(dfadfr)
883	dfbdx = np.zeros_like(dfadx)
884	dfbdui = np.zeros_like(dfadui)
885	dfbdua = np.zeros_like(dfadua)
886	dfbdba = np.zeros_like(dfadba)
887	dfadfl = 0.0
888	dfbdfl = 0.0
889	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!
890	if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
891	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)+ \
892	2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
893	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])+ \
894	2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
895	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])+ \
896	2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
897	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])+ \
898	2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
899	else:
900	SA = fas[0]+fas[1]
901	SB = fbs[0]+fbs[1]
902	if nTwin > 1:
903	dFdfr[iref] = [2.TwMask[it](SA[it]dfadfr[0][it]+SA[it]dfadfr[1][it]+SB[it]dfbdfr[0][it]+SB[it]dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
904	dFdx[iref] = [2.TwMask[it](SA[it]dfadx[it][0]+SA[it]dfadx[it][1]+SB[it]dfbdx[it][0]+SB[it]dfbdx[it][1]) for it in range(nTwin)]
905	dFdui[iref] = [2.TwMask[it](SA[it]dfadui[it][0]+SA[it]dfadui[it][1]+SB[it]dfbdui[it][0]+SB[it]dfbdui[it][1]) for it in range(nTwin)]
906	dFdua[iref] = [2.TwMask[it](SA[it]dfadua[it][0]+SA[it]dfadua[it][1]+SB[it]dfbdua[it][0]+SB[it]dfbdua[it][1]) for it in range(nTwin)]
907	dFdtw[iref] = np.sum(TwMaskfas,axis=0)2+np.sum(TwMaskfbs,axis=0)**2
908
909	else: #these are good for no twin single crystals
910	dFdfr[iref] = (2.SA(dfadfr[0]+dfadfr[1])+2.SB(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)
911	dFdx[iref] = 2.SA(dfadx[0]+dfadx[1])+2.SB(dfbdx[0]+dfbdx[1])
912	dFdui[iref] = 2.SA(dfadui[0]+dfadui[1])+2.SB(dfbdui[0]+dfbdui[1])
913	dFdua[iref] = 2.SA(dfadua[0]+dfadua[1])+2.SB(dfbdua[0]+dfbdua[1])
914	dFdfl[iref] = -SAdfadfl-SBdfbdfl #array(nRef,)
915	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T+ \
916	2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
917	# GSASIIpath.IPyBreak()
918	if not iref%100 :
919	print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
920	#loop over atoms - each dict entry is list of derivatives for all the reflections
921	if nTwin > 1:
922	for i in range(len(Mdata)): #these all OK?
923	dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0)
924	dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)
925	dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)
926	dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)
927	dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)
928	dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)
929	dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)
930	dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)
931	dFdvDict[pfx+'AU12:'+str(i)] = 2.np.sum(dFdua.T[3][i]TwinFr[:,np.newaxis],axis=0)
932	dFdvDict[pfx+'AU13:'+str(i)] = 2.np.sum(dFdua.T[4][i]TwinFr[:,np.newaxis],axis=0)
933	dFdvDict[pfx+'AU23:'+str(i)] = 2.np.sum(dFdua.T[5][i]TwinFr[:,np.newaxis],axis=0)
934	else:
935	for i in range(len(Mdata)):
936	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
937	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
938	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
939	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
940	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
941	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
942	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
943	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
944	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
945	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
946	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
947	dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
948	dFdvDict[phfx+'BabA'] = dFdbab.T[0]
949	dFdvDict[phfx+'BabU'] = dFdbab.T[1]
950	if nTwin > 1:
951	for i in range(nTwin):
952	dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
953	return dFdvDict
954
955	def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
956	'''
957	Compute super structure factors for all h,k,l,m for phase
958	puts the result, F^2, in each ref[8+im] in refList
959	input:
960
961	:param dict refDict: where
962	'RefList' list where each ref = h,k,l,m,d,...
963	'FF' dict of form factors - filed in below
964	:param np.array G: reciprocal metric tensor
965	:param str pfx: phase id string
966	:param dict SGData: space group info. dictionary output from SpcGroup
967	:param dict calcControls:
968	:param dict ParmDict:
969
970	'''
971	nxs = np.newaxis
972	phfx = pfx.split(':')[0]+hfx
973	ast = np.sqrt(np.diag(G))
974	Mast = twopisq*np.multiply.outer(ast,ast)
975	SGInv = SGData['SGInv']
976	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
977	SGT = np.array([ops[1] for ops in SGData['SGOps']])
978	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
979	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
980	# eps = SSGMT[:,3,3]
981	FFtables = calcControls['FFtables']
982	BLtables = calcControls['BLtables']
983	Flack = 1.0
984	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
985	Flack = 1.-2.*parmDict[phfx+'Flack']
986	TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],]) #4D?
987	TwDict = refDict.get('TwDict',{})
988	if 'S' in calcControls[hfx+'histType']:
989	NTL = calcControls[phfx+'NTL']
990	NM = calcControls[phfx+'TwinNMN']+1
991	TwinLaw = calcControls[phfx+'TwinLaw'] #this'll have to be 4D also...
992	TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
993	TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
994	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
995	waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
996	# Smult,TauT = GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata)
997	# GSASIIpath.IPyBreak()
998	# if SGInv:
999	# TauT = np.hstack((TauT,-TauT))
1000	# eps = np.concatenate((eps,-eps))
1001	modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1002	FF = np.zeros(len(Tdata))
1003	if 'NC' in calcControls[hfx+'histType']:
1004	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1005	else:
1006	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1007	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1008	Uij = np.array(G2lat.U6toUij(Uijdata)).T
1009	bij = Mast*Uij
1010	blkSize = 100 #no. of reflections in a block
1011	nRef = refDict['RefList'].shape[0]
1012	if not len(refDict['FF']):
1013	if 'N' in calcControls[hfx+'histType']:
1014	dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's
1015	refDict['FF']['El'] = dat.keys()
1016	refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
1017	else:
1018	dat = G2el.getFFvalues(FFtables,0.)
1019	refDict['FF']['El'] = dat.keys()
1020	refDict['FF']['FF'] = np.ones((nRef,len(dat)))
1021	for iref,ref in enumerate(refDict['RefList']):
1022	SQ = 1./(2.ref[5])*2
1023	dat = G2el.getFFvalues(FFtables,SQ)
1024	refDict['FF']['FF'][iref] *= dat.values()
1025	time0 = time.time()
1026	#reflection processing begins here - big arrays!
1027	iBeg = 0
1028	while iBeg < nRef:
1029	iFin = min(iBeg+blkSize,nRef)
1030	refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems)
1031	H = refl.T[:4] #array(blkSize,4)
1032	HP = H[:3]+modQ[:,nxs]*H[3:] #projected hklm to hkl
1033	# H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,nTwins,4) or (blkSize,4)
1034	# TwMask = np.any(H,axis=-1)
1035	# if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1036	# for ir in range(blkSize):
1037	# iref = ir+iBeg
1038	# if iref in TwDict:
1039	# for i in TwDict[iref]:
1040	# for n in range(NTL):
1041	# H[ir][i+nNM] = np.inner(TwinLaw[nNM],np.array(TwDict[iref][i])TwinInv[i+nNM])
1042	# TwMask = np.any(H,axis=-1)
1043	SQ = 1./(2.refl.T[5])*2 #array(blkSize)
1044	SQfactor = 4.0SQtwopisq #ditto prev.
1045	Uniq = np.inner(H.T,SSGMT)
1046	UniqP = np.inner(HP.T,SGMT)
1047	Phi = np.inner(H.T,SSGT)
1048	# PhiP = np.inner(HP.T,SGT)
1049	if SGInv: #if centro - expand HKL sets
1050	Uniq = np.hstack((Uniq,-Uniq))
1051	Phi = np.hstack((Phi,-Phi))
1052	UniqP = np.hstack((UniqP,-UniqP))
1053	# PhiP = np.hstack((PhiP,-PhiP))
1054	if 'T' in calcControls[hfx+'histType']:
1055	if 'P' in calcControls[hfx+'histType']:
1056	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1057	else:
1058	FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1059	FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1060	FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1061	Bab = np.repeat(parmDict[phfx+'BabA']np.exp(-parmDict[phfx+'BabU']SQfactor),Uniq.shape[1]*len(TwinLaw))
1062	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1063	FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1064	phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1065	sinp = np.sin(phase)
1066	cosp = np.cos(phase)
1067	biso = -SQfactor*Uisodata[:,nxs]
1068	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1069	HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1) #use hklt proj to hkl
1070	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1071	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/Uniq.shape[1] #refBlk x ops x atoms
1072	if 'T' in calcControls[hfx+'histType']:
1073	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-np.reshape(FlackFPP,sinp.shape)sinp*Tcorr])
1074	fb = np.array([np.reshape(FlackFPP,cosp.shape)cospTcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinp*Tcorr])
1075	else:
1076	fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)cospTcorr,-FlackFPPsinp*Tcorr])
1077	fb = np.array([FlackFPPcospTcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)sinp*Tcorr])
1078	GfpuA = G2mth.Modulation(waveTypes,Uniq,UniqP,FSSdata,XSSdata,USSdata,Mast) #2 x refBlk x sym X atoms
1079	fag = faGfpuA[0]-fbGfpuA[1] #real; 2 x refBlk x sym x atoms
1080	fbg = fbGfpuA[0]+faGfpuA[1]
1081	fas = np.sum(np.sum(fag,axis=-1),axis=-1) #2 x refBlk; sum sym & atoms
1082	fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1083	# GSASIIpath.IPyBreak()
1084	if 'P' in calcControls[hfx+'histType']:
1085	refl.T[10] = np.sum(fas2,axis=0)+np.sum(fbs2,axis=0)
1086	refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f"
1087	else:
1088	if len(TwinLaw) > 1:
1089	refl.T[10] = np.sum(fas[:,:,0],axis=0)2+np.sum(fbs[:,:,0],axis=0)2 #FcT from primary twin element
1090	refl.T[8] = np.sum(TwinFrnp.sum(TwMask[np.newaxis,:,:]fas,axis=0)**2,axis=-1)+ \
1091	np.sum(TwinFrnp.sum(TwMask[np.newaxis,:,:]fbs,axis=0)**2,axis=-1) #Fc sum over twins
1092	refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f"
1093	else:
1094	refl.T[10] = np.sum(fas,axis=0)2+np.sum(fbs,axis=0)2 #square of sums
1095	refl.T[8] = np.copy(refl.T[10])
1096	refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f"
1097	iBeg += blkSize
1098	print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
1099
1100	def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1101	'Needs a doc string'
1102	nxs = np.newaxis
1103	phfx = pfx.split(':')[0]+hfx
1104	ast = np.sqrt(np.diag(G))
1105	Mast = twopisq*np.multiply.outer(ast,ast)
1106	SGInv = SGData['SGInv']
1107	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1108	SGT = np.array([ops[1] for ops in SGData['SGOps']])
1109	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1110	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1111	FFtables = calcControls['FFtables']
1112	BLtables = calcControls['BLtables']
1113	TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
1114	TwDict = refDict.get('TwDict',{})
1115	if 'S' in calcControls[hfx+'histType']:
1116	NTL = calcControls[phfx+'NTL']
1117	NM = calcControls[phfx+'TwinNMN']+1
1118	TwinLaw = calcControls[phfx+'TwinLaw']
1119	TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1120	TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1121	nTwin = len(TwinLaw)
1122	nRef = len(refDict['RefList'])
1123	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1124	mSize = len(Mdata) #no. atoms
1125	waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1126	# Smult,TauT = GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata)
1127	# if SGInv:
1128	# TauT = np.hstack((TauT,-TauT))
1129	modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1130	FF = np.zeros(len(Tdata))
1131	if 'NC' in calcControls[hfx+'histType']:
1132	FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1133	elif 'X' in calcControls[hfx+'histType']:
1134	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1135	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1136	Uij = np.array(G2lat.U6toUij(Uijdata)).T
1137	bij = Mast*Uij
1138	if not len(refDict['FF']):
1139	if 'N' in calcControls[hfx+'histType']:
1140	dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's
1141	else:
1142	dat = G2el.getFFvalues(FFtables,0.)
1143	refDict['FF']['El'] = dat.keys()
1144	refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1145	dFdvDict = {}
1146	if nTwin > 1:
1147	dFdfr = np.zeros((nRef,nTwin,mSize))
1148	dFdx = np.zeros((nRef,nTwin,mSize,3))
1149	dFdui = np.zeros((nRef,nTwin,mSize))
1150	dFdua = np.zeros((nRef,nTwin,mSize,6))
1151	dFdbab = np.zeros((nRef,nTwin,2))
1152	dFdfl = np.zeros((nRef,nTwin))
1153	dFdtw = np.zeros((nRef,nTwin))
1154	dFdGfA = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
1155	dFdGfB = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
1156	dFdGxA = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
1157	dFdGxB = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
1158	dFdGuA = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
1159	dFdGuB = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
1160	else:
1161	dFdfr = np.zeros((nRef,mSize))
1162	dFdx = np.zeros((nRef,mSize,3))
1163	dFdui = np.zeros((nRef,mSize))
1164	dFdua = np.zeros((nRef,mSize,6))
1165	dFdbab = np.zeros((nRef,2))
1166	dFdfl = np.zeros((nRef))
1167	dFdtw = np.zeros((nRef))
1168	dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1169	dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1170	dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1171	Flack = 1.0
1172	if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1173	Flack = 1.-2.*parmDict[phfx+'Flack']
1174	time0 = time.time()
1175	nRef = len(refDict['RefList'])/100
1176	for iref,refl in enumerate(refDict['RefList']):
1177	if 'T' in calcControls[hfx+'histType']:
1178	FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1179	H = np.array(refl[:4])
1180	HP = H[:3]+modQ*H[3:] #projected hklm to hkl
1181	# H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(4,nTwins) or (4)
1182	# TwMask = np.any(H,axis=-1)
1183	# if TwinLaw.shape[0] > 1 and TwDict:
1184	# if iref in TwDict:
1185	# for i in TwDict[iref]:
1186	# for n in range(NTL):
1187	# H[i+nNM] = np.inner(TwinLaw[nNM],np.array(TwDict[iref][i])TwinInv[i+nNM])
1188	# TwMask = np.any(H,axis=-1)
1189	SQ = 1./(2.refl[4+im])2 # or (sin(theta)/lambda)*2
1190	SQfactor = 8.0SQnp.pi**2
1191	dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1192	Bab = parmDict[phfx+'BabA']*dBabdA
1193	Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1194	FF = refDict['FF']['FF'][iref].T[Tindx]
1195	Uniq = np.inner(H,SSGMT)
1196	Phi = np.inner(H,SSGT)
1197	UniqP = np.inner(HP,SGMT)
1198	if SGInv: #if centro - expand HKL sets
1199	Uniq = np.vstack((Uniq,-Uniq))
1200	Phi = np.hstack((Phi,-Phi))
1201	UniqP = np.vstack((UniqP,-UniqP))
1202	phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1203	sinp = np.sin(phase)
1204	cosp = np.cos(phase)
1205	occ = Mdata*Fdata/Uniq.shape[0]
1206	biso = -SQfactor*Uisodata[:,nxs]
1207	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T #ops x atoms
1208	HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1) #ops x atoms
1209	Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1210	Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
1211	Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms
1212	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/Uniq.shape[0] #ops x atoms
1213	fot = (FF+FP-Bab)*Tcorr #ops x atoms
1214	fotp = FPP*Tcorr #ops x atoms
1215	GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,UniqP,Hij,FSSdata,XSSdata,USSdata,Mast)
1216	# GfpuA is 2 x ops x atoms
1217	# derivs are: ops x atoms x waves x 2,6,or 12 parms as [real,imag] parts
1218	fa = np.array([((FF+FP).T-Bab).TcospTcorr,-FlackFPPsinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
1219	fb = np.array([((FF+FP).T-Bab).TsinpTcorr,FlackFPPcosp*Tcorr]) #or array(2,nEqv,nAtoms)
1220	fag = faGfpuA[0]-fbGfpuA[1]
1221	fbg = fbGfpuA[0]+faGfpuA[1]
1222
1223	fas = np.sum(np.sum(fag,axis=1),axis=1) # 2 x twin
1224	fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1225	fax = np.array([-fotsinp,-fotpcosp]) #positions; 2 x twin x ops x atoms
1226	fbx = np.array([fotcosp,-fotpsinp])
1227	fax = faxGfpuA[0]-fbxGfpuA[1]
1228	fbx = fbxGfpuA[0]+faxGfpuA[1]
1229	#sum below is over Uniq
1230	dfadfr = np.sum(fag/occ,axis=1) #Fdata != 0 ever avoids /0. problem
1231	dfbdfr = np.sum(fbg/occ,axis=1) #Fdata != 0 avoids /0. problem
1232	dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1233	dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1234	dfadui = np.sum(-SQfactor*fag,axis=1)
1235	dfbdui = np.sum(-SQfactor*fbg,axis=1)
1236	if nTwin > 1:
1237	dfadx = np.array([np.sum(twopiUniq[it,:,:3]np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1238	dfbdx = np.array([np.sum(twopiUniq[it,:,:3]np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1239	dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1240	dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1241	# array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
1242	dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]dGdf[1][nxs,nxs,:,:,:,:],axis=1)
1243	dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]dGdf[1][nxs,nxs,:,:,:,:],axis=1)
1244	dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]dGdx[1][nxs,nxs,:,:,:,:],axis=1)
1245	dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]dGdx[1][nxs,nxs,:,:,:,:],axis=1)
1246	dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]dGdu[1][nxs,nxs,:,:,:,:],axis=1)
1247	dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]dGdu[1][nxs,nxs,:,:,:,:],axis=1)
1248	else:
1249	dfadx = np.sum(twopiUniq[:,:3]np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2) #2 x nAtom x 3xyz; sum nOps
1250	dfbdx = np.sum(twopiUniq[:,:3]np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)
1251	dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2) #2 x nAtom x 6Uij; sum nOps
1252	dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2) #these are correct also for twins above
1253	# array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1254	dfadGf = np.sum(fa[:,:,:,nxs,nxs]dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]dGdf[1][nxs,:,:,:,:],axis=1)
1255	dfbdGf = np.sum(fb[:,:,:,nxs,nxs]dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]dGdf[1][nxs,:,:,:,:],axis=1)
1256	dfadGx = np.sum(fa[:,:,:,nxs,nxs]dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]dGdx[1][nxs,:,:,:,:],axis=1)
1257	dfbdGx = np.sum(fb[:,:,:,nxs,nxs]dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]dGdx[1][nxs,:,:,:,:],axis=1)
1258	dfadGu = np.sum(fa[:,:,:,nxs,nxs]dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]dGdu[1][nxs,:,:,:,:],axis=1)
1259	dfbdGu = np.sum(fb[:,:,:,nxs,nxs]dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]dGdu[1][nxs,:,:,:,:],axis=1)
1260	if not SGData['SGInv'] and len(TwinLaw) == 1: #Flack derivative
1261	dfadfl = np.sum(-FPPTcorrsinp)
1262	dfbdfl = np.sum(FPPTcorrcosp)
1263	else:
1264	dfadfl = 1.0
1265	dfbdfl = 1.0
1266	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!
1267	if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
1268	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)+ \
1269	2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
1270	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])+ \
1271	2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
1272	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])+ \
1273	2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
1274	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])+ \
1275	2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
1276	dFdGx[iref] = 2.(fas[0]dfadGx[0]+fas[1]*dfadGx[1])+ \
1277	2.(fbs[0]dfbdGx[0]+fbs[1]*dfbdGx[1])
1278	dFdGu[iref] = 2.(fas[0]dfadGu[0]+fas[1]*dfadGu[1])+ \
1279	2.(fbs[0]dfbdGu[0]+fbs[1]*dfbdGu[1])
1280	else:
1281	SA = fas[0]+fas[1] #float = A+A' (might be array[nTwin])
1282	SB = fbs[0]+fbs[1] #float = B+B' (might be array[nTwin])
1283	if nTwin > 1:
1284	dFdfr[iref] = [2.TwMask[it](SA[it]dfadfr[0][it]+SA[it]dfadfr[1][it]+SB[it]dfbdfr[0][it]+SB[it]dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
1285	dFdx[iref] = [2.TwMask[it](SA[it]dfadx[it][0]+SA[it]dfadx[it][1]+SB[it]dfbdx[it][0]+SB[it]dfbdx[it][1]) for it in range(nTwin)]
1286	dFdui[iref] = [2.TwMask[it](SA[it]dfadui[it][0]+SA[it]dfadui[it][1]+SB[it]dfbdui[it][0]+SB[it]dfbdui[it][1]) for it in range(nTwin)]
1287	dFdua[iref] = [2.TwMask[it](SA[it]dfadua[it][0]+SA[it]dfadua[it][1]+SB[it]dfbdua[it][0]+SB[it]dfbdua[it][1]) for it in range(nTwin)]
1288	dFdtw[iref] = np.sum(TwMaskfas,axis=0)2+np.sum(TwMaskfbs,axis=0)**2
1289
1290	dFdGf[iref] = [2.TwMask[it](SA[it]dfadGf[1]+SB[it]dfbdGf[1]) for it in range(nTwin)]
1291	dFdGx[iref] = [2.TwMask[it](SA[it]dfadGx[1]+SB[it]dfbdGx[1]) for it in range(nTwin)]
1292	dFdGu[iref] = [2.TwMask[it](SA[it]dfadGu[1]+SB[it]dfbdGu[1]) for it in range(nTwin)]
1293	else: #these are good for no twin single crystals
1294	dFdfr[iref] = 2.(SAdfadfr[0]+SAdfadfr[1]+SBdfbdfr[0]+SBdfbdfr[1])Mdata/len(Uniq) #array(nRef,nAtom)
1295	dFdx[iref] = 2.(SAdfadx[0]+SAdfadx[1]+SBdfbdx[0]+SB*dfbdx[1]) #array(nRef,nAtom,3)
1296	dFdui[iref] = 2.(SAdfadui[0]+SAdfadui[1]+SBdfbdui[0]+SB*dfbdui[1]) #array(nRef,nAtom)
1297	dFdua[iref] = 2.(SAdfadua[0]+SAdfadua[1]+SBdfbdua[0]+SB*dfbdua[1]) #array(nRef,nAtom,6)
1298	dFdfl[iref] = -SAdfadfl-SBdfbdfl #array(nRef,)
1299
1300	dFdGf[iref] = 2.(SAdfadGf[0]+SB*dfbdGf[1]) #array(nRef,natom,nwave,2)
1301	dFdGx[iref] = 2.(SAdfadGx[0]+SB*dfbdGx[1]) #array(nRef,natom,nwave,6)
1302	dFdGu[iref] = 2.(SAdfadGu[0]+SB*dfbdGu[1]) #array(nRef,natom,nwave,12)
1303	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T+ \
1304	2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
1305	#loop over atoms - each dict entry is list of derivatives for all the reflections
1306	# GSASIIpath.IPyBreak()
1307	if not iref%100 :
1308	print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
1309	for i in range(len(Mdata)): #loop over atoms
1310	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1311	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1312	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1313	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1314	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1315	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1316	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1317	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1318	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1319	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1320	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1321	for j in range(FSSdata.shape[1]): #loop over waves Fzero & Fwid?
1322	dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1323	dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1324	for j in range(XSSdata.shape[1]): #loop over waves 'Tzero','Xslope','Yslope','Zslope'?
1325	dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j)] = dFdGx.T[0][j][i]
1326	dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j)] = dFdGx.T[1][j][i]
1327	dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j)] = dFdGx.T[2][j][i]
1328	dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j)] = dFdGx.T[3][j][i]
1329	dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j)] = dFdGx.T[4][j][i]
1330	dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j)] = dFdGx.T[5][j][i]
1331	for j in range(USSdata.shape[1]): #loop over waves
1332	dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1333	dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1334	dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1335	dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
1336	dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
1337	dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
1338	dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1339	dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1340	dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1341	dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
1342	dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
1343	dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
1344
1345	dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1346	dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1347	return dFdvDict
1348
1349	def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1350	''' Single crystal extinction function; returns extinction & derivative
1351	'''
1352	extCor = 1.0
1353	dervDict = {}
1354	dervCor = 1.0
1355	if calcControls[phfx+'EType'] != 'None':
1356	SQ = 1/(4.ref[4+im]*2)
1357	if 'C' in parmDict[hfx+'Type']:
1358	cos2T = 1.0-2.SQparmDict[hfx+'Lam']**2 #cos(2theta)
1359	else: #'T'
1360	cos2T = 1.0-2.SQref[12+im]**2 #cos(2theta)
1361	if 'SXC' in parmDict[hfx+'Type']:
1362	AV = 7.9406e5/parmDict[pfx+'Vol']**2
1363	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1364	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
1365	PLZ = AVP12ref[9+im]parmDict[hfx+'Lam']*2
1366	elif 'SNT' in parmDict[hfx+'Type']:
1367	AV = 1.e7/parmDict[pfx+'Vol']**2
1368	PL = SQ
1369	PLZ = AVref[9+im]ref[12+im]**2
1370	elif 'SNC' in parmDict[hfx+'Type']:
1371	AV = 1.e7/parmDict[pfx+'Vol']**2
1372	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1373	PLZ = AVref[9+im]parmDict[hfx+'Lam']**2
1374
1375	if 'Primary' in calcControls[phfx+'EType']:
1376	PLZ *= 1.5
1377	else:
1378	if 'C' in parmDict[hfx+'Type']:
1379	PLZ *= calcControls[phfx+'Tbar']
1380	else: #'T'
1381	PLZ *= ref[13+im] #t-bar
1382	if 'Primary' in calcControls[phfx+'EType']:
1383	PLZ *= 1.5
1384	PSIG = parmDict[phfx+'Ep']
1385	elif 'I & II' in calcControls[phfx+'EType']:
1386	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
1387	elif 'Type II' in calcControls[phfx+'EType']:
1388	PSIG = parmDict[phfx+'Es']
1389	else: # 'Secondary Type I'
1390	PSIG = parmDict[phfx+'Eg']/PL
1391
1392	AG = 0.58+0.48cos2T+0.24cos2T**2
1393	AL = 0.025+0.285*cos2T
1394	BG = 0.02-0.025*cos2T
1395	BL = 0.15-0.2(0.75-cos2T)*2
1396	if cos2T < 0.:
1397	BL = -0.45*cos2T
1398	CG = 2.
1399	CL = 2.
1400	PF = PLZ*PSIG
1401
1402	if 'Gaussian' in calcControls[phfx+'EApprox']:
1403	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
1404	extCor = np.sqrt(PF4)
1405	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
1406	else:
1407	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
1408	extCor = np.sqrt(PF4)
1409	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
1410
1411	dervCor = (1.+PF)*PF3 #extinction corr for other derivatives
1412	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1413	dervDict[phfx+'Ep'] = -ref[7+im]PLZPF3
1414	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1415	dervDict[phfx+'Es'] = -ref[7+im]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
1416	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1417	dervDict[phfx+'Eg'] = -ref[7+im]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
1418
1419	return 1./extCor,dervDict,dervCor
1420
1421	def Dict2Values(parmdict, varylist):
1422	'''Use before call to leastsq to setup list of values for the parameters
1423	in parmdict, as selected by key in varylist'''
1424	return [parmdict[key] for key in varylist]
1425
1426	def Values2Dict(parmdict, varylist, values):
1427	''' Use after call to leastsq to update the parameter dictionary with
1428	values corresponding to keys in varylist'''
1429	parmdict.update(zip(varylist,values))
1430
1431	def GetNewCellParms(parmDict,varyList):
1432	'Needs a doc string'
1433	newCellDict = {}
1434	Anames = ['A'+str(i) for i in range(6)]
1435	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1436	for item in varyList:
1437	keys = item.split(':')
1438	if keys[2] in Ddict:
1439	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
1440	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
1441	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1442	return newCellDict # is e.g. {'0::D11':A0-D11}
1443
1444	def ApplyXYZshifts(parmDict,varyList):
1445	'''
1446	takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1447
1448	:param dict parmDict: parameter dictionary
1449	:param list varyList: list of variables (not used!)
1450	:returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1451
1452	'''
1453	newAtomDict = {}
1454	for item in parmDict:
1455	if 'dA' in item:
1456	parm = ''.join(item.split('d'))
1457	parmDict[parm] += parmDict[item]
1458	newAtomDict[item] = [parm,parmDict[parm]]
1459	return newAtomDict
1460
1461	def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1462	'Spherical harmonics texture'
1463	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1464	if 'T' in calcControls[hfx+'histType']:
1465	tth = parmDict[hfx+'2-theta']
1466	else:
1467	tth = refl[5+im]
1468	odfCor = 1.0
1469	H = refl[:3]
1470	cell = G2lat.Gmat2cell(g)
1471	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1472	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1473	phi,beta = G2lat.CrsAng(H,cell,SGData)
1474	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1475	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1476	for item in SHnames:
1477	L,M,N = eval(item.strip('C'))
1478	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1479	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1480	Lnorm = G2lat.Lnorm(L)
1481	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1482	return odfCor
1483
1484	def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1485	'Spherical harmonics texture derivatives'
1486	if 'T' in calcControls[hfx+'histType']:
1487	tth = parmDict[hfx+'2-theta']
1488	else:
1489	tth = refl[5+im]
1490	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1491	odfCor = 1.0
1492	dFdODF = {}
1493	dFdSA = [0,0,0]
1494	H = refl[:3]
1495	cell = G2lat.Gmat2cell(g)
1496	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1497	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1498	phi,beta = G2lat.CrsAng(H,cell,SGData)
1499	psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1500	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1501	for item in SHnames:
1502	L,M,N = eval(item.strip('C'))
1503	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1504	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1505	Lnorm = G2lat.Lnorm(L)
1506	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1507	dFdODF[pfx+item] = LnormKclKsl
1508	for i in range(3):
1509	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
1510	return odfCor,dFdODF,dFdSA
1511
1512	def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1513	'spherical harmonics preferred orientation (cylindrical symmetry only)'
1514	if 'T' in calcControls[hfx+'histType']:
1515	tth = parmDict[hfx+'2-theta']
1516	else:
1517	tth = refl[5+im]
1518	odfCor = 1.0
1519	H = refl[:3]
1520	cell = G2lat.Gmat2cell(g)
1521	Sangls = [0.,0.,0.]
1522	if 'Bragg' in calcControls[hfx+'instType']:
1523	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1524	IFCoup = True
1525	else:
1526	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1527	IFCoup = False
1528	phi,beta = G2lat.CrsAng(H,cell,SGData)
1529	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1530	SHnames = calcControls[phfx+'SHnames']
1531	for item in SHnames:
1532	L,N = eval(item.strip('C'))
1533	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1534	Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1535	Lnorm = G2lat.Lnorm(L)
1536	odfCor += parmDict[phfx+item]LnormKcl*Ksl
1537	return np.squeeze(odfCor)
1538
1539	def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1540	'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1541	if 'T' in calcControls[hfx+'histType']:
1542	tth = parmDict[hfx+'2-theta']
1543	else:
1544	tth = refl[5+im]
1545	odfCor = 1.0
1546	dFdODF = {}
1547	H = refl[:3]
1548	cell = G2lat.Gmat2cell(g)
1549	Sangls = [0.,0.,0.]
1550	if 'Bragg' in calcControls[hfx+'instType']:
1551	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1552	IFCoup = True
1553	else:
1554	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1555	IFCoup = False
1556	phi,beta = G2lat.CrsAng(H,cell,SGData)
1557	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1558	SHnames = calcControls[phfx+'SHnames']
1559	for item in SHnames:
1560	L,N = eval(item.strip('C'))
1561	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1562	Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1563	Lnorm = G2lat.Lnorm(L)
1564	odfCor += parmDict[phfx+item]LnormKcl*Ksl
1565	dFdODF[phfx+item] = KclKslLnorm
1566	return odfCor,dFdODF
1567
1568	def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1569	'March-Dollase preferred orientation correction'
1570	POcorr = 1.0
1571	MD = parmDict[phfx+'MD']
1572	if MD != 1.0:
1573	MDAxis = calcControls[phfx+'MDAxis']
1574	sumMD = 0
1575	for H in uniq:
1576	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1577	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1578	sumMD += A**3
1579	POcorr = sumMD/len(uniq)
1580	return POcorr
1581
1582	def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1583	'Needs a doc string'
1584	POcorr = 1.0
1585	POderv = {}
1586	if calcControls[phfx+'poType'] == 'MD':
1587	MD = parmDict[phfx+'MD']
1588	MDAxis = calcControls[phfx+'MDAxis']
1589	sumMD = 0
1590	sumdMD = 0
1591	for H in uniq:
1592	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1593	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1594	sumMD += A**3
1595	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
1596	POcorr = sumMD/len(uniq)
1597	POderv[phfx+'MD'] = sumdMD/len(uniq)
1598	else: #spherical harmonics
1599	if calcControls[phfx+'SHord']:
1600	POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1601	return POcorr,POderv
1602
1603	def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1604	'Needs a doc string'
1605	if 'Debye' in calcControls[hfx+'instType']:
1606	if 'T' in calcControls[hfx+'histType']:
1607	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1608	else:
1609	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1610	else:
1611	return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1612
1613	def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1614	'Needs a doc string'
1615	if 'Debye' in calcControls[hfx+'instType']:
1616	if 'T' in calcControls[hfx+'histType']:
1617	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1618	else:
1619	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1620	else:
1621	return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1622
1623	def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1624	'Needs a doc string'
1625	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1626	pi2 = np.sqrt(2./np.pi)
1627	if 'T' in calcControls[hfx+'histType']:
1628	sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1629	wave = refl[14+im]
1630	else: #'C'W
1631	sth2 = sind(refl[5+im]/2.)**2
1632	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1633	c2th = 1.-2.0*sth2
1634	flv2 = refl[9+im](wave/parmDict[pfx+'Vol'])*2
1635	if 'X' in calcControls[hfx+'histType']:
1636	flv2 = 0.079411(1.0+c2th**2)/2.0
1637	xfac = flv2*parmDict[phfx+'Extinction']
1638	exb = 1.0
1639	if xfac > -1.:
1640	exb = 1./np.sqrt(1.+xfac)
1641	exl = 1.0
1642	if 0 < xfac <= 1.:
1643	xn = np.array([xfac**(i+1) for i in range(6)])
1644	exl += np.sum(xn*coef)
1645	elif xfac > 1.:
1646	xfac2 = 1./np.sqrt(xfac)
1647	exl = pi2(1.-0.125/xfac)xfac2
1648	return exbsth2+exl(1.-sth2)
1649
1650	def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1651	'Needs a doc string'
1652	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1653	pi2 = np.sqrt(2./np.pi)
1654	if 'T' in calcControls[hfx+'histType']:
1655	sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1656	wave = refl[14+im]
1657	else: #'C'W
1658	sth2 = sind(refl[5+im]/2.)**2
1659	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1660	c2th = 1.-2.0*sth2
1661	flv2 = refl[9+im](wave/parmDict[pfx+'Vol'])*2
1662	if 'X' in calcControls[hfx+'histType']:
1663	flv2 = 0.079411(1.0+c2th**2)/2.0
1664	xfac = flv2*parmDict[phfx+'Extinction']
1665	dbde = -500.*flv2
1666	if xfac > -1.:
1667	dbde = -0.5flv2/np.sqrt(1.+xfac)*3
1668	dlde = 0.
1669	if 0 < xfac <= 1.:
1670	xn = np.array([iflv2xfac**i for i in [1,2,3,4,5,6]])
1671	dlde = np.sum(xn*coef)
1672	elif xfac > 1.:
1673	xfac2 = 1./np.sqrt(xfac)
1674	dlde = flv2pi2xfac2(-1./xfac+0.375/xfac*2)
1675
1676	return dbdesth2+dlde(1.-sth2)
1677
1678	def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1679	'Needs a doc string' #need powder extinction!
1680	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3+im] #scale*multiplicity
1681	if 'X' in parmDict[hfx+'Type']:
1682	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1683	POcorr = 1.0
1684	if pfx+'SHorder' in parmDict: #generalized spherical harmonics texture - takes precidence
1685	POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1686	elif calcControls[phfx+'poType'] == 'MD': #March-Dollase
1687	POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1688	elif calcControls[phfx+'SHord']: #cylindrical spherical harmonics
1689	POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1690	Icorr *= POcorr
1691	AbsCorr = 1.0
1692	AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1693	Icorr *= AbsCorr
1694	ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1695	Icorr *= ExtCorr
1696	return Icorr,POcorr,AbsCorr,ExtCorr
1697
1698	def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1699	'Needs a doc string' #need powder extinction derivs!
1700	dIdsh = 1./parmDict[hfx+'Scale']
1701	dIdsp = 1./parmDict[phfx+'Scale']
1702	if 'X' in parmDict[hfx+'Type']:
1703	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1704	dIdPola /= pola
1705	else: #'N'
1706	dIdPola = 0.0
1707	dFdODF = {}
1708	dFdSA = [0,0,0]
1709	dIdPO = {}
1710	if pfx+'SHorder' in parmDict:
1711	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1712	for iSH in dFdODF:
1713	dFdODF[iSH] /= odfCor
1714	for i in range(3):
1715	dFdSA[i] /= odfCor
1716	elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1717	POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1718	for iPO in dIdPO:
1719	dIdPO[iPO] /= POcorr
1720	if 'T' in parmDict[hfx+'Type']:
1721	dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1722	dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im] #/ext corr
1723	else:
1724	dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1725	dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im] #/ext corr
1726	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1727
1728	def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1729	'Needs a doc string'
1730	if 'C' in calcControls[hfx+'histType']: #All checked & OK
1731	costh = cosd(refl[5+im]/2.)
1732	#crystallite size
1733	if calcControls[phfx+'SizeType'] == 'isotropic':
1734	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1735	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1736	H = np.array(refl[:3])
1737	P = np.array(calcControls[phfx+'SizeAxis'])
1738	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1739	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
1740	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1741	else: #ellipsoidal crystallites
1742	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1743	H = np.array(refl[:3])
1744	lenR = G2pwd.ellipseSize(H,Sij,GB)
1745	Sgam = 1.8wave/(np.picosth*lenR)
1746	#microstrain
1747	if calcControls[phfx+'MustrainType'] == 'isotropic':
1748	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5+im]/2.)/np.pi
1749	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1750	H = np.array(refl[:3])
1751	P = np.array(calcControls[phfx+'MustrainAxis'])
1752	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1753	Si = parmDict[phfx+'Mustrain;i']
1754	Sa = parmDict[phfx+'Mustrain;a']
1755	Mgam = 0.018SiSatand(refl[5+im]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
1756	else: #generalized - P.W. Stephens model
1757	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1758	Sum = 0
1759	for i,strm in enumerate(Strms):
1760	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1761	Mgam = 0.018refl[4+im]2tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1762	elif 'T' in calcControls[hfx+'histType']: #All checked & OK
1763	#crystallite size
1764	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
1765	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/parmDict[phfx+'Size;i']
1766	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
1767	H = np.array(refl[:3])
1768	P = np.array(calcControls[phfx+'SizeAxis'])
1769	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1770	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]*2/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a'])
1771	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1772	else: #ellipsoidal crystallites #OK
1773	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1774	H = np.array(refl[:3])
1775	lenR = G2pwd.ellipseSize(H,Sij,GB)
1776	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/lenR
1777	#microstrain
1778	if calcControls[phfx+'MustrainType'] == 'isotropic': #OK
1779	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]*parmDict[phfx+'Mustrain;i']
1780	elif calcControls[phfx+'MustrainType'] == 'uniaxial': #OK
1781	H = np.array(refl[:3])
1782	P = np.array(calcControls[phfx+'MustrainAxis'])
1783	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1784	Si = parmDict[phfx+'Mustrain;i']
1785	Sa = parmDict[phfx+'Mustrain;a']
1786	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]SiSa/np.sqrt((SicosP)2+(SasinP)**2)
1787	else: #generalized - P.W. Stephens model OK
1788	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1789	Sum = 0
1790	for i,strm in enumerate(Strms):
1791	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1792	Mgam = 1.e-6parmDict[hfx+'difC']np.sqrt(Sum)refl[4+im]*3
1793
1794	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
1795	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
1796	sig /= ateln2
1797	return sig,gam
1798
1799	def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1800	'Needs a doc string'
1801	gamDict = {}
1802	sigDict = {}
1803	if 'C' in calcControls[hfx+'histType']: #All checked & OK
1804	costh = cosd(refl[5+im]/2.)
1805	tanth = tand(refl[5+im]/2.)
1806	#crystallite size derivatives
1807	if calcControls[phfx+'SizeType'] == 'isotropic':
1808	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1809	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
1810	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
1811	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1812	H = np.array(refl[:3])
1813	P = np.array(calcControls[phfx+'SizeAxis'])
1814	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1815	Si = parmDict[phfx+'Size;i']
1816	Sa = parmDict[phfx+'Size;a']
1817	gami = 1.8wave/(costhnp.piSiSa)
1818	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1819	Sgam = gami*sqtrm
1820	dsi = gamiSicosP**2/sqtrm-Sgam/Si
1821	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
1822	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1823	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1824	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1825	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1826	else: #ellipsoidal crystallites
1827	const = 1.8wave/(np.picosth)
1828	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1829	H = np.array(refl[:3])
1830	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1831	Sgam = const/lenR
1832	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1833	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1834	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1835	gamDict[phfx+'Size;mx'] = Sgam
1836	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
1837
1838	#microstrain derivatives
1839	if calcControls[phfx+'MustrainType'] == 'isotropic':
1840	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5+im]/2.)/np.pi
1841	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
1842	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
1843	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1844	H = np.array(refl[:3])
1845	P = np.array(calcControls[phfx+'MustrainAxis'])
1846	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1847	Si = parmDict[phfx+'Mustrain;i']
1848	Sa = parmDict[phfx+'Mustrain;a']
1849	gami = 0.018SiSa*tanth/np.pi
1850	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1851	Mgam = gami/sqtrm
1852	dsi = -gamiSicosP2/sqtrm3
1853	dsa = -gamiSasinP2/sqtrm3
1854	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1855	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1856	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1857	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1858	else: #generalized - P.W. Stephens model
1859	const = 0.018refl[4+im]2tanth/np.pi
1860	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1861	Sum = 0
1862	for i,strm in enumerate(Strms):
1863	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1864	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1865	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
1866	Mgam = const*np.sqrt(Sum)
1867	for i in range(len(Strms)):
1868	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1869	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
1870	gamDict[phfx+'Mustrain;mx'] = Mgam
1871	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1872	else: #'T'OF - All checked & OK
1873	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
1874	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/parmDict[phfx+'Size;i']
1875	gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1876	sigDict[phfx+'Size;i'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])*2/(ateln2parmDict[phfx+'Size;i'])
1877	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
1878	const = 1.e-4parmDict[hfx+'difC']refl[4+im]**2
1879	H = np.array(refl[:3])
1880	P = np.array(calcControls[phfx+'SizeAxis'])
1881	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1882	Si = parmDict[phfx+'Size;i']
1883	Sa = parmDict[phfx+'Size;a']
1884	gami = const/(Si*Sa)
1885	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1886	Sgam = gami*sqtrm
1887	dsi = gamiSicosP**2/sqtrm-Sgam/Si
1888	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
1889	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1890	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1891	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1892	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1893	else: #OK ellipsoidal crystallites
1894	const = 1.e-4parmDict[hfx+'difC']refl[4+im]**2
1895	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1896	H = np.array(refl[:3])
1897	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1898	Sgam = const/lenR
1899	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1900	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1901	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1902	gamDict[phfx+'Size;mx'] = Sgam #OK
1903	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2 #OK
1904
1905	#microstrain derivatives
1906	if calcControls[phfx+'MustrainType'] == 'isotropic':
1907	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]*parmDict[phfx+'Mustrain;i']
1908	gamDict[phfx+'Mustrain;i'] = 1.e-6refl[4+im]parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx'] #OK
1909	sigDict[phfx+'Mustrain;i'] = 2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])*2/(ateln2parmDict[phfx+'Mustrain;i'])
1910	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1911	H = np.array(refl[:3])
1912	P = np.array(calcControls[phfx+'MustrainAxis'])
1913	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1914	Si = parmDict[phfx+'Mustrain;i']
1915	Sa = parmDict[phfx+'Mustrain;a']
1916	gami = 1.e-6parmDict[hfx+'difC']refl[4+im]SiSa
1917	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1918	Mgam = gami/sqtrm
1919	dsi = -gamiSicosP2/sqtrm3
1920	dsa = -gamiSasinP2/sqtrm3
1921	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1922	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1923	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1924	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1925	else: #generalized - P.W. Stephens model OK
1926	pwrs = calcControls[phfx+'MuPwrs']
1927	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1928	const = 1.e-6parmDict[hfx+'difC']refl[4+im]**3
1929	Sum = 0
1930	for i,strm in enumerate(Strms):
1931	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1932	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1933	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
1934	Mgam = const*np.sqrt(Sum)
1935	for i in range(len(Strms)):
1936	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1937	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
1938	gamDict[phfx+'Mustrain;mx'] = Mgam
1939	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1940
1941	return sigDict,gamDict
1942
1943	def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1944	'Needs a doc string'
1945	if im:
1946	h,k,l,m = refl[:4]
1947	vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1948	d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1949	else:
1950	h,k,l = refl[:3]
1951	d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1952	refl[4+im] = d
1953	if 'C' in calcControls[hfx+'histType']:
1954	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
1955	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1956	if 'Bragg' in calcControls[hfx+'instType']:
1957	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1958	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
1959	else: #Debye-Scherrer - simple but maybe not right
1960	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1961	elif 'T' in calcControls[hfx+'histType']:
1962	pos = parmDict[hfx+'difC']d+parmDict[hfx+'difA']d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1963	#do I need sample position effects - maybe?
1964	return pos
1965
1966	def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1967	'Needs a doc string'
1968	dpr = 180./np.pi
1969	if im:
1970	h,k,l,m = refl[:4]
1971	vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1972	dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1973	else:
1974	m = 0
1975	h,k,l = refl[:3]
1976	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1977	dst = np.sqrt(dstsq)
1978	dsp = 1./dst
1979	if 'C' in calcControls[hfx+'histType']:
1980	pos = refl[5+im]-parmDict[hfx+'Zero']
1981	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
1982	dpdw = const*dst
1983	dpdA = np.array([h2,k2,l*2,hk,hl,kl])constwave/(2.0*dst)
1984	dpdZ = 1.0
1985	dpdV = np.array([2.hA[0]+kA[3]+lA[4],2kA[1]+hA[3]+lA[5],
1986	2lA[2]+hA[4]+kA[5]])mconstwave/(2.0dst)
1987	shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1988	if 'Bragg' in calcControls[hfx+'instType']:
1989	dpdSh = -4.shftcosd(pos/2.0)
1990	dpdTr = -shftsind(pos)100.0
1991	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1992	else: #Debye-Scherrer - simple but maybe not right
1993	dpdXd = -shft*cosd(pos)
1994	dpdYd = -shft*sind(pos)
1995	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1996	elif 'T' in calcControls[hfx+'histType']:
1997	dpdA = -np.array([h2,k2,l*2,hk,hl,kl])parmDict[hfx+'difC']dsp**3/2.
1998	dpdZ = 1.0
1999	dpdDC = dsp
2000	dpdDA = dsp**2
2001	dpdDB = 1./dsp
2002	dpdV = np.array([2.hA[0]+kA[3]+lA[4],2kA[1]+hA[3]+lA[5],
2003	2lA[2]+hA[4]+kA[5]])mparmDict[hfx+'difC']dsp**3/2.
2004	return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2005
2006	def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2007	'Needs a doc string'
2008	laue = SGData['SGLaue']
2009	uniq = SGData['SGUniq']
2010	h,k,l = refl[:3]
2011	if laue in ['m3','m3m']:
2012	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
2013	refl[4+im]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
2014	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2015	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
2016	elif laue in ['3R','3mR']:
2017	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
2018	elif laue in ['4/m','4/mmm']:
2019	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
2020	elif laue in ['mmm']:
2021	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2022	elif laue in ['2/m']:
2023	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2024	if uniq == 'a':
2025	Dij += parmDict[phfx+'D23']kl
2026	elif uniq == 'b':
2027	Dij += parmDict[phfx+'D13']hl
2028	elif uniq == 'c':
2029	Dij += parmDict[phfx+'D12']hk
2030	else:
2031	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
2032	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
2033	if 'C' in calcControls[hfx+'histType']:
2034	return -180.Dijrefl[4+im]*2tand(refl[5+im]/2.0)/np.pi
2035	else:
2036	return -DijparmDict[hfx+'difC']refl[4+im]**2/2.
2037
2038	def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2039	'Needs a doc string'
2040	laue = SGData['SGLaue']
2041	uniq = SGData['SGUniq']
2042	h,k,l = refl[:3]
2043	if laue in ['m3','m3m']:
2044	dDijDict = {phfx+'D11':h2+k2+l**2,
2045	phfx+'eA':refl[4+im]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
2046	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2047	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
2048	elif laue in ['3R','3mR']:
2049	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
2050	elif laue in ['4/m','4/mmm']:
2051	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
2052	elif laue in ['mmm']:
2053	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2054	elif laue in ['2/m']:
2055	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2056	if uniq == 'a':
2057	dDijDict[phfx+'D23'] = k*l
2058	elif uniq == 'b':
2059	dDijDict[phfx+'D13'] = h*l
2060	elif uniq == 'c':
2061	dDijDict[phfx+'D12'] = h*k
2062	else:
2063	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
2064	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
2065	if 'C' in calcControls[hfx+'histType']:
2066	for item in dDijDict:
2067	dDijDict[item] = 180.0refl[4+im]*2tand(refl[5+im]/2.0)/np.pi
2068	else:
2069	for item in dDijDict:
2070	dDijDict[item] = -parmDict[hfx+'difC']refl[4+im]**3/2.
2071	return dDijDict
2072
2073	def GetDij(phfx,SGData,parmDict):
2074	HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2075	return G2spc.HStrainVals(HSvals,SGData)
2076
2077	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2078	'Needs a doc string'
2079	histoList = Histograms.keys()
2080	histoList.sort()
2081	for histogram in histoList:
2082	if 'PWDR' in histogram[:4]:
2083	Histogram = Histograms[histogram]
2084	hId = Histogram['hId']
2085	hfx = ':%d:'%(hId)
2086	Limits = calcControls[hfx+'Limits']
2087	if 'C' in calcControls[hfx+'histType']:
2088	shl = max(parmDict[hfx+'SH/L'],0.0005)
2089	Ka2 = False
2090	kRatio = 0.0
2091	if hfx+'Lam1' in parmDict.keys():
2092	Ka2 = True
2093	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2094	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2095	x,y,w,yc,yb,yd = Histogram['Data']
2096	xB = np.searchsorted(x,Limits[0])
2097	xF = np.searchsorted(x,Limits[1])
2098	ymb = np.array(y-yb)
2099	ymb = np.where(ymb,ymb,1.0)
2100	ycmb = np.array(yc-yb)
2101	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
2102	refLists = Histogram['Reflection Lists']
2103	for phase in refLists:
2104	if phase not in Phases: #skips deleted or renamed phases silently!
2105	continue
2106	Phase = Phases[phase]
2107	im = 0
2108	if Phase['General']['Type'] in ['modulated','magnetic']:
2109	im = 1
2110	pId = Phase['pId']
2111	phfx = '%d:%d:'%(pId,hId)
2112	refDict = refLists[phase]
2113	sumFo = 0.0
2114	sumdF = 0.0
2115	sumFosq = 0.0
2116	sumdFsq = 0.0
2117	sumInt = 0.0
2118	for refl in refDict['RefList']:
2119	if 'C' in calcControls[hfx+'histType']:
2120	yp = np.zeros_like(yb)
2121	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2122	iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2123	iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2124	iFin2 = iFin
2125	if not iBeg+iFin: #peak below low limit - skip peak
2126	continue
2127	elif not iBeg-iFin: #peak above high limit - done
2128	break
2129	elif iBeg < iFin:
2130	yp[iBeg:iFin] = refl[11+im]refl[9+im]G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
2131	sumInt += refl[11+im]*refl[9+im]
2132	if Ka2:
2133	pos2 = refl[5+im]+lamRatiotand(refl[5+im]/2.0) # + 360/pi Dlam/lam * tan(th)
2134	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2135	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
2136	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
2137	if iFin2 > iBeg2:
2138	yp[iBeg2:iFin2] += refl[11+im]refl[9+im]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2])) #and here
2139	sumInt += refl[11+im]refl[9+im]kRatio
2140	refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[11+im](1.+kRatio)),0.0))
2141
2142	elif 'T' in calcControls[hfx+'histType']:
2143	yp = np.zeros_like(yb)
2144	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2145	iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2146	iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2147	if iBeg < iFin:
2148	yp[iBeg:iFin] = refl[11+im]refl[9+im]G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin])) #>90% of time spent here
2149	refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
2150	sumInt += refl[11+im]*refl[9+im]
2151	Fo = np.sqrt(np.abs(refl[8+im]))
2152	Fc = np.sqrt(np.abs(refl[9]+im))
2153	sumFo += Fo
2154	sumFosq += refl[8+im]**2
2155	sumdF += np.abs(Fo-Fc)
2156	sumdFsq += (refl[8+im]-refl[9+im])**2
2157	if sumFo:
2158	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2159	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2160	else:
2161	Histogram['Residuals'][phfx+'Rf'] = 100.
2162	Histogram['Residuals'][phfx+'Rf^2'] = 100.
2163	Histogram['Residuals'][phfx+'sumInt'] = sumInt
2164	Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
2165	Histogram['Residuals']['hId'] = hId
2166	elif 'HKLF' in histogram[:4]:
2167	Histogram = Histograms[histogram]
2168	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
2169
2170	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2171	'Needs a doc string'
2172
2173	def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
2174	U = parmDict[hfx+'U']
2175	V = parmDict[hfx+'V']
2176	W = parmDict[hfx+'W']
2177	X = parmDict[hfx+'X']
2178	Y = parmDict[hfx+'Y']
2179	tanPos = tand(refl[5+im]/2.0)
2180	Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2181	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
2182	sig = max(0.001,sig)
2183	gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam #save peak gamma
2184	gam = max(0.001,gam)
2185	return sig,gam
2186
2187	def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
2188	sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']refl[4+im]*2+ \
2189	parmDict[hfx+'sig-2']refl[4+im]4+parmDict[hfx+'sig-q']/refl[4+im]*2
2190	gam = parmDict[hfx+'X']refl[4+im]+parmDict[hfx+'Y']refl[4+im]**2
2191	Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2192	sig += Ssig
2193	gam += Sgam
2194	return sig,gam
2195
2196	def GetReflAlpBet(refl,im,hfx,parmDict):
2197	alp = parmDict[hfx+'alpha']/refl[4+im]
2198	bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]4+parmDict[hfx+'beta-q']/refl[4+im]2
2199	return alp,bet
2200
2201	hId = Histogram['hId']
2202	hfx = ':%d:'%(hId)
2203	bakType = calcControls[hfx+'bakType']
2204	yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2205	yc = np.zeros_like(yb)
2206	cw = np.diff(x)
2207	cw = np.append(cw,cw[-1])
2208
2209	if 'C' in calcControls[hfx+'histType']:
2210	shl = max(parmDict[hfx+'SH/L'],0.002)
2211	Ka2 = False
2212	if hfx+'Lam1' in parmDict.keys():
2213	wave = parmDict[hfx+'Lam1']
2214	Ka2 = True
2215	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2216	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2217	else:
2218	wave = parmDict[hfx+'Lam']
2219	for phase in Histogram['Reflection Lists']:
2220	refDict = Histogram['Reflection Lists'][phase]
2221	if phase not in Phases: #skips deleted or renamed phases silently!
2222	continue
2223	Phase = Phases[phase]
2224	pId = Phase['pId']
2225	pfx = '%d::'%(pId)
2226	phfx = '%d:%d:'%(pId,hId)
2227	hfx = ':%d:'%(hId)
2228	SGData = Phase['General']['SGData']
2229	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2230	im = 0
2231	if Phase['General']['Type'] in ['modulated','magnetic']:
2232	SSGData = Phase['General']['SSGData']
2233	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2234	im = 1 #offset in SS reflection list
2235	#??
2236	Dij = GetDij(phfx,SGData,parmDict)
2237	A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2238	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2239	if np.any(np.diag(G)<0.):
2240	raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
2241	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2242	Vst = np.sqrt(nl.det(G)) #V*
2243	if not Phase['General'].get('doPawley'):
2244	time0 = time.time()
2245	if im:
2246	SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2247	else:
2248	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2249	# print 'sf calc time: %.3fs'%(time.time()-time0)
2250	time0 = time.time()
2251	badPeak = False
2252	for iref,refl in enumerate(refDict['RefList']):
2253	if 'C' in calcControls[hfx+'histType']:
2254	if im:
2255	h,k,l,m = refl[:4]
2256	else:
2257	h,k,l = refl[:3]
2258	Uniq = np.inner(refl[:3],SGMT)
2259	refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict) #corrected reflection position
2260	Lorenz = 1./(2.sind(refl[5+im]/2.)2cosd(refl[5+im]/2.)) #Lorentz correction
2261	refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict) #peak sig & gam
2262	refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2263	refl[11+im] = VstLorenz
2264
2265	if Phase['General'].get('doPawley'):
2266	try:
2267	if im:
2268	pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2269	else:
2270	pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2271	refl[9+im] = parmDict[pInd]
2272	except KeyError:
2273	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2274	continue
2275	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2276	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2277	iFin = np.searchsorted(x,refl[5+im]+fmax)
2278	if not iBeg+iFin: #peak below low limit - skip peak
2279	continue
2280	elif not iBeg-iFin: #peak above high limit - done
2281	break
2282	elif iBeg > iFin: #bad peak coeff - skip
2283	badPeak = True
2284	continue
2285	yc[iBeg:iFin] += refl[11+im]refl[9+im]G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here
2286	if Ka2:
2287	pos2 = refl[5+im]+lamRatiotand(refl[5+im]/2.0) # + 360/pi Dlam/lam * tan(th)
2288	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2289	iBeg = np.searchsorted(x,pos2-fmin)
2290	iFin = np.searchsorted(x,pos2+fmax)
2291	if not iBeg+iFin: #peak below low limit - skip peak
2292	continue
2293	elif not iBeg-iFin: #peak above high limit - done
2294	return yc,yb
2295	elif iBeg > iFin: #bad peak coeff - skip
2296	continue
2297	yc[iBeg:iFin] += refl[11+im]refl[9+im]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin])) #and here
2298	elif 'T' in calcControls[hfx+'histType']:
2299	h,k,l = refl[:3]
2300	Uniq = np.inner(refl[:3],SGMT)
2301	refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict) #corrected reflection position
2302	Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)refl[4+im]*4 #TOF Lorentz correction
2303	# refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict) #apply hydrostatic strain shift
2304	refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict) #peak sig & gam
2305	refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2306	refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2307	refl[11+im] = VstLorenz
2308	if Phase['General'].get('doPawley'):
2309	try:
2310	if im:
2311	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2312	else:
2313	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2314	refl[9+im] = parmDict[pInd]
2315	except KeyError:
2316	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2317	continue
2318	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2319	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2320	iFin = np.searchsorted(x,refl[5+im]+fmax)
2321	if not iBeg+iFin: #peak below low limit - skip peak
2322	continue
2323	elif not iBeg-iFin: #peak above high limit - done
2324	break
2325	elif iBeg > iFin: #bad peak coeff - skip
2326	badPeak = True
2327	continue
2328	yc[iBeg:iFin] += refl[11+im]refl[9+im]G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
2329	# print 'profile calc time: %.3fs'%(time.time()-time0)
2330	if badPeak:
2331	print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped'
2332	return yc,yb
2333
2334	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2335	'Needs a doc string'
2336
2337	def cellVaryDerv(pfx,SGData,dpdA):
2338	if SGData['SGLaue'] in ['-1',]:
2339	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2340	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2341	elif SGData['SGLaue'] in ['2/m',]:
2342	if SGData['SGUniq'] == 'a':
2343	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2344	elif SGData['SGUniq'] == 'b':
2345	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2346	else:
2347	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2348	elif SGData['SGLaue'] in ['mmm',]:
2349	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2350	elif SGData['SGLaue'] in ['4/m','4/mmm']:
2351	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2352	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2353	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2354	elif SGData['SGLaue'] in ['3R', '3mR']:
2355	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
2356	elif SGData['SGLaue'] in ['m3m','m3']:
2357	return [[pfx+'A0',dpdA[0]]]
2358
2359	# create a list of dependent variables and set up a dictionary to hold their derivatives
2360	dependentVars = G2mv.GetDependentVars()
2361	depDerivDict = {}
2362	for j in dependentVars:
2363	depDerivDict[j] = np.zeros(shape=(len(x)))
2364	#print 'dependent vars',dependentVars
2365	lenX = len(x)
2366	hId = Histogram['hId']
2367	hfx = ':%d:'%(hId)
2368	bakType = calcControls[hfx+'bakType']
2369	dMdv = np.zeros(shape=(len(varylist),len(x)))
2370	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2371	if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2372	bBpos =varylist.index(hfx+'Back;0')
2373	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2374	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2375	for name in varylist:
2376	if 'Debye' in name:
2377	id = int(name.split(';')[-1])
2378	parm = name[:int(name.rindex(';'))]
2379	ip = names.index(parm)
2380	dMdv[varylist.index(name)] = dMddb[3*id+ip]
2381	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2382	for name in varylist:
2383	if 'BkPk' in name:
2384	parm,id = name.split(';')
2385	id = int(id)
2386	if parm in names:
2387	ip = names.index(parm)
2388	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2389	cw = np.diff(x)
2390	cw = np.append(cw,cw[-1])
2391	Ka2 = False #also for TOF!
2392	if 'C' in calcControls[hfx+'histType']:
2393	shl = max(parmDict[hfx+'SH/L'],0.002)
2394	if hfx+'Lam1' in parmDict.keys():
2395	wave = parmDict[hfx+'Lam1']
2396	Ka2 = True
2397	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2398	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2399	else:
2400	wave = parmDict[hfx+'Lam']
2401	for phase in Histogram['Reflection Lists']:
2402	refDict = Histogram['Reflection Lists'][phase]
2403	if phase not in Phases: #skips deleted or renamed phases silently!
2404	continue
2405	Phase = Phases[phase]
2406	SGData = Phase['General']['SGData']
2407	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2408	im = 0
2409	if Phase['General']['Type'] in ['modulated','magnetic']:
2410	SSGData = Phase['General']['SSGData']
2411	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2412	im = 1 #offset in SS reflection list
2413	#??
2414	pId = Phase['pId']
2415	pfx = '%d::'%(pId)
2416	phfx = '%d:%d:'%(pId,hId)
2417	Dij = GetDij(phfx,SGData,parmDict)
2418	A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2419	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2420	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2421	if not Phase['General'].get('doPawley'):
2422	time0 = time.time()
2423	if im:
2424	dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2425	else:
2426	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2427	# print 'sf-derv time %.3fs'%(time.time()-time0)
2428	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2429	time0 = time.time()
2430	for iref,refl in enumerate(refDict['RefList']):
2431	if im:
2432	h,k,l,m = refl[:4]
2433	else:
2434	h,k,l = refl[:3]
2435	Uniq = np.inner(refl[:3],SGMT)
2436	if 'T' in calcControls[hfx+'histType']:
2437	wave = refl[14+im]
2438	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2439	if 'C' in calcControls[hfx+'histType']: #CW powder
2440	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2441	else: #'T'OF
2442	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2443	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2444	iFin = np.searchsorted(x,refl[5+im]+fmax)
2445	if not iBeg+iFin: #peak below low limit - skip peak
2446	continue
2447	elif not iBeg-iFin: #peak above high limit - done
2448	break
2449	pos = refl[5+im]
2450	if 'C' in calcControls[hfx+'histType']:
2451	tanth = tand(pos/2.0)
2452	costh = cosd(pos/2.0)
2453	lenBF = iFin-iBeg
2454	dMdpk = np.zeros(shape=(6,lenBF))
2455	dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2456	for i in range(5):
2457	dMdpk[i] += 100.cw[iBeg:iFin]refl[11+im]refl[9+im]dMdipk[i]
2458	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2459	if Ka2:
2460	pos2 = refl[5+im]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
2461	iBeg2 = np.searchsorted(x,pos2-fmin)
2462	iFin2 = np.searchsorted(x,pos2+fmax)
2463	if iBeg2-iFin2:
2464	lenBF2 = iFin2-iBeg2
2465	dMdpk2 = np.zeros(shape=(6,lenBF2))
2466	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2467	for i in range(5):
2468	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[11+im]refl[9+im]kRatio*dMdipk2[i]
2469	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[11+im]*dMdipk2[0]
2470	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2471	else: #'T'OF
2472	lenBF = iFin-iBeg
2473	if lenBF < 0: #bad peak coeff
2474	break
2475	dMdpk = np.zeros(shape=(6,lenBF))
2476	dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2477	for i in range(6):
2478	dMdpk[i] += refl[11+im]refl[9+im]dMdipk[i] #cw[iBeg:iFin]*
2479	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
2480	if Phase['General'].get('doPawley'):
2481	dMdpw = np.zeros(len(x))
2482	try:
2483	if im:
2484	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2485	else:
2486	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2487	idx = varylist.index(pIdx)
2488	dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2489	if Ka2: #not for TOF either
2490	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2491	dMdv[idx] = dMdpw
2492	except: # ValueError:
2493	pass
2494	if 'C' in calcControls[hfx+'histType']:
2495	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2496	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2497	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2498	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2499	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2500	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2501	hfx+'DisplaceY':[dpdY,'pos'],}
2502	if 'Bragg' in calcControls[hfx+'instType']:
2503	names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2504	hfx+'SurfRoughB':[dFdAb[1],'int'],})
2505	else:
2506	names.update({hfx+'Absorption':[dFdAb,'int'],})
2507	else: #'T'OF
2508	dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2509	names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2510	hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2511	hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2512	hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2513	hfx+'beta-q':[1./refl[4+im]2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]2,'sig'],
2514	hfx+'sig-2':[refl[4+im]4,'sig'],hfx+'sig-q':[1./refl[4+im]2,'sig'],
2515	hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2516	for name in names:
2517	item = names[name]
2518	if name in varylist:
2519	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2520	if Ka2 and iFin2-iBeg2:
2521	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2522	elif name in dependentVars:
2523	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2524	if Ka2 and iFin2-iBeg2:
2525	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2526	for iPO in dIdPO:
2527	if iPO in varylist:
2528	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2529	if Ka2 and iFin2-iBeg2:
2530	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2531	elif iPO in dependentVars:
2532	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2533	if Ka2 and iFin2-iBeg2:
2534	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2535	for i,name in enumerate(['omega','chi','phi']):
2536	aname = pfx+'SH '+name
2537	if aname in varylist:
2538	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2539	if Ka2 and iFin2-iBeg2:
2540	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2541	elif aname in dependentVars:
2542	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2543	if Ka2 and iFin2-iBeg2:
2544	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2545	for iSH in dFdODF:
2546	if iSH in varylist:
2547	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2548	if Ka2 and iFin2-iBeg2:
2549	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2550	elif iSH in dependentVars:
2551	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2552	if Ka2 and iFin2-iBeg2:
2553	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2554	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2555	for name,dpdA in cellDervNames:
2556	if name in varylist:
2557	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2558	if Ka2 and iFin2-iBeg2:
2559	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2560	elif name in dependentVars:
2561	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2562	if Ka2 and iFin2-iBeg2:
2563	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2564	dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2565	for name in dDijDict:
2566	if name in varylist:
2567	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2568	if Ka2 and iFin2-iBeg2:
2569	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2570	elif name in dependentVars:
2571	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2572	if Ka2 and iFin2-iBeg2:
2573	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2574	for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2575	if name in varylist:
2576	dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2577	if Ka2 and iFin2-iBeg2:
2578	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2579	elif name in dependentVars:
2580	depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2581	if Ka2 and iFin2-iBeg2:
2582	depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2583	if 'C' in calcControls[hfx+'histType']:
2584	sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2585	else: #'T'OF
2586	sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2587	for name in gamDict:
2588	if name in varylist:
2589	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2590	if Ka2 and iFin2-iBeg2:
2591	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2592	elif name in dependentVars:
2593	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2594	if Ka2 and iFin2-iBeg2:
2595	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2596	for name in sigDict:
2597	if name in varylist:
2598	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2599	if Ka2 and iFin2-iBeg2:
2600	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2601	elif name in dependentVars:
2602	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2603	if Ka2 and iFin2-iBeg2:
2604	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2605	for name in ['BabA','BabU']:
2606	if refl[9+im]:
2607	if phfx+name in varylist:
2608	dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict['int']/refl[9+im]
2609	if Ka2 and iFin2-iBeg2:
2610	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict2['int']/refl[9+im]
2611	elif phfx+name in dependentVars:
2612	depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict['int']/refl[9+im]
2613	if Ka2 and iFin2-iBeg2:
2614	depDerivDict[phfx+name][iBeg2:iFin2] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict2['int']/refl[9+im]
2615	if not Phase['General'].get('doPawley'):
2616	#do atom derivatives - for RB,F,X & U so far
2617	corr = dervDict['int']/refl[9+im]
2618	if Ka2 and iFin2-iBeg2:
2619	corr2 = dervDict2['int']/refl[9+im]
2620	for name in varylist+dependentVars:
2621	if '::RBV;' in name:
2622	pass
2623	else:
2624	try:
2625	aname = name.split(pfx)[1][:2]
2626	if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2627	except IndexError:
2628	continue
2629	if name in varylist:
2630	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2631	if Ka2 and iFin2-iBeg2:
2632	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2633	elif name in dependentVars:
2634	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2635	if Ka2 and iFin2-iBeg2:
2636	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2637	# print 'profile derv time: %.3fs'%(time.time()-time0)
2638	# now process derivatives in constraints
2639	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2640	return dMdv
2641
2642	def UserRejectHKL(ref,im,userReject):
2643	if ref[5+im]/ref[6+im] < userReject['minF/sig']:
2644	return False
2645	elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
2646	return False
2647	elif ref[11+im] < userReject['MinExt']:
2648	return False
2649	elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
2650	return False
2651	return True
2652
2653	def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2654	'''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2655	model (M) with respect to all parameters. Independent and dependant dM/dp arrays
2656	are returned to either dervRefine or HessRefine.
2657
2658	:returns:
2659	'''
2660	nobs = Histogram['Residuals']['Nobs']
2661	hId = Histogram['hId']
2662	hfx = ':%d:'%(hId)
2663	pfx = '%d::'%(Phase['pId'])
2664	phfx = '%d:%d:'%(Phase['pId'],hId)
2665	SGData = Phase['General']['SGData']
2666	im = 0
2667	if Phase['General']['Type'] in ['modulated','magnetic']:
2668	SSGData = Phase['General']['SSGData']
2669	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2670	im = 1 #offset in SS reflection list
2671	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2672	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2673	refDict = Histogram['Data']
2674	if parmDict[phfx+'Scale'] < 0.:
2675	parmDict[phfx+'Scale'] = .001
2676	if im:
2677	dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2678	else:
2679	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2680	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2681	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2682	dependentVars = G2mv.GetDependentVars()
2683	depDerivDict = {}
2684	for j in dependentVars:
2685	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2686	wdf = np.zeros(len(refDict['RefList']))
2687	if calcControls['F**2']:
2688	for iref,ref in enumerate(refDict['RefList']):
2689	if ref[6+im] > 0:
2690	dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2691	w = 1.0/ref[6+im]
2692	if ref[3+im] > 0:
2693	wdf[iref] = w*(ref[5+im]-ref[7+im])
2694	for j,var in enumerate(varylist):
2695	if var in dFdvDict:
2696	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2697	for var in dependentVars:
2698	if var in dFdvDict:
2699	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2700	if phfx+'Scale' in varylist:
2701	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2702	elif phfx+'Scale' in dependentVars:
2703	depDerivDict[phfx+'Scale'][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2704	for item in ['Ep','Es','Eg']:
2705	if phfx+item in varylist and phfx+item in dervDict:
2706	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] #OK
2707	elif phfx+item in dependentVars and phfx+item in dervDict:
2708	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im] #OK
2709	for item in ['BabA','BabU']:
2710	if phfx+item in varylist:
2711	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2712	elif phfx+item in dependentVars:
2713	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2714	else: #F refinement
2715	for iref,ref in enumerate(refDict['RefList']):
2716	if ref[5+im] > 0.:
2717	dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2718	Fo = np.sqrt(ref[5+im])
2719	Fc = np.sqrt(ref[7+im])
2720	w = 1.0/ref[6+im]
2721	if ref[3+im] > 0:
2722	wdf[iref] = 2.0Fcw*(Fo-Fc)
2723	for j,var in enumerate(varylist):
2724	if var in dFdvDict:
2725	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2726	for var in dependentVars:
2727	if var in dFdvDict:
2728	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2729	if phfx+'Scale' in varylist:
2730	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2731	elif phfx+'Scale' in dependentVars:
2732	depDerivDict[phfx+'Scale'][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2733	for item in ['Ep','Es','Eg']: #OK!
2734	if phfx+item in varylist and phfx+item in dervDict:
2735	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]
2736	elif phfx+item in dependentVars and phfx+item in dervDict:
2737	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2738	for item in ['BabA','BabU']:
2739	if phfx+item in varylist:
2740	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2741	elif phfx+item in dependentVars:
2742	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2743	return dMdvh,depDerivDict,wdf
2744
2745
2746	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2747	'''Loop over histograms and compute derivatives of the fitting
2748	model (M) with respect to all parameters. Results are returned in
2749	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2750	n is the number of parameters and m is the number of data
2751	points. This can exceed memory when m gets large. This routine is
2752	used when refinement derivatives are selected as "analtytic
2753	Jacobian" in Controls.
2754
2755	:returns: Jacobian numpy.array dMdv for all histograms concatinated
2756	'''
2757	parmDict.update(zip(varylist,values))
2758	G2mv.Dict2Map(parmDict,varylist)
2759	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2760	nvar = len(varylist)
2761	dMdv = np.empty(0)
2762	histoList = Histograms.keys()
2763	histoList.sort()
2764	for histogram in histoList:
2765	if 'PWDR' in histogram[:4]:
2766	Histogram = Histograms[histogram]
2767	hId = Histogram['hId']
2768	hfx = ':%d:'%(hId)
2769	wtFactor = calcControls[hfx+'wtFactor']
2770	Limits = calcControls[hfx+'Limits']
2771	x,y,w,yc,yb,yd = Histogram['Data']
2772	xB = np.searchsorted(x,Limits[0])
2773	xF = np.searchsorted(x,Limits[1])
2774	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2775	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2776	elif 'HKLF' in histogram[:4]:
2777	Histogram = Histograms[histogram]
2778	phase = Histogram['Reflection Lists']
2779	Phase = Phases[phase]
2780	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2781	hfx = ':%d:'%(Histogram['hId'])
2782	wtFactor = calcControls[hfx+'wtFactor']
2783	# now process derivatives in constraints
2784	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2785	else:
2786	continue #skip non-histogram entries
2787	if len(dMdv):
2788	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2789	else:
2790	dMdv = np.sqrt(wtFactor)*dMdvh
2791
2792	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2793	if np.any(pVals):
2794	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2795	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2796
2797	return dMdv
2798
2799	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2800	'''Loop over histograms and compute derivatives of the fitting
2801	model (M) with respect to all parameters. For each histogram, the
2802	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2803	number of parameters and m is the number of data points *in the
2804	histogram*. The (n by n) Hessian is computed from each Jacobian
2805	and it is returned. This routine is used when refinement
2806	derivatives are selected as "analtytic Hessian" in Controls.
2807
2808	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2809	'''
2810	parmDict.update(zip(varylist,values))
2811	G2mv.Dict2Map(parmDict,varylist)
2812	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2813	#fixup H atom positions here?
2814	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
2815	nvar = len(varylist)
2816	Hess = np.empty(0)
2817	histoList = Histograms.keys()
2818	histoList.sort()
2819	for histogram in histoList:
2820	if 'PWDR' in histogram[:4]:
2821	Histogram = Histograms[histogram]
2822	hId = Histogram['hId']
2823	hfx = ':%d:'%(hId)
2824	wtFactor = calcControls[hfx+'wtFactor']
2825	Limits = calcControls[hfx+'Limits']
2826	x,y,w,yc,yb,yd = Histogram['Data']
2827	W = wtFactor*w
2828	dy = y-yc
2829	xB = np.searchsorted(x,Limits[0])
2830	xF = np.searchsorted(x,Limits[1])
2831	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2832	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2833	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2834	Dy = dy[xB:xF][np.newaxis,:]
2835	dMdvh *= Wt
2836	if dlg:
2837	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2838	if len(Hess):
2839	Hess += np.inner(dMdvh,dMdvh)
2840	dMdvh = WtDy
2841	Vec += np.sum(dMdvh,axis=1)
2842	else:
2843	Hess = np.inner(dMdvh,dMdvh)
2844	dMdvh = WtDy
2845	Vec = np.sum(dMdvh,axis=1)
2846	elif 'HKLF' in histogram[:4]:
2847	Histogram = Histograms[histogram]
2848	phase = Histogram['Reflection Lists']
2849	Phase = Phases[phase]
2850	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2851	hId = Histogram['hId']
2852	hfx = ':%d:'%(Histogram['hId'])
2853	wtFactor = calcControls[hfx+'wtFactor']
2854	# now process derivatives in constraints
2855	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2856	# print 'matrix build time: %.3f'%(time.time()-time0)
2857
2858	if dlg:
2859	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2860	if len(Hess):
2861	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
2862	Hess += wtFactor*np.inner(dMdvh,dMdvh)
2863	else:
2864	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
2865	Hess = wtFactor*np.inner(dMdvh,dMdvh)
2866	else:
2867	continue #skip non-histogram entries
2868	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2869	if np.any(pVals):
2870	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2871	Vec += np.sum(dpdvpWtpVals,axis=1)
2872	Hess += np.inner(dpdv*pWt,dpdv)
2873	return Vec,Hess
2874
2875	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):
2876	'''Computes the point-by-point discrepancies between every data point in every histogram
2877	and the observed value
2878
2879	:returns: an np array of differences between observed and computed diffraction values.
2880	'''
2881	Values2Dict(parmDict, varylist, values)
2882	G2mv.Dict2Map(parmDict,varylist)
2883	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2884	M = np.empty(0)
2885	SumwYo = 0
2886	Nobs = 0
2887	Nrej = 0
2888	Next = 0
2889	ApplyRBModels(parmDict,Phases,rigidbodyDict)
2890	#fixup Hatom positions here....
2891	histoList = Histograms.keys()
2892	histoList.sort()
2893	for histogram in histoList:
2894	if 'PWDR' in histogram[:4]:
2895	Histogram = Histograms[histogram]
2896	hId = Histogram['hId']
2897	hfx = ':%d:'%(hId)
2898	wtFactor = calcControls[hfx+'wtFactor']
2899	Limits = calcControls[hfx+'Limits']
2900	x,y,w,yc,yb,yd = Histogram['Data']
2901	yc *= 0.0 #zero full calcd profiles
2902	yb *= 0.0
2903	yd *= 0.0
2904	xB = np.searchsorted(x,Limits[0])
2905	xF = np.searchsorted(x,Limits[1])
2906	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2907	varylist,Histogram,Phases,calcControls,pawleyLookup)
2908	yc[xB:xF] += yb[xB:xF]
2909	if not np.any(y): #fill dummy data
2910	rv = st.poisson(yc[xB:xF])
2911	y[xB:xF] = rv.rvs()
2912	Z = np.ones_like(yc[xB:xF])
2913	Z[1::2] *= -1
2914	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2915	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2916	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2917	W = wtFactor*w
2918	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2919	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2920	Nobs += Histogram['Residuals']['Nobs']
2921	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
2922	SumwYo += Histogram['Residuals']['sumwYo']
2923	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2924	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
2925	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2926	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
2927	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.))
2928	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.))
2929	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2930	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2931	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2932	if dlg:
2933	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2934	M = np.concatenate((M,wdy))
2935	#end of PWDR processing
2936	elif 'HKLF' in histogram[:4]:
2937	Histogram = Histograms[histogram]
2938	Histogram['Residuals'] = {}
2939	phase = Histogram['Reflection Lists']
2940	Phase = Phases[phase]
2941	hId = Histogram['hId']
2942	hfx = ':%d:'%(hId)
2943	wtFactor = calcControls[hfx+'wtFactor']
2944	pfx = '%d::'%(Phase['pId'])
2945	phfx = '%d:%d:'%(Phase['pId'],hId)
2946	SGData = Phase['General']['SGData']
2947	im = 0
2948	if parmDict[phfx+'Scale'] < 0.:
2949	parmDict[phfx+'Scale'] = .001
2950	if Phase['General']['Type'] in ['modulated','magnetic']:
2951	SSGData = Phase['General']['SSGData']
2952	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2953	im = 1 #offset in SS reflection list
2954	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2955	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2956	refDict = Histogram['Data']
2957	time0 = time.time()
2958	if im:
2959	SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2960	else:
2961	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2962	# print 'sf-calc time: %.3f'%(time.time()-time0)
2963	df = np.zeros(len(refDict['RefList']))
2964	sumwYo = 0
2965	sumFo = 0
2966	sumFo2 = 0
2967	sumdF = 0
2968	sumdF2 = 0
2969	if im:
2970	sumSSFo = np.zeros(10)
2971	sumSSFo2 = np.zeros(10)
2972	sumSSdF = np.zeros(10)
2973	sumSSdF2 = np.zeros(10)
2974	sumSSwYo = np.zeros(10)
2975	sumSSwdf2 = np.zeros(10)
2976	SSnobs = np.zeros(10)
2977	nobs = 0
2978	nrej = 0
2979	next = 0
2980	maxH = 0
2981	if calcControls['F**2']:
2982	for i,ref in enumerate(refDict['RefList']):
2983	if ref[6+im] > 0:
2984	ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2985	w = 1.0/ref[6+im] # 1/sig(F^2)
2986	ref[7+im] = parmDict[phfx+'Scale']ref[11+im] #correct Fc^2 for extinction
2987	ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2988	if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0)
2989	ref[3+im] = abs(ref[3+im]) #mark as allowed
2990	Fo = np.sqrt(ref[5+im])
2991	sumFo += Fo
2992	sumFo2 += ref[5+im]
2993	sumdF += abs(Fo-np.sqrt(ref[7+im]))
2994	sumdF2 += abs(ref[5+im]-ref[7+im])
2995	nobs += 1
2996	df[i] = -w*(ref[5+im]-ref[7+im])
2997	sumwYo += (wref[5+im])2 #wFo^2
2998	if im: #accumulate super lattice sums
2999	ind = int(abs(ref[3]))
3000	sumSSFo[ind] += Fo
3001	sumSSFo2[ind] += ref[5+im]
3002	sumSSdF[ind] += abs(Fo-np.sqrt(ref[7+im]))
3003	sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
3004	sumSSwYo[ind] += (wref[5+im])2 #wFo^2
3005	sumSSwdf2[ind] += df[i]**2
3006	SSnobs[ind] += 1
3007	maxH = max(maxH,ind)
3008	else:
3009	if ref[3+im]:
3010	ref[3+im] = -abs(ref[3+im]) #mark as rejected
3011	nrej += 1
3012	else: #sp.gp.extinct
3013	next += 1
3014	else:
3015	for i,ref in enumerate(refDict['RefList']):
3016	if ref[5+im] > 0.:
3017	ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
3018	ref[7+im] = parmDict[phfx+'Scale']ref[11+im] #correct Fc^2 for extinction
3019	ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
3020	Fo = np.sqrt(ref[5+im])
3021	Fc = np.sqrt(ref[7+im])
3022	w = 2.0*Fo/ref[6+im] # 1/sig(F)?
3023	if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0)
3024	ref[3+im] = abs(ref[3+im]) #mark as allowed
3025	sumFo += Fo
3026	sumFo2 += ref[5+im]
3027	sumdF += abs(Fo-Fc)
3028	sumdF2 += abs(ref[5+im]-ref[7+im])
3029	nobs += 1
3030	df[i] = -w*(Fo-Fc)
3031	sumwYo += (wFo)*2
3032	if im:
3033	ind = int(abs(ref[3]))
3034	sumSSFo[ind] += Fo
3035	sumSSFo2[ind] += ref[5+im]
3036	sumSSdF[ind] += abs(Fo-Fc)
3037	sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
3038	sumSSwYo[ind] += (wFo)*2
3039	sumSSwdf2[ind] += df[i]**2
3040	SSnobs[ind] += 1
3041	maxH = max(maxH,ind)
3042	else:
3043	if ref[3+im]:
3044	ref[3+im] = -abs(ref[3+im]) #mark as rejected
3045	nrej += 1
3046	else: #sp.gp.extinct
3047	next += 1
3048	Histogram['Residuals']['Nobs'] = nobs
3049	Histogram['Residuals']['sumwYo'] = sumwYo
3050	SumwYo += sumwYo
3051	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/sumwYo)100.)
3052	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
3053	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
3054	Histogram['Residuals'][phfx+'Nref'] = nobs
3055	Histogram['Residuals'][phfx+'Nrej'] = nrej
3056	Histogram['Residuals'][phfx+'Next'] = next
3057	if im:
3058	Histogram['Residuals'][phfx+'SSRf'] = 100.*sumSSdF[:maxH+1]/sumSSFo[:maxH+1]
3059	Histogram['Residuals'][phfx+'SSRf^2'] = 100.*sumSSdF2[:maxH+1]/sumSSFo2[:maxH+1]
3060	Histogram['Residuals'][phfx+'SSNref'] = SSnobs[:maxH+1]
3061	Histogram['Residuals']['SSwR'] = np.sqrt(sumSSwdf2[:maxH+1]/sumSSwYo[:maxH+1])*100.
3062	Nobs += nobs
3063	Nrej += nrej
3064	Next += next
3065	if dlg:
3066	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
3067	M = np.concatenate((M,wtFactor*df))
3068	# end of HKLF processing
3069	Histograms['sumwYo'] = SumwYo
3070	Histograms['Nobs'] = Nobs
3071	Histograms['Nrej'] = Nrej
3072	Histograms['Next'] = Next
3073	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
3074	if dlg:
3075	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
3076	if not GoOn:
3077	parmDict['saved values'] = values
3078	dlg.Destroy()
3079	raise G2obj.G2Exception('User abort') #Abort!!
3080	pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
3081	if len(pVals):
3082	pSum = np.sum(pWtpVals*2)
3083	for name in pWsum:
3084	if pWsum[name]:
3085	print ' Penalty function for %8s = %12.5g'%(name,pWsum[name])
3086	print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
3087	Nobs += len(pVals)
3088	M = np.concatenate((M,np.sqrt(pWt)*pVals))
3089	return M
3090

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

Download in other formats: