source: trunk/GSASIIstrMath.py @ 1175

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

add surface roughness (Surotti model) to Bragg-Brentano sample parameters
seems to work OK but I don't have a good test data set.

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