source: trunk/GSASIIstrMath.py @ 2096

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

a bit more blocked derivative work - Babinet parms now OK.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 176.0 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2015-12-16 21:31:54 +0000 (Wed, 16 Dec 2015) $
8# $Author: vondreele $
9# $Revision: 2096 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 2096 2015-12-16 21:31:54Z 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: 2096 $")
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    input:
651   
652    :param dict refDict: where
653        'RefList' list where each ref = h,k,l,m,d,...
654        'FF' dict of form factors - filed in below
655    :param np.array G:      reciprocal metric tensor
656    :param str pfx:    phase id string
657    :param dict SGData: space group info. dictionary output from SpcGroup
658    :param dict calcControls:
659    :param dict ParmDict:
660
661    '''       
662    phfx = pfx.split(':')[0]+hfx
663    ast = np.sqrt(np.diag(G))
664    Mast = twopisq*np.multiply.outer(ast,ast)
665    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
666    SGT = np.array([ops[1] for ops in SGData['SGOps']])
667    FFtables = calcControls['FFtables']
668    BLtables = calcControls['BLtables']
669    Flack = 1.0
670    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
671        Flack = 1.-2.*parmDict[phfx+'Flack']
672    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
673    TwDict = refDict.get('TwDict',{})           
674    if 'S' in calcControls[hfx+'histType']:
675        NTL = calcControls[phfx+'NTL']
676        NM = calcControls[phfx+'TwinNMN']+1
677        TwinLaw = calcControls[phfx+'TwinLaw']
678        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
679        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
680    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
681    FF = np.zeros(len(Tdata))
682    if 'NC' in calcControls[hfx+'histType']:
683        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
684    elif 'X' in calcControls[hfx+'histType']:
685        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
686        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
687    Uij = np.array(G2lat.U6toUij(Uijdata))
688    bij = Mast*Uij.T
689    blkSize = 100       #no. of reflections in a block - size seems optimal
690    nRef = refDict['RefList'].shape[0]
691    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
692        if 'N' in calcControls[hfx+'histType']:
693            dat = G2el.getBLvalues(BLtables)
694            refDict['FF']['El'] = dat.keys()
695            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
696        else:       #'X'
697            dat = G2el.getFFvalues(FFtables,0.)
698            refDict['FF']['El'] = dat.keys()
699            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
700            for iref,ref in enumerate(refDict['RefList']):
701                SQ = 1./(2.*ref[4])**2
702                dat = G2el.getFFvalues(FFtables,SQ)
703                refDict['FF']['FF'][iref] *= dat.values()
704#reflection processing begins here - big arrays!
705    iBeg = 0
706    time0 = time.time()
707    while iBeg < nRef:
708        iFin = min(iBeg+blkSize,nRef)
709        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
710        H = refl.T[:3]                          #array(blkSize,3)
711        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
712        TwMask = np.any(H,axis=-1)
713        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
714            for ir in range(blkSize):
715                iref = ir+iBeg
716                if iref in TwDict:
717                    for i in TwDict[iref]:
718                        for n in range(NTL):
719                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
720            TwMask = np.any(H,axis=-1)
721        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
722        SQfactor = 4.0*SQ*twopisq               #ditto prev.
723        if 'T' in calcControls[hfx+'histType']:
724            if 'P' in calcControls[hfx+'histType']:
725                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
726            else:
727                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
728            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
729            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
730        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
731        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
732        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
733        Uniq = np.inner(H,SGMT)
734        Phi = np.inner(H,SGT)
735        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
736        sinp = np.sin(phase)
737        cosp = np.cos(phase)
738        biso = -SQfactor*Uisodata[:,nxs]
739        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
740        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
741        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
742        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
743        if 'T' in calcControls[hfx+'histType']:
744            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
745            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
746        else:
747            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
748            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
749        fas = np.sum(np.sum(fa,axis=-1),axis=-1)       #real sum over atoms & unique hkl
750        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag sum over atoms & uniq hkl
751        if SGData['SGInv']: #centrosymmetric; B=0
752            fbs[0] *= 0.
753        if 'P' in calcControls[hfx+'histType']:
754            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
755            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
756        else:
757            if len(TwinLaw) > 1:
758                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
759                refl.T[7] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
760                    np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
761                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
762            else:
763                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
764                refl.T[7] = np.copy(refl.T[9])               
765                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
766#        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
767        iBeg += blkSize
768#    print ' %d sf time %.4f\r'%(nRef,time.time()-time0)
769   
770def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
771    '''Compute structure factor derivatives on blocks of reflections - keep as it works for twins
772    but is slower for powders/nontwins
773    '''
774    phfx = pfx.split(':')[0]+hfx
775    ast = np.sqrt(np.diag(G))
776    Mast = twopisq*np.multiply.outer(ast,ast)
777    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
778    SGT = np.array([ops[1] for ops in SGData['SGOps']])
779    FFtables = calcControls['FFtables']
780    BLtables = calcControls['BLtables']
781    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
782    TwDict = refDict.get('TwDict',{})           
783    if 'S' in calcControls[hfx+'histType']:
784        NTL = calcControls[phfx+'NTL']
785        NM = calcControls[phfx+'TwinNMN']+1
786        TwinLaw = calcControls[phfx+'TwinLaw']
787        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
788        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
789    nTwin = len(TwinLaw)       
790    nRef = len(refDict['RefList'])
791    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
792    mSize = len(Mdata)
793    FF = np.zeros(len(Tdata))
794    if 'NC' in calcControls[hfx+'histType']:
795        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
796    elif 'X' in calcControls[hfx+'histType']:
797        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
798        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
799    Uij = np.array(G2lat.U6toUij(Uijdata))
800    bij = Mast*Uij.T
801    dFdvDict = {}
802    if nTwin > 1:
803        dFdfr = np.zeros((nRef,nTwin,mSize))
804        dFdx = np.zeros((nRef,nTwin,mSize,3))
805        dFdui = np.zeros((nRef,nTwin,mSize))
806        dFdua = np.zeros((nRef,nTwin,mSize,6))
807        dFdbab = np.zeros((nRef,nTwin,2))
808        dFdfl = np.zeros((nRef,nTwin))
809        dFdtw = np.zeros((nRef,nTwin))
810    else:
811        dFdfr = np.zeros((nRef,mSize))
812        dFdx = np.zeros((nRef,mSize,3))
813        dFdui = np.zeros((nRef,mSize))
814        dFdua = np.zeros((nRef,mSize,6))
815        dFdbab = np.zeros((nRef,2))
816        dFdfl = np.zeros((nRef))
817        dFdtw = np.zeros((nRef))
818    Flack = 1.0
819    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
820        Flack = 1.-2.*parmDict[phfx+'Flack']
821    time0 = time.time()
822    nref = len(refDict['RefList'])/100   
823    for iref,refl in enumerate(refDict['RefList']):
824        if 'T' in calcControls[hfx+'histType']:
825            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
826        H = np.array(refl[:3])
827        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
828        TwMask = np.any(H,axis=-1)
829        if TwinLaw.shape[0] > 1 and TwDict:
830            if iref in TwDict:
831                for i in TwDict[iref]:
832                    for n in range(NTL):
833                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
834            TwMask = np.any(H,axis=-1)
835        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
836        SQfactor = 8.0*SQ*np.pi**2
837        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
838        Bab = parmDict[phfx+'BabA']*dBabdA
839        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
840        FF = refDict['FF']['FF'][iref].T[Tindx].T
841        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
842        Phi = np.inner(H,SGT)
843        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
844        sinp = np.sin(phase)
845        cosp = np.cos(phase)
846        occ = Mdata*Fdata/len(SGT)
847        biso = -SQfactor*Uisodata[:,nxs]
848        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
849        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
850        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
851        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
852        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
853        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
854        fot = (FF+FP-Bab)*Tcorr
855        fotp = FPP*Tcorr       
856        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
857        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
858#        GSASIIpath.IPyBreak()
859        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
860        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
861        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
862        fbx = np.array([fot*cosp,-fotp*sinp])
863        #sum below is over Uniq
864        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
865        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
866        dfadui = np.sum(-SQfactor*fa,axis=-2)
867        if nTwin > 1:
868            dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
869            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
870            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
871        else:
872            dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)
873            dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,nxs],axis=-2)
874            # array(2,nAtom,3) & array(2,nAtom,6)
875        if not SGData['SGInv']:
876            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
877            dfadba /= 2.
878            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
879            dfbdui = np.sum(-SQfactor*fb,axis=-2)
880            if len(TwinLaw) > 1:
881                dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
882                dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
883            else:
884                dfadfl = np.sum(-FPP*Tcorr*sinp)
885                dfbdfl = np.sum(FPP*Tcorr*cosp)
886                dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=2)           
887                dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,nxs],axis=2)
888        else:
889            dfbdfr = np.zeros_like(dfadfr)
890            dfbdx = np.zeros_like(dfadx)
891            dfbdui = np.zeros_like(dfadui)
892            dfbdua = np.zeros_like(dfadua)
893            dfbdba = np.zeros_like(dfadba)
894            dfadfl = 0.0
895            dfbdfl = 0.0
896        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
897        SA = fas[0]+fas[1]
898        SB = fbs[0]+fbs[1]
899        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
900            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
901                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
902            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
903                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
904            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
905                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
906            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
907                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
908        else:
909            if nTwin > 1:
910                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)]
911                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)]
912                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)]
913                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)]
914                dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
915               
916            else:   #these are good for no twin single crystals
917                dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)
918                dFdx[iref] = 2.*SA*(dfadx[0]+dfadx[1])+2.*SB*(dfbdx[0]+dfbdx[1])
919                dFdui[iref] = 2.*SA*(dfadui[0]+dfadui[1])+2.*SB*(dfbdui[0]+dfbdui[1])
920                dFdua[iref] = 2.*SA*(dfadua[0]+dfadua[1])+2.*SB*(dfbdua[0]+dfbdua[1])
921                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
922        dFdbab[iref] = fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
923            fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
924#        GSASIIpath.IPyBreak()
925        if not iref%100 :
926            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
927    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
928        #loop over atoms - each dict entry is list of derivatives for all the reflections
929    if nTwin > 1:
930        for i in range(len(Mdata)):     #these all OK?
931            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
932            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
933            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
934            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
935            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
936            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
937            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
938            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
939            dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
940            dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
941            dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
942    else:
943        for i in range(len(Mdata)):
944            dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
945            dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
946            dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
947            dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
948            dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
949            dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
950            dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
951            dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
952            dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
953            dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
954            dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
955        dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
956    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
957    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
958    if nTwin > 1:
959        for i in range(nTwin):
960            dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
961    return dFdvDict
962   
963def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
964    '''Compute structure factor derivatives on blocks of reflections- works for powders/nontwins
965    faster than StructureFactorDerv
966    '''
967    phfx = pfx.split(':')[0]+hfx
968    ast = np.sqrt(np.diag(G))
969    Mast = twopisq*np.multiply.outer(ast,ast)
970    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
971    SGT = np.array([ops[1] for ops in SGData['SGOps']])
972    FFtables = calcControls['FFtables']
973    BLtables = calcControls['BLtables']
974    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
975    TwDict = refDict.get('TwDict',{})           
976    if 'S' in calcControls[hfx+'histType']:
977        NTL = calcControls[phfx+'NTL']
978        NM = calcControls[phfx+'TwinNMN']+1
979        TwinLaw = calcControls[phfx+'TwinLaw']
980        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
981        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
982    nTwin = len(TwinLaw)       
983    nRef = len(refDict['RefList'])
984    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
985    mSize = len(Mdata)
986    FF = np.zeros(len(Tdata))
987    if 'NC' in calcControls[hfx+'histType']:
988        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
989    elif 'X' in calcControls[hfx+'histType']:
990        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
991        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
992    Uij = np.array(G2lat.U6toUij(Uijdata))
993    bij = Mast*Uij.T
994    dFdvDict = {}
995    if nTwin > 1:
996        dFdfr = np.zeros((nRef,nTwin,mSize))
997        dFdx = np.zeros((nRef,nTwin,mSize,3))
998        dFdui = np.zeros((nRef,nTwin,mSize))
999        dFdua = np.zeros((nRef,nTwin,mSize,6))
1000        dFdbab = np.zeros((nRef,nTwin,2))
1001        dFdfl = np.zeros((nRef,nTwin))
1002        dFdtw = np.zeros((nRef,nTwin))
1003    else:
1004        dFdfr = np.zeros((nRef,mSize))
1005        dFdx = np.zeros((nRef,mSize,3))
1006        dFdui = np.zeros((nRef,mSize))
1007        dFdua = np.zeros((nRef,mSize,6))
1008        dFdbab = np.zeros((nRef,2))
1009        dFdfl = np.zeros((nRef))
1010        dFdtw = np.zeros((nRef))
1011    Flack = 1.0
1012    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1013        Flack = 1.-2.*parmDict[phfx+'Flack']
1014    time0 = time.time()
1015    nref = len(refDict['RefList'])/100   
1016#reflection processing begins here - big arrays!
1017    iBeg = 0
1018    blkSize = 32       #no. of reflections in a block - optimized for speed
1019    while iBeg < nRef:
1020        iFin = min(iBeg+blkSize,nRef)
1021        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1022        H = refl.T[:3]
1023        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
1024        TwMask = np.any(H,axis=-1)
1025        if TwinLaw.shape[0] > 1 and TwDict:
1026            for ir in range(blkSize):
1027                iref = ir+iBeg
1028                if iref in TwDict:
1029                    for i in TwDict[iref]:
1030                        for n in range(NTL):
1031                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1032            TwMask = np.any(H,axis=-1)
1033        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1034        SQfactor = 8.0*SQ*np.pi**2
1035        if 'T' in calcControls[hfx+'histType']:
1036            if 'P' in calcControls[hfx+'histType']:
1037                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1038            else:
1039                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1040            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
1041            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
1042        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1043        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
1044        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1045        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
1046        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
1047        Phi = np.inner(H,SGT)
1048        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1049        sinp = np.sin(phase)        #refBlk x nOps x nAtoms
1050        cosp = np.cos(phase)
1051        occ = Mdata*Fdata/len(SGT)
1052        biso = -SQfactor*Uisodata[:,nxs]
1053        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
1054        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1055        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1056        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
1057        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1058        if nTwin > 1:
1059            Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6))
1060        else:
1061            Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6))
1062        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
1063        fotp = FPP*Tcorr       
1064        if 'T' in calcControls[hfx+'histType']:
1065            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1066            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
1067        else:
1068            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1069            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
1070#        GSASIIpath.IPyBreak()
1071        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
1072        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1073        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,ntwi,nEqv,nAtoms)
1074        fbx = np.array([fot*cosp,-fotp*sinp])
1075        #sum below is over Uniq
1076        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,ntwin,nAtom) Fdata != 0 avoids /0. problem
1077        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
1078        if nTwin > 1:
1079            dfadfr = np.swapaxes(dfadfr,0,1)
1080            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)
1081            dfadui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms)
1082            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)
1083            # array(nTwin,2,refBlk,nAtom,3) & array(nTwin,2,refBlk,nAtom,6)
1084        else:
1085            dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
1086            dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nTw,nAtoms)
1087            dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
1088            # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
1089        if not SGData['SGInv']:
1090            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
1091            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
1092            if len(TwinLaw) > 1:
1093                dfbdfr = np.swapaxes(dfbdfr,0,1)
1094                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)         
1095                dfbdui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),0,1)      #array(Ops,refBlk,nTw,nAtoms)
1096                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)
1097            else:
1098                dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1)
1099                dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1)
1100                dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
1101                dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
1102                dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
1103        else:
1104            dfbdfr = np.zeros_like(dfadfr)
1105            dfbdx = np.zeros_like(dfadx)
1106            dfbdui = np.zeros_like(dfadui)
1107            dfbdua = np.zeros_like(dfadua)
1108            dfbdba = np.zeros_like(dfadba)
1109            dfadfl = 0.0
1110            dfbdfl = 0.0
1111        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1112        SA = fas[0]+fas[1]
1113        SB = fbs[0]+fbs[1]
1114        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
1115            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+   \
1116                2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT)
1117            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
1118                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
1119            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
1120                2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
1121            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
1122                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
1123        else:
1124            if nTwin > 1:
1125                dFdfr[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*(np.sum(dfadfr[:,:,it],axis=1))+ \
1126                    SB[:,it,nxs]*(np.sum(dfbdfr[:,:,it],axis=1)))*Mdata/len(SGMT) for it in range(nTwin)]),0,1)
1127                dFdx[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadx[:,:,it],axis=1)+    \
1128                    SB[:,it,nxs,nxs]*np.sum(dfbdx[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
1129                dFdui[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*np.sum(dfadui[:,:,it],axis=1)+  \
1130                    SB[:,it,nxs]*np.sum(dfbdui[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
1131                dFdua[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadua[:,:,it],axis=1)+  \
1132                    SB[:,it,nxs,nxs]*np.sum(dfbdua[:,:,it],axis=1)) for it in range(nTwin)]),0,1)
1133                dFdtw[iBeg:iFin] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
1134               
1135            else:   #these are good for no twin single crystals
1136                dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
1137                dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
1138                dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
1139                dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
1140                dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
1141        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)])+ \
1142                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
1143#        GSASIIpath.IPyBreak()
1144        iBeg += blkSize
1145    print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
1146        #loop over atoms - each dict entry is list of derivatives for all the reflections
1147    if nTwin > 1:
1148        for i in range(len(Mdata)):     #these all OK?
1149            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1150            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1151            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1152            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1153            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1154            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1155            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1156            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1157            dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1158            dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1159            dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1160    else:
1161        for i in range(len(Mdata)):
1162            dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1163            dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1164            dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1165            dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1166            dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1167            dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1168            dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1169            dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1170            dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1171            dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1172            dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1173        dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
1174    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1175    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1176    if nTwin > 1:
1177        for i in range(nTwin):
1178            dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1179    return dFdvDict
1180   
1181def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1182    '''
1183    Compute super structure factors for all h,k,l,m for phase
1184    puts the result, F^2, in each ref[8+im] in refList
1185    input:
1186   
1187    :param dict refDict: where
1188        'RefList' list where each ref = h,k,l,m,d,...
1189        'FF' dict of form factors - filed in below
1190    :param np.array G:      reciprocal metric tensor
1191    :param str pfx:    phase id string
1192    :param dict SGData: space group info. dictionary output from SpcGroup
1193    :param dict calcControls:
1194    :param dict ParmDict:
1195
1196    '''
1197    phfx = pfx.split(':')[0]+hfx
1198    ast = np.sqrt(np.diag(G))
1199    Mast = twopisq*np.multiply.outer(ast,ast)   
1200    SGInv = SGData['SGInv']
1201    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1202    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1203    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1204    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1205    FFtables = calcControls['FFtables']
1206    BLtables = calcControls['BLtables']
1207    Flack = 1.0
1208    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1209        Flack = 1.-2.*parmDict[phfx+'Flack']
1210    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1211    TwDict = refDict.get('TwDict',{})           
1212    if 'S' in calcControls[hfx+'histType']:
1213        NTL = calcControls[phfx+'NTL']
1214        NM = calcControls[phfx+'TwinNMN']+1
1215        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1216        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1217        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1218    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1219    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1220    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1221    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1222    FF = np.zeros(len(Tdata))
1223    if 'NC' in calcControls[hfx+'histType']:
1224        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1225    elif 'X' in calcControls[hfx+'histType']:
1226        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1227        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1228    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1229    bij = Mast*Uij
1230    blkSize = 100       #no. of reflections in a block
1231    nRef = refDict['RefList'].shape[0]
1232    if not len(refDict['FF']):
1233        if 'N' in calcControls[hfx+'histType']:
1234            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1235            refDict['FF']['El'] = dat.keys()
1236            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
1237        else:
1238            dat = G2el.getFFvalues(FFtables,0.)       
1239            refDict['FF']['El'] = dat.keys()
1240            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
1241            for iref,ref in enumerate(refDict['RefList']):
1242                SQ = 1./(2.*ref[5])**2
1243                dat = G2el.getFFvalues(FFtables,SQ)
1244                refDict['FF']['FF'][iref] *= dat.values()
1245    time0 = time.time()
1246#reflection processing begins here - big arrays!
1247    iBeg = 0
1248    while iBeg < nRef:
1249        iFin = min(iBeg+blkSize,nRef)
1250        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1251        H = refl.T[:4]                          #array(blkSize,4)
1252        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1253#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,4) or (blkSize,4)
1254#        TwMask = np.any(H,axis=-1)
1255#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1256#            for ir in range(blkSize):
1257#                iref = ir+iBeg
1258#                if iref in TwDict:
1259#                    for i in TwDict[iref]:
1260#                        for n in range(NTL):
1261#                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1262#            TwMask = np.any(H,axis=-1)
1263        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1264        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1265        Uniq = np.inner(H.T,SSGMT)
1266        UniqP = np.inner(HP.T,SGMT)
1267        Phi = np.inner(H.T,SSGT)
1268        if SGInv:   #if centro - expand HKL sets
1269            Uniq = np.hstack((Uniq,-Uniq))
1270            Phi = np.hstack((Phi,-Phi))
1271            UniqP = np.hstack((UniqP,-UniqP))
1272        if 'T' in calcControls[hfx+'histType']:
1273            if 'P' in calcControls[hfx+'histType']:
1274                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1275            else:
1276                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1277            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1278            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1279        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1280        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1281        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1282        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1283        sinp = np.sin(phase)
1284        cosp = np.cos(phase)
1285        biso = -SQfactor*Uisodata[:,nxs]
1286        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1287        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1288        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1289        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1290        if 'T' in calcControls[hfx+'histType']:
1291            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1292            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1293        else:
1294            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1295            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1296        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1297        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1298        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1299        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1300        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1301#        GSASIIpath.IPyBreak()
1302        if 'P' in calcControls[hfx+'histType']:
1303            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1304#            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
1305            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1306        else:
1307            if len(TwinLaw) > 1:
1308                refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1309                refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1310                    np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1311                refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1312            else:
1313                refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1314                refl.T[8] = np.copy(refl.T[10])               
1315                refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1316        iBeg += blkSize
1317    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
1318
1319def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1320    'Needs a doc string'
1321    phfx = pfx.split(':')[0]+hfx
1322    ast = np.sqrt(np.diag(G))
1323    Mast = twopisq*np.multiply.outer(ast,ast)
1324    SGInv = SGData['SGInv']
1325    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1326    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1327    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1328    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1329    FFtables = calcControls['FFtables']
1330    BLtables = calcControls['BLtables']
1331    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
1332    TwDict = refDict.get('TwDict',{})           
1333    if 'S' in calcControls[hfx+'histType']:
1334        NTL = calcControls[phfx+'NTL']
1335        NM = calcControls[phfx+'TwinNMN']+1
1336        TwinLaw = calcControls[phfx+'TwinLaw']
1337        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1338        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1339    nTwin = len(TwinLaw)       
1340    nRef = len(refDict['RefList'])
1341    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1342    mSize = len(Mdata)  #no. atoms
1343    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1344    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1345    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1346    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1347    FF = np.zeros(len(Tdata))
1348    if 'NC' in calcControls[hfx+'histType']:
1349        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1350    elif 'X' in calcControls[hfx+'histType']:
1351        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1352        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1353    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1354    bij = Mast*Uij
1355    if not len(refDict['FF']):
1356        if 'N' in calcControls[hfx+'histType']:
1357            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1358        else:
1359            dat = G2el.getFFvalues(FFtables,0.)       
1360        refDict['FF']['El'] = dat.keys()
1361        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1362    dFdvDict = {}
1363    if nTwin > 1:
1364        dFdfr = np.zeros((nRef,nTwin,mSize))
1365        dFdx = np.zeros((nRef,nTwin,mSize,3))
1366        dFdui = np.zeros((nRef,nTwin,mSize))
1367        dFdua = np.zeros((nRef,nTwin,mSize,6))
1368        dFdbab = np.zeros((nRef,nTwin,2))
1369        dFdfl = np.zeros((nRef,nTwin))
1370        dFdtw = np.zeros((nRef,nTwin))
1371        dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
1372        dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
1373        dFdGz = np.zeros((nRef,nTwin,mSize,5))
1374        dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
1375    else:
1376        dFdfr = np.zeros((nRef,mSize))
1377        dFdx = np.zeros((nRef,mSize,3))
1378        dFdui = np.zeros((nRef,mSize))
1379        dFdua = np.zeros((nRef,mSize,6))
1380        dFdbab = np.zeros((nRef,2))
1381        dFdfl = np.zeros((nRef))
1382        dFdtw = np.zeros((nRef))
1383        dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1384        dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1385        dFdGz = np.zeros((nRef,mSize,5))
1386        dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1387    Flack = 1.0
1388    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1389        Flack = 1.-2.*parmDict[phfx+'Flack']
1390    time0 = time.time()
1391    nRef = len(refDict['RefList'])/100
1392    for iref,refl in enumerate(refDict['RefList']):
1393        if 'T' in calcControls[hfx+'histType']:
1394            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1395        H = np.array(refl[:4])
1396        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1397#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(4,nTwins) or (4)
1398#        TwMask = np.any(H,axis=-1)
1399#        if TwinLaw.shape[0] > 1 and TwDict:
1400#            if iref in TwDict:
1401#                for i in TwDict[iref]:
1402#                    for n in range(NTL):
1403#                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1404#            TwMask = np.any(H,axis=-1)
1405        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1406        SQfactor = 8.0*SQ*np.pi**2
1407        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1408        Bab = parmDict[phfx+'BabA']*dBabdA
1409        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1410        FF = refDict['FF']['FF'][iref].T[Tindx]
1411        Uniq = np.inner(H,SSGMT)
1412        Phi = np.inner(H,SSGT)
1413        UniqP = np.inner(HP,SGMT)
1414        if SGInv:   #if centro - expand HKL sets
1415            Uniq = np.vstack((Uniq,-Uniq))
1416            Phi = np.hstack((Phi,-Phi))
1417            UniqP = np.vstack((UniqP,-UniqP))
1418        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1419        sinp = np.sin(phase)
1420        cosp = np.cos(phase)
1421        occ = Mdata*Fdata/Uniq.shape[0]
1422        biso = -SQfactor*Uisodata[:,nxs]
1423        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
1424        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1425        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1426        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
1427        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1428        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1429        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1430        fotp = FPP*Tcorr            #ops x atoms
1431        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1432        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1433        # GfpuA is 2 x ops x atoms
1434        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1435        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
1436        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1437        fag = fa*GfpuA[0]-fb*GfpuA[1]
1438        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1439       
1440        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1441        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1442        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
1443        fbx = np.array([fot*cosp,-fotp*sinp])
1444        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1445        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1446        #sum below is over Uniq
1447        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1448        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1449        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1450        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1451        dfadui = np.sum(-SQfactor*fag,axis=1)
1452        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1453        if nTwin > 1:
1454            dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1455            dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
1456            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1457            dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1458            # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
1459            dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
1460            dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
1461            dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
1462            dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
1463            dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
1464            dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
1465            dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
1466            dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
1467        else:
1468            dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1469            dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1470            dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1471            dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1472            # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1473            dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1474            dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1475            dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1476            dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1477            dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
1478            dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
1479            dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1480            dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1481#        GSASIIpath.IPyBreak()
1482        if not SGData['SGInv'] and len(TwinLaw) == 1:   #Flack derivative
1483            dfadfl = np.sum(-FPP*Tcorr*sinp)
1484            dfbdfl = np.sum(FPP*Tcorr*cosp)
1485        else:
1486            dfadfl = 1.0
1487            dfbdfl = 1.0
1488        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1489        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
1490        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
1491        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
1492            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1493            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1494                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1495            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1496                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1497            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1498                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1499            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1500                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1501            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1502                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1503            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1504                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1505            dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+  \
1506                2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1])
1507            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1508                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1509        else:
1510            if nTwin > 1:
1511                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)]
1512                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)]
1513                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)]
1514                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)]
1515                dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
1516
1517                dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
1518                dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
1519                dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
1520                dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
1521            else:   #these are good for no twin single crystals
1522                dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1523                dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1524                dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1525                dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1526                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1527                               
1528                dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1529                dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1530                dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
1531                dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1532#            GSASIIpath.IPyBreak()
1533        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1534            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1535        #loop over atoms - each dict entry is list of derivatives for all the reflections
1536        if not iref%100 :
1537            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
1538    for i in range(len(Mdata)):     #loop over atoms
1539        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1540        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1541        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1542        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1543        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1544        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1545        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1546        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1547        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1548        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1549        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1550        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
1551            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1552            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1553        nx = 0
1554        if waveTypes[i] in ['Block','ZigZag']:
1555            nx = 1 
1556            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
1557            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
1558            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
1559            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
1560            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
1561        for j in range(XSSdata.shape[1]-nx):       #loop over waves
1562            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
1563            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
1564            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
1565            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
1566            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
1567            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
1568        for j in range(USSdata.shape[1]):       #loop over waves
1569            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1570            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1571            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1572            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
1573            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
1574            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
1575            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1576            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1577            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1578            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
1579            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
1580            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
1581           
1582#        GSASIIpath.IPyBreak()
1583    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1584    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1585    return dFdvDict
1586   
1587def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1588    ''' Single crystal extinction function; returns extinction & derivative
1589    '''
1590    extCor = 1.0
1591    dervDict = {}
1592    dervCor = 1.0
1593    if calcControls[phfx+'EType'] != 'None':
1594        SQ = 1/(4.*ref[4+im]**2)
1595        if 'C' in parmDict[hfx+'Type']:           
1596            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
1597        else:   #'T'
1598            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
1599        if 'SXC' in parmDict[hfx+'Type']:
1600            AV = 7.9406e5/parmDict[pfx+'Vol']**2
1601            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1602            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
1603            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
1604        elif 'SNT' in parmDict[hfx+'Type']:
1605            AV = 1.e7/parmDict[pfx+'Vol']**2
1606            PL = SQ
1607            PLZ = AV*ref[9+im]*ref[12+im]**2
1608        elif 'SNC' in parmDict[hfx+'Type']:
1609            AV = 1.e7/parmDict[pfx+'Vol']**2
1610            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1611            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
1612           
1613        if 'Primary' in calcControls[phfx+'EType']:
1614            PLZ *= 1.5
1615        else:
1616            if 'C' in parmDict[hfx+'Type']:
1617                PLZ *= calcControls[phfx+'Tbar']
1618            else: #'T'
1619                PLZ *= ref[13+im]      #t-bar
1620        if 'Primary' in calcControls[phfx+'EType']:
1621            PLZ *= 1.5
1622            PSIG = parmDict[phfx+'Ep']
1623        elif 'I & II' in calcControls[phfx+'EType']:
1624            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
1625        elif 'Type II' in calcControls[phfx+'EType']:
1626            PSIG = parmDict[phfx+'Es']
1627        else:       # 'Secondary Type I'
1628            PSIG = parmDict[phfx+'Eg']/PL
1629           
1630        AG = 0.58+0.48*cos2T+0.24*cos2T**2
1631        AL = 0.025+0.285*cos2T
1632        BG = 0.02-0.025*cos2T
1633        BL = 0.15-0.2*(0.75-cos2T)**2
1634        if cos2T < 0.:
1635            BL = -0.45*cos2T
1636        CG = 2.
1637        CL = 2.
1638        PF = PLZ*PSIG
1639       
1640        if 'Gaussian' in calcControls[phfx+'EApprox']:
1641            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
1642            extCor = np.sqrt(PF4)
1643            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
1644        else:
1645            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
1646            extCor = np.sqrt(PF4)
1647            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
1648
1649        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
1650        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1651            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
1652        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1653            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
1654        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1655            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
1656               
1657    return 1./extCor,dervDict,dervCor
1658   
1659def Dict2Values(parmdict, varylist):
1660    '''Use before call to leastsq to setup list of values for the parameters
1661    in parmdict, as selected by key in varylist'''
1662    return [parmdict[key] for key in varylist] 
1663   
1664def Values2Dict(parmdict, varylist, values):
1665    ''' Use after call to leastsq to update the parameter dictionary with
1666    values corresponding to keys in varylist'''
1667    parmdict.update(zip(varylist,values))
1668   
1669def GetNewCellParms(parmDict,varyList):
1670    'Needs a doc string'
1671    newCellDict = {}
1672    Anames = ['A'+str(i) for i in range(6)]
1673    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1674    for item in varyList:
1675        keys = item.split(':')
1676        if keys[2] in Ddict:
1677            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
1678            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
1679            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1680    return newCellDict          # is e.g. {'0::D11':A0-D11}
1681   
1682def ApplyXYZshifts(parmDict,varyList):
1683    '''
1684    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1685   
1686    :param dict parmDict: parameter dictionary
1687    :param list varyList: list of variables (not used!)
1688    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1689
1690    '''
1691    newAtomDict = {}
1692    for item in parmDict:
1693        if 'dA' in item:
1694            parm = ''.join(item.split('d'))
1695            parmDict[parm] += parmDict[item]
1696            newAtomDict[item] = [parm,parmDict[parm]]
1697    return newAtomDict
1698   
1699def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1700    'Spherical harmonics texture'
1701    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1702    if 'T' in calcControls[hfx+'histType']:
1703        tth = parmDict[hfx+'2-theta']
1704    else:
1705        tth = refl[5+im]
1706    odfCor = 1.0
1707    H = refl[:3]
1708    cell = G2lat.Gmat2cell(g)
1709    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1710    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1711    phi,beta = G2lat.CrsAng(H,cell,SGData)
1712    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1713    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1714    for item in SHnames:
1715        L,M,N = eval(item.strip('C'))
1716        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1717        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1718        Lnorm = G2lat.Lnorm(L)
1719        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1720    return odfCor
1721   
1722def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1723    'Spherical harmonics texture derivatives'
1724    if 'T' in calcControls[hfx+'histType']:
1725        tth = parmDict[hfx+'2-theta']
1726    else:
1727        tth = refl[5+im]
1728    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1729    odfCor = 1.0
1730    dFdODF = {}
1731    dFdSA = [0,0,0]
1732    H = refl[:3]
1733    cell = G2lat.Gmat2cell(g)
1734    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1735    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1736    phi,beta = G2lat.CrsAng(H,cell,SGData)
1737    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1738    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1739    for item in SHnames:
1740        L,M,N = eval(item.strip('C'))
1741        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1742        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1743        Lnorm = G2lat.Lnorm(L)
1744        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1745        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1746        for i in range(3):
1747            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1748    return odfCor,dFdODF,dFdSA
1749   
1750def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1751    'spherical harmonics preferred orientation (cylindrical symmetry only)'
1752    if 'T' in calcControls[hfx+'histType']:
1753        tth = parmDict[hfx+'2-theta']
1754    else:
1755        tth = refl[5+im]
1756    odfCor = 1.0
1757    H = refl[:3]
1758    cell = G2lat.Gmat2cell(g)
1759    Sangls = [0.,0.,0.]
1760    if 'Bragg' in calcControls[hfx+'instType']:
1761        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1762        IFCoup = True
1763    else:
1764        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1765        IFCoup = False
1766    phi,beta = G2lat.CrsAng(H,cell,SGData)
1767    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1768    SHnames = calcControls[phfx+'SHnames']
1769    for item in SHnames:
1770        L,N = eval(item.strip('C'))
1771        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1772        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1773        Lnorm = G2lat.Lnorm(L)
1774        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1775    return np.squeeze(odfCor)
1776   
1777def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1778    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1779    if 'T' in calcControls[hfx+'histType']:
1780        tth = parmDict[hfx+'2-theta']
1781    else:
1782        tth = refl[5+im]
1783    odfCor = 1.0
1784    dFdODF = {}
1785    H = refl[:3]
1786    cell = G2lat.Gmat2cell(g)
1787    Sangls = [0.,0.,0.]
1788    if 'Bragg' in calcControls[hfx+'instType']:
1789        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1790        IFCoup = True
1791    else:
1792        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1793        IFCoup = False
1794    phi,beta = G2lat.CrsAng(H,cell,SGData)
1795    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1796    SHnames = calcControls[phfx+'SHnames']
1797    for item in SHnames:
1798        L,N = eval(item.strip('C'))
1799        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1800        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1801        Lnorm = G2lat.Lnorm(L)
1802        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1803        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
1804    return odfCor,dFdODF
1805   
1806def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1807    'March-Dollase preferred orientation correction'
1808    POcorr = 1.0
1809    MD = parmDict[phfx+'MD']
1810    if MD != 1.0:
1811        MDAxis = calcControls[phfx+'MDAxis']
1812        sumMD = 0
1813        for H in uniq:           
1814            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1815            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1816            sumMD += A**3
1817        POcorr = sumMD/len(uniq)
1818    return POcorr
1819   
1820def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1821    'Needs a doc string'
1822    POcorr = 1.0
1823    POderv = {}
1824    if calcControls[phfx+'poType'] == 'MD':
1825        MD = parmDict[phfx+'MD']
1826        MDAxis = calcControls[phfx+'MDAxis']
1827        sumMD = 0
1828        sumdMD = 0
1829        for H in uniq:           
1830            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1831            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1832            sumMD += A**3
1833            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1834        POcorr = sumMD/len(uniq)
1835        POderv[phfx+'MD'] = sumdMD/len(uniq)
1836    else:   #spherical harmonics
1837        if calcControls[phfx+'SHord']:
1838            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1839    return POcorr,POderv
1840   
1841def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1842    'Needs a doc string'
1843    if 'Debye' in calcControls[hfx+'instType']:
1844        if 'T' in calcControls[hfx+'histType']:
1845            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1846        else:
1847            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1848    else:
1849        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1850   
1851def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1852    'Needs a doc string'
1853    if 'Debye' in calcControls[hfx+'instType']:
1854        if 'T' in calcControls[hfx+'histType']:
1855            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1856        else:
1857            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1858    else:
1859        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1860       
1861def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1862    'Needs a doc string'
1863    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1864    pi2 = np.sqrt(2./np.pi)
1865    if 'T' in calcControls[hfx+'histType']:
1866        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1867        wave = refl[14+im]
1868    else:   #'C'W
1869        sth2 = sind(refl[5+im]/2.)**2
1870        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1871    c2th = 1.-2.0*sth2
1872    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1873    if 'X' in calcControls[hfx+'histType']:
1874        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1875    xfac = flv2*parmDict[phfx+'Extinction']
1876    exb = 1.0
1877    if xfac > -1.:
1878        exb = 1./np.sqrt(1.+xfac)
1879    exl = 1.0
1880    if 0 < xfac <= 1.:
1881        xn = np.array([xfac**(i+1) for i in range(6)])
1882        exl += np.sum(xn*coef)
1883    elif xfac > 1.:
1884        xfac2 = 1./np.sqrt(xfac)
1885        exl = pi2*(1.-0.125/xfac)*xfac2
1886    return exb*sth2+exl*(1.-sth2)
1887   
1888def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1889    'Needs a doc string'
1890    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1891    pi2 = np.sqrt(2./np.pi)
1892    if 'T' in calcControls[hfx+'histType']:
1893        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1894        wave = refl[14+im]
1895    else:   #'C'W
1896        sth2 = sind(refl[5+im]/2.)**2
1897        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1898    c2th = 1.-2.0*sth2
1899    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1900    if 'X' in calcControls[hfx+'histType']:
1901        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1902    xfac = flv2*parmDict[phfx+'Extinction']
1903    dbde = -500.*flv2
1904    if xfac > -1.:
1905        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
1906    dlde = 0.
1907    if 0 < xfac <= 1.:
1908        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
1909        dlde = np.sum(xn*coef)
1910    elif xfac > 1.:
1911        xfac2 = 1./np.sqrt(xfac)
1912        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
1913       
1914    return dbde*sth2+dlde*(1.-sth2)
1915   
1916def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1917    'Needs a doc string'    #need powder extinction!
1918    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
1919    if 'X' in parmDict[hfx+'Type']:
1920        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1921    POcorr = 1.0
1922    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
1923        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1924    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1925        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1926    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1927        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1928    Icorr *= POcorr
1929    AbsCorr = 1.0
1930    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1931    Icorr *= AbsCorr
1932    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1933    Icorr *= ExtCorr
1934    return Icorr,POcorr,AbsCorr,ExtCorr
1935   
1936def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1937    'Needs a doc string'    #need powder extinction derivs!
1938    dIdsh = 1./parmDict[hfx+'Scale']
1939    dIdsp = 1./parmDict[phfx+'Scale']
1940    if 'X' in parmDict[hfx+'Type']:
1941        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1942        dIdPola /= pola
1943    else:       #'N'
1944        dIdPola = 0.0
1945    dFdODF = {}
1946    dFdSA = [0,0,0]
1947    dIdPO = {}
1948    if pfx+'SHorder' in parmDict:
1949        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1950        for iSH in dFdODF:
1951            dFdODF[iSH] /= odfCor
1952        for i in range(3):
1953            dFdSA[i] /= odfCor
1954    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1955        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1956        for iPO in dIdPO:
1957            dIdPO[iPO] /= POcorr
1958    if 'T' in parmDict[hfx+'Type']:
1959        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1960        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
1961    else:
1962        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1963        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
1964    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1965       
1966def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1967    'Needs a doc string'
1968    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1969        costh = cosd(refl[5+im]/2.)
1970        #crystallite size
1971        if calcControls[phfx+'SizeType'] == 'isotropic':
1972            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1973        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1974            H = np.array(refl[:3])
1975            P = np.array(calcControls[phfx+'SizeAxis'])
1976            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1977            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1978            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1979        else:           #ellipsoidal crystallites
1980            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1981            H = np.array(refl[:3])
1982            lenR = G2pwd.ellipseSize(H,Sij,GB)
1983            Sgam = 1.8*wave/(np.pi*costh*lenR)
1984        #microstrain               
1985        if calcControls[phfx+'MustrainType'] == 'isotropic':
1986            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1987        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1988            H = np.array(refl[:3])
1989            P = np.array(calcControls[phfx+'MustrainAxis'])
1990            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1991            Si = parmDict[phfx+'Mustrain;i']
1992            Sa = parmDict[phfx+'Mustrain;a']
1993            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1994        else:       #generalized - P.W. Stephens model
1995            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1996            Sum = 0
1997            for i,strm in enumerate(Strms):
1998                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1999            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2000    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2001        #crystallite size
2002        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2003            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2004        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2005            H = np.array(refl[:3])
2006            P = np.array(calcControls[phfx+'SizeAxis'])
2007            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2008            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2009            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2010        else:           #ellipsoidal crystallites   #OK
2011            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2012            H = np.array(refl[:3])
2013            lenR = G2pwd.ellipseSize(H,Sij,GB)
2014            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2015        #microstrain               
2016        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2017            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2018        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2019            H = np.array(refl[:3])
2020            P = np.array(calcControls[phfx+'MustrainAxis'])
2021            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2022            Si = parmDict[phfx+'Mustrain;i']
2023            Sa = parmDict[phfx+'Mustrain;a']
2024            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2025        else:       #generalized - P.W. Stephens model  OK
2026            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2027            Sum = 0
2028            for i,strm in enumerate(Strms):
2029                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
2030            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2031           
2032    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2033    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2034    sig /= ateln2
2035    return sig,gam
2036       
2037def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2038    'Needs a doc string'
2039    gamDict = {}
2040    sigDict = {}
2041    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
2042        costh = cosd(refl[5+im]/2.)
2043        tanth = tand(refl[5+im]/2.)
2044        #crystallite size derivatives
2045        if calcControls[phfx+'SizeType'] == 'isotropic':
2046            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2047            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
2048            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2049        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2050            H = np.array(refl[:3])
2051            P = np.array(calcControls[phfx+'SizeAxis'])
2052            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2053            Si = parmDict[phfx+'Size;i']
2054            Sa = parmDict[phfx+'Size;a']
2055            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2056            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2057            Sgam = gami*sqtrm
2058            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2059            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2060            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2061            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2062            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2063            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2064        else:           #ellipsoidal crystallites
2065            const = 1.8*wave/(np.pi*costh)
2066            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2067            H = np.array(refl[:3])
2068            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2069            Sgam = const/lenR
2070            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2071                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2072                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2073        gamDict[phfx+'Size;mx'] = Sgam
2074        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2075               
2076        #microstrain derivatives               
2077        if calcControls[phfx+'MustrainType'] == 'isotropic':
2078            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2079            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2080            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2081        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2082            H = np.array(refl[:3])
2083            P = np.array(calcControls[phfx+'MustrainAxis'])
2084            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2085            Si = parmDict[phfx+'Mustrain;i']
2086            Sa = parmDict[phfx+'Mustrain;a']
2087            gami = 0.018*Si*Sa*tanth/np.pi
2088            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2089            Mgam = gami/sqtrm
2090            dsi = -gami*Si*cosP**2/sqtrm**3
2091            dsa = -gami*Sa*sinP**2/sqtrm**3
2092            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2093            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2094            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2095            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2096        else:       #generalized - P.W. Stephens model
2097            const = 0.018*refl[4+im]**2*tanth/np.pi
2098            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2099            Sum = 0
2100            for i,strm in enumerate(Strms):
2101                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
2102                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2103                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2104            Mgam = const*np.sqrt(Sum)
2105            for i in range(len(Strms)):
2106                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
2107                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
2108        gamDict[phfx+'Mustrain;mx'] = Mgam
2109        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2110    else:   #'T'OF - All checked & OK
2111        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2112            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2113            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2114            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2115        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2116            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2117            H = np.array(refl[:3])
2118            P = np.array(calcControls[phfx+'SizeAxis'])
2119            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2120            Si = parmDict[phfx+'Size;i']
2121            Sa = parmDict[phfx+'Size;a']
2122            gami = const/(Si*Sa)
2123            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2124            Sgam = gami*sqtrm
2125            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2126            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2127            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2128            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2129            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2130            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2131        else:           #OK  ellipsoidal crystallites
2132            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2133            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2134            H = np.array(refl[:3])
2135            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2136            Sgam = const/lenR
2137            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2138                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2139                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2140        gamDict[phfx+'Size;mx'] = Sgam  #OK
2141        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2142               
2143        #microstrain derivatives               
2144        if calcControls[phfx+'MustrainType'] == 'isotropic':
2145            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2146            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2147            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2148        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2149            H = np.array(refl[:3])
2150            P = np.array(calcControls[phfx+'MustrainAxis'])
2151            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2152            Si = parmDict[phfx+'Mustrain;i']
2153            Sa = parmDict[phfx+'Mustrain;a']
2154            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2155            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2156            Mgam = gami/sqtrm
2157            dsi = -gami*Si*cosP**2/sqtrm**3
2158            dsa = -gami*Sa*sinP**2/sqtrm**3
2159            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2160            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2161            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2162            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2163        else:       #generalized - P.W. Stephens model OK
2164            pwrs = calcControls[phfx+'MuPwrs']
2165            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2166            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2167            Sum = 0
2168            for i,strm in enumerate(Strms):
2169                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
2170                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2171                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2172            Mgam = const*np.sqrt(Sum)
2173            for i in range(len(Strms)):
2174                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
2175                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
2176        gamDict[phfx+'Mustrain;mx'] = Mgam
2177        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2178       
2179    return sigDict,gamDict
2180       
2181def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2182    'Needs a doc string'
2183    if im:
2184        h,k,l,m = refl[:4]
2185        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2186        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2187    else:
2188        h,k,l = refl[:3]
2189        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2190    refl[4+im] = d
2191    if 'C' in calcControls[hfx+'histType']:
2192        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2193        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2194        if 'Bragg' in calcControls[hfx+'instType']:
2195            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2196                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2197        else:               #Debye-Scherrer - simple but maybe not right
2198            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2199    elif 'T' in calcControls[hfx+'histType']:
2200        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2201        #do I need sample position effects - maybe?
2202    return pos
2203
2204def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2205    'Needs a doc string'
2206    dpr = 180./np.pi
2207    if im:
2208        h,k,l,m = refl[:4]
2209        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2210        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2211        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
2212    else:
2213        m = 0
2214        h,k,l = refl[:3]       
2215        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2216    dst = np.sqrt(dstsq)
2217    dsp = 1./dst
2218    if 'C' in calcControls[hfx+'histType']:
2219        pos = refl[5+im]-parmDict[hfx+'Zero']
2220        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2221        dpdw = const*dst
2222        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2223        dpdZ = 1.0
2224        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2225            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2226        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2227        if 'Bragg' in calcControls[hfx+'instType']:
2228            dpdSh = -4.*shft*cosd(pos/2.0)
2229            dpdTr = -shft*sind(pos)*100.0
2230            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2231        else:               #Debye-Scherrer - simple but maybe not right
2232            dpdXd = -shft*cosd(pos)
2233            dpdYd = -shft*sind(pos)
2234            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2235    elif 'T' in calcControls[hfx+'histType']:
2236        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2237        dpdZ = 1.0
2238        dpdDC = dsp
2239        dpdDA = dsp**2
2240        dpdDB = 1./dsp
2241        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2242            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2243        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2244           
2245def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2246    'Needs a doc string'
2247    laue = SGData['SGLaue']
2248    uniq = SGData['SGUniq']
2249    h,k,l = refl[:3]
2250    if laue in ['m3','m3m']:
2251        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2252            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2253    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2254        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2255    elif laue in ['3R','3mR']:
2256        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2257    elif laue in ['4/m','4/mmm']:
2258        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2259    elif laue in ['mmm']:
2260        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2261    elif laue in ['2/m']:
2262        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2263        if uniq == 'a':
2264            Dij += parmDict[phfx+'D23']*k*l
2265        elif uniq == 'b':
2266            Dij += parmDict[phfx+'D13']*h*l
2267        elif uniq == 'c':
2268            Dij += parmDict[phfx+'D12']*h*k
2269    else:
2270        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2271            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2272    if 'C' in calcControls[hfx+'histType']:
2273        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2274    else:
2275        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2276           
2277def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2278    'Needs a doc string'
2279    laue = SGData['SGLaue']
2280    uniq = SGData['SGUniq']
2281    h,k,l = refl[:3]
2282    if laue in ['m3','m3m']:
2283        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2284            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2285    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2286        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2287    elif laue in ['3R','3mR']:
2288        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2289    elif laue in ['4/m','4/mmm']:
2290        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2291    elif laue in ['mmm']:
2292        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2293    elif laue in ['2/m']:
2294        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2295        if uniq == 'a':
2296            dDijDict[phfx+'D23'] = k*l
2297        elif uniq == 'b':
2298            dDijDict[phfx+'D13'] = h*l
2299        elif uniq == 'c':
2300            dDijDict[phfx+'D12'] = h*k
2301    else:
2302        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2303            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2304    if 'C' in calcControls[hfx+'histType']:
2305        for item in dDijDict:
2306            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2307    else:
2308        for item in dDijDict:
2309            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
2310    return dDijDict
2311   
2312def GetDij(phfx,SGData,parmDict):
2313    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2314    return G2spc.HStrainVals(HSvals,SGData)
2315               
2316def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2317    'Needs a doc string'
2318    histoList = Histograms.keys()
2319    histoList.sort()
2320    for histogram in histoList:
2321        if 'PWDR' in histogram[:4]:
2322            Histogram = Histograms[histogram]
2323            hId = Histogram['hId']
2324            hfx = ':%d:'%(hId)
2325            Limits = calcControls[hfx+'Limits']
2326            if 'C' in calcControls[hfx+'histType']:
2327                shl = max(parmDict[hfx+'SH/L'],0.0005)
2328                Ka2 = False
2329                kRatio = 0.0
2330                if hfx+'Lam1' in parmDict.keys():
2331                    Ka2 = True
2332                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2333                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
2334            x,y,w,yc,yb,yd = Histogram['Data']
2335            xB = np.searchsorted(x,Limits[0])
2336            xF = np.searchsorted(x,Limits[1])
2337            ymb = np.array(y-yb)
2338            ymb = np.where(ymb,ymb,1.0)
2339            ycmb = np.array(yc-yb)
2340            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
2341            refLists = Histogram['Reflection Lists']
2342            for phase in refLists:
2343                if phase not in Phases:     #skips deleted or renamed phases silently!
2344                    continue
2345                Phase = Phases[phase]
2346                im = 0
2347                if Phase['General']['Type'] in ['modulated','magnetic']:
2348                    im = 1
2349                pId = Phase['pId']
2350                phfx = '%d:%d:'%(pId,hId)
2351                refDict = refLists[phase]
2352                sumFo = 0.0
2353                sumdF = 0.0
2354                sumFosq = 0.0
2355                sumdFsq = 0.0
2356                sumInt = 0.0
2357                for refl in refDict['RefList']:
2358                    if 'C' in calcControls[hfx+'histType']:
2359                        yp = np.zeros_like(yb)
2360                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2361                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2362                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2363                        iFin2 = iFin
2364                        if not iBeg+iFin:       #peak below low limit - skip peak
2365                            continue
2366                        elif not iBeg-iFin:     #peak above high limit - done
2367                            break
2368                        elif iBeg < iFin:
2369                            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
2370                            sumInt += refl[11+im]*refl[9+im]
2371                            if Ka2:
2372                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2373                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2374                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
2375                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
2376                                if iFin2 > iBeg2: 
2377                                    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
2378                                    sumInt += refl[11+im]*refl[9+im]*kRatio
2379                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
2380                               
2381                    elif 'T' in calcControls[hfx+'histType']:
2382                        yp = np.zeros_like(yb)
2383                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2384                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2385                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2386                        if iBeg < iFin:
2387                            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
2388                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
2389                            sumInt += refl[11+im]*refl[9+im]
2390                    Fo = np.sqrt(np.abs(refl[8+im]))
2391                    Fc = np.sqrt(np.abs(refl[9]+im))
2392                    sumFo += Fo
2393                    sumFosq += refl[8+im]**2
2394                    sumdF += np.abs(Fo-Fc)
2395                    sumdFsq += (refl[8+im]-refl[9+im])**2
2396                if sumFo:
2397                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2398                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2399                else:
2400                    Histogram['Residuals'][phfx+'Rf'] = 100.
2401                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
2402                Histogram['Residuals'][phfx+'sumInt'] = sumInt
2403                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
2404                Histogram['Residuals']['hId'] = hId
2405        elif 'HKLF' in histogram[:4]:
2406            Histogram = Histograms[histogram]
2407            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
2408               
2409def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2410    'Needs a doc string'
2411   
2412    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
2413        U = parmDict[hfx+'U']
2414        V = parmDict[hfx+'V']
2415        W = parmDict[hfx+'W']
2416        X = parmDict[hfx+'X']
2417        Y = parmDict[hfx+'Y']
2418        tanPos = tand(refl[5+im]/2.0)
2419        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2420        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
2421        sig = max(0.001,sig)
2422        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
2423        gam = max(0.001,gam)
2424        return sig,gam
2425               
2426    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
2427        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
2428            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
2429        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
2430        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2431        sig += Ssig
2432        gam += Sgam
2433        return sig,gam
2434       
2435    def GetReflAlpBet(refl,im,hfx,parmDict):
2436        alp = parmDict[hfx+'alpha']/refl[4+im]
2437        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
2438        return alp,bet
2439       
2440    hId = Histogram['hId']
2441    hfx = ':%d:'%(hId)
2442    bakType = calcControls[hfx+'bakType']
2443    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2444    yc = np.zeros_like(yb)
2445    cw = np.diff(x)
2446    cw = np.append(cw,cw[-1])
2447       
2448    if 'C' in calcControls[hfx+'histType']:   
2449        shl = max(parmDict[hfx+'SH/L'],0.002)
2450        Ka2 = False
2451        if hfx+'Lam1' in parmDict.keys():
2452            wave = parmDict[hfx+'Lam1']
2453            Ka2 = True
2454            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2455            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2456        else:
2457            wave = parmDict[hfx+'Lam']
2458    for phase in Histogram['Reflection Lists']:
2459        refDict = Histogram['Reflection Lists'][phase]
2460        if phase not in Phases:     #skips deleted or renamed phases silently!
2461            continue
2462        Phase = Phases[phase]
2463        pId = Phase['pId']
2464        pfx = '%d::'%(pId)
2465        phfx = '%d:%d:'%(pId,hId)
2466        hfx = ':%d:'%(hId)
2467        SGData = Phase['General']['SGData']
2468        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2469        im = 0
2470        if Phase['General']['Type'] in ['modulated','magnetic']:
2471            SSGData = Phase['General']['SSGData']
2472            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2473            im = 1  #offset in SS reflection list
2474            #??
2475        Dij = GetDij(phfx,SGData,parmDict)
2476        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2477        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2478        if np.any(np.diag(G)<0.):
2479            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
2480        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2481        Vst = np.sqrt(nl.det(G))    #V*
2482        if not Phase['General'].get('doPawley'):
2483            time0 = time.time()
2484            if im:
2485                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2486            else:
2487                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2488#            print 'sf calc time: %.3fs'%(time.time()-time0)
2489        time0 = time.time()
2490        badPeak = False
2491        for iref,refl in enumerate(refDict['RefList']):
2492            if 'C' in calcControls[hfx+'histType']:
2493                if im:
2494                    h,k,l,m = refl[:4]
2495                else:
2496                    h,k,l = refl[:3]
2497                Uniq = np.inner(refl[:3],SGMT)
2498                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2499                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
2500                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2501                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2502                refl[11+im] *= Vst*Lorenz
2503                 
2504                if Phase['General'].get('doPawley'):
2505                    try:
2506                        if im:
2507                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2508                        else:
2509                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2510                        refl[9+im] = parmDict[pInd]
2511                    except KeyError:
2512#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2513                        continue
2514                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2515                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2516                iFin = np.searchsorted(x,refl[5+im]+fmax)
2517                if not iBeg+iFin:       #peak below low limit - skip peak
2518                    continue
2519                elif not iBeg-iFin:     #peak above high limit - done
2520                    break
2521                elif iBeg > iFin:   #bad peak coeff - skip
2522                    badPeak = True
2523                    continue
2524                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
2525                if Ka2:
2526                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2527                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2528                    iBeg = np.searchsorted(x,pos2-fmin)
2529                    iFin = np.searchsorted(x,pos2+fmax)
2530                    if not iBeg+iFin:       #peak below low limit - skip peak
2531                        continue
2532                    elif not iBeg-iFin:     #peak above high limit - done
2533                        return yc,yb
2534                    elif iBeg > iFin:   #bad peak coeff - skip
2535                        continue
2536                    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
2537            elif 'T' in calcControls[hfx+'histType']:
2538                h,k,l = refl[:3]
2539                Uniq = np.inner(refl[:3],SGMT)
2540                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2541                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
2542#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
2543                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2544                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2545                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2546                refl[11+im] *= Vst*Lorenz
2547                if Phase['General'].get('doPawley'):
2548                    try:
2549                        if im:
2550                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2551                        else:
2552                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2553                        refl[9+im] = parmDict[pInd]
2554                    except KeyError:
2555#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2556                        continue
2557                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2558                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2559                iFin = np.searchsorted(x,refl[5+im]+fmax)
2560                if not iBeg+iFin:       #peak below low limit - skip peak
2561                    continue
2562                elif not iBeg-iFin:     #peak above high limit - done
2563                    break
2564                elif iBeg > iFin:   #bad peak coeff - skip
2565                    badPeak = True
2566                    continue
2567                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]
2568#        print 'profile calc time: %.3fs'%(time.time()-time0)
2569    if badPeak:
2570        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
2571    return yc,yb
2572   
2573def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2574    'Needs a doc string'
2575   
2576    def cellVaryDerv(pfx,SGData,dpdA): 
2577        if SGData['SGLaue'] in ['-1',]:
2578            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2579                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2580        elif SGData['SGLaue'] in ['2/m',]:
2581            if SGData['SGUniq'] == 'a':
2582                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2583            elif SGData['SGUniq'] == 'b':
2584                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2585            else:
2586                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2587        elif SGData['SGLaue'] in ['mmm',]:
2588            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2589        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2590            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2591        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2592            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2593        elif SGData['SGLaue'] in ['3R', '3mR']:
2594            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2595        elif SGData['SGLaue'] in ['m3m','m3']:
2596            return [[pfx+'A0',dpdA[0]]]
2597           
2598    # create a list of dependent variables and set up a dictionary to hold their derivatives
2599    dependentVars = G2mv.GetDependentVars()
2600    depDerivDict = {}
2601    for j in dependentVars:
2602        depDerivDict[j] = np.zeros(shape=(len(x)))
2603    #print 'dependent vars',dependentVars
2604    lenX = len(x)               
2605    hId = Histogram['hId']
2606    hfx = ':%d:'%(hId)
2607    bakType = calcControls[hfx+'bakType']
2608    dMdv = np.zeros(shape=(len(varylist),len(x)))
2609    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2610    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2611        bBpos =varylist.index(hfx+'Back;0')
2612        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2613    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2614    for name in varylist:
2615        if 'Debye' in name:
2616            id = int(name.split(';')[-1])
2617            parm = name[:int(name.rindex(';'))]
2618            ip = names.index(parm)
2619            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2620    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2621    for name in varylist:
2622        if 'BkPk' in name:
2623            parm,id = name.split(';')
2624            id = int(id)
2625            if parm in names:
2626                ip = names.index(parm)
2627                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2628    cw = np.diff(x)
2629    cw = np.append(cw,cw[-1])
2630    Ka2 = False #also for TOF!
2631    if 'C' in calcControls[hfx+'histType']:   
2632        shl = max(parmDict[hfx+'SH/L'],0.002)
2633        if hfx+'Lam1' in parmDict.keys():
2634            wave = parmDict[hfx+'Lam1']
2635            Ka2 = True
2636            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2637            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2638        else:
2639            wave = parmDict[hfx+'Lam']
2640    for phase in Histogram['Reflection Lists']:
2641        refDict = Histogram['Reflection Lists'][phase]
2642        if phase not in Phases:     #skips deleted or renamed phases silently!
2643            continue
2644        Phase = Phases[phase]
2645        SGData = Phase['General']['SGData']
2646        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2647        im = 0
2648        if Phase['General']['Type'] in ['modulated','magnetic']:
2649            SSGData = Phase['General']['SSGData']
2650            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2651            im = 1  #offset in SS reflection list
2652            #??
2653        pId = Phase['pId']
2654        pfx = '%d::'%(pId)
2655        phfx = '%d:%d:'%(pId,hId)
2656        Dij = GetDij(phfx,SGData,parmDict)
2657        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2658        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2659        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2660        if not Phase['General'].get('doPawley'):
2661            time0 = time.time()
2662            if im:
2663                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2664            else:
2665                dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2666#            print 'sf-derv time %.3fs'%(time.time()-time0)
2667            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2668        time0 = time.time()
2669        for iref,refl in enumerate(refDict['RefList']):
2670            if im:
2671                h,k,l,m = refl[:4]
2672            else:
2673                h,k,l = refl[:3]
2674            Uniq = np.inner(refl[:3],SGMT)
2675            if 'T' in calcControls[hfx+'histType']:
2676                wave = refl[14+im]
2677            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2678            if 'C' in calcControls[hfx+'histType']:        #CW powder
2679                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2680            else: #'T'OF
2681                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2682            iBeg = np.searchsorted(x,refl[5+im]-fmin)
2683            iFin = np.searchsorted(x,refl[5+im]+fmax)
2684            if not iBeg+iFin:       #peak below low limit - skip peak
2685                continue
2686            elif not iBeg-iFin:     #peak above high limit - done
2687                break
2688            pos = refl[5+im]
2689            if 'C' in calcControls[hfx+'histType']:
2690                tanth = tand(pos/2.0)
2691                costh = cosd(pos/2.0)
2692                lenBF = iFin-iBeg
2693                dMdpk = np.zeros(shape=(6,lenBF))
2694                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2695                for i in range(5):
2696                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
2697                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2698                if Ka2:
2699                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2700                    iBeg2 = np.searchsorted(x,pos2-fmin)
2701                    iFin2 = np.searchsorted(x,pos2+fmax)
2702                    if iBeg2-iFin2:
2703                        lenBF2 = iFin2-iBeg2
2704                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2705                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2706                        for i in range(5):
2707                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
2708                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
2709                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2710            else:   #'T'OF
2711                lenBF = iFin-iBeg
2712                if lenBF < 0:   #bad peak coeff
2713                    break
2714                dMdpk = np.zeros(shape=(6,lenBF))
2715                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2716                for i in range(6):
2717                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
2718                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
2719            if Phase['General'].get('doPawley'):
2720                dMdpw = np.zeros(len(x))
2721                try:
2722                    if im:
2723                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2724                    else:
2725                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2726                    idx = varylist.index(pIdx)
2727                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2728                    if Ka2: #not for TOF either
2729                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2730                    dMdv[idx] = dMdpw
2731                except: # ValueError:
2732                    pass
2733            if 'C' in calcControls[hfx+'histType']:
2734                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2735                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2736                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2737                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2738                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2739                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2740                    hfx+'DisplaceY':[dpdY,'pos'],}
2741                if 'Bragg' in calcControls[hfx+'instType']:
2742                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2743                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
2744                else:
2745                    names.update({hfx+'Absorption':[dFdAb,'int'],})
2746            else:   #'T'OF
2747                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2748                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2749                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2750                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2751                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2752                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
2753                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
2754                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2755            for name in names:
2756                item = names[name]
2757                if name in varylist:
2758                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2759                    if Ka2 and iFin2-iBeg2:
2760                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2761                elif name in dependentVars:
2762                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2763                    if Ka2 and iFin2-iBeg2:
2764                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2765            for iPO in dIdPO:
2766                if iPO in varylist:
2767                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2768                    if Ka2 and iFin2-iBeg2:
2769                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2770                elif iPO in dependentVars:
2771                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2772                    if Ka2 and iFin2-iBeg2:
2773                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2774            for i,name in enumerate(['omega','chi','phi']):
2775                aname = pfx+'SH '+name
2776                if aname in varylist:
2777                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2778                    if Ka2 and iFin2-iBeg2:
2779                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2780                elif aname in dependentVars:
2781                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2782                    if Ka2 and iFin2-iBeg2:
2783                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2784            for iSH in dFdODF:
2785                if iSH in varylist:
2786                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2787                    if Ka2 and iFin2-iBeg2:
2788                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2789                elif iSH in dependentVars:
2790                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2791                    if Ka2 and iFin2-iBeg2:
2792                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2793            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2794            for name,dpdA in cellDervNames:
2795                if name in varylist:
2796                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2797                    if Ka2 and iFin2-iBeg2:
2798                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2799                elif name in dependentVars:
2800                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2801                    if Ka2 and iFin2-iBeg2:
2802                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2803            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2804            for name in dDijDict:
2805                if name in varylist:
2806                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2807                    if Ka2 and iFin2-iBeg2:
2808                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2809                elif name in dependentVars:
2810                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2811                    if Ka2 and iFin2-iBeg2:
2812                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2813            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2814                if name in varylist:
2815                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2816                    if Ka2 and iFin2-iBeg2:
2817                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2818                elif name in dependentVars:
2819                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2820                    if Ka2 and iFin2-iBeg2:
2821                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2822            if 'C' in calcControls[hfx+'histType']:
2823                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2824            else:   #'T'OF
2825                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2826            for name in gamDict:
2827                if name in varylist:
2828                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2829                    if Ka2 and iFin2-iBeg2:
2830                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2831                elif name in dependentVars:
2832                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2833                    if Ka2 and iFin2-iBeg2:
2834                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2835            for name in sigDict:
2836                if name in varylist:
2837                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2838                    if Ka2 and iFin2-iBeg2:
2839                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2840                elif name in dependentVars:
2841                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2842                    if Ka2 and iFin2-iBeg2:
2843                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2844            for name in ['BabA','BabU']:
2845                if refl[9+im]:
2846                    if phfx+name in varylist:
2847                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
2848                        if Ka2 and iFin2-iBeg2:
2849                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]
2850                    elif phfx+name in dependentVars:                   
2851                        depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
2852                        if Ka2 and iFin2-iBeg2:
2853                            depDerivDict[phfx+name][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]                 
2854            if not Phase['General'].get('doPawley'):
2855                #do atom derivatives -  for RB,F,X & U so far             
2856                corr = dervDict['int']/refl[9+im]
2857                if Ka2 and iFin2-iBeg2:
2858                    corr2 = dervDict2['int']/refl[9+im]
2859                for name in varylist+dependentVars:
2860                    if '::RBV;' in name:
2861                        pass
2862                    else:
2863                        try:
2864                            aname = name.split(pfx)[1][:2]
2865                            if aname not in ['Af','dA','AU','RB','Xs','Xc','Ys','Yc','Zs','Zc','Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
2866                        except IndexError:
2867                            continue
2868                    if name in varylist:
2869                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2870                        if Ka2 and iFin2-iBeg2:
2871                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2872                    elif name in dependentVars:
2873                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2874                        if Ka2 and iFin2-iBeg2:
2875                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2876    #        print 'profile derv time: %.3fs'%(time.time()-time0)
2877    # now process derivatives in constraints
2878    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2879    return dMdv
2880   
2881def UserRejectHKL(ref,im,userReject):
2882    if ref[5+im]/ref[6+im] < userReject['minF/sig']:
2883        return False
2884    elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
2885        return False
2886    elif ref[11+im] < userReject['MinExt']:
2887        return False
2888    elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
2889        return False
2890    return True
2891   
2892def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2893    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2894    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
2895    are returned to either dervRefine or HessRefine.
2896
2897    :returns:
2898    '''
2899    nobs = Histogram['Residuals']['Nobs']
2900    hId = Histogram['hId']
2901    hfx = ':%d:'%(hId)
2902    pfx = '%d::'%(Phase['pId'])
2903    phfx = '%d:%d:'%(Phase['pId'],hId)
2904    SGData = Phase['General']['SGData']
2905    im = 0
2906    if Phase['General']['Type'] in ['modulated','magnetic']:
2907        SSGData = Phase['General']['SSGData']
2908        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2909        im = 1  #offset in SS reflection list
2910    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2911    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2912    TwinLaw = calcControls[phfx+'TwinLaw']
2913    refDict = Histogram['Data']
2914    if parmDict[phfx+'Scale'] < 0.:
2915        parmDict[phfx+'Scale'] = .001
2916    if im:
2917        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2918    else:
2919        if len(TwinLaw) > 1:
2920            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2921        else:
2922            dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2923    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2924    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2925    dependentVars = G2mv.GetDependentVars()
2926    depDerivDict = {}
2927    for j in dependentVars:
2928        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2929    wdf = np.zeros(len(refDict['RefList']))
2930    if calcControls['F**2']:
2931        for iref,ref in enumerate(refDict['RefList']):
2932            if ref[6+im] > 0:
2933                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2934                w = 1.0/ref[6+im]
2935                if ref[3+im] > 0:
2936                    wdf[iref] = w*(ref[5+im]-ref[7+im])
2937                    for j,var in enumerate(varylist):
2938                        if var in dFdvDict:
2939                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2940                    for var in dependentVars:
2941                        if var in dFdvDict:
2942                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2943                    if phfx+'Scale' in varylist:
2944                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[7+im]*ref[11+im]/parmDict[phfx+'Scale']  #OK
2945                    elif phfx+'Scale' in dependentVars:
2946                        depDerivDict[phfx+'Scale'][iref] = w*ref[7+im]*ref[11+im]/parmDict[phfx+'Scale']   #OK
2947                    for item in ['Ep','Es','Eg']:
2948                        if phfx+item in varylist and phfx+item in dervDict:
2949                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2950                        elif phfx+item in dependentVars and phfx+item in dervDict:
2951                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2952                    for item in ['BabA','BabU']:
2953                        if phfx+item in varylist:
2954                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2955                        elif phfx+item in dependentVars:
2956                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2957    else:   #F refinement
2958        for iref,ref in enumerate(refDict['RefList']):
2959            if ref[5+im] > 0.:
2960                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2961                Fo = np.sqrt(ref[5+im])
2962                Fc = np.sqrt(ref[7+im])
2963                w = 1.0/ref[6+im]
2964                if ref[3+im] > 0:
2965                    wdf[iref] = 2.0*Fc*w*(Fo-Fc)
2966                    for j,var in enumerate(varylist):
2967                        if var in dFdvDict:
2968                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2969                    for var in dependentVars:
2970                        if var in dFdvDict:
2971                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2972                    if phfx+'Scale' in varylist:
2973                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[7+im]*ref[11+im]/parmDict[phfx+'Scale']  #OK
2974                    elif phfx+'Scale' in dependentVars:
2975                        depDerivDict[phfx+'Scale'][iref] = w*ref[7+im]*ref[11+im]/parmDict[phfx+'Scale']   #OK                   
2976                    for item in ['Ep','Es','Eg']:   #OK!
2977                        if phfx+item in varylist and phfx+item in dervDict:
2978                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] 
2979                        elif phfx+item in dependentVars and phfx+item in dervDict:
2980                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2981                    for item in ['BabA','BabU']:
2982                        if phfx+item in varylist:
2983                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[phfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2984                        elif phfx+item in dependentVars:
2985                            depDerivDict[phfx+item][iref] = w*dFdvDict[phfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2986    return dMdvh,depDerivDict,wdf
2987   
2988
2989def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2990    '''Loop over histograms and compute derivatives of the fitting
2991    model (M) with respect to all parameters.  Results are returned in
2992    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2993    n is the number of parameters and m is the number of data
2994    points. This can exceed memory when m gets large. This routine is
2995    used when refinement derivatives are selected as "analtytic
2996    Jacobian" in Controls.
2997
2998    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2999    '''
3000    parmDict.update(zip(varylist,values))
3001    G2mv.Dict2Map(parmDict,varylist)
3002    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
3003    nvar = len(varylist)
3004    dMdv = np.empty(0)
3005    histoList = Histograms.keys()
3006    histoList.sort()
3007    for histogram in histoList:
3008        if 'PWDR' in histogram[:4]:
3009            Histogram = Histograms[histogram]
3010            hId = Histogram['hId']
3011            hfx = ':%d:'%(hId)
3012            wtFactor = calcControls[hfx+'wtFactor']
3013            Limits = calcControls[hfx+'Limits']
3014            x,y,w,yc,yb,yd = Histogram['Data']
3015            xB = np.searchsorted(x,Limits[0])
3016            xF = np.searchsorted(x,Limits[1])
3017            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
3018                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
3019        elif 'HKLF' in histogram[:4]:
3020            Histogram = Histograms[histogram]
3021            phase = Histogram['Reflection Lists']
3022            Phase = Phases[phase]
3023            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
3024            hfx = ':%d:'%(Histogram['hId'])
3025            wtFactor = calcControls[hfx+'wtFactor']
3026            # now process derivatives in constraints
3027            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
3028        else:
3029            continue        #skip non-histogram entries
3030        if len(dMdv):
3031            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
3032        else:
3033            dMdv = np.sqrt(wtFactor)*dMdvh
3034           
3035    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
3036    if np.any(pVals):
3037        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
3038        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
3039       
3040    return dMdv
3041
3042def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
3043    '''Loop over histograms and compute derivatives of the fitting
3044    model (M) with respect to all parameters.  For each histogram, the
3045    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
3046    number of parameters and m is the number of data points *in the
3047    histogram*. The (n by n) Hessian is computed from each Jacobian
3048    and it is returned.  This routine is used when refinement
3049    derivatives are selected as "analtytic Hessian" in Controls.
3050
3051    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
3052    '''
3053    parmDict.update(zip(varylist,values))
3054    G2mv.Dict2Map(parmDict,varylist)
3055    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
3056    #fixup H atom positions here?
3057    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
3058    nvar = len(varylist)
3059    Hess = np.empty(0)
3060    histoList = Histograms.keys()
3061    histoList.sort()
3062    for histogram in histoList:
3063        if 'PWDR' in histogram[:4]:
3064            Histogram = Histograms[histogram]
3065            hId = Histogram['hId']
3066            hfx = ':%d:'%(hId)
3067            wtFactor = calcControls[hfx+'wtFactor']
3068            Limits = calcControls[hfx+'Limits']
3069            x,y,w,yc,yb,yd = Histogram['Data']
3070            W = wtFactor*w
3071            dy = y-yc
3072            xB = np.searchsorted(x,Limits[0])
3073            xF = np.searchsorted(x,Limits[1])
3074            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
3075                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
3076            Wt = ma.sqrt(W[xB:xF])[nxs,:]
3077            Dy = dy[xB:xF][nxs,:]
3078            dMdvh *= Wt
3079            if dlg:
3080                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
3081            if len(Hess):
3082                Hess += np.inner(dMdvh,dMdvh)
3083                dMdvh *= Wt*Dy
3084                Vec += np.sum(dMdvh,axis=1)
3085            else:
3086                Hess = np.inner(dMdvh,dMdvh)
3087                dMdvh *= Wt*Dy
3088                Vec = np.sum(dMdvh,axis=1)
3089        elif 'HKLF' in histogram[:4]:
3090            Histogram = Histograms[histogram]
3091            phase = Histogram['Reflection Lists']
3092            Phase = Phases[phase]
3093            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
3094            hId = Histogram['hId']
3095            hfx = ':%d:'%(Histogram['hId'])
3096            wtFactor = calcControls[hfx+'wtFactor']
3097            # now process derivatives in constraints
3098            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
3099#            print 'matrix build time: %.3f'%(time.time()-time0)
3100
3101            if dlg:
3102                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
3103            if len(Hess):
3104                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
3105                Hess += wtFactor*np.inner(dMdvh,dMdvh)
3106            else:
3107                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
3108                Hess = wtFactor*np.inner(dMdvh,dMdvh)
3109        else:
3110            continue        #skip non-histogram entries
3111    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
3112    if np.any(pVals):
3113        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
3114        Vec += np.sum(dpdv*pWt*pVals,axis=1)
3115        Hess += np.inner(dpdv*pWt,dpdv)
3116    return Vec,Hess
3117
3118def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
3119    '''Computes the point-by-point discrepancies between every data point in every histogram
3120    and the observed value
3121   
3122    :returns: an np array of differences between observed and computed diffraction values.
3123    '''
3124    Values2Dict(parmDict, varylist, values)
3125    G2mv.Dict2Map(parmDict,varylist)
3126    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
3127    M = np.empty(0)
3128    SumwYo = 0
3129    Nobs = 0
3130    Nrej = 0
3131    Next = 0
3132    ApplyRBModels(parmDict,Phases,rigidbodyDict)
3133    #fixup Hatom positions here....
3134    histoList = Histograms.keys()
3135    histoList.sort()
3136    for histogram in histoList:
3137        if 'PWDR' in histogram[:4]:
3138            Histogram = Histograms[histogram]
3139            hId = Histogram['hId']
3140            hfx = ':%d:'%(hId)
3141            wtFactor = calcControls[hfx+'wtFactor']
3142            Limits = calcControls[hfx+'Limits']
3143            x,y,w,yc,yb,yd = Histogram['Data']
3144            yc *= 0.0                           #zero full calcd profiles
3145            yb *= 0.0
3146            yd *= 0.0
3147            xB = np.searchsorted(x,Limits[0])
3148            xF = np.searchsorted(x,Limits[1])
3149            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
3150                varylist,Histogram,Phases,calcControls,pawleyLookup)
3151            yc[xB:xF] += yb[xB:xF]
3152            if not np.any(y):                   #fill dummy data
3153                rv = st.poisson(yc[xB:xF])
3154                y[xB:xF] = rv.rvs()
3155                Z = np.ones_like(yc[xB:xF])
3156                Z[1::2] *= -1
3157                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
3158                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
3159            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
3160            W = wtFactor*w
3161            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
3162            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
3163            Nobs += Histogram['Residuals']['Nobs']
3164            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
3165            SumwYo += Histogram['Residuals']['sumwYo']
3166            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
3167            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
3168            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
3169            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
3170            sumYB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yc[xB:xF])*ma.abs(y[xB:xF]-yb[xB:xF])/y[xB:xF],0.))
3171            sumwYB2 = ma.sum(ma.where(yc[xB:xF]!=