source: trunk/GSASIIstrMath.py @ 2512

Last change on this file since 2512 was 2512, checked in by vondreele, 5 years ago

remove 'constraint' and 'restraint' from Consraints & Restraints tabs
work on automatic nucl-mag constraints
allow selection of possible mag atoms to use in new mag phase
fix issues with chemical composition restraints

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