source: trunk/GSASIIstrMath.py @ 1111

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

small change to FF calc in GSASIIstrMath.py
add a bit of doc test to testDeriv.py

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