source: trunk/GSASIIstrMath.py @ 1299

Last change on this file since 1299 was 1299, checked in by toby, 9 years ago

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