source: trunk/GSASIIstrMath.py @ 3383

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