source: trunk/GSASIIstrMath.py @ 2097

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

name fix in AtomTypeSelect? - remove ( & )
split StructureFactorDeriv? (SS & non SS versions) into twin & nontwin/PWDR versions
non twin versions all test OK for refl blocks, twin versions need work (use single refl version for now). SS versions do single refl for now.

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