source: trunk/GSASIIstrMath.py @ 954

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

fix to Hessian LSQ to stop if lam > 10e5!
fix to rigid body quaternion derivatives

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