source: trunk/GSASIIstrMath.py @ 1322

Last change on this file since 1322 was 1322, checked in by toby, 8 years ago

rework pseudovars and parametric fitting; remove step size for derivatives from expression objects

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 105.5 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2014-05-05 03:30:16 +0000 (Mon, 05 May 2014) $
8# $Author: toby $
9# $Revision: 1322 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1322 2014-05-05 03:30:16Z 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: 1322 $")
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    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.?
530    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
531   
532def StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
533    ''' Compute structure factors for all h,k,l for phase
534    puts the result, F^2, in each ref[8] in refList
535    input:
536   
537    :param dict refDict: where
538        'RefList' list where each ref = h,k,l,m,d,...
539        'FF' dict of form factors - filed in below
540    :param np.array G:      reciprocal metric tensor
541    :param str pfx:    phase id string
542    :param dict SGData: space group info. dictionary output from SpcGroup
543    :param dict calcControls:
544    :param dict ParmDict:
545
546    '''       
547    twopi = 2.0*np.pi
548    twopisq = 2.0*np.pi**2
549    phfx = pfx.split(':')[0]+hfx
550    ast = np.sqrt(np.diag(G))
551    Mast = twopisq*np.multiply.outer(ast,ast)
552    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
553    SGT = np.array([ops[1] for ops in SGData['SGOps']])
554    FFtables = calcControls['FFtables']
555    BLtables = calcControls['BLtables']
556    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
557    FF = np.zeros(len(Tdata))
558    if 'N' in calcControls[hfx+'histType']:
559        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
560    else:
561        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
562        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
563    Uij = np.array(G2lat.U6toUij(Uijdata))
564    bij = Mast*Uij.T
565    if not len(refDict['FF']):
566        if 'N' in calcControls[hfx+'histType']:
567            dat = G2el.getBLvalues(BLtables)
568        else:
569            dat = G2el.getFFvalues(FFtables,0.)       
570        refDict['FF']['El'] = dat.keys()
571        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
572    for iref,refl in enumerate(refDict['RefList']):
573        fbs = np.array([0,0])
574        H = refl[:3]
575        SQ = 1./(2.*refl[4])**2
576        SQfactor = 4.0*SQ*twopisq
577        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
578        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
579            if 'N' in calcControls[hfx+'histType']:
580                dat = G2el.getBLvalues(BLtables)
581                refDict['FF']['FF'][iref] = dat.values()
582            else:       #'X'
583                dat = G2el.getFFvalues(FFtables,SQ)
584                refDict['FF']['FF'][iref] = dat.values()
585        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
586        FF = refDict['FF']['FF'][iref][Tindx]
587        Uniq = np.inner(H,SGMT)
588        Phi = np.inner(H,SGT)
589        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
590        sinp = np.sin(phase)
591        cosp = np.cos(phase)
592        biso = -SQfactor*Uisodata
593        Tiso = np.where(biso<1.,np.exp(biso),1.0)
594        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
595        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
596        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
597        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
598        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
599        if not SGData['SGInv']:
600            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
601            fbs = np.sum(np.sum(fb,axis=1),axis=1)
602        fasq = fas**2
603        fbsq = fbs**2        #imaginary
604        refl[9] = np.sum(fasq)+np.sum(fbsq)
605        refl[10] = atan2d(fbs[0],fas[0])
606   
607def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
608    ''' Compute structure factors for all h,k,l for phase
609    puts the result, F^2, in each ref[8] in refList
610    input:
611   
612    :param dict refDict: where
613        'RefList' list where each ref = h,k,l,m,d,...
614        'FF' dict of form factors - filed in below
615    :param np.array G:      reciprocal metric tensor
616    :param str pfx:    phase id string
617    :param dict SGData: space group info. dictionary output from SpcGroup
618    :param dict calcControls:
619    :param dict ParmDict:
620
621    '''       
622    twopi = 2.0*np.pi
623    twopisq = 2.0*np.pi**2
624    phfx = pfx.split(':')[0]+hfx
625    ast = np.sqrt(np.diag(G))
626    Mast = twopisq*np.multiply.outer(ast,ast)
627    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
628    SGT = np.array([ops[1] for ops in SGData['SGOps']])
629    FFtables = calcControls['FFtables']
630    BLtables = calcControls['BLtables']
631    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
632    FF = np.zeros(len(Tdata))
633    if 'N' in calcControls[hfx+'histType']:
634        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
635    else:
636        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
637        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
638    Uij = np.array(G2lat.U6toUij(Uijdata))
639    bij = Mast*Uij.T
640    blkSize = 100       #no. of reflections in a block
641    nRef = refDict['RefList'].shape[0]
642    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
643        if 'N' in calcControls[hfx+'histType']:
644            dat = G2el.getBLvalues(BLtables)
645            refDict['FF']['El'] = dat.keys()
646            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
647        else:       #'X'
648            dat = G2el.getFFvalues(FFtables,0.)
649            refDict['FF']['El'] = dat.keys()
650            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
651            for iref,ref in enumerate(refDict['RefList']):
652                SQ = 1./(2.*ref[4])**2
653                dat = G2el.getFFvalues(FFtables,SQ)
654                refDict['FF']['FF'][iref] *= dat.values()
655#reflection processing begins here - big arrays!
656    iBeg = 0           
657    while iBeg < nRef:
658        iFin = min(iBeg+blkSize,nRef)
659        refl = refDict['RefList'][iBeg:iFin]
660        H = refl.T[:3]
661        SQ = 1./(2.*refl.T[4])**2
662        SQfactor = 4.0*SQ*twopisq
663        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
664        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
665        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
666        Uniq = np.reshape(np.inner(H.T,SGMT),(-1,3))
667        Phi = np.inner(H.T,SGT).flatten()
668        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T)+Phi[:,np.newaxis])
669        sinp = np.sin(phase)
670        cosp = np.cos(phase)
671        biso = -SQfactor*Uisodata[:,np.newaxis]
672        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
673        HbH = -np.sum(Uniq.T*np.inner(bij,Uniq),axis=1)
674        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
675        Tcorr = Tiso*Tuij*Mdata*Fdata/len(SGMT)
676        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])
677        fa = np.reshape(fa,(2,len(refl),len(SGT),len(Mdata)))
678        fas = np.sum(np.sum(fa,axis=2),axis=2)        #real
679        fbs = np.zeros_like(fas)
680        if not SGData['SGInv']:
681            fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])
682            fb = np.reshape(fb,(2,len(refl),len(SGT),len(Mdata)))
683            fbs = np.sum(np.sum(fb,axis=2),axis=2)
684        fasq = fas**2
685        fbsq = fbs**2        #imaginary
686        refl.T[9] = np.sum(fasq,axis=0)+np.sum(fbsq,axis=0)
687        refl.T[10] = atan2d(fbs[0],fas[0])
688        iBeg += blkSize
689   
690def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
691    'Needs a doc string'
692    twopi = 2.0*np.pi
693    twopisq = 2.0*np.pi**2
694    phfx = pfx.split(':')[0]+hfx
695    ast = np.sqrt(np.diag(G))
696    Mast = twopisq*np.multiply.outer(ast,ast)
697    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
698    SGT = np.array([ops[1] for ops in SGData['SGOps']])
699    FFtables = calcControls['FFtables']
700    BLtables = calcControls['BLtables']
701    nRef = len(refDict['RefList'])
702    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
703    mSize = len(Mdata)
704    FF = np.zeros(len(Tdata))
705    if 'N' in calcControls[hfx+'histType']:
706        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
707    else:
708        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
709        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
710    Uij = np.array(G2lat.U6toUij(Uijdata))
711    bij = Mast*Uij.T
712    dFdvDict = {}
713    dFdfr = np.zeros((nRef,mSize))
714    dFdx = np.zeros((nRef,mSize,3))
715    dFdui = np.zeros((nRef,mSize))
716    dFdua = np.zeros((nRef,mSize,6))
717    dFdbab = np.zeros((nRef,2))
718    for iref,refl in enumerate(refDict['RefList']):
719        H = np.array(refl[:3])
720        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
721        SQfactor = 8.0*SQ*np.pi**2
722        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
723        Bab = parmDict[phfx+'BabA']*dBabdA
724        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
725        FF = refDict['FF']['FF'][iref].T[Tindx]
726#        FF = [refDict['FF'][iref][El] for El in Tdata]         
727        Uniq = np.inner(H,SGMT)
728        Phi = np.inner(H,SGT)
729        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
730        sinp = np.sin(phase)
731        cosp = np.cos(phase)
732        occ = Mdata*Fdata/len(Uniq)
733        biso = -SQfactor*Uisodata
734        Tiso = np.where(biso<1.,np.exp(biso),1.0)
735        HbH = -np.inner(H,np.inner(bij,H))
736        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
737        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
738        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
739        Tcorr = Tiso*Tuij
740        fot = (FF+FP-Bab)*occ*Tcorr
741        fotp = FPP*occ*Tcorr
742        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
743        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
744       
745        fas = np.sum(np.sum(fa,axis=1),axis=1)
746        fbs = np.sum(np.sum(fb,axis=1),axis=1)
747        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
748        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
749        #sum below is over Uniq
750        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
751        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
752        dfadui = np.sum(-SQfactor*fa,axis=2)
753        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
754        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
755        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
756        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
757        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
758        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
759        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
760        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
761        if not SGData['SGInv']:
762            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
763            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
764            dfbdui = np.sum(-SQfactor*fb,axis=2)
765            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
766            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
767            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
768            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
769            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
770            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
771            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
772        #loop over atoms - each dict entry is list of derivatives for all the reflections
773    for i in range(len(Mdata)):     
774        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
775        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
776        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
777        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
778        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
779        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
780        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
781        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
782        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
783        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
784        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
785    dFdvDict[pfx+'BabA'] = dFdbab.T[0]
786    dFdvDict[pfx+'BabU'] = dFdbab.T[1]
787    return dFdvDict
788   
789def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
790    ''' Single crystal extinction function; puts correction in ref[13] and returns
791    corrections needed for derivatives
792    '''
793    ref[11] = 1.0
794    dervCor = 1.0
795    dervDict = {}
796    if calcControls[phfx+'EType'] != 'None':
797        cos2T = 1.0-0.5*(parmDict[hfx+'Lam']/ref[4])**2         #cos(2theta)
798        if 'SXC' in parmDict[hfx+'Type']:
799            AV = 7.9406e5/parmDict[pfx+'Vol']**2
800            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
801            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
802        elif 'SNT' in parmDict[hfx+'Type']:
803            AV = 1.e7/parmDict[pfx+'Vol']**2
804            PL = 1./(4.*refl[4]**2)
805            P12 = 1.0
806        elif 'SNC' in parmDict[hfx+'Type']:
807            AV = 1.e7/parmDict[pfx+'Vol']**2
808            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
809            P12 = 1.0
810           
811        PLZ = AV*P12*parmDict[hfx+'Lam']**2*ref[7]
812        if 'Primary' in calcControls[phfx+'EType']:
813            PLZ *= 1.5
814        else:
815            PLZ *= calcControls[phfx+'Tbar']
816                       
817        if 'Primary' in calcControls[phfx+'EType']:
818            PSIG = parmDict[phfx+'Ep']
819        elif 'I & II' in calcControls[phfx+'EType']:
820            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
821        elif 'Type II' in calcControls[phfx+'EType']:
822            PSIG = parmDict[phfx+'Es']
823        else:       # 'Secondary Type I'
824            PSIG = parmDict[phfx+'Eg']/PL
825           
826        AG = 0.58+0.48*cos2T+0.24*cos2T**2
827        AL = 0.025+0.285*cos2T
828        BG = 0.02-0.025*cos2T
829        BL = 0.15-0.2*(0.75-cos2T)**2
830        if cos2T < 0.:
831            BL = -0.45*cos2T
832        CG = 2.
833        CL = 2.
834        PF = PLZ*PSIG
835       
836        if 'Gaussian' in calcControls[phfx+'EApprox']:
837            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
838            extCor = np.sqrt(PF4)
839            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
840        else:
841            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
842            extCor = np.sqrt(PF4)
843            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
844
845        dervCor = (1.+PF)*PF3
846        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
847            dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
848        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
849            dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
850        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
851            dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
852               
853        ref[11] = 1./extCor
854    return dervCor,dervDict
855   
856def Dict2Values(parmdict, varylist):
857    '''Use before call to leastsq to setup list of values for the parameters
858    in parmdict, as selected by key in varylist'''
859    return [parmdict[key] for key in varylist] 
860   
861def Values2Dict(parmdict, varylist, values):
862    ''' Use after call to leastsq to update the parameter dictionary with
863    values corresponding to keys in varylist'''
864    parmdict.update(zip(varylist,values))
865   
866def GetNewCellParms(parmDict,varyList):
867    'Needs a doc string'
868    newCellDict = {}
869    Anames = ['A'+str(i) for i in range(6)]
870    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
871    for item in varyList:
872        keys = item.split(':')
873        if keys[2] in Ddict:
874            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
875            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
876            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
877    return newCellDict          # is e.g. {'0::D11':A0+D11}
878   
879def ApplyXYZshifts(parmDict,varyList):
880    '''
881    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
882   
883    :param dict parmDict: parameter dictionary
884    :param list varyList: list of variables
885    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
886
887    '''
888    newAtomDict = {}
889    for item in parmDict:
890        if 'dA' in item:
891            parm = ''.join(item.split('d'))
892            parmDict[parm] += parmDict[item]
893            newAtomDict[item] = [parm,parmDict[parm]]
894    return newAtomDict
895   
896def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
897    'Spherical harmonics texture'
898    IFCoup = 'Bragg' in calcControls[hfx+'instType']
899    odfCor = 1.0
900    H = refl[:3]
901    cell = G2lat.Gmat2cell(g)
902    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
903    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
904    phi,beta = G2lat.CrsAng(H,cell,SGData)
905    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
906    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
907    for item in SHnames:
908        L,M,N = eval(item.strip('C'))
909        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
910        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
911        Lnorm = G2lat.Lnorm(L)
912        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
913    return odfCor
914   
915def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
916    'Spherical harmonics texture derivatives'
917    FORPI = 4.0*np.pi
918    IFCoup = 'Bragg' in calcControls[hfx+'instType']
919    odfCor = 1.0
920    dFdODF = {}
921    dFdSA = [0,0,0]
922    H = refl[:3]
923    cell = G2lat.Gmat2cell(g)
924    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
925    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
926    phi,beta = G2lat.CrsAng(H,cell,SGData)
927    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
928    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
929    for item in SHnames:
930        L,M,N = eval(item.strip('C'))
931        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
932        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
933        Lnorm = G2lat.Lnorm(L)
934        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
935        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
936        for i in range(3):
937            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
938    return odfCor,dFdODF,dFdSA
939   
940def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
941    'spherical harmonics preferred orientation (cylindrical symmetry only)'
942    odfCor = 1.0
943    H = refl[:3]
944    cell = G2lat.Gmat2cell(g)
945    Sangl = [0.,0.,0.]
946    if 'Bragg' in calcControls[hfx+'instType']:
947        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
948        IFCoup = True
949    else:
950        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
951        IFCoup = False
952    phi,beta = G2lat.CrsAng(H,cell,SGData)
953    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
954    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
955    for item in SHnames:
956        L,N = eval(item.strip('C'))
957        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
958        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
959    return np.squeeze(odfCor)
960   
961def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
962    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
963    FORPI = 12.5663706143592
964    odfCor = 1.0
965    dFdODF = {}
966    H = refl[:3]
967    cell = G2lat.Gmat2cell(g)
968    Sangl = [0.,0.,0.]
969    if 'Bragg' in calcControls[hfx+'instType']:
970        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
971        IFCoup = True
972    else:
973        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
974        IFCoup = False
975    phi,beta = G2lat.CrsAng(H,cell,SGData)
976    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
977    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
978    for item in SHnames:
979        L,N = eval(item.strip('C'))
980        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
981        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
982        dFdODF[phfx+item] = Kcsl*Lnorm
983    return odfCor,dFdODF
984   
985def GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
986    'Needs a doc string'
987    POcorr = 1.0
988    if calcControls[phfx+'poType'] == 'MD':
989        MD = parmDict[phfx+'MD']
990        if MD != 1.0:
991            MDAxis = calcControls[phfx+'MDAxis']
992            sumMD = 0
993            for H in uniq:           
994                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
995                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
996                sumMD += A**3
997            POcorr = sumMD/len(uniq)
998    else:   #spherical harmonics
999        if calcControls[phfx+'SHord']:
1000            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1001    return POcorr
1002   
1003def GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1004    'Needs a doc string'
1005    POcorr = 1.0
1006    POderv = {}
1007    if calcControls[phfx+'poType'] == 'MD':
1008        MD = parmDict[phfx+'MD']
1009        MDAxis = calcControls[phfx+'MDAxis']
1010        sumMD = 0
1011        sumdMD = 0
1012        for H in uniq:           
1013            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1014            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1015            sumMD += A**3
1016            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1017        POcorr = sumMD/len(uniq)
1018        POderv[phfx+'MD'] = sumdMD/len(uniq)
1019    else:   #spherical harmonics
1020        if calcControls[phfx+'SHord']:
1021            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1022    return POcorr,POderv
1023   
1024def GetAbsorb(refl,hfx,calcControls,parmDict):
1025    'Needs a doc string'
1026    if 'Debye' in calcControls[hfx+'instType']:
1027        return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1028    else:
1029        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1030   
1031def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
1032    'Needs a doc string'
1033    if 'Debye' in calcControls[hfx+'instType']:
1034        return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
1035    else:
1036        return G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
1037   
1038def GetIntensityCorr(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1039    'Needs a doc string'    #need powder extinction!
1040    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1041    if 'X' in parmDict[hfx+'Type']:
1042        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1043    Icorr *= GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1044    if pfx+'SHorder' in parmDict:
1045        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1046    Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
1047    refl[11] = Icorr       
1048   
1049def GetIntensityDerv(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1050    'Needs a doc string'    #need powder extinction derivs!
1051    dIdsh = 1./parmDict[hfx+'Scale']
1052    dIdsp = 1./parmDict[phfx+'Scale']
1053    if 'X' in parmDict[hfx+'Type']:
1054        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1055        dIdPola /= pola
1056    else:       #'N'
1057        dIdPola = 0.0
1058    POcorr,dIdPO = GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1059    for iPO in dIdPO:
1060        dIdPO[iPO] /= POcorr
1061    dFdODF = {}
1062    dFdSA = [0,0,0]
1063    if pfx+'SHorder' in parmDict:
1064        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1065        for iSH in dFdODF:
1066            dFdODF[iSH] /= odfCor
1067        for i in range(3):
1068            dFdSA[i] /= odfCor
1069    dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
1070    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
1071       
1072def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
1073    'Needs a doc string'
1074    costh = cosd(refl[5]/2.)
1075    #crystallite size
1076    if calcControls[phfx+'SizeType'] == 'isotropic':
1077        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1078    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1079        H = np.array(refl[:3])
1080        P = np.array(calcControls[phfx+'SizeAxis'])
1081        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1082        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1083        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1084    else:           #ellipsoidal crystallites
1085        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1086        H = np.array(refl[:3])
1087        lenR = G2pwd.ellipseSize(H,Sij,GB)
1088        Sgam = 1.8*wave/(np.pi*costh*lenR)
1089    #microstrain               
1090    if calcControls[phfx+'MustrainType'] == 'isotropic':
1091        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1092    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1093        H = np.array(refl[:3])
1094        P = np.array(calcControls[phfx+'MustrainAxis'])
1095        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1096        Si = parmDict[phfx+'Mustrain;i']
1097        Sa = parmDict[phfx+'Mustrain;a']
1098        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1099    else:       #generalized - P.W. Stephens model
1100        pwrs = calcControls[phfx+'MuPwrs']
1101        sum = 0
1102        for i,pwr in enumerate(pwrs):
1103            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1104        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
1105    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1106    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1107    sig /= ateln2
1108    return sig,gam
1109       
1110def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
1111    'Needs a doc string'
1112    gamDict = {}
1113    sigDict = {}
1114    costh = cosd(refl[5]/2.)
1115    tanth = tand(refl[5]/2.)
1116    #crystallite size derivatives
1117    if calcControls[phfx+'SizeType'] == 'isotropic':
1118        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1119        gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1120        sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1121    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1122        H = np.array(refl[:3])
1123        P = np.array(calcControls[phfx+'SizeAxis'])
1124        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1125        Si = parmDict[phfx+'Size;i']
1126        Sa = parmDict[phfx+'Size;a']
1127        gami = (1.8*wave/np.pi)/(Si*Sa)
1128        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1129        Sgam = gami*sqtrm
1130        gam = Sgam/costh
1131        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
1132        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
1133        gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1134        gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1135        sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1136        sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1137    else:           #ellipsoidal crystallites
1138        const = 1.8*wave/(np.pi*costh)
1139        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1140        H = np.array(refl[:3])
1141        lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1142        Sgam = 1.8*wave/(np.pi*costh*lenR)
1143        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1144            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1145            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1146    gamDict[phfx+'Size;mx'] = Sgam
1147    sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1148           
1149    #microstrain derivatives               
1150    if calcControls[phfx+'MustrainType'] == 'isotropic':
1151        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
1152        gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1153        sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1154    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1155        H = np.array(refl[:3])
1156        P = np.array(calcControls[phfx+'MustrainAxis'])
1157        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1158        Si = parmDict[phfx+'Mustrain;i']
1159        Sa = parmDict[phfx+'Mustrain;a']
1160        gami = 0.018*Si*Sa*tanth/np.pi
1161        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1162        Mgam = gami/sqtrm
1163        dsi = -gami*Si*cosP**2/sqtrm**3
1164        dsa = -gami*Sa*sinP**2/sqtrm**3
1165        gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1166        gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1167        sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1168        sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1169    else:       #generalized - P.W. Stephens model
1170        pwrs = calcControls[phfx+'MuPwrs']
1171        const = 0.018*refl[4]**2*tanth
1172        sum = 0
1173        for i,pwr in enumerate(pwrs):
1174            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1175            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
1176            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
1177            sigDict[phfx+'Mustrain:'+str(i)] = \
1178                2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1179        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
1180        for i in range(len(pwrs)):
1181            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
1182    gamDict[phfx+'Mustrain;mx'] = Mgam
1183    sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1184    return sigDict,gamDict
1185       
1186def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1187    'Needs a doc string'
1188    h,k,l = refl[:3]
1189    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1190    d = np.sqrt(dsq)
1191
1192    refl[4] = d
1193    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1194    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1195    if 'Bragg' in calcControls[hfx+'instType']:
1196        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1197            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1198    else:               #Debye-Scherrer - simple but maybe not right
1199        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1200    return pos
1201
1202def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1203    'Needs a doc string'
1204    dpr = 180./np.pi
1205    h,k,l = refl[:3]
1206    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1207    dst = np.sqrt(dstsq)
1208    pos = refl[5]-parmDict[hfx+'Zero']
1209    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1210    dpdw = const*dst
1211    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1212    dpdA *= const*wave/(2.0*dst)
1213    dpdZ = 1.0
1214    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1215    if 'Bragg' in calcControls[hfx+'instType']:
1216        dpdSh = -4.*const*cosd(pos/2.0)
1217        dpdTr = -const*sind(pos)*100.0
1218        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1219    else:               #Debye-Scherrer - simple but maybe not right
1220        dpdXd = -const*cosd(pos)
1221        dpdYd = -const*sind(pos)
1222        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1223           
1224def GetHStrainShift(refl,SGData,phfx,parmDict):
1225    'Needs a doc string'
1226    laue = SGData['SGLaue']
1227    uniq = SGData['SGUniq']
1228    h,k,l = refl[:3]
1229    if laue in ['m3','m3m']:
1230        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1231            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1232    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1233        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1234    elif laue in ['3R','3mR']:
1235        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1236    elif laue in ['4/m','4/mmm']:
1237        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1238    elif laue in ['mmm']:
1239        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1240    elif laue in ['2/m']:
1241        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1242        if uniq == 'a':
1243            Dij += parmDict[phfx+'D23']*k*l
1244        elif uniq == 'b':
1245            Dij += parmDict[phfx+'D13']*h*l
1246        elif uniq == 'c':
1247            Dij += parmDict[phfx+'D12']*h*k
1248    else:
1249        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1250            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1251    return -Dij*refl[4]**2*tand(refl[5]/2.0)
1252           
1253def GetHStrainShiftDerv(refl,SGData,phfx):
1254    'Needs a doc string'
1255    laue = SGData['SGLaue']
1256    uniq = SGData['SGUniq']
1257    h,k,l = refl[:3]
1258    if laue in ['m3','m3m']:
1259        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1260            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1261    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1262        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1263    elif laue in ['3R','3mR']:
1264        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1265    elif laue in ['4/m','4/mmm']:
1266        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1267    elif laue in ['mmm']:
1268        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1269    elif laue in ['2/m']:
1270        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1271        if uniq == 'a':
1272            dDijDict[phfx+'D23'] = k*l
1273        elif uniq == 'b':
1274            dDijDict[phfx+'D13'] = h*l
1275        elif uniq == 'c':
1276            dDijDict[phfx+'D12'] = h*k
1277            names.append()
1278    else:
1279        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1280            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1281    for item in dDijDict:
1282        dDijDict[item] *= -refl[4]**2*tand(refl[5]/2.0)
1283    return dDijDict
1284   
1285def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1286    'Needs a doc string'
1287    histoList = Histograms.keys()
1288    histoList.sort()
1289    for histogram in histoList:
1290        if 'PWDR' in histogram[:4]:
1291            Histogram = Histograms[histogram]
1292            hId = Histogram['hId']
1293            hfx = ':%d:'%(hId)
1294            Limits = calcControls[hfx+'Limits']
1295            shl = max(parmDict[hfx+'SH/L'],0.0005)
1296            Ka2 = False
1297            kRatio = 0.0
1298            if hfx+'Lam1' in parmDict.keys():
1299                Ka2 = True
1300                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1301                kRatio = parmDict[hfx+'I(L2)/I(L1)']
1302            x,y,w,yc,yb,yd = Histogram['Data']
1303            xB = np.searchsorted(x,Limits[0])
1304            xF = np.searchsorted(x,Limits[1])
1305            ymb = np.array(y-yb)
1306            ymb = np.where(ymb,ymb,1.0)
1307            ycmb = np.array(yc-yb)
1308            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
1309            refLists = Histogram['Reflection Lists']
1310            for phase in refLists:
1311                Phase = Phases[phase]
1312                pId = Phase['pId']
1313                phfx = '%d:%d:'%(pId,hId)
1314                refDict = refLists[phase]
1315                sumFo = 0.0
1316                sumdF = 0.0
1317                sumFosq = 0.0
1318                sumdFsq = 0.0
1319                for refl in refDict['RefList']:
1320                    if 'C' in calcControls[hfx+'histType']:
1321                        yp = np.zeros_like(yb)
1322                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1323                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
1324                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
1325                        iFin2 = iFin
1326                        yp[iBeg:iFin] = refl[11]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1327                        if Ka2:
1328                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1329                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1330                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
1331                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
1332                            if not iBeg2+iFin2:       #peak below low limit - skip peak
1333                                continue
1334                            elif not iBeg2-iFin2:     #peak above high limit - done
1335                                break
1336                            yp[iBeg2:iFin2] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
1337                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11]*(1.+kRatio)),0.0))
1338                    elif 'T' in calcControls[hfx+'histType']:
1339                        print 'TOF Undefined at present'
1340                        raise Exception    #no TOF yet
1341                    Fo = np.sqrt(np.abs(refl[8]))
1342                    Fc = np.sqrt(np.abs(refl[9]))
1343                    sumFo += Fo
1344                    sumFosq += refl[8]**2
1345                    sumdF += np.abs(Fo-Fc)
1346                    sumdFsq += (refl[8]-refl[9])**2
1347                Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1348                Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1349                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
1350                Histogram['Residuals']['hId'] = hId
1351        elif 'HKLF' in histogram[:4]:
1352            Histogram = Histograms[histogram]
1353            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
1354               
1355def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1356    'Needs a doc string'
1357   
1358    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1359        U = parmDict[hfx+'U']
1360        V = parmDict[hfx+'V']
1361        W = parmDict[hfx+'W']
1362        X = parmDict[hfx+'X']
1363        Y = parmDict[hfx+'Y']
1364        tanPos = tand(refl[5]/2.0)
1365        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
1366        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1367        sig = max(0.001,sig)
1368        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
1369        gam = max(0.001,gam)
1370        return sig,gam
1371               
1372    hId = Histogram['hId']
1373    hfx = ':%d:'%(hId)
1374    bakType = calcControls[hfx+'bakType']
1375    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1376    yc = np.zeros_like(yb)
1377       
1378    if 'C' in calcControls[hfx+'histType']:   
1379        shl = max(parmDict[hfx+'SH/L'],0.002)
1380        Ka2 = False
1381        if hfx+'Lam1' in parmDict.keys():
1382            wave = parmDict[hfx+'Lam1']
1383            Ka2 = True
1384            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1385            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1386        else:
1387            wave = parmDict[hfx+'Lam']
1388    else:
1389        print 'TOF Undefined at present'
1390        raise ValueError
1391    for phase in Histogram['Reflection Lists']:
1392        refDict = Histogram['Reflection Lists'][phase]
1393        Phase = Phases[phase]
1394        pId = Phase['pId']
1395        pfx = '%d::'%(pId)
1396        phfx = '%d:%d:'%(pId,hId)
1397        hfx = ':%d:'%(hId)
1398        SGData = Phase['General']['SGData']
1399        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1400        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]     #Do I want to modify by Dij?
1401        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1402        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1403        Vst = np.sqrt(nl.det(G))    #V*
1404        if not Phase['General'].get('doPawley'):
1405            time0 = time.time()
1406            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1407#            print 'sf calc time: %.3fs'%(time.time()-time0)
1408        time0 = time.time()
1409        for iref,refl in enumerate(refDict['RefList']):
1410            if 'C' in calcControls[hfx+'histType']:
1411                h,k,l = refl[:3]
1412                Uniq = np.inner(refl[:3],SGMT)
1413                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1414                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1415                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1416                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
1417                GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[11]
1418                refl[11] *= Vst*Lorenz
1419                if Phase['General'].get('doPawley'):
1420                    try:
1421                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1422                        refl[9] = parmDict[pInd]
1423                    except KeyError:
1424#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1425                        continue
1426                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1427                iBeg = np.searchsorted(x,refl[5]-fmin)
1428                iFin = np.searchsorted(x,refl[5]+fmax)
1429                if not iBeg+iFin:       #peak below low limit - skip peak
1430                    continue
1431                elif not iBeg-iFin:     #peak above high limit - done
1432                    break
1433                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
1434                if Ka2:
1435                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1436                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1437                    iBeg = np.searchsorted(x,pos2-fmin)
1438                    iFin = np.searchsorted(x,pos2+fmax)
1439                    if not iBeg+iFin:       #peak below low limit - skip peak
1440                        continue
1441                    elif not iBeg-iFin:     #peak above high limit - done
1442                        return yc,yb
1443                    yc[iBeg:iFin] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))        #and here
1444            elif 'T' in calcControls[hfx+'histType']:
1445                print 'TOF Undefined at present'
1446                raise Exception    #no TOF yet
1447#        print 'profile calc time: %.3fs'%(time.time()-time0)
1448    return yc,yb
1449   
1450def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1451    'Needs a doc string'
1452   
1453    def cellVaryDerv(pfx,SGData,dpdA): 
1454        if SGData['SGLaue'] in ['-1',]:
1455            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1456                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1457        elif SGData['SGLaue'] in ['2/m',]:
1458            if SGData['SGUniq'] == 'a':
1459                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1460            elif SGData['SGUniq'] == 'b':
1461                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1462            else:
1463                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1464        elif SGData['SGLaue'] in ['mmm',]:
1465            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1466        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1467            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1468        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1469            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1470        elif SGData['SGLaue'] in ['3R', '3mR']:
1471            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1472        elif SGData['SGLaue'] in ['m3m','m3']:
1473            return [[pfx+'A0',dpdA[0]]]
1474           
1475    # create a list of dependent variables and set up a dictionary to hold their derivatives
1476    dependentVars = G2mv.GetDependentVars()
1477    depDerivDict = {}
1478    for j in dependentVars:
1479        depDerivDict[j] = np.zeros(shape=(len(x)))
1480    #print 'dependent vars',dependentVars
1481    lenX = len(x)               
1482    hId = Histogram['hId']
1483    hfx = ':%d:'%(hId)
1484    bakType = calcControls[hfx+'bakType']
1485    dMdv = np.zeros(shape=(len(varylist),len(x)))
1486    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1487    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1488        bBpos =varylist.index(hfx+'Back:0')
1489        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1490    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1491    for name in varylist:
1492        if 'Debye' in name:
1493            id = int(name.split(':')[-1])
1494            parm = name[:int(name.rindex(':'))]
1495            ip = names.index(parm)
1496            dMdv[varylist.index(name)] = dMddb[3*id+ip]
1497    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1498    for name in varylist:
1499        if 'BkPk' in name:
1500            parm,id = name.split(';')
1501            id = int(id)
1502            if parm in names:
1503                ip = names.index(parm)
1504                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1505    cw = np.diff(x)
1506    cw = np.append(cw,cw[-1])
1507    if 'C' in calcControls[hfx+'histType']:   
1508        shl = max(parmDict[hfx+'SH/L'],0.002)
1509        Ka2 = False
1510        if hfx+'Lam1' in parmDict.keys():
1511            wave = parmDict[hfx+'Lam1']
1512            Ka2 = True
1513            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1514            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1515        else:
1516            wave = parmDict[hfx+'Lam']
1517    else:
1518        print 'TOF Undefined at present'
1519        raise ValueError
1520    for phase in Histogram['Reflection Lists']:
1521        refDict = Histogram['Reflection Lists'][phase]
1522        Phase = Phases[phase]
1523        SGData = Phase['General']['SGData']
1524        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1525        pId = Phase['pId']
1526        pfx = '%d::'%(pId)
1527        phfx = '%d:%d:'%(pId,hId)
1528        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]     #And modify here by Dij?
1529        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1530        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1531        if not Phase['General'].get('doPawley'):
1532            time0 = time.time()
1533            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1534#            print 'sf-derv time %.3fs'%(time.time()-time0)
1535            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1536        time0 = time.time()
1537        for iref,refl in enumerate(refDict['RefList']):
1538            if 'C' in calcControls[hfx+'histType']:        #CW powder
1539                h,k,l = refl[:3]
1540                Uniq = np.inner(refl[:3],SGMT)
1541                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1542                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1543                iBeg = np.searchsorted(x,refl[5]-fmin)
1544                iFin = np.searchsorted(x,refl[5]+fmax)
1545                if not iBeg+iFin:       #peak below low limit - skip peak
1546                    continue
1547                elif not iBeg-iFin:     #peak above high limit - done
1548                    break
1549                pos = refl[5]
1550                tanth = tand(pos/2.0)
1551                costh = cosd(pos/2.0)
1552                lenBF = iFin-iBeg
1553                dMdpk = np.zeros(shape=(6,lenBF))
1554                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1555                for i in range(5):
1556                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11]*refl[9]*dMdipk[i]
1557                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1558                if Ka2:
1559                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1560                    iBeg2 = np.searchsorted(x,pos2-fmin)
1561                    iFin2 = np.searchsorted(x,pos2+fmax)
1562                    if iBeg2-iFin2:
1563                        lenBF2 = iFin2-iBeg2
1564                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1565                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1566                        for i in range(5):
1567                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11]*refl[9]*kRatio*dMdipk2[i]
1568                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11]*dMdipk2[0]
1569                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1570                if Phase['General'].get('doPawley'):
1571                    dMdpw = np.zeros(len(x))
1572                    try:
1573                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1574                        idx = varylist.index(pIdx)
1575                        dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1576                        if Ka2:
1577                            dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1578                        dMdv[idx] = dMdpw
1579                    except: # ValueError:
1580                        pass
1581                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1582                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1583                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1584                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1585                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1586                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1587                    hfx+'DisplaceY':[dpdY,'pos'],}
1588                if 'Bragg' in calcControls[hfx+'instType']:
1589                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
1590                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
1591                else:
1592                    names.update({hfx+'Absorption':[dFdAb,'int'],})
1593                for name in names:
1594                    item = names[name]
1595                    if name in varylist:
1596                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1597                        if Ka2:
1598                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1599                    elif name in dependentVars:
1600                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1601                        if Ka2:
1602                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1603                for iPO in dIdPO:
1604                    if iPO in varylist:
1605                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1606                        if Ka2:
1607                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1608                    elif iPO in dependentVars:
1609                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1610                        if Ka2:
1611                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1612                for i,name in enumerate(['omega','chi','phi']):
1613                    aname = pfx+'SH '+name
1614                    if aname in varylist:
1615                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1616                        if Ka2:
1617                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1618                    elif aname in dependentVars:
1619                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1620                        if Ka2:
1621                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1622                for iSH in dFdODF:
1623                    if iSH in varylist:
1624                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1625                        if Ka2:
1626                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1627                    elif iSH in dependentVars:
1628                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1629                        if Ka2:
1630                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1631                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1632                for name,dpdA in cellDervNames:
1633                    if name in varylist:
1634                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1635                        if Ka2:
1636                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1637                    elif name in dependentVars:
1638                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1639                        if Ka2:
1640                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1641                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1642                for name in dDijDict:
1643                    if name in varylist:
1644                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1645                        if Ka2:
1646                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1647                    elif name in dependentVars:
1648                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1649                        if Ka2:
1650                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1651                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
1652                for name in gamDict:
1653                    if name in varylist:
1654                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1655                        if Ka2:
1656                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1657                    elif name in dependentVars:
1658                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1659                        if Ka2:
1660                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1661                for name in sigDict:
1662                    if name in varylist:
1663                        dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1664                        if Ka2:
1665                            dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1666                    elif name in dependentVars:
1667                        depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1668                        if Ka2:
1669                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1670                for name in ['BabA','BabU']:
1671                    if refl[9]:
1672                        if phfx+name in varylist:
1673                            dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1674                            if Ka2:
1675                                dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
1676                        elif phfx+name in dependentVars:                   
1677                            depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
1678                            if Ka2:
1679                                depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]                 
1680            elif 'T' in calcControls[hfx+'histType']:
1681                print 'TOF Undefined at present'
1682                raise Exception    #no TOF yet
1683            if not Phase['General'].get('doPawley'):
1684                #do atom derivatives -  for RB,F,X & U so far             
1685                corr = dervDict['int']/refl[9]
1686                if Ka2:
1687                    corr2 = dervDict2['int']/refl[9]
1688                for name in varylist+dependentVars:
1689                    if '::RBV;' in name:
1690                        pass
1691                    else:
1692                        try:
1693                            aname = name.split(pfx)[1][:2]
1694                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1695                        except IndexError:
1696                            continue
1697                    if name in varylist:
1698                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1699                        if Ka2:
1700                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1701                    elif name in dependentVars:
1702                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1703                        if Ka2:
1704                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1705#        print 'profile derv time: %.3fs'%(time.time()-time0)
1706    # now process derivatives in constraints
1707    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1708    return dMdv
1709
1710def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1711    '''Loop over histograms and compute derivatives of the fitting
1712    model (M) with respect to all parameters.  Results are returned in
1713    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
1714    n is the number of parameters and m is the number of data
1715    points. This can exceed memory when m gets large. This routine is
1716    used when refinement derivatives are selected as "analtytic
1717    Jacobian" in Controls.
1718
1719    :returns: Jacobian numpy.array dMdv for all histograms concatinated
1720    '''
1721    parmDict.update(zip(varylist,values))
1722    G2mv.Dict2Map(parmDict,varylist)
1723    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1724    nvar = len(varylist)
1725    dMdv = np.empty(0)
1726    histoList = Histograms.keys()
1727    histoList.sort()
1728    for histogram in histoList:
1729        if 'PWDR' in histogram[:4]:
1730            Histogram = Histograms[histogram]
1731            hId = Histogram['hId']
1732            hfx = ':%d:'%(hId)
1733            wtFactor = calcControls[hfx+'wtFactor']
1734            Limits = calcControls[hfx+'Limits']
1735            x,y,w,yc,yb,yd = Histogram['Data']
1736            W = wtFactor*w
1737            xB = np.searchsorted(x,Limits[0])
1738            xF = np.searchsorted(x,Limits[1])
1739            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
1740                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1741        elif 'HKLF' in histogram[:4]:
1742            Histogram = Histograms[histogram]
1743            nobs = Histogram['Residuals']['Nobs']
1744            phase = Histogram['Reflection Lists']
1745            Phase = Phases[phase]
1746            hId = Histogram['hId']
1747            hfx = ':%d:'%(hId)
1748            wtFactor = calcControls[hfx+'wtFactor']
1749            pfx = '%d::'%(Phase['pId'])
1750            phfx = '%d:%d:'%(Phase['pId'],hId)
1751            SGData = Phase['General']['SGData']
1752            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1753            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1754            refDict = Histogram['Data']
1755            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1756            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1757            dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1758            dependentVars = G2mv.GetDependentVars()
1759            depDerivDict = {}
1760            for j in dependentVars:
1761                depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
1762            if calcControls['F**2']:
1763                for iref,ref in enumerate(refDict['RefList']):
1764                    if ref[6] > 0:
1765                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1766                        w = 1.0/ref[6]
1767                        if w*ref[5] >= calcControls['minF/sig']:
1768                            for j,var in enumerate(varylist):
1769                                if var in dFdvDict:
1770                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1771                            for var in dependentVars:
1772                                if var in dFdvDict:
1773                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1774                            if phfx+'Scale' in varylist:
1775                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1776                            elif phfx+'Scale' in dependentVars:
1777                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
1778                            for item in ['Ep','Es','Eg']:
1779                                if phfx+item in varylist:
1780                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1781                                elif phfx+item in dependentVars:
1782                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1783                            for item in ['BabA','BabU']:
1784                                if phfx+item in varylist:
1785                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1786                                elif phfx+item in dependentVars:
1787                                    depDerivDict[phfx+item][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1788            else:
1789                for iref,ref in enumerate(refDict['RefList']):
1790                    if ref[5] > 0.:
1791                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1792                        Fo = np.sqrt(ref[5])
1793                        Fc = np.sqrt(ref[7])
1794                        w = 1.0/ref[6]
1795                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1796                            for j,var in enumerate(varylist):
1797                                if var in dFdvDict:
1798                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1799                            for var in dependentVars:
1800                                if var in dFdvDict:
1801                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1802                            if phfx+'Scale' in varylist:
1803                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1804                            elif phfx+'Scale' in dependentVars:
1805                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor                           
1806                            for item in ['Ep','Es','Eg']:
1807                                if phfx+item in varylist:
1808                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1809                                elif phfx+item in dependentVars:
1810                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1811                            for item in ['BabA','BabU']:
1812                                if phfx+item in varylist:
1813                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1814                                elif phfx+item in dependentVars:
1815                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1816            # now process derivatives in constraints
1817            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1818        else:
1819            continue        #skip non-histogram entries
1820        if len(dMdv):
1821            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
1822        else:
1823            dMdv = np.sqrt(wtFactor)*dMdvh
1824           
1825    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
1826    if np.any(pVals):
1827        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1828        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
1829       
1830    return dMdv
1831
1832def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1833    '''Loop over histograms and compute derivatives of the fitting
1834    model (M) with respect to all parameters.  For each histogram, the
1835    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
1836    number of parameters and m is the number of data points *in the
1837    histogram*. The (n by n) Hessian is computed from each Jacobian
1838    and it is returned.  This routine is used when refinement
1839    derivatives are selected as "analtytic Hessian" in Controls.
1840
1841    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
1842    '''
1843    parmDict.update(zip(varylist,values))
1844    G2mv.Dict2Map(parmDict,varylist)
1845    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1846    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
1847    nvar = len(varylist)
1848    Hess = np.empty(0)
1849    histoList = Histograms.keys()
1850    histoList.sort()
1851    for histogram in histoList:
1852        if 'PWDR' in histogram[:4]:
1853            Histogram = Histograms[histogram]
1854            hId = Histogram['hId']
1855            hfx = ':%d:'%(hId)
1856            wtFactor = calcControls[hfx+'wtFactor']
1857            Limits = calcControls[hfx+'Limits']
1858            x,y,w,yc,yb,yd = Histogram['Data']
1859            W = wtFactor*w
1860            dy = y-yc
1861            xB = np.searchsorted(x,Limits[0])
1862            xF = np.searchsorted(x,Limits[1])
1863            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
1864                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1865            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
1866            Dy = dy[xB:xF][np.newaxis,:]
1867            dMdvh *= Wt
1868            if dlg:
1869                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1870            if len(Hess):
1871                Hess += np.inner(dMdvh,dMdvh)
1872                dMdvh *= Wt*Dy
1873                Vec += np.sum(dMdvh,axis=1)
1874            else:
1875                Hess = np.inner(dMdvh,dMdvh)
1876                dMdvh *= Wt*Dy
1877                Vec = np.sum(dMdvh,axis=1)
1878        elif 'HKLF' in histogram[:4]:
1879            Histogram = Histograms[histogram]
1880            nobs = Histogram['Residuals']['Nobs']
1881            phase = Histogram['Reflection Lists']
1882            Phase = Phases[phase]
1883            hId = Histogram['hId']
1884            hfx = ':%d:'%(hId)
1885            wtFactor = calcControls[hfx+'wtFactor']
1886            pfx = '%d::'%(Phase['pId'])
1887            phfx = '%d:%d:'%(Phase['pId'],hId)
1888            SGData = Phase['General']['SGData']
1889            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1890            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1891            refDict = Histogram['Data']
1892            time0 = time.time()
1893            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
1894            print 'sf-deriv time: %.3f'%(time.time()-time0)
1895            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1896            dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
1897            dependentVars = G2mv.GetDependentVars()
1898            depDerivDict = {}
1899            for j in dependentVars:
1900                depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
1901            wdf = np.zeros(len(refDict['RefList']))
1902            time0 = time.time()
1903            if calcControls['F**2']:
1904                for iref,ref in enumerate(refDict['RefList']):
1905                    if ref[6] > 0:
1906                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1907                        w =  1.0/ref[6]
1908                        if w*ref[5] >= calcControls['minF/sig']:
1909                            wdf[iref] = w*(ref[5]-ref[7])
1910                            for j,var in enumerate(varylist):
1911                                if var in dFdvDict:
1912                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1913                            for var in dependentVars:
1914                                if var in dFdvDict:
1915                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1916                            if phfx+'Scale' in varylist:
1917                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1918                            elif phfx+'Scale' in dependentVars:
1919                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
1920                            for item in ['Ep','Es','Eg']:
1921                                if phfx+item in varylist:
1922                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1923                                elif phfx+item in dependentVars:
1924                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1925                            for item in ['BabA','BabU']:
1926                                if phfx+item in varylist:
1927                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1928                                elif phfx+item in dependentVars:
1929                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1930            else:
1931                for iref,ref in enumerate(refDict['RefList']):
1932                    if ref[5] > 0.:
1933                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1934                        Fo = np.sqrt(ref[5])
1935                        Fc = np.sqrt(ref[7])
1936                        w = 1.0/ref[6]
1937                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1938                            wdf[iref] = 2.0*Fo*w*(Fo-Fc)
1939                            for j,var in enumerate(varylist):
1940                                if var in dFdvDict:
1941                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1942                            for var in dependentVars:
1943                                if var in dFdvDict:
1944                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1945                            if phfx+'Scale' in varylist:
1946                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1947                            elif phfx+'Scale' in dependentVars:
1948                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor                           
1949                            for item in ['Ep','Es','Eg']:
1950                                if phfx+item in varylist:
1951                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1952                                elif phfx+item in dependentVars:
1953                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1954                            for item in ['BabA','BabU']:
1955                                if phfx+item in varylist:
1956                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1957                                elif phfx+item in dependentVars:
1958                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1959            # now process derivatives in constraints
1960            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1961            print 'matrix build time: %.3f'%(time.time()-time0)
1962
1963            if dlg:
1964                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1965            if len(Hess):
1966                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
1967                Hess += wtFactor*np.inner(dMdvh,dMdvh)
1968            else:
1969                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
1970                Hess = wtFactor*np.inner(dMdvh,dMdvh)
1971        else:
1972            continue        #skip non-histogram entries
1973    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
1974    if np.any(pVals):
1975        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1976        Vec += np.sum(dpdv*pWt*pVals,axis=1)
1977        Hess += np.inner(dpdv*pWt,dpdv)
1978    return Vec,Hess
1979
1980def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
1981    'Needs a doc string'
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.