source: trunk/GSASIIstrMath.py @ 2481

Last change on this file since 2481 was 2481, checked in by vondreele, 7 years ago

force Transform to delete nonmagnetic atoms if phase made magnetic & add 'mag' to new phase name
fix TOF cosine background function
magnetic structure refinement works (in numeric mode only)
magnetic structure factor calculation correct

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