From ec53eea4bccc61f163be5703b87c30eb16e14754 Mon Sep 17 00:00:00 2001 From: RemiHelleboid Date: Mon, 19 Feb 2024 10:40:08 +0100 Subject: [PATCH] up --- apps/elph.cpp | 2 + apps/mpi_DOS_MeshBZ.cpp | 2 +- src/BZ_MESH/bz_mesh.cpp | 1 - src/BZ_MESH/electron_phonon.cpp | 125 ++++++++++++++++++++++++++++++-- src/BZ_MESH/electron_phonon.hpp | 1 + src/BZ_MESH/mesh_tetra.cpp | 87 ++++++++++++++-------- src/BZ_MESH/mesh_tetra.hpp | 10 +-- 7 files changed, 184 insertions(+), 44 deletions(-) diff --git a/apps/elph.cpp b/apps/elph.cpp index 7ad2099..322cc0c 100644 --- a/apps/elph.cpp +++ b/apps/elph.cpp @@ -81,6 +81,8 @@ int main(int argc, char const *argv[]) ElectronPhonon.compute_plot_electron_phonon_rates_vs_energy_over_mesh(my_options.nrLevels, 10.0, 0.01, "rates_vs_energy.csv"); + ElectronPhonon.add_electron_phonon_rates_to_mesh(mesh_band_input_file, "rates.msh"); + auto stop = std::chrono::high_resolution_clock::now(); std::cout << "Time taken: " << std::chrono::duration_cast(stop - start).count() << " seconds" << std::endl; diff --git a/apps/mpi_DOS_MeshBZ.cpp b/apps/mpi_DOS_MeshBZ.cpp index 8ef5e47..b5c31f9 100644 --- a/apps/mpi_DOS_MeshBZ.cpp +++ b/apps/mpi_DOS_MeshBZ.cpp @@ -261,7 +261,7 @@ int main(int argc, char *argv[]) { const std::string python_plot_dos = std::string(CMAKE_SOURCE_DIR) + "/python/plots/plot_density_of_states.py"; bool call_python_plot = plot_with_python.isSet(); if (call_python_plot) { - std::string python_call = "python3 " + python_plot_dos + " --file " + out_file_bands + ".csv"; + std::string python_call = "nohup python3 " + python_plot_dos + " --file " + out_file_bands + ".csv &"; std::cout << "Executing: " << python_call << std::endl; int succes_plot = system(python_call.c_str()); if (succes_plot != 0) { diff --git a/src/BZ_MESH/bz_mesh.cpp b/src/BZ_MESH/bz_mesh.cpp index 14ef31d..0d80ba1 100644 --- a/src/BZ_MESH/bz_mesh.cpp +++ b/src/BZ_MESH/bz_mesh.cpp @@ -221,7 +221,6 @@ double MeshBZ::compute_dos_at_energy_and_band(double iso_energy, int band_index) for (auto&& tetra : m_list_tetrahedra) { total_dos += tetra.compute_tetra_dos_energy_band(iso_energy, band_index); } - total_dos /= this->m_total_volume; return total_dos; } diff --git a/src/BZ_MESH/electron_phonon.cpp b/src/BZ_MESH/electron_phonon.cpp index 2337b17..7738129 100644 --- a/src/BZ_MESH/electron_phonon.cpp +++ b/src/BZ_MESH/electron_phonon.cpp @@ -23,6 +23,8 @@ #include "Vector3D.h" #include "bz_states.hpp" +#include "gmsh.h" + #include "omp.h" namespace bz_mesh { @@ -103,7 +105,8 @@ RateValues ElectronPhonon::compute_electron_phonon_rate(int idx_n1, std::size_t } // std::cout << "Energy phonon: " << e_ph << std::endl; for (auto sign_phonon : {-1.0, 1.0}) { - double bose_part = (bose_einstein_distribution(e_ph, m_temperature) + 0.5 + 0.5 * std::pow(-1.0, sign_phonon)); + double add_phonon = (sign_phonon == 1.0) ? 1.0 : 0.0; + double bose_part = (bose_einstein_distribution(e_ph, m_temperature) + add_phonon); double energy_final = energy_n1_k1 + e_ph * sign_phonon; // std::cout << "Energy (n1, k1): " << energy_n1_k1 << " Energy phonon: " << e_ph << " Energy final: " << energy_final // << std::endl; @@ -125,7 +128,7 @@ RateValues ElectronPhonon::compute_electron_phonon_rate(int idx_n1, std::size_t deformation_potential_value * overlap_integral * overlap_integral * bose_part * dos_tetra; // std::cout << this->get_volume() << std::endl; rate_value /= this->get_volume(); - rate_value *= EmpiricalPseudopotential::Constants::q; + // rate_value *= EmpiricalPseudopotential::Constants::q; if (rate_value < 0.0 || rate_value > 1e50 || std::isnan(rate_value) || std::isinf(rate_value)) { std::cout << "Rate value: " << rate_value << std::endl; @@ -140,7 +143,7 @@ RateValues ElectronPhonon::compute_electron_phonon_rate(int idx_n1, std::size_t std::cout << "Energy (n1, k1): " << energy_n1_k1 << std::endl; std::cout << "--------------------------------------------------------------------------------" << std::endl; } - // std::cout << "Rate value: " << (deformation_potential_value) << std::endl; + // std::cout << "Rate value: " << rate_value << std::endl; RateValue rate(phonon_mode, phonon_direction, phonon_event, rate_value); rates_k1_n1.add_rate(rate); } @@ -212,8 +215,8 @@ RateValues ElectronPhonon::compute_hole_phonon_rate(int idx_n1, std::size_t idx_ } // std::cout << "Volume: " << this->get_volume() << std::endl; double dos_tetra = tetra.compute_tetra_dos_energy_band(energy_final, idx_n2); - // std::cout << "Energy final: " << energy_final << " dos tetra: " << dos_tetra - // << std::endl; + // double dos_tetra = 1.0; + // std::cout << "DOS: " << dos_tetra << std::endl; DeformationPotential deformation_potential = (mode_direction.first == PhononMode::acoustic) ? m_acoustic_deformation_potential_h : m_optical_deformation_potential_h; @@ -263,7 +266,7 @@ void ElectronPhonon::compute_electron_phonon_rates_over_mesh() { std::random_device rd; std::mt19937 gen(rd()); std::uniform_real_distribution dis(0.0, 1.0); - double p_compute_rate = 1.0; + double p_compute_rate = 0.1; #pragma omp parallel for schedule(dynamic) for (std::size_t idx_k1 = 0; idx_k1 < m_list_vertices.size(); ++idx_k1) { @@ -291,7 +294,8 @@ void ElectronPhonon::compute_electron_phonon_rates_over_mesh() { } } if (omp_get_thread_num() == 0){ - std::cout << "\rComputing rates for k-point " << idx_k1 << std::flush;} + std::cout << "\rComputing rates for k-point " << idx_k1 << " / " << m_list_vertices.size() << std::flush; + } } std::cout << std::endl; @@ -359,6 +363,113 @@ void ElectronPhonon::plot_phonon_dispersion(const std::string& filename) const { } } +void ElectronPhonon::add_electron_phonon_rates_to_mesh(const std::string& initial_filename, + const std::string& final_filename) { + gmsh::initialize(); + gmsh::option::setNumber("General.Verbosity", 0); + gmsh::model::add("bz_mesh"); + gmsh::open(initial_filename); + + std::string model_file_name; + gmsh::model::getCurrent(model_file_name); + + + std::vector node_tags; + std::vector nodeCoords; + std::vector nodeParams; + gmsh::model::mesh::reclassifyNodes(); + gmsh::model::mesh::getNodes(node_tags, nodeCoords, nodeParams, -1, -1, false, false); + + + auto indices_conduction_bands = m_indices_conduction_bands; + auto min_idx_conduction_band = *std::min_element(indices_conduction_bands.begin(), indices_conduction_bands.end()); + auto max_idx_conduction_band = *std::max_element(indices_conduction_bands.begin(), indices_conduction_bands.end()); + std::cout << "Min index conduction band: " << min_idx_conduction_band << std::endl; + int nb_bands = m_indices_conduction_bands.size() + m_indices_valence_bands.size(); + + for (int idx_val_band=0; idx_val_band rates_hole_lo_em(m_list_vertices.size()); + std::vector rates_hole_lo_ab(m_list_vertices.size()); + std::vector rates_hole_tr_em(m_list_vertices.size()); + std::vector rates_hole_tr_ab(m_list_vertices.size()); + std::vector rates_elec_lo_em(m_list_vertices.size()); + std::vector rates_elec_lo_ab(m_list_vertices.size()); + std::vector rates_elec_tr_em(m_list_vertices.size()); + std::vector rates_elec_tr_ab(m_list_vertices.size()); + + for (std::size_t idx_k1 = 0; idx_k1 < m_list_vertices.size(); ++idx_k1) { + auto rates = m_list_vertices[idx_k1].get_electron_phonon_rates(idx_val_band); + rates_hole_lo_em[idx_k1] = rates[0]; + rates_hole_lo_ab[idx_k1] = rates[1]; + rates_hole_tr_em[idx_k1] = rates[2]; + rates_hole_tr_ab[idx_k1] = rates[3]; + rates_elec_lo_em[idx_k1] = rates[4]; + rates_elec_lo_ab[idx_k1] = rates[5]; + rates_elec_tr_em[idx_k1] = rates[6]; + rates_elec_tr_ab[idx_k1] = rates[7]; + } + std::string name_rate_hole_lo_em = "hole_lo_em_" + std::to_string(idx_val_band); + std::string name_rate_hole_lo_ab = "hole_lo_ab_" + std::to_string(idx_val_band); + std::string name_rate_hole_tr_em = "hole_tr_em_" + std::to_string(idx_val_band); + std::string name_rate_hole_tr_ab = "hole_tr_ab_" + std::to_string(idx_val_band); + std::string name_rate_elec_lo_em = "elec_lo_em_" + std::to_string(idx_val_band); + std::string name_rate_elec_lo_ab = "elec_lo_ab_" + std::to_string(idx_val_band); + std::string name_rate_elec_tr_em = "elec_tr_em_" + std::to_string(idx_val_band); + std::string name_rate_elec_tr_ab = "elec_tr_ab_" + std::to_string(idx_val_band); + + int data_tag_hole_lo_em = gmsh::view::add(name_rate_hole_lo_em); + int data_tag_hole_lo_ab = gmsh::view::add(name_rate_hole_lo_ab); + int data_tag_hole_tr_em = gmsh::view::add(name_rate_hole_tr_em); + int data_tag_hole_tr_ab = gmsh::view::add(name_rate_hole_tr_ab); + int data_tag_elec_lo_em = gmsh::view::add(name_rate_elec_lo_em); + int data_tag_elec_lo_ab = gmsh::view::add(name_rate_elec_lo_ab); + int data_tag_elec_tr_em = gmsh::view::add(name_rate_elec_tr_em); + int data_tag_elec_tr_ab = gmsh::view::add(name_rate_elec_tr_ab); + + gmsh::view::addHomogeneousModelData(data_tag_hole_lo_em, 0, model_file_name, "NodeData", node_tags, rates_hole_lo_em); + gmsh::view::addHomogeneousModelData(data_tag_hole_lo_ab, 0, model_file_name, "NodeData", node_tags, rates_hole_lo_ab); + gmsh::view::addHomogeneousModelData(data_tag_hole_tr_em, 0, model_file_name, "NodeData", node_tags, rates_hole_tr_em); + gmsh::view::addHomogeneousModelData(data_tag_hole_tr_ab, 0, model_file_name, "NodeData", node_tags, rates_hole_tr_ab); + gmsh::view::addHomogeneousModelData(data_tag_elec_lo_em, 0, model_file_name, "NodeData", node_tags, rates_elec_lo_em); + gmsh::view::addHomogeneousModelData(data_tag_elec_lo_ab, 0, model_file_name, "NodeData", node_tags, rates_elec_lo_ab); + gmsh::view::addHomogeneousModelData(data_tag_elec_tr_em, 0, model_file_name, "NodeData", node_tags, rates_elec_tr_em); + gmsh::view::addHomogeneousModelData(data_tag_elec_tr_ab, 0, model_file_name, "NodeData", node_tags, rates_elec_tr_ab); + + gmsh::view::write(data_tag_hole_lo_em, final_filename, true); + gmsh::view::write(data_tag_hole_lo_ab, final_filename, true); + gmsh::view::write(data_tag_hole_tr_em, final_filename, true); + gmsh::view::write(data_tag_hole_tr_ab, final_filename, true); + gmsh::view::write(data_tag_elec_lo_em, final_filename, true); + gmsh::view::write(data_tag_elec_lo_ab, final_filename, true); + gmsh::view::write(data_tag_elec_tr_em, final_filename, true); + gmsh::view::write(data_tag_elec_tr_ab, final_filename, true); + } + + + // for (int index_band = 0; index_band < number_bands; ++index_band) { + // std::string band_name = "band_" + std::to_string(index_band); + // std::vector current_band_values(m_node_tags.size()); + // for (std::size_t index_node = 0; index_node < m_node_tags.size(); ++index_node) { + // current_band_values[index_node] = band_values[index_node * number_bands + index_band]; + // // std::cout << "band_values[" << index_node << "]: " << band_values[index_node * number_bands + index_band] << std::endl; + // } + // int data_tag = gmsh::view::add(band_name); + // if (m_node_tags.size() != current_band_values.size()) { + // std::cout << "number of nodes: " << m_node_tags.size() << std::endl; + // std::cout << "number of values: " << current_band_values.size() << std::endl; + // throw std::runtime_error("Number of nodes and number of values are not the same. Abort."); + // } + // gmsh::view::addHomogeneousModelData(data_tag, 0, model_file_name, "NodeData", m_node_tags, current_band_values); + // const int index_view = gmsh::view::getIndex(data_tag); + // std::string name_object_visibility = "View[" + std::to_string(index_view) + "].Visible"; + // gmsh::option::setNumber(name_object_visibility, 0); + // gmsh::view::write(data_tag, out_filename, true); + // } + + +} + + void ElectronPhonon::load_phonon_parameters(const std::string& filename) { YAML::Node config = YAML::LoadFile(filename); diff --git a/src/BZ_MESH/electron_phonon.hpp b/src/BZ_MESH/electron_phonon.hpp index 5011020..5d9f7ee 100644 --- a/src/BZ_MESH/electron_phonon.hpp +++ b/src/BZ_MESH/electron_phonon.hpp @@ -233,6 +233,7 @@ class ElectronPhonon : public BZ_States { void set_density(double rho) { m_rho = rho; } void compute_electron_phonon_rates_over_mesh(); + void add_electron_phonon_rates_to_mesh(const std::string& initial_filename, const std::string& final_filename); void export_rate_values(const std::string& filename) const; diff --git a/src/BZ_MESH/mesh_tetra.cpp b/src/BZ_MESH/mesh_tetra.cpp index 62f5986..2325869 100644 --- a/src/BZ_MESH/mesh_tetra.cpp +++ b/src/BZ_MESH/mesh_tetra.cpp @@ -16,8 +16,9 @@ #include #include +#include "Constants.hpp" #include "iso_triangle.hpp" -#include "Constants.hpp" + namespace bz_mesh { @@ -47,33 +48,55 @@ vector3 Tetra::compute_barycenter() const { 4.0; } - -std::vector Tetra::compute_gradient_scalar_field(std::array values){ - double v0 = values[0]; - double v1 = values[1]; - double v2 = values[2]; - double v3 = values[3]; - double x_0 = m_list_vertices[0]->get_position().x(); - double y_0 = m_list_vertices[0]->get_position().y(); - double z_0 = m_list_vertices[0]->get_position().z(); - double x_1 = m_list_vertices[1]->get_position().x(); - double y_1 = m_list_vertices[1]->get_position().y(); - double z_1 = m_list_vertices[1]->get_position().z(); - double x_2 = m_list_vertices[2]->get_position().x(); - double y_2 = m_list_vertices[2]->get_position().y(); - double z_2 = m_list_vertices[2]->get_position().z(); - double x_3 = m_list_vertices[3]->get_position().x(); - double y_3 = m_list_vertices[3]->get_position().y(); - double z_3 = m_list_vertices[3]->get_position().z(); - double det = 6.0 * m_signed_volume; - - - +vector3 Tetra::compute_gradient_scalar_field(std::array values) { + double value_0 = values[0]; + double value_1 = values[1]; + double value_2 = values[2]; + double value_3 = values[3]; + double x_0 = m_list_vertices[0]->get_position().x(); + double y_0 = m_list_vertices[0]->get_position().y(); + double z_0 = m_list_vertices[0]->get_position().z(); + double x_1 = m_list_vertices[1]->get_position().x(); + double y_1 = m_list_vertices[1]->get_position().y(); + double z_1 = m_list_vertices[1]->get_position().z(); + double x_2 = m_list_vertices[2]->get_position().x(); + double y_2 = m_list_vertices[2]->get_position().y(); + double z_2 = m_list_vertices[2]->get_position().z(); + double x_3 = m_list_vertices[3]->get_position().x(); + double y_3 = m_list_vertices[3]->get_position().y(); + double z_3 = m_list_vertices[3]->get_position().z(); + constexpr double one_half = 0.5; + + const double volABC = + one_half * (x_0 * y_1 * z_2 - x_0 * y_1 * z_3 - x_0 * y_2 * z_1 + x_0 * y_2 * z_3 + x_0 * y_3 * z_1 - x_0 * y_3 * z_2 - + x_1 * y_0 * z_2 + x_1 * y_0 * z_3 + x_1 * y_2 * z_0 - x_1 * y_2 * z_3 - x_1 * y_3 * z_0 + x_1 * y_3 * z_2 + + x_2 * y_0 * z_1 - x_2 * y_0 * z_3 - x_2 * y_1 * z_0 + x_2 * y_1 * z_3 + x_2 * y_3 * z_0 - x_2 * y_3 * z_1 - + x_3 * y_0 * z_1 + x_3 * y_0 * z_2 + x_3 * y_1 * z_0 - x_3 * y_1 * z_2 - x_3 * y_2 * z_0 + x_3 * y_2 * z_1); + const double grad_x = -one_half * + (value_0 * y_1 * z_2 - value_0 * y_1 * z_3 - value_0 * y_2 * z_1 + value_0 * y_2 * z_3 + value_0 * y_3 * z_1 - + value_0 * y_3 * z_2 - value_1 * y_0 * z_2 + value_1 * y_0 * z_3 + value_1 * y_2 * z_0 - value_1 * y_2 * z_3 - + value_1 * y_3 * z_0 + value_1 * y_3 * z_2 + value_2 * y_0 * z_1 - value_2 * y_0 * z_3 - value_2 * y_1 * z_0 + + value_2 * y_1 * z_3 + value_2 * y_3 * z_0 - value_2 * y_3 * z_1 - value_3 * y_0 * z_1 + value_3 * y_0 * z_2 + + value_3 * y_1 * z_0 - value_3 * y_1 * z_2 - value_3 * y_2 * z_0 + value_3 * y_2 * z_1) / + volABC; + const double grad_y = -one_half * + (-value_0 * x_1 * z_2 + value_0 * x_1 * z_3 + value_0 * x_2 * z_1 - value_0 * x_2 * z_3 - value_0 * x_3 * z_1 + + value_0 * x_3 * z_2 + value_1 * x_0 * z_2 - value_1 * x_0 * z_3 - value_1 * x_2 * z_0 + value_1 * x_2 * z_3 + + value_1 * x_3 * z_0 - value_1 * x_3 * z_2 - value_2 * x_0 * z_1 + value_2 * x_0 * z_3 + value_2 * x_1 * z_0 - + value_2 * x_1 * z_3 - value_2 * x_3 * z_0 + value_2 * x_3 * z_1 + value_3 * x_0 * z_1 - value_3 * x_0 * z_2 - + value_3 * x_1 * z_0 + value_3 * x_1 * z_2 + value_3 * x_2 * z_0 - value_3 * x_2 * z_1) / + volABC; + const double grad_z = -one_half * + (value_0 * x_1 * y_2 - value_0 * x_1 * y_3 - value_0 * x_2 * y_1 + value_0 * x_2 * y_3 + value_0 * x_3 * y_1 - + value_0 * x_3 * y_2 - value_1 * x_0 * y_2 + value_1 * x_0 * y_3 + value_1 * x_2 * y_0 - value_1 * x_2 * y_3 - + value_1 * x_3 * y_0 + value_1 * x_3 * y_2 + value_2 * x_0 * y_1 - value_2 * x_0 * y_3 - value_2 * x_1 * y_0 + + value_2 * x_1 * y_3 + value_2 * x_3 * y_0 - value_2 * x_3 * y_1 - value_3 * x_0 * y_1 + value_3 * x_0 * y_2 + + value_3 * x_1 * y_0 - value_3 * x_1 * y_2 - value_3 * x_2 * y_0 + value_3 * x_2 * y_1) / + volABC; + + return {grad_x, grad_y, grad_z}; } - - - 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(); @@ -85,8 +108,12 @@ void Tetra::compute_gradient_energy_at_bands() { const double eps_12 = (e_0 - energies_at_vertices[indices_sort[1]]); const double eps_13 = (e_0 - energies_at_vertices[indices_sort[2]]); const double eps_14 = (e_0 - energies_at_vertices[indices_sort[3]]); - const double gradient_energy = sqrt((eps_12 * eps_12 + eps_13 * eps_13 + eps_14 * eps_14)); - if (gradient_energy == 0) { + // const double gradient_energy = sqrt((eps_12 * eps_12 + eps_13 * eps_13 + eps_14 * eps_14)); + const vector3 gradient_energy = compute_gradient_scalar_field(energies_at_vertices); + double norm_grad_energy = gradient_energy.norm(); + // std::cout << "norm_grad_energy: " << norm_grad_energy << std::endl; + // std::cout << "norm_grad_energy quick: " << gradient_energy * (6 * m_signed_volume) << std::endl; + if (norm_grad_energy == 0) { std::cout << "Gradient energy is zero for tetra " << m_index << " at band " << band_index << std::endl; std::cout << "Energies: " << energies_at_vertices[0] << " " << energies_at_vertices[1] << " " << energies_at_vertices[2] << " " << energies_at_vertices[3] << std::endl; @@ -95,7 +122,7 @@ void Tetra::compute_gradient_energy_at_bands() { << m_list_vertices[2]->get_position() << "\n" << m_list_vertices[3]->get_position() << std::endl; } - m_gradient_energy_per_band.push_back(gradient_energy); + m_gradient_energy_per_band.push_back(norm_grad_energy); } } @@ -366,7 +393,7 @@ double Tetra::compute_tetra_dos_energy_band(double energy, std::size_t band_inde if (energy < m_min_energy_per_band[band_index] || energy > m_max_energy_per_band[band_index]) { return 0.0; } - const double renormalization = (2.0 / pow(M_2_PI, 3)) / EmpiricalPseudopotential::Constants::h_bar; + const double renormalization = 1.0; return renormalization * (1.0 / m_gradient_energy_per_band[band_index]) * this->compute_tetra_iso_surface_energy_band(energy, band_index); } diff --git a/src/BZ_MESH/mesh_tetra.hpp b/src/BZ_MESH/mesh_tetra.hpp index 32f7b83..647a8d9 100644 --- a/src/BZ_MESH/mesh_tetra.hpp +++ b/src/BZ_MESH/mesh_tetra.hpp @@ -102,11 +102,11 @@ class Tetra { std::array get_band_energies_at_vertices(std::size_t index_band) const; - double compute_signed_volume() const; - double get_signed_volume() const { return m_signed_volume; } - vector3 compute_edge(std::size_t index_vtx_1, std::size_t index_vtx_2) const; - void compute_gradient_energy_at_bands(); - std::vector Tetra::compute_gradient_scalar_field(std::array values); + double compute_signed_volume() const; + double get_signed_volume() const { return m_signed_volume; } + vector3 compute_edge(std::size_t index_vtx_1, std::size_t index_vtx_2) const; + void compute_gradient_energy_at_bands(); + vector3 compute_gradient_scalar_field(std::array values); vector3 compute_barycenter() const; bool is_location_inside(const vector3& location) const;