source: trunk/GSASIIstrMath.py @ 989

Last change on this file since 989 was 989, checked in by toby, 9 years ago

add constraint derivs for single xtal; allow cellFill to work without esds; update docs; fix controls copy bug

File size: 98.3 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[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1239                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1240                Histogram[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,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,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,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,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    'Needs a doc string'
1588    parmDict.update(zip(varylist,values))
1589    G2mv.Dict2Map(parmDict,varylist)
1590    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1591    nvar = len(varylist)
1592    dMdv = np.empty(0)
1593    histoList = Histograms.keys()
1594    histoList.sort()
1595    for histogram in histoList:
1596        if 'PWDR' in histogram[:4]:
1597            Histogram = Histograms[histogram]
1598            hId = Histogram['hId']
1599            hfx = ':%d:'%(hId)
1600            wtFactor = calcControls[hfx+'wtFactor']
1601            Limits = calcControls[hfx+'Limits']
1602            x,y,w,yc,yb,yd = Histogram['Data']
1603            W = wtFactor*w
1604            xB = np.searchsorted(x,Limits[0])
1605            xF = np.searchsorted(x,Limits[1])
1606            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
1607                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1608        elif 'HKLF' in histogram[:4]:
1609            Histogram = Histograms[histogram]
1610            nobs = Histogram['Nobs']
1611            phase = Histogram['Reflection Lists']
1612            Phase = Phases[phase]
1613            hId = Histogram['hId']
1614            hfx = ':%d:'%(hId)
1615            wtFactor = calcControls[hfx+'wtFactor']
1616            pfx = '%d::'%(Phase['pId'])
1617            phfx = '%d:%d:'%(Phase['pId'],hId)
1618            SGData = Phase['General']['SGData']
1619            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1620            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1621            refList = Histogram['Data']
1622            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1623            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1624            dMdvh = np.zeros((len(varylist),len(refList)))
1625            dependentVars = G2mv.GetDependentVars()
1626            depDerivDict = {}
1627            for j in dependentVars:
1628                depDerivDict[j] = np.zeros(shape=(len(refList)))
1629            if calcControls['F**2']:
1630                for iref,ref in enumerate(refList):
1631                    if ref[6] > 0:
1632                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1633                        w = 1.0/ref[6]
1634                        if w*ref[5] >= calcControls['minF/sig']:
1635                            for j,var in enumerate(varylist):
1636                                if var in dFdvDict:
1637                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1638                            for var in dependentVars:
1639                                if var in dFdvDict:
1640                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*dervCor
1641                            if phfx+'Scale' in varylist:
1642                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1643                            elif phfx+'Scale' in dependentVars:
1644                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
1645                            for item in ['Ep','Es','Eg']:
1646                                if phfx+item in varylist:
1647                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1648                                elif phfx+item in dependentVars:
1649                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1650                            for item in ['BabA','BabU']:
1651                                if phfx+item in varylist:
1652                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1653                                elif phfx+item in dependentVars:
1654                                    depDerivDict[phfx+item][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1655            else:
1656                for iref,ref in enumerate(refList):
1657                    if ref[5] > 0.:
1658                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1659                        Fo = np.sqrt(ref[5])
1660                        Fc = np.sqrt(ref[7])
1661                        w = 1.0/ref[6]
1662                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1663                            for j,var in enumerate(varylist):
1664                                if var in dFdvDict:
1665                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1666                            for var in dependentVars:
1667                                if var in dFdvDict:
1668                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1669                            if phfx+'Scale' in varylist:
1670                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1671                            elif phfx+'Scale' in dependentVars:
1672                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor                           
1673                            for item in ['Ep','Es','Eg']:
1674                                if phfx+item in varylist:
1675                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1676                                elif phfx+item in dependentVars:
1677                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1678                            for item in ['BabA','BabU']:
1679                                if phfx+item in varylist:
1680                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']
1681                                elif phfx+item in dependentVars:
1682                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1683            # now process derivatives in constraints
1684            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1685        else:
1686            continue        #skip non-histogram entries
1687        if len(dMdv):
1688            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
1689        else:
1690            dMdv = np.sqrt(wtFactor)*dMdvh
1691           
1692    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1693    if np.any(pVals):
1694        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1695        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
1696       
1697    return dMdv
1698
1699def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
1700    'Needs a doc string'
1701    parmDict.update(zip(varylist,values))
1702    G2mv.Dict2Map(parmDict,varylist)
1703    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1704    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
1705    nvar = len(varylist)
1706    Hess = np.empty(0)
1707    histoList = Histograms.keys()
1708    histoList.sort()
1709    for histogram in histoList:
1710        if 'PWDR' in histogram[:4]:
1711            Histogram = Histograms[histogram]
1712            hId = Histogram['hId']
1713            hfx = ':%d:'%(hId)
1714            wtFactor = calcControls[hfx+'wtFactor']
1715            Limits = calcControls[hfx+'Limits']
1716            x,y,w,yc,yb,yd = Histogram['Data']
1717            W = wtFactor*w
1718            dy = y-yc
1719            xB = np.searchsorted(x,Limits[0])
1720            xF = np.searchsorted(x,Limits[1])
1721            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
1722                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
1723            Wt = np.sqrt(W[xB:xF])[np.newaxis,:]
1724            Dy = dy[xB:xF][np.newaxis,:]
1725            dMdvh *= Wt
1726            if dlg:
1727                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1728            if len(Hess):
1729                Hess += np.inner(dMdvh,dMdvh)
1730                dMdvh *= Wt*Dy
1731                Vec += np.sum(dMdvh,axis=1)
1732            else:
1733                Hess = np.inner(dMdvh,dMdvh)
1734                dMdvh *= Wt*Dy
1735                Vec = np.sum(dMdvh,axis=1)
1736        elif 'HKLF' in histogram[:4]:
1737            Histogram = Histograms[histogram]
1738            nobs = Histogram['Nobs']
1739            phase = Histogram['Reflection Lists']
1740            Phase = Phases[phase]
1741            hId = Histogram['hId']
1742            hfx = ':%d:'%(hId)
1743            wtFactor = calcControls[hfx+'wtFactor']
1744            pfx = '%d::'%(Phase['pId'])
1745            phfx = '%d:%d:'%(Phase['pId'],hId)
1746            SGData = Phase['General']['SGData']
1747            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1748            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1749            refList = Histogram['Data']
1750            time0 = time.time()
1751            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)  #accurate for powders!
1752            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
1753            dMdvh = np.zeros((len(varylist),len(refList)))
1754            dependentVars = G2mv.GetDependentVars()
1755            depDerivDict = {}
1756            for j in dependentVars:
1757                depDerivDict[j] = np.zeros(shape=(len(refList)))
1758            wdf = np.zeros(len(refList))
1759            if calcControls['F**2']:
1760                for iref,ref in enumerate(refList):
1761                    if ref[6] > 0:
1762                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1763                        w =  1.0/ref[6]
1764                        if w*ref[5] >= calcControls['minF/sig']:
1765                            wdf[iref] = w*(ref[5]-ref[7])
1766                            for j,var in enumerate(varylist):
1767                                if var in dFdvDict:
1768                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1769                            for var in dependentVars:
1770                                if var in dFdvDict:
1771                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1772                            if phfx+'Scale' in varylist:
1773                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1774                            elif phfx+'Scale' in dependentVars:
1775                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor
1776                            for item in ['Ep','Es','Eg']:
1777                                if phfx+item in varylist:
1778                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1779                                elif phfx+item in dependentVars:
1780                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1781                            for item in ['BabA','BabU']:
1782                                if phfx+item in varylist:
1783                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1784                                elif phfx+item in dependentVars:
1785                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1786            else:
1787                for iref,ref in enumerate(refList):
1788                    if ref[5] > 0.:
1789                        dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1790                        Fo = np.sqrt(ref[5])
1791                        Fc = np.sqrt(ref[7])
1792                        w = 1.0/ref[6]
1793                        if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
1794                            wdf[iref] = 2.0*Fo*w*(Fo-Fc)
1795                            for j,var in enumerate(varylist):
1796                                if var in dFdvDict:
1797                                    dMdvh[j][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1798                            for var in dependentVars:
1799                                if var in dFdvDict:
1800                                    depDerivDict[var][iref] = w*dFdvDict[var][iref]*dervCor*parmDict[phfx+'Scale']
1801                            if phfx+'Scale' in varylist:
1802                                dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*dervCor
1803                            elif phfx+'Scale' in dependentVars:
1804                                depDerivDict[phfx+'Scale'][iref] = w*ref[9]*dervCor                           
1805                            for item in ['Ep','Es','Eg']:
1806                                if phfx+item in varylist:
1807                                    dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1808                                elif phfx+item in dependentVars:
1809                                    depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]*parmDict[phfx+'Scale']
1810                            for item in ['BabA','BabU']:
1811                                if phfx+item in varylist:
1812                                    dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1813                                elif phfx+item in dependentVars:
1814                                    depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor
1815            # now process derivatives in constraints
1816            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
1817
1818            if dlg:
1819                dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1820            if len(Hess):
1821                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
1822                Hess += wtFactor*np.inner(dMdvh,dMdvh)
1823            else:
1824                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
1825                Hess = wtFactor*np.inner(dMdvh,dMdvh)
1826        else:
1827            continue        #skip non-histogram entries
1828    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1829    if np.any(pVals):
1830        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
1831        Vec += np.sum(dpdv*pWt*pVals,axis=1)
1832        Hess += np.inner(dpdv*pWt,dpdv)
1833    return Vec,Hess
1834
1835def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
1836    'Needs a doc string'
1837    parmDict.update(zip(varylist,values))
1838    Values2Dict(parmDict, varylist, values)
1839    G2mv.Dict2Map(parmDict,varylist)
1840    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
1841    M = np.empty(0)
1842    SumwYo = 0
1843    Nobs = 0
1844    ApplyRBModels(parmDict,Phases,rigidbodyDict)
1845    histoList = Histograms.keys()
1846    histoList.sort()
1847    for histogram in histoList:
1848        if 'PWDR' in histogram[:4]:
1849            Histogram = Histograms[histogram]
1850            hId = Histogram['hId']
1851            hfx = ':%d:'%(hId)
1852            wtFactor = calcControls[hfx+'wtFactor']
1853            Limits = calcControls[hfx+'Limits']
1854            x,y,w,yc,yb,yd = Histogram['Data']
1855            W = wtFactor*w
1856            yc *= 0.0                           #zero full calcd profiles
1857            yb *= 0.0
1858            yd *= 0.0
1859            xB = np.searchsorted(x,Limits[0])
1860            xF = np.searchsorted(x,Limits[1])
1861            Histogram['Nobs'] = xF-xB
1862            Nobs += Histogram['Nobs']
1863            Histogram['sumwYo'] = np.sum(W[xB:xF]*y[xB:xF]**2)
1864            SumwYo += Histogram['sumwYo']
1865            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
1866                varylist,Histogram,Phases,calcControls,pawleyLookup)
1867            yc[xB:xF] += yb[xB:xF]
1868            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
1869            Histogram['sumwYd'] = np.sum(np.sqrt(W[xB:xF])*(yd[xB:xF]))
1870            wdy = -np.sqrt(W[xB:xF])*(yd[xB:xF])
1871            Histogram['wR'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
1872            if dlg:
1873                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1874            M = np.concatenate((M,wdy))
1875#end of PWDR processing
1876        elif 'HKLF' in histogram[:4]:
1877            Histogram = Histograms[histogram]
1878            phase = Histogram['Reflection Lists']
1879            Phase = Phases[phase]
1880            hId = Histogram['hId']
1881            hfx = ':%d:'%(hId)
1882            wtFactor = calcControls[hfx+'wtFactor']
1883            pfx = '%d::'%(Phase['pId'])
1884            phfx = '%d:%d:'%(Phase['pId'],hId)
1885            SGData = Phase['General']['SGData']
1886            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1887            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1888            refList = Histogram['Data']
1889            time0 = time.time()
1890            StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1891#            print 'sf-calc time: %.3f'%(time.time()-time0)
1892            df = np.zeros(len(refList))
1893            sumwYo = 0
1894            sumFo = 0
1895            sumFo2 = 0
1896            sumdF = 0
1897            sumdF2 = 0
1898            nobs = 0
1899            if calcControls['F**2']:
1900                for i,ref in enumerate(refList):
1901                    if ref[6] > 0:
1902                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1903                        w = 1.0/ref[6]
1904                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1905                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1906                        ref[8] = ref[5]/parmDict[phfx+'Scale']
1907                        if w*ref[5] >= calcControls['minF/sig']:
1908                            sumFo2 += ref[5]
1909                            Fo = np.sqrt(ref[5])
1910                            sumFo += Fo
1911                            sumFo2 += ref[5]
1912                            sumdF += abs(Fo-np.sqrt(ref[7]))
1913                            sumdF2 += abs(ref[5]-ref[7])
1914                            nobs += 1
1915                            df[i] = -w*(ref[5]-ref[7])
1916                            sumwYo += (w*ref[5])**2
1917            else:
1918                for i,ref in enumerate(refList):
1919                    if ref[5] > 0.:
1920                        SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
1921                        ref[7] = parmDict[phfx+'Scale']*ref[9]
1922                        ref[7] *= ref[13]                       #correct Fc^2 for extinction
1923                        Fo = np.sqrt(ref[5])
1924                        Fc = np.sqrt(ref[7])
1925                        w = 2.0*Fo/ref[6]
1926                        if w*Fo >= calcControls['minF/sig']:
1927                            sumFo += Fo
1928                            sumFo2 += ref[5]
1929                            sumdF += abs(Fo-Fc)
1930                            sumdF2 += abs(ref[5]-ref[7])
1931                            nobs += 1
1932                            df[i] = -w*(Fo-Fc)
1933                            sumwYo += (w*Fo)**2
1934            Histogram['Nobs'] = nobs
1935            Histogram['sumwYo'] = sumwYo
1936            SumwYo += sumwYo
1937            Histogram['wR'] = min(100.,np.sqrt(np.sum(df**2)/Histogram['sumwYo'])*100.)
1938            Histogram[phfx+'Rf'] = 100.*sumdF/sumFo
1939            Histogram[phfx+'Rf^2'] = 100.*sumdF2/sumFo2
1940            Histogram[phfx+'Nref'] = nobs
1941            Nobs += nobs
1942            if dlg:
1943                dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
1944            M = np.concatenate((M,wtFactor*df))
1945# end of HKLF processing
1946    Histograms['sumwYo'] = SumwYo
1947    Histograms['Nobs'] = Nobs
1948    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
1949    if dlg:
1950        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
1951        if not GoOn:
1952            parmDict['saved values'] = values
1953            dlg.Destroy()
1954            raise Exception         #Abort!!
1955    pDict,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
1956    if np.any(pVals):
1957        pSum = np.sum(pWt*pVals**2)
1958        print 'Penalty function: %.3f on %d terms'%(pSum,len(pVals))
1959        Nobs += len(pVals)
1960        M = np.concatenate((M,np.sqrt(pWt)*pVals))
1961    return M
1962                       
Note: See TracBrowser for help on using the repository browser.