source: trunk/GSASIIstrMath.py @ 3340

Last change on this file since 3340 was 3340, checked in by vondreele, 6 years ago

make testDeriv py3 compatible
fix problem with Reload draw atoms - mag moments
add |Mag| to lst output of magnetic moments
correct errors in mag moment structure factor & derivatives

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