source: trunk/GSASIIElem.py @ 881

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

more RB editing & fixes
'+0' & '-0' (vice versa) removed from any cif import atom names
fix spelling of "Sequential"
begin implementation of 'Omit map"

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