Tmfp

README.rst

@@ -27,7 +27,7 @@ clone the needed repositories::
     ...
     > git clone http://github.com/eScatter/pyelsepa.git
     ...
-    > git clone http://github.com/eScatter/cstool.git
+    > git clone -b tmfp http://github.com/eScatter/cstool.git
     ...
 
 create a virtual environment::

cstool/compile.py

@@ -147,3 +147,11 @@ def compute_tcs_icdf(f, a, b, P, sampling=100000):
     if cf[-1]==0:
         return 0 * x.units*y.units, np.zeros_like(P) * x.units
     return cf[-1] * x.units*y.units, np.interp(P, cf/cf[-1], x) * x.units
+
+def compute_tcs(f, a, b, P, sampling=100000):
+    x = np.linspace(a, b.to(a.units), sampling) * a.units
+    y = f(x)
+    cf = np.r_[0, np.cumsum((x[1:] - x[:-1]) * (y[:-1] + y[1:]) / 2.0)]
+    if cf[-1]==0:
+        return 0 * x.units*y.units
+    return cf[-1] * x.units*y.units

cstool/inelastic.py

@@ -39,6 +39,37 @@ def L_Kieft(K, w0, F):
     return np.maximum(0, (w0 < 50 * units.eV) * L1
                       + (w0 > 50 * units.eV) * L2)
 
+def L_dv1(K, w0, F):
+    """Computes electron cross-sections for inelastic scattering from
+    optical data. Model is conceptually somewhere between L_Kieft and L_Ashley:
+    L1 is a Fermi-corrected version of Ashley without the factor 3/2
+    rescale by Kieft; L2 is the same as in Kieft.
+
+    :param K:
+        Kinetic energy of electron.
+    :param w0:
+        ω₀ - transition energy
+    :param F:
+        Fermi energy
+    """
+
+    a = (w0 / K).magnitude
+    b = (F / K).magnitude
+    s = sqrt(1 - 2*a, where = (a <= .5), out = np.zeros(a.shape))
+
+    L1_range = (a > 0) * (a < .5) * (a - s < 1 - 2*b)
+    L2_range = (a > 0) * (a < 1 - b)
+
+    # Calculate L1
+    wm = (1 + a - s)/2
+    wp = np.minimum((1 + a + s)/2, 1 - b)
+    L1 = log((wp - a) * wm / (wp * (wm - a)), where = L1_range, out = np.zeros(a.shape))
+
+    # Calculate L2
+    L2 = -log(a, where = L2_range, out = np.zeros(a.shape))
+
+    return np.maximum(0, (w0 < 50 * units.eV) * L1
+                      + (w0 > 50 * units.eV) * L2)
 
 def L_Ashley_w_ex(K, w0, _):
     a = w0 / K
@@ -54,7 +85,8 @@ def L_Ashley_wo_ex(K, w0, _):
 methods = {
     'Kieft': L_Kieft,
     'Ashley_w_ex': L_Ashley_w_ex,
-    'Ashley_wo_ex': L_Ashley_wo_ex
+    'Ashley_wo_ex': L_Ashley_wo_ex,
+    'dv1': L_dv1
 }
 
 

cstool/phonon.py

@@ -103,7 +103,7 @@ def dcs(theta, E):
             lambda E: 1, partial(dcs_hi, m),
             h, E_BZ / 4, E_BZ)
 
-        return g(E) * norm(m, E)
+        return g(E) * norm(m, E) * 3.0
 
     # should have units of m²/sr
     return dcs

data/materials/alumina.yaml

@@ -0,0 +1,55 @@
+name:           alumina                 # http://www.webelements.com (also for density,
+                                        # atomic mass values for all elements)
+rho_m:          3.98 g/cm³              # http://accuratus.com/alumox.html this
+                                        # link shows the density of alumina increases
+                                        # with the purity level
+                                        # http://www.azom.com/properties.aspx?ArticleID=52
+                                        # this link gives a range of densities for alumina
+                                        # from 3 to 3.98 g/cm³
+                                        # http://www-ferp.ucsd.edu/LIB/PROPS/PANOS/al2o3.html
+                                        # this site has a density of 3.9 g/cm³
+                                        # and list a reference: Goodfellow Cambridge Ltd.,
+                                        # "Metals, Alloys, Compounds, Ceramics,
+                                        # Polymers, Composites", Catalogue 1993/94.
+elf_file:       data/elf/df_Al2O3.dat
+
+band_structure:
+    model:      insulator
+    fermi:      0.0 eV                  # E.O. Filatova, A.S. Konashuk,
+                                        # DOI: 10.1021/acs.jpcc.5b06843 J.Phys.
+                                        # Chem.C 2015, 119, 20755 − 20761
+                                        # gives the valence band max as
+                                        # 3.64 +- 0.04 eV for am-Al2O3
+    band_gap:   7.0 eV                  # E.O. Filatova, A.S. Konashuk,
+                                        # DOI: 10.1021/acs.jpcc.5b06843 J.Phys.
+                                        # Chem.C 2015, 119, 20755 − 20761
+                                        # gives the band gap as
+                                        # 7.0 +- 0.1 eV for am-Al2O3
+                                        # and 7.6 +- 0.1 eV for gamma-Al2O3
+    affinity:   1.0 eV                  # (unrealistic value?) D.V. Morgan et al.,
+                                        # J.Phys.D 13, 307 (1980).
+
+phonon:
+    model:      single                  # dual parameters for ac. def. potentail are yet unknown
+    lattice:    4.76 Å                  # http://www.ceramics.nist.gov (does not work anymore)
+                                        # https://srdata.nist.gov/CeramicDataPortal/Pds/Scdaos (for sintered alumina)
+                                        # hexagonal: a = 4.761 Å and c = 12.991 Å
+                                        # Landolt-Bornstein: 5.140 Å and alpha = 55”16’
+    m_dos:      1.0 m_e                 # Density of state mass [] (unknown in KB)
+    m_eff:      1.0 m_e                 # Effective mass (unknown in KB)
+    single:
+        c_s:    8009 m/s                # http://www.ceramics.nist.gov (does not work anymore)
+                                        # https://srdata.nist.gov/CeramicDataPortal/Pds/Scdaos (for sintered alumina)
+        eps_ac: 13.0 eV                 # J. Shan et al., Phys.Rev.Lett. 90(24),
+                                        # 247401 (2003), using model of L.P.
+                                        # Kadanoff, Phys.Rev. 130, 1364 (1963).
+    longitudinal:                       # idem dito for longitudinal,
+        c_s:    11003 m/s               # https://srdata.nist.gov/CeramicDataPortal/Pds/Scdaos (for sintered alumina)
+        eps_ac: 6.39 eV                 # value from silicon
+    transversal:                        # and transversal modes
+        c_s:    6512 m/s                # https://srdata.nist.gov/CeramicDataPortal/Pds/Scdaos (for sintered alumina)
+        eps_ac: 3.01 eV                 # value from silicon
+
+elements:
+    Al: { count: 2, Z: 13, M: 26.982 g/mol }
+    O:  { count: 3, Z:  8, M: 15.999 g/mol }

data/materials/gold.yaml

@@ -1,18 +1,32 @@
 name:           gold
-rho_m:          19.30 g/cm³
+rho_m:          19.30 g/cm³           # http://www.webelements.com (also for density,
+                                      # atomic mass values for all elements)
+                                      # also resistivity for gold (2.2 e-8 Ohm m)
 elf_file:       data/elf/df_Au.dat
 
 band_structure:
     model:      metal
-    fermi:      5.53 eV
-    work_func:  5.4 eV
+    fermi:      5.53 eV               # http://hyperphysics.phy-astr.gsu.edu/hbase/hframe.html,
+                                      # quoting Ashcroft, N. W. and Mermin, N.
+                                      # D., Solid State Physics, Saunders, 1976. (5.53 eV)
+                                      # and Ohya et al., report NIFS-DATA-84, Japan. (9.11 eV)
+    work_func:  5.38 eV               # Ohya et al., report NIFS-DATA-84, Japan. (5.38 eV)
 
 phonon:
-    model:      single
-    lattice:    4.079 Å
+    model:      dual
+    lattice:    4.0782 Å              # https://www.webelements.com/gold/crystal_structure.html
+                                      # Landolt-Bornstein: 4.0786 Å
+    m_dos:      1.0 m_e               # Density of state mass [] (unknown in KB)
+    m_eff:      1.0 m_e               # Effective mass (unknown in KB)
     single:
