source: trunk/GSASIIstrMath.py @ 3595

Last change on this file since 3595 was 3595, checked in by vondreele, 5 years ago

fix singular matrix problem for magnetic atoms with zero moment.

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