source: trunk/GSASIIstrMath.py @ 3320

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

implement new delt/sig subplot for PWDR & SASD plots - nvoked with 'w' key
make new MagStructureFactor2 routine - removed magnetism stuff from StructureFactor2

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