-        c_s:    2.03 km/s
-        eps_ac: 8.92 eV
+        c_s:    1880 m/s             # speed of sound
+        eps_ac: 2.82 eV               # acoustic deformation potential
+    longitudinal:                     # idem dito for longitudinal,
+        c_s:    3240 m/s              # https://en.wikipedia.org/wiki/Speeds_of_sound_of_the_elements_(data_page)
+        eps_ac: 4.86 eV               # calculated
+    transversal:                      # and transversal modes
+        c_s:    1200 m/s              # https://en.wikipedia.org/wiki/Speeds_of_sound_of_the_elements_(data_page)
+        eps_ac: 1.80 eV               # calculated
 
 elements:
     Au: { count: 1, Z: 79, M: 196.97 g/mol }

data/materials/silicon.yaml

@@ -1,29 +1,45 @@
 name:           silicon
-rho_m:          2.3290 g/cm³
+rho_m:          2.3290 g/cm³            # http://www.webelements.com (also for density,
+                                        # atomic mass values for all elements)
+                                        # and http://www.ioffe.ru/SVA/NSM/Semicond/Si/basic.html
 elf_file:       data/elf/df_Si.dat
 
 band_structure:
     model:      semiconductor
-    fermi:      7.83 eV
-    band_gap:   1.12 eV
-    affinity:   4.05 eV
+    fermi:      7.83 eV                 # Ohya et al., report NIFS-DATA-84, Japan. (and workfunction? 4.79 eV)
+    band_gap:   1.12 eV                 # http://www.ioffe.rssi.ru/SVA/NSM/Semicond/
+                                        # http://www.ioffe.ru/SVA/NSM/Semicond/Si/bandstr.html
+    affinity:   4.05 eV                 # http://www.ioffe.ru/SVA/NSM/Semicond/Si/basic.html
 
 phonon:
     model:      dual
-    lattice:    5.430710 Å
-    m_dos:      1.09 m_e                # Density of state mass []
+    lattice:    5.430710 Å              # Several volumes in Landolt-Börnstein,
+                                        # book series, Group III Condensed Matter,
+                                        # Springer-Verlag. (5.43072 Å)
+    m_dos:      1.08 m_e                # Density of state mass []
+                                        # https://ecee.colorado.edu/~bart/book/effmass.htm
+                                        # http://onlinelibrary.wiley.com/doi/10.1002/9780470769522.app2/pdf
     m_eff:      0.26 m_e                # Effective mass
+                                        # https://ecee.colorado.edu/~bart/book/effmass.htm
     single:
-        c_s:    6490 m/s
-        eps_ac: 4.14 eV
-    longitudinal:                       # idem dito for longitudinal,
-        alpha:  2.00e-7 m²/s
-        c_s:    9010 m/s
-        eps_ac: 6.39 eV
-    transversal:                        # and transversal modes
-        alpha:  2.26e-7 m²/s
-        c_s:    5230 m/s
-        eps_ac: 3.01 eV
+        c_s:    6938 m/s                # Speed of sound
+                                        # weighted average from dual c_s's
+        eps_ac: 6.4 eV                  # Acoustic deformation potential
+                                        # weighted average from dual eps_ac's
+    longitudinal:                       # parameters for longitudinal branch,
+        alpha:  2.00e-7 m²/s            # thesis T. Verduin
+        c_s:    9130 m/s                # Landolt-Bornstein Vol. III/41A1a 872
+                                        # Smirnov: 9033 m/s
+                                        # http://www.iue.tuwien.ac.at/phd/smirnov/
+        eps_ac: 9.2 eV                  # Landolt-Bornstein Vol. III/41A1a 648
+                                        # Smirnov (see link above): 11.0 eV
+                                        # see also http://onlinelibrary.wiley.com/doi/10.1002/pssb.2220430167/pdf
+    transversal:                        # and transversal branch parameters
+        alpha:  2.26e-7 m²/s            # thesis T. Verduin
+        c_s:    5842 m/s                # Landolt-Bornstein Vol. III/41A1a 872
+                                        # Smirnov (see link above): 5410 m/s
+        eps_ac: 5.0 eV                  # Landolt-Bornstein Vol. III/41A1a 648
+                                        # Smirnov (see link above): 7.2 eV
 
 elements:
     Si: { count: 1, Z: 14, M: 28.0855 g/mol }

examples/cs.py

@@ -6,6 +6,7 @@
 from cstool.phonon import phonon_cs_fn
 from cstool.inelastic import inelastic_cs_fn
 from cstool.compile import compute_tcs_icdf
+from cstool.compile import compute_tcs
 from cstool.ionization import ionization_shells, outer_shell_energies, \
                               loglog_interpolate as ion_loglog_interp
 from cslib.noodles import registry
