Context Navigation

source: trunk/GSASIIstrMath.py @ 2024

Last change on this file since 2024 was 2024, checked in by vondreele, 8 years ago
dmin in indexing of unit cells is now taken directly from the powder pattern limits rather than some saved value which caused issues with TOF data. minor revisions to SS structure factor
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 157.2 KB

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

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

Download in other formats: