source: trunk/GSASIIstrMath.py @ 2837

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

fix aniso Uij derivatives - were messing up Uii on trig/hex 3-fold axes

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