Context Navigation

source: trunk/GSASIIstrMath.py @ 3383

Last change on this file since 3383 was 3383, checked in by vondreele, 4 years ago

Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 245.1 KB

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