source: trunk/GSASIIstrMath.py @ 927

Last change on this file since 927 was 927, checked in by vondreele, 9 years ago

more G2str fixes

File size: 94.2 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIstrMath - structure computation math routines
3########### SVN repository information ###################
4# $Date: 2013-05-17 12:13:15 -0500 (Fri, 17 May 2013) $
5# $Author: vondreele $
6# $Revision: 920 $
7# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIstrMath.py $
8# $Id: GSASIIstrMath.py 920 2013-05-17 17:13:15Z vondreele $
9########### SVN repository information ###################
10import time
11import math
12import copy
13import numpy as np
14import numpy.ma as ma
15import numpy.linalg as nl
16import scipy.optimize as so
17import GSASIIpath
18GSASIIpath.SetVersionNumber("$Revision: 920 $")
19import GSASIIElem as G2el
20import GSASIIlattice as G2lat
21import GSASIIspc as G2spc
22import GSASIIpwd as G2pwd
23import GSASIImapvars as G2mv
24import GSASIImath as G2mth
25
26sind = lambda x: np.sin(x*np.pi/180.)
27cosd = lambda x: np.cos(x*np.pi/180.)
28tand = lambda x: np.tan(x*np.pi/180.)
29asind = lambda x: 180.*np.arcsin(x)/np.pi
30acosd = lambda x: 180.*np.arccos(x)/np.pi
31atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
32   
33ateln2 = 8.0*math.log(2.0)
34
35################################################################################
36##### Rigid Body Models
37################################################################################
38       
39def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
40    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
41    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
42    '''
43    atxIds = ['Ax:','Ay:','Az:']
44    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
45    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
46    if not RBIds['Vector'] and not RBIds['Residue']:
47        return
48    VRBIds = RBIds['Vector']
49    RRBIds = RBIds['Residue']
50    if Update:
51        RBData = rigidbodyDict
52    else:
53        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
54    if RBIds['Vector']:                       # first update the vector magnitudes
55        VRBData = RBData['Vector']
56        for i,rbId in enumerate(VRBIds):
57            for j in range(len(VRBData[rbId]['VectMag'])):
58                name = '::RBV;'+str(j)+':'+str(i)
59                VRBData[rbId]['VectMag'][j] = parmDict[name]
60    for phase in Phases:
61        Phase = Phases[phase]
62        General = Phase['General']
63        cell = General['Cell'][1:7]
64        Amat,Bmat = G2lat.cell2AB(cell)
65        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
66        pfx = str(Phase['pId'])+'::'
67        if Update:
68            RBModels = Phase['RBModels']
69        else:
70            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
71        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
72            jrb = VRBIds.index(RBObj['RBId'])
73            rbsx = str(irb)+':'+str(jrb)
74            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
75                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
76            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
77                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
78            TLS = RBObj['ThermalMotion']
79            if 'T' in TLS[0]:
80                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
81                    TLS[1][i] = parmDict[pfx+pt+rbsx]
82            if 'L' in TLS[0]:
83                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
84                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
85            if 'S' in TLS[0]:
86                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
87                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
88            if 'U' in TLS[0]:
89                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
90            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
91            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
92            for i,x in enumerate(XYZ):
93                atId = RBObj['Ids'][i]
94                for j in [0,1,2]:
95                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
96                if UIJ[i][0] == 'A':
97                    for j in range(6):
98                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
99                elif UIJ[i][0] == 'I':
100                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
101           
102        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
103            jrb = RRBIds.index(RBObj['RBId'])
104            rbsx = str(irb)+':'+str(jrb)
105            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
106                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
107            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
108                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
109            TLS = RBObj['ThermalMotion']
110            if 'T' in TLS[0]:
111                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
112                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
113            if 'L' in TLS[0]:
114                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
115                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
116            if 'S' in TLS[0]:
117                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
118                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
119            if 'U' in TLS[0]:
120                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
121            for itors,tors in enumerate(RBObj['Torsions']):
122                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
123            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
124            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
125            for i,x in enumerate(XYZ):
126                atId = RBObj['Ids'][i]
127                for j in [0,1,2]:
128                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
129                if UIJ[i][0] == 'A':
130                    for j in range(6):
131                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
132                elif UIJ[i][0] == 'I':
133                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
134                   
135def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
136    atxIds = ['dAx:','dAy:','dAz:']
137    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
138    TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
139    LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
140    SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
141    PIds = ['Px:','Py:','Pz:']
142    OIds = ['Oa:','Oi:','Oj:','Ok:']
143    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
144    if not RBIds['Vector'] and not RBIds['Residue']:
145        return
146    VRBIds = RBIds['Vector']
147    RRBIds = RBIds['Residue']
148    RBData = rigidbodyDict
149    for item in parmDict:
150        if 'RB' in item:
151            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
152    General = Phase['General']
153    cell = General['Cell'][1:7]
154    Amat,Bmat = G2lat.cell2AB(cell)
155    rpd = np.pi/180.
156    rpd2 = rpd**2
157    g = nl.inv(np.inner(Bmat,Bmat))
158    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
159        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
160    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
161    pfx = str(Phase['pId'])+'::'
162    RBModels =  Phase['RBModels']
163    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
164        VModel = RBData['Vector'][RBObj['RBId']]
165        Q = RBObj['Orient'][0]
166        Pos = RBObj['Orig'][0]
167        jrb = VRBIds.index(RBObj['RBId'])
168        rbsx = str(irb)+':'+str(jrb)
169        dXdv = []
170        for iv in range(len(VModel['VectMag'])):
171            dCdv = []
172            for vec in VModel['rbVect'][iv]:
173                dCdv.append(G2mth.prodQVQ(Q,vec))
174            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
175        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
176        for ia,atId in enumerate(RBObj['Ids']):
177            atNum = AtLookup[atId]
178            dx = 0.0001
179            for iv in range(len(VModel['VectMag'])):
180                for ix in [0,1,2]:
181                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
182            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
183                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
184            for iv in range(4):
185                Q[iv] -= dx
186                XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
187                Q[iv] += 2.*dx
188                XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
189                Q[iv] -= dx
190                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
191                for ix in [0,1,2]:
192                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
193            X = G2mth.prodQVQ(Q,Cart[ia])
194            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
195            dFdu = G2lat.U6toUij(dFdu.T)
196            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
197            dFdu = G2lat.UijtoU6(dFdu)
198            atNum = AtLookup[atId]
199            if 'T' in RBObj['ThermalMotion'][0]:
200                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
201                    dFdvDict[pfx+name+rbsx] += dFdu[i]
202            if 'L' in RBObj['ThermalMotion'][0]:
203                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
204                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
205                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
206                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
207                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
208                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
209                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
210                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
211                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
212            if 'S' in RBObj['ThermalMotion'][0]:
213                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
214                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
215                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
216                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
217                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
218                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
219                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
220                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
221            if 'U' in RBObj['ThermalMotion'][0]:
222                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
223
224
225    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
226        Q = RBObj['Orient'][0]
227        Pos = RBObj['Orig'][0]
228        jrb = RRBIds.index(RBObj['RBId'])
229        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
230        rbsx = str(irb)+':'+str(jrb)
231        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
232        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
233            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
234            orId,pvId = torData[itors][:2]
235            pivotVec = Cart[orId]-Cart[pvId]
236            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
237            QB = G2mth.AVdeg2Q(0.001,pivotVec)
238            for ir in torData[itors][3]:
239                atNum = AtLookup[RBObj['Ids'][ir]]
240                rVec = Cart[ir]-Cart[pvId]
241                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
242                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
243                for ix in [0,1,2]:
244                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
245        for ia,atId in enumerate(RBObj['Ids']):
246            atNum = AtLookup[atId]
247            dx = 0.0001
248            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
249                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
250            for iv in range(4):
251                Q[iv] -= dx
252                XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
253                Q[iv] += 2.*dx
254                XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
255                Q[iv] -= dx
256                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
257                for ix in [0,1,2]:
258                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
259            X = G2mth.prodQVQ(Q,Cart[ia])
260            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
261            dFdu = G2lat.U6toUij(dFdu.T)
262            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
263            dFdu = G2lat.UijtoU6(dFdu)
264            atNum = AtLookup[atId]
265            if 'T' in RBObj['ThermalMotion'][0]:
266                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
267                    dFdvDict[pfx+name+rbsx] += dFdu[i]
268            if 'L' in RBObj['ThermalMotion'][0]:
269                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
270                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
271                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
272                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
273                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
274                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
275                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
276                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
277                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
278            if 'S' in RBObj['ThermalMotion'][0]:
279                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
280                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
281                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
282                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
283                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
284                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
285                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
286                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
287            if 'U' in RBObj['ThermalMotion'][0]:
288                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
289   
290################################################################################
291##### Penalty & restraint functions
292################################################################################
293
294def penaltyFxn(HistoPhases,parmDict,varyList):
295    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
296    pNames = []
297    pVals = []
298    pWt = []
299    negWt = {}
300    for phase in Phases:
301        pId = Phases[phase]['pId']
302        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
303        General = Phases[phase]['General']
304        textureData = General['SH Texture']
305        SGData = General['SGData']
306        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
307        cell = General['Cell'][1:7]
308        Amat,Bmat = G2lat.cell2AB(cell)
309        if phase not in restraintDict:
310            continue
311        phaseRest = restraintDict[phase]
312        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
313            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
314            ['ChemComp','Sites'],['Texture','HKLs']]
315        for name,rest in names:
316            itemRest = phaseRest[name]
317            if itemRest[rest] and itemRest['Use']:
318                wt = itemRest['wtFactor']
319                if name in ['Bond','Angle','Plane','Chiral']:
320                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
321                        pNames.append(str(pId)+':'+name+':'+str(i))
322                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
323                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
324                        if name == 'Bond':
325                            calc = G2mth.getRestDist(XYZ,Amat)
326                        elif name == 'Angle':
327                            calc = G2mth.getRestAngle(XYZ,Amat)
328                        elif name == 'Plane':
329                            calc = G2mth.getRestPlane(XYZ,Amat)
330                        elif name == 'Chiral':
331                            calc = G2mth.getRestChiral(XYZ,Amat)
332                        pVals.append(obs-calc)
333                        pWt.append(wt/esd**2)
334                elif name in ['Torsion','Rama']:
335                    coeffDict = itemRest['Coeff']
336                    for i,[indx,ops,cofName,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 == 'Torsion':
341                            tor = G2mth.getRestTorsion(XYZ,Amat)
342                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
343                        else:
344                            phi,psi = G2mth.getRestRama(XYZ,Amat)
345                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
346                        pVals.append(obs-calc)
347                        pWt.append(wt/esd**2)
348                elif name == 'ChemComp':
349                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
350                        pNames.append(str(pId)+':'+name+':'+str(i))
351                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
352                        frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
353                        calc = np.sum(mul*frac*factors)
354                        pVals.append(obs-calc)
355                        pWt.append(wt/esd**2)                   
356                elif name == 'Texture':
357                    SHkeys = textureData['SH Coeff'][1].keys()
358                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
359                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
360                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
361                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
362                        PH = np.array(hkl)
363                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
364                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
365                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
366                        Z1 = -ma.masked_greater(Z,0.0)
367                        IndZ1 = np.array(ma.nonzero(Z1))
368                        for ind in IndZ1.T:
369                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
370                            pVals.append(Z1[ind[0]][ind[1]])
371                            pWt.append(wt/esd1**2)
372                        if ifesd2:
373                            Z2 = 1.-Z
374                            for ind in np.ndindex(grid,grid):
375                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
376                                pVals.append(Z1[ind[0]][ind[1]])
377                                pWt.append(wt/esd2**2)
378         
379    for item in varyList:
380        if 'PWLref' in item and parmDict[item] < 0.:
381            pId = int(item.split(':')[0])
382            if negWt[pId]:
383                pNames.append(item)
384                pVals.append(-parmDict[item])
385                pWt.append(negWt[pId])
386    pVals = np.array(pVals)
387    pWt = np.array(pWt)         #should this be np.sqrt?
388    return pNames,pVals,pWt
389   
390def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
391    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
392    pDerv = np.zeros((len(varyList),len(pVal)))
393    for phase in Phases:
394#        if phase not in restraintDict:
395#            continue
396        pId = Phases[phase]['pId']
397        General = Phases[phase]['General']
398        SGData = General['SGData']
399        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
400        cell = General['Cell'][1:7]
401        Amat,Bmat = G2lat.cell2AB(cell)
402        textureData = General['SH Texture']
403
404        SHkeys = textureData['SH Coeff'][1].keys()
405        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
406        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
407        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
408        sam = SamSym[textureData['Model']]
409        phaseRest = restraintDict.get(phase,{})
410        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
411            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
412            'ChemComp':'Sites','Texture':'HKLs'}
413        lasthkl = np.array([0,0,0])
414        for ip,pName in enumerate(pNames):
415            pnames = pName.split(':')
416            if pId == int(pnames[0]):
417                name = pnames[1]
418                if 'PWL' in pName:
419                    pDerv[varyList.index(pName)][ip] += 1.
420                    continue
421                id = int(pnames[2]) 
422                itemRest = phaseRest[name]
423                if name in ['Bond','Angle','Plane','Chiral']:
424                    indx,ops,obs,esd = itemRest[names[name]][id]
425                    dNames = []
426                    for ind in indx:
427                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
428                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
429                    if name == 'Bond':
430                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
431                    elif name == 'Angle':
432                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
433                    elif name == 'Plane':
434                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
435                    elif name == 'Chiral':
436                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
437                elif name in ['Torsion','Rama']:
438                    coffDict = itemRest['Coeff']
439                    indx,ops,cofName,esd = itemRest[names[name]][id]
440                    dNames = []
441                    for ind in indx:
442                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
443                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
444                    if name == 'Torsion':
445                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
446                    else:
447                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
448                elif name == 'ChemComp':
449                    indx,factors,obs,esd = itemRest[names[name]][id]
450                    dNames = []
451                    for ind in indx:
452                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
453                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
454                        deriv = mul*factors
455                elif 'Texture' in name:
456                    deriv = []
457                    dNames = []
458                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
459                    hkl = np.array(hkl)
460                    if np.any(lasthkl-hkl):
461                        PH = np.array(hkl)
462                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
463                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
464                        lasthkl = copy.copy(hkl)                       
465                    if 'unit' in name:
466                        pass
467                    else:
468                        gam = float(pnames[3])
469                        psi = float(pnames[4])
470                        for SHname in ODFln:
471                            l,m,n = eval(SHname[1:])
472                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
473                            dNames += [str(pId)+'::'+SHname]
474                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
475                for dName,drv in zip(dNames,deriv):
476                    try:
477                        ind = varyList.index(dName)
478                        pDerv[ind][ip] += drv
479                    except ValueError:
480                        pass
481    return pDerv
482
483################################################################################
484##### Function & derivative calculations
485################################################################################       
486                   
487def GetAtomFXU(pfx,calcControls,parmDict):
488    Natoms = calcControls['Natoms'][pfx]
489    Tdata = Natoms*[' ',]
490    Mdata = np.zeros(Natoms)
491    IAdata = Natoms*[' ',]
492    Fdata = np.zeros(Natoms)
493    FFdata = []
494    BLdata = []
495    Xdata = np.zeros((3,Natoms))
496    dXdata = np.zeros((3,Natoms))
497    Uisodata = np.zeros(Natoms)
498    Uijdata = np.zeros((6,Natoms))
499    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
500        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
501        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
502        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
503        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
504    for iatm in range(Natoms):
505        for key in keys:
506            parm = pfx+key+str(iatm)
507            if parm in parmDict:
508                keys[key][iatm] = parmDict[parm]
509    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
510   
511def getFFvalues(FFtables,SQ):
512    FFvals = {}
513    for El in FFtables:
514        FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
515    return FFvals
516   
517def getBLvalues(BLtables):
518    BLvals = {}
519    for El in BLtables:
520        BLvals[El] = BLtables[El][1][1]
521    return BLvals
522       
523def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
524    ''' Compute structure factors for all h,k,l for phase
525    input:
526        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
527        G:      reciprocal metric tensor
528        pfx:    phase id string
529        SGData: space group info. dictionary output from SpcGroup
530        calcControls:
531        ParmDict:
532    puts result F^2 in each ref[8] in refList
533    '''       
534    twopi = 2.0*np.pi
535    twopisq = 2.0*np.pi**2
536    phfx = pfx.split(':')[0]+hfx
537    ast = np.sqrt(np.diag(G))
538    Mast = twopisq*np.multiply.outer(ast,ast)
539    FFtables = calcControls['FFtables']
540    BLtables = calcControls['BLtables']
541    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
542    FF = np.zeros(len(Tdata))
543    if 'N' in parmDict[hfx+'Type']:
544        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
545    else:
546        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
547        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
548    maxPos = len(SGData['SGOps'])
549    Uij = np.array(G2lat.U6toUij(Uijdata))
550    bij = Mast*Uij.T
551    for refl in refList:
552        fbs = np.array([0,0])
553        H = refl[:3]
554        SQ = 1./(2.*refl[4])**2
555        SQfactor = 4.0*SQ*twopisq
556        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
557        if not len(refl[-1]):                #no form factors
558            if 'N' in parmDict[hfx+'Type']:
559                refl[-1] = getBLvalues(BLtables)
560            else:       #'X'
561                refl[-1] = getFFvalues(FFtables,SQ)
562        for i,El in enumerate(Tdata):
563            FF[i] = refl[-1][El]           
564        Uniq = refl[11]
565        phi = refl[12]
566        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
567        sinp = np.sin(phase)
568        cosp = np.cos(phase)
569        occ = Mdata*Fdata/len(Uniq)
570        biso = -SQfactor*Uisodata
571        Tiso = np.where(biso<1.,np.exp(biso),1.0)
572        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
573        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
574        Tcorr = Tiso*Tuij
575        fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
576        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
577        if not SGData['SGInv']:
578            fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
579            fbs = np.sum(np.sum(fb,axis=1),axis=1)
580        fasq = fas**2
581        fbsq = fbs**2        #imaginary
582        refl[9] = np.sum(fasq)+np.sum(fbsq)
583        refl[10] = atan2d(fbs[0],fas[0])
584   
585def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
586    twopi = 2.0*np.pi
587    twopisq = 2.0*np.pi**2
588    phfx = pfx.split(':')[0]+hfx
589    ast = np.sqrt(np.diag(G))
590    Mast = twopisq*np.multiply.outer(ast,ast)
591    FFtables = calcControls['FFtables']
592    BLtables = calcControls['BLtables']
593    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
594    FF = np.zeros(len(Tdata))
595    if 'N' in parmDict[hfx+'Type']:
596        FP = 0.
597        FPP = 0.
598    else:
599        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
600        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
601    maxPos = len(SGData['SGOps'])       
602    Uij = np.array(G2lat.U6toUij(Uijdata))
603    bij = Mast*Uij.T
604    dFdvDict = {}
605    dFdfr = np.zeros((len(refList),len(Mdata)))
606    dFdx = np.zeros((len(refList),len(Mdata),3))
607    dFdui = np.zeros((len(refList),len(Mdata)))
608    dFdua = np.zeros((len(refList),len(Mdata),6))
609    dFdbab = np.zeros((len(refList),2))
610    for iref,refl in enumerate(refList):
611        H = np.array(refl[:3])
612        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
613        SQfactor = 8.0*SQ*np.pi**2
614        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
615        Bab = parmDict[phfx+'BabA']*dBabdA
616        for i,El in enumerate(Tdata):           
617            FF[i] = refl[-1][El]           
618        Uniq = refl[11]
619        phi = refl[12]
620        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
621        sinp = np.sin(phase)
622        cosp = np.cos(phase)
623        occ = Mdata*Fdata/len(Uniq)
624        biso = -SQfactor*Uisodata
625        Tiso = np.where(biso<1.,np.exp(biso),1.0)
626        HbH = -np.inner(H,np.inner(bij,H))
627        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
628        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
629        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
630        Tcorr = Tiso*Tuij
631        fot = (FF+FP-Bab)*occ*Tcorr
632        fotp = FPP*occ*Tcorr
633        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
634        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
635       
636        fas = np.sum(np.sum(fa,axis=1),axis=1)
637        fbs = np.sum(np.sum(fb,axis=1),axis=1)
638        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
639        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
640        #sum below is over Uniq
641        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
642        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
643        dfadui = np.sum(-SQfactor*fa,axis=2)
644        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
645        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
646        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
647        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
648        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
649        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
650        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
651        dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
652        if not SGData['SGInv']:
653            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
654            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
655            dfbdui = np.sum(-SQfactor*fb,axis=2)
656            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
657            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
658            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
659            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
660            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
661            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
662            dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
663        #loop over atoms - each dict entry is list of derivatives for all the reflections
664    for i in range(len(Mdata)):     
665        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
666        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
667        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
668        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
669        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
670        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
671        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
672        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
673        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
674        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
675        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
676        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
677        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
678    return dFdvDict
679   
680def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
681    ''' Single crystal extinction function; puts correction in ref[13] and returns
682    corrections needed for derivatives
683    '''
684    ref[13] = 1.0
685    dervCor = 1.0
686    dervDict = {}
687    if calcControls[phfx+'EType'] != 'None':
688        cos2T = 1.0-0.5*(parmDict[hfx+'Lam']/ref[4])**2         #cos(2theta)
689        if 'SXC' in parmDict[hfx+'Type']:
690            AV = 7.9406e5/parmDict[pfx+'Vol']**2
691            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
692            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
693        elif 'SNT' in parmDict[hfx+'Type']:
694            AV = 1.e7/parmDict[pfx+'Vol']**2
695            PL = 1./(4.*refl[4]**2)
696            P12 = 1.0
697        elif 'SNC' in parmDict[hfx+'Type']:
698            AV = 1.e7/parmDict[pfx+'Vol']**2
699            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
700            P12 = 1.0
701           
702        PLZ = AV*P12*parmDict[hfx+'Lam']**2*ref[7]
703        if 'Primary' in calcControls[phfx+'EType']:
704            PLZ *= 1.5
705        else:
706            PLZ *= calcControls[phfx+'Tbar']
707                       
708        if 'Primary' in calcControls[phfx+'EType']:
709            PSIG = parmDict[phfx+'Ep']
710        elif 'I & II' in calcControls[phfx+'EType']:
711            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
712        elif 'Type II' in calcControls[phfx+'EType']:
713            PSIG = parmDict[phfx+'Es']
714        else:       # 'Secondary Type I'
715            PSIG = parmDict[phfx+'Eg']/PL
716           
717        AG = 0.58+0.48*cos2T+0.24*cos2T**2
718        AL = 0.025+0.285*cos2T
719        BG = 0.02-0.025*cos2T
720        BL = 0.15-0.2*(0.75-cos2T)**2
721        if cos2T < 0.:
722            BL = -0.45*cos2T
723        CG = 2.
724        CL = 2.
725        PF = PLZ*PSIG
726       
727        if 'Gaussian' in calcControls[phfx+'EApprox']:
728            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
729            extCor = np.sqrt(PF4)
730            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
731        else:
732            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
733            extCor = np.sqrt(PF4)
734            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
735
736        dervCor = (1.+PF)*PF3
737        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
738            dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
739        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
740            dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
741        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
742            dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
743               
744        ref[13] = 1./extCor
745    return dervCor,dervDict
746       
747   
748def Dict2Values(parmdict, varylist):
749    '''Use before call to leastsq to setup list of values for the parameters
750    in parmdict, as selected by key in varylist'''
751    return [parmdict[key] for key in varylist] 
752   
753def Values2Dict(parmdict, varylist, values):
754    ''' Use after call to leastsq to update the parameter dictionary with
755    values corresponding to keys in varylist'''
756    parmdict.update(zip(varylist,values))
757   
758def GetNewCellParms(parmDict,varyList):
759    newCellDict = {}
760    Anames = ['A'+str(i) for i in range(6)]
761    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
762    for item in varyList:
763        keys = item.split(':')
764        if keys[2] in Ddict:
765            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
766            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
767            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
768    return newCellDict          # is e.g. {'0::D11':A0+D11}
769   
770def ApplyXYZshifts(parmDict,varyList):
771    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
772        input:
773            parmDict - parameter dictionary
774            varyList - list of variables
775        returns:
776            newAtomDict - dictionary of new atomic coordinate names & values;
777                key is parameter shift name
778    '''
779    newAtomDict = {}
780    for item in parmDict:
781        if 'dA' in item:
782            parm = ''.join(item.split('d'))
783            parmDict[parm] += parmDict[item]
784            newAtomDict[item] = [parm,parmDict[parm]]
785    return newAtomDict
786   
787def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
788    #Spherical harmonics texture
789    IFCoup = 'Bragg' in calcControls[hfx+'instType']
790    odfCor = 1.0
791    H = refl[:3]
792    cell = G2lat.Gmat2cell(g)
793    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
794    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
795    phi,beta = G2lat.CrsAng(H,cell,SGData)
796    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
797    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
798    for item in SHnames:
799        L,M,N = eval(item.strip('C'))
800        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
801        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
802        Lnorm = G2lat.Lnorm(L)
803        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
804    return odfCor
805   
806def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
807    #Spherical harmonics texture derivatives
808    FORPI = 4.0*np.pi
809    IFCoup = 'Bragg' in calcControls[hfx+'instType']
810    odfCor = 1.0
811    dFdODF = {}
812    dFdSA = [0,0,0]
813    H = refl[:3]
814    cell = G2lat.Gmat2cell(g)
815    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
816    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
817    phi,beta = G2lat.CrsAng(H,cell,SGData)
818    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
819    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
820    for item in SHnames:
821        L,M,N = eval(item.strip('C'))
822        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
823        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
824        Lnorm = G2lat.Lnorm(L)
825        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
826        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
827        for i in range(3):
828            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
829    return odfCor,dFdODF,dFdSA
830   
831def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
832    #sphericaql harmonics preferred orientation (cylindrical symmetry only)
833    odfCor = 1.0
834    H = refl[:3]
835    cell = G2lat.Gmat2cell(g)
836    Sangl = [0.,0.,0.]
837    if 'Bragg' in calcControls[hfx+'instType']:
838        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
839        IFCoup = True
840    else:
841        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
842        IFCoup = False
843    phi,beta = G2lat.CrsAng(H,cell,SGData)
844    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
845    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
846    for item in SHnames:
847        L,N = eval(item.strip('C'))
848        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
849        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
850    return odfCor
851   
852def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
853    #spherical harmonics preferred orientation derivatives (cylindrical symmetry only)
854    FORPI = 12.5663706143592
855    odfCor = 1.0
856    dFdODF = {}
857    H = refl[:3]
858    cell = G2lat.Gmat2cell(g)
859    Sangl = [0.,0.,0.]
860    if 'Bragg' in calcControls[hfx+'instType']:
861        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
862        IFCoup = True
863    else:
864        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
865        IFCoup = False
866    phi,beta = G2lat.CrsAng(H,cell,SGData)
867    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
868    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
869    for item in SHnames:
870        L,N = eval(item.strip('C'))
871        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
872        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
873        dFdODF[phfx+item] = Kcsl*Lnorm
874    return odfCor,dFdODF
875   
876def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
877    POcorr = 1.0
878    if calcControls[phfx+'poType'] == 'MD':
879        MD = parmDict[phfx+'MD']
880        if MD != 1.0:
881            MDAxis = calcControls[phfx+'MDAxis']
882            sumMD = 0
883            for H in refl[11]:           
884                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
885                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
886                sumMD += A**3
887            POcorr = sumMD/len(refl[11])
888    else:   #spherical harmonics
889        if calcControls[phfx+'SHord']:
890            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
891    return POcorr
892   
893def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
894    POcorr = 1.0
895    POderv = {}
896    if calcControls[phfx+'poType'] == 'MD':
897        MD = parmDict[phfx+'MD']
898        MDAxis = calcControls[phfx+'MDAxis']
899        sumMD = 0
900        sumdMD = 0
901        for H in refl[11]:           
902            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
903            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
904            sumMD += A**3
905            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
906        POcorr = sumMD/len(refl[11])
907        POderv[phfx+'MD'] = sumdMD/len(refl[11])
908    else:   #spherical harmonics
909        if calcControls[phfx+'SHord']:
910            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
911    return POcorr,POderv
912   
913def GetAbsorb(refl,hfx,calcControls,parmDict):
914    if 'Debye' in calcControls[hfx+'instType']:
915        return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
916    else:
917        return 1.0
918   
919def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
920    if 'Debye' in calcControls[hfx+'instType']:
921        return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
922    else:
923        return 0.0
924   
925def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
926    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
927    if 'X' in parmDict[hfx+'Type']:
928        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
929    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
930    if pfx+'SHorder' in parmDict:
931        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
932    Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
933    refl[13] = Icorr       
934   
935def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
936    dIdsh = 1./parmDict[hfx+'Scale']
937    dIdsp = 1./parmDict[phfx+'Scale']
938    if 'X' in parmDict[hfx+'Type']:
939        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
940        dIdPola /= pola
941    else:       #'N'
942        dIdPola = 0.0
943    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
944    for iPO in dIdPO:
945        dIdPO[iPO] /= POcorr
946    dFdODF = {}
947    dFdSA = [0,0,0]
948    if pfx+'SHorder' in parmDict:
949        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
950        for iSH in dFdODF:
951            dFdODF[iSH] /= odfCor
952        for i in range(3):
953            dFdSA[i] /= odfCor
954    dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
955    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
956       
957def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
958    costh = cosd(refl[5]/2.)
959    #crystallite size
960    if calcControls[phfx+'SizeType'] == 'isotropic':
961        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
962    elif calcControls[phfx+'SizeType'] == 'uniaxial':
963        H = np.array(refl[:3])
964        P = np.array(calcControls[phfx+'SizeAxis'])
965        cosP,sinP = G2lat.CosSinAngle(H,P,G)
966        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
967        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
968    else:           #ellipsoidal crystallites
969        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
970        H = np.array(refl[:3])
971        lenR = G2pwd.ellipseSize(H,Sij,GB)
972        Sgam = 1.8*wave/(np.pi*costh*lenR)
973    #microstrain               
974    if calcControls[phfx+'MustrainType'] == 'isotropic':
975        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
976    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
977        H = np.array(refl[:3])
978        P = np.array(calcControls[phfx+'MustrainAxis'])
979        cosP,sinP = G2lat.CosSinAngle(H,P,G)
980        Si = parmDict[phfx+'Mustrain;i']
981        Sa = parmDict[phfx+'Mustrain;a']
982        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
983    else:       #generalized - P.W. Stephens model
984        pwrs = calcControls[phfx+'MuPwrs']
985        sum = 0
986        for i,pwr in enumerate(pwrs):
987            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
988        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
989    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
990    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
991    sig /= ateln2
992    return sig,gam
993       
994def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
995    gamDict = {}
996    sigDict = {}
997    costh = cosd(refl[5]/2.)
998    tanth = tand(refl[5]/2.)
999    #crystallite size derivatives
1000    if calcControls[phfx+'SizeType'] == 'isotropic':
1001        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1002        gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1003        sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1004    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1005        H = np.array(refl[:3])
1006        P = np.array(calcControls[phfx+'SizeAxis'])
1007        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1008        Si = parmDict[phfx+'Size;i']
1009        Sa = parmDict[phfx+'Size;a']
1010        gami = (1.8*wave/np.pi)/(Si*Sa)
1011        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1012        Sgam = gami*sqtrm
1013        gam = Sgam/costh
1014        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
1015        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
1016        gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1017        gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1018        sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1019        sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1020    else:           #ellipsoidal crystallites
1021        const = 1.8*wave/(np.pi*costh)
1022        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1023        H = np.array(refl[:3])
1024        lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1025        Sgam = 1.8*wave/(np.pi*costh*lenR)
1026        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1027            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1028            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1029    gamDict[phfx+'Size;mx'] = Sgam
1030    sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1031           
1032    #microstrain derivatives               
1033    if calcControls[phfx+'MustrainType'] == 'isotropic':
1034        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1035        gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1036        sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1037    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1038        H = np.array(refl[:3])
1039        P = np.array(calcControls[phfx+'MustrainAxis'])
1040        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1041        Si = parmDict[phfx+'Mustrain;i']
1042        Sa = parmDict[phfx+'Mustrain;a']
1043        gami = 0.018*Si*Sa*tanth/np.pi
1044        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1045        Mgam = gami/sqtrm
1046        dsi = -gami*Si*cosP**2/sqtrm**3
1047        dsa = -gami*Sa*sinP**2/sqtrm**3
1048        gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1049        gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1050        sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1051        sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1052    else:       #generalized - P.W. Stephens model
1053        pwrs = calcControls[phfx+'MuPwrs']
1054        const = 0.018*refl[4]**2*tanth
1055        sum = 0
1056        for i,pwr in enumerate(pwrs):
1057            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1058            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1059            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
1060            sigDict[phfx+'Mustrain:'+str(i)] = \
1061                2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1062        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
1063        for i in range(len(pwrs)):
1064            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1065    gamDict[phfx+'Mustrain;mx'] = Mgam
1066    sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1067    return sigDict,gamDict
1068       
1069def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1070    h,k,l = refl[:3]
1071    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1072    d = np.sqrt(dsq)
1073
1074    refl[4] = d
1075    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1076    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1077    if 'Bragg' in calcControls[hfx+'instType']:
1078        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1079            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1080    else:               #Debye-Scherrer - simple but maybe not right
1081        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1082    return pos
1083
1084def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1085    dpr = 180./np.pi
1086    h,k,l = refl[:3]
1087    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1088    dst = np.sqrt(dstsq)
1089    pos = refl[5]-parmDict[hfx+'Zero']
1090    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1091    dpdw = const*dst
1092    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1093    dpdA *= const*wave/(2.0*dst)
1094    dpdZ = 1.0
1095    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1096    if 'Bragg' in calcControls[hfx+'instType']:
1097        dpdSh = -4.*const*cosd(pos/2.0)
1098        dpdTr = -const*sind(pos)*100.0
1099        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1100    else:               #Debye-Scherrer - simple but maybe not right
1101        dpdXd = -const*cosd(pos)
1102        dpdYd = -const*sind(pos)
1103        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1104           
1105def GetHStrainShift(refl,SGData,phfx,parmDict):
1106    laue = SGData['SGLaue']
1107    uniq = SGData['SGUniq']
1108    h,k,l = refl[:3]
1109    if laue in ['m3','m3m']:
1110        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1111            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1112    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1113        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1114    elif laue in ['3R','3mR']:
1115        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1116    elif laue in ['4/m','4/mmm']:
1117        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1118    elif laue in ['mmm']:
1119        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1120    elif laue in ['2/m']:
1121        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1122        if uniq == 'a':
1123            Dij += parmDict[phfx+'D23']*k*l
1124        elif uniq == 'b':
1125            Dij += parmDict[phfx+'D13']*h*l
1126        elif uniq == 'c':
1127            Dij += parmDict[phfx+'D12']*h*k
1128    else:
1129        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1130            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1131    return -Dij*refl[4]**2*tand(refl[5]/2.0)
1132           
1133def GetHStrainShiftDerv(refl,SGData,phfx):
1134    laue = SGData['SGLaue']
1135    uniq = SGData['SGUniq']
1136    h,k,l = refl[:3]
1137    if laue in ['m3','m3m']:
1138        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1139            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1140    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1141        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1142    elif laue in ['3R','3mR']:
1143        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1144    elif laue in ['4/m','4/mmm']:
1145        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1146    elif laue in ['mmm']:
1147        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1148    elif laue in ['2/m']:
1149        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1150        if uniq == 'a':
1151            dDijDict[phfx+'D23'] = k*l
1152        elif uniq == 'b':
1153            dDijDict[phfx+'D13'] = h*l
1154        elif uniq == 'c':
1155            dDijDict[phfx+'D12'] = h*k
1156            names.append()
1157    else:
1158        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1159            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1160    for item in dDijDict:
1161        dDijDict[item] *= -refl[4]**2*tand(refl[5]/2.0)
1162    return dDijDict
1163   
1164def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1165    histoList = Histograms.keys()
1166    histoList.sort()
1167    for histogram in histoList:
1168        if 'PWDR' in histogram[:4]:
1169            Histogram = Histograms[histogram]
1170            hId = Histogram['hId']
1171            hfx = ':%d:'%(hId)
1172            Limits = calcControls[hfx+'Limits']
1173            shl = max(parmDict[hfx+'SH/L'],0.0005)
1174            Ka2 = False
1175            kRatio = 0.0
1176            if hfx+'Lam1' in parmDict.keys():
1177                Ka2 = True
1178                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1179                kRatio = parmDict[hfx+'I(L2)/I(L1)']
1180            x,y,w,yc,yb,yd = Histogram['Data']
1181            xB = np.searchsorted(x,Limits[0])
1182            xF = np.searchsorted(x,Limits[1])
1183            ymb = np.array(y-yb)
1184            ymb = np.where(ymb,ymb,1.0)
1185            ycmb = np.array(yc-yb)
1186            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1187            refLists = Histogram['Reflection Lists']
1188            for phase in refLists:
1189                Phase = Phases[phase]
1190                pId = Phase['pId']
1191                phfx = '%d:%d:'%(pId,hId)
1192                refList = refLists[phase]
1193                sumFo = 0.0
1194                sumdF = 0.0
1195                sumFosq = 0.0
1196                sumdFsq = 0.0
1197                for refl in refList:
1198                    if 'C' in calcControls[hfx+'histType']:
1199                        yp = np.zeros_like(yb)
1200                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1201                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1202                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1203                        iFin2 = iFin
1204                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1205                        if Ka2:
1206                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1207                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1208                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1209                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1210                            if not iBeg2+iFin2:       #peak below low limit - skip peak
1211                                continue
1212                            elif not iBeg2-iFin2:     #peak above high limit - done
1213                                break
1214                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
1215                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
1216                    elif 'T' in calcControls[hfx+'histType']:
1217                        print 'TOF Undefined at present'
1218                        raise Exception    #no TOF yet
1219                    Fo = np.sqrt(np.abs(refl[8]))
1220                    Fc = np.sqrt(np.abs(refl[9]))
1221                    sumFo += Fo
1222                    sumFosq += refl[8]**2
1223                    sumdF += np.abs(Fo-Fc)
1224                    sumdFsq += (refl[8]-refl[9])**2
1225                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1226                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1227                Histogram[phfx+'Nref'] = len(refList)
1228               
1229def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1230   
1231    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1232        U = parmDict[hfx+'U']
1233        V = parmDict[hfx+'V']
1234        W = parmDict[hfx+'W']
1235        X = parmDict[hfx+'X']
1236        Y = parmDict[hfx+'Y']
1237        tanPos = tand(refl[5]/2.0)
1238        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
1239        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1240        sig = max(0.001,sig)
1241        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
1242        gam = max(0.001,gam)
1243        return sig,gam
1244               
1245    hId = Histogram['hId']
1246    hfx = ':%d:'%(hId)
1247    bakType = calcControls[hfx+'bakType']
1248    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1249    yc = np.zeros_like(yb)
1250       
1251    if 'C' in calcControls[hfx+'histType']:   
1252        shl = max(parmDict[hfx+'SH/L'],0.002)
1253        Ka2 = False
1254        if hfx+'Lam1' in parmDict.keys():
1255            wave = parmDict[hfx+'Lam1']
1256            Ka2 = True
1257            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1258            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1259        else:
1260            wave = parmDict[hfx+'Lam']
1261    else:
1262        print 'TOF Undefined at present'
1263        raise ValueError
1264    for phase in Histogram['Reflection Lists']:
1265        refList = Histogram['Reflection Lists'][phase]
1266        Phase = Phases[phase]
1267        pId = Phase['pId']
1268        pfx = '%d::'%(pId)
1269        phfx = '%d:%d:'%(pId,hId)
1270        hfx = ':%d:'%(hId)
1271        SGData = Phase['General']['SGData']
1272        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1273        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1274        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1275        Vst = np.sqrt(nl.det(G))    #V*
1276        if not Phase['General'].get('doPawley'):
1277            time0 = time.time()
1278            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1279#            print 'sf calc time: %.3fs'%(time.time()-time0)
1280        time0 = time.time()
1281        for refl in refList:
1282            if 'C' in calcControls[hfx+'histType']:
1283                h,k,l = refl[:3]
1284                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1285                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1286                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1287                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
1288                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
1289                refl[13] *= Vst*Lorenz
1290                if Phase['General'].get('doPawley'):
1291                    try:
1292                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1293                        refl[9] = parmDict[pInd]
1294                    except KeyError:
1295#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1296                        continue
1297                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1298                iBeg = np.searchsorted(x,refl[5]-fmin)
1299                iFin = np.searchsorted(x,refl[5]+fmax)
1300                if not iBeg+iFin:       #peak below low limit - skip peak
1301                    continue
1302                elif not iBeg-iFin:     #peak above high limit - done
1303                    break
1304                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1305                if Ka2:
1306                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1307                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1308                    iBeg = np.searchsorted(x,pos2-fmin)
1309                    iFin = np.searchsorted(x,pos2+fmax)
1310                    if not iBeg+iFin:       #peak below low limit - skip peak
1311                        continue
1312                    elif not iBeg-iFin:     #peak above high limit - done
1313                        return yc,yb
1314                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
1315            elif 'T' in calcControls[hfx+'histType']:
1316                print 'TOF Undefined at present'
1317                raise Exception    #no TOF yet
1318#        print 'profile calc time: %.3fs'%(time.time()-time0)
1319    return yc,yb
1320   
1321def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1322   
1323    def cellVaryDerv(pfx,SGData,dpdA): 
1324        if SGData['SGLaue'] in ['-1',]:
1325            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1326                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1327        elif SGData['SGLaue'] in ['2/m',]:
1328            if SGData['SGUniq'] == 'a':
1329                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1330            elif SGData['SGUniq'] == 'b':
1331                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1332            else:
1333                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1334        elif SGData['SGLaue'] in ['mmm',]:
1335            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1336        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1337            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1338        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1339            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1340        elif SGData['SGLaue'] in ['3R', '3mR']:
1341            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1342        elif SGData['SGLaue'] in ['m3m','m3']:
1343            return [[pfx+'A0',dpdA[0]]]
1344           
1345    # create a list of dependent variables and set up a dictionary to hold their derivatives
1346    dependentVars = G2mv.GetDependentVars()
1347    depDerivDict = {}
1348    for j in dependentVars:
1349        depDerivDict[j] = np.zeros(shape=(len(x)))
1350    #print 'dependent vars',dependentVars
1351    lenX = len(x)               
1352    hId = Histogram['hId']
1353    hfx = ':%d:'%(hId)
1354    bakType = calcControls[hfx+'bakType']
1355    dMdv = np.zeros(shape=(len(varylist),len(x)))
1356    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1357    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1358        bBpos =varylist.index(hfx+'Back:0')
1359        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1360    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1361    for name in varylist:
1362        if 'Debye' in name:
1363            id = int(name.split(':')[-1])
1364            parm = name[:int(name.rindex(':'))]
1365            ip = names.index(parm)
1366            dMdv[varylist.index(name)] = dMddb[3*id+ip]
1367    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1368    for name in varylist:
1369        if 'BkPk' in name:
1370            id = int(name.split(':')[-1])
1371            parm = name[:int(name.rindex(':'))]
1372            ip = names.index(parm)
1373            dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1374    cw = np.diff(x)
1375    cw = np.append(cw,cw[-1])
1376    if 'C' in calcControls[hfx+'histType']:   
1377        shl = max(parmDict[hfx+'SH/L'],0.002)
1378        Ka2 = False
1379        if hfx+'Lam1' in parmDict.keys():
1380            wave = parmDict[hfx+'Lam1']
1381            Ka2 = True
1382            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1383            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1384        else:
1385            wave = parmDict[hfx+'Lam']
1386    else:
1387        print 'TOF Undefined at present'
1388        raise ValueError
1389    for phase in Histogram['Reflection Lists']:
1390        refList = Histogram['Reflection Lists'][phase]
1391        Phase = Phases[phase]
1392        SGData = Phase['General']['SGData']
1393        pId = Phase['pId']
1394        pfx = '%d::'%(pId)
1395        phfx = '%d:%d:'%(pId,hId)
1396        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1397        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1398        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1399        if not Phase['General'].get('doPawley'):
1400            time0 = time.time()
1401            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1402#            print 'sf-derv time %.3fs'%(time.time()-time0)
1403            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1404        time0 = time.time()
1405        for iref,refl in enumerate(refList):
1406            if 'C' in calcControls[hfx+'histType']:        #CW powder
1407                h,k,l = refl[:3]
1408                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1409                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1410                iBeg = np.searchsorted(x,refl[5]-fmin)
1411                iFin = np.searchsorted(x,refl[5]+fmax)
1412                if not iBeg+iFin:       #peak below low limit - skip peak
1413                    continue
1414                elif not iBeg-iFin:     #peak above high limit - done
1415                    break
1416                pos = refl[5]
1417                tanth = tand(pos/2.0)
1418                costh = cosd(pos/2.0)
1419                lenBF = iFin-iBeg
1420                dMdpk = np.zeros(shape=(6,lenBF))
1421                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
1422                for i in range(5):
1423                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[i]
1424                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1425                if Ka2:
1426                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1427                    iBeg2 = np.searchsorted(x,pos2-fmin)
1428                    iFin2 = np.searchsorted(x,pos2+fmax)
1429                    if iBeg2-iFin2:
1430                        lenBF2 = iFin2-iBeg2
1431                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1432                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
1433                        for i in range(5):
1434                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[i]
1435                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[13]*dMdipk2[0]
1436                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1437                if Phase['General'].get('doPawley'):
1438                    dMdpw = np.zeros(len(x))
1439                    try:
1440                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1441                        idx = varylist.index(pIdx)
1442                        dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1443                        if Ka2:
1444                            dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1445                        dMdv[idx] = dMdpw
1446                    except: # ValueError:
1447                        pass
1448                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1449                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1450                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1451                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1452                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1453                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1454                    hfx+'DisplaceY':[dpdY,'pos'],hfx+'Absorption':[dFdAb,'int'],}
1455                for name in names:
1456                    item = names[name]
1457                    if name in varylist:
1458                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1459                        if Ka2:
1460                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1461                    elif name in dependentVars:
1462                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1463                        if Ka2:
1464                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1465                for iPO in dIdPO:
1466                    if iPO in varylist:
1467                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1468                        if Ka2:
1469                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1470                    elif iPO in dependentVars:
1471                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1472                        if Ka2:
1473                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1474                for i,name in enumerate(['omega','chi','phi']):
1475                    aname = pfx+'SH '+name
1476                    if aname in varylist:
1477                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1478                        if Ka2:
1479                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1480                    elif aname in dependentVars:
1481                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1482                        if Ka2:
1483                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1484                for iSH in dFdODF:
1485                    if iSH in varylist:
1486                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1487                        if Ka2:
1488                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1489                    elif iSH in dependentVars:
1490                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1491                        if Ka2:
1492                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1493                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1494                for name,dpdA in cellDervNames:
1495                    if name in varylist:
1496                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1497                        if Ka2:
1498                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1499                    elif name in dependentVars:
1500                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1501                        if Ka2:
1502                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1503                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1504                for name in dDijDict:
1505                    if name in varylist:
1506                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1507                        if Ka2:
1508                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1509                    elif name in dependentVars:
1510                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1511                        if Ka2:
1512                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1513                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
1514                for name in gamDict:
1515                    if name in varylist:
1516                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1517                        if Ka2:
1518                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1519                    elif name in dependentVars:
1520                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1521                        if Ka2:
1522                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1523                for name in sigDict:
1524                    if name in varylist:
1525                        dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1526                        if Ka2:
1527                            dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1528                    elif name in dependentVars:
1529                        depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1530                        if Ka2:
1531                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1532                for name in ['BabA','BabU']:
1533                    if phfx+name in varylist:
1534                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
1535                        if Ka2:
1536                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]
1537                    elif phfx+name in dependentVars:                   
1538                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
1539                        if Ka2:
1540                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]                 
1541            elif 'T' in calcControls[hfx+'histType']:
1542                print 'TOF Undefined at present'
1543                raise Exception    #no TOF yet
1544            #do atom derivatives -  for RB,F,X & U so far             
1545            corr = dervDict['int']/refl[9]
1546            if Ka2:
1547                corr2 = dervDict2['int']/refl[9]
1548            for name in varylist+dependentVars:
1549                if '::RBV;' in name:
1550                    pass
1551                else:
1552                    try:
1553                        aname = name.split(pfx)[1][:2]
1554                        if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1555                    except IndexError:
1556                        continue
1557                if name in varylist:
1558                    dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1559                    if Ka2:
1560                        dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1561                elif name in dependentVars:
1562                    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1563                    if Ka2:
1564                        depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1565#        print 'profile derv time: %.3fs'%(time.time()-time0)
1566    # now process derivatives in constraints
1567    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1568    return dMdv
1569
1570def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1571    parmDict.update(zip(varylist,values))
1572    G2mv.Dict2Map(parmDict,varylist)
1573    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1574    nvar = len(varylist)
1575    dMdv = np.empty(0)
1576    histoList = Histograms.keys()
1577    histoList.sort()
1578    for histogram in histoList:
1579        if 'PWDR' in histogram[:4]:
1580            Histogram = Histograms[histogram]
1581            hId = Histogram['hId']
1582            hfx = ':%d:'%(hId)
1583            wtFactor = calcControls[hfx+'wtFactor']
1584            Limits = calcControls[hfx+'Limits']
1585            x,y,w,yc,yb,yd = Histogram['Data']
1586            W = wtFactor*w
1587            xB = np.searchsorted(x,Limits[0])
1588            xF = np.searchsorted(x,Limits[1])
1589            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
1590                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1591        elif 'HKLF' in histogram[:4]:
1592            Histogram = Histograms[histogram]
1593            nobs = Histogram['Nobs']
1594            phase = Histogram['Reflection Lists']
1595            Phase = Phases[phase]
1596            hId = Histogram['hId']
1597            hfx = ':%d:'%(hId)
1598            wtFactor = calcControls[hfx+'wtFactor']
1599            pfx = '%d::'%(Phase['pId'])
1600            phfx = '%d:%d:'%(Phase['pId'],hId)
1601            SGData = Phase['General']['SGData']
1602            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1603            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1604            refList = Histogram['Data']
1605            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1606            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1607            dMdvh = np.zeros((len(varylist),len(refList)))
1608            if calcControls['F**2']:
1609                for iref,ref in enumerate(refList):
1610                    if ref[6] > 0:
1611                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1612                        w = 1.0/ref[6]
1613                        if w*ref[5] >= calcControls['minF/sig']:
1614                            for j,var in enumerate(varylist):
1615                                if var in dFdvDict:
1616                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1617                            if phfx+'Scale' in varylist:
1618                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1619                            for item in ['Ep','Es','Eg']:
1620                                if phfx+item in varylist:
1621                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1622                            for item in ['BabA','BabU']:
1623                                if phfx+item in varylist:
1624                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1625            else:
1626                for iref,ref in enumerate(refList):
1627                    if ref[5] > 0.:
1628                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1629                        Fo = np.sqrt(ref[5])
1630                        Fc = np.sqrt(ref[7])
1631                        w = 1.0/ref[6]
1632                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1633                            for j,var in enumerate(varylist):
1634                                if var in dFdvDict:
1635                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1636                            if phfx+'Scale' in varylist:
1637                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor                           
1638                            for item in ['Ep','Es','Eg']:
1639                                if phfx+item in varylist:
1640                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1641                            for item in ['BabA','BabU']:
1642                                if phfx+item in varylist:
1643                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1644        else:
1645            continue        #skip non-histogram entries
1646        if len(dMdv):
1647            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
1648        else:
1649            dMdv = np.sqrt(wtFactor)*dMdvh
1650           
1651    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1652    if np.any(pVals):
1653        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1654        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
1655       
1656    return dMdv
1657
1658def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1659    parmDict.update(zip(varylist,values))
1660    G2mv.Dict2Map(parmDict,varylist)
1661    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1662    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
1663    nvar = len(varylist)
1664    Hess = np.empty(0)
1665    histoList = Histograms.keys()
1666    histoList.sort()
1667    for histogram in histoList:
1668        if 'PWDR' in histogram[:4]:
1669            Histogram = Histograms[histogram]
1670            hId = Histogram['hId']
1671            hfx = ':%d:'%(hId)
1672            wtFactor = calcControls[hfx+'wtFactor']
1673            Limits = calcControls[hfx+'Limits']
1674            x,y,w,yc,yb,yd = Histogram['Data']
1675            W = wtFactor*w
1676            dy = y-yc
1677            xB = np.searchsorted(x,Limits[0])
1678            xF = np.searchsorted(x,Limits[1])
1679            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
1680                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1681            Wt = np.sqrt(W[xB:xF])[np.newaxis,:]
1682            Dy = dy[xB:xF][np.newaxis,:]
1683            dMdvh *= Wt
1684            if dlg:
1685                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1686            if len(Hess):
1687                Hess += np.inner(dMdvh,dMdvh)
1688                dMdvh *= Wt*Dy
1689                Vec += np.sum(dMdvh,axis=1)
1690            else:
1691                Hess = np.inner(dMdvh,dMdvh)
1692                dMdvh *= Wt*Dy
1693                Vec = np.sum(dMdvh,axis=1)
1694        elif 'HKLF' in histogram[:4]:
1695            Histogram = Histograms[histogram]
1696            nobs = Histogram['Nobs']
1697            phase = Histogram['Reflection Lists']
1698            Phase = Phases[phase]
1699            hId = Histogram['hId']
1700            hfx = ':%d:'%(hId)
1701            wtFactor = calcControls[hfx+'wtFactor']
1702            pfx = '%d::'%(Phase['pId'])
1703            phfx = '%d:%d:'%(Phase['pId'],hId)
1704            SGData = Phase['General']['SGData']
1705            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1706            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1707            refList = Histogram['Data']
1708            time0 = time.time()
1709            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)  #accurate for powders!
1710            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1711            dMdvh = np.zeros((len(varylist),len(refList)))
1712            wdf = np.zeros(len(refList))
1713            if calcControls['F**2']:
1714                for iref,ref in enumerate(refList):
1715                    if ref[6] > 0:
1716                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1717                        w =  1.0/ref[6]
1718                        if w*ref[5] >= calcControls['minF/sig']:
1719                            wdf[iref] = w*(ref[5]-ref[7])
1720                            for j,var in enumerate(varylist):
1721                                if var in dFdvDict:
1722                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1723                            if phfx+'Scale' in varylist:
1724                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1725                            for item in ['Ep','Es','Eg']:
1726                                if phfx+item in varylist:
1727                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1728                            for item in ['BabA','BabU']:
1729                                if phfx+item in varylist:
1730                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1731            else:
1732                for iref,ref in enumerate(refList):
1733                    if ref[5] > 0.:
1734                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1735                        Fo = np.sqrt(ref[5])
1736                        Fc = np.sqrt(ref[7])
1737                        w = 1.0/ref[6]
1738                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1739                            wdf[iref] = 2.0*Fo*w*(Fo-Fc)
1740                            for j,var in enumerate(varylist):
1741                                if var in dFdvDict:
1742                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1743                            if phfx+'Scale' in varylist:
1744                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor                           
1745                            for item in ['Ep','Es','Eg']:
1746                                if phfx+item in varylist:
1747                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1748                            for item in ['BabA','BabU']:
1749                                if phfx+item in varylist:
1750                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1751                       
1752            if dlg:
1753                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1754            if len(Hess):
1755                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
1756                Hess += wtFactor*np.inner(dMdvh,dMdvh)
1757            else:
1758                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
1759                Hess = wtFactor*np.inner(dMdvh,dMdvh)
1760        else:
1761            continue        #skip non-histogram entries
1762    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1763    if np.any(pVals):
1764        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1765        Vec += np.sum(dpdv*pWt*pVals,axis=1)
1766        Hess += np.inner(dpdv*pWt,dpdv)
1767    return Vec,Hess
1768
1769def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
1770    parmDict.update(zip(varylist,values))
1771    Values2Dict(parmDict, varylist, values)
1772    G2mv.Dict2Map(parmDict,varylist)
1773    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1774    M = np.empty(0)
1775    SumwYo = 0
1776    Nobs = 0
1777    ApplyRBModels(parmDict,Phases,rigidbodyDict)
1778    histoList = Histograms.keys()
1779    histoList.sort()
1780    for histogram in histoList:
1781        if 'PWDR' in histogram[:4]:
1782            Histogram = Histograms[histogram]
1783            hId = Histogram['hId']
1784            hfx = ':%d:'%(hId)
1785            wtFactor = calcControls[hfx+'wtFactor']
1786            Limits = calcControls[hfx+'Limits']
1787            x,y,w,yc,yb,yd = Histogram['Data']
1788            W = wtFactor*w
1789            yc *= 0.0                           #zero full calcd profiles
1790            yb *= 0.0
1791            yd *= 0.0
1792            xB = np.searchsorted(x,Limits[0])
1793            xF = np.searchsorted(x,Limits[1])
1794            Histogram['Nobs'] = xF-xB
1795            Nobs += Histogram['Nobs']
1796            Histogram['sumwYo'] = np.sum(W[xB:xF]*y[xB:xF]**2)
1797            SumwYo += Histogram['sumwYo']
1798            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
1799                varylist,Histogram,Phases,calcControls,pawleyLookup)
1800            yc[xB:xF] += yb[xB:xF]
1801            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
1802            Histogram['sumwYd'] = np.sum(np.sqrt(W[xB:xF])*(yd[xB:xF]))
1803            wdy = -np.sqrt(W[xB:xF])*(yd[xB:xF])
1804            Histogram['wR'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
1805            if dlg:
1806                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1807            M = np.concatenate((M,wdy))
1808#end of PWDR processing
1809        elif 'HKLF' in histogram[:4]:
1810            Histogram = Histograms[histogram]
1811            phase = Histogram['Reflection Lists']
1812            Phase = Phases[phase]
1813            hId = Histogram['hId']
1814            hfx = ':%d:'%(hId)
1815            wtFactor = calcControls[hfx+'wtFactor']
1816            pfx = '%d::'%(Phase['pId'])
1817            phfx = '%d:%d:'%(Phase['pId'],hId)
1818            SGData = Phase['General']['SGData']
1819            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1820            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1821            refList = Histogram['Data']
1822            time0 = time.time()
1823            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1824#            print 'sf-calc time: %.3f'%(time.time()-time0)
1825            df = np.zeros(len(refList))
1826            sumwYo = 0
1827            sumFo = 0
1828            sumFo2 = 0
1829            sumdF = 0
1830            sumdF2 = 0
1831            nobs = 0
1832            if calcControls['F**2']:
1833                for i,ref in enumerate(refList):
1834                    if ref[6] > 0:
1835                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1836                        w = 1.0/ref[6]
1837                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1838                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1839                        ref[8] = ref[5]/parmDict[phfx+'Scale']
1840                        if w*ref[5] >= calcControls['minF/sig']:
1841                            sumFo2 += ref[5]
1842                            Fo = np.sqrt(ref[5])
1843                            sumFo += Fo
1844                            sumFo2 += ref[5]
1845                            sumdF += abs(Fo-np.sqrt(ref[7]))
1846                            sumdF2 += abs(ref[5]-ref[7])
1847                            nobs += 1
1848                            df[i] = -w*(ref[5]-ref[7])
1849                            sumwYo += (w*ref[5])**2
1850            else:
1851                for i,ref in enumerate(refList):
1852                    if ref[5] > 0.:
1853                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1854                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1855                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1856                        Fo = np.sqrt(ref[5])
1857                        Fc = np.sqrt(ref[7])
1858                        w = 2.0*Fo/ref[6]
1859                        if w*Fo >= calcControls['minF/sig']:
1860                            sumFo += Fo
1861                            sumFo2 += ref[5]
1862                            sumdF += abs(Fo-Fc)
1863                            sumdF2 += abs(ref[5]-ref[7])
1864                            nobs += 1
1865                            df[i] = -w*(Fo-Fc)
1866                            sumwYo += (w*Fo)**2
1867            Histogram['Nobs'] = nobs
1868            Histogram['sumwYo'] = sumwYo
1869            SumwYo += sumwYo
1870            Histogram['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['sumwYo'])*100.)
1871            Histogram[phfx+'Rf'] = 100.*sumdF/sumFo
1872            Histogram[phfx+'Rf^2'] = 100.*sumdF2/sumFo2
1873            Histogram[phfx+'Nref'] = nobs
1874            Nobs += nobs
1875            if dlg:
1876                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1877            M = np.concatenate((M,wtFactor*df))
1878# end of HKLF processing
1879    Histograms['sumwYo'] = SumwYo
1880    Histograms['Nobs'] = Nobs
1881    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
1882    if dlg:
1883        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
1884        if not GoOn:
1885            parmDict['saved values'] = values
1886            dlg.Destroy()
1887            raise Exception         #Abort!!
1888    pDict,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1889    if np.any(pVals):
1890        pSum = np.sum(pWt*pVals**2)
1891        print 'Penalty function: %.3f on %d terms'%(pSum,len(pVals))
1892        Nobs += len(pVals)
1893        M = np.concatenate((M,np.sqrt(pWt)*pVals))
1894    return M
1895                       
Note: See TracBrowser for help on using the repository browser.