source: trunk/GSASIIstrMath.py @ 1992

Last change on this file since 1992 was 1992, checked in by vondreele, 6 years ago

fix SS position derivative error

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 155.1 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2015-10-07 20:35:11 +0000 (Wed, 07 Oct 2015) $
8# $Author: vondreele $
9# $Revision: 1992 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 1992 2015-10-07 20:35:11Z vondreele $
12########### SVN repository information ###################
13import time
14import copy
15import numpy as np
16import numpy.ma as ma
17import numpy.linalg as nl
18import scipy.optimize as so
19import scipy.stats as st
20import GSASIIpath
21GSASIIpath.SetVersionNumber("$Revision: 1992 $")
22import GSASIIElem as G2el
23import GSASIIlattice as G2lat
24import GSASIIspc as G2spc
25import GSASIIpwd as G2pwd
26import GSASIImapvars as G2mv
27import GSASIImath as G2mth
28import GSASIIobj as G2obj
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*np.log(2.0)
38twopi = 2.0*np.pi
39twopisq = 2.0*np.pi**2
40
41################################################################################
42##### Rigid Body Models
43################################################################################
44       
45def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
46    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
47    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
48    '''
49    atxIds = ['Ax:','Ay:','Az:']
50    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
51    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
52    if not RBIds['Vector'] and not RBIds['Residue']:
53        return
54    VRBIds = RBIds['Vector']
55    RRBIds = RBIds['Residue']
56    if Update:
57        RBData = rigidbodyDict
58    else:
59        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
60    if RBIds['Vector']:                       # first update the vector magnitudes
61        VRBData = RBData['Vector']
62        for i,rbId in enumerate(VRBIds):
63            if VRBData[rbId]['useCount']:
64                for j in range(len(VRBData[rbId]['VectMag'])):
65                    name = '::RBV;'+str(j)+':'+str(i)
66                    VRBData[rbId]['VectMag'][j] = parmDict[name]
67    for phase in Phases:
68        Phase = Phases[phase]
69        General = Phase['General']
70        cx,ct,cs,cia = General['AtomPtrs']
71        cell = General['Cell'][1:7]
72        Amat,Bmat = G2lat.cell2AB(cell)
73        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
74        pfx = str(Phase['pId'])+'::'
75        if Update:
76            RBModels = Phase['RBModels']
77        else:
78            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
79        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
80            jrb = VRBIds.index(RBObj['RBId'])
81            rbsx = str(irb)+':'+str(jrb)
82            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
83                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
84            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
85                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
86            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
87            TLS = RBObj['ThermalMotion']
88            if 'T' in TLS[0]:
89                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
90                    TLS[1][i] = parmDict[pfx+pt+rbsx]
91            if 'L' in TLS[0]:
92                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
93                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
94            if 'S' in TLS[0]:
95                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
96                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
97            if 'U' in TLS[0]:
98                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
99            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
100            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
101            for i,x in enumerate(XYZ):
102                atId = RBObj['Ids'][i]
103                for j in [0,1,2]:
104                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
105                if UIJ[i][0] == 'A':
106                    for j in range(6):
107                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
108                elif UIJ[i][0] == 'I':
109                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
110           
111        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
112            jrb = RRBIds.index(RBObj['RBId'])
113            rbsx = str(irb)+':'+str(jrb)
114            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
115                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
116            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
117                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
118            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
119            TLS = RBObj['ThermalMotion']
120            if 'T' in TLS[0]:
121                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
122                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
123            if 'L' in TLS[0]:
124                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
125                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
126            if 'S' in TLS[0]:
127                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
128                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
129            if 'U' in TLS[0]:
130                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
131            for itors,tors in enumerate(RBObj['Torsions']):
132                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
133            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
134            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
135            for i,x in enumerate(XYZ):
136                atId = RBObj['Ids'][i]
137                for j in [0,1,2]:
138                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
139                if UIJ[i][0] == 'A':
140                    for j in range(6):
141                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
142                elif UIJ[i][0] == 'I':
143                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
144                   
145def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
146    'Needs a doc string'
147    atxIds = ['dAx:','dAy:','dAz:']
148    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
149    TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
150    LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
151    SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
152    PIds = ['Px:','Py:','Pz:']
153    OIds = ['Oa:','Oi:','Oj:','Ok:']
154    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
155    if not RBIds['Vector'] and not RBIds['Residue']:
156        return
157    VRBIds = RBIds['Vector']
158    RRBIds = RBIds['Residue']
159    RBData = rigidbodyDict
160    for item in parmDict:
161        if 'RB' in item:
162            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
163    General = Phase['General']
164    cx,ct,cs,cia = General['AtomPtrs']
165    cell = General['Cell'][1:7]
166    Amat,Bmat = G2lat.cell2AB(cell)
167    rpd = np.pi/180.
168    rpd2 = rpd**2
169    g = nl.inv(np.inner(Bmat,Bmat))
170    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
171        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
172    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
173    pfx = str(Phase['pId'])+'::'
174    RBModels =  Phase['RBModels']
175    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
176        VModel = RBData['Vector'][RBObj['RBId']]
177        Q = RBObj['Orient'][0]
178        Pos = RBObj['Orig'][0]
179        jrb = VRBIds.index(RBObj['RBId'])
180        rbsx = str(irb)+':'+str(jrb)
181        dXdv = []
182        for iv in range(len(VModel['VectMag'])):
183            dCdv = []
184            for vec in VModel['rbVect'][iv]:
185                dCdv.append(G2mth.prodQVQ(Q,vec))
186            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
187        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
188        for ia,atId in enumerate(RBObj['Ids']):
189            atNum = AtLookup[atId]
190            dx = 0.00001
191            for iv in range(len(VModel['VectMag'])):
192                for ix in [0,1,2]:
193                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
194            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
195                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
196            for iv in range(4):
197                Q[iv] -= dx
198                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
199                Q[iv] += 2.*dx
200                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
201                Q[iv] -= dx
202                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
203                for ix in [0,1,2]:
204                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
205            X = G2mth.prodQVQ(Q,Cart[ia])
206            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
207            dFdu = G2lat.U6toUij(dFdu.T)
208            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
209            dFdu = G2lat.UijtoU6(dFdu)
210            atNum = AtLookup[atId]
211            if 'T' in RBObj['ThermalMotion'][0]:
212                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
213                    dFdvDict[pfx+name+rbsx] += dFdu[i]
214            if 'L' in RBObj['ThermalMotion'][0]:
215                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
216                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
217                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
218                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
219                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
220                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
221                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
222                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
223                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
224            if 'S' in RBObj['ThermalMotion'][0]:
225                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
226                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
227                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
228                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
229                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
230                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
231                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
232                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
233            if 'U' in RBObj['ThermalMotion'][0]:
234                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
235
236
237    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
238        Q = RBObj['Orient'][0]
239        Pos = RBObj['Orig'][0]
240        jrb = RRBIds.index(RBObj['RBId'])
241        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
242        rbsx = str(irb)+':'+str(jrb)
243        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
244        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
245            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
246            orId,pvId = torData[itors][:2]
247            pivotVec = Cart[orId]-Cart[pvId]
248            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
249            QB = G2mth.AVdeg2Q(0.001,pivotVec)
250            for ir in torData[itors][3]:
251                atNum = AtLookup[RBObj['Ids'][ir]]
252                rVec = Cart[ir]-Cart[pvId]
253                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
254                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
255                for ix in [0,1,2]:
256                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
257        for ia,atId in enumerate(RBObj['Ids']):
258            atNum = AtLookup[atId]
259            dx = 0.00001
260            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
261                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
262            for iv in range(4):
263                Q[iv] -= dx
264                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
265                Q[iv] += 2.*dx
266                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
267                Q[iv] -= dx
268                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
269                for ix in [0,1,2]:
270                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
271            X = G2mth.prodQVQ(Q,Cart[ia])
272            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
273            dFdu = G2lat.U6toUij(dFdu.T)
274            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
275            dFdu = G2lat.UijtoU6(dFdu)
276            atNum = AtLookup[atId]
277            if 'T' in RBObj['ThermalMotion'][0]:
278                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
279                    dFdvDict[pfx+name+rbsx] += dFdu[i]
280            if 'L' in RBObj['ThermalMotion'][0]:
281                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
282                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
283                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
284                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
285                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
286                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
287                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
288                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
289                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
290            if 'S' in RBObj['ThermalMotion'][0]:
291                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
292                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
293                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
294                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
295                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
296                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
297                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
298                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
299            if 'U' in RBObj['ThermalMotion'][0]:
300                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
301   
302################################################################################
303##### Penalty & restraint functions
304################################################################################
305
306def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
307    'Needs a doc string'
308    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
309    pNames = []
310    pVals = []
311    pWt = []
312    negWt = {}
313    pWsum = {}
314    for phase in Phases:
315        pId = Phases[phase]['pId']
316        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
317        General = Phases[phase]['General']
318        cx,ct,cs,cia = General['AtomPtrs']
319        textureData = General['SH Texture']
320        SGData = General['SGData']
321        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
322        cell = General['Cell'][1:7]
323        Amat,Bmat = G2lat.cell2AB(cell)
324        if phase not in restraintDict:
325            continue
326        phaseRest = restraintDict[phase]
327        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
328            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
329            ['ChemComp','Sites'],['Texture','HKLs'],]
330        for name,rest in names:
331            pWsum[name] = 0.
332            itemRest = phaseRest[name]
333            if itemRest[rest] and itemRest['Use']:
334                wt = itemRest['wtFactor']
335                if name in ['Bond','Angle','Plane','Chiral']:
336                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
337                        pNames.append(str(pId)+':'+name+':'+str(i))
338                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
339                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
340                        if name == 'Bond':
341                            calc = G2mth.getRestDist(XYZ,Amat)
342                        elif name == 'Angle':
343                            calc = G2mth.getRestAngle(XYZ,Amat)
344                        elif name == 'Plane':
345                            calc = G2mth.getRestPlane(XYZ,Amat)
346                        elif name == 'Chiral':
347                            calc = G2mth.getRestChiral(XYZ,Amat)
348                        pVals.append(obs-calc)
349                        pWt.append(wt/esd**2)
350                        pWsum[name] += wt*((obs-calc)/esd)**2
351                elif name in ['Torsion','Rama']:
352                    coeffDict = itemRest['Coeff']
353                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
354                        pNames.append(str(pId)+':'+name+':'+str(i))
355                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
356                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
357                        if name == 'Torsion':
358                            tor = G2mth.getRestTorsion(XYZ,Amat)
359                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
360                        else:
361                            phi,psi = G2mth.getRestRama(XYZ,Amat)
362                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
363                        pVals.append(restr)
364                        pWt.append(wt/esd**2)
365                        pWsum[name] += wt*(restr/esd)**2
366                elif name == 'ChemComp':
367                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
368                        pNames.append(str(pId)+':'+name+':'+str(i))
369                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
370                        frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
371                        calc = np.sum(mul*frac*factors)
372                        pVals.append(obs-calc)
373                        pWt.append(wt/esd**2)                   
374                        pWsum[name] += wt*((obs-calc)/esd)**2
375                elif name == 'Texture':
376                    SHkeys = textureData['SH Coeff'][1].keys()
377                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
378                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
379                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
380                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
381                        PH = np.array(hkl)
382                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
383                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
384                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
385                        Z1 = -ma.masked_greater(Z,0.0)
386                        IndZ1 = np.array(ma.nonzero(Z1))
387                        for ind in IndZ1.T:
388                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
389                            pVals.append(Z1[ind[0]][ind[1]])
390                            pWt.append(wt/esd1**2)
391                            pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd1)**2
392                        if ifesd2:
393                            Z2 = 1.-Z
394                            for ind in np.ndindex(grid,grid):
395                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
396                                pVals.append(Z1[ind[0]][ind[1]])
397                                pWt.append(wt/esd2**2)
398                                pWsum[name] += wt*(Z2/esd2)**2
399       
400    for phase in Phases:
401        name = 'SH-Pref.Ori.'
402        pId = Phases[phase]['pId']
403        General = Phases[phase]['General']
404        SGData = General['SGData']
405        cell = General['Cell'][1:7]
406        pWsum[name] = 0.0
407        for hist in Phases[phase]['Histograms']:
408            if hist in Histograms and 'PWDR' in hist:
409                hId = Histograms[hist]['hId']
410                phfx = '%d:%d:'%(pId,hId)
411                if calcControls[phfx+'poType'] == 'SH':
412                    toler = calcControls[phfx+'SHtoler']
413                    wt = 1./toler**2
414                    HKLs = np.array(calcControls[phfx+'SHhkl'])
415                    SHnames = calcControls[phfx+'SHnames']
416                    SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
417                    for i,PH in enumerate(HKLs):
418                        phi,beta = G2lat.CrsAng(PH,cell,SGData)
419                        SH3Coef = {}
420                        for item in SHcof:
421                            L,N = eval(item.strip('C'))
422                            SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
423                        ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
424                        X = np.linspace(0,90.0,26)
425                        Y = -ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)
426                        IndY = ma.nonzero(Y)
427                        for ind in IndY[0]:
428                            pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
429                            pVals.append(Y[ind])
430                            pWt.append(wt)
431                            pWsum[name] += wt*(Y[ind])**2
432    pWsum['PWLref'] = 0.
433    for item in varyList:
434        if 'PWLref' in item and parmDict[item] < 0.:
435            pId = int(item.split(':')[0])
436            if negWt[pId]:
437                pNames.append(item)
438                pVals.append(-parmDict[item])
439                pWt.append(negWt[pId])
440                pWsum['PWLref'] += negWt[pId]*(-parmDict[item])**2
441    pVals = np.array(pVals)
442    pWt = np.array(pWt)         #should this be np.sqrt?
443    return pNames,pVals,pWt,pWsum
444   
445def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
446    'Needs a doc string'
447    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
448    pDerv = np.zeros((len(varyList),len(pVal)))
449    for phase in Phases:
450#        if phase not in restraintDict:
451#            continue
452        pId = Phases[phase]['pId']
453        General = Phases[phase]['General']
454        cx,ct,cs,cia = General['AtomPtrs']
455        SGData = General['SGData']
456        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
457        cell = General['Cell'][1:7]
458        Amat,Bmat = G2lat.cell2AB(cell)
459        textureData = General['SH Texture']
460
461        SHkeys = textureData['SH Coeff'][1].keys()
462        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
463        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
464        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
465        sam = SamSym[textureData['Model']]
466        phaseRest = restraintDict.get(phase,{})
467        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
468            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
469            'ChemComp':'Sites','Texture':'HKLs'}
470        lasthkl = np.array([0,0,0])
471        for ip,pName in enumerate(pNames):
472            pnames = pName.split(':')
473            if pId == int(pnames[0]):
474                name = pnames[1]
475                if 'PWL' in pName:
476                    pDerv[varyList.index(pName)][ip] += 1.
477                    continue
478                elif 'SH-' in pName:
479                    continue
480                id = int(pnames[2]) 
481                itemRest = phaseRest[name]
482                if name in ['Bond','Angle','Plane','Chiral']:
483                    indx,ops,obs,esd = itemRest[names[name]][id]
484                    dNames = []
485                    for ind in indx:
486                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
487                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
488                    if name == 'Bond':
489                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
490                    elif name == 'Angle':
491                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
492                    elif name == 'Plane':
493                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
494                    elif name == 'Chiral':
495                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
496                elif name in ['Torsion','Rama']:
497                    coffDict = itemRest['Coeff']
498                    indx,ops,cofName,esd = itemRest[names[name]][id]
499                    dNames = []
500                    for ind in indx:
501                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
502                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
503                    if name == 'Torsion':
504                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
505                    else:
506                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
507                elif name == 'ChemComp':
508                    indx,factors,obs,esd = itemRest[names[name]][id]
509                    dNames = []
510                    for ind in indx:
511                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
512                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
513                        deriv = mul*factors
514                elif 'Texture' in name:
515                    deriv = []
516                    dNames = []
517                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
518                    hkl = np.array(hkl)
519                    if np.any(lasthkl-hkl):
520                        PH = np.array(hkl)
521                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
522                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
523                        lasthkl = copy.copy(hkl)                       
524                    if 'unit' in name:
525                        pass
526                    else:
527                        gam = float(pnames[3])
528                        psi = float(pnames[4])
529                        for SHname in ODFln:
530                            l,m,n = eval(SHname[1:])
531                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
532                            dNames += [str(pId)+'::'+SHname]
533                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
534                for dName,drv in zip(dNames,deriv):
535                    try:
536                        ind = varyList.index(dName)
537                        pDerv[ind][ip] += drv
538                    except ValueError:
539                        pass
540       
541        lasthkl = np.array([0,0,0])
542        for ip,pName in enumerate(pNames):
543            deriv = []
544            dNames = []
545            pnames = pName.split(':')
546            if 'SH-' in pName and pId == int(pnames[0]):
547                hId = int(pnames[1])
548                phfx = '%d:%d:'%(pId,hId)
549                psi = float(pnames[4])
550                HKLs = calcControls[phfx+'SHhkl']
551                SHnames = calcControls[phfx+'SHnames']
552                SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
553                hkl = np.array(HKLs[int(pnames[3])])     
554                if np.any(lasthkl-hkl):
555                    PH = np.array(hkl)
556                    phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
557                    SH3Coef = {}
558                    for item in SHcof:
559                        L,N = eval(item.strip('C'))
560                        SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
561                    ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
562                    lasthkl = copy.copy(hkl)                       
563                for SHname in SHnames:
564                    l,n = eval(SHname[1:])
565                    SH3name = 'C%d,0,%d'%(l,n)
566                    Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
567                    dNames += [phfx+SHname]
568                    deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
569            for dName,drv in zip(dNames,deriv):
570                try:
571                    ind = varyList.index(dName)
572                    pDerv[ind][ip] += drv
573                except ValueError:
574                    pass
575    return pDerv
576
577################################################################################
578##### Function & derivative calculations
579################################################################################       
580                   
581def GetAtomFXU(pfx,calcControls,parmDict):
582    'Needs a doc string'
583    Natoms = calcControls['Natoms'][pfx]
584    Tdata = Natoms*[' ',]
585    Mdata = np.zeros(Natoms)
586    IAdata = Natoms*[' ',]
587    Fdata = np.zeros(Natoms)
588    FFdata = []
589    BLdata = []
590    Xdata = np.zeros((3,Natoms))
591    dXdata = np.zeros((3,Natoms))
592    Uisodata = np.zeros(Natoms)
593    Uijdata = np.zeros((6,Natoms))
594    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
595        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
596        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
597        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
598        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
599    for iatm in range(Natoms):
600        for key in keys:
601            parm = pfx+key+str(iatm)
602            if parm in parmDict:
603                keys[key][iatm] = parmDict[parm]
604    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.
605    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
606   
607def GetAtomSSFXU(pfx,calcControls,parmDict):
608    'Needs a doc string'
609    Natoms = calcControls['Natoms'][pfx]
610    maxSSwave = calcControls['maxSSwave'][pfx]
611    Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
612        'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
613    XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
614    FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
615    USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
616    MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
617    waveTypes = []
618    keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
619        'Tzero:':XSSdata[0],'Xslope:':XSSdata[1],'Yslope:':XSSdata[2],'Zslope:':XSSdata[3],
620        'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
621        'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
622        'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
623        'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
624    for iatm in range(Natoms):
625        waveTypes.append(parmDict[pfx+'waveType:'+str(iatm)])
626        for key in keys:
627            for m in range(Nwave[key[0]]):
628                parm = pfx+key+str(iatm)+':%d'%(m)
629                if parm in parmDict:
630                    keys[key][m][iatm] = parmDict[parm]
631    return np.array(waveTypes),FSSdata,XSSdata,USSdata,MSSdata
632   
633def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):
634   
635    Natoms = calcControls['Natoms'][pfx]
636    maxSSwave = calcControls['maxSSwave'][pfx]
637    Smult = np.zeros((Natoms,len(SGOps)))
638    TauT = np.zeros((Natoms,len(SGOps)))
639    for ix,xyz in enumerate(XData.T):
640        for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):
641            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)
642            Smult[ix][isym] = sdet
643            TauT[ix][isym] = tauT
644    return Smult,TauT
645   
646def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
647    ''' Compute structure factors for all h,k,l for phase
648    puts the result, F^2, in each ref[8] in refList
649    input:
650   
651    :param dict refDict: where
652        'RefList' list where each ref = h,k,l,m,d,...
653        'FF' dict of form factors - filed in below
654    :param np.array G:      reciprocal metric tensor
655    :param str pfx:    phase id string
656    :param dict SGData: space group info. dictionary output from SpcGroup
657    :param dict calcControls:
658    :param dict ParmDict:
659
660    '''       
661    phfx = pfx.split(':')[0]+hfx
662    ast = np.sqrt(np.diag(G))
663    Mast = twopisq*np.multiply.outer(ast,ast)
664    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
665    SGT = np.array([ops[1] for ops in SGData['SGOps']])
666    FFtables = calcControls['FFtables']
667    BLtables = calcControls['BLtables']
668    Flack = 1.0
669    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
670        Flack = 1.-2.*parmDict[phfx+'Flack']
671    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
672    TwDict = refDict.get('TwDict',{})           
673    if 'S' in calcControls[hfx+'histType']:
674        NTL = calcControls[phfx+'NTL']
675        NM = calcControls[phfx+'TwinNMN']+1
676        TwinLaw = calcControls[phfx+'TwinLaw']
677        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
678        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
679    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
680    FF = np.zeros(len(Tdata))
681    if 'NC' in calcControls[hfx+'histType']:
682        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
683    elif 'X' in calcControls[hfx+'histType']:
684        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
685        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
686    Uij = np.array(G2lat.U6toUij(Uijdata))
687    bij = Mast*Uij.T
688    blkSize = 100       #no. of reflections in a block
689    nRef = refDict['RefList'].shape[0]
690    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
691        if 'N' in calcControls[hfx+'histType']:
692            dat = G2el.getBLvalues(BLtables)
693            refDict['FF']['El'] = dat.keys()
694            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
695        else:       #'X'
696            dat = G2el.getFFvalues(FFtables,0.)
697            refDict['FF']['El'] = dat.keys()
698            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
699            for iref,ref in enumerate(refDict['RefList']):
700                SQ = 1./(2.*ref[4])**2
701                dat = G2el.getFFvalues(FFtables,SQ)
702                refDict['FF']['FF'][iref] *= dat.values()
703#reflection processing begins here - big arrays!
704    iBeg = 0
705    while iBeg < nRef:
706        iFin = min(iBeg+blkSize,nRef)
707        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
708        H = refl.T[:3]                          #array(blkSize,3)
709        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
710        TwMask = np.any(H,axis=-1)
711        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
712            for ir in range(blkSize):
713                iref = ir+iBeg
714                if iref in TwDict:
715                    for i in TwDict[iref]:
716                        for n in range(NTL):
717                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
718            TwMask = np.any(H,axis=-1)
719        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
720        SQfactor = 4.0*SQ*twopisq               #ditto prev.
721        if 'T' in calcControls[hfx+'histType']:
722            if 'P' in calcControls[hfx+'histType']:
723                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
724            else:
725                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
726            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
727            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
728        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
729        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
730        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
731        Uniq = np.inner(H,SGMT)
732        Phi = np.inner(H,SGT)
733        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
734        sinp = np.sin(phase)
735        cosp = np.cos(phase)
736        biso = -SQfactor*Uisodata[:,np.newaxis]
737        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
738        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
739        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
740        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
741        if 'T' in calcControls[hfx+'histType']:
742            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
743            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
744        else:
745            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
746            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
747        fas = np.sum(np.sum(fa,axis=-1),axis=-1)       #real sum over atoms & unique hkl
748        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag sum over atoms & uniq hkl
749        if SGData['SGInv']: #centrosymmetric; B=0
750            fbs[0] *= 0.
751        if 'P' in calcControls[hfx+'histType']:
752            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
753            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
754        else:
755            if len(TwinLaw) > 1:
756                refl.T[9] = np.sum(fas[:,:,0]**2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0)   #FcT from primary twin element
757                refl.T[7] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+   \
758                    np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
759                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
760            else:
761                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
762                refl.T[7] = np.copy(refl.T[9])               
763                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
764#        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
765        iBeg += blkSize
766   
767def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
768    'Needs a doc string'
769    phfx = pfx.split(':')[0]+hfx
770    ast = np.sqrt(np.diag(G))
771    Mast = twopisq*np.multiply.outer(ast,ast)
772    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
773    SGT = np.array([ops[1] for ops in SGData['SGOps']])
774    FFtables = calcControls['FFtables']
775    BLtables = calcControls['BLtables']
776    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
777    TwDict = refDict.get('TwDict',{})           
778    if 'S' in calcControls[hfx+'histType']:
779        NTL = calcControls[phfx+'NTL']
780        NM = calcControls[phfx+'TwinNMN']+1
781        TwinLaw = calcControls[phfx+'TwinLaw']
782        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
783        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
784    nTwin = len(TwinLaw)       
785    nRef = len(refDict['RefList'])
786    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
787    mSize = len(Mdata)
788    FF = np.zeros(len(Tdata))
789    if 'NC' in calcControls[hfx+'histType']:
790        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
791    elif 'X' in calcControls[hfx+'histType']:
792        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
793        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
794    Uij = np.array(G2lat.U6toUij(Uijdata))
795    bij = Mast*Uij.T
796    dFdvDict = {}
797    if nTwin > 1:
798        dFdfr = np.zeros((nRef,nTwin,mSize))
799        dFdx = np.zeros((nRef,nTwin,mSize,3))
800        dFdui = np.zeros((nRef,nTwin,mSize))
801        dFdua = np.zeros((nRef,nTwin,mSize,6))
802        dFdbab = np.zeros((nRef,nTwin,2))
803        dFdfl = np.zeros((nRef,nTwin))
804        dFdtw = np.zeros((nRef,nTwin))
805    else:
806        dFdfr = np.zeros((nRef,mSize))
807        dFdx = np.zeros((nRef,mSize,3))
808        dFdui = np.zeros((nRef,mSize))
809        dFdua = np.zeros((nRef,mSize,6))
810        dFdbab = np.zeros((nRef,2))
811        dFdfl = np.zeros((nRef))
812        dFdtw = np.zeros((nRef))
813    Flack = 1.0
814    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
815        Flack = 1.-2.*parmDict[phfx+'Flack']
816    time0 = time.time()
817    nref = len(refDict['RefList'])/100   
818    for iref,refl in enumerate(refDict['RefList']):
819        if 'T' in calcControls[hfx+'histType']:
820            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
821        H = np.array(refl[:3])
822        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
823        TwMask = np.any(H,axis=-1)
824        if TwinLaw.shape[0] > 1 and TwDict:
825            if iref in TwDict:
826                for i in TwDict[iref]:
827                    for n in range(NTL):
828                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
829            TwMask = np.any(H,axis=-1)
830        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
831        SQfactor = 8.0*SQ*np.pi**2
832        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
833        Bab = parmDict[phfx+'BabA']*dBabdA
834        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
835        FF = refDict['FF']['FF'][iref].T[Tindx].T
836        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
837        Phi = np.inner(H,SGT)
838        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
839        sinp = np.sin(phase)
840        cosp = np.cos(phase)
841        occ = Mdata*Fdata/len(SGT)
842        biso = -SQfactor*Uisodata[:,np.newaxis]
843        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
844        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
845        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
846        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
847        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
848        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
849        fot = (FF+FP-Bab)*Tcorr
850        fotp = FPP*Tcorr       
851        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
852        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
853        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
854        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
855        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
856        fbx = np.array([fot*cosp,-fotp*sinp])
857        #sum below is over Uniq
858        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
859        dfadba = np.sum(-cosp*Tcorr[:,np.newaxis],axis=1)
860        dfadui = np.sum(-SQfactor*fa,axis=-2)
861        if nTwin > 1:
862            dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
863            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
864            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
865        else:
866            dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
867            dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
868            # array(2,nAtom,3) & array(2,nAtom,6)
869        if not SGData['SGInv']:
870            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
871            dfbdba = np.sum(-sinp*Tcorr[:,np.newaxis],axis=1)
872            dfbdui = np.sum(-SQfactor*fb,axis=-2)
873            if len(TwinLaw) > 1:
874                dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])           
875                dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
876            else:
877                dfadfl = np.sum(-FPP*Tcorr*sinp)
878                dfbdfl = np.sum(FPP*Tcorr*cosp)
879                dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)           
880                dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)
881        else:
882            dfbdfr = np.zeros_like(dfadfr)
883            dfbdx = np.zeros_like(dfadx)
884            dfbdui = np.zeros_like(dfadui)
885            dfbdua = np.zeros_like(dfadua)
886            dfbdba = np.zeros_like(dfadba)
887            dfadfl = 0.0
888            dfbdfl = 0.0
889        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
890        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
891            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
892                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
893            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
894                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
895            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
896                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
897            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
898                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
899        else:
900            SA = fas[0]-fbs[1]
901            SB = fbs[0]+fas[1]
902            if nTwin > 1:
903                dFdfr[iref] = [2.*TwMask[it]*SA[it]*(dfadfr[0][it]+dfbdfr[1][it])*Mdata/len(Uniq[it])+ \
904                    2.*TwMask[it]*SB[it]*(dfbdfr[0][it]+dfadfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
905                dFdx[iref] = [2.*TwMask[it]*SA[it]*(dfadx[it][0]+dfbdx[it][1])+2.*TwMask[it]*SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)]
906                dFdui[iref] = [2.*TwMask[it]*SA[it]*(dfadui[0][it]+dfbdui[1][it])+2.*TwMask[it]*SB[it]*(dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)]
907                dFdua[iref] = [2.*TwMask[it]*SA[it]*(dfadua[it][0]+dfbdua[it][1])+2.*TwMask[it]*SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)]
908                dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
909               
910            else:   #these are good for no twin single crystals
911                dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
912                    2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
913                dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])    #array(nRef,nAtom,3)
914                dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])   #array(nRef,nAtom)
915                dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])    #array(nRef,nAtom,6)
916                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
917        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
918            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
919#        GSASIIpath.IPyBreak()
920    print ' derivative time %.4f\r'%(time.time()-time0)
921        #loop over atoms - each dict entry is list of derivatives for all the reflections
922    if nTwin > 1:
923        for i in range(len(Mdata)):     #these all OK?
924            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0)
925            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0)
926            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0)
927            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0)
928            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0)
929            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0)
930            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0)
931            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0)
932            dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5*dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0)
933            dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5*dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0)
934            dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5*dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0)
935    else:
936        for i in range(len(Mdata)):
937            dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
938            dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
939            dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
940            dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
941            dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
942            dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
943            dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
944            dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
945            dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i]
946            dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i]
947            dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i]
948        dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
949    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
950    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
951    if nTwin > 1:
952        for i in range(nTwin):
953            dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
954    return dFdvDict
955   
956def SStructureFactor2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
957    '''
958    Compute super structure factors for all h,k,l,m for phase
959    puts the result, F^2, in each ref[8+im] in refList
960    input:
961   
962    :param dict refDict: where
963        'RefList' list where each ref = h,k,l,m,d,...
964        'FF' dict of form factors - filed in below
965    :param np.array G:      reciprocal metric tensor
966    :param str pfx:    phase id string
967    :param dict SGData: space group info. dictionary output from SpcGroup
968    :param dict calcControls:
969    :param dict ParmDict:
970
971    '''
972    nxs = np.newaxis       
973    phfx = pfx.split(':')[0]+hfx
974    ast = np.sqrt(np.diag(G))
975    Mast = twopisq*np.multiply.outer(ast,ast)   
976    SGInv = SGData['SGInv']
977    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
978    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
979    FFtables = calcControls['FFtables']
980    BLtables = calcControls['BLtables']
981    Flack = 1.0
982    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
983        Flack = 1.-2.*parmDict[phfx+'Flack']
984    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
985    TwDict = refDict.get('TwDict',{})           
986    if 'S' in calcControls[hfx+'histType']:
987        NTL = calcControls[phfx+'NTL']
988        NM = calcControls[phfx+'TwinNMN']+1
989        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
990        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
991        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
992    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
993    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
994    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
995    FF = np.zeros(len(Tdata))
996    if 'NC' in calcControls[hfx+'histType']:
997        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
998    else:
999        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1000        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1001    Uij = np.array(G2lat.U6toUij(Uijdata))
1002    bij = Mast*Uij.T
1003    blkSize = 100       #no. of reflections in a block
1004    nRef = refDict['RefList'].shape[0]
1005    if not len(refDict['FF']):
1006        if 'N' in calcControls[hfx+'histType']:
1007            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1008            refDict['FF']['El'] = dat.keys()
1009            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
1010        else:
1011            dat = G2el.getFFvalues(FFtables,0.)       
1012            refDict['FF']['El'] = dat.keys()
1013            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
1014            for iref,ref in enumerate(refDict['RefList']):
1015                SQ = 1./(2.*ref[5])**2
1016                dat = G2el.getFFvalues(FFtables,SQ)
1017                refDict['FF']['FF'][iref] *= dat.values()
1018    time0 = time.time()
1019#reflection processing begins here - big arrays!
1020    iBeg = 0
1021    while iBeg < nRef:
1022        iFin = min(iBeg+blkSize,nRef)
1023        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1024        H = refl.T[:4]                          #array(blkSize,4)
1025#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,4) or (blkSize,4)
1026#        TwMask = np.any(H,axis=-1)
1027#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1028#            for ir in range(blkSize):
1029#                iref = ir+iBeg
1030#                if iref in TwDict:
1031#                    for i in TwDict[iref]:
1032#                        for n in range(NTL):
1033#                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1034#            TwMask = np.any(H,axis=-1)
1035        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1036        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1037        Uniq = np.inner(H.T,SSGMT)
1038        Phi = np.inner(H.T,SSGT)
1039        if SGInv:   #if centro - expand HKL sets
1040            Uniq = np.hstack((Uniq,-Uniq))
1041            Phi = np.hstack((Phi,-Phi))
1042        if 'T' in calcControls[hfx+'histType']:
1043            if 'P' in calcControls[hfx+'histType']:
1044                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1045            else:
1046                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1047            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1048            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1049        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1050        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1051        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1052        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))+Phi[:,:,nxs])
1053        sinp = np.sin(phase)
1054        cosp = np.cos(phase)
1055        biso = -SQfactor*Uisodata[:,nxs]
1056        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1057        HbH = -np.sum(Uniq[:,:,nxs,:3]*np.inner(Uniq[:,:,:3],bij),axis=-1)
1058        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1059        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]
1060        if 'T' in calcControls[hfx+'histType']:
1061            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1062            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1063        else:
1064            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1065            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1066        GfpuA = G2mth.Modulation(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) #2 x refBlk x sym X atoms
1067        fag = fa*GfpuA[0]-fb*GfpuA[1]
1068        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1069#        GSASIIpath.IPyBreak()
1070        fas = np.sum(np.sum(fag,axis=-1),axis=-1)
1071        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1072        if 'P' in calcControls[hfx+'histType']:
1073            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
1074            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1075        else:
1076            if len(TwinLaw) > 1:
1077                refl.T[10] = np.sum(fas[:,:,0]**2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0)   #FcT from primary twin element
1078                refl.T[8] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+   \
1079                    np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
1080                refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1081            else:
1082                refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
1083                refl.T[8] = np.copy(refl.T[10])               
1084                refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1085        iBeg += blkSize
1086    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
1087
1088def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1089    'Needs a doc string'
1090    nxs = np.newaxis       
1091    phfx = pfx.split(':')[0]+hfx
1092    ast = np.sqrt(np.diag(G))
1093    Mast = twopisq*np.multiply.outer(ast,ast)
1094    SGInv = SGData['SGInv']
1095    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1096    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1097    FFtables = calcControls['FFtables']
1098    BLtables = calcControls['BLtables']
1099    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
1100    TwDict = refDict.get('TwDict',{})           
1101    if 'S' in calcControls[hfx+'histType']:
1102        NTL = calcControls[phfx+'NTL']
1103        NM = calcControls[phfx+'TwinNMN']+1
1104        TwinLaw = calcControls[phfx+'TwinLaw']
1105        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1106        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1107    nTwin = len(TwinLaw)       
1108    nRef = len(refDict['RefList'])
1109    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1110    mSize = len(Mdata)  #no. atoms
1111    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1112    FF = np.zeros(len(Tdata))
1113    if 'NC' in calcControls[hfx+'histType']:
1114        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1115    elif 'X' in calcControls[hfx+'histType']:
1116        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1117        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1118    Uij = np.array(G2lat.U6toUij(Uijdata))
1119    bij = Mast*Uij.T
1120    if not len(refDict['FF']):
1121        if 'N' in calcControls[hfx+'histType']:
1122            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1123        else:
1124            dat = G2el.getFFvalues(FFtables,0.)       
1125        refDict['FF']['El'] = dat.keys()
1126        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1127    dFdvDict = {}
1128    if nTwin > 1:
1129        dFdfr = np.zeros((nRef,nTwin,mSize))
1130        dFdx = np.zeros((nRef,nTwin,mSize,3))
1131        dFdui = np.zeros((nRef,nTwin,mSize))
1132        dFdua = np.zeros((nRef,nTwin,mSize,6))
1133        dFdbab = np.zeros((nRef,nTwin,2))
1134        dFdfl = np.zeros((nRef,nTwin))
1135        dFdtw = np.zeros((nRef,nTwin))
1136        dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1],2))
1137        dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],6))
1138        dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],12))
1139    else:
1140        dFdfr = np.zeros((nRef,mSize))
1141        dFdx = np.zeros((nRef,mSize,3))
1142        dFdui = np.zeros((nRef,mSize))
1143        dFdua = np.zeros((nRef,mSize,6))
1144        dFdbab = np.zeros((nRef,2))
1145        dFdfl = np.zeros((nRef))
1146        dFdtw = np.zeros((nRef))
1147        dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1148        dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1149        dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1150    Flack = 1.0
1151    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1152        Flack = 1.-2.*parmDict[phfx+'Flack']
1153    time0 = time.time()
1154    nRef = len(refDict['RefList'])/100
1155    for iref,refl in enumerate(refDict['RefList']):
1156        if 'T' in calcControls[hfx+'histType']:
1157            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1158        H = np.array(refl[:4])
1159#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(4,nTwins) or (4)
1160#        TwMask = np.any(H,axis=-1)
1161#        if TwinLaw.shape[0] > 1 and TwDict:
1162#            if iref in TwDict:
1163#                for i in TwDict[iref]:
1164#                    for n in range(NTL):
1165#                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1166#            TwMask = np.any(H,axis=-1)
1167        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1168        SQfactor = 8.0*SQ*np.pi**2
1169        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1170        Bab = parmDict[phfx+'BabA']*dBabdA
1171        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1172        FF = refDict['FF']['FF'][iref].T[Tindx]
1173        Uniq = np.inner(H,SSGMT)
1174        Phi = np.inner(H,SSGT)
1175        if SGInv:   #if centro - expand HKL sets
1176            Uniq = np.vstack((Uniq,-Uniq))
1177            Phi = np.hstack((Phi,-Phi))
1178        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1179        sinp = np.sin(phase)
1180        cosp = np.cos(phase)
1181        occ = Mdata*Fdata/Uniq.shape[0]
1182        biso = -SQfactor*Uisodata[:,nxs]
1183        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
1184        HbH = -np.sum(Uniq[:,nxs,:3]*np.inner(Uniq[:,:3],bij),axis=-1)  #ops x atoms
1185        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in Uniq]) #atoms x 3x3
1186        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
1187        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1188        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1189        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1190        fotp = FPP*Tcorr            #ops x atoms
1191        GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast)
1192        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
1193        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
1194        fag = fa*GfpuA[0]-fb*GfpuA[1]
1195        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1196       
1197        fas = np.sum(np.sum(fag,axis=1),axis=1)
1198        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1199        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
1200        fbx = np.array([fot*cosp,-fotp*sinp])
1201        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1202        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1203        #sum below is over Uniq
1204        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1205        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1206        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1207        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1208        dfadui = np.sum(-SQfactor*fag,axis=1)
1209        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1210        if nTwin > 1:
1211            dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
1212            dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])           
1213            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
1214            dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
1215            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
1216        else:
1217            dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
1218            dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=-2)           
1219            dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,np.newaxis],axis=-2)
1220            dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,np.newaxis],axis=-2)
1221            # array(2,nAtom,3) & array(2,nAtom,6)
1222        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
1223#        GSASIIpath.IPyBreak()
1224        if not SGData['SGInv'] and len(TwinLaw) == 1:   #Flack derivative
1225            dfadfl = np.sum(-FPP*Tcorr*sinp)
1226            dfbdfl = np.sum(FPP*Tcorr*cosp)
1227        else:
1228            dfadfl = 1.0
1229            dfbdfl = 1.0
1230        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1231        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
1232            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1233                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1234            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1235                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1236            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1237                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1238            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1239                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1240        else:
1241            SA = fas[0]-fbs[1]
1242            SB = fbs[0]+fas[1]
1243            if nTwin > 1:
1244                dFdfr[iref] = [2.*TwMask[it]*SA[it]*(dfadfr[0][it]+dfbdfr[1][it])*Mdata/len(Uniq[it])+ \
1245                    2.*TwMask[it]*SB[it]*(dfbdfr[0][it]+dfadfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
1246                dFdx[iref] = [2.*TwMask[it]*SA[it]*(dfadx[it][0]+dfbdx[it][1])+2.*TwMask[it]*SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)]
1247                dFdui[iref] = [2.*TwMask[it]*SA[it]*(dfadui[0][it]+dfbdui[1][it])+2.*TwMask[it]*SB[it]*(dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)]
1248                dFdua[iref] = [2.*TwMask[it]*SA[it]*(dfadua[it][0]+dfbdua[it][1])+2.*TwMask[it]*SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)]
1249                dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
1250               
1251            else:   #these are good for no twin single crystals
1252                dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
1253                    2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1254                dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1])    #array(nRef,nAtom,3)
1255                dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])   #array(nRef,nAtom)
1256                dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])    #array(nRef,nAtom,6)
1257                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
1258        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1259            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1260        #loop over atoms - each dict entry is list of derivatives for all the reflections
1261    for i in range(len(Mdata)):     
1262        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1263        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1264        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1265        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1266        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1267        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1268        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1269        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1270        dFdvDict[pfx+'AU12:'+str(i)] = .5*dFdua.T[3][i]
1271        dFdvDict[pfx+'AU13:'+str(i)] = .5*dFdua.T[4][i]
1272        dFdvDict[pfx+'AU23:'+str(i)] = .5*dFdua.T[5][i]
1273        #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
1274    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1275    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1276    return dFdvDict
1277   
1278def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
1279    ''' Single crystal extinction function; returns extinction & derivative
1280    '''
1281    extCor = 1.0
1282    dervDict = {}
1283    dervCor = 1.0
1284    if calcControls[phfx+'EType'] != 'None':
1285        SQ = 1/(4.*ref[4+im]**2)
1286        if 'C' in parmDict[hfx+'Type']:           
1287            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
1288        else:   #'T'
1289            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
1290        if 'SXC' in parmDict[hfx+'Type']:
1291            AV = 7.9406e5/parmDict[pfx+'Vol']**2
1292            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1293            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
1294            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
1295        elif 'SNT' in parmDict[hfx+'Type']:
1296            AV = 1.e7/parmDict[pfx+'Vol']**2
1297            PL = SQ
1298            PLZ = AV*ref[9+im]*ref[12+im]**2
1299        elif 'SNC' in parmDict[hfx+'Type']:
1300            AV = 1.e7/parmDict[pfx+'Vol']**2
1301            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
1302            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
1303           
1304        if 'Primary' in calcControls[phfx+'EType']:
1305            PLZ *= 1.5
1306        else:
1307            if 'C' in parmDict[hfx+'Type']:
1308                PLZ *= calcControls[phfx+'Tbar']
1309            else: #'T'
1310                PLZ *= ref[13+im]      #t-bar
1311        if 'Primary' in calcControls[phfx+'EType']:
1312            PLZ *= 1.5
1313            PSIG = parmDict[phfx+'Ep']
1314        elif 'I & II' in calcControls[phfx+'EType']:
1315            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
1316        elif 'Type II' in calcControls[phfx+'EType']:
1317            PSIG = parmDict[phfx+'Es']
1318        else:       # 'Secondary Type I'
1319            PSIG = parmDict[phfx+'Eg']/PL
1320           
1321        AG = 0.58+0.48*cos2T+0.24*cos2T**2
1322        AL = 0.025+0.285*cos2T
1323        BG = 0.02-0.025*cos2T
1324        BL = 0.15-0.2*(0.75-cos2T)**2
1325        if cos2T < 0.:
1326            BL = -0.45*cos2T
1327        CG = 2.
1328        CL = 2.
1329        PF = PLZ*PSIG
1330       
1331        if 'Gaussian' in calcControls[phfx+'EApprox']:
1332            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
1333            extCor = np.sqrt(PF4)
1334            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
1335        else:
1336            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
1337            extCor = np.sqrt(PF4)
1338            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
1339
1340        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
1341        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
1342            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
1343        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
1344            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
1345        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
1346            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
1347               
1348    return 1./extCor,dervDict,dervCor
1349   
1350def Dict2Values(parmdict, varylist):
1351    '''Use before call to leastsq to setup list of values for the parameters
1352    in parmdict, as selected by key in varylist'''
1353    return [parmdict[key] for key in varylist] 
1354   
1355def Values2Dict(parmdict, varylist, values):
1356    ''' Use after call to leastsq to update the parameter dictionary with
1357    values corresponding to keys in varylist'''
1358    parmdict.update(zip(varylist,values))
1359   
1360def GetNewCellParms(parmDict,varyList):
1361    'Needs a doc string'
1362    newCellDict = {}
1363    Anames = ['A'+str(i) for i in range(6)]
1364    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
1365    for item in varyList:
1366        keys = item.split(':')
1367        if keys[2] in Ddict:
1368            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
1369            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
1370            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1371    return newCellDict          # is e.g. {'0::D11':A0-D11}
1372   
1373def ApplyXYZshifts(parmDict,varyList):
1374    '''
1375    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1376   
1377    :param dict parmDict: parameter dictionary
1378    :param list varyList: list of variables (not used!)
1379    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
1380
1381    '''
1382    newAtomDict = {}
1383    for item in parmDict:
1384        if 'dA' in item:
1385            parm = ''.join(item.split('d'))
1386            parmDict[parm] += parmDict[item]
1387            newAtomDict[item] = [parm,parmDict[parm]]
1388    return newAtomDict
1389   
1390def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1391    'Spherical harmonics texture'
1392    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1393    if 'T' in calcControls[hfx+'histType']:
1394        tth = parmDict[hfx+'2-theta']
1395    else:
1396        tth = refl[5+im]
1397    odfCor = 1.0
1398    H = refl[:3]
1399    cell = G2lat.Gmat2cell(g)
1400    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1401    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1402    phi,beta = G2lat.CrsAng(H,cell,SGData)
1403    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1404    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1405    for item in SHnames:
1406        L,M,N = eval(item.strip('C'))
1407        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1408        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1409        Lnorm = G2lat.Lnorm(L)
1410        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1411    return odfCor
1412   
1413def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
1414    'Spherical harmonics texture derivatives'
1415    if 'T' in calcControls[hfx+'histType']:
1416        tth = parmDict[hfx+'2-theta']
1417    else:
1418        tth = refl[5+im]
1419    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1420    odfCor = 1.0
1421    dFdODF = {}
1422    dFdSA = [0,0,0]
1423    H = refl[:3]
1424    cell = G2lat.Gmat2cell(g)
1425    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1426    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1427    phi,beta = G2lat.CrsAng(H,cell,SGData)
1428    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
1429    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1430    for item in SHnames:
1431        L,M,N = eval(item.strip('C'))
1432        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1433        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1434        Lnorm = G2lat.Lnorm(L)
1435        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1436        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1437        for i in range(3):
1438            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1439    return odfCor,dFdODF,dFdSA
1440   
1441def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1442    'spherical harmonics preferred orientation (cylindrical symmetry only)'
1443    if 'T' in calcControls[hfx+'histType']:
1444        tth = parmDict[hfx+'2-theta']
1445    else:
1446        tth = refl[5+im]
1447    odfCor = 1.0
1448    H = refl[:3]
1449    cell = G2lat.Gmat2cell(g)
1450    Sangls = [0.,0.,0.]
1451    if 'Bragg' in calcControls[hfx+'instType']:
1452        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1453        IFCoup = True
1454    else:
1455        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1456        IFCoup = False
1457    phi,beta = G2lat.CrsAng(H,cell,SGData)
1458    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1459    SHnames = calcControls[phfx+'SHnames']
1460    for item in SHnames:
1461        L,N = eval(item.strip('C'))
1462        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1463        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1464        Lnorm = G2lat.Lnorm(L)
1465        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1466    return np.squeeze(odfCor)
1467   
1468def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
1469    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
1470    if 'T' in calcControls[hfx+'histType']:
1471        tth = parmDict[hfx+'2-theta']
1472    else:
1473        tth = refl[5+im]
1474    odfCor = 1.0
1475    dFdODF = {}
1476    H = refl[:3]
1477    cell = G2lat.Gmat2cell(g)
1478    Sangls = [0.,0.,0.]
1479    if 'Bragg' in calcControls[hfx+'instType']:
1480        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1481        IFCoup = True
1482    else:
1483        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
1484        IFCoup = False
1485    phi,beta = G2lat.CrsAng(H,cell,SGData)
1486    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1487    SHnames = calcControls[phfx+'SHnames']
1488    for item in SHnames:
1489        L,N = eval(item.strip('C'))
1490        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1491        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
1492        Lnorm = G2lat.Lnorm(L)
1493        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
1494        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
1495    return odfCor,dFdODF
1496   
1497def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1498    'March-Dollase preferred orientation correction'
1499    POcorr = 1.0
1500    MD = parmDict[phfx+'MD']
1501    if MD != 1.0:
1502        MDAxis = calcControls[phfx+'MDAxis']
1503        sumMD = 0
1504        for H in uniq:           
1505            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1506            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1507            sumMD += A**3
1508        POcorr = sumMD/len(uniq)
1509    return POcorr
1510   
1511def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
1512    'Needs a doc string'
1513    POcorr = 1.0
1514    POderv = {}
1515    if calcControls[phfx+'poType'] == 'MD':
1516        MD = parmDict[phfx+'MD']
1517        MDAxis = calcControls[phfx+'MDAxis']
1518        sumMD = 0
1519        sumdMD = 0
1520        for H in uniq:           
1521            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1522            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1523            sumMD += A**3
1524            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1525        POcorr = sumMD/len(uniq)
1526        POderv[phfx+'MD'] = sumdMD/len(uniq)
1527    else:   #spherical harmonics
1528        if calcControls[phfx+'SHord']:
1529            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1530    return POcorr,POderv
1531   
1532def GetAbsorb(refl,im,hfx,calcControls,parmDict):
1533    'Needs a doc string'
1534    if 'Debye' in calcControls[hfx+'instType']:
1535        if 'T' in calcControls[hfx+'histType']:
1536            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1537        else:
1538            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1539    else:
1540        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
1541   
1542def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
1543    'Needs a doc string'
1544    if 'Debye' in calcControls[hfx+'instType']:
1545        if 'T' in calcControls[hfx+'histType']:
1546            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
1547        else:
1548            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
1549    else:
1550        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
1551       
1552def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1553    'Needs a doc string'
1554    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1555    pi2 = np.sqrt(2./np.pi)
1556    if 'T' in calcControls[hfx+'histType']:
1557        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1558        wave = refl[14+im]
1559    else:   #'C'W
1560        sth2 = sind(refl[5+im]/2.)**2
1561        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1562    c2th = 1.-2.0*sth2
1563    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1564    if 'X' in calcControls[hfx+'histType']:
1565        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1566    xfac = flv2*parmDict[phfx+'Extinction']
1567    exb = 1.0
1568    if xfac > -1.:
1569        exb = 1./np.sqrt(1.+xfac)
1570    exl = 1.0
1571    if 0 < xfac <= 1.:
1572        xn = np.array([xfac**(i+1) for i in range(6)])
1573        exl += np.sum(xn*coef)
1574    elif xfac > 1.:
1575        xfac2 = 1./np.sqrt(xfac)
1576        exl = pi2*(1.-0.125/xfac)*xfac2
1577    return exb*sth2+exl*(1.-sth2)
1578   
1579def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
1580    'Needs a doc string'
1581    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
1582    pi2 = np.sqrt(2./np.pi)
1583    if 'T' in calcControls[hfx+'histType']:
1584        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
1585        wave = refl[14+im]
1586    else:   #'C'W
1587        sth2 = sind(refl[5+im]/2.)**2
1588        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1589    c2th = 1.-2.0*sth2
1590    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
1591    if 'X' in calcControls[hfx+'histType']:
1592        flv2 *= 0.079411*(1.0+c2th**2)/2.0
1593    xfac = flv2*parmDict[phfx+'Extinction']
1594    dbde = -500.*flv2
1595    if xfac > -1.:
1596        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
1597    dlde = 0.
1598    if 0 < xfac <= 1.:
1599        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
1600        dlde = np.sum(xn*coef)
1601    elif xfac > 1.:
1602        xfac2 = 1./np.sqrt(xfac)
1603        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
1604       
1605    return dbde*sth2+dlde*(1.-sth2)
1606   
1607def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1608    'Needs a doc string'    #need powder extinction!
1609    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
1610    if 'X' in parmDict[hfx+'Type']:
1611        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
1612    POcorr = 1.0
1613    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
1614        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1615    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
1616        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
1617    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
1618        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
1619    Icorr *= POcorr
1620    AbsCorr = 1.0
1621    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
1622    Icorr *= AbsCorr
1623    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
1624    Icorr *= ExtCorr
1625    return Icorr,POcorr,AbsCorr,ExtCorr
1626   
1627def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1628    'Needs a doc string'    #need powder extinction derivs!
1629    dIdsh = 1./parmDict[hfx+'Scale']
1630    dIdsp = 1./parmDict[phfx+'Scale']
1631    if 'X' in parmDict[hfx+'Type']:
1632        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
1633        dIdPola /= pola
1634    else:       #'N'
1635        dIdPola = 0.0
1636    dFdODF = {}
1637    dFdSA = [0,0,0]
1638    dIdPO = {}
1639    if pfx+'SHorder' in parmDict:
1640        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
1641        for iSH in dFdODF:
1642            dFdODF[iSH] /= odfCor
1643        for i in range(3):
1644            dFdSA[i] /= odfCor
1645    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
1646        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
1647        for iPO in dIdPO:
1648            dIdPO[iPO] /= POcorr
1649    if 'T' in parmDict[hfx+'Type']:
1650        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
1651        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
1652    else:
1653        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
1654        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
1655    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
1656       
1657def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1658    'Needs a doc string'
1659    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
1660        costh = cosd(refl[5+im]/2.)
1661        #crystallite size
1662        if calcControls[phfx+'SizeType'] == 'isotropic':
1663            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1664        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1665            H = np.array(refl[:3])
1666            P = np.array(calcControls[phfx+'SizeAxis'])
1667            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1668            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
1669            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1670        else:           #ellipsoidal crystallites
1671            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1672            H = np.array(refl[:3])
1673            lenR = G2pwd.ellipseSize(H,Sij,GB)
1674            Sgam = 1.8*wave/(np.pi*costh*lenR)
1675        #microstrain               
1676        if calcControls[phfx+'MustrainType'] == 'isotropic':
1677            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1678        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1679            H = np.array(refl[:3])
1680            P = np.array(calcControls[phfx+'MustrainAxis'])
1681            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1682            Si = parmDict[phfx+'Mustrain;i']
1683            Sa = parmDict[phfx+'Mustrain;a']
1684            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1685        else:       #generalized - P.W. Stephens model
1686            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1687            Sum = 0
1688            for i,strm in enumerate(Strms):
1689                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1690            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
1691    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
1692        #crystallite size
1693        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1694            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1695        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1696            H = np.array(refl[:3])
1697            P = np.array(calcControls[phfx+'SizeAxis'])
1698            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1699            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
1700            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
1701        else:           #ellipsoidal crystallites   #OK
1702            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1703            H = np.array(refl[:3])
1704            lenR = G2pwd.ellipseSize(H,Sij,GB)
1705            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
1706        #microstrain               
1707        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
1708            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1709        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
1710            H = np.array(refl[:3])
1711            P = np.array(calcControls[phfx+'MustrainAxis'])
1712            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1713            Si = parmDict[phfx+'Mustrain;i']
1714            Sa = parmDict[phfx+'Mustrain;a']
1715            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1716        else:       #generalized - P.W. Stephens model  OK
1717            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1718            Sum = 0
1719            for i,strm in enumerate(Strms):
1720                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1721            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
1722           
1723    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
1724    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
1725    sig /= ateln2
1726    return sig,gam
1727       
1728def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
1729    'Needs a doc string'
1730    gamDict = {}
1731    sigDict = {}
1732    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
1733        costh = cosd(refl[5+im]/2.)
1734        tanth = tand(refl[5+im]/2.)
1735        #crystallite size derivatives
1736        if calcControls[phfx+'SizeType'] == 'isotropic':
1737            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
1738            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
1739            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
1740        elif calcControls[phfx+'SizeType'] == 'uniaxial':
1741            H = np.array(refl[:3])
1742            P = np.array(calcControls[phfx+'SizeAxis'])
1743            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1744            Si = parmDict[phfx+'Size;i']
1745            Sa = parmDict[phfx+'Size;a']
1746            gami = 1.8*wave/(costh*np.pi*Si*Sa)
1747            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1748            Sgam = gami*sqtrm
1749            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1750            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1751            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1752            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1753            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1754            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1755        else:           #ellipsoidal crystallites
1756            const = 1.8*wave/(np.pi*costh)
1757            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1758            H = np.array(refl[:3])
1759            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1760            Sgam = const/lenR
1761            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1762                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1763                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1764        gamDict[phfx+'Size;mx'] = Sgam
1765        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
1766               
1767        #microstrain derivatives               
1768        if calcControls[phfx+'MustrainType'] == 'isotropic':
1769            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
1770            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
1771            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
1772        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1773            H = np.array(refl[:3])
1774            P = np.array(calcControls[phfx+'MustrainAxis'])
1775            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1776            Si = parmDict[phfx+'Mustrain;i']
1777            Sa = parmDict[phfx+'Mustrain;a']
1778            gami = 0.018*Si*Sa*tanth/np.pi
1779            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1780            Mgam = gami/sqtrm
1781            dsi = -gami*Si*cosP**2/sqtrm**3
1782            dsa = -gami*Sa*sinP**2/sqtrm**3
1783            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1784            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1785            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1786            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1787        else:       #generalized - P.W. Stephens model
1788            const = 0.018*refl[4+im]**2*tanth/np.pi
1789            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1790            Sum = 0
1791            for i,strm in enumerate(Strms):
1792                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1793                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1794                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1795            Mgam = const*np.sqrt(Sum)
1796            for i in range(len(Strms)):
1797                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1798                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2
1799        gamDict[phfx+'Mustrain;mx'] = Mgam
1800        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1801    else:   #'T'OF - All checked & OK
1802        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
1803            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
1804            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
1805            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
1806        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
1807            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1808            H = np.array(refl[:3])
1809            P = np.array(calcControls[phfx+'SizeAxis'])
1810            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1811            Si = parmDict[phfx+'Size;i']
1812            Sa = parmDict[phfx+'Size;a']
1813            gami = const/(Si*Sa)
1814            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1815            Sgam = gami*sqtrm
1816            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
1817            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
1818            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
1819            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
1820            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1821            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1822        else:           #OK  ellipsoidal crystallites
1823            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
1824            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1825            H = np.array(refl[:3])
1826            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1827            Sgam = const/lenR
1828            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1829                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
1830                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
1831        gamDict[phfx+'Size;mx'] = Sgam  #OK
1832        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
1833               
1834        #microstrain derivatives               
1835        if calcControls[phfx+'MustrainType'] == 'isotropic':
1836            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
1837            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
1838            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
1839        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1840            H = np.array(refl[:3])
1841            P = np.array(calcControls[phfx+'MustrainAxis'])
1842            cosP,sinP = G2lat.CosSinAngle(H,P,G)
1843            Si = parmDict[phfx+'Mustrain;i']
1844            Sa = parmDict[phfx+'Mustrain;a']
1845            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
1846            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1847            Mgam = gami/sqtrm
1848            dsi = -gami*Si*cosP**2/sqtrm**3
1849            dsa = -gami*Sa*sinP**2/sqtrm**3
1850            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
1851            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
1852            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
1853            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
1854        else:       #generalized - P.W. Stephens model OK
1855            pwrs = calcControls[phfx+'MuPwrs']
1856            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
1857            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
1858            Sum = 0
1859            for i,strm in enumerate(Strms):
1860                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
1861                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
1862                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
1863            Mgam = const*np.sqrt(Sum)
1864            for i in range(len(Strms)):
1865                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
1866                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
1867        gamDict[phfx+'Mustrain;mx'] = Mgam
1868        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
1869       
1870    return sigDict,gamDict
1871       
1872def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1873    'Needs a doc string'
1874    if im:
1875        h,k,l,m = refl[:4]
1876        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1877        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
1878    else:
1879        h,k,l = refl[:3]
1880        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
1881    refl[4+im] = d
1882    if 'C' in calcControls[hfx+'histType']:
1883        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1884        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1885        if 'Bragg' in calcControls[hfx+'instType']:
1886            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1887                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1888        else:               #Debye-Scherrer - simple but maybe not right
1889            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1890    elif 'T' in calcControls[hfx+'histType']:
1891        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
1892        #do I need sample position effects - maybe?
1893    return pos
1894
1895def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
1896    'Needs a doc string'
1897    dpr = 180./np.pi
1898    if im:
1899        h,k,l,m = refl[:4]
1900        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1901        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
1902    else:
1903        m = 0
1904        h,k,l = refl[:3]       
1905        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1906    dst = np.sqrt(dstsq)
1907    dsp = 1./dst
1908    if 'C' in calcControls[hfx+'histType']:
1909        pos = refl[5+im]-parmDict[hfx+'Zero']
1910        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1911        dpdw = const*dst
1912        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
1913        dpdZ = 1.0
1914        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1915            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
1916        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1917        if 'Bragg' in calcControls[hfx+'instType']:
1918            dpdSh = -4.*shft*cosd(pos/2.0)
1919            dpdTr = -shft*sind(pos)*100.0
1920            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
1921        else:               #Debye-Scherrer - simple but maybe not right
1922            dpdXd = -shft*cosd(pos)
1923            dpdYd = -shft*sind(pos)
1924            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
1925    elif 'T' in calcControls[hfx+'histType']:
1926        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
1927        dpdZ = 1.0
1928        dpdDC = dsp
1929        dpdDA = dsp**2
1930        dpdDB = 1./dsp
1931        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
1932            2*l*A[2]+h*A[4]+k*A[5]])*m**parmDict[hfx+'difC']*dsp**3/2.
1933        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
1934           
1935def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1936    'Needs a doc string'
1937    laue = SGData['SGLaue']
1938    uniq = SGData['SGUniq']
1939    h,k,l = refl[:3]
1940    if laue in ['m3','m3m']:
1941        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1942            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1943    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1944        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1945    elif laue in ['3R','3mR']:
1946        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1947    elif laue in ['4/m','4/mmm']:
1948        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1949    elif laue in ['mmm']:
1950        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1951    elif laue in ['2/m']:
1952        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1953        if uniq == 'a':
1954            Dij += parmDict[phfx+'D23']*k*l
1955        elif uniq == 'b':
1956            Dij += parmDict[phfx+'D13']*h*l
1957        elif uniq == 'c':
1958            Dij += parmDict[phfx+'D12']*h*k
1959    else:
1960        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1961            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1962    if 'C' in calcControls[hfx+'histType']:
1963        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1964    else:
1965        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
1966           
1967def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
1968    'Needs a doc string'
1969    laue = SGData['SGLaue']
1970    uniq = SGData['SGUniq']
1971    h,k,l = refl[:3]
1972    if laue in ['m3','m3m']:
1973        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1974            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1975    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1976        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1977    elif laue in ['3R','3mR']:
1978        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1979    elif laue in ['4/m','4/mmm']:
1980        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1981    elif laue in ['mmm']:
1982        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1983    elif laue in ['2/m']:
1984        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1985        if uniq == 'a':
1986            dDijDict[phfx+'D23'] = k*l
1987        elif uniq == 'b':
1988            dDijDict[phfx+'D13'] = h*l
1989        elif uniq == 'c':
1990            dDijDict[phfx+'D12'] = h*k
1991    else:
1992        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1993            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1994    if 'C' in calcControls[hfx+'histType']:
1995        for item in dDijDict:
1996            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
1997    else:
1998        for item in dDijDict:
1999            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
2000    return dDijDict
2001   
2002def GetDij(phfx,SGData,parmDict):
2003    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2004    return G2spc.HStrainVals(HSvals,SGData)
2005               
2006def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2007    'Needs a doc string'
2008    histoList = Histograms.keys()
2009    histoList.sort()
2010    for histogram in histoList:
2011        if 'PWDR' in histogram[:4]:
2012            Histogram = Histograms[histogram]
2013            hId = Histogram['hId']
2014            hfx = ':%d:'%(hId)
2015            Limits = calcControls[hfx+'Limits']
2016            if 'C' in calcControls[hfx+'histType']:
2017                shl = max(parmDict[hfx+'SH/L'],0.0005)
2018                Ka2 = False
2019                kRatio = 0.0
2020                if hfx+'Lam1' in parmDict.keys():
2021                    Ka2 = True
2022                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2023                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
2024            x,y,w,yc,yb,yd = Histogram['Data']
2025            xB = np.searchsorted(x,Limits[0])
2026            xF = np.searchsorted(x,Limits[1])
2027            ymb = np.array(y-yb)
2028            ymb = np.where(ymb,ymb,1.0)
2029            ycmb = np.array(yc-yb)
2030            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
2031            refLists = Histogram['Reflection Lists']
2032            for phase in refLists:
2033                if phase not in Phases:     #skips deleted or renamed phases silently!
2034                    continue
2035                Phase = Phases[phase]
2036                im = 0
2037                if Phase['General']['Type'] in ['modulated','magnetic']:
2038                    im = 1
2039                pId = Phase['pId']
2040                phfx = '%d:%d:'%(pId,hId)
2041                refDict = refLists[phase]
2042                sumFo = 0.0
2043                sumdF = 0.0
2044                sumFosq = 0.0
2045                sumdFsq = 0.0
2046                sumInt = 0.0
2047                for refl in refDict['RefList']:
2048                    if 'C' in calcControls[hfx+'histType']:
2049                        yp = np.zeros_like(yb)
2050                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2051                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2052                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2053                        iFin2 = iFin
2054                        if not iBeg+iFin:       #peak below low limit - skip peak
2055                            continue
2056                        elif not iBeg-iFin:     #peak above high limit - done
2057                            break
2058                        elif iBeg < iFin:
2059                            yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
2060                            sumInt += refl[11+im]*refl[9+im]
2061                            if Ka2:
2062                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2063                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2064                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
2065                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
2066                                if iFin2 > iBeg2: 
2067                                    yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
2068                                    sumInt += refl[11+im]*refl[9+im]*kRatio
2069                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
2070                               
2071                    elif 'T' in calcControls[hfx+'histType']:
2072                        yp = np.zeros_like(yb)
2073                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2074                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
2075                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
2076                        if iBeg < iFin:
2077                            yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))  #>90% of time spent here
2078                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
2079                            sumInt += refl[11+im]*refl[9+im]
2080                    Fo = np.sqrt(np.abs(refl[8+im]))
2081                    Fc = np.sqrt(np.abs(refl[9]+im))
2082                    sumFo += Fo
2083                    sumFosq += refl[8+im]**2
2084                    sumdF += np.abs(Fo-Fc)
2085                    sumdFsq += (refl[8+im]-refl[9+im])**2
2086                if sumFo:
2087                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2088                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2089                else:
2090                    Histogram['Residuals'][phfx+'Rf'] = 100.
2091                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
2092                Histogram['Residuals'][phfx+'sumInt'] = sumInt
2093                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])
2094                Histogram['Residuals']['hId'] = hId
2095        elif 'HKLF' in histogram[:4]:
2096            Histogram = Histograms[histogram]
2097            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
2098               
2099def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2100    'Needs a doc string'
2101   
2102    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
2103        U = parmDict[hfx+'U']
2104        V = parmDict[hfx+'V']
2105        W = parmDict[hfx+'W']
2106        X = parmDict[hfx+'X']
2107        Y = parmDict[hfx+'Y']
2108        tanPos = tand(refl[5+im]/2.0)
2109        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2110        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
2111        sig = max(0.001,sig)
2112        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
2113        gam = max(0.001,gam)
2114        return sig,gam
2115               
2116    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
2117        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
2118            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
2119        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
2120        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2121        sig += Ssig
2122        gam += Sgam
2123        return sig,gam
2124       
2125    def GetReflAlpBet(refl,im,hfx,parmDict):
2126        alp = parmDict[hfx+'alpha']/refl[4+im]
2127        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
2128        return alp,bet
2129       
2130    hId = Histogram['hId']
2131    hfx = ':%d:'%(hId)
2132    bakType = calcControls[hfx+'bakType']
2133    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2134    yc = np.zeros_like(yb)
2135    cw = np.diff(x)
2136    cw = np.append(cw,cw[-1])
2137       
2138    if 'C' in calcControls[hfx+'histType']:   
2139        shl = max(parmDict[hfx+'SH/L'],0.002)
2140        Ka2 = False
2141        if hfx+'Lam1' in parmDict.keys():
2142            wave = parmDict[hfx+'Lam1']
2143            Ka2 = True
2144            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2145            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2146        else:
2147            wave = parmDict[hfx+'Lam']
2148    for phase in Histogram['Reflection Lists']:
2149        refDict = Histogram['Reflection Lists'][phase]
2150        if phase not in Phases:     #skips deleted or renamed phases silently!
2151            continue
2152        Phase = Phases[phase]
2153        pId = Phase['pId']
2154        pfx = '%d::'%(pId)
2155        phfx = '%d:%d:'%(pId,hId)
2156        hfx = ':%d:'%(hId)
2157        SGData = Phase['General']['SGData']
2158        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2159        im = 0
2160        if Phase['General']['Type'] in ['modulated','magnetic']:
2161            SSGData = Phase['General']['SSGData']
2162            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2163            im = 1  #offset in SS reflection list
2164            #??
2165        Dij = GetDij(phfx,SGData,parmDict)
2166        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2167        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2168        if np.any(np.diag(G)<0.):
2169            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
2170        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2171        Vst = np.sqrt(nl.det(G))    #V*
2172        if not Phase['General'].get('doPawley'):
2173            time0 = time.time()
2174            if im:
2175                SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2176            else:
2177                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2178#            print 'sf calc time: %.3fs'%(time.time()-time0)
2179        time0 = time.time()
2180        badPeak = False
2181        for iref,refl in enumerate(refDict['RefList']):
2182            if 'C' in calcControls[hfx+'histType']:
2183                if im:
2184                    h,k,l,m = refl[:4]
2185                else:
2186                    h,k,l = refl[:3]
2187                Uniq = np.inner(refl[:3],SGMT)
2188                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2189                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
2190                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2191                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2192                refl[11+im] *= Vst*Lorenz
2193                 
2194                if Phase['General'].get('doPawley'):
2195                    try:
2196                        if im:
2197                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2198                        else:
2199                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2200                        refl[9+im] = parmDict[pInd]
2201                    except KeyError:
2202#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2203                        continue
2204                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2205                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2206                iFin = np.searchsorted(x,refl[5+im]+fmax)
2207                if not iBeg+iFin:       #peak below low limit - skip peak
2208                    continue
2209                elif not iBeg-iFin:     #peak above high limit - done
2210                    break
2211                elif iBeg > iFin:   #bad peak coeff - skip
2212                    badPeak = True
2213                    continue
2214                yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
2215                if Ka2:
2216                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2217                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
2218                    iBeg = np.searchsorted(x,pos2-fmin)
2219                    iFin = np.searchsorted(x,pos2+fmax)
2220                    if not iBeg+iFin:       #peak below low limit - skip peak
2221                        continue
2222                    elif not iBeg-iFin:     #peak above high limit - done
2223                        return yc,yb
2224                    elif iBeg > iFin:   #bad peak coeff - skip
2225                        continue
2226                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))        #and here
2227            elif 'T' in calcControls[hfx+'histType']:
2228                h,k,l = refl[:3]
2229                Uniq = np.inner(refl[:3],SGMT)
2230                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
2231                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
2232#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
2233                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
2234                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
2235                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2236                refl[11+im] *= Vst*Lorenz
2237                if Phase['General'].get('doPawley'):
2238                    try:
2239                        if im:
2240                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2241                        else:
2242                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2243                        refl[9+im] = parmDict[pInd]
2244                    except KeyError:
2245#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2246                        continue
2247                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2248                iBeg = np.searchsorted(x,refl[5+im]-fmin)
2249                iFin = np.searchsorted(x,refl[5+im]+fmax)
2250                if not iBeg+iFin:       #peak below low limit - skip peak
2251                    continue
2252                elif not iBeg-iFin:     #peak above high limit - done
2253                    break
2254                elif iBeg > iFin:   #bad peak coeff - skip
2255                    badPeak = True
2256                    continue
2257                yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
2258#        print 'profile calc time: %.3fs'%(time.time()-time0)
2259    if badPeak:
2260        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped' 
2261    return yc,yb
2262   
2263def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
2264    'Needs a doc string'
2265   
2266    def cellVaryDerv(pfx,SGData,dpdA): 
2267        if SGData['SGLaue'] in ['-1',]:
2268            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2269                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2270        elif SGData['SGLaue'] in ['2/m',]:
2271            if SGData['SGUniq'] == 'a':
2272                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2273            elif SGData['SGUniq'] == 'b':
2274                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2275            else:
2276                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2277        elif SGData['SGLaue'] in ['mmm',]:
2278            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2279        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2280            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2281        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2282            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2283        elif SGData['SGLaue'] in ['3R', '3mR']:
2284            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2285        elif SGData['SGLaue'] in ['m3m','m3']:
2286            return [[pfx+'A0',dpdA[0]]]
2287           
2288    # create a list of dependent variables and set up a dictionary to hold their derivatives
2289    dependentVars = G2mv.GetDependentVars()
2290    depDerivDict = {}
2291    for j in dependentVars:
2292        depDerivDict[j] = np.zeros(shape=(len(x)))
2293    #print 'dependent vars',dependentVars
2294    lenX = len(x)               
2295    hId = Histogram['hId']
2296    hfx = ':%d:'%(hId)
2297    bakType = calcControls[hfx+'bakType']
2298    dMdv = np.zeros(shape=(len(varylist),len(x)))
2299    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
2300    if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
2301        bBpos =varylist.index(hfx+'Back;0')
2302        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2303    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2304    for name in varylist:
2305        if 'Debye' in name:
2306            id = int(name.split(';')[-1])
2307            parm = name[:int(name.rindex(';'))]
2308            ip = names.index(parm)
2309            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2310    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2311    for name in varylist:
2312        if 'BkPk' in name:
2313            parm,id = name.split(';')
2314            id = int(id)
2315            if parm in names:
2316                ip = names.index(parm)
2317                dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2318    cw = np.diff(x)
2319    cw = np.append(cw,cw[-1])
2320    Ka2 = False #also for TOF!
2321    if 'C' in calcControls[hfx+'histType']:   
2322        shl = max(parmDict[hfx+'SH/L'],0.002)
2323        if hfx+'Lam1' in parmDict.keys():
2324            wave = parmDict[hfx+'Lam1']
2325            Ka2 = True
2326            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2327            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2328        else:
2329            wave = parmDict[hfx+'Lam']
2330    for phase in Histogram['Reflection Lists']:
2331        refDict = Histogram['Reflection Lists'][phase]
2332        if phase not in Phases:     #skips deleted or renamed phases silently!
2333            continue
2334        Phase = Phases[phase]
2335        SGData = Phase['General']['SGData']
2336        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2337        im = 0
2338        if Phase['General']['Type'] in ['modulated','magnetic']:
2339            SSGData = Phase['General']['SSGData']
2340            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2341            im = 1  #offset in SS reflection list
2342            #??
2343        pId = Phase['pId']
2344        pfx = '%d::'%(pId)
2345        phfx = '%d:%d:'%(pId,hId)
2346        Dij = GetDij(phfx,SGData,parmDict)
2347        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
2348        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2349        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2350        if not Phase['General'].get('doPawley'):
2351            time0 = time.time()
2352            if im:
2353                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2354            else:
2355                dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2356#            print 'sf-derv time %.3fs'%(time.time()-time0)
2357            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2358        time0 = time.time()
2359        for iref,refl in enumerate(refDict['RefList']):
2360            if im:
2361                h,k,l,m = refl[:4]
2362            else:
2363                h,k,l = refl[:3]
2364            Uniq = np.inner(refl[:3],SGMT)
2365            if 'T' in calcControls[hfx+'histType']:
2366                wave = refl[14+im]
2367            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2368            if 'C' in calcControls[hfx+'histType']:        #CW powder
2369                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
2370            else: #'T'OF
2371                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
2372            iBeg = np.searchsorted(x,refl[5+im]-fmin)
2373            iFin = np.searchsorted(x,refl[5+im]+fmax)
2374            if not iBeg+iFin:       #peak below low limit - skip peak
2375                continue
2376            elif not iBeg-iFin:     #peak above high limit - done
2377                break
2378            pos = refl[5+im]
2379            if 'C' in calcControls[hfx+'histType']:
2380                tanth = tand(pos/2.0)
2381                costh = cosd(pos/2.0)
2382                lenBF = iFin-iBeg
2383                dMdpk = np.zeros(shape=(6,lenBF))
2384                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
2385                for i in range(5):
2386                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
2387                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2388                if Ka2:
2389                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2390                    iBeg2 = np.searchsorted(x,pos2-fmin)
2391                    iFin2 = np.searchsorted(x,pos2+fmax)
2392                    if iBeg2-iFin2:
2393                        lenBF2 = iFin2-iBeg2
2394                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2395                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
2396                        for i in range(5):
2397                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
2398                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
2399                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2400            else:   #'T'OF
2401                lenBF = iFin-iBeg
2402                if lenBF < 0:   #bad peak coeff
2403                    break
2404                dMdpk = np.zeros(shape=(6,lenBF))
2405                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
2406                for i in range(6):
2407                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
2408                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
2409            if Phase['General'].get('doPawley'):
2410                dMdpw = np.zeros(len(x))
2411                try:
2412                    if im:
2413                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
2414                    else:
2415                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2416                    idx = varylist.index(pIdx)
2417                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
2418                    if Ka2: #not for TOF either
2419                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
2420                    dMdv[idx] = dMdpw
2421                except: # ValueError:
2422                    pass
2423            if 'C' in calcControls[hfx+'histType']:
2424                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
2425                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2426                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2427                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2428                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2429                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2430                    hfx+'DisplaceY':[dpdY,'pos'],}
2431                if 'Bragg' in calcControls[hfx+'instType']:
2432                    names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
2433                        hfx+'SurfRoughB':[dFdAb[1],'int'],})
2434                else:
2435                    names.update({hfx+'Absorption':[dFdAb,'int'],})
2436            else:   #'T'OF
2437                dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
2438                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
2439                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
2440                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
2441                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
2442                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
2443                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
2444                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
2445            for name in names:
2446                item = names[name]
2447                if name in varylist:
2448                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2449                    if Ka2 and iFin2-iBeg2:
2450                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2451                elif name in dependentVars:
2452                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2453                    if Ka2 and iFin2-iBeg2:
2454                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2455            for iPO in dIdPO:
2456                if iPO in varylist:
2457                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2458                    if Ka2 and iFin2-iBeg2:
2459                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2460                elif iPO in dependentVars:
2461                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2462                    if Ka2 and iFin2-iBeg2:
2463                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2464            for i,name in enumerate(['omega','chi','phi']):
2465                aname = pfx+'SH '+name
2466                if aname in varylist:
2467                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2468                    if Ka2 and iFin2-iBeg2:
2469                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2470                elif aname in dependentVars:
2471                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2472                    if Ka2 and iFin2-iBeg2:
2473                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2474            for iSH in dFdODF:
2475                if iSH in varylist:
2476                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2477                    if Ka2 and iFin2-iBeg2:
2478                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2479                elif iSH in dependentVars:
2480                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2481                    if Ka2 and iFin2-iBeg2:
2482                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2483            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2484            for name,dpdA in cellDervNames:
2485                if name in varylist:
2486                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2487                    if Ka2 and iFin2-iBeg2:
2488                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2489                elif name in dependentVars:
2490                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2491                    if Ka2 and iFin2-iBeg2:
2492                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2493            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
2494            for name in dDijDict:
2495                if name in varylist:
2496                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2497                    if Ka2 and iFin2-iBeg2:
2498                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2499                elif name in dependentVars:
2500                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2501                    if Ka2 and iFin2-iBeg2:
2502                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2503            for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
2504                if name in varylist:
2505                    dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
2506                    if Ka2 and iFin2-iBeg2:
2507                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2508                elif name in dependentVars:
2509                    depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
2510                    if Ka2 and iFin2-iBeg2:
2511                        depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
2512            if 'C' in calcControls[hfx+'histType']:
2513                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2514            else:   #'T'OF
2515                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
2516            for name in gamDict:
2517                if name in varylist:
2518                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2519                    if Ka2 and iFin2-iBeg2:
2520                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2521                elif name in dependentVars:
2522                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2523                    if Ka2 and iFin2-iBeg2:
2524                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2525            for name in sigDict:
2526                if name in varylist:
2527                    dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2528                    if Ka2 and iFin2-iBeg2:
2529                        dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2530                elif name in dependentVars:
2531                    depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2532                    if Ka2 and iFin2-iBeg2:
2533                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2534            for name in ['BabA','BabU']:
2535                if refl[9+im]:
2536                    if phfx+name in varylist:
2537                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
2538                        if Ka2 and iFin2-iBeg2:
2539                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]
2540                    elif phfx+name in dependentVars:                   
2541                        depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
2542                        if Ka2 and iFin2-iBeg2:
2543                            depDerivDict[phfx+name][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]                 
2544            if not Phase['General'].get('doPawley'):
2545                #do atom derivatives -  for RB,F,X & U so far             
2546                corr = dervDict['int']/refl[9+im]
2547                if Ka2 and iFin2-iBeg2:
2548                    corr2 = dervDict2['int']/refl[9+im]
2549                for name in varylist+dependentVars:
2550                    if '::RBV;' in name:
2551                        pass
2552                    else:
2553                        try:
2554                            aname = name.split(pfx)[1][:2]
2555                            if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
2556                        except IndexError:
2557                            continue
2558                    if name in varylist:
2559                        dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2560                        if Ka2 and iFin2-iBeg2:
2561                            dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2562                    elif name in dependentVars:
2563                        depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2564                        if Ka2 and iFin2-iBeg2:
2565                            depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2566    #        print 'profile derv time: %.3fs'%(time.time()-time0)
2567    # now process derivatives in constraints
2568    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2569    return dMdv
2570   
2571def UserRejectHKL(ref,im,userReject):
2572    if ref[5+im]/ref[6+im] < userReject['minF/sig']:
2573        return False
2574    elif userReject['MaxD'] < ref[4+im] > userReject['MinD']:
2575        return False
2576    elif ref[11+im] < userReject['MinExt']:
2577        return False
2578    elif abs(ref[5+im]-ref[7+im])/ref[6+im] > userReject['MaxDF/F']:
2579        return False
2580    return True
2581   
2582def dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict):
2583    '''Loop over reflections in a HKLF histogram and compute derivatives of the fitting
2584    model (M) with respect to all parameters.  Independent and dependant dM/dp arrays
2585    are returned to either dervRefine or HessRefine.
2586
2587    :returns:
2588    '''
2589    nobs = Histogram['Residuals']['Nobs']
2590    hId = Histogram['hId']
2591    hfx = ':%d:'%(hId)
2592    pfx = '%d::'%(Phase['pId'])
2593    phfx = '%d:%d:'%(Phase['pId'],hId)
2594    SGData = Phase['General']['SGData']
2595    im = 0
2596    if Phase['General']['Type'] in ['modulated','magnetic']:
2597        SSGData = Phase['General']['SSGData']
2598        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2599        im = 1  #offset in SS reflection list
2600    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2601    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2602    refDict = Histogram['Data']
2603    if parmDict[phfx+'Scale'] < 0.:
2604        parmDict[phfx+'Scale'] = .001
2605    if im:
2606        dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2607    else:
2608        dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2609    ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
2610    dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
2611    dependentVars = G2mv.GetDependentVars()
2612    depDerivDict = {}
2613    for j in dependentVars:
2614        depDerivDict[j] = np.zeros(shape=(len(refDict['RefList'])))
2615    wdf = np.zeros(len(refDict['RefList']))
2616    if calcControls['F**2']:
2617        for iref,ref in enumerate(refDict['RefList']):
2618            if ref[6+im] > 0:
2619                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2620                w = 1.0/ref[6+im]
2621                if ref[3+im] > 0:
2622                    wdf[iref] = w*(ref[5+im]-ref[7+im])
2623                    for j,var in enumerate(varylist):
2624                        if var in dFdvDict:
2625                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2626                    for var in dependentVars:
2627                        if var in dFdvDict:
2628                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2629                    if phfx+'Scale' in varylist:
2630                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[7+im]*ref[11+im]/parmDict[phfx+'Scale']  #OK
2631                    elif phfx+'Scale' in dependentVars:
2632                        depDerivDict[phfx+'Scale'][iref] = w*ref[7+im]*ref[11+im]/parmDict[phfx+'Scale']   #OK
2633                    for item in ['Ep','Es','Eg']:
2634                        if phfx+item in varylist and phfx+item in dervDict:
2635                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2636                        elif phfx+item in dependentVars and phfx+item in dervDict:
2637                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
2638                    for item in ['BabA','BabU']:
2639                        if phfx+item in varylist:
2640                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2641                        elif phfx+item in dependentVars:
2642                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2643    else:   #F refinement
2644        for iref,ref in enumerate(refDict['RefList']):
2645            if ref[5+im] > 0.:
2646                dervDict,dervCor = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1:]
2647                Fo = np.sqrt(ref[5+im])
2648                Fc = np.sqrt(ref[7+im])
2649                w = 1.0/ref[6+im]
2650                if ref[3+im] > 0:
2651                    wdf[iref] = 2.0*Fc*w*(Fo-Fc)
2652                    for j,var in enumerate(varylist):
2653                        if var in dFdvDict:
2654                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2655                    for var in dependentVars:
2656                        if var in dFdvDict:
2657                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
2658                    if phfx+'Scale' in varylist:
2659                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[7+im]*ref[11+im]/parmDict[phfx+'Scale']  #OK
2660                    elif phfx+'Scale' in dependentVars:
2661                        depDerivDict[phfx+'Scale'][iref] = w*ref[7+im]*ref[11+im]/parmDict[phfx+'Scale']   #OK                   
2662                    for item in ['Ep','Es','Eg']:   #OK!
2663                        if phfx+item in varylist and phfx+item in dervDict:
2664                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] 
2665                        elif phfx+item in dependentVars and phfx+item in dervDict:
2666                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
2667                    for item in ['BabA','BabU']:
2668                        if phfx+item in varylist:
2669                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2670                        elif phfx+item in dependentVars:
2671                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
2672    return dMdvh,depDerivDict,wdf
2673   
2674
2675def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2676    '''Loop over histograms and compute derivatives of the fitting
2677    model (M) with respect to all parameters.  Results are returned in
2678    a Jacobian matrix (aka design matrix) of dimensions (n by m) where
2679    n is the number of parameters and m is the number of data
2680    points. This can exceed memory when m gets large. This routine is
2681    used when refinement derivatives are selected as "analtytic
2682    Jacobian" in Controls.
2683
2684    :returns: Jacobian numpy.array dMdv for all histograms concatinated
2685    '''
2686    parmDict.update(zip(varylist,values))
2687    G2mv.Dict2Map(parmDict,varylist)
2688    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2689    nvar = len(varylist)
2690    dMdv = np.empty(0)
2691    histoList = Histograms.keys()
2692    histoList.sort()
2693    for histogram in histoList:
2694        if 'PWDR' in histogram[:4]:
2695            Histogram = Histograms[histogram]
2696            hId = Histogram['hId']
2697            hfx = ':%d:'%(hId)
2698            wtFactor = calcControls[hfx+'wtFactor']
2699            Limits = calcControls[hfx+'Limits']
2700            x,y,w,yc,yb,yd = Histogram['Data']
2701            xB = np.searchsorted(x,Limits[0])
2702            xF = np.searchsorted(x,Limits[1])
2703            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
2704                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2705        elif 'HKLF' in histogram[:4]:
2706            Histogram = Histograms[histogram]
2707            phase = Histogram['Reflection Lists']
2708            Phase = Phases[phase]
2709            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2710            hfx = ':%d:'%(Histogram['hId'])
2711            wtFactor = calcControls[hfx+'wtFactor']
2712            # now process derivatives in constraints
2713            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2714        else:
2715            continue        #skip non-histogram entries
2716        if len(dMdv):
2717            dMdv = np.concatenate((dMdv.T,np.sqrt(wtFactor)*dMdvh.T)).T
2718        else:
2719            dMdv = np.sqrt(wtFactor)*dMdvh
2720           
2721    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2722    if np.any(pVals):
2723        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2724        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
2725       
2726    return dMdv
2727
2728def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
2729    '''Loop over histograms and compute derivatives of the fitting
2730    model (M) with respect to all parameters.  For each histogram, the
2731    Jacobian matrix, dMdv, with dimensions (n by m) where n is the
2732    number of parameters and m is the number of data points *in the
2733    histogram*. The (n by n) Hessian is computed from each Jacobian
2734    and it is returned.  This routine is used when refinement
2735    derivatives are selected as "analtytic Hessian" in Controls.
2736
2737    :returns: Vec,Hess where Vec is the least-squares vector and Hess is the Hessian
2738    '''
2739    parmDict.update(zip(varylist,values))
2740    G2mv.Dict2Map(parmDict,varylist)
2741    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2742    #fixup H atom positions here?
2743    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
2744    nvar = len(varylist)
2745    Hess = np.empty(0)
2746    histoList = Histograms.keys()
2747    histoList.sort()
2748    for histogram in histoList:
2749        if 'PWDR' in histogram[:4]:
2750            Histogram = Histograms[histogram]
2751            hId = Histogram['hId']
2752            hfx = ':%d:'%(hId)
2753            wtFactor = calcControls[hfx+'wtFactor']
2754            Limits = calcControls[hfx+'Limits']
2755            x,y,w,yc,yb,yd = Histogram['Data']
2756            W = wtFactor*w
2757            dy = y-yc
2758            xB = np.searchsorted(x,Limits[0])
2759            xF = np.searchsorted(x,Limits[1])
2760            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
2761                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
2762            Wt = ma.sqrt(W[xB:xF])[np.newaxis,:]
2763            Dy = dy[xB:xF][np.newaxis,:]
2764            dMdvh *= Wt
2765            if dlg:
2766                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))
2767            if len(Hess):
2768                Hess += np.inner(dMdvh,dMdvh)
2769                dMdvh *= Wt*Dy
2770                Vec += np.sum(dMdvh,axis=1)
2771            else:
2772                Hess = np.inner(dMdvh,dMdvh)
2773                dMdvh *= Wt*Dy
2774                Vec = np.sum(dMdvh,axis=1)
2775        elif 'HKLF' in histogram[:4]:
2776            Histogram = Histograms[histogram]
2777            phase = Histogram['Reflection Lists']
2778            Phase = Phases[phase]
2779            dMdvh,depDerivDict,wdf = dervHKLF(Histogram,Phase,calcControls,varylist,parmDict,rigidbodyDict)
2780            hId = Histogram['hId']
2781            hfx = ':%d:'%(Histogram['hId'])
2782            wtFactor = calcControls[hfx+'wtFactor']
2783            # now process derivatives in constraints
2784            G2mv.Dict2Deriv(varylist,depDerivDict,dMdvh)
2785#            print 'matrix build time: %.3f'%(time.time()-time0)
2786
2787            if dlg:
2788                dlg.Update(Histogram['Residuals']['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2789            if len(Hess):
2790                Vec += wtFactor*np.sum(dMdvh*wdf,axis=1)
2791                Hess += wtFactor*np.inner(dMdvh,dMdvh)
2792            else:
2793                Vec = wtFactor*np.sum(dMdvh*wdf,axis=1)
2794                Hess = wtFactor*np.inner(dMdvh,dMdvh)
2795        else:
2796            continue        #skip non-histogram entries
2797    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
2798    if np.any(pVals):
2799        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
2800        Vec += np.sum(dpdv*pWt*pVals,axis=1)
2801        Hess += np.inner(dpdv*pWt,dpdv)
2802    return Vec,Hess
2803
2804def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None):       
2805    '''Computes the point-by-point discrepancies between every data point in every histogram
2806    and the observed value
2807   
2808    :returns: an np array of differences between observed and computed diffraction values.
2809    '''
2810    Values2Dict(parmDict, varylist, values)
2811    G2mv.Dict2Map(parmDict,varylist)
2812    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2813    M = np.empty(0)
2814    SumwYo = 0
2815    Nobs = 0
2816    Nrej = 0
2817    Next = 0
2818    ApplyRBModels(parmDict,Phases,rigidbodyDict)
2819    #fixup Hatom positions here....
2820    histoList = Histograms.keys()
2821    histoList.sort()
2822    for histogram in histoList:
2823        if 'PWDR' in histogram[:4]:
2824            Histogram = Histograms[histogram]
2825            hId = Histogram['hId']
2826            hfx = ':%d:'%(hId)
2827            wtFactor = calcControls[hfx+'wtFactor']
2828            Limits = calcControls[hfx+'Limits']
2829            x,y,w,yc,yb,yd = Histogram['Data']
2830            yc *= 0.0                           #zero full calcd profiles
2831            yb *= 0.0
2832            yd *= 0.0
2833            xB = np.searchsorted(x,Limits[0])
2834            xF = np.searchsorted(x,Limits[1])
2835            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
2836                varylist,Histogram,Phases,calcControls,pawleyLookup)
2837            yc[xB:xF] += yb[xB:xF]
2838            if not np.any(y):                   #fill dummy data
2839                rv = st.poisson(yc[xB:xF])
2840                y[xB:xF] = rv.rvs()
2841                Z = np.ones_like(yc[xB:xF])
2842                Z[1::2] *= -1
2843                y[xB:xF] = yc[xB:xF]+np.abs(y[xB:xF]-yc[xB:xF])*Z
2844                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
2845            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2846            W = wtFactor*w
2847            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
2848            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
2849            Nobs += Histogram['Residuals']['Nobs']
2850            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
2851            SumwYo += Histogram['Residuals']['sumwYo']
2852            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
2853            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
2854            sumYmB = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],ma.abs(y[xB:xF]-yb[xB:xF]),0.))
2855            sumwYmB2 = ma.sum(ma.where(yc[xB:xF]!=yb[xB:xF],W[xB:xF]*(y[xB:xF]-yb[xB:xF])**2,0.))
2856            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.))
2857            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.))
2858            Histogram['Residuals']['Rb'] = min(100.,100.*sumYB/sumYmB)
2859            Histogram['Residuals']['wRb'] = min(100.,100.*ma.sqrt(sumwYB2/sumwYmB2))
2860            Histogram['Residuals']['wRmin'] = min(100.,100.*ma.sqrt(Histogram['Residuals']['Nobs']/Histogram['Residuals']['sumwYo']))
2861            if dlg:
2862                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2863            M = np.concatenate((M,wdy))
2864#end of PWDR processing
2865        elif 'HKLF' in histogram[:4]:
2866            Histogram = Histograms[histogram]
2867            Histogram['Residuals'] = {}
2868            phase = Histogram['Reflection Lists']
2869            Phase = Phases[phase]
2870            hId = Histogram['hId']
2871            hfx = ':%d:'%(hId)
2872            wtFactor = calcControls[hfx+'wtFactor']
2873            pfx = '%d::'%(Phase['pId'])
2874            phfx = '%d:%d:'%(Phase['pId'],hId)
2875            SGData = Phase['General']['SGData']
2876            im = 0
2877            if parmDict[phfx+'Scale'] < 0.:
2878                parmDict[phfx+'Scale'] = .001               
2879            if Phase['General']['Type'] in ['modulated','magnetic']:
2880                SSGData = Phase['General']['SSGData']
2881                SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2882                im = 1  #offset in SS reflection list
2883            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2884            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2885            refDict = Histogram['Data']
2886            time0 = time.time()
2887            if im:
2888                SStructureFactor2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2889            else:
2890                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
2891#            print 'sf-calc time: %.3f'%(time.time()-time0)
2892            df = np.zeros(len(refDict['RefList']))
2893            sumwYo = 0
2894            sumFo = 0
2895            sumFo2 = 0
2896            sumdF = 0
2897            sumdF2 = 0
2898            if im:
2899                sumSSFo = np.zeros(10)
2900                sumSSFo2 = np.zeros(10)
2901                sumSSdF = np.zeros(10)
2902                sumSSdF2 = np.zeros(10)
2903                sumSSwYo = np.zeros(10)
2904                sumSSwdf2 = np.zeros(10)
2905                SSnobs = np.zeros(10)
2906            nobs = 0
2907            nrej = 0
2908            next = 0
2909            maxH = 0
2910            if calcControls['F**2']:
2911                for i,ref in enumerate(refDict['RefList']):
2912                    if ref[6+im] > 0:
2913                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2914                        w = 1.0/ref[6+im]   # 1/sig(F^2)
2915                        ref[7+im] *= parmDict[phfx+'Scale']*ref[11+im]  #correct Fc^2 for extinction
2916                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2917                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
2918                            ref[3+im] = abs(ref[3+im])      #mark as allowed
2919                            Fo = np.sqrt(ref[5+im])
2920                            sumFo += Fo
2921                            sumFo2 += ref[5+im]
2922                            sumdF += abs(Fo-np.sqrt(ref[7+im]))
2923                            sumdF2 += abs(ref[5+im]-ref[7+im])
2924                            nobs += 1
2925                            df[i] = -w*(ref[5+im]-ref[7+im])
2926                            sumwYo += (w*ref[5+im])**2      #w*Fo^2
2927                            if im:  #accumulate super lattice sums
2928                                ind = int(abs(ref[3]))
2929                                sumSSFo[ind] += Fo
2930                                sumSSFo2[ind] += ref[5+im]
2931                                sumSSdF[ind] += abs(Fo-np.sqrt(ref[7+im]))
2932                                sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
2933                                sumSSwYo[ind] += (w*ref[5+im])**2      #w*Fo^2
2934                                sumSSwdf2[ind] +=  df[i]**2
2935                                SSnobs[ind] += 1
2936                                maxH = max(maxH,ind)                           
2937                        else:
2938                            if ref[3+im]:
2939                                ref[3+im] = -abs(ref[3+im])      #mark as rejected
2940                                nrej += 1
2941                            else:   #sp.gp.extinct
2942                                next += 1
2943            else:
2944                for i,ref in enumerate(refDict['RefList']):
2945                    if ref[5+im] > 0.:
2946                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
2947                        ref[7+im] *= parmDict[phfx+'Scale']*ref[11+im]    #correct Fc^2 for extinction
2948                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
2949                        Fo = np.sqrt(ref[5+im])
2950                        Fc = np.sqrt(ref[7+im])
2951                        w = 2.0*Fo/ref[6+im]    # 1/sig(F)?
2952                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
2953                            ref[3+im] = abs(ref[3+im])      #mark as allowed
2954                            sumFo += Fo
2955                            sumFo2 += ref[5+im]
2956                            sumdF += abs(Fo-Fc)
2957                            sumdF2 += abs(ref[5+im]-ref[7+im])
2958                            nobs += 1
2959                            df[i] = -w*(Fo-Fc)
2960                            sumwYo += (w*Fo)**2
2961                            if im:
2962                                ind = int(abs(ref[3]))
2963                                sumSSFo[ind] += Fo
2964                                sumSSFo2[ind] += ref[5+im]
2965                                sumSSdF[ind] += abs(Fo-Fc)
2966                                sumSSdF2[ind] += abs(ref[5+im]-ref[7+im])
2967                                sumSSwYo[ind] += (w*Fo)**2                                                           
2968                                sumSSwdf2[ind] +=  df[i]**2
2969                                SSnobs[ind] += 1                           
2970                                maxH = max(maxH,ind)                           
2971                        else:
2972                            if ref[3+im]:
2973                                ref[3+im] = -abs(ref[3+im])      #mark as rejected
2974                                nrej += 1
2975                            else:   #sp.gp.extinct
2976                                next += 1
2977            Histogram['Residuals']['Nobs'] = nobs
2978            Histogram['Residuals']['sumwYo'] = sumwYo
2979            SumwYo += sumwYo
2980            Histogram['Residuals']['wR'] = min(100.,np.sqrt(np.sum(df**2)/sumwYo)*100.)
2981            Histogram['Residuals'][phfx+'Rf'] = 100.*sumdF/sumFo
2982            Histogram['Residuals'][phfx+'Rf^2'] = 100.*sumdF2/sumFo2
2983            Histogram['Residuals'][phfx+'Nref'] = nobs
2984            Histogram['Residuals'][phfx+'Nrej'] = nrej
2985            Histogram['Residuals'][phfx+'Next'] = next
2986            if im:
2987                Histogram['Residuals'][phfx+'SSRf'] = 100.*sumSSdF[:maxH+1]/sumSSFo[:maxH+1]
2988                Histogram['Residuals'][phfx+'SSRf^2'] = 100.*sumSSdF2[:maxH+1]/sumSSFo2[:maxH+1]
2989                Histogram['Residuals'][phfx+'SSNref'] = SSnobs[:maxH+1]
2990                Histogram['Residuals']['SSwR'] = np.sqrt(sumSSwdf2[:maxH+1]/sumSSwYo[:maxH+1])*100.               
2991            Nobs += nobs
2992            Nrej += nrej
2993            Next += next
2994            if dlg:
2995                dlg.Update(Histogram['Residuals']['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['Residuals']['wR'],'%'))[0]
2996            M = np.concatenate((M,wtFactor*df))
2997# end of HKLF processing
2998    Histograms['sumwYo'] = SumwYo
2999    Histograms['Nobs'] = Nobs
3000    Histograms['Nrej'] = Nrej
3001    Histograms['Next'] = Next
3002    Rw = min(100.,np.sqrt(np.sum(M**2)/SumwYo)*100.)
3003    if dlg:
3004        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
3005        if not GoOn:
3006            parmDict['saved values'] = values
3007            dlg.Destroy()
3008            raise G2obj.G2Exception('User abort')         #Abort!!
3009    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
3010    if len(pVals):
3011        pSum = np.sum(pWt*pVals**2)
3012        for name in pWsum:
3013            if pWsum[name]:
3014                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
3015        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
3016        Nobs += len(pVals)
3017        M = np.concatenate((M,np.sqrt(pWt)*pVals))
3018    return M
3019                       
Note: See TracBrowser for help on using the repository browser.