source: trunk/GSASIIstrMath.py @ 1113

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

fix Babinet derivatives for SC & powder

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