source: trunk/GSASIIstrMath.py @ 3374

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

further mag derivatives - getting closer

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