source: trunk/GSASIIstrMath.py @ 3423

Last change on this file since 3423 was 3423, checked in by vondreele, 5 years ago

fix chemical composition restraint. Wasn't updating with new fracs after each cycle. Now OK.
new routine GetAtomFracById? in G2math called in G2strMath
fix a wx.wx.... bug in G2ctrlGUI!
fix G2testplot to use wx-Phoenix; ditto testSSymbols

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