source: trunk/GSASIIstrMath.py @ 4519

Last change on this file since 4519 was 4519, checked in by vondreele, 17 months ago

G2fpaGUI - minor cleanup
G2lattice - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2math - add getPinkalpha, getPinkbeta % deriv routines for pink beam peak shapes

modify setPeakparams to use them

G2plot - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2pwd - add pink beam peak shape function - same as TOF peak shape; scale sig & gam to be in centidegrees

  • add wtFactor to FitPeaks? fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2pwdGUI - test on 'T' instead of 'C' so that 'B' PWDR types properly handled

  • add wtFactor to DoPeakFit? to fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2strMath - latest in incommensurate mag str factors - no real improvement

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