source: trunk/GSASIIstrMath.py @ 1134

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

add MemoryError? check to StructureFactor2
allow inversion of x & y axes on image plot
these are saved
fix missing stuff in save/restore image controls

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