From ec98b529c4697ca6d7f4211c2bb45540cffc73ae Mon Sep 17 00:00:00 2001 From: RemiHelleboid Date: Sun, 11 Feb 2024 15:58:01 +0100 Subject: [PATCH] Start implementing el-ph --- apps/CMakeLists.txt | 6 + apps/elph.cpp | 70 ++++++++ apps/epsilon.cpp | 1 + parameter_files/phonon_michaillat.yaml | 42 +++++ src/BZ_MESH/CMakeLists.txt | 3 + src/BZ_MESH/electron_phonon.cpp | 205 ++++++++++++++++++++++ src/BZ_MESH/electron_phonon.hpp | 233 +++++++++++++++++++++++++ src/BZ_MESH/mesh_tetra.cpp | 6 + src/BZ_MESH/mesh_tetra.hpp | 20 ++- src/BZ_MESH/mesh_vertex.hpp | 26 +++ src/BZ_MESH/vector.hpp | 2 + src/EPP/Constants.hpp | 14 +- src/EPP/Material.h | 2 +- src/EPP/Vector3.cpp | 0 src/EPP/Vector3D.h | 9 +- 15 files changed, 626 insertions(+), 13 deletions(-) create mode 100644 apps/elph.cpp create mode 100644 parameter_files/phonon_michaillat.yaml create mode 100644 src/BZ_MESH/electron_phonon.cpp create mode 100644 src/BZ_MESH/electron_phonon.hpp create mode 100644 src/EPP/Vector3.cpp diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt index b140722..81d8a6f 100644 --- a/apps/CMakeLists.txt +++ b/apps/CMakeLists.txt @@ -21,6 +21,11 @@ target_link_libraries(fullstates.epm PUBLIC libepp lib_bzmesh) target_compile_features(fullstates.epm PRIVATE cxx_std_20) target_link_libraries(fullstates.epm PUBLIC MPI::MPI_CXX) +add_executable(elph.epm elph.cpp) +target_link_libraries(elph.epm PUBLIC libepp lib_bzmesh) +target_compile_features(elph.epm PRIVATE cxx_std_20) +target_link_libraries(elph.epm PUBLIC MPI::MPI_CXX) + if(USE_MPI_ACCELERATION) add_executable(mpiBandsOnBZ mpi_BandsOnBZ.cpp) target_link_libraries(mpiBandsOnBZ PUBLIC libepp Eigen3::Eigen) @@ -32,3 +37,4 @@ if(USE_MPI_ACCELERATION) target_link_libraries(mpiDOS_MeshBZ PUBLIC MPI::MPI_CXX) target_compile_features(mpiDOS_MeshBZ PRIVATE cxx_std_20) endif(USE_MPI_ACCELERATION) + diff --git a/apps/elph.cpp b/apps/elph.cpp new file mode 100644 index 0000000..4e46ba3 --- /dev/null +++ b/apps/elph.cpp @@ -0,0 +1,70 @@ +/** + * @file electron_phonon.cpp + * @author your name (you@domain.com) + * @brief + * @version 0.1 + * @date 2024-02-10 + * + * @copyright Copyright (c) 2024 + * + */ + + +#include + +#include +#include +#include +#include + +#include "BandStructure.h" +#include "Material.h" +#include "Options.h" +#include "bz_mesh.hpp" +#include "bz_meshfile.hpp" +#include "bz_states.hpp" +#include "electron_phonon.hpp" + + +int main(int argc, char const *argv[]) +{ + TCLAP::CmdLine cmd("EPP PROGRAM. COMPUTE BAND STRUCTURE ON A BZ MESH.", ' ', "1.0"); + TCLAP::ValueArg arg_mesh_file("f", "meshbandfile", "File with BZ mesh and bands energy.", true, "bz.msh", "string"); + TCLAP::ValueArg arg_material("m", "material", "Symbol of the material to use (Si, Ge, GaAs, ...)", true, "Si", "string"); + TCLAP::ValueArg arg_nb_energies("e", "nenergy", "Number of energies to compute", false, 250, "int"); + TCLAP::ValueArg arg_nb_bands("b", "nbands", "Number of bands to consider", false, 12, "int"); + TCLAP::ValueArg arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int"); + TCLAP::SwitchArg plot_with_python("P", "plot", "Call a python script after the computation to plot the band structure.", false); + cmd.add(plot_with_python); + cmd.add(arg_mesh_file); + cmd.add(arg_material); + cmd.add(arg_nb_bands); + cmd.add(arg_nb_energies); + cmd.add(arg_nb_threads); + + cmd.parse(argc, argv); + + EmpiricalPseudopotential::Materials materials; + const std::string file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials.yaml"; + materials.load_material_parameters(file_material_parameters); + + Options my_options; + my_options.materialName = arg_material.getValue(); + my_options.nrLevels = arg_nb_bands.getValue(); + my_options.nrThreads = arg_nb_threads.getValue(); + int number_energies = arg_nb_energies.getValue(); + + EmpiricalPseudopotential::Material current_material = materials.materials.at(arg_material.getValue()); + + const std::string mesh_band_input_file = arg_mesh_file.getValue(); + bz_mesh::ElectronPhonon ElectronPhonon{current_material}; + const std::string phonon_file = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/phonon_michaillat.yaml"; + ElectronPhonon.load_phonon_parameters(phonon_file); + + + ElectronPhonon.read_mesh_geometry_from_msh_file(mesh_band_input_file); + ElectronPhonon.read_mesh_bands_from_msh_file(mesh_band_input_file); + + // ElectronPhonon.compute_electron_phonon_rates_over_mesh(); + +} diff --git a/apps/epsilon.cpp b/apps/epsilon.cpp index 67e6607..616c097 100644 --- a/apps/epsilon.cpp +++ b/apps/epsilon.cpp @@ -26,6 +26,7 @@ #define MASTER 0 #include + #include #include diff --git a/parameter_files/phonon_michaillat.yaml b/parameter_files/phonon_michaillat.yaml new file mode 100644 index 0000000..22d9d56 --- /dev/null +++ b/parameter_files/phonon_michaillat.yaml @@ -0,0 +1,42 @@ +# Materials parameters for the electron-phonon interaction. +# All the parameters are from M. Michaillat PhD Manuscript. + +materials: + - name: Si + Radius-WS: 2.122 + dispersion: + longitudinal: + acoustic: + w0: 0.0 + vs: 9.01E3 + c: -2.00e-7 + optic: + w0: 9.88e13 + vs: 0.0 + c: -1.60e-7 + transverse: + acoustic: + w0: 0.0 + vs: 5.23e3 + c: -2.26e-7 + optic: + w0: 1.02e14 + vs: -2.57e3 + c: 1.11e-7 + deformation-potential: + electron: + acoustic: + A: 6.0 + B: -2.0 + optic: + A: 2.5e17 + B: -8.0e16 + energy-threshold: 2.0 + hole: + acoustic: + A: 8.0 + B: -3.0 + optic: + A: 2.0e17 + B: -9.0e16 + energy-threshold: 1.5 diff --git a/src/BZ_MESH/CMakeLists.txt b/src/BZ_MESH/CMakeLists.txt index 4a5289f..ab46b9e 100644 --- a/src/BZ_MESH/CMakeLists.txt +++ b/src/BZ_MESH/CMakeLists.txt @@ -5,12 +5,15 @@ set(HEADER_FILES_LIBBZMESH vector.hpp iso_triangle.hpp bz_states.hpp + electron_phonon.hpp ) + set(SOURCE_FILES_LIBBZMESH mesh_tetra.cpp bz_mesh.cpp bz_states.cpp + electron_phonon.cpp ) add_library(lib_bzmesh STATIC ${SOURCE_FILES_LIBBZMESH} ${HEADER_FILES_LIBBZMESH}) diff --git a/src/BZ_MESH/electron_phonon.cpp b/src/BZ_MESH/electron_phonon.cpp new file mode 100644 index 0000000..81b6f05 --- /dev/null +++ b/src/BZ_MESH/electron_phonon.cpp @@ -0,0 +1,205 @@ +/** + * @file elelectron_phonon.cpp + * @author your name (you@domain.com) + * @brief + * @version 0.1 + * @date 2024-02-09 + * + * @copyright Copyright (c) 2024 + * + */ + +#include "electron_phonon.hpp" + +#include +#include +#include +#include +#include + +#include "Constants.hpp" +#include "Vector3D.h" +#include "bz_states.hpp" + +namespace bz_mesh { + +double ElectronPhonon::bose_einstein_distribution(double energy, double temperature) { + double f = 1.0 / (std::expm1(energy / (EmpiricalPseudopotential::Constants::k_b_eV * temperature))); + return f; +} + +double ElectronPhonon::electron_overlap_integral(const Vector3D& k1, const Vector3D& k2) { + const double R_Wigner_Seitz = 2.122e-10; + const double qRws = (k1 - k2).Length() * R_Wigner_Seitz; + double integral = 3.0 * (std::sin(qRws) - qRws * std::cos(qRws)) / (qRws * qRws * qRws); + return integral; +} + +double ElectronPhonon::hole_overlap_integral(int n1, const Vector3D& k1, int n2, const Vector3D& k2) { + const double cos_angle_k1_k2_2 = compte_cos_angle(k1, k2) * compte_cos_angle(k1, k2); + auto A_B_params = m_hole_overlap_int_params.get_params(n1, n2); + double integral = (1.0 / 2.0) * std::sqrt(A_B_params[0] + A_B_params[1] * cos_angle_k1_k2_2); + return integral; +} + +RateValues ElectronPhonon::compute_electron_phonon_rate(int idx_n1, std::size_t idx_k1) { + RateValues rates_k1_n1; + const std::vector& list_tetrahedra = m_list_tetrahedra; + auto indices_conduction_bands = m_indices_conduction_bands; + double energy_n1_k1 = m_list_vertices[idx_k1].get_energy_at_band(idx_n1); + + for (auto&& idx_n2 : indices_conduction_bands) { + for (auto&& tetra : list_tetrahedra) { + auto k_1 = m_list_vertices[idx_k1].get_position(); + auto k_2 = tetra.compute_barycenter(); + Vector3D k1{k_1.x(), k_1.y(), k_1.z()}; + Vector3D k2{k_2.x(), k_2.y(), k_2.z()}; + + double overlap_integral = electron_overlap_integral(k1, k2); + auto q = k1 - k2; + Vector3D q_ph{q.X, q.X, q.X}; + // for (auto&& mode_phonon : m_phonon_dispersion) { + // double e_ph = mode_phonon.get_phonon_dispersion(q_ph.Length()); + // for (auto sign_phonon : {-1.0, 1.0}) { + // double bose_part = (bose_einstein_distribution(e_ph, m_temperature) + sign_phonon); + // double energy_final = energy_n1_k1 + e_ph * sign_phonon; + // double dos_tetra = tetra.compute_tetra_dos_energy_band(energy_final, idx_n2); + // DeformationPotential deformation_potential = + // (mode_phonon.m_mode == PhononMode::acoustic) ? m_acoustic_deformation_potential : + // m_optical_deformation_potential; + // double deformation_potential_value = deformation_potential.get_deformation_potential(q_ph, energy_final); + // PhononMode phonon_mode = mode_phonon.m_mode; + // PhononDirection phonon_direction = mode_phonon.m_direction; + // PhononEvent phonon_event = (sign_phonon == 1.0) ? PhononEvent::absorption : PhononEvent::emission; + // double rate_value = (EmpiricalPseudopotential::Constants::pi / (m_rho * e_ph)) * deformation_potential_value + // * + // deformation_potential_value * overlap_integral * overlap_integral * bose_part * dos_tetra; + // RateValue rate(phonon_mode, phonon_direction, phonon_event, rate_value); + // rates_k1_n1.add_rate(rate); + // } + // } + for (auto&& mode_phonon : m_phonon_dispersion) { + PhononModeDirection mode_direction = mode_phonon.first; + double e_ph = mode_phonon.second.get_phonon_dispersion(q_ph.Length()); + for (auto sign_phonon : {-1.0, 1.0}) { + double bose_part = (bose_einstein_distribution(e_ph, m_temperature) + sign_phonon); + double energy_final = energy_n1_k1 + e_ph * sign_phonon; + double dos_tetra = tetra.compute_tetra_dos_energy_band(energy_final, idx_n2); + DeformationPotential deformation_potential = + (mode_direction.first == PhononMode::acoustic) ? m_acoustic_deformation_potential : m_optical_deformation_potential; + double deformation_potential_value = deformation_potential.get_deformation_potential(q_ph, energy_final); + PhononMode phonon_mode = mode_direction.first; + PhononDirection phonon_direction = mode_direction.second; + PhononEvent phonon_event = (sign_phonon == 1.0) ? PhononEvent::absorption : PhononEvent::emission; + double rate_value = (EmpiricalPseudopotential::Constants::pi / (m_rho * e_ph)) * deformation_potential_value * + deformation_potential_value * overlap_integral * overlap_integral * bose_part * dos_tetra; + RateValue rate(phonon_mode, phonon_direction, phonon_event, rate_value); + rates_k1_n1.add_rate(rate); + } + } + } + } + rates_k1_n1.print_rates(); + return rates_k1_n1; +} + +void ElectronPhonon::compute_electron_phonon_rates_over_mesh() { + for (std::size_t idx_k1 = 0; idx_k1 < m_list_vertices.size(); ++idx_k1) { + std::cout << "Computing rates for k-point " << idx_k1 << std::endl; + for (std::size_t idx_n1 = 0; idx_n1 < m_indices_valence_bands.size(); ++idx_n1) { + auto rate = compute_electron_phonon_rate(idx_n1, idx_k1); + m_list_vertices[idx_k1].add_electron_phonon_rates(rate.to_array()); + } + } + + // Add rate to the mesh +} + +void ElectronPhonon::compute_plot_electron_phonon_rates_vs_energy_over_mesh(int nb_bands, + double max_energy, + double energy_step, + const std::string& filename) { + std::vector energies; + std::vector> rates; + for (double energy = 0.0; energy < max_energy; energy += energy_step) { + energies.push_back(energy); + for (int idx_band = 0; idx_band < nb_bands; ++idx_band) { + double rate = 0.0; + for (auto&& vertex : m_list_vertices) { + for (std::size_t idx_band = 0; idx_band < vertex.get_number_bands(); ++idx_band) { + rate += vertex.get_energy_at_band(idx_band); + } + } + // rates[idx_band].push_back(rate); + } + } +} + +void ElectronPhonon::load_phonon_parameters(const std::string& filename) { + YAML::Node config = YAML::LoadFile(filename); + + // Check if the file is empty + if (config.IsNull()) { + throw std::runtime_error("File " + filename + " is empty"); + } + // Print the file + std::cout << "File " << filename << " contains: " << std::endl; + std::cout << config << std::endl; + + auto list_materials = config["materials"]; + + const std::string& my_material = m_material.get_name(); + + auto same_material = [my_material](const YAML::Node& node) { return node["name"].as() == my_material; }; + auto it_material = std::find_if(list_materials.begin(), list_materials.end(), same_material); + + + + if (it_material == list_materials.end()) { + throw std::runtime_error("Material " + my_material + " not found in file " + filename); + } + auto material = *it_material; + + double radiusWS = material["Radius-WS"].as(); + m_radii_wigner_seitz = radiusWS; + + // Example for reading nested dispersion properties + auto dispersion = material["dispersion"]; + for (const auto& type : {"longitudinal", "transverse"}) { + auto dispType = dispersion[type]; + for (const auto& wave : {"acoustic", "optic"}) { + auto waveType = dispType[wave]; + double w0 = waveType["w0"].as(); + double vs = waveType["vs"].as(); + double c = waveType["c"].as(); + std::cout << "w0: " << w0 << " vs: " << vs << " c: " << c << std::endl; + + PhononDirection direction = (type == "longitudinal") ? PhononDirection::longitudinal : PhononDirection::transverse; + PhononMode mode = (wave == "acoustic") ? PhononMode::acoustic : PhononMode::optical; + PhononDispersion phononDispersion(mode, direction, w0, vs, c); + m_phonon_dispersion[std::make_pair(mode, direction)] = phononDispersion; + } + } + + // Example for reading deformation potential + auto node_deformationPotential = material["deformation-potential"]; + for (const auto& carrierType : {"electron", "hole"}) { + auto carrier = node_deformationPotential[carrierType]; + for (const auto& wave : {"acoustic", "optic"}) { + auto waveType = carrier[wave]; + double A = waveType["A"].as(); + double B = waveType["B"].as(); + std::cout << "A: " << A << " B: " << B << std::endl; + + PhononMode mode = (wave == "acoustic") ? PhononMode::acoustic : PhononMode::optical; + DeformationPotential deformationPotential(mode, A, B); + if (mode == PhononMode::acoustic) { + m_acoustic_deformation_potential = deformationPotential; + } else { + m_optical_deformation_potential = deformationPotential; + } + } + } +} + +} // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/electron_phonon.hpp b/src/BZ_MESH/electron_phonon.hpp new file mode 100644 index 0000000..a3a5d64 --- /dev/null +++ b/src/BZ_MESH/electron_phonon.hpp @@ -0,0 +1,233 @@ +/** + * @file elelectron_phonon.hpp + * @author your name (you@domain.com) + * @brief + * @version 0.1 + * @date 2024-02-09 + * + * @copyright Copyright (c) 2024 + * + * Phonon are always stored in the same order: + * 0: acoustic longitudinal absorption + * 1: acoustic transverse absorption + * 2: optical longitudinal absorption + * 3: optical transverse absorption + * + */ + +#pragma once + +#include + +#include +#include +#include +#include +#include + +#include "bz_states.hpp" +#include "yaml-cpp/yaml.h" + +namespace bz_mesh { + +enum class PhononMode { acoustic, optical, none }; +enum class PhononDirection { longitudinal, transverse, none }; +enum class PhononEvent { absorption, emission, none }; + +struct PhononDispersion { + PhononMode m_mode = PhononMode::none; + PhononDirection m_direction = PhononDirection::none; + double m_omega = 0.0; + double m_vs = 0.0; + double m_c = 0.0; + + PhononDispersion() = default; + + PhononDispersion(PhononMode mode, PhononDirection direction, double omega, double vs, double c) + : m_mode(mode), + m_direction(direction), + m_omega(omega), + m_vs(vs), + m_c(c) {} + PhononDispersion(PhononMode mode, PhononDirection direction) + : m_mode(mode), + m_direction(direction), + m_omega(0.0), + m_vs(0.0), + m_c(0.0) {} + + double get_phonon_dispersion(double q) const { + double e_ph = m_omega + m_vs * q + m_c * q * q; + return e_ph; + } +}; + +struct DeformationPotential { + PhononMode m_mode = PhononMode::none; + double m_A = 0.0; + double m_B = 0.0; + double m_energy_threshold = 0.0; + + DeformationPotential() = default; + DeformationPotential(PhononMode type, double A, double B) : m_mode(type), m_A(A), m_B(B) {} + + double get_deformation_potential(Vector3D q, double energy) const { + double clamp_energy = std::max(energy, m_energy_threshold); + if (m_mode == PhononMode::acoustic) { + return std::sqrt(m_A + clamp_energy * m_B) * q.Length(); + } else { + return std::sqrt(m_A + clamp_energy * m_B); + } + } +}; + +struct HoleOverlapIntParams { + std::vector, double>> A_params = + {{{1, 1}, 1.0}, {{2, 2}, 1.0}, {{3, 3}, 5.0 / 8.0}, {{1, 2}, 3.0}, {{1, 3}, 3.0 / 8.0}, {{2, 3}, 3.0 / 8.0}}; + + std::vector, double>> B_params = + {{{1, 1}, 3.0}, {{2, 2}, 3.0}, {{3, 3}, 0.0}, {{1, 2}, -3.0}, {{1, 3}, 0.0}, {{2, 3}, 0.0}}; + + /** + * @brief Return the parameters for the overlap integral, given the indices of the bands. + * + * @param n1 + * @param n2 + * @return std::array + */ + std::array get_params(int n1, int n2) const { + for (const auto& p : A_params) { + if (p.first.first == n1 && p.first.second == n2) { + return {p.second, 0.0}; + } + } + for (const auto& p : B_params) { + if (p.first.first == n1 && p.first.second == n2) { + return {0.0, p.second}; + } + } + return {0.0, 0.0}; + } + + HoleOverlapIntParams() = default; +}; + +struct RateValue { + PhononMode m_mode = PhononMode::none; + PhononDirection m_direction = PhononDirection::none; + PhononEvent m_type = PhononEvent::none; + double m_value = 0.0; + + RateValue(PhononMode mode, PhononDirection direction, PhononEvent type, double value) + : m_mode(mode), + m_direction(direction), + m_type(type), + m_value(value) {} +}; + +struct RateValues { + RateValue m_ac_long_absorption = RateValue(PhononMode::acoustic, PhononDirection::longitudinal, PhononEvent::absorption, 0.0); + RateValue m_ac_trans_absorption = RateValue(PhononMode::acoustic, PhononDirection::transverse, PhononEvent::absorption, 0.0); + RateValue m_opt_long_absorption = RateValue(PhononMode::optical, PhononDirection::longitudinal, PhononEvent::absorption, 0.0); + RateValue m_opt_trans_absorption = RateValue(PhononMode::optical, PhononDirection::transverse, PhononEvent::absorption, 0.0); + + RateValue m_ac_long_emission = RateValue(PhononMode::acoustic, PhononDirection::longitudinal, PhononEvent::emission, 0.0); + RateValue m_ac_trans_emission = RateValue(PhononMode::acoustic, PhononDirection::transverse, PhononEvent::emission, 0.0); + RateValue m_opt_long_emission = RateValue(PhononMode::optical, PhononDirection::longitudinal, PhononEvent::emission, 0.0); + RateValue m_opt_trans_emission = RateValue(PhononMode::optical, PhononDirection::transverse, PhononEvent::emission, 0.0); + + /** + * @brief Print the rates to the standard output. + * + */ + void print_rates() const { + std::cout << "Acoustic longitudinal absorption: " << m_ac_long_absorption.m_value << std::endl; + std::cout << "Acoustic transverse absorption: " << m_ac_trans_absorption.m_value << std::endl; + std::cout << "Optical longitudinal absorption: " << m_opt_long_absorption.m_value << std::endl; + std::cout << "Optical transverse absorption: " << m_opt_trans_absorption.m_value << std::endl; + std::cout << "Acoustic longitudinal emission: " << m_ac_long_emission.m_value << std::endl; + std::cout << "Acoustic transverse emission: " << m_ac_trans_emission.m_value << std::endl; + std::cout << "optical longitudinal emission: " << m_opt_long_emission.m_value << std::endl; + std::cout << "Optical transverse emission: " << m_opt_trans_emission.m_value << std::endl; + } + + std::array to_array() const { + return {m_ac_long_absorption.m_value, + m_ac_trans_absorption.m_value, + m_opt_long_absorption.m_value, + m_opt_trans_absorption.m_value, + m_ac_long_emission.m_value, + m_ac_trans_emission.m_value, + m_opt_long_emission.m_value, + m_opt_trans_emission.m_value}; + } + + void add_rate(RateValue rate) { + if (rate.m_mode == PhononMode::acoustic) { + if (rate.m_direction == PhononDirection::longitudinal) { + if (rate.m_type == PhononEvent::absorption) { + m_ac_long_absorption.m_value += rate.m_value; + } else { + m_ac_long_emission.m_value += rate.m_value; + } + } else { + if (rate.m_type == PhononEvent::absorption) { + m_ac_trans_absorption.m_value += rate.m_value; + } else { + m_ac_trans_emission.m_value += rate.m_value; + } + } + } else { + if (rate.m_direction == PhononDirection::longitudinal) { + if (rate.m_type == PhononEvent::absorption) { + m_opt_long_absorption.m_value += rate.m_value; + } else { + m_opt_long_emission.m_value += rate.m_value; + } + } else { + if (rate.m_type == PhononEvent::absorption) { + m_opt_trans_absorption.m_value += rate.m_value; + } else { + m_opt_trans_emission.m_value += rate.m_value; + } + } + } + } +}; + +typedef std::pair PhononModeDirection; + +class ElectronPhonon : public BZ_States { + private: + double m_temperature = 0.0; + double m_rho = 0.0; + double m_radii_wigner_seitz = 0.0; + HoleOverlapIntParams m_hole_overlap_int_params; + DeformationPotential m_acoustic_deformation_potential; + DeformationPotential m_optical_deformation_potential; + + std::map m_phonon_dispersion; + + public: + explicit ElectronPhonon(const EmpiricalPseudopotential::Material& material) : BZ_States(material) {} + + void load_phonon_parameters(const std::string& filename); + + double bose_einstein_distribution(double energy, double temperature); + double electron_overlap_integral(const Vector3D& k1, const Vector3D& k2); + double hole_overlap_integral(int n1, const Vector3D& k1, int n2, const Vector3D& k2); + + RateValues compute_electron_phonon_rate(int idx_n1, std::size_t idx_k1); + + void set_temperature(double temperature) { m_temperature = temperature; } + void set_density(double rho) { m_rho = rho; } + + void compute_electron_phonon_rates_over_mesh(); + + void compute_plot_electron_phonon_rates_vs_energy_over_mesh(int nb_bands, + double max_energy, + double energy_step, + const std::string& filename); +}; + +} // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/mesh_tetra.cpp b/src/BZ_MESH/mesh_tetra.cpp index 45f5c24..eecd908 100644 --- a/src/BZ_MESH/mesh_tetra.cpp +++ b/src/BZ_MESH/mesh_tetra.cpp @@ -40,6 +40,12 @@ Tetra::Tetra(std::size_t index, const std::array& list_vertices) m_signed_volume = compute_signed_volume(); } +vector3 Tetra::compute_barycenter() const { + return (m_list_vertices[0]->get_position() + m_list_vertices[1]->get_position() + m_list_vertices[2]->get_position() + + m_list_vertices[3]->get_position()) / + 4.0; +} + void Tetra::compute_gradient_energy_at_bands() { m_gradient_energy_per_band.clear(); std::size_t m_nb_bands = m_list_vertices[0]->get_number_bands(); diff --git a/src/BZ_MESH/mesh_tetra.hpp b/src/BZ_MESH/mesh_tetra.hpp index 38b13dd..b65854c 100644 --- a/src/BZ_MESH/mesh_tetra.hpp +++ b/src/BZ_MESH/mesh_tetra.hpp @@ -102,6 +102,7 @@ class Tetra { vector3 compute_edge(std::size_t index_vtx_1, std::size_t index_vtx_2) const; void compute_gradient_energy_at_bands(); + vector3 compute_barycenter() const; bool is_location_inside(const vector3& location) const; std::array compute_barycentric_coordinates(const vector3& location) const; vector3 compute_euclidean_coordinates(const std::array& barycentric_coordinates) const; @@ -113,19 +114,20 @@ class Tetra { double compute_tetra_iso_surface_energy_band(double energy, std::size_t band_index) const; double compute_tetra_dos_energy_band(double energy, std::size_t band_index) const; + + double interpolate_scalar_at_position(const std::array& barycentric_coordinates, const std::vector& scalar_field) const; vector3 interpolate_vector_at_position(const std::array& barycentric_coordinates, const std::vector& vector_field) const; -template - T interpolate_at_position(const std::array& barycentric_coordinates, - const std::vector& field) const { - T interpolated_value = T::Zero(); - for (std::size_t idx_vtx = 0; idx_vtx < 4; ++idx_vtx) { - interpolated_value += barycentric_coordinates[idx_vtx] * field[idx_vtx]; - } - return interpolated_value; - } + template + T interpolate_at_position(const std::array& barycentric_coordinates, const std::vector& field) const { + T interpolated_value = T::Zero(); + for (std::size_t idx_vtx = 0; idx_vtx < 4; ++idx_vtx) { + interpolated_value += barycentric_coordinates[idx_vtx] * field[idx_vtx]; + } + return interpolated_value; + } static void reset_stat_iso_computing() { ms_case_stats = std::vector(5, 0.0); } diff --git a/src/BZ_MESH/mesh_vertex.hpp b/src/BZ_MESH/mesh_vertex.hpp index 31a06b2..8681015 100644 --- a/src/BZ_MESH/mesh_vertex.hpp +++ b/src/BZ_MESH/mesh_vertex.hpp @@ -43,6 +43,14 @@ class Vertex { */ std::vector m_band_energies; + /** + * @brief Electron-phonon rates for each band. + * The rates are stored as an array of 8 values: + * [ALO, ALA, ATO, ATA, ELO, ELA, ETO, ETA] (ALO: absorption longitudinal optical, ALA: absorption longitudinal acoustic, ...) + * + */ + std::vector> m_electron_phonon_rates; + public: /** * @brief Default constructor of a new Vertex object. @@ -124,6 +132,24 @@ class Vertex { * @return double */ double get_energy_at_band(std::size_t band_index) const { return m_band_energies[band_index]; } + + const std::vector& get_band_energies() const { return m_band_energies; } + + /** + * @brief Add the electron-phonon rates for a given band. + * + * @param rates + */ + void add_electron_phonon_rates(const std::array& rates) { m_electron_phonon_rates.push_back(rates); } + + /** + * @brief Get the electron-phonon rates for a given band. + * + * @param band_index + * @return const std::array& + */ + const std::array& get_electron_phonon_rates(std::size_t band_index) const { return m_electron_phonon_rates[band_index]; } + }; } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/vector.hpp b/src/BZ_MESH/vector.hpp index 09f079d..6c12acf 100644 --- a/src/BZ_MESH/vector.hpp +++ b/src/BZ_MESH/vector.hpp @@ -18,6 +18,8 @@ #include "Vector3D.h" + + namespace bz_mesh { enum class permutaion_type { XY, XZ, YZ, XYZ, YZX, ZXY }; diff --git a/src/EPP/Constants.hpp b/src/EPP/Constants.hpp index b822653..0cdce8b 100644 --- a/src/EPP/Constants.hpp +++ b/src/EPP/Constants.hpp @@ -21,13 +21,23 @@ constexpr double Ryd_to_eV = 13.6; constexpr double h = 6.62606957e-34; constexpr double h_bar = 6.62606957e-34 / (2 * M_PI); constexpr double m0 = 9.10938356e-31; -// constexpr double q = 1.602176634e-19; -constexpr double q = 1.6e-19; +constexpr double q = 1.602176634e-19; constexpr double eps_zero = 8.854187e-12; constexpr double k_b = 1.38064852e-23; +constexpr double k_b_eV = 8.617333262145e-5; constexpr double bohr_radius = 5.291e-11; constexpr double angstrom_to_m = 1e-10; constexpr double ryd_to_hartree = 1.0 / 2.0; +constexpr double hartree_to_eV = 27.21138602; +constexpr double hartree_to_ryd = 2.0; +constexpr double hartree_to_J = 4.3597447222071e-18; +constexpr double eV_to_J = 1.602176634e-19; +constexpr double eV_to_ryd = 1.0 / 13.6056980659; +constexpr double eV_to_hartree = 1.0 / 27.21138602; +constexpr double eV_to_kg_m2 = 1.783e-36; +constexpr double eV_to_cm_1 = 8065.54429; +constexpr double eV_to_mol = 9.64853399e4; +constexpr double pi = M_PI; } // namespace Constants diff --git a/src/EPP/Material.h b/src/EPP/Material.h index 9c1d99c..e9f1cdf 100644 --- a/src/EPP/Material.h +++ b/src/EPP/Material.h @@ -14,7 +14,7 @@ namespace EmpiricalPseudopotential { class Material { protected: /** - * @brief Name of the material. (Actually it is the symbole, eg.g., Si) + * @brief Name of the material. (Actually it is the symbol, eg.g., Si) * */ std::string m_name; diff --git a/src/EPP/Vector3.cpp b/src/EPP/Vector3.cpp new file mode 100644 index 0000000..e69de29 diff --git a/src/EPP/Vector3D.h b/src/EPP/Vector3D.h index dc8fd41..afa54c2 100644 --- a/src/EPP/Vector3D.h +++ b/src/EPP/Vector3D.h @@ -4,6 +4,7 @@ #include #include + template class Vector3D { public: @@ -77,6 +78,13 @@ class Vector3D { return (norm_product < 1.0e-13) ? 1.0 : dot_product / (V1.Length() * V2.Length()); } + /** + * @brief print the vector + * + * @param os + * @param v + * @return std::ostream& + */ friend std::ostream& operator<<(std::ostream& os, const Vector3D& v) { os << "(" << v.X << ", " << v.Y << ", " << v.Z << ")"; return os; @@ -103,7 +111,6 @@ bool operator==(const Vector3D& f, const Vector3D& t) { return f.X == t.X && f.Y == t.Y && f.Z == t.Z; } - template bool operator<(const Vector3D& lhs, const Vector3D& rhs) { return std::tie(lhs.X, lhs.Y, lhs.Z) < std::tie(rhs.X, rhs.Y, rhs.Z);