source: trunk/GSASIIstrMath.py @ 3341

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

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