source: trunk/GSASIIstrMath.py @ 1107

Last change on this file since 1107 was 1107, checked in by vondreele, 8 years ago

make reflist np.array
fix bug - hfx+'Type' needs to be in parmDict

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