source: trunk/GSASIIElem.py @ 939

Last change on this file since 939 was 939, checked in by toby, 8 years ago

fix & cleanup unit tests; add/change doc strings for sphinx; add all G2 py files to sphinx

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