source: trunk/GSASIIstrMath.py @ 1455

Last change on this file since 1455 was 1455, checked in by vondreele, 7 years ago

fix ddata copy flags so it copies extinction models as well as refine flags
more fixes to derivatives,

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