@@ -23,6 +24,12 @@ def integrant(theta):
 
     return compute_tcs_icdf(integrant, 0*units('rad'), np.pi*units('rad'), P)
 
+def compute_elastic_tcs(dcs, P):
+    def integrant(theta):
+        return dcs(theta) * 2 * np.pi * np.sin(theta)
+
+    return compute_tcs(integrant, 0*units('rad'), np.pi*units('rad'), P)
+
 
 def compute_inelastic_tcs_icdf(dcs, P, K0, K, max_interval):
     def integrant(w):
@@ -33,6 +40,15 @@ def integrant(w):
             int(np.ceil((K - K0) / max_interval))
             ]))
 
+def compute_inelastic_tcs(dcs, P, K0, K, max_interval):
+    def integrant(w):
+        return dcs(w)
+
+    return compute_tcs(integrant, K0, K, P,
+        sampling = np.min([100000,
+            int(np.ceil((K - K0) / max_interval))
+            ]))
+
 
 if __name__ == "__main__":
     import argparse
@@ -193,4 +209,74 @@ def dcs(w):
     ionization_osi = ionization_grp.create_dataset("outer_shells", data=s.elf_file.get_outer_shells().to('eV'))
     ionization_osi.attrs['units'] = 'eV'
 
+    # electron range
+    e_ran = np.logspace(-2, 3, 129) * units.eV
+    p_ran = np.linspace(0.0, 1.0, 1024)
+
+    electron_range_grp = outfile.create_group("electron_range")
+    ran_energies = electron_range_grp.create_dataset("energy", data=e_ran.to('eV'))
+    ran_energies.attrs['units'] = 'eV'
+    ran_range = electron_range_grp.create_dataset("range", e_ran.shape)
+    ran_range.attrs['units'] = 'm'
+
+    ran_tmfp = np.zeros(e_ran.shape) * units.m;
+    tmfpmax = 0. * units.m
+    print("# Computing transport mean free path.")
+    for i, K in enumerate(e_ran):
+        def dcs(theta):
+            return elastic_cs_fn(theta, K) * (1 - np.cos(theta)) * s.rho_n
+        inv_tmfp = compute_elastic_tcs(dcs, p_ran)
+        tmfp = 1./inv_tmfp
+        if ( K >= s.band_structure.barrier):
+            # make sure the tmfp is monotonously increasing for energies above
+            # the vacuum potential barrier
+            if (tmfp < tmfpmax):
+                tmfp = tmfpmax
+            else:
+                tmfpmax = tmfp
+        else:
+            tmfpmax = tmfp
+        ran_tmfp[i] = tmfp
+        print('.', end='', flush=True)
+    print()
+
+    rangemax = 0. * units.m
+    tmprange = np.zeros(e_ran.shape) * units.m
+    print("# Computing inelastic electron range.")
+    for i, K in enumerate(e_ran):
+        ran_energies[i] = K.to('eV')
+        if (K < s.band_structure.barrier):
+            ran_range[i] = 0. * units.m
+        else:
+            w0_max = K-s.band_structure.fermi # it is not possible to lose so
+            # much energy that the primary electron ends up below the Fermi
+            # level in an inelastic scattering event
+
+            def dcs(w):
+                return inelastic_cs_fn(s)(K, w) * s.rho_n
+            tcs = compute_inelastic_tcs(dcs, p_inel,
+                s.elf_file.get_min_energy(), w0_max,
+                s.elf_file.get_min_energy_interval())
+            tmprange[i] = 1/tcs
+            for j, omega in enumerate(e_ran):
+                if (omega < K/2 and omega < K-s.band_structure.fermi):
+                    if j==0:
+                        omega_old = 0. * units.eV
+                    else:
+                        omega_old = e_ran[j-1]
+                    cdcsvalue = inelastic_cs_fn(s)(K,omega) * s.rho_n * (omega - omega_old)
+                    index = sum(e_ran <= (K - omega)) - 1 # get last energy index smaller
+                    # than the energy remaining after one energy loss event
+                    tmprange[i] += tmprange[index] * cdcsvalue / tcs
+                else:
+                    break
+            rangevalue = tmprange[i]
+            if (2 * ran_tmfp[i] < tmprange[i]):
+                rangevalue = np.sqrt(2 * ran_tmfp[i] * tmprange[i])
+            if (rangevalue > rangemax):
+                rangemax = rangevalue
+            ran_range[i] = rangemax.to('m')
+        print('.', end='', flush=True)
+    print()
+
     outfile.close()