source: trunk/GSASIIstrMath.py @ 3353

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

remove call to paint after selecting all atoms. (paint cleared selection!)
fix perspective issue in powder waterfall plots

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