source: trunk/GSASIIstrMath.py @ 1244

Last change on this file since 1244 was 1244, checked in by vondreele, 9 years ago

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