Skip to content

Commit

Permalink
up
Browse files Browse the repository at this point in the history
  • Loading branch information
@RemiHelleboid
RemiHelleboid committed Feb 19, 2024
1 parent 0a9054c commit ec53eea
Show file tree
Hide file tree
Showing 7 changed files with 184 additions and 44 deletions.
2 changes: 2 additions & 0 deletions apps/elph.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::chrono::seconds>(stop - start).count() << " seconds" << std::endl;
Expand Down
2 changes: 1 addition & 1 deletion apps/mpi_DOS_MeshBZ.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
1 change: 0 additions & 1 deletion src/BZ_MESH/bz_mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

Expand Down
125 changes: 118 additions & 7 deletions src/BZ_MESH/electron_phonon.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@
#include "Vector3D.h"
#include "bz_states.hpp"

#include "gmsh.h"

#include "omp.h"

namespace bz_mesh {
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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);
}
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -263,7 +266,7 @@ void ElectronPhonon::compute_electron_phonon_rates_over_mesh() {
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<double> 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) {
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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<std::size_t> node_tags;
std::vector<double> nodeCoords;
std::vector<double> 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<max_idx_conduction_band; ++idx_val_band) {
std::vector<double> rates_hole_lo_em(m_list_vertices.size());
std::vector<double> rates_hole_lo_ab(m_list_vertices.size());
std::vector<double> rates_hole_tr_em(m_list_vertices.size());
std::vector<double> rates_hole_tr_ab(m_list_vertices.size());
std::vector<double> rates_elec_lo_em(m_list_vertices.size());
std::vector<double> rates_elec_lo_ab(m_list_vertices.size());
std::vector<double> rates_elec_tr_em(m_list_vertices.size());
std::vector<double> 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<double> 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);

Expand Down
1 change: 1 addition & 0 deletions src/BZ_MESH/electron_phonon.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
Loading

0 comments on commit ec53eea

Please sign in to comment.