source: trunk/GSASIIstrMath.py @ 1030

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

mods to CIF export
add residuals to histogram GUI
save histogram residual info in data[0] dictionary

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