source: trunk/GSASIIElem.py @ 642

Last change on this file since 642 was 642, checked in by vondreele, 9 years ago

add CheckElement?

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