source: trunk/GSASIIstrMath.py @ 2499

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

move setting of G2frame.HKL from MakeReflectionTable? to ShowReflTable?; makes reflection tic mark tool tip match selected reflection list
work on Mag moment derivs.

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