Context Navigation

source: trunk/GSASIIstrMath.py @ 2038

Last change on this file since 2038 was 2038, checked in by vondreele, 8 years ago
add betaij2Uij to G2lattice revisions to SS structure factor calcs. Use hklt proj to hkl is several places fxn works for thiourea derivs OK except X,Y,Zcos modulations; no Uijsin/cos derivatives yet adj scaling of 4D charge flip maps convert betaij vals from Jana2K files to Uij start on SS read phase from cif added a hklt F import (might vanish)
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 159.9 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMath - structure math routines
4	-----------------------------------------
5	'''
6	########### SVN repository information ###################
7	# $Date: 2015-10-30 21:02:02 +0000 (Fri, 30 Oct 2015) $
8	# $Author: vondreele $
9	# $Revision: 2038 $
10	# $URL: trunk/GSASIIstrMath.py $
11	# $Id: GSASIIstrMath.py 2038 2015-10-30 21:02:02Z vondreele $
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: 2038 $")
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]2,axis=0)+np.sum(fbs[:,:,0]2,axis=0) #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(fas2,axis=0)+np.sum(fbs2,axis=0)
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]+fbs[1]
901	SB = fbs[0]+fas[1]
902	if nTwin > 1:
903	dFdfr[iref] = [2.TwMask[it](fas[0][it]dfadfr[0][it]+fas[1][it]dfadfr[1][it]+fbs[0][it]dfbdfr[0][it]+fbs[1][it]dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
904	dFdx[iref] = [2.TwMask[it](fas[0][it]dfadx[it][0]+fas[1][it]dfadx[it][1]+fbs[0][it]dfbdx[it][0]+fbs[1][it]dfbdx[it][1]) for it in range(nTwin)]
905	dFdui[iref] = [2.TwMask[it](fas[0][it]dfadui[it][0]+fas[1][it]dfadui[it][1]+fbs[0][it]dfbdui[it][0]+fbs[1][it]dfbdui[it][1]) for it in range(nTwin)]
906	dFdua[iref] = [2.TwMask[it](fas[0][it]dfadua[it][0]+fas[1][it]dfadua[it][1]+fbs[0][it]dfbdua[it][0]+fbs[1][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.(fas[0]dfadfr[0]+fas[1]dfadfr[1]+fbs[0]dfbdfr[0]+fbs[1]dfbdfr[1])Mdata/len(Uniq)
911	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]dfadx[1]+fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
912	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]dfadui[1]+fbs[0]dfbdui[0]+fbs[1]*dfbdui[1])
913	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]dfadua[1]+fbs[0]dfbdua[0]+fbs[1]*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]2,axis=0)+np.sum(fbs[:,:,0]2,axis=0) #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(fas2,axis=0)+np.sum(fbs2,axis=0) #sum of squares
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	PhiP = np.inner(HP,SGT)
1199	if SGInv: #if centro - expand HKL sets
1200	Uniq = np.vstack((Uniq,-Uniq))
1201	Phi = np.hstack((Phi,-Phi))
1202	UniqP = np.vstack((UniqP,-UniqP))
1203	PhiP = np.hstack((PhiP,-PhiP))
1204	# GSASIIpath.IPyBreak()
1205	phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1206	sinp = np.sin(phase)
1207	cosp = np.cos(phase)
1208	occ = Mdata*Fdata/Uniq.shape[0]
1209	biso = -SQfactor*Uisodata[:,nxs]
1210	Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T #ops x atoms
1211	HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1) #ops x atoms
1212	Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1213	Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
1214	Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms
1215	Tcorr = np.reshape(Tiso,Tuij.shape)TuijMdata*Fdata/Uniq.shape[0] #ops x atoms
1216	fot = (FF+FP-Bab)*Tcorr #ops x atoms
1217	fotp = FPP*Tcorr #ops x atoms
1218	GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,UniqP,Hij,FSSdata,XSSdata,USSdata,Mast)
1219	# derivs are: ops x atoms x waves x 1,3,or 6 parms as [real,imag] parts
1220	fa = np.array([((FF+FP).T-Bab).TcospTcorr,-FlackFPPsinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
1221	fb = np.array([((FF+FP).T-Bab).TsinpTcorr,FlackFPPcosp*Tcorr])
1222	fag = faGfpuA[0]-fbGfpuA[1]
1223	fbg = fbGfpuA[0]+faGfpuA[1]
1224
1225	fas = np.sum(np.sum(fag,axis=1),axis=1)
1226	fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1227	fax = np.array([-fotsinp,-fotpcosp]) #positions; 2 x twin x ops x atoms
1228	fbx = np.array([fotcosp,-fotpsinp])
1229	fax = faxGfpuA[0]-fbxGfpuA[1]
1230	fbx = fbxGfpuA[0]+faxGfpuA[1]
1231	#sum below is over Uniq
1232	dfadfr = np.sum(fag/occ,axis=1) #Fdata != 0 ever avoids /0. problem
1233	dfbdfr = np.sum(fbg/occ,axis=1) #Fdata != 0 avoids /0. problem
1234	dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1235	dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1236	dfadui = np.sum(-SQfactor*fag,axis=1)
1237	dfbdui = np.sum(-SQfactor*fbg,axis=1)
1238	if nTwin > 1:
1239	dfadx = np.array([np.sum(twopiUniq[it,:,:3]np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1240	dfbdx = np.array([np.sum(twopiUniq[it,:,:3]np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1241	dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1242	dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1243	# array(nTwin,2,nAtom,3) & array(2,nTwin,nAtom,6)
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	dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]dGdu[1][nxs,nxs,:,:,:,:],axis=1)
1247	dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]dGdu[1][nxs,nxs,:,:,:,:],axis=1)
1248	# array (2,nAtom,wave,3)
1249	else:
1250	dfadx = np.sum(twopiUniq[:,:3]np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)
1251	dfbdx = np.sum(twopiUniq[:,:3]np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)
1252	dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)
1253	dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)
1254	# array(2,nAtom,3) & array(2,nAtom,6)
1255	dfadGx = np.sum(fa[:,:,:,nxs,nxs]dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]dGdx[1][nxs,:,:,:,:],axis=1)
1256	dfbdGx = np.sum(fb[:,:,:,nxs,nxs]dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]dGdx[1][nxs,:,:,:,:],axis=1)
1257	dfadGu = np.sum(fa[:,:,:,nxs,nxs]dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]dGdu[1][nxs,:,:,:,:],axis=1)
1258	dfbdGu = np.sum(fb[:,:,:,nxs,nxs]dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]dGdu[1][nxs,:,:,:,:],axis=1)
1259	# array (2,nAtom,wave,3)
1260	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!
1261	# GSASIIpath.IPyBreak()
1262	if not SGData['SGInv'] and len(TwinLaw) == 1: #Flack derivative
1263	dfadfl = np.sum(-FPPTcorrsinp)
1264	dfbdfl = np.sum(FPPTcorrcosp)
1265	else:
1266	dfadfl = 1.0
1267	dfbdfl = 1.0
1268	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!
1269	if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
1270	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)+ \
1271	2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
1272	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])+ \
1273	2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
1274	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])+ \
1275	2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
1276	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])+ \
1277	2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
1278	dFdGx[iref] = 2.(fas[0]dfadGx[0]+fas[1]*dfadGx[1])+ \
1279	2.(fbs[0]dfbdGx[0]+fbs[1]*dfbdGx[1])
1280	dFdGu[iref] = 2.(fas[0]dfadGu[0]+fas[1]*dfadGu[1])+ \
1281	2.(fbs[0]dfbdGu[0]+fbs[1]*dfbdGu[1])
1282	else:
1283	SA = fas[0]-fbs[1]
1284	SB = fbs[0]+fas[1]
1285	if nTwin > 1:
1286	dFdfr[iref] = [2.TwMask[it](fas[0][it]dfadfr[0][it]+fas[1][it]dfadfr[1][it]+fbs[0][it]dfbdfr[0][it]+fbs[1][it]dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
1287	dFdx[iref] = [2.TwMask[it](fas[0][it]dfadx[it][0]+fas[1][it]dfadx[it][1]+fbs[0][it]dfbdx[it][0]+fbs[1][it]dfbdx[it][1]) for it in range(nTwin)]
1288	dFdui[iref] = [2.TwMask[it](fas[0][it]dfadui[it][0]+fas[1][it]dfadui[it][1]+fbs[0][it]dfbdui[it][0]+fbs[1][it]dfbdui[it][1]) for it in range(nTwin)]
1289	dFdua[iref] = [2.TwMask[it](fas[0][it]dfadua[it][0]+fas[1][it]dfadua[it][1]+fbs[0][it]dfbdua[it][0]+fbs[1][it]dfbdua[it][1]) for it in range(nTwin)]
1290	dFdGx[iref] = [2.TwMask[it](fas[0][it]dfadGx[0]+fas[1][it]dfadGx[1]+fbs[0][it]dfbdGx[0]+fbs[1][it]dfbdGx[1]) for it in range(nTwin)]
1291	dFdGu[iref] = [2.TwMask[it](fas[0][it]dfadGu[0]+fas[1][it]dfadGu[1]+fbs[0][it]dfbdGu[0]+fbs[1][it]dfbdGu[1]) for it in range(nTwin)]
1292	dFdtw[iref] = np.sum(TwMaskfas,axis=0)2+np.sum(TwMaskfbs,axis=0)**2
1293
1294	else: #these are good for no twin single crystals
1295	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1]+fbs[0]dfbdfr[0]+fbs[1]dfbdfr[1])Mdata/len(Uniq) #array(nRef,nAtom)
1296	dFdx[iref] = 2.(fas[0]dfadx[0]+fbs[1]dfbdx[1]+fbs[0]dfbdx[0]+fas[1]*dfadx[1]) #array(nRef,nAtom,3)
1297	dFdui[iref] = 2.(fas[0]dfadui[0]+fbs[1]dfbdui[1]+fbs[0]dfbdui[0]+fas[1]*dfadui[1]) #array(nRef,nAtom)
1298	dFdua[iref] = 2.(fas[0]dfadua[0]+fbs[1]dfbdua[1]+fbs[0]dfbdua[0]+fas[1]*dfadua[1]) #array(nRef,nAtom,6)
1299	dFdfl[iref] = -SAdfadfl-SBdfbdfl #array(nRef,)
1300	dFdGx[iref] = 2.(fas[0]dfadGx[0]+fas[1]dfadGx[1]+fbs[0]dfbdGx[0]+fbs[1]*dfbdGx[1]) #array(nRef,natom,nwave,6)
1301	dFdGu[iref] = 2.(fas[0]dfadGu[0]+fas[1]dfadGu[1]+fbs[0]dfbdGu[0]+fbs[1]*dfbdGu[1]) #array(nRef,natom,nwave,6)
1302	dFdbab[iref] = 2.fas[0]np.array([np.sum(dfadbadBabdA),np.sum(-dfadbaparmDict[phfx+'BabA']SQfactordBabdA)]).T+ \
1303	2.fbs[0]np.array([np.sum(dfbdbadBabdA),np.sum(-dfbdbaparmDict[phfx+'BabA']SQfactordBabdA)]).T
1304	#loop over atoms - each dict entry is list of derivatives for all the reflections
1305	if not iref%100 :
1306	print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
1307	for i in range(len(Mdata)): #loop over atoms
1308	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1309	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1310	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1311	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1312	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1313	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1314	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1315	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1316	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1317	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1318	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1319	for j in range(XSSdata.shape[1]):
1320	dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j)] = dFdGx.T[0][j][i]
1321	dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j)] = dFdGx.T[1][j][i]
1322	dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j)] = dFdGx.T[2][j][i]
1323	dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j)] = dFdGx.T[3][j][i]
1324	dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j)] = dFdGx.T[4][j][i]
1325	dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j)] = dFdGx.T[5][j][i]
1326	for j in range(USSdata.shape[1]):
1327	dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1328	dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1329	dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1330	dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
1331	dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
1332	dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
1333	dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1334	dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1335	dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1336	dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
1337	dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
1338	dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
1339
1340	dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1341	dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1342	return dFdvDict
1343
1344	def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1345	''' Single crystal extinction function; returns extinction & derivative
1346	'''
1347	extCor = 1.0
1348	dervDict = {}
1349	dervCor = 1.0
1350	if calcControls[phfx+'EType'] != 'None':
1351	SQ = 1/(4.ref[4+im]*2)
1352	if 'C' in parmDict[hfx+'Type']:
1353	cos2T = 1.0-2.SQparmDict[hfx+'Lam']**2 #cos(2theta)
1354	else: #'T'
1355	cos2T = 1.0-2.SQref[12+im]**2 #cos(2theta)
1356	if 'SXC' in parmDict[hfx+'Type']:
1357	AV = 7.9406e5/parmDict[pfx+'Vol']**2
1358	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1359	P12 = (calcControls[phfx+'Cos2TM']+cos2T4)/(calcControls[phfx+'Cos2TM']+cos2T2)
1360	PLZ = AVP12ref[9+im]parmDict[hfx+'Lam']*2
1361	elif 'SNT' in parmDict[hfx+'Type']:
1362	AV = 1.e7/parmDict[pfx+'Vol']**2
1363	PL = SQ
1364	PLZ = AVref[9+im]ref[12+im]**2
1365	elif 'SNC' in parmDict[hfx+'Type']:
1366	AV = 1.e7/parmDict[pfx+'Vol']**2
1367	PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1368	PLZ = AVref[9+im]parmDict[hfx+'Lam']**2
1369
1370	if 'Primary' in calcControls[phfx+'EType']:
1371	PLZ *= 1.5
1372	else:
1373	if 'C' in parmDict[hfx+'Type']:
1374	PLZ *= calcControls[phfx+'Tbar']
1375	else: #'T'
1376	PLZ *= ref[13+im] #t-bar
1377	if 'Primary' in calcControls[phfx+'EType']:
1378	PLZ *= 1.5
1379	PSIG = parmDict[phfx+'Ep']
1380	elif 'I & II' in calcControls[phfx+'EType']:
1381	PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']PL/parmDict[phfx+'Eg'])*2)
1382	elif 'Type II' in calcControls[phfx+'EType']:
1383	PSIG = parmDict[phfx+'Es']
1384	else: # 'Secondary Type I'
1385	PSIG = parmDict[phfx+'Eg']/PL
1386
1387	AG = 0.58+0.48cos2T+0.24cos2T**2
1388	AL = 0.025+0.285*cos2T
1389	BG = 0.02-0.025*cos2T
1390	BL = 0.15-0.2(0.75-cos2T)*2
1391	if cos2T < 0.:
1392	BL = -0.45*cos2T
1393	CG = 2.
1394	CL = 2.
1395	PF = PLZ*PSIG
1396
1397	if 'Gaussian' in calcControls[phfx+'EApprox']:
1398	PF4 = 1.+CGPF+AGPF*2/(1.+BGPF)
1399	extCor = np.sqrt(PF4)
1400	PF3 = 0.5(CG+2.AGPF/(1.+BGPF)-AGPF2BG/(1.+BGPF)2)/(PF4extCor)
1401	else:
1402	PF4 = 1.+CLPF+ALPF*2/(1.+BLPF)
1403	extCor = np.sqrt(PF4)
1404	PF3 = 0.5(CL+2.ALPF/(1.+BLPF)-ALPF2BL/(1.+BLPF)2)/(PF4extCor)
1405
1406	dervCor = (1.+PF)*PF3 #extinction corr for other derivatives
1407	if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1408	dervDict[phfx+'Ep'] = -ref[7+im]PLZPF3
1409	if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1410	dervDict[phfx+'Es'] = -ref[7+im]PLZPF3(PSIG/parmDict[phfx+'Es'])*3
1411	if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1412	dervDict[phfx+'Eg'] = -ref[7+im]PLZPF3(PSIG/parmDict[phfx+'Eg'])3PL**2
1413
1414	return 1./extCor,dervDict,dervCor
1415
1416	def Dict2Values(parmdict, varylist):
1417	'''Use before call to leastsq to setup list of values for the parameters
1418	in parmdict, as selected by key in varylist'''
1419	return [parmdict[key] for key in varylist]
1420
1421	def Values2Dict(parmdict, varylist, values):
1422	''' Use after call to leastsq to update the parameter dictionary with
1423	values corresponding to keys in varylist'''
1424	parmdict.update(zip(varylist,values))
1425
1426	def GetNewCellParms(parmDict,varyList):
1427	'Needs a doc string'
1428	newCellDict = {}
1429	Anames = ['A'+str(i) for i in range(6)]
1430	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1431	for item in varyList:
1432	keys = item.split(':')
1433	if keys[2] in Ddict:
1434	key = keys[0]+'::'+Ddict[keys[2]] #key is e.g. '0::A0'
1435	parm = keys[0]+'::'+keys[2] #parm is e.g. '0::D11'
1436	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1437	return newCellDict # is e.g. {'0::D11':A0-D11}
1438
1439	def ApplyXYZshifts(parmDict,varyList):
1440	'''
1441	takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1442
1443	:param dict parmDict: parameter dictionary
1444	:param list varyList: list of variables (not used!)
1445	:returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1446
1447	'''
1448	newAtomDict = {}
1449	for item in parmDict:
1450	if 'dA' in item:
1451	parm = ''.join(item.split('d'))
1452	parmDict[parm] += parmDict[item]
1453	newAtomDict[item] = [parm,parmDict[parm]]
1454	return newAtomDict
1455
1456	def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1457	'Spherical harmonics texture'
1458	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1459	if 'T' in calcControls[hfx+'histType']:
1460	tth = parmDict[hfx+'2-theta']
1461	else:
1462	tth = refl[5+im]
1463	odfCor = 1.0
1464	H = refl[:3]
1465	cell = G2lat.Gmat2cell(g)
1466	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1467	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1468	phi,beta = G2lat.CrsAng(H,cell,SGData)
1469	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1470	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1471	for item in SHnames:
1472	L,M,N = eval(item.strip('C'))
1473	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1474	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1475	Lnorm = G2lat.Lnorm(L)
1476	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1477	return odfCor
1478
1479	def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1480	'Spherical harmonics texture derivatives'
1481	if 'T' in calcControls[hfx+'histType']:
1482	tth = parmDict[hfx+'2-theta']
1483	else:
1484	tth = refl[5+im]
1485	IFCoup = 'Bragg' in calcControls[hfx+'instType']
1486	odfCor = 1.0
1487	dFdODF = {}
1488	dFdSA = [0,0,0]
1489	H = refl[:3]
1490	cell = G2lat.Gmat2cell(g)
1491	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1492	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1493	phi,beta = G2lat.CrsAng(H,cell,SGData)
1494	psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1495	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1496	for item in SHnames:
1497	L,M,N = eval(item.strip('C'))
1498	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1499	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1500	Lnorm = G2lat.Lnorm(L)
1501	odfCor += parmDict[pfx+item]LnormKcl*Ksl
1502	dFdODF[pfx+item] = LnormKclKsl
1503	for i in range(3):
1504	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
1505	return odfCor,dFdODF,dFdSA
1506
1507	def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1508	'spherical harmonics preferred orientation (cylindrical symmetry only)'
1509	if 'T' in calcControls[hfx+'histType']:
1510	tth = parmDict[hfx+'2-theta']
1511	else:
1512	tth = refl[5+im]
1513	odfCor = 1.0
1514	H = refl[:3]
1515	cell = G2lat.Gmat2cell(g)
1516	Sangls = [0.,0.,0.]
1517	if 'Bragg' in calcControls[hfx+'instType']:
1518	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1519	IFCoup = True
1520	else:
1521	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1522	IFCoup = False
1523	phi,beta = G2lat.CrsAng(H,cell,SGData)
1524	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1525	SHnames = calcControls[phfx+'SHnames']
1526	for item in SHnames:
1527	L,N = eval(item.strip('C'))
1528	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1529	Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1530	Lnorm = G2lat.Lnorm(L)
1531	odfCor += parmDict[phfx+item]LnormKcl*Ksl
1532	return np.squeeze(odfCor)
1533
1534	def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1535	'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1536	if 'T' in calcControls[hfx+'histType']:
1537	tth = parmDict[hfx+'2-theta']
1538	else:
1539	tth = refl[5+im]
1540	odfCor = 1.0
1541	dFdODF = {}
1542	H = refl[:3]
1543	cell = G2lat.Gmat2cell(g)
1544	Sangls = [0.,0.,0.]
1545	if 'Bragg' in calcControls[hfx+'instType']:
1546	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1547	IFCoup = True
1548	else:
1549	Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1550	IFCoup = False
1551	phi,beta = G2lat.CrsAng(H,cell,SGData)
1552	psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1553	SHnames = calcControls[phfx+'SHnames']
1554	for item in SHnames:
1555	L,N = eval(item.strip('C'))
1556	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1557	Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1558	Lnorm = G2lat.Lnorm(L)
1559	odfCor += parmDict[phfx+item]LnormKcl*Ksl
1560	dFdODF[phfx+item] = KclKslLnorm
1561	return odfCor,dFdODF
1562
1563	def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1564	'March-Dollase preferred orientation correction'
1565	POcorr = 1.0
1566	MD = parmDict[phfx+'MD']
1567	if MD != 1.0:
1568	MDAxis = calcControls[phfx+'MDAxis']
1569	sumMD = 0
1570	for H in uniq:
1571	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1572	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1573	sumMD += A**3
1574	POcorr = sumMD/len(uniq)
1575	return POcorr
1576
1577	def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1578	'Needs a doc string'
1579	POcorr = 1.0
1580	POderv = {}
1581	if calcControls[phfx+'poType'] == 'MD':
1582	MD = parmDict[phfx+'MD']
1583	MDAxis = calcControls[phfx+'MDAxis']
1584	sumMD = 0
1585	sumdMD = 0
1586	for H in uniq:
1587	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1588	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
1589	sumMD += A**3
1590	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
1591	POcorr = sumMD/len(uniq)
1592	POderv[phfx+'MD'] = sumdMD/len(uniq)
1593	else: #spherical harmonics
1594	if calcControls[phfx+'SHord']:
1595	POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1596	return POcorr,POderv
1597
1598	def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1599	'Needs a doc string'
1600	if 'Debye' in calcControls[hfx+'instType']:
1601	if 'T' in calcControls[hfx+'histType']:
1602	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1603	else:
1604	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1605	else:
1606	return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1607
1608	def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1609	'Needs a doc string'
1610	if 'Debye' in calcControls[hfx+'instType']:
1611	if 'T' in calcControls[hfx+'histType']:
1612	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1613	else:
1614	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1615	else:
1616	return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1617
1618	def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1619	'Needs a doc string'
1620	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1621	pi2 = np.sqrt(2./np.pi)
1622	if 'T' in calcControls[hfx+'histType']:
1623	sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1624	wave = refl[14+im]
1625	else: #'C'W
1626	sth2 = sind(refl[5+im]/2.)**2
1627	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1628	c2th = 1.-2.0*sth2
1629	flv2 = refl[9+im](wave/parmDict[pfx+'Vol'])*2
1630	if 'X' in calcControls[hfx+'histType']:
1631	flv2 = 0.079411(1.0+c2th**2)/2.0
1632	xfac = flv2*parmDict[phfx+'Extinction']
1633	exb = 1.0
1634	if xfac > -1.:
1635	exb = 1./np.sqrt(1.+xfac)
1636	exl = 1.0
1637	if 0 < xfac <= 1.:
1638	xn = np.array([xfac**(i+1) for i in range(6)])
1639	exl += np.sum(xn*coef)
1640	elif xfac > 1.:
1641	xfac2 = 1./np.sqrt(xfac)
1642	exl = pi2(1.-0.125/xfac)xfac2
1643	return exbsth2+exl(1.-sth2)
1644
1645	def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1646	'Needs a doc string'
1647	coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1648	pi2 = np.sqrt(2./np.pi)
1649	if 'T' in calcControls[hfx+'histType']:
1650	sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1651	wave = refl[14+im]
1652	else: #'C'W
1653	sth2 = sind(refl[5+im]/2.)**2
1654	wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1655	c2th = 1.-2.0*sth2
1656	flv2 = refl[9+im](wave/parmDict[pfx+'Vol'])*2
1657	if 'X' in calcControls[hfx+'histType']:
1658	flv2 = 0.079411(1.0+c2th**2)/2.0
1659	xfac = flv2*parmDict[phfx+'Extinction']
1660	dbde = -500.*flv2
1661	if xfac > -1.:
1662	dbde = -0.5flv2/np.sqrt(1.+xfac)*3
1663	dlde = 0.
1664	if 0 < xfac <= 1.:
1665	xn = np.array([iflv2xfac**i for i in [1,2,3,4,5,6]])
1666	dlde = np.sum(xn*coef)
1667	elif xfac > 1.:
1668	xfac2 = 1./np.sqrt(xfac)
1669	dlde = flv2pi2xfac2(-1./xfac+0.375/xfac*2)
1670
1671	return dbdesth2+dlde(1.-sth2)
1672
1673	def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1674	'Needs a doc string' #need powder extinction!
1675	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3+im] #scale*multiplicity
1676	if 'X' in parmDict[hfx+'Type']:
1677	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1678	POcorr = 1.0
1679	if pfx+'SHorder' in parmDict: #generalized spherical harmonics texture - takes precidence
1680	POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1681	elif calcControls[phfx+'poType'] == 'MD': #March-Dollase
1682	POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1683	elif calcControls[phfx+'SHord']: #cylindrical spherical harmonics
1684	POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1685	Icorr *= POcorr
1686	AbsCorr = 1.0
1687	AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1688	Icorr *= AbsCorr
1689	ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1690	Icorr *= ExtCorr
1691	return Icorr,POcorr,AbsCorr,ExtCorr
1692
1693	def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1694	'Needs a doc string' #need powder extinction derivs!
1695	dIdsh = 1./parmDict[hfx+'Scale']
1696	dIdsp = 1./parmDict[phfx+'Scale']
1697	if 'X' in parmDict[hfx+'Type']:
1698	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1699	dIdPola /= pola
1700	else: #'N'
1701	dIdPola = 0.0
1702	dFdODF = {}
1703	dFdSA = [0,0,0]
1704	dIdPO = {}
1705	if pfx+'SHorder' in parmDict:
1706	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1707	for iSH in dFdODF:
1708	dFdODF[iSH] /= odfCor
1709	for i in range(3):
1710	dFdSA[i] /= odfCor
1711	elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1712	POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1713	for iPO in dIdPO:
1714	dIdPO[iPO] /= POcorr
1715	if 'T' in parmDict[hfx+'Type']:
1716	dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1717	dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im] #/ext corr
1718	else:
1719	dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1720	dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im] #/ext corr
1721	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1722
1723	def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1724	'Needs a doc string'
1725	if 'C' in calcControls[hfx+'histType']: #All checked & OK
1726	costh = cosd(refl[5+im]/2.)
1727	#crystallite size
1728	if calcControls[phfx+'SizeType'] == 'isotropic':
1729	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1730	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1731	H = np.array(refl[:3])
1732	P = np.array(calcControls[phfx+'SizeAxis'])
1733	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1734	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a']*costh)
1735	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1736	else: #ellipsoidal crystallites
1737	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1738	H = np.array(refl[:3])
1739	lenR = G2pwd.ellipseSize(H,Sij,GB)
1740	Sgam = 1.8wave/(np.picosth*lenR)
1741	#microstrain
1742	if calcControls[phfx+'MustrainType'] == 'isotropic':
1743	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5+im]/2.)/np.pi
1744	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1745	H = np.array(refl[:3])
1746	P = np.array(calcControls[phfx+'MustrainAxis'])
1747	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1748	Si = parmDict[phfx+'Mustrain;i']
1749	Sa = parmDict[phfx+'Mustrain;a']
1750	Mgam = 0.018SiSatand(refl[5+im]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
1751	else: #generalized - P.W. Stephens model
1752	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1753	Sum = 0
1754	for i,strm in enumerate(Strms):
1755	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1756	Mgam = 0.018refl[4+im]2tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1757	elif 'T' in calcControls[hfx+'histType']: #All checked & OK
1758	#crystallite size
1759	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
1760	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/parmDict[phfx+'Size;i']
1761	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
1762	H = np.array(refl[:3])
1763	P = np.array(calcControls[phfx+'SizeAxis'])
1764	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1765	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]*2/(parmDict[phfx+'Size;i']parmDict[phfx+'Size;a'])
1766	Sgam = np.sqrt((sinPparmDict[phfx+'Size;a'])*2+(cosPparmDict[phfx+'Size;i'])**2)
1767	else: #ellipsoidal crystallites #OK
1768	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1769	H = np.array(refl[:3])
1770	lenR = G2pwd.ellipseSize(H,Sij,GB)
1771	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/lenR
1772	#microstrain
1773	if calcControls[phfx+'MustrainType'] == 'isotropic': #OK
1774	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]*parmDict[phfx+'Mustrain;i']
1775	elif calcControls[phfx+'MustrainType'] == 'uniaxial': #OK
1776	H = np.array(refl[:3])
1777	P = np.array(calcControls[phfx+'MustrainAxis'])
1778	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1779	Si = parmDict[phfx+'Mustrain;i']
1780	Sa = parmDict[phfx+'Mustrain;a']
1781	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]SiSa/np.sqrt((SicosP)2+(SasinP)**2)
1782	else: #generalized - P.W. Stephens model OK
1783	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1784	Sum = 0
1785	for i,strm in enumerate(Strms):
1786	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1787	Mgam = 1.e-6parmDict[hfx+'difC']np.sqrt(Sum)refl[4+im]*3
1788
1789	gam = SgamparmDict[phfx+'Size;mx']+MgamparmDict[phfx+'Mustrain;mx']
1790	sig = (Sgam(1.-parmDict[phfx+'Size;mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain;mx']))**2
1791	sig /= ateln2
1792	return sig,gam
1793
1794	def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1795	'Needs a doc string'
1796	gamDict = {}
1797	sigDict = {}
1798	if 'C' in calcControls[hfx+'histType']: #All checked & OK
1799	costh = cosd(refl[5+im]/2.)
1800	tanth = tand(refl[5+im]/2.)
1801	#crystallite size derivatives
1802	if calcControls[phfx+'SizeType'] == 'isotropic':
1803	Sgam = 1.8wave/(np.piparmDict[phfx+'Size;i']*costh)
1804	gamDict[phfx+'Size;i'] = -1.8waveparmDict[phfx+'Size;mx']/(np.pi*costh)
1805	sigDict[phfx+'Size;i'] = -3.6Sgamwave(1.-parmDict[phfx+'Size;mx'])2/(np.picosth*ateln2)
1806	elif calcControls[phfx+'SizeType'] == 'uniaxial':
1807	H = np.array(refl[:3])
1808	P = np.array(calcControls[phfx+'SizeAxis'])
1809	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1810	Si = parmDict[phfx+'Size;i']
1811	Sa = parmDict[phfx+'Size;a']
1812	gami = 1.8wave/(costhnp.piSiSa)
1813	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1814	Sgam = gami*sqtrm
1815	dsi = gamiSicosP**2/sqtrm-Sgam/Si
1816	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
1817	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1818	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1819	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1820	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1821	else: #ellipsoidal crystallites
1822	const = 1.8wave/(np.picosth)
1823	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1824	H = np.array(refl[:3])
1825	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1826	Sgam = const/lenR
1827	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1828	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1829	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1830	gamDict[phfx+'Size;mx'] = Sgam
1831	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2
1832
1833	#microstrain derivatives
1834	if calcControls[phfx+'MustrainType'] == 'isotropic':
1835	Mgam = 0.018parmDict[phfx+'Mustrain;i']tand(refl[5+im]/2.)/np.pi
1836	gamDict[phfx+'Mustrain;i'] = 0.018tanthparmDict[phfx+'Mustrain;mx']/np.pi
1837	sigDict[phfx+'Mustrain;i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain;mx'])2/(np.piateln2)
1838	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1839	H = np.array(refl[:3])
1840	P = np.array(calcControls[phfx+'MustrainAxis'])
1841	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1842	Si = parmDict[phfx+'Mustrain;i']
1843	Sa = parmDict[phfx+'Mustrain;a']
1844	gami = 0.018SiSa*tanth/np.pi
1845	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1846	Mgam = gami/sqtrm
1847	dsi = -gamiSicosP2/sqtrm3
1848	dsa = -gamiSasinP2/sqtrm3
1849	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1850	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1851	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1852	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1853	else: #generalized - P.W. Stephens model
1854	const = 0.018refl[4+im]2tanth/np.pi
1855	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1856	Sum = 0
1857	for i,strm in enumerate(Strms):
1858	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1859	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1860	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
1861	Mgam = const*np.sqrt(Sum)
1862	for i in range(len(Strms)):
1863	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1864	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
1865	gamDict[phfx+'Mustrain;mx'] = Mgam
1866	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1867	else: #'T'OF - All checked & OK
1868	if calcControls[phfx+'SizeType'] == 'isotropic': #OK
1869	Sgam = 1.e-4parmDict[hfx+'difC']refl[4+im]**2/parmDict[phfx+'Size;i']
1870	gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1871	sigDict[phfx+'Size;i'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])*2/(ateln2parmDict[phfx+'Size;i'])
1872	elif calcControls[phfx+'SizeType'] == 'uniaxial': #OK
1873	const = 1.e-4parmDict[hfx+'difC']refl[4+im]**2
1874	H = np.array(refl[:3])
1875	P = np.array(calcControls[phfx+'SizeAxis'])
1876	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1877	Si = parmDict[phfx+'Size;i']
1878	Sa = parmDict[phfx+'Size;a']
1879	gami = const/(Si*Sa)
1880	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
1881	Sgam = gami*sqtrm
1882	dsi = gamiSicosP**2/sqtrm-Sgam/Si
1883	dsa = gamiSasinP**2/sqtrm-Sgam/Sa
1884	gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1885	gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1886	sigDict[phfx+'Size;i'] = 2.dsiSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1887	sigDict[phfx+'Size;a'] = 2.dsaSgam(1.-parmDict[phfx+'Size;mx'])*2/ateln2
1888	else: #OK ellipsoidal crystallites
1889	const = 1.e-4parmDict[hfx+'difC']refl[4+im]**2
1890	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1891	H = np.array(refl[:3])
1892	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1893	Sgam = const/lenR
1894	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1895	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size;mx']
1896	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size;mx'])*2/ateln2
1897	gamDict[phfx+'Size;mx'] = Sgam #OK
1898	sigDict[phfx+'Size;mx'] = -2.Sgam2(1.-parmDict[phfx+'Size;mx'])/ateln2 #OK
1899
1900	#microstrain derivatives
1901	if calcControls[phfx+'MustrainType'] == 'isotropic':
1902	Mgam = 1.e-6parmDict[hfx+'difC']refl[4+im]*parmDict[phfx+'Mustrain;i']
1903	gamDict[phfx+'Mustrain;i'] = 1.e-6refl[4+im]parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx'] #OK
1904	sigDict[phfx+'Mustrain;i'] = 2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])*2/(ateln2parmDict[phfx+'Mustrain;i'])
1905	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1906	H = np.array(refl[:3])
1907	P = np.array(calcControls[phfx+'MustrainAxis'])
1908	cosP,sinP = G2lat.CosSinAngle(H,P,G)
1909	Si = parmDict[phfx+'Mustrain;i']
1910	Sa = parmDict[phfx+'Mustrain;a']
1911	gami = 1.e-6parmDict[hfx+'difC']refl[4+im]SiSa
1912	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
1913	Mgam = gami/sqtrm
1914	dsi = -gamiSicosP2/sqtrm3
1915	dsa = -gamiSasinP2/sqtrm3
1916	gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1917	gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1918	sigDict[phfx+'Mustrain;i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1919	sigDict[phfx+'Mustrain;a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain;mx'])*2/ateln2
1920	else: #generalized - P.W. Stephens model OK
1921	pwrs = calcControls[phfx+'MuPwrs']
1922	Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1923	const = 1.e-6parmDict[hfx+'difC']refl[4+im]**3
1924	Sum = 0
1925	for i,strm in enumerate(Strms):
1926	Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1927	gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1928	sigDict[phfx+'Mustrain:'+str(i)] = strm(1.-parmDict[phfx+'Mustrain;mx'])*2
1929	Mgam = const*np.sqrt(Sum)
1930	for i in range(len(Strms)):
1931	gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1932	sigDict[phfx+'Mustrain:'+str(i)] = const*2/ateln2
1933	gamDict[phfx+'Mustrain;mx'] = Mgam
1934	sigDict[phfx+'Mustrain;mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1935
1936	return sigDict,gamDict
1937
1938	def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1939	'Needs a doc string'
1940	if im:
1941	h,k,l,m = refl[:4]
1942	vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1943	d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1944	else:
1945	h,k,l = refl[:3]
1946	d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1947	refl[4+im] = d
1948	if 'C' in calcControls[hfx+'histType']:
1949	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
1950	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1951	if 'Bragg' in calcControls[hfx+'instType']:
1952	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1953	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
1954	else: #Debye-Scherrer - simple but maybe not right
1955	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1956	elif 'T' in calcControls[hfx+'histType']:
1957	pos = parmDict[hfx+'difC']d+parmDict[hfx+'difA']d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1958	#do I need sample position effects - maybe?
1959	return pos
1960
1961	def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1962	'Needs a doc string'
1963	dpr = 180./np.pi
1964	if im:
1965	h,k,l,m = refl[:4]
1966	vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1967	dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1968	else:
1969	m = 0
1970	h,k,l = refl[:3]
1971	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1972	dst = np.sqrt(dstsq)
1973	dsp = 1./dst
1974	if 'C' in calcControls[hfx+'histType']:
1975	pos = refl[5+im]-parmDict[hfx+'Zero']
1976	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
1977	dpdw = const*dst
1978	dpdA = np.array([h2,k2,l*2,hk,hl,kl])constwave/(2.0*dst)
1979	dpdZ = 1.0
1980	dpdV = np.array([2.hA[0]+kA[3]+lA[4],2kA[1]+hA[3]+lA[5],
1981	2lA[2]+hA[4]+kA[5]])mconstwave/(2.0dst)
1982	shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
1983	if 'Bragg' in calcControls[hfx+'instType']:
1984	dpdSh = -4.shftcosd(pos/2.0)
1985	dpdTr = -shftsind(pos)100.0
1986	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1987	else: #Debye-Scherrer - simple but maybe not right
1988	dpdXd = -shft*cosd(pos)
1989	dpdYd = -shft*sind(pos)
1990	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1991	elif 'T' in calcControls[hfx+'histType']:
1992	dpdA = -np.array([h2,k2,l*2,hk,hl,kl])parmDict[hfx+'difC']dsp**3/2.
1993	dpdZ = 1.0
1994	dpdDC = dsp
1995	dpdDA = dsp**2
1996	dpdDB = 1./dsp
1997	dpdV = np.array([2.hA[0]+kA[3]+lA[4],2kA[1]+hA[3]+lA[5],
1998	2lA[2]+hA[4]+kA[5]])mparmDict[hfx+'difC']dsp**3/2.
1999	return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2000
2001	def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2002	'Needs a doc string'
2003	laue = SGData['SGLaue']
2004	uniq = SGData['SGUniq']
2005	h,k,l = refl[:3]
2006	if laue in ['m3','m3m']:
2007	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
2008	refl[4+im]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
2009	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2010	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
2011	elif laue in ['3R','3mR']:
2012	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
2013	elif laue in ['4/m','4/mmm']:
2014	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
2015	elif laue in ['mmm']:
2016	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2017	elif laue in ['2/m']:
2018	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2019	if uniq == 'a':
2020	Dij += parmDict[phfx+'D23']kl
2021	elif uniq == 'b':
2022	Dij += parmDict[phfx+'D13']hl
2023	elif uniq == 'c':
2024	Dij += parmDict[phfx+'D12']hk
2025	else:
2026	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
2027	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
2028	if 'C' in calcControls[hfx+'histType']:
2029	return -180.Dijrefl[4+im]*2tand(refl[5+im]/2.0)/np.pi
2030	else:
2031	return -DijparmDict[hfx+'difC']refl[4+im]**2/2.
2032
2033	def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2034	'Needs a doc string'
2035	laue = SGData['SGLaue']
2036	uniq = SGData['SGUniq']
2037	h,k,l = refl[:3]
2038	if laue in ['m3','m3m']:
2039	dDijDict = {phfx+'D11':h2+k2+l**2,
2040	phfx+'eA':refl[4+im]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
2041	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2042	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
2043	elif laue in ['3R','3mR']:
2044	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
2045	elif laue in ['4/m','4/mmm']:
2046	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
2047	elif laue in ['mmm']:
2048	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2049	elif laue in ['2/m']:
2050	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2051	if uniq == 'a':
2052	dDijDict[phfx+'D23'] = k*l
2053	elif uniq == 'b':
2054	dDijDict[phfx+'D13'] = h*l
2055	elif uniq == 'c':
2056	dDijDict[phfx+'D12'] = h*k
2057	else:
2058	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
2059	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
2060	if 'C' in calcControls[hfx+'histType']:
2061	for item in dDijDict:
2062	dDijDict[item] = 180.0refl[4+im]*2tand(refl[5+im]/2.0)/np.pi
2063	else:
2064	for item in dDijDict:
2065	dDijDict[item] = -parmDict[hfx+'difC']refl[4+im]**3/2.
2066	return dDijDict
2067
2068	def GetDij(phfx,SGData,parmDict):
2069	HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2070	return G2spc.HStrainVals(HSvals,SGData)
2071
2072	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2073	'Needs a doc string'
2074	histoList = Histograms.keys()
2075	histoList.sort()
2076	for histogram in histoList:
2077	if 'PWDR' in histogram[:4]:
2078	Histogram = Histograms[histogram]
2079	hId = Histogram['hId']
2080	hfx = ':%d:'%(hId)
2081	Limits = calcControls[hfx+'Limits']
2082	if 'C' in calcControls[hfx+'histType']:
2083	shl = max(parmDict[hfx+'SH/L'],0.0005)
2084	Ka2 = False
2085	kRatio = 0.0
2086	if hfx+'Lam1' in parmDict.keys():
2087	Ka2 = True
2088	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2089	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2090	x,y,w,yc,yb,yd = Histogram['Data']
2091	xB = np.searchsorted(x,Limits[0])
2092	xF = np.searchsorted(x,Limits[1])
2093	ymb = np.array(y-yb)
2094	ymb = np.where(ymb,ymb,1.0)
2095	ycmb = np.array(yc-yb)
2096	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
2097	refLists = Histogram['Reflection Lists']
2098	for phase in refLists:
2099	if phase not in Phases: #skips deleted or renamed phases silently!
2100	continue
2101	Phase = Phases[phase]
2102	im = 0
2103	if Phase['General']['Type'] in ['modulated','magnetic']:
2104	im = 1
2105	pId = Phase['pId']
2106	phfx = '%d:%d:'%(pId,hId)
2107	refDict = refLists[phase]
2108	sumFo = 0.0
2109	sumdF = 0.0
2110	sumFosq = 0.0
2111	sumdFsq = 0.0
2112	sumInt = 0.0
2113	for refl in refDict['RefList']:
2114	if 'C' in calcControls[hfx+'histType']:
2115	yp = np.zeros_like(yb)
2116	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2117	iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2118	iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2119	iFin2 = iFin
2120	if not iBeg+iFin: #peak below low limit - skip peak
2121	continue
2122	elif not iBeg-iFin: #peak above high limit - done
2123	break
2124	elif iBeg < iFin:
2125	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
2126	sumInt += refl[11+im]*refl[9+im]
2127	if Ka2:
2128	pos2 = refl[5+im]+lamRatiotand(refl[5+im]/2.0) # + 360/pi Dlam/lam * tan(th)
2129	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2130	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
2131	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
2132	if iFin2 > iBeg2:
2133	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
2134	sumInt += refl[11+im]refl[9+im]kRatio
2135	refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[11+im](1.+kRatio)),0.0))
2136
2137	elif 'T' in calcControls[hfx+'histType']:
2138	yp = np.zeros_like(yb)
2139	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2140	iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2141	iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2142	if iBeg < iFin:
2143	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
2144	refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
2145	sumInt += refl[11+im]*refl[9+im]
2146	Fo = np.sqrt(np.abs(refl[8+im]))
2147	Fc = np.sqrt(np.abs(refl[9]+im))
2148	sumFo += Fo
2149	sumFosq += refl[8+im]**2
2150	sumdF += np.abs(Fo-Fc)
2151	sumdFsq += (refl[8+im]-refl[9+im])**2
2152	if sumFo:
2153	Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2154	Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2155	else:
2156	Histogram['Residuals'][phfx+'Rf'] = 100.
2157	Histogram['Residuals'][phfx+'Rf^2'] = 100.
2158	Histogram['Residuals'][phfx+'sumInt'] = sumInt
2159	Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
2160	Histogram['Residuals']['hId'] = hId
2161	elif 'HKLF' in histogram[:4]:
2162	Histogram = Histograms[histogram]
2163	Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
2164
2165	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2166	'Needs a doc string'
2167
2168	def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
2169	U = parmDict[hfx+'U']
2170	V = parmDict[hfx+'V']
2171	W = parmDict[hfx+'W']
2172	X = parmDict[hfx+'X']
2173	Y = parmDict[hfx+'Y']
2174	tanPos = tand(refl[5+im]/2.0)
2175	Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2176	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
2177	sig = max(0.001,sig)
2178	gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam #save peak gamma
2179	gam = max(0.001,gam)
2180	return sig,gam
2181
2182	def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
2183	sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']refl[4+im]*2+ \
2184	parmDict[hfx+'sig-2']refl[4+im]4+parmDict[hfx+'sig-q']/refl[4+im]*2
2185	gam = parmDict[hfx+'X']refl[4+im]+parmDict[hfx+'Y']refl[4+im]**2
2186	Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2187	sig += Ssig
2188	gam += Sgam
2189	return sig,gam
2190
2191	def GetReflAlpBet(refl,im,hfx,parmDict):
2192	alp = parmDict[hfx+'alpha']/refl[4+im]
2193	bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]4+parmDict[hfx+'beta-q']/refl[4+im]2
2194	return alp,bet
2195
2196	hId = Histogram['hId']
2197	hfx = ':%d:'%(hId)
2198	bakType = calcControls[hfx+'bakType']
2199	yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2200	yc = np.zeros_like(yb)
2201	cw = np.diff(x)
2202	cw = np.append(cw,cw[-1])
2203
2204	if 'C' in calcControls[hfx+'histType']:
2205	shl = max(parmDict[hfx+'SH/L'],0.002)
2206	Ka2 = False
2207	if hfx+'Lam1' in parmDict.keys():
2208	wave = parmDict[hfx+'Lam1']
2209	Ka2 = True
2210	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2211	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2212	else:
2213	wave = parmDict[hfx+'Lam']
2214	for phase in Histogram['Reflection Lists']:
2215	refDict = Histogram['Reflection Lists'][phase]
2216	if phase not in Phases: #skips deleted or renamed phases silently!
2217	continue
2218	Phase = Phases[phase]
2219	pId = Phase['pId']
2220	pfx = '%d::'%(pId)
2221	phfx = '%d:%d:'%(pId,hId)
2222	hfx = ':%d:'%(hId)
2223	SGData = Phase['General']['SGData']
2224	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2225	im = 0
2226	if Phase['General']['Type'] in ['modulated','magnetic']:
2227	SSGData = Phase['General']['SSGData']
2228	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2229	im = 1 #offset in SS reflection list
2230	#??
2231	Dij = GetDij(phfx,SGData,parmDict)
2232	A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2233	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2234	if np.any(np.diag(G)<0.):
2235	raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
2236	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2237	Vst = np.sqrt(nl.det(G)) #V*
2238	if not Phase['General'].get('doPawley'):
2239	time0 = time.time()
2240	if im:
2241	SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2242	else:
2243	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2244	# print 'sf calc time: %.3fs'%(time.time()-time0)
2245	time0 = time.time()
2246	badPeak = False
2247	for iref,refl in enumerate(refDict['RefList']):
2248	if 'C' in calcControls[hfx+'histType']:
2249	if im:
2250	h,k,l,m = refl[:4]
2251	else:
2252	h,k,l = refl[:3]
2253	Uniq = np.inner(refl[:3],SGMT)
2254	refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict) #corrected reflection position
2255	Lorenz = 1./(2.sind(refl[5+im]/2.)2cosd(refl[5+im]/2.)) #Lorentz correction
2256	refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict) #peak sig & gam
2257	refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2258	refl[11+im] = VstLorenz
2259
2260	if Phase['General'].get('doPawley'):
2261	try:
2262	if im:
2263	pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2264	else:
2265	pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2266	refl[9+im] = parmDict[pInd]
2267	except KeyError:
2268	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2269	continue
2270	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2271	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2272	iFin = np.searchsorted(x,refl[5+im]+fmax)
2273	if not iBeg+iFin: #peak below low limit - skip peak
2274	continue
2275	elif not iBeg-iFin: #peak above high limit - done
2276	break
2277	elif iBeg > iFin: #bad peak coeff - skip
2278	badPeak = True
2279	continue
2280	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
2281	if Ka2:
2282	pos2 = refl[5+im]+lamRatiotand(refl[5+im]/2.0) # + 360/pi Dlam/lam * tan(th)
2283	Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2284	iBeg = np.searchsorted(x,pos2-fmin)
2285	iFin = np.searchsorted(x,pos2+fmax)
2286	if not iBeg+iFin: #peak below low limit - skip peak
2287	continue
2288	elif not iBeg-iFin: #peak above high limit - done
2289	return yc,yb
2290	elif iBeg > iFin: #bad peak coeff - skip
2291	continue
2292	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
2293	elif 'T' in calcControls[hfx+'histType']:
2294	h,k,l = refl[:3]
2295	Uniq = np.inner(refl[:3],SGMT)
2296	refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict) #corrected reflection position
2297	Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)refl[4+im]*4 #TOF Lorentz correction
2298	# refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict) #apply hydrostatic strain shift
2299	refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict) #peak sig & gam
2300	refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2301	refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2302	refl[11+im] = VstLorenz
2303	if Phase['General'].get('doPawley'):
2304	try:
2305	if im:
2306	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2307	else:
2308	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2309	refl[9+im] = parmDict[pInd]
2310	except KeyError:
2311	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2312	continue
2313	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2314	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2315	iFin = np.searchsorted(x,refl[5+im]+fmax)
2316	if not iBeg+iFin: #peak below low limit - skip peak
2317	continue
2318	elif not iBeg-iFin: #peak above high limit - done
2319	break
2320	elif iBeg > iFin: #bad peak coeff - skip
2321	badPeak = True
2322	continue
2323	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]
2324	# print 'profile calc time: %.3fs'%(time.time()-time0)
2325	if badPeak:
2326	print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped'
2327	return yc,yb
2328
2329	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2330	'Needs a doc string'
2331
2332	def cellVaryDerv(pfx,SGData,dpdA):
2333	if SGData['SGLaue'] in ['-1',]:
2334	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2335	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2336	elif SGData['SGLaue'] in ['2/m',]:
2337	if SGData['SGUniq'] == 'a':
2338	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2339	elif SGData['SGUniq'] == 'b':
2340	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2341	else:
2342	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2343	elif SGData['SGLaue'] in ['mmm',]:
2344	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2345	elif SGData['SGLaue'] in ['4/m','4/mmm']:
2346	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2347	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2348	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2349	elif SGData['SGLaue'] in ['3R', '3mR']:
2350	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
2351	elif SGData['SGLaue'] in ['m3m','m3']:
2352	return [[pfx+'A0',dpdA[0]]]
2353
2354	# create a list of dependent variables and set up a dictionary to hold their derivatives
2355	dependentVars = G2mv.GetDependentVars()
2356	depDerivDict = {}
2357	for j in dependentVars:
2358	depDerivDict[j] = np.zeros(shape=(len(x)))
2359	#print 'dependent vars',dependentVars
2360	lenX = len(x)
2361	hId = Histogram['hId']
2362	hfx = ':%d:'%(hId)
2363	bakType = calcControls[hfx+'bakType']
2364	dMdv = np.zeros(shape=(len(varylist),len(x)))
2365	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2366	if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2367	bBpos =varylist.index(hfx+'Back;0')
2368	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2369	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2370	for name in varylist:
2371	if 'Debye' in name:
2372	id = int(name.split(';')[-1])
2373	parm = name[:int(name.rindex(';'))]
2374	ip = names.index(parm)
2375	dMdv[varylist.index(name)] = dMddb[3*id+ip]
2376	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2377	for name in varylist:
2378	if 'BkPk' in name:
2379	parm,id = name.split(';')
2380	id = int(id)
2381	if parm in names:
2382	ip = names.index(parm)
2383	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2384	cw = np.diff(x)
2385	cw = np.append(cw,cw[-1])
2386	Ka2 = False #also for TOF!
2387	if 'C' in calcControls[hfx+'histType']:
2388	shl = max(parmDict[hfx+'SH/L'],0.002)
2389	if hfx+'Lam1' in parmDict.keys():
2390	wave = parmDict[hfx+'Lam1']
2391	Ka2 = True
2392	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2393	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2394	else:
2395	wave = parmDict[hfx+'Lam']
2396	for phase in Histogram['Reflection Lists']:
2397	refDict = Histogram['Reflection Lists'][phase]
2398	if phase not in Phases: #skips deleted or renamed phases silently!
2399	continue
2400	Phase = Phases[phase]
2401	SGData = Phase['General']['SGData']
2402	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2403	im = 0
2404	if Phase['General']['Type'] in ['modulated','magnetic']:
2405	SSGData = Phase['General']['SSGData']
2406	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2407	im = 1 #offset in SS reflection list
2408	#??
2409	pId = Phase['pId']
2410	pfx = '%d::'%(pId)
2411	phfx = '%d:%d:'%(pId,hId)
2412	Dij = GetDij(phfx,SGData,parmDict)
2413	A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2414	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2415	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2416	if not Phase['General'].get('doPawley'):
2417	time0 = time.time()
2418	if im:
2419	dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2420	else:
2421	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2422	# print 'sf-derv time %.3fs'%(time.time()-time0)
2423	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2424	time0 = time.time()
2425	for iref,refl in enumerate(refDict['RefList']):
2426	if im:
2427	h,k,l,m = refl[:4]
2428	else:
2429	h,k,l = refl[:3]
2430	Uniq = np.inner(refl[:3],SGMT)
2431	if 'T' in calcControls[hfx+'histType']:
2432	wave = refl[14+im]
2433	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2434	if 'C' in calcControls[hfx+'histType']: #CW powder
2435	Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2436	else: #'T'OF
2437	Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2438	iBeg = np.searchsorted(x,refl[5+im]-fmin)
2439	iFin = np.searchsorted(x,refl[5+im]+fmax)
2440	if not iBeg+iFin: #peak below low limit - skip peak
2441	continue
2442	elif not iBeg-iFin: #peak above high limit - done
2443	break
2444	pos = refl[5+im]
2445	if 'C' in calcControls[hfx+'histType']:
2446	tanth = tand(pos/2.0)
2447	costh = cosd(pos/2.0)
2448	lenBF = iFin-iBeg
2449	dMdpk = np.zeros(shape=(6,lenBF))
2450	dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2451	for i in range(5):
2452	dMdpk[i] += 100.cw[iBeg:iFin]refl[11+im]refl[9+im]dMdipk[i]
2453	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2454	if Ka2:
2455	pos2 = refl[5+im]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
2456	iBeg2 = np.searchsorted(x,pos2-fmin)
2457	iFin2 = np.searchsorted(x,pos2+fmax)
2458	if iBeg2-iFin2:
2459	lenBF2 = iFin2-iBeg2
2460	dMdpk2 = np.zeros(shape=(6,lenBF2))
2461	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2462	for i in range(5):
2463	dMdpk2[i] = 100.cw[iBeg2:iFin2]refl[11+im]refl[9+im]kRatio*dMdipk2[i]
2464	dMdpk2[5] = 100.cw[iBeg2:iFin2]refl[11+im]*dMdipk2[0]
2465	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2466	else: #'T'OF
2467	lenBF = iFin-iBeg
2468	if lenBF < 0: #bad peak coeff
2469	break
2470	dMdpk = np.zeros(shape=(6,lenBF))
2471	dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2472	for i in range(6):
2473	dMdpk[i] += refl[11+im]refl[9+im]dMdipk[i] #cw[iBeg:iFin]*
2474	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
2475	if Phase['General'].get('doPawley'):
2476	dMdpw = np.zeros(len(x))
2477	try:
2478	if im:
2479	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2480	else:
2481	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2482	idx = varylist.index(pIdx)
2483	dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2484	if Ka2: #not for TOF either
2485	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2486	dMdv[idx] = dMdpw
2487	except: # ValueError:
2488	pass
2489	if 'C' in calcControls[hfx+'histType']:
2490	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2491	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2492	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2493	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2494	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2495	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2496	hfx+'DisplaceY':[dpdY,'pos'],}
2497	if 'Bragg' in calcControls[hfx+'instType']:
2498	names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2499	hfx+'SurfRoughB':[dFdAb[1],'int'],})
2500	else:
2501	names.update({hfx+'Absorption':[dFdAb,'int'],})
2502	else: #'T'OF
2503	dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2504	names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2505	hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2506	hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2507	hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2508	hfx+'beta-q':[1./refl[4+im]2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]2,'sig'],
2509	hfx+'sig-2':[refl[4+im]4,'sig'],hfx+'sig-q':[1./refl[4+im]2,'sig'],
2510	hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2511	for name in names:
2512	item = names[name]
2513	if name in varylist:
2514	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2515	if Ka2 and iFin2-iBeg2:
2516	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2517	elif name in dependentVars:
2518	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2519	if Ka2 and iFin2-iBeg2:
2520	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2521	for iPO in dIdPO:
2522	if iPO in varylist:
2523	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2524	if Ka2 and iFin2-iBeg2:
2525	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2526	elif iPO in dependentVars:
2527	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2528	if Ka2 and iFin2-iBeg2:
2529	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2530	for i,name in enumerate(['omega','chi','phi']):
2531	aname = pfx+'SH '+name
2532	if aname in varylist:
2533	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2534	if Ka2 and iFin2-iBeg2:
2535	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2536	elif aname in dependentVars:
2537	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2538	if Ka2 and iFin2-iBeg2:
2539	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2540	for iSH in dFdODF:
2541	if iSH in varylist:
2542	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2543	if Ka2 and iFin2-iBeg2:
2544	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2545	elif iSH in dependentVars:
2546	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2547	if Ka2 and iFin2-iBeg2:
2548	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2549	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2550	for name,dpdA in cellDervNames:
2551	if name in varylist:
2552	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2553	if Ka2 and iFin2-iBeg2:
2554	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2555	elif name in dependentVars:
2556	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2557	if Ka2 and iFin2-iBeg2:
2558	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2559	dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2560	for name in dDijDict:
2561	if name in varylist:
2562	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2563	if Ka2 and iFin2-iBeg2:
2564	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2565	elif name in dependentVars:
2566	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2567	if Ka2 and iFin2-iBeg2:
2568	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2569	for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2570	if name in varylist:
2571	dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2572	if Ka2 and iFin2-iBeg2:
2573	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2574	elif name in dependentVars:
2575	depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2576	if Ka2 and iFin2-iBeg2:
2577	depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2578	if 'C' in calcControls[hfx+'histType']:
2579	sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2580	else: #'T'OF
2581	sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2582	for name in gamDict:
2583	if name in varylist:
2584	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2585	if Ka2 and iFin2-iBeg2:
2586	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2587	elif name in dependentVars:
2588	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2589	if Ka2 and iFin2-iBeg2:
2590	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2591	for name in sigDict:
2592	if name in varylist:
2593	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2594	if Ka2 and iFin2-iBeg2:
2595	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2596	elif name in dependentVars:
2597	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2598	if Ka2 and iFin2-iBeg2:
2599	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2600	for name in ['BabA','BabU']:
2601	if refl[9+im]:
2602	if phfx+name in varylist:
2603	dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict['int']/refl[9+im]
2604	if Ka2 and iFin2-iBeg2:
2605	dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict2['int']/refl[9+im]
2606	elif phfx+name in dependentVars:
2607	depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict['int']/refl[9+im]
2608	if Ka2 and iFin2-iBeg2:
2609	depDerivDict[phfx+name][iBeg2:iFin2] += parmDict[phfx+'Scale']dFdvDict[phfx+name][iref]dervDict2['int']/refl[9+im]
2610	if not Phase['General'].get('doPawley'):
2611	#do atom derivatives - for RB,F,X & U so far
2612	corr = dervDict['int']/refl[9+im]
2613	if Ka2 and iFin2-iBeg2:
2614	corr2 = dervDict2['int']/refl[9+im]
2615	for name in varylist+dependentVars:
2616	if '::RBV;' in name:
2617	pass
2618	else:
2619	try:
2620	aname = name.split(pfx)[1][:2]
2621	if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2622	except IndexError:
2623	continue
2624	if name in varylist:
2625	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2626	if Ka2 and iFin2-iBeg2:
2627	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2628	elif name in dependentVars:
2629	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2630	if Ka2 and iFin2-iBeg2:
2631	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2632	# print 'profile derv time: %.3fs'%(time.time()-time0)
2633	# now process derivatives in constraints
2634	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2635	return dMdv
2636
2637	def UserRejectHKL(ref,im,userReject):
2638	if ref[5+im]/ref[6+im] < userReject['minF/sig']:
2639	return False
2640	elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
2641	return False
2642	elif ref[11+im] < userReject['MinExt']:
2643	return False
2644	elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
2645	return False
2646	return True
2647
2648	def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2649	'''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2650	model (M) with respect to all parameters. Independent and dependant dM/dp arrays
2651	are returned to either dervRefine or HessRefine.
2652
2653	:returns:
2654	'''
2655	nobs = Histogram['Residuals']['Nobs']
2656	hId = Histogram['hId']
2657	hfx = ':%d:'%(hId)
2658	pfx = '%d::'%(Phase['pId'])
2659	phfx = '%d:%d:'%(Phase['pId'],hId)
2660	SGData = Phase['General']['SGData']
2661	im = 0
2662	if Phase['General']['Type'] in ['modulated','magnetic']:
2663	SSGData = Phase['General']['SSGData']
2664	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2665	im = 1 #offset in SS reflection list
2666	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2667	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2668	refDict = Histogram['Data']
2669	if parmDict[phfx+'Scale'] < 0.:
2670	parmDict[phfx+'Scale'] = .001
2671	if im:
2672	dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2673	else:
2674	dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2675	ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2676	dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2677	dependentVars = G2mv.GetDependentVars()
2678	depDerivDict = {}
2679	for j in dependentVars:
2680	depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2681	wdf = np.zeros(len(refDict['RefList']))
2682	if calcControls['F**2']:
2683	for iref,ref in enumerate(refDict['RefList']):
2684	if ref[6+im] > 0:
2685	dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2686	w = 1.0/ref[6+im]
2687	if ref[3+im] > 0:
2688	wdf[iref] = w*(ref[5+im]-ref[7+im])
2689	for j,var in enumerate(varylist):
2690	if var in dFdvDict:
2691	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2692	for var in dependentVars:
2693	if var in dFdvDict:
2694	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2695	if phfx+'Scale' in varylist:
2696	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2697	elif phfx+'Scale' in dependentVars:
2698	depDerivDict[phfx+'Scale'][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2699	for item in ['Ep','Es','Eg']:
2700	if phfx+item in varylist and phfx+item in dervDict:
2701	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] #OK
2702	elif phfx+item in dependentVars and phfx+item in dervDict:
2703	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im] #OK
2704	for item in ['BabA','BabU']:
2705	if phfx+item in varylist:
2706	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2707	elif phfx+item in dependentVars:
2708	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2709	else: #F refinement
2710	for iref,ref in enumerate(refDict['RefList']):
2711	if ref[5+im] > 0.:
2712	dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2713	Fo = np.sqrt(ref[5+im])
2714	Fc = np.sqrt(ref[7+im])
2715	w = 1.0/ref[6+im]
2716	if ref[3+im] > 0:
2717	wdf[iref] = 2.0Fcw*(Fo-Fc)
2718	for j,var in enumerate(varylist):
2719	if var in dFdvDict:
2720	dMdvh[j][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2721	for var in dependentVars:
2722	if var in dFdvDict:
2723	depDerivDict[var][iref] = wdFdvDict[var][iref]parmDict[phfx+'Scale']*ref[11+im]
2724	if phfx+'Scale' in varylist:
2725	dMdvh[varylist.index(phfx+'Scale')][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2726	elif phfx+'Scale' in dependentVars:
2727	depDerivDict[phfx+'Scale'][iref] = wref[7+im]ref[11+im]/parmDict[phfx+'Scale'] #OK
2728	for item in ['Ep','Es','Eg']: #OK!
2729	if phfx+item in varylist and phfx+item in dervDict:
2730	dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]
2731	elif phfx+item in dependentVars and phfx+item in dervDict:
2732	depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2733	for item in ['BabA','BabU']:
2734	if phfx+item in varylist:
2735	dMdvh[varylist.index(phfx+item)][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2736	elif phfx+item in dependentVars:
2737	depDerivDict[phfx+item][iref] = wdFdvDict[pfx+item][iref]parmDict[phfx+'Scale']*ref[11+im]
2738	return dMdvh,depDerivDict,wdf
2739
2740
2741	def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2742	'''Loop over histograms and compute derivatives of the fitting
2743	model (M) with respect to all parameters. Results are returned in
2744	a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2745	n is the number of parameters and m is the number of data
2746	points. This can exceed memory when m gets large. This routine is
2747	used when refinement derivatives are selected as "analtytic
2748	Jacobian" in Controls.
2749
2750	:returns: Jacobian numpy.array dMdv for all histograms concatinated
2751	'''
2752	parmDict.update(zip(varylist,values))
2753	G2mv.Dict2Map(parmDict,varylist)
2754	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2755	nvar = len(varylist)
2756	dMdv = np.empty(0)
2757	histoList = Histograms.keys()
2758	histoList.sort()
2759	for histogram in histoList:
2760	if 'PWDR' in histogram[:4]:
2761	Histogram = Histograms[histogram]
2762	hId = Histogram['hId']
2763	hfx = ':%d:'%(hId)
2764	wtFactor = calcControls[hfx+'wtFactor']
2765	Limits = calcControls[hfx+'Limits']
2766	x,y,w,yc,yb,yd = Histogram['Data']
2767	xB = np.searchsorted(x,Limits[0])
2768	xF = np.searchsorted(x,Limits[1])
2769	dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2770	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2771	elif 'HKLF' in histogram[:4]:
2772	Histogram = Histograms[histogram]
2773	phase = Histogram['Reflection Lists']
2774	Phase = Phases[phase]
2775	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2776	hfx = ':%d:'%(Histogram['hId'])
2777	wtFactor = calcControls[hfx+'wtFactor']
2778	# now process derivatives in constraints
2779	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2780	else:
2781	continue #skip non-histogram entries
2782	if len(dMdv):
2783	dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2784	else:
2785	dMdv = np.sqrt(wtFactor)*dMdvh
2786
2787	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2788	if np.any(pVals):
2789	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2790	dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2791
2792	return dMdv
2793
2794	def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2795	'''Loop over histograms and compute derivatives of the fitting
2796	model (M) with respect to all parameters. For each histogram, the
2797	Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2798	number of parameters and m is the number of data points *in the
2799	histogram*. The (n by n) Hessian is computed from each Jacobian
2800	and it is returned. This routine is used when refinement
2801	derivatives are selected as "analtytic Hessian" in Controls.
2802
2803	:returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2804	'''
2805	parmDict.update(zip(varylist,values))
2806	G2mv.Dict2Map(parmDict,varylist)
2807	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2808	#fixup H atom positions here?
2809	ApplyRBModels(parmDict,Phases,rigidbodyDict) #,Update=True??
2810	nvar = len(varylist)
2811	Hess = np.empty(0)
2812	histoList = Histograms.keys()
2813	histoList.sort()
2814	for histogram in histoList:
2815	if 'PWDR' in histogram[:4]:
2816	Histogram = Histograms[histogram]
2817	hId = Histogram['hId']
2818	hfx = ':%d:'%(hId)
2819	wtFactor = calcControls[hfx+'wtFactor']
2820	Limits = calcControls[hfx+'Limits']
2821	x,y,w,yc,yb,yd = Histogram['Data']
2822	W = wtFactor*w
2823	dy = y-yc
2824	xB = np.searchsorted(x,Limits[0])
2825	xF = np.searchsorted(x,Limits[1])
2826	dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2827	varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2828	Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2829	Dy = dy[xB:xF][np.newaxis,:]
2830	dMdvh *= Wt
2831	if dlg:
2832	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2833	if len(Hess):
2834	Hess += np.inner(dMdvh,dMdvh)
2835	dMdvh = WtDy
2836	Vec += np.sum(dMdvh,axis=1)
2837	else:
2838	Hess = np.inner(dMdvh,dMdvh)
2839	dMdvh = WtDy
2840	Vec = np.sum(dMdvh,axis=1)
2841	elif 'HKLF' in histogram[:4]:
2842	Histogram = Histograms[histogram]
2843	phase = Histogram['Reflection Lists']
2844	Phase = Phases[phase]
2845	dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2846	hId = Histogram['hId']
2847	hfx = ':%d:'%(Histogram['hId'])
2848	wtFactor = calcControls[hfx+'wtFactor']
2849	# now process derivatives in constraints
2850	G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2851	# print 'matrix build time: %.3f'%(time.time()-time0)
2852
2853	if dlg:
2854	dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2855	if len(Hess):
2856	Vec += wtFactornp.sum(dMdvhwdf,axis=1)
2857	Hess += wtFactor*np.inner(dMdvh,dMdvh)
2858	else:
2859	Vec = wtFactornp.sum(dMdvhwdf,axis=1)
2860	Hess = wtFactor*np.inner(dMdvh,dMdvh)
2861	else:
2862	continue #skip non-histogram entries
2863	pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2864	if np.any(pVals):
2865	dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2866	Vec += np.sum(dpdvpWtpVals,axis=1)
2867	Hess += np.inner(dpdv*pWt,dpdv)
2868	return Vec,Hess
2869
2870	def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):
2871	'''Computes the point-by-point discrepancies between every data point in every histogram
2872	and the observed value
2873
2874	:returns: an np array of differences between observed and computed diffraction values.
2875	'''
2876	Values2Dict(parmDict, varylist, values)
2877	G2mv.Dict2Map(parmDict,varylist)
2878	Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2879	M = np.empty(0)
2880	SumwYo = 0
2881	Nobs = 0
2882	Nrej = 0
2883	Next = 0
2884	ApplyRBModels(parmDict,Phases,rigidbodyDict)
2885	#fixup Hatom positions here....
2886	histoList = Histograms.keys()
2887	histoList.sort()
2888	for histogram in histoList:
2889	if 'PWDR' in histogram[:4]:
2890	Histogram = Histograms[histogram]
2891	hId = Histogram['hId']
2892	hfx = ':%d:'%(hId)
2893	wtFactor = calcControls[hfx+'wtFactor']
2894	Limits = calcControls[hfx+'Limits']
2895	x,y,w,yc,yb,yd = Histogram['Data']
2896	yc *= 0.0 #zero full calcd profiles
2897	yb *= 0.0
2898	yd *= 0.0
2899	xB = np.searchsorted(x,Limits[0])
2900	xF = np.searchsorted(x,Limits[1])
2901	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2902	varylist,Histogram,Phases,calcControls,pawleyLookup)
2903	yc[xB:xF] += yb[xB:xF]
2904	if not np.any(y): #fill dummy data
2905	rv = st.poisson(yc[xB:xF])
2906	y[xB:xF] = rv.rvs()
2907	Z = np.ones_like(yc[xB:xF])
2908	Z[1::2] *= -1
2909	y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2910	w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2911	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2912	W = wtFactor*w
2913	wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2914	Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2915	Nobs += Histogram['Residuals']['Nobs']
2916	Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]y[xB:xF]*2)
2917	SumwYo += Histogram['Residuals']['sumwYo']
2918	Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2919	Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy*2)/Histogram['Residuals']['sumwYo'])100.)
2920	sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2921	sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF](y[xB:xF]-yb[xB:xF])*2,0.))
2922	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.))
2923	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.))
2924	Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2925	Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2926	Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2927	if dlg:
2928	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2929	M = np.concatenate((M,wdy))
2930	#end of PWDR processing
2931	elif 'HKLF' in histogram[:4]:
2932	Histogram = Histograms[histogram]
2933	Histogram['Residuals'] = {}
2934	phase = Histogram['Reflection Lists']
2935	Phase = Phases[phase]
2936	hId = Histogram['hId']
2937	hfx = ':%d:'%(hId)
2938	wtFactor = calcControls[hfx+'wtFactor']
2939	pfx = '%d::'%(Phase['pId'])
2940	phfx = '%d:%d:'%(Phase['pId'],hId)
2941	SGData = Phase['General']['SGData']
2942	im = 0
2943	if parmDict[phfx+'Scale'] < 0.:
2944	parmDict[phfx+'Scale'] = .001
2945	if Phase['General']['Type'] in ['modulated','magnetic']:
2946	SSGData = Phase['General']['SSGData']
2947	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2948	im = 1 #offset in SS reflection list
2949	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2950	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2951	refDict = Histogram['Data']
2952	time0 = time.time()
2953	if im:
2954	SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2955	else:
2956	StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2957	# print 'sf-calc time: %.3f'%(time.time()-time0)
2958	df = np.zeros(len(refDict['RefList']))
2959	sumwYo = 0
2960	sumFo = 0
2961	sumFo2 = 0
2962	sumdF = 0
2963	sumdF2 = 0
2964	if im:
2965	sumSSFo = np.zeros(10)
2966	sumSSFo2 = np.zeros(10)
2967	sumSSdF = np.zeros(10)
2968	sumSSdF2 = np.zeros(10)
2969	sumSSwYo = np.zeros(10)
2970	sumSSwdf2 = np.zeros(10)
2971	SSnobs = np.zeros(10)
2972	nobs = 0
2973	nrej = 0
2974	next = 0
2975	maxH = 0
2976	if calcControls['F**2']:
2977	for i,ref in enumerate(refDict['RefList']):
2978	if ref[6+im] > 0:
2979	ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2980	w = 1.0/ref[6+im] # 1/sig(F^2)
2981	ref[7+im] = parmDict[phfx+'Scale']ref[11+im] #correct Fc^2 for extinction
2982	ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2983	if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0)
2984	ref[3+im] = abs(ref[3+im]) #mark as allowed
2985	Fo = np.sqrt(ref[5+im])
2986	sumFo += Fo
2987	sumFo2 += ref[5+im]
2988	sumdF += abs(Fo-np.sqrt(ref[7+im]))
2989	sumdF2 += abs(ref[5+im]-ref[7+im])
2990	nobs += 1
2991	df[i] = -w*(ref[5+im]-ref[7+im])
2992	sumwYo += (wref[5+im])2 #wFo^2
2993	if im: #accumulate super lattice sums
2994	ind = int(abs(ref[3]))
2995	sumSSFo[ind] += Fo
2996	sumSSFo2[ind] += ref[5+im]
2997	sumSSdF[ind] += abs(Fo-np.sqrt(ref[7+im]))
2998	sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
2999	sumSSwYo[ind] += (wref[5+im])2 #wFo^2
3000	sumSSwdf2[ind] += df[i]**2
3001	SSnobs[ind] += 1
3002	maxH = max(maxH,ind)
3003	else:
3004	if ref[3+im]:
3005	ref[3+im] = -abs(ref[3+im]) #mark as rejected
3006	nrej += 1
3007	else: #sp.gp.extinct
3008	next += 1
3009	else:
3010	for i,ref in enumerate(refDict['RefList']):
3011	if ref[5+im] > 0.:
3012	ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
3013	ref[7+im] = parmDict[phfx+'Scale']ref[11+im] #correct Fc^2 for extinction
3014	ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
3015	Fo = np.sqrt(ref[5+im])
3016	Fc = np.sqrt(ref[7+im])
3017	w = 2.0*Fo/ref[6+im] # 1/sig(F)?
3018	if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0)
3019	ref[3+im] = abs(ref[3+im]) #mark as allowed
3020	sumFo += Fo
3021	sumFo2 += ref[5+im]
3022	sumdF += abs(Fo-Fc)
3023	sumdF2 += abs(ref[5+im]-ref[7+im])
3024	nobs += 1
3025	df[i] = -w*(Fo-Fc)
3026	sumwYo += (wFo)*2
3027	if im:
3028	ind = int(abs(ref[3]))
3029	sumSSFo[ind] += Fo
3030	sumSSFo2[ind] += ref[5+im]
3031	sumSSdF[ind] += abs(Fo-Fc)
3032	sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
3033	sumSSwYo[ind] += (wFo)*2
3034	sumSSwdf2[ind] += df[i]**2
3035	SSnobs[ind] += 1
3036	maxH = max(maxH,ind)
3037	else:
3038	if ref[3+im]:
3039	ref[3+im] = -abs(ref[3+im]) #mark as rejected
3040	nrej += 1
3041	else: #sp.gp.extinct
3042	next += 1
3043	Histogram['Residuals']['Nobs'] = nobs
3044	Histogram['Residuals']['sumwYo'] = sumwYo
3045	SumwYo += sumwYo
3046	Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df*2)/sumwYo)100.)
3047	Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
3048	Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
3049	Histogram['Residuals'][phfx+'Nref'] = nobs
3050	Histogram['Residuals'][phfx+'Nrej'] = nrej
3051	Histogram['Residuals'][phfx+'Next'] = next
3052	if im:
3053	Histogram['Residuals'][phfx+'SSRf'] = 100.*sumSSdF[:maxH+1]/sumSSFo[:maxH+1]
3054	Histogram['Residuals'][phfx+'SSRf^2'] = 100.*sumSSdF2[:maxH+1]/sumSSFo2[:maxH+1]
3055	Histogram['Residuals'][phfx+'SSNref'] = SSnobs[:maxH+1]
3056	Histogram['Residuals']['SSwR'] = np.sqrt(sumSSwdf2[:maxH+1]/sumSSwYo[:maxH+1])*100.
3057	Nobs += nobs
3058	Nrej += nrej
3059	Next += next
3060	if dlg:
3061	dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
3062	M = np.concatenate((M,wtFactor*df))
3063	# end of HKLF processing
3064	Histograms['sumwYo'] = SumwYo
3065	Histograms['Nobs'] = Nobs
3066	Histograms['Nrej'] = Nrej
3067	Histograms['Next'] = Next
3068	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
3069	if dlg:
3070	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
3071	if not GoOn:
3072	parmDict['saved values'] = values
3073	dlg.Destroy()
3074	raise G2obj.G2Exception('User abort') #Abort!!
3075	pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
3076	if len(pVals):
3077	pSum = np.sum(pWtpVals*2)
3078	for name in pWsum:
3079	if pWsum[name]:
3080	print ' Penalty function for %8s = %12.5g'%(name,pWsum[name])
3081	print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
3082	Nobs += len(pVals)
3083	M = np.concatenate((M,np.sqrt(pWt)*pVals))
3084	return M
3085

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

Download in other formats: