source: trunk/GSASIIstrMath.py @ 1136

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

block StructureFactor2 calcs - less memory

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