source: trunk/GSASIIstrMath.py @ 2040

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

make sure XSC is reflection import default Type
convert Biso to Uiso for J2K atoms import
get plots to show hklm for SS 2D structure factor plots on status bar & tool tip
make wave coeff 5 places not 4
now get correct Fcalc2 for SS reflections; proper derivatives for x & Uij
(not for x waves or U waves)
correct error in single crystal Fcalc
2 & it's derivatives

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