source: trunk/GSASIIstrMath.py @ 2534

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

fix exp importer for magnetic structures - mishandled U's

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