source: trunk/GSASIIstrMath.py @ 2111

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

fix to mod vec edit & search
comment out unused StructureFactorDerv? versions

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