source: trunk/GSASIIstrMath.py @ 2501

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

fix (again!) mag moment drawings
clean up testDeriv

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