source: trunk/GSASIIstrMath.py @ 1087

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

simulated CW powder data now complete with Poisson errors

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