source: trunk/GSASIIstrMath.py @ 2094

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

implement reflection block processing for sf derivatives for non twinned data
block size optimized at 32 for sf derivatives & 100 for sf calcs

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