source: trunk/GSASIIstrMath.py @ 926

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

split GSASIIstruct.py into 3 pieces

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