source: trunk/GSASIIstrMath.py @ 2401

Last change on this file since 2401 was 2401, checked in by vondreele, 5 years ago

begin adding magnetism; change space group fortran code
modify tutorial menus

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