source: trunk/GSASIIstrMath.py @ 1036

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

fix a couple of mistakes & get SC to cif OK

File size: 99.9 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['Residuals']['hId'] = Histograms[histogram]['hId']
1244               
1245def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1246    'Needs a doc string'
1247   
1248    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1249        U = parmDict[hfx+'U']
1250        V = parmDict[hfx+'V']
1251        W = parmDict[hfx+'W']
1252        X = parmDict[hfx+'X']
1253        Y = parmDict[hfx+'Y']
1254        tanPos = tand(refl[5]/2.0)
1255        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
1256        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
1257        sig = max(0.001,sig)
1258        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
1259        gam = max(0.001,gam)
1260        return sig,gam
1261               
1262    hId = Histogram['hId']
1263    hfx = ':%d:'%(hId)
1264    bakType = calcControls[hfx+'bakType']
1265    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1266    yc = np.zeros_like(yb)
1267       
1268    if 'C' in calcControls[hfx+'histType']:   
1269        shl = max(parmDict[hfx+'SH/L'],0.002)
1270        Ka2 = False
1271        if hfx+'Lam1' in parmDict.keys():
1272            wave = parmDict[hfx+'Lam1']
1273            Ka2 = True
1274            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1275            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1276        else:
1277            wave = parmDict[hfx+'Lam']
1278    else:
1279        print 'TOF Undefined at present'
1280        raise ValueError
1281    for phase in Histogram['Reflection Lists']:
1282        refList = Histogram['Reflection Lists'][phase]
1283        Phase = Phases[phase]
1284        pId = Phase['pId']
1285        pfx = '%d::'%(pId)
1286        phfx = '%d:%d:'%(pId,hId)
1287        hfx = ':%d:'%(hId)
1288        SGData = Phase['General']['SGData']
1289        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1290        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1291        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1292        Vst = np.sqrt(nl.det(G))    #V*
1293        if not Phase['General'].get('doPawley'):
1294            time0 = time.time()
1295            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1296#            print 'sf calc time: %.3fs'%(time.time()-time0)
1297        time0 = time.time()
1298        for refl in refList:
1299            if 'C' in calcControls[hfx+'histType']:
1300                h,k,l = refl[:3]
1301                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1302                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1303                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1304                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
1305                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
1306                refl[13] *= Vst*Lorenz
1307                if Phase['General'].get('doPawley'):
1308                    try:
1309                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1310                        refl[9] = parmDict[pInd]
1311                    except KeyError:
1312#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1313                        continue
1314                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1315                iBeg = np.searchsorted(x,refl[5]-fmin)
1316                iFin = np.searchsorted(x,refl[5]+fmax)
1317                if not iBeg+iFin:       #peak below low limit - skip peak
1318                    continue
1319                elif not iBeg-iFin:     #peak above high limit - done
1320                    break
1321                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
1322                if Ka2:
1323                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1324                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
1325                    iBeg = np.searchsorted(x,pos2-fmin)
1326                    iFin = np.searchsorted(x,pos2+fmax)
1327                    if not iBeg+iFin:       #peak below low limit - skip peak
1328                        continue
1329                    elif not iBeg-iFin:     #peak above high limit - done
1330                        return yc,yb
1331                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))        #and here
1332            elif 'T' in calcControls[hfx+'histType']:
1333                print 'TOF Undefined at present'
1334                raise Exception    #no TOF yet
1335#        print 'profile calc time: %.3fs'%(time.time()-time0)
1336    return yc,yb
1337   
1338def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
1339    'Needs a doc string'
1340   
1341    def cellVaryDerv(pfx,SGData,dpdA): 
1342        if SGData['SGLaue'] in ['-1',]:
1343            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1344                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1345        elif SGData['SGLaue'] in ['2/m',]:
1346            if SGData['SGUniq'] == 'a':
1347                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1348            elif SGData['SGUniq'] == 'b':
1349                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1350            else:
1351                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1352        elif SGData['SGLaue'] in ['mmm',]:
1353            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1354        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1355            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1356        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1357            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1358        elif SGData['SGLaue'] in ['3R', '3mR']:
1359            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1360        elif SGData['SGLaue'] in ['m3m','m3']:
1361            return [[pfx+'A0',dpdA[0]]]
1362           
1363    # create a list of dependent variables and set up a dictionary to hold their derivatives
1364    dependentVars = G2mv.GetDependentVars()
1365    depDerivDict = {}
1366    for j in dependentVars:
1367        depDerivDict[j] = np.zeros(shape=(len(x)))
1368    #print 'dependent vars',dependentVars
1369    lenX = len(x)               
1370    hId = Histogram['hId']
1371    hfx = ':%d:'%(hId)
1372    bakType = calcControls[hfx+'bakType']
1373    dMdv = np.zeros(shape=(len(varylist),len(x)))
1374    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1375    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
1376        bBpos =varylist.index(hfx+'Back:0')
1377        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1378    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
1379    for name in varylist:
1380        if 'Debye' in name:
1381            id = int(name.split(':')[-1])
1382            parm = name[:int(name.rindex(':'))]
1383            ip = names.index(parm)
1384            dMdv[varylist.index(name)] = dMddb[3*id+ip]
1385    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
1386    for name in varylist:
1387        if 'BkPk' in name:
1388            id = int(name.split(':')[-1])
1389            parm = name[:int(name.rindex(':'))]
1390            ip = names.index(parm)
1391            dMdv[varylist.index(name)] = dMdpk[4*id+ip]
1392    cw = np.diff(x)
1393    cw = np.append(cw,cw[-1])
1394    if 'C' in calcControls[hfx+'histType']:   
1395        shl = max(parmDict[hfx+'SH/L'],0.002)
1396        Ka2 = False
1397        if hfx+'Lam1' in parmDict.keys():
1398            wave = parmDict[hfx+'Lam1']
1399            Ka2 = True
1400            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1401            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1402        else:
1403            wave = parmDict[hfx+'Lam']
1404    else:
1405        print 'TOF Undefined at present'
1406        raise ValueError
1407    for phase in Histogram['Reflection Lists']:
1408        refList = Histogram['Reflection Lists'][phase]
1409        Phase = Phases[phase]
1410        SGData = Phase['General']['SGData']
1411        pId = Phase['pId']
1412        pfx = '%d::'%(pId)
1413        phfx = '%d:%d:'%(pId,hId)
1414        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1415        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1416        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1417        if not Phase['General'].get('doPawley'):
1418            time0 = time.time()
1419            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1420#            print 'sf-derv time %.3fs'%(time.time()-time0)
1421            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1422        time0 = time.time()
1423        for iref,refl in enumerate(refList):
1424            if 'C' in calcControls[hfx+'histType']:        #CW powder
1425                h,k,l = refl[:3]
1426                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
1427                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
1428                iBeg = np.searchsorted(x,refl[5]-fmin)
1429                iFin = np.searchsorted(x,refl[5]+fmax)
1430                if not iBeg+iFin:       #peak below low limit - skip peak
1431                    continue
1432                elif not iBeg-iFin:     #peak above high limit - done
1433                    break
1434                pos = refl[5]
1435                tanth = tand(pos/2.0)
1436                costh = cosd(pos/2.0)
1437                lenBF = iFin-iBeg
1438                dMdpk = np.zeros(shape=(6,lenBF))
1439                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
1440                for i in range(5):
1441                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[i]
1442                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
1443                if Ka2:
1444                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1445                    iBeg2 = np.searchsorted(x,pos2-fmin)
1446                    iFin2 = np.searchsorted(x,pos2+fmax)
1447                    if iBeg2-iFin2:
1448                        lenBF2 = iFin2-iBeg2
1449                        dMdpk2 = np.zeros(shape=(6,lenBF2))
1450                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
1451                        for i in range(5):
1452                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[i]
1453                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[13]*dMdipk2[0]
1454                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
1455                if Phase['General'].get('doPawley'):
1456                    dMdpw = np.zeros(len(x))
1457                    try:
1458                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
1459                        idx = varylist.index(pIdx)
1460                        dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
1461                        if Ka2:
1462                            dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
1463                        dMdv[idx] = dMdpw
1464                    except: # ValueError:
1465                        pass
1466                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1467                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1468                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1469                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1470                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1471                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1472                    hfx+'DisplaceY':[dpdY,'pos'],hfx+'Absorption':[dFdAb,'int'],}
1473                for name in names:
1474                    item = names[name]
1475                    if name in varylist:
1476                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
1477                        if Ka2:
1478                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1479                    elif name in dependentVars:
1480                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
1481                        if Ka2:
1482                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
1483                for iPO in dIdPO:
1484                    if iPO in varylist:
1485                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1486                        if Ka2:
1487                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1488                    elif iPO in dependentVars:
1489                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
1490                        if Ka2:
1491                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
1492                for i,name in enumerate(['omega','chi','phi']):
1493                    aname = pfx+'SH '+name
1494                    if aname in varylist:
1495                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
1496                        if Ka2:
1497                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1498                    elif aname in dependentVars:
1499                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
1500                        if Ka2:
1501                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
1502                for iSH in dFdODF:
1503                    if iSH in varylist:
1504                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1505                        if Ka2:
1506                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1507                    elif iSH in dependentVars:
1508                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
1509                        if Ka2:
1510                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
1511                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1512                for name,dpdA in cellDervNames:
1513                    if name in varylist:
1514                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
1515                        if Ka2:
1516                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
1517                    elif name in dependentVars:
1518                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
1519                        if Ka2:
1520                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
1521                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1522                for name in dDijDict:
1523                    if name in varylist:
1524                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1525                        if Ka2:
1526                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1527                    elif name in dependentVars:
1528                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
1529                        if Ka2:
1530                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
1531                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
1532                for name in gamDict:
1533                    if name in varylist:
1534                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
1535                        if Ka2:
1536                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1537                    elif name in dependentVars:
1538                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
1539                        if Ka2:
1540                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
1541                for name in sigDict:
1542                    if name in varylist:
1543                        dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
1544                        if Ka2:
1545                            dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1546                    elif name in dependentVars:
1547                        depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
1548                        if Ka2:
1549                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
1550                for name in ['BabA','BabU']:
1551                    if phfx+name in varylist:
1552                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
1553                        if Ka2:
1554                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]
1555                    elif phfx+name in dependentVars:                   
1556                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
1557                        if Ka2:
1558                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]                 
1559            elif 'T' in calcControls[hfx+'histType']:
1560                print 'TOF Undefined at present'
1561                raise Exception    #no TOF yet
1562            if not Phase['General'].get('doPawley'):
1563                #do atom derivatives -  for RB,F,X & U so far             
1564                corr = dervDict['int']/refl[9]
1565                if Ka2:
1566                    corr2 = dervDict2['int']/refl[9]
1567                for name in varylist+dependentVars:
1568                    if '::RBV;' in name:
1569                        pass
1570                    else:
1571                        try:
1572                            aname = name.split(pfx)[1][:2]
1573                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
1574                        except IndexError:
1575                            continue
1576                    if name in varylist:
1577                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
1578                        if Ka2:
1579                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1580                    elif name in dependentVars:
1581                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
1582                        if Ka2:
1583                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
1584#        print 'profile derv time: %.3fs'%(time.time()-time0)
1585    # now process derivatives in constraints
1586    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
1587    return dMdv
1588
1589def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1590    '''Loop over histograms and compute derivatives of the fitting
1591    model (M) with respect to all parameters.  Results are returned in
1592    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
1593    n is the number of parameters and m is the number of data
1594    points. This can exceed memory when m gets large. This routine is
1595    used when refinement derivatives are selected as "analtytic
1596    Jacobian" in Controls.
1597
1598    :returns: Jacobian numpy.array dMdv for all histograms concatinated
1599    '''
1600    parmDict.update(zip(varylist,values))
1601    G2mv.Dict2Map(parmDict,varylist)
1602    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1603    nvar = len(varylist)
1604    dMdv = np.empty(0)
1605    histoList = Histograms.keys()
1606    histoList.sort()
1607    for histogram in histoList:
1608        if 'PWDR' in histogram[:4]:
1609            Histogram = Histograms[histogram]
1610            hId = Histogram['hId']
1611            hfx = ':%d:'%(hId)
1612            wtFactor = calcControls[hfx+'wtFactor']
1613            Limits = calcControls[hfx+'Limits']
1614            x,y,w,yc,yb,yd = Histogram['Data']
1615            W = wtFactor*w
1616            xB = np.searchsorted(x,Limits[0])
1617            xF = np.searchsorted(x,Limits[1])
1618            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
1619                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1620        elif 'HKLF' in histogram[:4]:
1621            Histogram = Histograms[histogram]
1622            nobs = Histogram['Residuals']['Nobs']
1623            phase = Histogram['Reflection Lists']
1624            Phase = Phases[phase]
1625            hId = Histogram['hId']
1626            hfx = ':%d:'%(hId)
1627            wtFactor = calcControls[hfx+'wtFactor']
1628            pfx = '%d::'%(Phase['pId'])
1629            phfx = '%d:%d:'%(Phase['pId'],hId)
1630            SGData = Phase['General']['SGData']
1631            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1632            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1633            refList = Histogram['Data']
1634            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1635            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1636            dMdvh = np.zeros((len(varylist),len(refList)))
1637            dependentVars = G2mv.GetDependentVars()
1638            depDerivDict = {}
1639            for j in dependentVars:
1640                depDerivDict[j] = np.zeros(shape=(len(refList)))
1641            if calcControls['F**2']:
1642                for iref,ref in enumerate(refList):
1643                    if ref[6] > 0:
1644                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1645                        w = 1.0/ref[6]
1646                        if w*ref[5] >= calcControls['minF/sig']:
1647                            for j,var in enumerate(varylist):
1648                                if var in dFdvDict:
1649                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1650                            for var in dependentVars:
1651                                if var in dFdvDict:
1652                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1653                            if phfx+'Scale' in varylist:
1654                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1655                            elif phfx+'Scale' in dependentVars:
1656                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
1657                            for item in ['Ep','Es','Eg']:
1658                                if phfx+item in varylist:
1659                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1660                                elif phfx+item in dependentVars:
1661                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1662                            for item in ['BabA','BabU']:
1663                                if phfx+item in varylist:
1664                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1665                                elif phfx+item in dependentVars:
1666                                    depDerivDict[phfx+item][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1667            else:
1668                for iref,ref in enumerate(refList):
1669                    if ref[5] > 0.:
1670                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1671                        Fo = np.sqrt(ref[5])
1672                        Fc = np.sqrt(ref[7])
1673                        w = 1.0/ref[6]
1674                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1675                            for j,var in enumerate(varylist):
1676                                if var in dFdvDict:
1677                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1678                            for var in dependentVars:
1679                                if var in dFdvDict:
1680                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1681                            if phfx+'Scale' in varylist:
1682                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1683                            elif phfx+'Scale' in dependentVars:
1684                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor                           
1685                            for item in ['Ep','Es','Eg']:
1686                                if phfx+item in varylist:
1687                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1688                                elif phfx+item in dependentVars:
1689                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1690                            for item in ['BabA','BabU']:
1691                                if phfx+item in varylist:
1692                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1693                                elif phfx+item in dependentVars:
1694                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1695            # now process derivatives in constraints
1696            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1697        else:
1698            continue        #skip non-histogram entries
1699        if len(dMdv):
1700            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
1701        else:
1702            dMdv = np.sqrt(wtFactor)*dMdvh
1703           
1704    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1705    if np.any(pVals):
1706        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1707        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
1708       
1709    return dMdv
1710
1711def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1712    '''Loop over histograms and compute derivatives of the fitting
1713    model (M) with respect to all parameters.  For each histogram, the
1714    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
1715    number of parameters and m is the number of data points *in the
1716    histogram*. The (n by n) Hessian is computed from each Jacobian
1717    and it is returned.  This routine is used when refinement
1718    derivatives are selected as "analtytic Hessian" in Controls.
1719
1720    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
1721    '''
1722    'Needs a doc string'
1723    parmDict.update(zip(varylist,values))
1724    G2mv.Dict2Map(parmDict,varylist)
1725    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1726    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
1727    nvar = len(varylist)
1728    Hess = np.empty(0)
1729    histoList = Histograms.keys()
1730    histoList.sort()
1731    for histogram in histoList:
1732        if 'PWDR' in histogram[:4]:
1733            Histogram = Histograms[histogram]
1734            hId = Histogram['hId']
1735            hfx = ':%d:'%(hId)
1736            wtFactor = calcControls[hfx+'wtFactor']
1737            Limits = calcControls[hfx+'Limits']
1738            x,y,w,yc,yb,yd = Histogram['Data']
1739            W = wtFactor*w
1740            dy = y-yc
1741            xB = np.searchsorted(x,Limits[0])
1742            xF = np.searchsorted(x,Limits[1])
1743            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
1744                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1745            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
1746            Dy = dy[xB:xF][np.newaxis,:]
1747            dMdvh *= Wt
1748            if dlg:
1749                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1750            if len(Hess):
1751                Hess += np.inner(dMdvh,dMdvh)
1752                dMdvh *= Wt*Dy
1753                Vec += np.sum(dMdvh,axis=1)
1754            else:
1755                Hess = np.inner(dMdvh,dMdvh)
1756                dMdvh *= Wt*Dy
1757                Vec = np.sum(dMdvh,axis=1)
1758        elif 'HKLF' in histogram[:4]:
1759            Histogram = Histograms[histogram]
1760            nobs = Histogram['Residuals']['Nobs']
1761            phase = Histogram['Reflection Lists']
1762            Phase = Phases[phase]
1763            hId = Histogram['hId']
1764            hfx = ':%d:'%(hId)
1765            wtFactor = calcControls[hfx+'wtFactor']
1766            pfx = '%d::'%(Phase['pId'])
1767            phfx = '%d:%d:'%(Phase['pId'],hId)
1768            SGData = Phase['General']['SGData']
1769            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1770            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1771            refList = Histogram['Data']
1772            time0 = time.time()
1773            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)  #accurate for powders!
1774            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1775            dMdvh = np.zeros((len(varylist),len(refList)))
1776            dependentVars = G2mv.GetDependentVars()
1777            depDerivDict = {}
1778            for j in dependentVars:
1779                depDerivDict[j] = np.zeros(shape=(len(refList)))
1780            wdf = np.zeros(len(refList))
1781            if calcControls['F**2']:
1782                for iref,ref in enumerate(refList):
1783                    if ref[6] > 0:
1784                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1785                        w =  1.0/ref[6]
1786                        if w*ref[5] >= calcControls['minF/sig']:
1787                            wdf[iref] = w*(ref[5]-ref[7])
1788                            for j,var in enumerate(varylist):
1789                                if var in dFdvDict:
1790                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1791                            for var in dependentVars:
1792                                if var in dFdvDict:
1793                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1794                            if phfx+'Scale' in varylist:
1795                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1796                            elif phfx+'Scale' in dependentVars:
1797                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
1798                            for item in ['Ep','Es','Eg']:
1799                                if phfx+item in varylist:
1800                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1801                                elif phfx+item in dependentVars:
1802                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1803                            for item in ['BabA','BabU']:
1804                                if phfx+item in varylist:
1805                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1806                                elif phfx+item in dependentVars:
1807                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1808            else:
1809                for iref,ref in enumerate(refList):
1810                    if ref[5] > 0.:
1811                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1812                        Fo = np.sqrt(ref[5])
1813                        Fc = np.sqrt(ref[7])
1814                        w = 1.0/ref[6]
1815                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1816                            wdf[iref] = 2.0*Fo*w*(Fo-Fc)
1817                            for j,var in enumerate(varylist):
1818                                if var in dFdvDict:
1819                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1820                            for var in dependentVars:
1821                                if var in dFdvDict:
1822                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1823                            if phfx+'Scale' in varylist:
1824                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1825                            elif phfx+'Scale' in dependentVars:
1826                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor                           
1827                            for item in ['Ep','Es','Eg']:
1828                                if phfx+item in varylist:
1829                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1830                                elif phfx+item in dependentVars:
1831                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1832                            for item in ['BabA','BabU']:
1833                                if phfx+item in varylist:
1834                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1835                                elif phfx+item in dependentVars:
1836                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1837            # now process derivatives in constraints
1838            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1839
1840            if dlg:
1841                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1842            if len(Hess):
1843                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
1844                Hess += wtFactor*np.inner(dMdvh,dMdvh)
1845            else:
1846                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
1847                Hess = wtFactor*np.inner(dMdvh,dMdvh)
1848        else:
1849            continue        #skip non-histogram entries
1850    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1851    if np.any(pVals):
1852        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1853        Vec += np.sum(dpdv*pWt*pVals,axis=1)
1854        Hess += np.inner(dpdv*pWt,dpdv)
1855    return Vec,Hess
1856
1857def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
1858    'Needs a doc string'
1859    parmDict.update(zip(varylist,values))
1860    Values2Dict(parmDict, varylist, values)
1861    G2mv.Dict2Map(parmDict,varylist)
1862    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1863    M = np.empty(0)
1864    SumwYo = 0
1865    Nobs = 0
1866    ApplyRBModels(parmDict,Phases,rigidbodyDict)
1867    histoList = Histograms.keys()
1868    histoList.sort()
1869    for histogram in histoList:
1870        if 'PWDR' in histogram[:4]:
1871            Histogram = Histograms[histogram]
1872            hId = Histogram['hId']
1873            hfx = ':%d:'%(hId)
1874            wtFactor = calcControls[hfx+'wtFactor']
1875            Limits = calcControls[hfx+'Limits']
1876            x,y,w,yc,yb,yd = Histogram['Data']
1877            W = wtFactor*w
1878            yc *= 0.0                           #zero full calcd profiles
1879            yb *= 0.0
1880            yd *= 0.0
1881            xB = np.searchsorted(x,Limits[0])
1882            xF = np.searchsorted(x,Limits[1])
1883            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
1884            Nobs += Histogram['Residuals']['Nobs']
1885            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
1886            SumwYo += Histogram['Residuals']['sumwYo']
1887            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
1888                varylist,Histogram,Phases,calcControls,pawleyLookup)
1889            yc[xB:xF] += yb[xB:xF]
1890            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
1891            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
1892            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
1893            if dlg:
1894                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1895            M = np.concatenate((M,wdy))
1896#end of PWDR processing
1897        elif 'HKLF' in histogram[:4]:
1898            Histogram = Histograms[histogram]
1899            phase = Histogram['Reflection Lists']
1900            Phase = Phases[phase]
1901            hId = Histogram['hId']
1902            hfx = ':%d:'%(hId)
1903            wtFactor = calcControls[hfx+'wtFactor']
1904            pfx = '%d::'%(Phase['pId'])
1905            phfx = '%d:%d:'%(Phase['pId'],hId)
1906            SGData = Phase['General']['SGData']
1907            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1908            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1909            refList = Histogram['Data']
1910            time0 = time.time()
1911            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1912#            print 'sf-calc time: %.3f'%(time.time()-time0)
1913            df = np.zeros(len(refList))
1914            sumwYo = 0
1915            sumFo = 0
1916            sumFo2 = 0
1917            sumdF = 0
1918            sumdF2 = 0
1919            nobs = 0
1920            if calcControls['F**2']:
1921                for i,ref in enumerate(refList):
1922                    if ref[6] > 0:
1923                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1924                        w = 1.0/ref[6]
1925                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1926                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1927                        ref[8] = ref[5]/parmDict[phfx+'Scale']
1928                        if w*ref[5] >= calcControls['minF/sig']:
1929                            sumFo2 += ref[5]
1930                            Fo = np.sqrt(ref[5])
1931                            sumFo += Fo
1932                            sumFo2 += ref[5]
1933                            sumdF += abs(Fo-np.sqrt(ref[7]))
1934                            sumdF2 += abs(ref[5]-ref[7])
1935                            nobs += 1
1936                            df[i] = -w*(ref[5]-ref[7])
1937                            sumwYo += (w*ref[5])**2
1938            else:
1939                for i,ref in enumerate(refList):
1940                    if ref[5] > 0.:
1941                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1942                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1943                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1944                        Fo = np.sqrt(ref[5])
1945                        Fc = np.sqrt(ref[7])
1946                        w = 2.0*Fo/ref[6]
1947                        if w*Fo >= calcControls['minF/sig']:
1948                            sumFo += Fo
1949                            sumFo2 += ref[5]
1950                            sumdF += abs(Fo-Fc)
1951                            sumdF2 += abs(ref[5]-ref[7])
1952                            nobs += 1
1953                            df[i] = -w*(Fo-Fc)
1954                            sumwYo += (w*Fo)**2
1955            Histogram['Residuals']['Nobs'] = nobs
1956            Histogram['Residuals']['sumwYo'] = sumwYo
1957            SumwYo += sumwYo
1958            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['Residuals']['sumwYo'])*100.)
1959            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
1960            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
1961            Histogram['Residuals'][phfx+'Nref'] = nobs
1962            Nobs += nobs
1963            if dlg:
1964                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
1965            M = np.concatenate((M,wtFactor*df))
1966# end of HKLF processing
1967    Histograms['sumwYo'] = SumwYo
1968    Histograms['Nobs'] = Nobs
1969    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
1970    if dlg:
1971        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
1972        if not GoOn:
1973            parmDict['saved values'] = values
1974            dlg.Destroy()
1975            raise Exception         #Abort!!
1976    pDict,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1977    if np.any(pVals):
1978        pSum = np.sum(pWt*pVals**2)
1979        print 'Penalty function: %.3f on %d terms'%(pSum,len(pVals))
1980        Nobs += len(pVals)
1981        M = np.concatenate((M,np.sqrt(pWt)*pVals))
1982    return M
1983                       
Note: See TracBrowser for help on using the repository browser.