source: trunk/GSASIIstrMath.py @ 1106

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

modify reflection array to be a dictionary with 4 items: RefList?, Uniq, Phi & FF. Each is list. Use patches to convert old format to new in various places.
Fix bug in simulate PWDR code.

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