source: branch/2frame/GSASIIstrMath.py @ 2962

Last change on this file since 2962 was 2962, checked in by toby, 4 years ago

fix bug in derivative for excluded regions

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