source: trunk/GSASIIstrMath.py @ 1391

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

remove stray prints
Breit-Wigner & Lynn & Seeger neutron anomalous form factors for CW data
Add 'Ext' to HKLF reflection lists
make sure wave & tbar are in TOF HKLF reflection lists
add '_refln_f_sigma' & '_refln_f_squared_sigma' to possible cif reflection columns

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