source: trunk/GSASIIstrMath.py @ 1456

Last change on this file since 1456 was 1456, checked in by vondreele, 7 years ago

consolidate HKLF parts of dervRefine & HessRefine? into new routine dervHKLF

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