source: branch/MPbranch/GSASIIElem.py @ 3377

Last change on this file since 3377 was 489, checked in by vondreele, 13 years ago

fix to speed up scat fac calculations & skip texture calcs if order=0
  • Property svn:keywords set to Date Author Revision URL Id
File size: 13.9 KB
Line 
1"""Element: functions for element types
2   Copyright: 2008, Robert B. Von Dreele & Brian H. Toby (Argonne National Laboratory)
3"""
4########### SVN repository information ###################
5# $Date: 2012-02-22 22:08:40 +0000 (Wed, 22 Feb 2012) $
6# $Author: vondreele $
7# $Revision: 489 $
8# $URL: branch/MPbranch/GSASIIElem.py $
9# $Id: GSASIIElem.py 489 2012-02-22 22:08:40Z vondreele $
10########### SVN repository information ###################
11
12import math
13import os.path
14import GSASIIpath
15import numpy as np
16
17def GetFormFactorCoeff(El):
18    """Read X-ray form factor coefficients from `atomdata.asc` file
19
20    :param El: element 1-2 character symbol case irrevelant
21    :return: `FormFactors`: list of form factor dictionaries
22   
23    Each X-ray form factor dictionary is:
24   
25    * `Symbol`: 4 character element symbol with valence (e.g. 'NI+2')
26    * `Z`: atomic number
27    * `fa`: 4 A coefficients
28    * `fb`: 4 B coefficients
29    * `fc`: C coefficient
30   
31    """
32    ElS = El.upper()
33    ElS = ElS.rjust(2)
34    filename = os.path.join(os.path.split(__file__)[0],'atmdata.dat')
35    try:
36        FFdata = open(filename,'Ur')
37    except:
38        print "**** ERROR - File atmdata.dat not found in directory %s" % sys.path[0]
39        sys.exit()
40    S = '1'
41    FormFactors = []
42    while S:
43        S = FFdata.readline()
44        if S[3:5] == ElS:
45            if S[5:6] != '_' and S[8] not in ['N','M']:
46                Z=int(S[:2])
47                Symbol = S[3:7].strip()
48                S = S[12:]
49                fa = (float(S[:7]),float(S[14:21]),float(S[28:35]),float(S[42:49]))
50                fb = (float(S[7:14]),float(S[21:28]),float(S[35:42]),float(S[49:56]))
51                FormFac = {'Symbol':Symbol,'Z':Z,'fa':fa,'fb':fb,'fc':float(S[56:63])}
52                FormFactors.append(FormFac)               
53    FFdata.close()
54    return FormFactors
55   
56def GetFFC5(ElSym):
57    '''Get 5 term form factor and Compton scattering data
58    @param ElSym: str(1-2 character element symbol with proper case);
59    @return El: dictionary with 5 term form factor & compton coefficients
60    '''
61    import FormFactors as FF
62    El = {}
63    FF5 = FF.FFac5term[ElSym]
64    El['fa'] = FF5[:5]
65    El['fc'] = FF5[5]
66    El['fb'] = FF5[6:]
67    Cmp5 = FF.Compton[ElSym]
68    El['cmpz'] = Cmp5[0]
69    El['cmpa'] = Cmp5[1:6]
70    El['cmpb'] = Cmp5[6:]
71    return El
72   
73def GetAtomInfo(El):
74   
75    import ElementTable as ET
76    Elements = []
77    for elem in ET.ElTable:
78        Elements.append(elem[0][0])
79    if len(El) in [2,4]:
80        ElS = El.upper()[:2].rjust(2)
81    else:
82        ElS = El.upper()[:1].rjust(2)
83    filename = os.path.join(os.path.split(__file__)[0],'atmdata.dat')
84    try:
85        FFdata = open(filename,'Ur')
86    except:
87        print '**** ERROR - File atmdata.dat not found in directory %s' % sys.path[0]
88        sys.exit()
89    S = '1'
90    AtomInfo = {}
91    Isotopes = {}
92    Mass = []
93    while S:
94        S = FFdata.readline()
95        if S[3:5] == ElS:
96            if S[5:6] == '_':
97                if not Mass:                                 #picks 1st one; natural abundance or 1st isotope
98                    Mass = float(S[10:19])
99                if S[6] in [' ','1','2','3','4','5','6','7','8','9']:                       
100                    isoName = S[6:9]
101                    if isoName == '   ':
102                        isoName = 'Nat. Abund.'              #natural abundance
103                    if S[76:78] in ['LS','BW']:     #special anomalous scattering length info
104                        St = [S[10:19],S[19:25],S[25:31],S[31:38],S[38:44],S[44:50],
105                            S[50:56],S[56:62],S[62:68],S[68:74],]
106                        Vals = []
107                        for item in St:
108                            if item.strip():
109                                Vals.append(float(item.strip()))
110                        Isotopes[isoName.rstrip()] = Vals                       
111                    else:
112                        Isotopes[isoName.rstrip()] = [float(S[10:19]),float(S[19:25])]
113                elif S[5:9] == '_SIZ':
114                    Z=int(S[:2])
115                    Symbol = S[3:5].strip().lower().capitalize()
116                    Drad = float(S[12:22])
117                    Arad = float(S[22:32])
118                    Vdrad = float(S[32:38])
119                    Color = ET.ElTable[Elements.index(Symbol)][6]
120    FFdata.close()
121    AtomInfo={'Symbol':Symbol,'Isotopes':Isotopes,'Mass':Mass,'Z':Z,'Drad':Drad,'Arad':Arad,'Vdrad':Vdrad,'Color':Color}   
122    return AtomInfo
123     
124def GetXsectionCoeff(El):
125    """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm
126    @param El: 2 character element symbol
127    @return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc
128    each dictionary is:
129    'OrbName': Orbital name read from file
130    'IfBe' 0/2 depending on orbital
131    'BindEn': binding energy
132    'BB': BindEn/0.02721
133    'XSectIP': 5 cross section inflection points
134    'ElEterm': energy correction term
135    'SEdge': absorption edge for orbital
136    'Nval': 10/11 depending on IfBe
137    'LEner': 10/11 values of log(energy)
138    'LXSect': 10/11 values of log(cross section)
139    """
140    AU = 2.80022e+7
141    C1 = 0.02721
142    ElS = El.upper()
143    ElS = ElS.ljust(2)
144    filename = os.path.join(os.path.split(__file__)[0],'Xsect.dat')
145    try:
146        xsec = open(filename,'Ur')
147    except:
148        print '**** ERROR - File Xsect.dat not found in directory %s' % sys.path[0]
149        sys.exit()
150    S = '1'
151    Orbs = []
152    while S:
153        S = xsec.readline()
154        if S[:2] == ElS:
155            S = S[:-1]+xsec.readline()[:-1]+xsec.readline()
156            OrbName = S[9:14]
157            S = S[14:]
158            IfBe = int(S[0])
159            S = S[1:]
160            val = S.split()
161            BindEn = float(val[0])
162            BB = BindEn/C1
163            Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB}
164            Energy = []
165            XSect = []
166            for i in range(11):
167                Energy.append(float(val[2*i+1]))
168                XSect.append(float(val[2*i+2]))
169            XSecIP = []
170            for i in range(5): XSecIP.append(XSect[i+5]/AU)
171            Orb['XSecIP'] = XSecIP
172            if IfBe == 0:
173                Orb['SEdge'] = XSect[10]/AU
174                Nval = 11
175            else:
176                Orb['ElEterm'] = XSect[10]
177                del Energy[10]
178                del XSect[10]
179                Nval = 10
180                Orb['SEdge'] = 0.0
181            Orb['Nval'] = Nval
182            D = dict(zip(Energy,XSect))
183            Energy.sort()
184            X = []
185            for key in Energy:
186                X.append(D[key])
187            XSect = X
188            LEner = []
189            LXSect = []
190            for i in range(Nval):
191                LEner.append(math.log(Energy[i]))
192                if XSect[i] > 0.0:
193                    LXSect.append(math.log(XSect[i]))
194                else:
195                    LXSect.append(0.0)
196            Orb['LEner'] = LEner
197            Orb['LXSect'] = LXSect
198            Orbs.append(Orb)
199    xsec.close()
200    return Orbs
201   
202def GetMagFormFacCoeff(El):
203    """Read magnetic form factor data from atomdata.asc file
204    @param El: 2 character element symbol
205    @return: MagFormFactors: list of all magnetic form factors dictionaries for element El.
206    each dictionary contains:
207    'Symbol':Symbol
208    'Z':Z
209    'mfa': 4 MA coefficients
210    'nfa': 4 NA coefficients
211    'mfb': 4 MB coefficients
212    'nfb': 4 NB coefficients
213    'mfc': MC coefficient
214    'nfc': NC coefficient
215    """
216    ElS = El.upper()
217    ElS = ElS.rjust(2)
218    filename = os.path.join(os.path.split(__file__)[0],'atmdata.dat')
219    try:
220        FFdata = open(filename,'Ur')
221    except:
222        print '**** ERROR - File atmdata.dat not found in directory %s' % sys.path[0]
223        sys.exit()
224    S = '1'
225    MagFormFactors = []
226    while S:
227        S = FFdata.readline()
228        if S[3:5] == ElS:
229            if S[8:9] == 'M':
230                SN = FFdata.readline()               #'N' is assumed to follow 'M' in Atomdata.asc
231                Z=int(S[:2])
232                Symbol = S[3:7]
233                S = S[12:]
234                SN = SN[12:]
235                mfa = (float(S[:7]),float(S[14:21]),float(S[28:35]),float(S[42:49]))
236                mfb = (float(S[7:14]),float(S[21:28]),float(S[35:42]),float(S[49:56]))
237                nfa = (float(SN[:7]),float(SN[14:21]),float(SN[28:35]),float(SN[42:49]))
238                nfb = (float(SN[7:14]),float(SN[21:28]),float(SN[35:42]),float(SN[49:56]))
239                FormFac = {'Symbol':Symbol,'Z':Z,'mfa':mfa,'nfa':nfa,'mfb':mfb,'nfb':nfb,
240                    'mfc':float(S[56:63]),'nfc':float(SN[56:63])}
241                MagFormFactors.append(FormFac)
242    FFdata.close()
243    return MagFormFactors
244
245def ScatFac(El, SQ):
246    """compute value of form factor
247    @param El: element dictionary defined in GetFormFactorCoeff
248    @param SQ: (sin-theta/lambda)**2
249    @return: real part of form factor
250    """
251    fa = np.array(El['fa'])
252    fb = np.array(El['fb'])
253    t = -fb[:,np.newaxis]*SQ
254    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
255       
256def BlenRes(Elist,BLtables,wave):
257    FP = np.zeros(len(Elist))
258    FPP = np.zeros(len(Elist))
259    Emev = 81.80703/wave**2
260    for i,El in enumerate(Elist):
261        BL = BLtables[El]
262        if len(BL) >= 6:
263            Emev = 81.80703/wave**2
264            G2 = BL[5]**2
265            T = [Emev-BL[4],0,0]
266            D = [T**2+G2,0,0]
267            fp = T/D
268            fpp = 1.0/D
269            if len(BL) == 8:
270                T = Emev-BL[7]
271                D = T**2+G2
272                fp += BL[6]*T/D
273                fpp += BL[6]/D
274            if len(BL) == 10:
275                T = Emev-BL[9]
276                D = T**2+G2
277                fp += BL[8]*T/D
278                fpp += BL[8]/D
279            FP[i] = (BL[2]*fp)
280            FPP[i] = (-BL[3]*fpp)
281        else:
282            FP[i] = 0.0
283            FPP[i] = 0.0
284    return FP,FPP
285   
286#def BlenRes(BLdata,wave):
287#    FP = []
288#    FPP = []
289#    Emev = 81.80703/wave**2
290#    for BL in BLdata:
291#        if len(BL) >= 6:
292#            Emev = 81.80703/wave**2
293#            G2 = BL[5]**2
294#            T = [Emev-BL[4],0,0]
295#            D = [T**2+G2,0,0]
296#            fp = T/D
297#            fpp = 1.0/D
298#            if len(BL) == 8:
299#                T = Emev-BL[7]
300#                D = T**2+G2
301#                fp += BL[6]*T/D
302#                fpp += BL[6]/D
303#            if len(BL) == 10:
304#                T = Emev-BL[9]
305#                D = T**2+G2
306#                fp += BL[8]*T/D
307#                fpp += BL[8]/D
308#            FP.append(BL[2]*fp)
309#            FPP.append(-BL[3]*fpp)
310#        else:
311#            FP.append(0.0)
312#            FPP.append(0.0)
313#    return np.array(FP),np.array(FPP)
314   
315def ComptonFac(El,SQ):
316    """compute Compton scattering factor
317    @param El: element dictionary
318    @param SQ: (sin-theta/lambda)**2
319    @return: compton scattering factor
320    """   
321    ca = np.array(El['cmpa'])
322    cb = np.array(El['cmpb'])
323    t = -cb[:,np.newaxis]*SQ       
324    return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) 
325           
326def FPcalc(Orbs, KEv):
327    """Compute real & imaginary resonant X-ray scattering factors
328    @param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
329    @param KEv: x-ray energy in keV
330    @return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
331    """
332    def Aitken(Orb, LKev):
333        Nval = Orb['Nval']
334        j = Nval-1
335        LEner = Orb['LEner']
336        for i in range(Nval):
337            if LEner[i] <= LKev: j = i
338        if j > Nval-3: j= Nval-3
339        T = [0,0,0,0,0,0]
340        LXSect = Orb['LXSect']
341        for i in range(3):
342           T[i] = LXSect[i+j]
343           T[i+3] = LEner[i+j]-LKev
344        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
345        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
346        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
347        C = T[2]
348        return C
349   
350    def DGauss(Orb,CX,RX,ISig):
351        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
352        0.23931433524968,0.11846344252810)
353        XLG = (0.04691007703067,0.23076534494716,0.5,
354        0.76923465505284,0.95308992296933)
355       
356        D = 0.0
357        B2 = Orb['BB']**2
358        R2 = RX**2
359        XSecIP = Orb['XSecIP']
360        for i in range(5):
361            X = XLG[i]
362            X2 = X**2
363            XS = XSecIP[i]
364            if ISig == 0:
365                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
366            elif ISig == 1:
367                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
368            elif ISig == 2:
369                T = X*X2*R2-B2/X
370                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
371            else:
372                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
373            A = ALG[i]
374            D += A*S
375        return D
376   
377    AU = 2.80022e+7
378    C1 = 0.02721
379    C = 137.0367
380    FP = 0.0
381    FPP = 0.0
382    Mu = 0.0
383    LKev = math.log(KEv)
384    RX = KEv/C1
385    if Orbs:
386        for Orb in Orbs:
387            CX = 0.0
388            BB = Orb['BB']
389            BindEn = Orb['BindEn']
390            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
391            if BindEn <= KEv:
392                CX = math.exp(Aitken(Orb,LKev))
393                Mu += CX
394                CX /= AU
395            Corr = 0.0
396            if Orb['IfBe'] == 0 and BindEn >= KEv:
397                CX = 0.0
398                FPI = DGauss(Orb,CX,RX,3)
399                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
400            else:
401                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
402                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
403            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
404            FPPI = C*CX*RX/(4.0*math.pi)
405            FP += FPI
406            FPP += FPPI
407        FP -= ElEterm
408   
409    return (FP, FPP, Mu)
410   
411
Note: See TracBrowser for help on using the repository browser.