source: trunk/GSASIIstrMath.py @ 4533

Last change on this file since 4533 was 4533, checked in by vondreele, 3 years ago

fix math error in incomm mag str fctr stuff $ expand equations

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