source: trunk/GSASIIstrMath.py @ 1038

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

more work on residuals & cif files

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