source: trunk/GSASIIstrMath.py @ 2480

Last change on this file since 2480 was 2480, checked in by vondreele, 7 years ago

add make magnetic phase to General/Transform? option
trap missing rigid bodies to define torsion seq option
fix lighting issues for polygons
fix xye importer to stop on trailing blank lines rather than crashing
fix error in powder structure factor calc
add magnetic structure factor calc. (some error still & no derivatives yet)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 223.1 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2016-09-28 17:16:00 +0000 (Wed, 28 Sep 2016) $
8# $Author: vondreele $
9# $Revision: 2480 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 2480 2016-09-28 17:16:00Z 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: 2480 $")
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    Gdata = np.zeros((3,Natoms))
596    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
597        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
598        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
599        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
600        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5],
601        'AMx:':Gdata[0],'AMy:':Gdata[1],'AMz:':Gdata[2],}
602    for iatm in range(Natoms):
603        for key in keys:
604            parm = pfx+key+str(iatm)
605            if parm in parmDict:
606                keys[key][iatm] = parmDict[parm]
607    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.
608    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata
609   
610def GetAtomSSFXU(pfx,calcControls,parmDict):
611    'Needs a doc string'
612    Natoms = calcControls['Natoms'][pfx]
613    maxSSwave = calcControls['maxSSwave'][pfx]
614    Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
615        'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
616    XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
617    FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
618    USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
619    MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
620    waveTypes = []
621    keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
622        'Tmin:':XSSdata[0],'Tmax:':XSSdata[1],'Xmax:':XSSdata[2],'Ymax:':XSSdata[3],'Zmax:':XSSdata[4],
623        'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
624        'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
625        'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
626        'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
627    for iatm in range(Natoms):
628        waveTypes.append(parmDict[pfx+'waveType:'+str(iatm)])
629        for key in keys:
630            for m in range(Nwave[key[0]]):
631                parm = pfx+key+str(iatm)+':%d'%(m)
632                if parm in parmDict:
633                    keys[key][m][iatm] = parmDict[parm]
634    return np.array(waveTypes),FSSdata,XSSdata,USSdata,MSSdata
635   
636#def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):
637#   
638#    Natoms = calcControls['Natoms'][pfx]
639#    maxSSwave = calcControls['maxSSwave'][pfx]
640#    Smult = np.zeros((Natoms,len(SGOps)))
641#    TauT = np.zeros((Natoms,len(SGOps)))
642#    for ix,xyz in enumerate(XData.T):
643#        for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):
644#            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)
645#            Smult[ix][isym] = sdet
646#            TauT[ix][isym] = tauT
647#    return Smult,TauT
648#   
649def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
650    ''' Compute structure factors for all h,k,l for phase
651    puts the result, F^2, in each ref[8] in refList
652    operates on blocks of 100 reflections for speed
653    input:
654   
655    :param dict refDict: where
656        'RefList' list where each ref = h,k,l,it,d,...
657        'FF' dict of form factors - filed in below
658    :param np.array G:      reciprocal metric tensor
659    :param str pfx:    phase id string
660    :param dict SGData: space group info. dictionary output from SpcGroup
661    :param dict calcControls:
662    :param dict ParmDict:
663
664    '''       
665    phfx = pfx.split(':')[0]+hfx
666    ast = np.sqrt(np.diag(G))
667    Mast = twopisq*np.multiply.outer(ast,ast)
668    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
669    SGT = np.array([ops[1] for ops in SGData['SGOps']])
670    Ncen = len(SGData['SGCen'])
671    FFtables = calcControls['FFtables']
672    BLtables = calcControls['BLtables']
673    MFtables = calcControls['MFtables']
674    Amat,Bmat = G2lat.Gmat2AB(G)
675    Flack = 1.0
676    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
677        Flack = 1.-2.*parmDict[phfx+'Flack']
678    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
679    TwDict = refDict.get('TwDict',{})           
680    if 'S' in calcControls[hfx+'histType']:
681        NTL = calcControls[phfx+'NTL']
682        NM = calcControls[phfx+'TwinNMN']+1
683        TwinLaw = calcControls[phfx+'TwinLaw']
684        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
685        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
686    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
687        GetAtomFXU(pfx,calcControls,parmDict)
688    if parmDict[pfx+'isMag']:
689        Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
690        Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
691        Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
692        Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
693        if SGData['SGInv']:
694            Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
695        Gdata = np.repeat(Gdata,Ncen,axis=1)    #dup over cell centering
696        Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
697        Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps*Inv*Ncen
698        Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last
699#        GSASIIpath.IPyBreak()
700        Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
701        if SGData['SGInv']:
702            Mag = np.repeat(Mag,2,axis=0)                  #Mag same length as Gdata
703    if 'NC' in calcControls[hfx+'histType']:
704        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
705    elif 'X' in calcControls[hfx+'histType']:
706        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
707        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
708    Uij = np.array(G2lat.U6toUij(Uijdata))
709    bij = Mast*Uij.T
710    blkSize = 100       #no. of reflections in a block - size seems optimal
711    nRef = refDict['RefList'].shape[0]
712    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
713        SQ = 1./(2.*refDict['RefList'].T[4])**2
714        if 'N' in calcControls[hfx+'histType']:
715            dat = G2el.getBLvalues(BLtables)
716            refDict['FF']['El'] = dat.keys()
717            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
718            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
719            for iel,El in enumerate(refDict['FF']['El']):
720                if El in MFtables:
721                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
722        else:       #'X'
723            dat = G2el.getFFvalues(FFtables,0.)
724            refDict['FF']['El'] = dat.keys()
725            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
726            for iel,El in enumerate(refDict['FF']['El']):
727                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
728#reflection processing begins here - big arrays!
729    iBeg = 0
730    time0 = time.time()
731    while iBeg < nRef:
732        iFin = min(iBeg+blkSize,nRef)
733        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
734        H = refl.T[:3]                          #array(blkSize,3)
735        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
736        TwMask = np.any(H,axis=-1)
737        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
738            for ir in range(blkSize):
739                iref = ir+iBeg
740                if iref in TwDict:
741                    for i in TwDict[iref]:
742                        for n in range(NTL):
743                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
744            TwMask = np.any(H,axis=-1)
745        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
746        SQfactor = 4.0*SQ*twopisq               #ditto prev.
747        if 'T' in calcControls[hfx+'histType']:
748            if 'P' in calcControls[hfx+'histType']:
749                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
750            else:
751                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
752            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
753            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
754        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
755        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
756        Uniq = np.inner(H,SGMT)
757        Phi = np.inner(H,SGT)
758        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
759        sinp = np.sin(phase)
760        cosp = np.cos(phase)
761        biso = -SQfactor*Uisodata[:,nxs]
762        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
763        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
764        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
765        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
766        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
767        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
768            MF = np.repeat(refDict['FF']['MF'][iBeg:iFin].T[Tindx].T,len(TwinLaw),axis=0)   #Nref,Natm
769            TMcorr = 0.539*Tcorr[:,0,:]*MF                                                  #Nref,Natm
770            if SGData['SGInv']:
771                mphase = np.hstack((phase,-phase))
772            else:
773                mphase = phase
774            mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
775            mphase = np.concatenate(mphase,axis=1)              #Nref,Nop,Natm
776            sinm = np.sin(mphase)                               #ditto - match magstrfc.for
777            cosm = np.cos(mphase)                               #ditto
778            HM = np.inner(Bmat.T,H)                             #put into cartesian space
779            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
780            eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
781            Q = -HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]+Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
782            fam = Q*TMcorr[nxs,:,nxs,:]*SGData['MagMom'][nxs,nxs,:,nxs]*cosm[nxs,:,:,:]    #ditto
783            fbm = Q*TMcorr[nxs,:,nxs,:]*SGData['MagMom'][nxs,nxs,:,nxs]*sinm[nxs,:,:,:]    #ditto
784            fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
785            fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
786            GSASIIpath.IPyBreak()
787        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
788            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
789            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
790        else:
791            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
792            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
793        fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
794        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
795        if SGData['SGInv']: #centrosymmetric; B=0
796            fbs[0] *= 0.
797            fas[1] *= 0.
798        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
799            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
800            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
801            if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
802                refl.T[9] += np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
803        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
804            if len(TwinLaw) > 1:
805                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
806                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
807                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
808                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
809            else:   # checked correct!!
810                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
811                refl.T[7] = np.copy(refl.T[9])               
812                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
813#        GSASIIpath.IPyBreak()
814#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
815        iBeg += blkSize
816#    print ' %d sf time %.4f\r'%(nRef,time.time()-time0)
817   
818#def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
819#    '''Compute structure factor derivatives on single reflections - keep as it works for twins
820#    but is slower for powders/nontwins
821#    input:
822#   
823#    :param dict refDict: where
824#        'RefList' list where each ref = h,k,l,it,d,...
825#        'FF' dict of form factors - filled in below
826#    :param np.array G:      reciprocal metric tensor
827#    :param str hfx:    histogram id string
828#    :param str pfx:    phase id string
829#    :param dict SGData: space group info. dictionary output from SpcGroup
830#    :param dict calcControls:
831#    :param dict ParmDict:
832#   
833#    :returns: dict dFdvDict: dictionary of derivatives
834#    '''
835#    phfx = pfx.split(':')[0]+hfx
836#    ast = np.sqrt(np.diag(G))
837#    Mast = twopisq*np.multiply.outer(ast,ast)
838#    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
839#    SGT = np.array([ops[1] for ops in SGData['SGOps']])
840#    FFtables = calcControls['FFtables']
841#    BLtables = calcControls['BLtables']
842#    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
843#    TwDict = refDict.get('TwDict',{})           
844#    if 'S' in calcControls[hfx+'histType']:
845#        NTL = calcControls[phfx+'NTL']
846#        NM = calcControls[phfx+'TwinNMN']+1
847#        TwinLaw = calcControls[phfx+'TwinLaw']
848#        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
849#        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
850#    nTwin = len(TwinLaw)       
851#    nRef = len(refDict['RefList'])
852#    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
853#        GetAtomFXU(pfx,calcControls,parmDict)
854#    mSize = len(Mdata)
855#    FF = np.zeros(len(Tdata))
856#    if 'NC' in calcControls[hfx+'histType']:
857#        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
858#    elif 'X' in calcControls[hfx+'histType']:
859#        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
860#        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
861#    Uij = np.array(G2lat.U6toUij(Uijdata))
862#    bij = Mast*Uij.T
863#    dFdvDict = {}
864#    if nTwin > 1:
865#        dFdfr = np.zeros((nRef,nTwin,mSize))
866#        dFdx = np.zeros((nRef,nTwin,mSize,3))
867#        dFdui = np.zeros((nRef,nTwin,mSize))
868#        dFdua = np.zeros((nRef,nTwin,mSize,6))
869#        dFdbab = np.zeros((nRef,nTwin,2))
870#        dFdfl = np.zeros((nRef,nTwin))
871#        dFdtw = np.zeros((nRef,nTwin))
872#    else:
873#        dFdfr = np.zeros((nRef,mSize))
874#        dFdx = np.zeros((nRef,mSize,3))
875#        dFdui = np.zeros((nRef,mSize))
876#        dFdua = np.zeros((nRef,mSize,6))
877#        dFdbab = np.zeros((nRef,2))
878#        dFdfl = np.zeros((nRef))
879#        dFdtw = np.zeros((nRef))
880#    Flack = 1.0
881#    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
882#        Flack = 1.-2.*parmDict[phfx+'Flack']
883#    time0 = time.time()
884#    nref = len(refDict['RefList'])/100   
885#    for iref,refl in enumerate(refDict['RefList']):
886#        if 'T' in calcControls[hfx+'histType']:
887#            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
888#        H = np.array(refl[:3])
889#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
890#        TwMask = np.any(H,axis=-1)
891#        if TwinLaw.shape[0] > 1 and TwDict:
892#            if iref in TwDict:
893#                for i in TwDict[iref]:
894#                    for n in range(NTL):
895#                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
896#            TwMask = np.any(H,axis=-1)
897#        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
898#        SQfactor = 8.0*SQ*np.pi**2
899#        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
900#        Bab = parmDict[phfx+'BabA']*dBabdA
901#        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
902#        FF = refDict['FF']['FF'][iref].T[Tindx].T
903#        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
904#        Phi = np.inner(H,SGT)
905#        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
906#        sinp = np.sin(phase)
907#        cosp = np.cos(phase)
908#        occ = Mdata*Fdata/len(SGT)
909#        biso = -SQfactor*Uisodata[:,nxs]
910#        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
911#        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
912#        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
913#        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
914#        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
915#        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
916#        fot = (FF+FP-Bab)*Tcorr
917#        fotp = FPP*Tcorr       
918#        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
919#        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
920##        GSASIIpath.IPyBreak()
921#        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
922#        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
923#        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
924#        fbx = np.array([fot*cosp,-fotp*sinp])
925#        #sum below is over Uniq
926#        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
927#        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
928#        dfadui = np.sum(-SQfactor*fa,axis=-2)
929#        if nTwin > 1:
930#            dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
931#            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
932#            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
933#        else:
934#            dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)
935#            dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,nxs],axis=-2)
936#            # array(2,nAtom,3) & array(2,nAtom,6)
937#        if not SGData['SGInv']:
938#            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
939#            dfadba /= 2.
940#            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
941#            dfbdui = np.sum(-SQfactor*fb,axis=-2)
942#            if len(TwinLaw) > 1:
943#                dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
944#                dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
945#            else:
946#                dfadfl = np.sum(-FPP*Tcorr*sinp)
947#                dfbdfl = np.sum(FPP*Tcorr*cosp)
948#                dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=2)           
949#                dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,nxs],axis=2)
950#        else:
951#            dfbdfr = np.zeros_like(dfadfr)
952#            dfbdx = np.zeros_like(dfadx)
953#            dfbdui = np.zeros_like(dfadui)
954#            dfbdua = np.zeros_like(dfadua)
955#            dfbdba = np.zeros_like(dfadba)
956#            dfadfl = 0.0
957#            dfbdfl = 0.0
958#        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
959#        SA = fas[0]+fas[1]
960#        SB = fbs[0]+fbs[1]
961#        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
962#            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
963#                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
964#            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
965#                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
966#            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
967#                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
968#            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
969#                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
970#        else:
971#            if nTwin > 1:
972#                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)]
973#                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)]
974#                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)]
975#                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)]
976#                dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
977#            else:   #these are good for no twin single crystals
978#                dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)
979#                dFdx[iref] = 2.*SA*(dfadx[0]+dfadx[1])+2.*SB*(dfbdx[0]+dfbdx[1])
980#                dFdui[iref] = 2.*SA*(dfadui[0]+dfadui[1])+2.*SB*(dfbdui[0]+dfbdui[1])
981#                dFdua[iref] = 2.*SA*(dfadua[0]+dfadua[1])+2.*SB*(dfbdua[0]+dfbdua[1])
982#                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
983#        dFdbab[iref] = fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
984#            fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
985##        GSASIIpath.IPyBreak()
986#        if not iref%100 :
987#            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
988#    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
989#        #loop over atoms - each dict entry is list of derivatives for all the reflections
990#    if nTwin > 1:
991#        for i in range(len(Mdata)):     #these all OK?
992#            dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
993#            dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
994#            dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
995#            dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
996#            dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
997#            dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
998#            dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
999#            dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1000#            dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1001#            dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1002#            dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1003#    else:
1004#        for i in range(len(Mdata)):
1005#            dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1006#            dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1007#            dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1008#            dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1009#            dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1010#            dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1011#            dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1012#            dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1013#            dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1014#            dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1015#            dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1016#        dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
1017#    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1018#    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1019#    if nTwin > 1:
1020#        for i in range(nTwin):
1021#            dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1022#    return dFdvDict
1023   
1024def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1025    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
1026    faster than StructureFactorDerv - correct for powders/nontwins!!
1027    input:
1028   
1029    :param dict refDict: where
1030        'RefList' list where each ref = h,k,l,it,d,...
1031        'FF' dict of form factors - filled in below
1032    :param np.array G:      reciprocal metric tensor
1033    :param str hfx:    histogram id string
1034    :param str pfx:    phase id string
1035    :param dict SGData: space group info. dictionary output from SpcGroup
1036    :param dict calcControls:
1037    :param dict parmDict:
1038   
1039    :returns: dict dFdvDict: dictionary of derivatives
1040    '''
1041    phfx = pfx.split(':')[0]+hfx
1042    ast = np.sqrt(np.diag(G))
1043    Mast = twopisq*np.multiply.outer(ast,ast)
1044    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1045    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1046    FFtables = calcControls['FFtables']
1047    BLtables = calcControls['BLtables']
1048    nRef = len(refDict['RefList'])
1049    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1050        GetAtomFXU(pfx,calcControls,parmDict)
1051    mSize = len(Mdata)
1052    FF = np.zeros(len(Tdata))
1053    if 'NC' in calcControls[hfx+'histType']:
1054        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1055    elif 'X' in calcControls[hfx+'histType']:
1056        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1057        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1058    Uij = np.array(G2lat.U6toUij(Uijdata))
1059    bij = Mast*Uij.T
1060    dFdvDict = {}
1061    dFdfr = np.zeros((nRef,mSize))
1062    dFdx = np.zeros((nRef,mSize,3))
1063    dFdui = np.zeros((nRef,mSize))
1064    dFdua = np.zeros((nRef,mSize,6))
1065    dFdbab = np.zeros((nRef,2))
1066    dFdfl = np.zeros((nRef))
1067    Flack = 1.0
1068    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1069        Flack = 1.-2.*parmDict[phfx+'Flack']
1070    time0 = time.time()
1071    nref = len(refDict['RefList'])/100   
1072#reflection processing begins here - big arrays!
1073    iBeg = 0
1074    blkSize = 32       #no. of reflections in a block - optimized for speed
1075    while iBeg < nRef:
1076        iFin = min(iBeg+blkSize,nRef)
1077        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1078        H = refl.T[:3]
1079        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1080        SQfactor = 8.0*SQ*np.pi**2
1081        if 'T' in calcControls[hfx+'histType']:
1082            if 'P' in calcControls[hfx+'histType']:
1083                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1084            else:
1085                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1086            FP = np.repeat(FP.T,len(SGT),axis=0)
1087            FPP = np.repeat(FPP.T,len(SGT),axis=0)
1088        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1089        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
1090        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1091        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
1092#        GSASIIpath.IPyBreak()
1093        Uniq = np.inner(H.T,SGMT)             # array(nSGOp,3)
1094        Phi = np.inner(H.T,SGT)
1095        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1096        sinp = np.sin(phase)        #refBlk x nOps x nAtoms
1097        cosp = np.cos(phase)
1098        occ = Mdata*Fdata/len(SGT)
1099        biso = -SQfactor*Uisodata[:,nxs]
1100        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
1101        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1102        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1103        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
1104        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1105        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6))
1106        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
1107        if len(FPP.shape) > 1:
1108            fotp = np.reshape(FPP,cosp.shape)*Tcorr
1109        else:
1110            fotp = FPP*Tcorr     
1111        if 'T' in calcControls[hfx+'histType']:
1112            fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1113            fb = np.array([fot*sinp,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
1114        else:
1115            fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr])
1116            fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])
1117        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
1118        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1119        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,nEqv,nAtoms)
1120        fbx = np.array([fot*cosp,-fotp*sinp])
1121        #sum below is over Uniq
1122        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
1123        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
1124        dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
1125        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
1126        dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
1127        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
1128        if not SGData['SGInv']:
1129            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
1130            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
1131            dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
1132            dfbdfl = np.sum(np.sum(fotp*cosp,axis=-1),axis=-1)
1133            dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
1134            dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
1135            dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
1136        else:
1137            dfbdfr = np.zeros_like(dfadfr)
1138            dfbdx = np.zeros_like(dfadx)
1139            dfbdui = np.zeros_like(dfadui)
1140            dfbdua = np.zeros_like(dfadua)
1141            dfbdba = np.zeros_like(dfadba)
1142            dfadfl = 0.0
1143            dfbdfl = 0.0
1144        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1145        SA = fas[0]+fas[1]
1146        SB = fbs[0]+fbs[1]
1147        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
1148            dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
1149#            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+   \
1150#                2.*(fbs[0,:,nxs]*dfbdfr[0]+fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT)
1151            dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
1152#            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
1153#                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
1154            dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
1155#            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
1156#                2.*(fbs[0,:,nxs]*dfbdui[0]+fbs[1,:,nxs]*dfbdui[1])
1157            dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
1158#            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
1159#                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
1160        else:
1161            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
1162            dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
1163            dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
1164            dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
1165            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
1166        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)])+ \
1167                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
1168#        GSASIIpath.IPyBreak()
1169        iBeg += blkSize
1170    print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
1171        #loop over atoms - each dict entry is list of derivatives for all the reflections
1172    for i in range(len(Mdata)):
1173        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1174        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1175        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1176        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1177        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1178        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1179        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1180        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1181        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1182        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1183        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1184    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
1185    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1186    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1187    return dFdvDict
1188   
1189#def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1190#    '''Compute structure factor derivatives on single reflections - for twins only
1191#    input:
1192#   
1193#    :param dict refDict: where
1194#        'RefList' list where each ref = h,k,l,it,d,...
1195#        'FF' dict of form factors - filled in below
1196#    :param np.array G:      reciprocal metric tensor
1197#    :param str hfx:    histogram id string
1198#    :param str pfx:    phase id string
1199#    :param dict SGData: space group info. dictionary output from SpcGroup
1200#    :param dict calcControls:
1201#    :param dict parmDict:
1202#   
1203#    :returns: dict dFdvDict: dictionary of derivatives
1204#    '''
1205#    phfx = pfx.split(':')[0]+hfx
1206#    ast = np.sqrt(np.diag(G))
1207#    Mast = twopisq*np.multiply.outer(ast,ast)
1208#    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1209#    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1210#    FFtables = calcControls['FFtables']
1211#    BLtables = calcControls['BLtables']
1212#    TwDict = refDict.get('TwDict',{})           
1213#    NTL = calcControls[phfx+'NTL']
1214#    NM = calcControls[phfx+'TwinNMN']+1
1215#    TwinLaw = calcControls[phfx+'TwinLaw']
1216#    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1217#    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1218#    nTwin = len(TwinLaw)       
1219#    nRef = len(refDict['RefList'])
1220#    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1221#        GetAtomFXU(pfx,calcControls,parmDict)
1222#    mSize = len(Mdata)
1223#    FF = np.zeros(len(Tdata))
1224#    if 'NC' in calcControls[hfx+'histType']:
1225#        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1226#    elif 'X' in calcControls[hfx+'histType']:
1227#        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1228#        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1229#    Uij = np.array(G2lat.U6toUij(Uijdata))
1230#    bij = Mast*Uij.T
1231#    dFdvDict = {}
1232#    dFdfr = np.zeros((nRef,nTwin,mSize))
1233#    dFdx = np.zeros((nRef,nTwin,mSize,3))
1234#    dFdui = np.zeros((nRef,nTwin,mSize))
1235#    dFdua = np.zeros((nRef,nTwin,mSize,6))
1236#    dFdbab = np.zeros((nRef,nTwin,2))
1237#    dFdtw = np.zeros((nRef,nTwin))
1238#    time0 = time.time()
1239#    nref = len(refDict['RefList'])/100   
1240#    for iref,refl in enumerate(refDict['RefList']):
1241#        if 'T' in calcControls[hfx+'histType']:
1242#            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
1243#        H = np.array(refl[:3])
1244#        H = np.inner(H.T,TwinLaw)   #maybe array(3,nTwins) or (3)
1245#        TwMask = np.any(H,axis=-1)
1246#        if iref in TwDict:
1247#            for i in TwDict[iref]:
1248#                for n in range(NTL):
1249#                    H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1250#            TwMask = np.any(H,axis=-1)
1251#        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
1252#        SQfactor = 8.0*SQ*np.pi**2
1253#        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1254#        Bab = parmDict[phfx+'BabA']*dBabdA
1255#        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1256#        FF = refDict['FF']['FF'][iref].T[Tindx].T
1257#        Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
1258#        Phi = np.inner(H,SGT)
1259#        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1260#        sinp = np.sin(phase)
1261#        cosp = np.cos(phase)
1262#        occ = Mdata*Fdata/len(SGT)
1263#        biso = -SQfactor*Uisodata[:,nxs]
1264#        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
1265#        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1266#        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1267#        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))
1268#        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1269#        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
1270#        fot = (FF+FP-Bab)*Tcorr
1271#        fotp = FPP*Tcorr       
1272#        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])
1273#        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])
1274##        GSASIIpath.IPyBreak()
1275#        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
1276#        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1277#        if SGData['SGInv']: #centrosymmetric; B=0
1278#            fbs[0] *= 0.
1279#            fas[1] *= 0.
1280#        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
1281#        fbx = np.array([fot*cosp,-fotp*sinp])
1282#        #sum below is over Uniq
1283#        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
1284#        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1285#        dfadui = np.sum(-SQfactor*fa,axis=-2)
1286#        dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1287#        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
1288#        # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
1289#        if not SGData['SGInv']:
1290#            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
1291#            dfadba /= 2.
1292#            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
1293#            dfbdui = np.sum(-SQfactor*fb,axis=-2)
1294#            dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
1295#            dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
1296#        else:
1297#            dfbdfr = np.zeros_like(dfadfr)
1298#            dfbdx = np.zeros_like(dfadx)
1299#            dfbdui = np.zeros_like(dfadui)
1300#            dfbdua = np.zeros_like(dfadua)
1301#            dfbdba = np.zeros_like(dfadba)
1302#        SA = fas[0]+fas[1]
1303#        SB = fbs[0]+fbs[1]
1304#        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)]
1305#        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)]
1306#        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)]
1307#        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)]
1308#        if SGData['SGInv']: #centrosymmetric; B=0
1309#            dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
1310#        else:               
1311#            dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
1312#        dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1313#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1314##        GSASIIpath.IPyBreak()
1315#    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
1316#    #loop over atoms - each dict entry is list of derivatives for all the reflections
1317#    for i in range(len(Mdata)):     #these all OK?
1318#        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1319#        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1320#        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1321#        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1322#        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1323#        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1324#        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1325#        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1326#        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1327#        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1328#        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1329#    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1330#    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1331#    for i in range(nTwin):
1332#        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1333#    return dFdvDict
1334   
1335def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1336    '''Compute structure factor derivatives on blocks of reflections - for twins only
1337    faster than StructureFactorDervTw
1338    input:
1339   
1340    :param dict refDict: where
1341        'RefList' list where each ref = h,k,l,it,d,...
1342        'FF' dict of form factors - filled in below
1343    :param np.array G:      reciprocal metric tensor
1344    :param str hfx:    histogram id string
1345    :param str pfx:    phase id string
1346    :param dict SGData: space group info. dictionary output from SpcGroup
1347    :param dict calcControls:
1348    :param dict parmDict:
1349   
1350    :returns: dict dFdvDict: dictionary of derivatives
1351    '''
1352    phfx = pfx.split(':')[0]+hfx
1353    ast = np.sqrt(np.diag(G))
1354    Mast = twopisq*np.multiply.outer(ast,ast)
1355    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1356    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1357    FFtables = calcControls['FFtables']
1358    BLtables = calcControls['BLtables']
1359    TwDict = refDict.get('TwDict',{})           
1360    NTL = calcControls[phfx+'NTL']
1361    NM = calcControls[phfx+'TwinNMN']+1
1362    TwinLaw = calcControls[phfx+'TwinLaw']
1363    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1364    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1365    nTwin = len(TwinLaw)       
1366    nRef = len(refDict['RefList'])
1367    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1368        GetAtomFXU(pfx,calcControls,parmDict)
1369    mSize = len(Mdata)
1370    FF = np.zeros(len(Tdata))
1371    if 'NC' in calcControls[hfx+'histType']:
1372        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1373    elif 'X' in calcControls[hfx+'histType']:
1374        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1375        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1376    Uij = np.array(G2lat.U6toUij(Uijdata))
1377    bij = Mast*Uij.T
1378    dFdvDict = {}
1379    dFdfr = np.zeros((nRef,nTwin,mSize))
1380    dFdx = np.zeros((nRef,nTwin,mSize,3))
1381    dFdui = np.zeros((nRef,nTwin,mSize))
1382    dFdua = np.zeros((nRef,nTwin,mSize,6))
1383    dFdbab = np.zeros((nRef,nTwin,2))
1384    dFdtw = np.zeros((nRef,nTwin))
1385    time0 = time.time()
1386#reflection processing begins here - big arrays!
1387    iBeg = 0
1388    blkSize = 16       #no. of reflections in a block - optimized for speed
1389    while iBeg < nRef:
1390        iFin = min(iBeg+blkSize,nRef)
1391        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1392        H = refl.T[:3]
1393        H = np.inner(H.T,TwinLaw)   #array(3,nTwins)
1394        TwMask = np.any(H,axis=-1)
1395        for ir in range(blkSize):
1396            iref = ir+iBeg
1397            if iref in TwDict:
1398                for i in TwDict[iref]:
1399                    for n in range(NTL):
1400                        H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1401        TwMask = np.any(H,axis=-1)
1402        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1403        SQfactor = 8.0*SQ*np.pi**2
1404        if 'T' in calcControls[hfx+'histType']:
1405            if 'P' in calcControls[hfx+'histType']:
1406                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1407            else:
1408                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1409            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
1410            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
1411        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1412        Bab = np.repeat(parmDict[phfx+'BabA']*dBabdA,len(SGT)*nTwin)
1413        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1414        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
1415        Uniq = np.inner(H,SGMT)             # (nTwin,nSGOp,3)
1416        Phi = np.inner(H,SGT)
1417        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1418        sinp = np.sin(phase)
1419        cosp = np.cos(phase)
1420        occ = Mdata*Fdata/len(SGT)
1421        biso = -SQfactor*Uisodata[:,nxs]
1422        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
1423        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1424        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1425        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6))
1426        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1427        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT)
1428        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
1429        fotp = FPP*Tcorr       
1430        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
1431            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(FPP,sinp.shape)*sinp*Tcorr])
1432            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(FPP,cosp.shape)*cosp*Tcorr])
1433        else:
1434            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-FPP*sinp*Tcorr])
1435            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,FPP*cosp*Tcorr])
1436        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
1437        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1438        if SGData['SGInv']: #centrosymmetric; B=0
1439            fbs[0] *= 0.
1440            fas[1] *= 0.
1441        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,nRef,ntwi,nEqv,nAtoms)
1442        fbx = np.array([fot*cosp,-fotp*sinp])
1443        #sum below is over Uniq
1444        dfadfr = np.sum(np.sum(fa/occ,axis=-2),axis=0)        #array(2,nRef,ntwin,nAtom) Fdata != 0 avoids /0. problem
1445        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1446        dfadui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),axis=0)           
1447        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'
1448        dfadua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fa[:,:,:,:,:,nxs],axis=-3),axis=0) 
1449        if not SGData['SGInv']:
1450            dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0)        #Fdata != 0 avoids /0. problem
1451            dfadba /= 2.
1452            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
1453            dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0)
1454            dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0)
1455            dfbdua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fb[:,:,:,:,:,nxs],axis=-3),axis=0)
1456        else:
1457            dfbdfr = np.zeros_like(dfadfr)
1458            dfbdx = np.zeros_like(dfadx)
1459            dfbdui = np.zeros_like(dfadui)
1460            dfbdua = np.zeros_like(dfadua)
1461            dfbdba = np.zeros_like(dfadba)
1462        SA = fas[0]+fas[1]
1463        SB = fbs[0]+fbs[1]
1464#        GSASIIpath.IPyBreak()
1465        dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT)
1466        dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx
1467        dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui
1468        dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua
1469        if SGData['SGInv']: #centrosymmetric; B=0
1470            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
1471        else:               
1472            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
1473#        dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1474#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1475        iBeg += blkSize
1476#        GSASIIpath.IPyBreak()
1477    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
1478    #loop over atoms - each dict entry is list of derivatives for all the reflections
1479    for i in range(len(Mdata)):     #these all OK
1480        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1481        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1482        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1483        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1484        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1485        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1486        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1487        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1488        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1489        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1490        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1491    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1492    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1493    for i in range(nTwin):
1494        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1495    return dFdvDict
1496   
1497def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1498    '''
1499    Compute super structure factors for all h,k,l,m for phase - no twins
1500    puts the result, F^2, in each ref[9] in refList
1501    works on blocks of 32 reflections for speed
1502    input:
1503   
1504    :param dict refDict: where
1505        'RefList' list where each ref = h,k,l,m,it,d,...
1506        'FF' dict of form factors - filed in below
1507    :param np.array G:      reciprocal metric tensor
1508    :param str pfx:    phase id string
1509    :param dict SGData: space group info. dictionary output from SpcGroup
1510    :param dict calcControls:
1511    :param dict ParmDict:
1512
1513    '''
1514    phfx = pfx.split(':')[0]+hfx
1515    ast = np.sqrt(np.diag(G))
1516    Mast = twopisq*np.multiply.outer(ast,ast)   
1517    SGInv = SGData['SGInv']
1518    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1519    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1520    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1521    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1522    FFtables = calcControls['FFtables']
1523    BLtables = calcControls['BLtables']
1524    Flack = 1.0
1525    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1526        Flack = 1.-2.*parmDict[phfx+'Flack']
1527    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1528        GetAtomFXU(pfx,calcControls,parmDict)
1529    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1530    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1531    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1532    FF = np.zeros(len(Tdata))
1533    if 'NC' in calcControls[hfx+'histType']:
1534        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1535    elif 'X' in calcControls[hfx+'histType']:
1536        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1537        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1538    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1539    bij = Mast*Uij
1540    blkSize = 32       #no. of reflections in a block
1541    nRef = refDict['RefList'].shape[0]
1542    if not len(refDict['FF']):
1543        SQ = 1./(2.*refDict['RefList'].T[5])**2
1544        if 'N' in calcControls[hfx+'histType']:
1545            dat = G2el.getBLvalues(BLtables)
1546            refDict['FF']['El'] = dat.keys()
1547            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
1548            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1549            for iel,El in enumerate(refDict['FF']['El']):
1550                if El in MFtables:
1551                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1552        else:
1553            dat = G2el.getFFvalues(FFtables,0.)
1554            refDict['FF']['El'] = dat.keys()
1555            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1556            for iel,El in enumerate(refDict['FF']['El']):
1557                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1558    time0 = time.time()
1559#reflection processing begins here - big arrays!
1560    iBeg = 0
1561    while iBeg < nRef:
1562        iFin = min(iBeg+blkSize,nRef)
1563        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1564        H = refl.T[:4]                          #array(blkSize,4)
1565        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1566        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1567        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1568        Uniq = np.inner(H.T,SSGMT)
1569        UniqP = np.inner(HP.T,SGMT)
1570        Phi = np.inner(H.T,SSGT)
1571        if SGInv:   #if centro - expand HKL sets
1572            Uniq = np.hstack((Uniq,-Uniq))
1573            Phi = np.hstack((Phi,-Phi))
1574            UniqP = np.hstack((UniqP,-UniqP))
1575        if 'T' in calcControls[hfx+'histType']:
1576            if 'P' in calcControls[hfx+'histType']:
1577                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1578            else:
1579                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1580            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1581            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1582        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1583        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1584        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1585        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1586        sinp = np.sin(phase)
1587        cosp = np.cos(phase)
1588        biso = -SQfactor*Uisodata[:,nxs]
1589        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1590        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1591        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1592        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1593        if 'T' in calcControls[hfx+'histType']:
1594            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1595            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1596        else:
1597            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1598            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1599        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1600        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1601        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1602        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1603        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1604#        GSASIIpath.IPyBreak()
1605        if 'P' in calcControls[hfx+'histType']:
1606            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1607#            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
1608            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1609        else:
1610            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1611            refl.T[8] = np.copy(refl.T[10])               
1612            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1613        iBeg += blkSize
1614    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
1615
1616def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1617    '''
1618    Compute super structure factors for all h,k,l,m for phase - twins only
1619    puts the result, F^2, in each ref[8+im] in refList
1620    works on blocks of 32 reflections for speed
1621    input:
1622   
1623    :param dict refDict: where
1624        'RefList' list where each ref = h,k,l,m,it,d,...
1625        'FF' dict of form factors - filed in below
1626    :param np.array G:      reciprocal metric tensor
1627    :param str pfx:    phase id string
1628    :param dict SGData: space group info. dictionary output from SpcGroup
1629    :param dict calcControls:
1630    :param dict ParmDict:
1631
1632    '''
1633    phfx = pfx.split(':')[0]+hfx
1634    ast = np.sqrt(np.diag(G))
1635    Mast = twopisq*np.multiply.outer(ast,ast)   
1636    SGInv = SGData['SGInv']
1637    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1638    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1639    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1640    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1641    FFtables = calcControls['FFtables']
1642    BLtables = calcControls['BLtables']
1643    Flack = 1.0
1644    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1645        Flack = 1.-2.*parmDict[phfx+'Flack']
1646    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1647    TwDict = refDict.get('TwDict',{})           
1648    if 'S' in calcControls[hfx+'histType']:
1649        NTL = calcControls[phfx+'NTL']
1650        NM = calcControls[phfx+'TwinNMN']+1
1651        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1652        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1653        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1654    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1655        GetAtomFXU(pfx,calcControls,parmDict)
1656    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1657    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1658    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1659    FF = np.zeros(len(Tdata))
1660    if 'NC' in calcControls[hfx+'histType']:
1661        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1662    elif 'X' in calcControls[hfx+'histType']:
1663        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1664        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1665    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1666    bij = Mast*Uij
1667    blkSize = 32       #no. of reflections in a block
1668    nRef = refDict['RefList'].shape[0]
1669    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1670        SQ = 1./(2.*refDict['RefList'].T[5])**2
1671        if 'N' in calcControls[hfx+'histType']:
1672            dat = G2el.getBLvalues(BLtables)
1673            refDict['FF']['El'] = dat.keys()
1674            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
1675            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1676            for iel,El in enumerate(refDict['FF']['El']):
1677                if El in MFtables:
1678                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1679        else:
1680            dat = G2el.getFFvalues(FFtables,0.)
1681            refDict['FF']['El'] = dat.keys()
1682            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1683            for iel,El in enumerate(refDict['FF']['El']):
1684                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1685    time0 = time.time()
1686#reflection processing begins here - big arrays!
1687    iBeg = 0
1688    while iBeg < nRef:
1689        iFin = min(iBeg+blkSize,nRef)
1690        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1691        H = refl[:,:4]                          #array(blkSize,4)
1692        H3 = refl[:,:3]
1693        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1694        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1695        H3 = np.inner(H3,TwinLaw)       
1696        TwMask = np.any(HP,axis=-1)
1697        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1698            for ir in range(blkSize):
1699                iref = ir+iBeg
1700                if iref in TwDict:
1701                    for i in TwDict[iref]:
1702                        for n in range(NTL):
1703                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1704                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1705            TwMask = np.any(HP,axis=-1)
1706        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1707        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1708        Uniq = np.inner(H,SSGMT)
1709        Uniq3 = np.inner(H3,SGMT)
1710        UniqP = np.inner(HP,SGMT)
1711        Phi = np.inner(H,SSGT)
1712        if SGInv:   #if centro - expand HKL sets
1713            Uniq = np.hstack((Uniq,-Uniq))
1714            Uniq3 = np.hstack((Uniq3,-Uniq3))
1715            Phi = np.hstack((Phi,-Phi))
1716            UniqP = np.hstack((UniqP,-UniqP))
1717        if 'T' in calcControls[hfx+'histType']:
1718            if 'P' in calcControls[hfx+'histType']:
1719                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1720            else:
1721                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1722            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1723            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1724        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1725        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1726        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1727        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1728        sinp = np.sin(phase)
1729        cosp = np.cos(phase)
1730        biso = -SQfactor*Uisodata[:,nxs]
1731        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1732        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1733        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1734        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1735#        GSASIIpath.IPyBreak()
1736        if 'T' in calcControls[hfx+'histType']:
1737            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1738            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1739        else:
1740            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1741            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1742        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1743        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1744        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1745        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1746        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1747        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1748        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1749            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1750        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1751        iBeg += blkSize
1752    print 'nRef %d time %.4f\r'%(nRef,time.time()-time0)
1753
1754def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1755    '''
1756    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1757    input:
1758   
1759    :param dict refDict: where
1760        'RefList' list where each ref = h,k,l,m,it,d,...
1761        'FF' dict of form factors - filled in below
1762    :param int im: = 1 (could be eliminated)
1763    :param np.array G:      reciprocal metric tensor
1764    :param str hfx:    histogram id string
1765    :param str pfx:    phase id string
1766    :param dict SGData: space group info. dictionary output from SpcGroup
1767    :param dict SSGData: super space group info.
1768    :param dict calcControls:
1769    :param dict ParmDict:
1770   
1771    :returns: dict dFdvDict: dictionary of derivatives
1772    '''
1773    phfx = pfx.split(':')[0]+hfx
1774    ast = np.sqrt(np.diag(G))
1775    Mast = twopisq*np.multiply.outer(ast,ast)
1776    SGInv = SGData['SGInv']
1777    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1778    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1779    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1780    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1781    FFtables = calcControls['FFtables']
1782    BLtables = calcControls['BLtables']
1783    nRef = len(refDict['RefList'])
1784    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1785        GetAtomFXU(pfx,calcControls,parmDict)
1786    mSize = len(Mdata)  #no. atoms
1787    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1788    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1789    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1790    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1791    FF = np.zeros(len(Tdata))
1792    if 'NC' in calcControls[hfx+'histType']:
1793        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1794    elif 'X' in calcControls[hfx+'histType']:
1795        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1796        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1797    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1798    bij = Mast*Uij
1799    if not len(refDict['FF']):
1800        if 'N' in calcControls[hfx+'histType']:
1801            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1802        else:
1803            dat = G2el.getFFvalues(FFtables,0.)       
1804        refDict['FF']['El'] = dat.keys()
1805        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1806    dFdvDict = {}
1807    dFdfr = np.zeros((nRef,mSize))
1808    dFdx = np.zeros((nRef,mSize,3))
1809    dFdui = np.zeros((nRef,mSize))
1810    dFdua = np.zeros((nRef,mSize,6))
1811    dFdbab = np.zeros((nRef,2))
1812    dFdfl = np.zeros((nRef))
1813    dFdtw = np.zeros((nRef))
1814    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1815    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1816    dFdGz = np.zeros((nRef,mSize,5))
1817    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1818    Flack = 1.0
1819    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1820        Flack = 1.-2.*parmDict[phfx+'Flack']
1821    time0 = time.time()
1822    nRef = len(refDict['RefList'])/100
1823    for iref,refl in enumerate(refDict['RefList']):
1824        if 'T' in calcControls[hfx+'histType']:
1825            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1826        H = np.array(refl[:4])
1827        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1828        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1829        SQfactor = 8.0*SQ*np.pi**2
1830        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1831        Bab = parmDict[phfx+'BabA']*dBabdA
1832        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1833        FF = refDict['FF']['FF'][iref].T[Tindx]
1834        Uniq = np.inner(H,SSGMT)
1835        Phi = np.inner(H,SSGT)
1836        UniqP = np.inner(HP,SGMT)
1837        if SGInv:   #if centro - expand HKL sets
1838            Uniq = np.vstack((Uniq,-Uniq))
1839            Phi = np.hstack((Phi,-Phi))
1840            UniqP = np.vstack((UniqP,-UniqP))
1841        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1842        sinp = np.sin(phase)
1843        cosp = np.cos(phase)
1844        occ = Mdata*Fdata/Uniq.shape[0]
1845        biso = -SQfactor*Uisodata[:,nxs]
1846        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1847        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1848        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1849        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])                     #atoms x 6
1850        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1851        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1852        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1853        fotp = FPP*Tcorr            #ops x atoms
1854        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1855        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1856        # GfpuA is 2 x ops x atoms
1857        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1858        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1859        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1860        fag = fa*GfpuA[0]-fb*GfpuA[1]
1861        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1862       
1863        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1864        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1865        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1866        fbx = np.array([fot*cosp,-fotp*sinp])
1867        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1868        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1869        #sum below is over Uniq
1870        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1871        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1872        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1873        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1874        dfadui = np.sum(-SQfactor*fag,axis=1)
1875        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1876        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1877        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1878        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1879        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1880        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1881        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1882        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1883        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1884        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1885        dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
1886        dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
1887        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1888        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1889        if not SGData['SGInv']:   #Flack derivative
1890            dfadfl = np.sum(-FPP*Tcorr*sinp)
1891            dfbdfl = np.sum(FPP*Tcorr*cosp)
1892        else:
1893            dfadfl = 1.0
1894            dfbdfl = 1.0
1895#        GSASIIpath.IPyBreak()
1896        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1897        SA = fas[0]+fas[1]      #float = A+A'
1898        SB = fbs[0]+fbs[1]      #float = B+B'
1899        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1900            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1901            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1902                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1903            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1904                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1905            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1906                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1907            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1908                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1909            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1910                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1911            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1912                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1913            dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+  \
1914                2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1])
1915            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1916                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1917        else:                       #OK, I think
1918            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1919            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1920            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1921            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1922            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1923                           
1924            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1925            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1926            dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
1927            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1928#            GSASIIpath.IPyBreak()
1929        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1930            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1931        #loop over atoms - each dict entry is list of derivatives for all the reflections
1932        if not iref%100 :
1933            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
1934    for i in range(len(Mdata)):     #loop over atoms
1935        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1936        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1937        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1938        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1939        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1940        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1941        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1942        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1943        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1944        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1945        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1946        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
1947            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1948            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1949        nx = 0
1950        if waveTypes[i] in ['Block','ZigZag']:
1951            nx = 1 
1952            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
1953            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
1954            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
1955            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
1956            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
1957        for j in range(XSSdata.shape[1]-nx):       #loop over waves
1958            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
1959            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
1960            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
1961            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
1962            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
1963            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
1964        for j in range(USSdata.shape[1]):       #loop over waves
1965            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1966            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1967            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1968            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
1969            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
1970            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
1971            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1972            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1973            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1974            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
1975            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
1976            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
1977           
1978#        GSASIIpath.IPyBreak()
1979    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
1980    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1981    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1982    return dFdvDict
1983
1984def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1985    'Needs a doc string - no twins'
1986    phfx = pfx.split(':')[0]+hfx
1987    ast = np.sqrt(np.diag(G))
1988    Mast = twopisq*np.multiply.outer(ast,ast)
1989    SGInv = SGData['SGInv']
1990    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1991    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1992    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1993    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1994    FFtables = calcControls['FFtables']
1995    BLtables = calcControls['BLtables']
1996    nRef = len(refDict['RefList'])
1997    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1998        GetAtomFXU(pfx,calcControls,parmDict)
1999    mSize = len(Mdata)  #no. atoms
2000    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2001    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
2002    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2003    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2004    FF = np.zeros(len(Tdata))
2005    if 'NC' in calcControls[hfx+'histType']:
2006        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2007    elif 'X' in calcControls[hfx+'histType']:
2008        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2009        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2010    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2011    bij = Mast*Uij
2012    if not len(refDict['FF']):
2013        if 'N' in calcControls[hfx+'histType']:
2014            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2015        else:
2016            dat = G2el.getFFvalues(FFtables,0.)       
2017        refDict['FF']['El'] = dat.keys()
2018        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2019    dFdvDict = {}
2020    dFdfr = np.zeros((nRef,mSize))
2021    dFdx = np.zeros((nRef,mSize,3))
2022    dFdui = np.zeros((nRef,mSize))
2023    dFdua = np.zeros((nRef,mSize,6))
2024    dFdbab = np.zeros((nRef,2))
2025    dFdfl = np.zeros((nRef))
2026    dFdtw = np.zeros((nRef))
2027    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
2028    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
2029    dFdGz = np.zeros((nRef,mSize,5))
2030    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
2031    Flack = 1.0
2032    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2033        Flack = 1.-2.*parmDict[phfx+'Flack']
2034    time0 = time.time()
2035    iBeg = 0
2036    blkSize = 4       #no. of reflections in a block - optimized for speed
2037    while iBeg < nRef:
2038        iFin = min(iBeg+blkSize,nRef)
2039        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
2040        H = refl.T[:4]
2041        HP = H[:3].T+modQ*H.T[:,3:]            #projected hklm to hkl
2042        SQ = 1./(2.*refl.T[4+im])**2             # or (sin(theta)/lambda)**2
2043        SQfactor = 8.0*SQ*np.pi**2
2044        if 'T' in calcControls[hfx+'histType']:
2045            if 'P' in calcControls[hfx+'histType']:
2046                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[15])
2047            else:
2048                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[13])
2049            FP = np.repeat(FP.T,len(SGT),axis=0)
2050            FPP = np.repeat(FPP.T,len(SGT),axis=0)
2051        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2052        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
2053        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2054        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
2055        Uniq = np.inner(H.T,SSGMT)
2056        Phi = np.inner(H.T,SSGT)
2057        UniqP = np.inner(HP,SGMT)
2058        if SGInv:   #if centro - expand HKL sets
2059            Uniq = np.hstack((Uniq,-Uniq))
2060            Phi = np.hstack((Phi,-Phi))
2061            UniqP = np.hstack((UniqP,-UniqP))
2062            FF = np.vstack((FF,FF))
2063            Bab = np.concatenate((Bab,Bab))
2064        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata+Xdata).T)+Phi[:,:,nxs])
2065        sinp = np.sin(phase)
2066        cosp = np.cos(phase)
2067        occ = Mdata*Fdata/Uniq.shape[1]
2068        biso = -SQfactor*Uisodata[:,nxs]
2069        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T    #ops x atoms
2070        HbH = -np.sum(UniqP[:,:,nxs,:3]*np.inner(UniqP[:,:,:3],bij),axis=-1)  #ops x atoms
2071        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in np.reshape(UniqP,(-1,3))]) #atoms x 3x3
2072        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(iFin-iBeg,-1,6))                     #atoms x 6
2073        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2074#        GSASIIpath.IPyBreak()
2075        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2076        fot = np.reshape(FF+FP[nxs,:]-Bab[:,nxs],cosp.shape)*Tcorr     #ops x atoms
2077        fotp = FPP*Tcorr            #ops x atoms
2078        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2079        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv2(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2080        # GfpuA is 2 x ops x atoms
2081        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2082        fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
2083        fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2084        fag = fa*GfpuA[0]-fb*GfpuA[1]
2085        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2086       
2087        fas = np.sum(np.sum(fag,axis=-1),axis=-1)     # 2 x refBlk
2088        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
2089        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
2090        fbx = np.array([fot*cosp,-fotp*sinp])
2091        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2092        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2093        #sum below is over Uniq
2094        dfadfr = np.sum(fag/occ,axis=-2)        #Fdata != 0 ever avoids /0. problem
2095        dfbdfr = np.sum(fbg/occ,axis=-2)        #Fdata != 0 avoids /0. problem
2096        dfadba = np.sum(-cosp*Tcorr,axis=-2)
2097        dfbdba = np.sum(-sinp*Tcorr,axis=-2)
2098        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fag,axis=-2)
2099        dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fbg,axis=-2)
2100        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fax[:,:,:,:,nxs],axis=-3)  #2 refBlk x x nAtom x 3xyz; sum nOps
2101        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fbx[:,:,:,:,nxs],axis=-3)  #2 refBlk x x nAtom x 3xyz; sum nOps
2102        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fag[:,:,:,:,nxs],axis=-3)             #2 x nAtom x 6Uij; sum nOps
2103        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbg[:,:,:,:,nxs],axis=-3)             #2 x nAtom x 6Uij; sum nOps
2104        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
2105        dfadGf = np.sum(fa[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2)
2106        dfbdGf = np.sum(fb[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2)
2107        dfadGx = np.sum(fa[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2)
2108        dfbdGx = np.sum(fb[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2)
2109        dfadGz = np.sum(fa[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2)
2110        dfbdGz = np.sum(fb[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2)
2111        dfadGu = np.sum(fa[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2)
2112        dfbdGu = np.sum(fb[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2)   
2113        if not SGData['SGInv']:   #Flack derivative
2114            dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1)
2115            dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1)
2116        else:
2117            dfadfl = 1.0
2118            dfbdfl = 1.0
2119        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2120        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2121        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2122        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
2123            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
2124            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(Uniq)+   \
2125                2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(Uniq)
2126            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
2127                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
2128            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
2129                2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
2130            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
2131                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
2132            dFdGf[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGf[0]+fas[1,:,nxs,nxs,nxs]*dfadGf[1])+  \
2133                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGf[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGf[1])
2134            dFdGx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGx[0]+fas[1,:,nxs,nxs,nxs]*dfadGx[1])+  \
2135                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGx[0]-fbs[1,:,nxs,nxs,nxs]*dfbdGx[1])
2136            dFdGz[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadGz[0]+fas[1,:,nxs,nxs]*dfadGz[1])+  \
2137                2.*(fbs[0,:,nxs,nxs]*dfbdGz[0]+fbs[1,:,nxs,nxs]*dfbdGz[1])
2138            dFdGu[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGu[0]+fas[1,:,nxs,nxs,nxs]*dfadGu[1])+  \
2139                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGu[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGu[1])
2140        else:                       #OK, I think
2141            dFdfr[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadfr[0]+dfadfr[1])+SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq) #array(nRef,nAtom)
2142            dFdx[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1]))    #array(nRef,nAtom,3)
2143            dFdui[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadui[0]+dfadui[1])+SB[:,nxs]*(dfbdui[0]+dfbdui[1]))   #array(nRef,nAtom)
2144            dFdua[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1]))    #array(nRef,nAtom,6)
2145            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
2146                           
2147            dFdGf[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGf[0]+SB[:,nxs,nxs,nxs]*dfbdGf[1])      #array(nRef,natom,nwave,2)
2148            dFdGx[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGx[0]+SB[:,nxs,nxs,nxs]*dfbdGx[1])      #array(nRef,natom,nwave,6)
2149            dFdGz[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*dfadGz[0]+SB[:,nxs,nxs]*dfbdGz[1])      #array(nRef,natom,5)
2150            dFdGu[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGu[0]+SB[:,nxs,nxs,nxs]*dfbdGu[1])      #array(nRef,natom,nwave,12)
2151#            GSASIIpath.IPyBreak()
2152#        dFdbab[iBeg:iFin] = 2.*fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2153#            2.*fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2154        #loop over atoms - each dict entry is list of derivatives for all the reflections
2155        print ' %d derivative time %.4f\r'%(iBeg,time.time()-time0),
2156        iBeg += blkSize
2157    for i in range(len(Mdata)):     #loop over atoms
2158        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2159        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2160        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2161        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2162        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2163        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2164        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2165        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2166        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2167        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2168        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2169        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2170            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2171            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2172        nx = 0
2173        if waveTypes[i] in ['Block','ZigZag']:
2174            nx = 1 
2175            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2176            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2177            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2178            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2179            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2180        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2181            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2182            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2183            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2184            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2185            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2186            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2187        for j in range(USSdata.shape[1]):       #loop over waves
2188            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2189            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2190            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2191            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2192            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2193            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2194            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2195            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2196            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2197            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2198            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2199            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2200           
2201#        GSASIIpath.IPyBreak()
2202    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2203    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2204    return dFdvDict
2205   
2206def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2207    'Needs a doc string'
2208    phfx = pfx.split(':')[0]+hfx
2209    ast = np.sqrt(np.diag(G))
2210    Mast = twopisq*np.multiply.outer(ast,ast)
2211    SGInv = SGData['SGInv']
2212    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2213    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2214    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2215    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2216    FFtables = calcControls['FFtables']
2217    BLtables = calcControls['BLtables']
2218    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2219    TwDict = refDict.get('TwDict',{})           
2220    if 'S' in calcControls[hfx+'histType']:
2221        NTL = calcControls[phfx+'NTL']
2222        NM = calcControls[phfx+'TwinNMN']+1
2223        TwinLaw = calcControls[phfx+'TwinLaw']
2224        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
2225        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2226    nTwin = len(TwinLaw)       
2227    nRef = len(refDict['RefList'])
2228    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2229        GetAtomFXU(pfx,calcControls,parmDict)
2230    mSize = len(Mdata)  #no. atoms
2231    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2232    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
2233    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2234    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2235    FF = np.zeros(len(Tdata))
2236    if 'NC' in calcControls[hfx+'histType']:
2237        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2238    elif 'X' in calcControls[hfx+'histType']:
2239        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2240        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2241    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2242    bij = Mast*Uij
2243    if not len(refDict['FF']):
2244        if 'N' in calcControls[hfx+'histType']:
2245            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2246        else:
2247            dat = G2el.getFFvalues(FFtables,0.)       
2248        refDict['FF']['El'] = dat.keys()
2249        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2250    dFdvDict = {}
2251    dFdfr = np.zeros((nRef,nTwin,mSize))
2252    dFdx = np.zeros((nRef,nTwin,mSize,3))
2253    dFdui = np.zeros((nRef,nTwin,mSize))
2254    dFdua = np.zeros((nRef,nTwin,mSize,6))
2255    dFdbab = np.zeros((nRef,nTwin,2))
2256    dFdfl = np.zeros((nRef,nTwin))
2257    dFdtw = np.zeros((nRef,nTwin))
2258    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2259    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2260    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2261    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2262    Flack = 1.0
2263    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2264        Flack = 1.-2.*parmDict[phfx+'Flack']
2265    time0 = time.time()
2266    nRef = len(refDict['RefList'])/100
2267    for iref,refl in enumerate(refDict['RefList']):
2268        if 'T' in calcControls[hfx+'histType']:
2269            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2270        H = np.array(refl[:4])
2271        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2272        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2273        TwMask = np.any(H,axis=-1)
2274        if TwinLaw.shape[0] > 1 and TwDict:
2275            if iref in TwDict:
2276                for i in TwDict[iref]:
2277                    for n in range(NTL):
2278                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2279            TwMask = np.any(H,axis=-1)
2280        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2281        SQfactor = 8.0*SQ*np.pi**2
2282        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2283        Bab = parmDict[phfx+'BabA']*dBabdA
2284        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2285        FF = refDict['FF']['FF'][iref].T[Tindx]
2286        Uniq = np.inner(H,SSGMT)
2287        Phi = np.inner(H,SSGT)
2288        UniqP = np.inner(HP,SGMT)
2289        if SGInv:   #if centro - expand HKL sets
2290            Uniq = np.vstack((Uniq,-Uniq))
2291            Phi = np.hstack((Phi,-Phi))
2292            UniqP = np.vstack((UniqP,-UniqP))
2293        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2294        sinp = np.sin(phase)
2295        cosp = np.cos(phase)
2296        occ = Mdata*Fdata/Uniq.shape[0]
2297        biso = -SQfactor*Uisodata[:,nxs]
2298        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2299        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2300        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2301        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
2302        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2303        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2304        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2305        fotp = FPP*Tcorr            #ops x atoms
2306        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2307        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2308        # GfpuA is 2 x ops x atoms
2309        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2310        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2311        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2312        fag = fa*GfpuA[0]-fb*GfpuA[1]
2313        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2314       
2315        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2316        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2317        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2318        fbx = np.array([fot*cosp,-fotp*sinp])
2319        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2320        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2321        #sum below is over Uniq
2322        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2323        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2324        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2325        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2326        dfadui = np.sum(-SQfactor*fag,axis=1)
2327        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2328        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2329        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2330        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2331        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2332        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2333        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2334        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2335        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2336        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2337        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2338        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2339        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2340        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2341#        GSASIIpath.IPyBreak()
2342        if not SGData['SGInv'] and len(TwinLaw) == 1:   #Flack derivative
2343            dfadfl = np.sum(-FPP*Tcorr*sinp)
2344            dfbdfl = np.sum(FPP*Tcorr*cosp)
2345        else:
2346            dfadfl = 1.0
2347            dfbdfl = 1.0
2348        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2349        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2350        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2351        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)]
2352        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)]
2353        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)]
2354        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)]
2355        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2356
2357        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2358        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2359        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2360        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2361#            GSASIIpath.IPyBreak()
2362        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2363            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2364        #loop over atoms - each dict entry is list of derivatives for all the reflections
2365        if not iref%100 :
2366            print ' %d derivative time %.4f\r'%(iref,time.time()-time0),
2367    for i in range(len(Mdata)):     #loop over atoms
2368        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2369        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2370        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2371        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2372        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2373        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2374        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2375        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2376        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2377        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2378        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2379        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2380            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2381            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2382        nx = 0
2383        if waveTypes[i] in ['Block','ZigZag']:
2384            nx = 1 
2385            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2386            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2387            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2388            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2389            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2390        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2391            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2392            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2393            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2394            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2395            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2396            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2397        for j in range(USSdata.shape[1]):       #loop over waves
2398            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2399            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2400            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2401            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2402            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2403            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2404            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2405            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2406            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2407            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2408            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2409            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2410           
2411#        GSASIIpath.IPyBreak()
2412    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2413    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2414    return dFdvDict
2415   
2416def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2417    ''' Single crystal extinction function; returns extinction & derivative
2418    '''
2419    extCor = 1.0
2420    dervDict = {}
2421    dervCor = 1.0
2422    if calcControls[phfx+'EType'] != 'None':
2423        SQ = 1/(4.*ref[4+im]**2)
2424        if 'C' in parmDict[hfx+'Type']:           
2425            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2426        else:   #'T'
2427            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2428        if 'SXC' in parmDict[hfx+'Type']:
2429            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2430            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2431            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2432            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2433        elif 'SNT' in parmDict[hfx+'Type']:
2434            AV = 1.e7/parmDict[pfx+'Vol']**2
2435            PL = SQ
2436            PLZ = AV*ref[9+im]*ref[12+im]**2
2437        elif 'SNC' in parmDict[hfx+'Type']:
2438            AV = 1.e7/parmDict[pfx+'Vol']**2
2439            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2440            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2441           
2442        if 'Primary' in calcControls[phfx+'EType']:
2443            PLZ *= 1.5
2444        else:
2445            if 'C' in parmDict[hfx+'Type']:
2446                PLZ *= calcControls[phfx+'Tbar']
2447            else: #'T'
2448                PLZ *= ref[13+im]      #t-bar
2449        if 'Primary' in calcControls[phfx+'EType']:
2450            PLZ *= 1.5
2451            PSIG = parmDict[phfx+'Ep']
2452        elif 'I & II' in calcControls[phfx+'EType']:
2453            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2454        elif 'Type II' in calcControls[phfx+'EType']:
2455            PSIG = parmDict[phfx+'Es']
2456        else:       # 'Secondary Type I'
2457            PSIG = parmDict[phfx+'Eg']/PL
2458           
2459        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2460        AL = 0.025+0.285*cos2T
2461        BG = 0.02-0.025*cos2T
2462        BL = 0.15-0.2*(0.75-cos2T)**2
2463        if cos2T < 0.:
2464            BL = -0.45*cos2T
2465        CG = 2.
2466        CL = 2.
2467        PF = PLZ*PSIG
2468       
2469        if 'Gaussian' in calcControls[phfx+'EApprox']:
2470            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2471            extCor = np.sqrt(PF4)
2472            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2473        else:
2474            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2475            extCor = np.sqrt(PF4)
2476            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2477
2478        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2479        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2480            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2481        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2482            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2483        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2484            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2485               
2486    return 1./extCor,dervDict,dervCor
2487   
2488def Dict2Values(parmdict, varylist):
2489    '''Use before call to leastsq to setup list of values for the parameters
2490    in parmdict, as selected by key in varylist'''
2491    return [parmdict[key] for key in varylist] 
2492   
2493def Values2Dict(parmdict, varylist, values):
2494    ''' Use after call to leastsq to update the parameter dictionary with
2495    values corresponding to keys in varylist'''
2496    parmdict.update(zip(varylist,values))
2497   
2498def GetNewCellParms(parmDict,varyList):
2499    'Needs a doc string'
2500    newCellDict = {}
2501    Anames = ['A'+str(i) for i in range(6)]
2502    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2503    for item in varyList:
2504        keys = item.split(':')
2505        if keys[2] in Ddict:
2506            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2507            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2508            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2509    return newCellDict          # is e.g. {'0::D11':A0-D11}
2510   
2511def ApplyXYZshifts(parmDict,varyList):
2512    '''
2513    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2514   
2515    :param dict parmDict: parameter dictionary
2516    :param list varyList: list of variables (not used!)
2517    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2518
2519    '''
2520    newAtomDict = {}
2521    for item in parmDict:
2522        if 'dA' in item:
2523            parm = ''.join(item.split('d'))
2524            parmDict[parm] += parmDict[item]
2525            newAtomDict[item] = [parm,parmDict[parm]]
2526    return newAtomDict
2527   
2528def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2529    'Spherical harmonics texture'
2530    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2531    if 'T' in calcControls[hfx+'histType']:
2532        tth = parmDict[hfx+'2-theta']
2533    else:
2534        tth = refl[5+im]
2535    odfCor = 1.0
2536    H = refl[:3]
2537    cell = G2lat.Gmat2cell(g)
2538    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2539    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2540    phi,beta = G2lat.CrsAng(H,cell,SGData)
2541    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2542    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2543    for item in SHnames:
2544        L,M,N = eval(item.strip('C'))
2545        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2546        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2547        Lnorm = G2lat.Lnorm(L)
2548        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2549    return odfCor
2550   
2551def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2552    'Spherical harmonics texture derivatives'
2553    if 'T' in calcControls[hfx+'histType']:
2554        tth = parmDict[hfx+'2-theta']
2555    else:
2556        tth = refl[5+im]
2557    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2558    odfCor = 1.0
2559    dFdODF = {}
2560    dFdSA = [0,0,0]
2561    H = refl[:3]
2562    cell = G2lat.Gmat2cell(g)
2563    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2564    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2565    phi,beta = G2lat.CrsAng(H,cell,SGData)
2566    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2567    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2568    for item in SHnames:
2569        L,M,N = eval(item.strip('C'))
2570        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2571        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2572        Lnorm = G2lat.Lnorm(L)
2573        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2574        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2575        for i in range(3):
2576            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2577    return odfCor,dFdODF,dFdSA
2578   
2579def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2580    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2581    if 'T' in calcControls[hfx+'histType']:
2582        tth = parmDict[hfx+'2-theta']
2583    else:
2584        tth = refl[5+im]
2585    odfCor = 1.0
2586    H = refl[:3]
2587    cell = G2lat.Gmat2cell(g)
2588    Sangls = [0.,0.,0.]
2589    if 'Bragg' in calcControls[hfx+'instType']:
2590        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2591        IFCoup = True
2592    else:
2593        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2594        IFCoup = False
2595    phi,beta = G2lat.CrsAng(H,cell,SGData)
2596    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2597    SHnames = calcControls[phfx+'SHnames']
2598    for item in SHnames:
2599        L,N = eval(item.strip('C'))
2600        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2601        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2602        Lnorm = G2lat.Lnorm(L)
2603        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2604    return np.squeeze(odfCor)
2605   
2606def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2607    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2608    if 'T' in calcControls[hfx+'histType']:
2609        tth = parmDict[hfx+'2-theta']
2610    else:
2611        tth = refl[5+im]
2612    odfCor = 1.0
2613    dFdODF = {}
2614    H = refl[:3]
2615    cell = G2lat.Gmat2cell(g)
2616    Sangls = [0.,0.,0.]
2617    if 'Bragg' in calcControls[hfx+'instType']:
2618        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2619        IFCoup = True
2620    else:
2621        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2622        IFCoup = False
2623    phi,beta = G2lat.CrsAng(H,cell,SGData)
2624    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2625    SHnames = calcControls[phfx+'SHnames']
2626    for item in SHnames:
2627        L,N = eval(item.strip('C'))
2628        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2629        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2630        Lnorm = G2lat.Lnorm(L)
2631        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2632        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2633    return odfCor,dFdODF
2634   
2635def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2636    'March-Dollase preferred orientation correction'
2637    POcorr = 1.0
2638    MD = parmDict[phfx+'MD']
2639    if MD != 1.0:
2640        MDAxis = calcControls[phfx+'MDAxis']
2641        sumMD = 0
2642        for H in uniq:           
2643            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2644            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2645            sumMD += A**3
2646        POcorr = sumMD/len(uniq)
2647    return POcorr
2648   
2649def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2650    'Needs a doc string'
2651    POcorr = 1.0
2652    POderv = {}
2653    if calcControls[phfx+'poType'] == 'MD':
2654        MD = parmDict[phfx+'MD']
2655        MDAxis = calcControls[phfx+'MDAxis']
2656        sumMD = 0
2657        sumdMD = 0
2658        for H in uniq:           
2659            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2660            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2661            sumMD += A**3
2662            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2663        POcorr = sumMD/len(uniq)
2664        POderv[phfx+'MD'] = sumdMD/len(uniq)
2665    else:   #spherical harmonics
2666        if calcControls[phfx+'SHord']:
2667            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2668    return POcorr,POderv
2669   
2670def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2671    'Needs a doc string'
2672    if 'Debye' in calcControls[hfx+'instType']:
2673        if 'T' in calcControls[hfx+'histType']:
2674            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2675        else:
2676            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2677    else:
2678        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2679   
2680def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2681    'Needs a doc string'
2682    if 'Debye' in calcControls[hfx+'instType']:
2683        if 'T' in calcControls[hfx+'histType']:
2684            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2685        else:
2686            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2687    else:
2688        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2689       
2690def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2691    'Needs a doc string'
2692    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2693    pi2 = np.sqrt(2./np.pi)
2694    if 'T' in calcControls[hfx+'histType']:
2695        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2696        wave = refl[14+im]
2697    else:   #'C'W
2698        sth2 = sind(refl[5+im]/2.)**2
2699        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2700    c2th = 1.-2.0*sth2
2701    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2702    if 'X' in calcControls[hfx+'histType']:
2703        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2704    xfac = flv2*parmDict[phfx+'Extinction']
2705    exb = 1.0
2706    if xfac > -1.:
2707        exb = 1./np.sqrt(1.+xfac)
2708    exl = 1.0
2709    if 0 < xfac <= 1.:
2710        xn = np.array([xfac**(i+1) for i in range(6)])
2711        exl += np.sum(xn*coef)
2712    elif xfac > 1.:
2713        xfac2 = 1./np.sqrt(xfac)
2714        exl = pi2*(1.-0.125/xfac)*xfac2
2715    return exb*sth2+exl*(1.-sth2)
2716   
2717def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2718    'Needs a doc string'
2719    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2720    pi2 = np.sqrt(2./np.pi)
2721    if 'T' in calcControls[hfx+'histType']:
2722        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2723        wave = refl[14+im]
2724    else:   #'C'W
2725        sth2 = sind(refl[5+im]/2.)**2
2726        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2727    c2th = 1.-2.0*sth2
2728    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2729    if 'X' in calcControls[hfx+'histType']:
2730        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2731    xfac = flv2*parmDict[phfx+'Extinction']
2732    dbde = -500.*flv2
2733    if xfac > -1.:
2734        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2735    dlde = 0.
2736    if 0 < xfac <= 1.:
2737        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2738        dlde = np.sum(xn*coef)
2739    elif xfac > 1.:
2740        xfac2 = 1./np.sqrt(xfac)
2741        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2742       
2743    return dbde*sth2+dlde*(1.-sth2)
2744   
2745def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2746    'Needs a doc string'    #need powder extinction!
2747    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2748    if 'X' in parmDict[hfx+'Type']:
2749        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2750    POcorr = 1.0
2751    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2752        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2753    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2754        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2755    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2756        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2757    Icorr *= POcorr
2758    AbsCorr = 1.0
2759    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2760    Icorr *= AbsCorr
2761    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2762    Icorr *= ExtCorr
2763    return Icorr,POcorr,AbsCorr,ExtCorr
2764   
2765def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2766    'Needs a doc string'    #need powder extinction derivs!
2767    dIdsh = 1./parmDict[hfx+'Scale']
2768    dIdsp = 1./parmDict[phfx+'Scale']
2769    if 'X' in parmDict[hfx+'Type']:
2770        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2771        dIdPola /= pola
2772    else:       #'N'
2773        dIdPola = 0.0
2774    dFdODF = {}
2775    dFdSA = [0,0,0]
2776    dIdPO = {}
2777    if pfx+'SHorder' in parmDict:
2778        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2779        for iSH in dFdODF:
2780            dFdODF[iSH] /= odfCor
2781        for i in range(3):
2782            dFdSA[i] /= odfCor
2783    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2784        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2785        for iPO in dIdPO:
2786            dIdPO[iPO] /= POcorr
2787    if 'T' in parmDict[hfx+'Type']:
2788        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2789        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2790    else:
2791        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2792        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2793    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2794       
2795def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2796    'Needs a doc string'
2797    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
2798        costh = cosd(refl[5+im]/2.)
2799        #crystallite size
2800        if calcControls[phfx+'SizeType'] == 'isotropic':
2801            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2802        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2803            H = np.array(refl[:3])
2804            P = np.array(calcControls[phfx+'SizeAxis'])
2805            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2806            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2807            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2808        else:           #ellipsoidal crystallites
2809            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2810            H = np.array(refl[:3])
2811            lenR = G2pwd.ellipseSize(H,Sij,GB)
2812            Sgam = 1.8*wave/(np.pi*costh*lenR)
2813        #microstrain               
2814        if calcControls[phfx+'MustrainType'] == 'isotropic':
2815            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2816        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2817            H = np.array(refl[:3])
2818            P = np.array(calcControls[phfx+'MustrainAxis'])
2819            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2820            Si = parmDict[phfx+'Mustrain;i']
2821            Sa = parmDict[phfx+'Mustrain;a']
2822            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2823        else:       #generalized - P.W. Stephens model
2824            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2825            Sum = 0
2826            for i,strm in enumerate(Strms):
2827                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2828            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2829    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2830        #crystallite size
2831        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2832            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2833        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2834            H = np.array(refl[:3])
2835            P = np.array(calcControls[phfx+'SizeAxis'])
2836            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2837            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2838            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2839        else:           #ellipsoidal crystallites   #OK
2840            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2841            H = np.array(refl[:3])
2842            lenR = G2pwd.ellipseSize(H,Sij,GB)
2843            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2844        #microstrain               
2845        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2846            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2847        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2848            H = np.array(refl[:3])
2849            P = np.array(calcControls[phfx+'MustrainAxis'])
2850            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2851            Si = parmDict[phfx+'Mustrain;i']
2852            Sa = parmDict[phfx+'Mustrain;a']
2853            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2854        else:       #generalized - P.W. Stephens model  OK
2855            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2856            Sum = 0
2857            for i,strm in enumerate(Strms):
2858                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2859            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2860           
2861    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2862    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2863    sig /= ateln2
2864    return sig,gam
2865       
2866def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2867    'Needs a doc string'
2868    gamDict = {}
2869    sigDict = {}
2870    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
2871        costh = cosd(refl[5+im]/2.)
2872        tanth = tand(refl[5+im]/2.)
2873        #crystallite size derivatives
2874        if calcControls[phfx+'SizeType'] == 'isotropic':
2875            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2876            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
2877            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2878        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2879            H = np.array(refl[:3])
2880            P = np.array(calcControls[phfx+'SizeAxis'])
2881            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2882            Si = parmDict[phfx+'Size;i']
2883            Sa = parmDict[phfx+'Size;a']
2884            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2885            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2886            Sgam = gami*sqtrm
2887            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2888            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2889            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2890            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2891            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2892            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2893        else:           #ellipsoidal crystallites
2894            const = 1.8*wave/(np.pi*costh)
2895            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2896            H = np.array(refl[:3])
2897            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2898            Sgam = const/lenR
2899            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2900                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2901                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2902        gamDict[phfx+'Size;mx'] = Sgam
2903        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2904               
2905        #microstrain derivatives               
2906        if calcControls[phfx+'MustrainType'] == 'isotropic':
2907            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2908            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2909            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2910        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2911            H = np.array(refl[:3])
2912            P = np.array(calcControls[phfx+'MustrainAxis'])
2913            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2914            Si = parmDict[phfx+'Mustrain;i']
2915            Sa = parmDict[phfx+'Mustrain;a']
2916            gami = 0.018*Si*Sa*tanth/np.pi
2917            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2918            Mgam = gami/sqtrm
2919            dsi = -gami*Si*cosP**2/sqtrm**3
2920            dsa = -gami*Sa*sinP**2/sqtrm**3
2921            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2922            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2923            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2924            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2925        else:       #generalized - P.W. Stephens model
2926            const = 0.018*refl[4+im]**2*tanth/np.pi
2927            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2928            Sum = 0
2929            for i,strm in enumerate(Strms):
2930                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2931                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2932                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2933            Mgam = const*np.sqrt(Sum)
2934            for i in range(len(Strms)):
2935                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2936                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2937        gamDict[phfx+'Mustrain;mx'] = Mgam
2938        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2939    else:   #'T'OF - All checked & OK
2940        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2941            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2942            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2943            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2944        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2945            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2946            H = np.array(refl[:3])
2947            P = np.array(calcControls[phfx+'SizeAxis'])
2948            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2949            Si = parmDict[phfx+'Size;i']
2950            Sa = parmDict[phfx+'Size;a']
2951            gami = const/(Si*Sa)
2952            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2953            Sgam = gami*sqtrm
2954            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2955            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2956            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2957            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2958            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2959            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2960        else:           #OK  ellipsoidal crystallites
2961            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2962            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2963            H = np.array(refl[:3])
2964            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2965            Sgam = const/lenR
2966            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2967                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2968                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2969        gamDict[phfx+'Size;mx'] = Sgam  #OK
2970        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2971               
2972        #microstrain derivatives               
2973        if calcControls[phfx+'MustrainType'] == 'isotropic':
2974            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2975            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2976            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2977        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2978            H = np.array(refl[:3])
2979            P = np.array(calcControls[phfx+'MustrainAxis'])
2980            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2981            Si = parmDict[phfx+'Mustrain;i']
2982            Sa = parmDict[phfx+'Mustrain;a']
2983            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2984            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2985            Mgam = gami/sqtrm
2986            dsi = -gami*Si*cosP**2/sqtrm**3
2987            dsa = -gami*Sa*sinP**2/sqtrm**3
2988            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2989            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2990            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2991            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2992        else:       #generalized - P.W. Stephens model OK
2993            pwrs = calcControls[phfx+'MuPwrs']
2994            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2995            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2996            Sum = 0
2997            for i,strm in enumerate(Strms):
2998                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2999                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
3000                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
3001            Mgam = const*np.sqrt(Sum)
3002            for i in range(len(Strms)):
3003                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
3004                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
3005        gamDict[phfx+'Mustrain;mx'] = Mgam
3006        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
3007       
3008    return sigDict,gamDict
3009       
3010def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
3011    'Needs a doc string'
3012    if im:
3013        h,k,l,m = refl[:4]
3014        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
3015        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
3016    else:
3017        h,k,l = refl[:3]
3018        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
3019    refl[4+im] = d
3020    if 'C' in calcControls[hfx+'histType']:
3021        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
3022        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
3023        if 'Bragg' in calcControls[hfx+'instType']:
3024            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
3025                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
3026        else:               #Debye-Scherrer - simple but maybe not right
3027            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
3028    elif 'T' in calcControls[hfx+'histType']:
3029        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
3030        #do I need sample position effects - maybe?
3031    return pos
3032
3033def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
3034    'Needs a doc string'
3035    dpr = 180./np.pi
3036    if im:
3037        h,k,l,m = refl[:4]
3038        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
3039        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
3040        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
3041    else:
3042        m = 0
3043        h,k,l = refl[:3]       
3044        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
3045    dst = np.sqrt(dstsq)
3046    dsp = 1./dst
3047    if 'C' in calcControls[hfx+'histType']:
3048        pos = refl[5+im]-parmDict[hfx+'Zero']
3049        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
3050        dpdw = const*dst
3051        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
3052        dpdZ = 1.0
3053        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
3054            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
3055        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
3056        if 'Bragg' in calcControls[hfx+'instType']:
3057            dpdSh = -4.*shft*cosd(pos/2.0)
3058            dpdTr = -shft*sind(pos)*100.0
3059            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
3060        else:               #Debye-Scherrer - simple but maybe not right
3061            dpdXd = -shft*cosd(pos)
3062            dpdYd = -shft*sind(pos)
3063            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
3064    elif 'T' in calcControls[hfx+'histType']:
3065        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
3066        dpdZ = 1.0
3067        dpdDC = dsp
3068        dpdDA = dsp**2
3069        dpdDB = 1./dsp
3070        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
3071            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
3072        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
3073           
3074def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
3075    'Needs a doc string'
3076    laue = SGData['SGLaue']
3077    uniq = SGData['SGUniq']
3078    h,k,l = refl[:3]
3079    if laue in ['m3','m3m']:
3080        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
3081            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
3082    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3083        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
3084    elif laue in ['3R','3mR']:
3085        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
3086    elif laue in ['4/m','4/mmm']:
3087        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
3088    elif laue in ['mmm']:
3089        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
3090    elif laue in ['2/m']:
3091        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
3092        if uniq == 'a':
3093            Dij += parmDict[phfx+'D23']*k*l
3094        elif uniq == 'b':
3095            Dij += parmDict[phfx+'D13']*h*l
3096        elif uniq == 'c':
3097            Dij += parmDict[phfx+'D12']*h*k
3098    else:
3099        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
3100            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
3101    if 'C' in calcControls[hfx+'histType']:
3102        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
3103    else:
3104        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
3105           
3106def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
3107    'Needs a doc string'
3108    laue = SGData['SGLaue']
3109    uniq = SGData['SGUniq']
3110    h,k,l = refl[:3]
3111    if laue in ['m3','m3m']:
3112        dDijDict = {phfx+'D11':h**2+k**2+l**2,
3113            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
3114    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3115        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
3116    elif laue in ['3R','3mR']:
3117        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
3118    elif laue in ['4/m','4/mmm']:
3119        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
3120    elif laue in ['mmm']:
3121        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
3122    elif laue in ['2/m']:
3123        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
3124        if uniq == 'a':
3125            dDijDict[phfx+'D23'] = k*l
3126        elif uniq == 'b':
3127            dDijDict[phfx+'D13'] = h*l
3128        elif uniq == 'c':
3129            dDijDict[phfx+'D12'] = h*k
3130    else:
3131        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
3132            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
3133    if 'C' in calcControls[hfx+'histType']:
3134        for item in dDijDict:
3135            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
3136    else:
3137        for item in dDijDict:
3138            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
3139    return dDijDict
3140   
3141def GetDij(phfx,SGData,parmDict):
3142    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
3143    return G2spc.HStrainVals(HSvals,SGData)
3144               
3145def GetFobsSq(Histograms,Phases,parmDict,calcControls):
3146    'Needs a doc string'
3147    histoList = Histograms.keys()
3148    histoList.sort()
3149    for histogram in histoList:
3150        if 'PWDR' in histogram[:4]:
3151            Histogram = Histograms[histogram]
3152            hId = Histogram['hId']
3153            hfx = ':%d:'%(hId)
3154            Limits = calcControls[hfx+'Limits']
3155            if 'C' in calcControls[hfx+'histType']:
3156                shl = max(parmDict[hfx+'SH/L'],0.0005)
3157                Ka2 = False
3158                kRatio = 0.0
3159                if hfx+'Lam1' in parmDict.keys():
3160                    Ka2 = True
3161                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3162                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
3163            x,y,w,yc,yb,yd = Histogram['Data']
3164            xMask = ma.getmaskarray(x)
3165            xB = np.searchsorted(x,Limits[0])
3166            xF = np.searchsorted(x,Limits[1])
3167            ymb = np.array(y-yb)
3168            ymb = np.where(ymb,ymb,1.0)
3169            ycmb = np.array(yc-yb)
3170            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3171            refLists = Histogram['Reflection Lists']
3172            for phase in refLists:
3173                if phase not in Phases:     #skips deleted or renamed phases silently!
3174                    continue
3175                Phase = Phases[phase]
3176                im = 0
3177                if Phase['General'].get('Modulated',False):
3178                    im = 1
3179                pId = Phase['pId']
3180                phfx = '%d:%d:'%(pId,hId)
3181                refDict = refLists[phase]
3182                sumFo = 0.0
3183                sumdF = 0.0
3184                sumFosq = 0.0
3185                sumdFsq = 0.0
3186                sumInt = 0.0
3187                nExcl = 0
3188                for refl in refDict['RefList']:
3189                    if 'C' in calcControls[hfx+'histType']:
3190                        yp = np.zeros_like(yb)
3191                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3192                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
3193                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
3194                        iFin2 = iFin
3195                        if not iBeg+iFin:       #peak below low limit - skip peak
3196                            continue
3197                        if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
3198                            refl[3+im] *= -1
3199                            nExcl += 1
3200                            continue
3201                        elif not iBeg-iFin:     #peak above high limit - done
3202                            break
3203                        elif iBeg < iFin:
3204                            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
3205                            sumInt += refl[11+im]*refl[9+im]
3206                            if Ka2:
3207                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3208                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3209                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
3210                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
3211                                if iFin2 > iBeg2: 
3212                                    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
3213                                    sumInt += refl[11+im]*refl[9+im]*kRatio
3214                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
3215                               
3216                    elif 'T' in calcControls[hfx+'histType']:
3217                        yp = np.zeros_like(yb)
3218                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
3219                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
3220                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
3221                        if not iBeg+iFin:       #peak below low limit - skip peak
3222                            continue
3223                        if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
3224                            refl[3+im] *= -1
3225                            nExcl += 1
3226                            continue
3227                        elif not iBeg-iFin:     #peak above high limit - done
3228                            break
3229                        if iBeg < iFin:
3230                            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
3231                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
3232                            sumInt += refl[11+im]*refl[9+im]
3233                    Fo = np.sqrt(np.abs(refl[8+im]))
3234                    Fc = np.sqrt(np.abs(refl[9]+im))
3235                    sumFo += Fo
3236                    sumFosq += refl[8+im]**2
3237                    sumdF += np.abs(Fo-Fc)
3238                    sumdFsq += (refl[8+im]-refl[9+im])**2
3239                if sumFo:
3240                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3241                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3242                else:
3243                    Histogram['Residuals'][phfx+'Rf'] = 100.
3244                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3245                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3246                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3247                Histogram['Residuals']['hId'] = hId
3248        elif 'HKLF' in histogram[:4]:
3249            Histogram = Histograms[histogram]
3250            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3251               
3252def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3253    'Needs a doc string'
3254   
3255    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3256        U = parmDict[hfx+'U']
3257        V = parmDict[hfx+'V']
3258        W = parmDict[hfx+'W']
3259        X = parmDict[hfx+'X']
3260        Y = parmDict[hfx+'Y']
3261        tanPos = tand(refl[5+im]/2.0)
3262        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3263        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3264        sig = max(0.001,sig)
3265        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
3266        gam = max(0.001,gam)
3267        return sig,gam
3268               
3269    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3270        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3271            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
3272        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
3273        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3274        sig += Ssig
3275        gam += Sgam
3276        return sig,gam
3277       
3278    def GetReflAlpBet(refl,im,hfx,parmDict):
3279        alp = parmDict[hfx+'alpha']/refl[4+im]
3280        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3281        return alp,bet
3282       
3283    hId = Histogram['hId']
3284    hfx = ':%d:'%(hId)
3285    bakType = calcControls[hfx+'bakType']
3286    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
3287    yc = np.zeros_like(yb)
3288    cw = np.diff(x)
3289    cw = np.append(cw,cw[-1])
3290       
3291    if 'C' in calcControls[hfx+'histType']:   
3292        shl = max(parmDict[hfx+'SH/L'],0.002)
3293        Ka2 = False
3294        if hfx+'Lam1' in parmDict.keys():
3295            wave = parmDict[hfx+'Lam1']
3296            Ka2 = True
3297            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3298            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3299        else:
3300            wave = parmDict[hfx+'Lam']
3301    for phase in Histogram['Reflection Lists']:
3302        refDict = Histogram['Reflection Lists'][phase]
3303        if phase not in Phases:     #skips deleted or renamed phases silently!
3304            continue
3305        Phase = Phases[phase]
3306        pId = Phase['pId']
3307        pfx = '%d::'%(pId)
3308        phfx = '%d:%d:'%(pId,hId)
3309        hfx = ':%d:'%(hId)
3310        SGData = Phase['General']['SGData']
3311        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3312        im = 0
3313        if Phase['General'].get('Modulated',False):
3314            SSGData = Phase['General']['SSGData']
3315            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3316            im = 1  #offset in SS reflection list
3317            #??
3318        Dij = GetDij(phfx,SGData,parmDict)
3319        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
3320        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3321        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3322            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3323        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3324        Vst = np.sqrt(nl.det(G))    #V*
3325        if not Phase['General'].get('doPawley'):
3326            time0 = time.time()
3327            if im:
3328                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3329            else:
3330                StructureFactor2(refDict,G,hfx,pfx,SGData