diff --git a/README.rst b/README.rst index 96c8627..235edb0 100644 --- a/README.rst +++ b/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:: diff --git a/cstool/compile.py b/cstool/compile.py index adf2a25..2f63c04 100644 --- a/cstool/compile.py +++ b/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 diff --git a/cstool/inelastic.py b/cstool/inelastic.py index 87085f8..7da242f 100644 --- a/cstool/inelastic.py +++ b/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 } diff --git a/cstool/phonon.py b/cstool/phonon.py index 41319d3..ff04608 100644 --- a/cstool/phonon.py +++ b/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 diff --git a/data/materials/alumina.yaml b/data/materials/alumina.yaml new file mode 100644 index 0000000..d742e09 --- /dev/null +++ b/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 } diff --git a/data/materials/gold.yaml b/data/materials/gold.yaml index d36cf94..7f655d3 100644 --- a/data/materials/gold.yaml +++ b/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 } diff --git a/data/materials/silicon.yaml b/data/materials/silicon.yaml index 9057fcb..a51af4f 100644 --- a/data/materials/silicon.yaml +++ b/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 } diff --git a/examples/cs.py b/examples/cs.py index 33f8f7c..dcf9d22 100644 --- a/examples/cs.py +++ b/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()