source: trunk/GSASIIstrMath.py @ 3754

Last change on this file since 3754 was 3754, checked in by vondreele, 3 years ago

begin work on incommensurate mag. structure factors
fix problem of magnetic reflection extinctions
replace analytical derivatives of magnetic moments with numerical ones using deltM = 0.0001

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