source: trunk/GSASIIstrMath.py @ 2506

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

some work on automatic mag constraints
move a dlg.Destroy() in OnSeqPeakFit? & OnIndexPeaks? - get rid of orphan wait cursor & progress dialog
more work on mag moment derivs - still wrong! Other derivs now all OK

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