Skip to content

Commit

Permalink
up
Browse files Browse the repository at this point in the history
  • Loading branch information
@RemiHelleboid
RemiHelleboid committed Feb 16, 2024
1 parent 280c378 commit 0a9054c
Show file tree
Hide file tree
Showing 8 changed files with 338 additions and 37 deletions.
8 changes: 8 additions & 0 deletions apps/elph.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,14 @@ int main(int argc, char const *argv[])
ElectronPhonon.read_mesh_geometry_from_msh_file(mesh_band_input_file);
ElectronPhonon.read_mesh_bands_from_msh_file(mesh_band_input_file);

unsigned int nb_bands = ElectronPhonon.get_number_bands();
std::cout << "Number of bands: " << nb_bands << std::endl;
if (my_options.nrLevels > nb_bands) {
std::cout << "Number of bands requested is greater than the number of bands in the mesh file. Resetting to " << nb_bands << std::endl;
my_options.nrLevels = nb_bands;
}


ElectronPhonon.compute_electron_phonon_rates_over_mesh();

ElectronPhonon.export_rate_values("rates_all.csv");
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 @@ -258,7 +258,7 @@ int main(int argc, char *argv[]) {
std::cout << "Exporting DOS values to file: " << out_file_bands << std::endl;
export_multiple_vector_to_csv(out_file_bands + ".csv", list_header, dos_at_bands);

const std::string python_plot_dos = std::string(CMAKE_SOURCE_DIR) + "/python/plot_density_of_states.py";
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";
Expand Down
33 changes: 26 additions & 7 deletions python/plots/plot_phonon_rate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ def scatter_plot_rates(filename):
data = np.loadtxt(filename)
energy = data[:,0]
energy -= np.min(energy)
for i in range(1, data.shape[1]):
for i in range(2, data.shape[1]):
ax.scatter(energy, data[:,i], label=f"Mode {i}", s=1)
ax.set_xlabel("Energy (eV)")
ax.set_ylabel("Rate (s$^-1$)")
Expand All @@ -34,24 +34,43 @@ def scatter_plot_rates(filename):
ax.legend()

def plot_rates(filename):
fig, ax = plt.subplots()
data = np.loadtxt(filename)
energy = data[:,0]
energy -= np.min(energy)
dos = data[:,1]
ax.plot(energy, dos, label="DOS")
nb_modes = data.shape[1] - 2
fig, ax = plt.subplots(nb_modes, figsize=(6, 4*nb_modes))
for i in range(1, data.shape[1]-3):
ax[i-2].plot(energy, data[:,i], label=f"Mode {i}", color="black")
ax[i-2].set_xlabel("Energy (eV)")
ax[i-2].set_ylabel("Rate (s$^-1$)")
ax[i-2].set_title(f"Mode {i}")

# max_rate = np.max(data[:,2:])
# dos /= np.max(dos)
# ax.plot(energy, dos, label="DOS", color="black", linestyle="--")

# ax.legend()
fig.tight_layout()


def plo_dos(filename):
data = np.loadtxt(filename)
energy = data[:,0]
energy -= np.min(energy)
dos = data[:,1]
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(energy, dos, label="DOS", color="black", linestyle="--")
ax.set_xlabel("Energy (eV)")
ax.set_ylabel("DOS")
fig.tight_layout()
ax.legend()




if __name__ == "__main__":
parser = ArgumentParser()
parser.add_argument("-f", "--filename", type=str, required=True, help="Filename to plot")
args = parser.parse_args()
# scatter_plot_rates(args.filename)
scatter_plot_rates(args.filename)
plot_rates(args.filename)
plo_dos(args.filename)
plt.show()
110 changes: 110 additions & 0 deletions python/plots/test_gradient.ipynb
View file

Large diffs are not rendered by default.

169 changes: 149 additions & 20 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 "omp.h"

namespace bz_mesh {

double ElectronPhonon::bose_einstein_distribution(double energy, double temperature) {
Expand Down Expand Up @@ -90,12 +92,12 @@ RateValues ElectronPhonon::compute_electron_phonon_rate(int idx_n1, std::size_t
throw std::runtime_error("Energy phonon too high");
}
if (e_ph <= 0.0) {
std::cout << "Energy phonon: " << e_ph << std::endl << std::endl;
std::cout << "Q: " << q_ph << std::endl << std::endl;
std::cout << "Q: " << q_ph.Length() / m_material.get_fourier_factor() << std::endl << std::endl;
std::cout << "Energy (n1, k1): " << energy_n1_k1 << std::endl << std::endl;
// std::cout << "Energy phonon: " << e_ph << std::endl << std::endl;
// std::cout << "Q: " << q_ph << std::endl << std::endl;
// std::cout << "Q: " << q_ph.Length() / m_material.get_fourier_factor() << std::endl << std::endl;
// std::cout << "Energy (n1, k1): " << energy_n1_k1 << std::endl << std::endl;
bool is_in_bz = is_inside_mesh_geometry(q_ph / m_material.get_fourier_factor());
std::cout << "Is in BZ: " << is_in_bz << std::endl << std::endl;
// std::cout << "Is in BZ: " << is_in_bz << std::endl << std::endl;
// throw std::runtime_error("Energy phonon negative");
continue;
}
Expand Down Expand Up @@ -151,6 +153,106 @@ RateValues ElectronPhonon::compute_electron_phonon_rate(int idx_n1, std::size_t
return rates_k1_n1;
}

RateValues ElectronPhonon::compute_hole_phonon_rate(int idx_n1, std::size_t idx_k1) {
RateValues rates_k1_n1;
const std::vector<Tetra>& list_tetrahedra = m_list_tetrahedra;
auto indices_valence_bands = m_indices_valence_bands;
double energy_n1_k1 = m_list_vertices[idx_k1].get_energy_at_band(idx_n1);

for (auto&& idx_n2 : indices_valence_bands) {
for (auto&& tetra : list_tetrahedra) {
auto k_1 = m_list_vertices[idx_k1].get_position();
auto k_2 = tetra.compute_barycenter();
Vector3D<double> k1{k_1.x(), k_1.y(), k_1.z()};
Vector3D<double> k2{k_2.x(), k_2.y(), k_2.z()};

double overlap_integral = hole_overlap_integral(idx_n1, k1, idx_n2, k2);
auto q = k2 - k1;
auto initial_q = q;
Vector3D<double> q_ph{q.X, q.X, q.X};
bool is_in_bz = is_inside_mesh_geometry(q_ph / m_material.get_fourier_factor());

// No Umklapp process for now.
if (!is_in_bz) {
auto new_q = retrieve_k_inside_mesh_geometry(q_ph / m_material.get_fourier_factor()) * m_material.get_fourier_factor();
q_ph = Vector3D<double>{new_q.x(), new_q.y(), new_q.z()};
// std::cout << "q: " << q_ph / m_material.get_fourier_factor() << " new_q: " << new_q / m_material.get_fourier_factor() <<
// std::endl;
}
is_in_bz = is_inside_mesh_geometry(q_ph / m_material.get_fourier_factor());
if (!is_in_bz) {
throw std::runtime_error("Q is not inside the BZ");
}

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()) * EmpiricalPseudopotential::Constants::h_bar_eV;
if (e_ph > 100e-3) {
std::cout << "Energy phonon: " << e_ph << std::endl << std::endl;
throw std::runtime_error("Energy phonon too high");
}
if (e_ph <= 0.0) {
// std::cout << "Energy phonon: " << e_ph << std::endl << std::endl;
// std::cout << "Q: " << q_ph << std::endl << std::endl;
// std::cout << "Q: " << q_ph.Length() / m_material.get_fourier_factor() << std::endl << std::endl;
// std::cout << "Energy (n1, k1): " << energy_n1_k1 << std::endl << std::endl;
// bool is_in_bz = is_inside_mesh_geometry(q_ph / m_material.get_fourier_factor());
// std::cout << "Is in BZ: " << is_in_bz << std::endl << std::endl;
// throw std::runtime_error("Energy phonon negative");
continue;
}
// 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 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;
if (!tetra.is_energy_inside_band(energy_final, idx_n2)) {
continue;
}
// 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;
DeformationPotential deformation_potential = (mode_direction.first == PhononMode::acoustic)
? m_acoustic_deformation_potential_h
: m_optical_deformation_potential_h;
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;
// std::cout << this->get_volume() << std::endl;
rate_value /= this->get_volume();
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;
std::cout << "Overlap integral: " << overlap_integral << std::endl;
std::cout << "Deformation potential: " << deformation_potential_value << std::endl;
std::cout << "Bose part: " << bose_part << std::endl;
std::cout << "DOS tetra: " << dos_tetra << std::endl;
std::cout << "Energy final: " << energy_final << std::endl;
std::cout << "q: " << q_ph << std::endl;
std::cout << "initial q: " << initial_q / m_material.get_fourier_factor() << std::endl;
std::cout << "Energy phonon: " << e_ph << std::endl;
std::cout << "Energy (n1, k1): " << energy_n1_k1 << std::endl;
std::cout << "--------------------------------------------------------------------------------" << std::endl;
}
// std::cout << "Rate value: " << (deformation_potential_value) << std::endl;
RateValue rate(phonon_mode, phonon_direction, phonon_event, rate_value);
rates_k1_n1.add_rate(rate);
}
}
}
}
// rates_k1_n1.print_rates();
// std::cout << "*************************************************************************************************************" <<
// std::endl;
return rates_k1_n1;
}

void ElectronPhonon::compute_electron_phonon_rates_over_mesh() {
auto indices_conduction_bands = m_indices_conduction_bands;
auto min_idx_conduction_band = *std::min_element(indices_conduction_bands.begin(), indices_conduction_bands.end());
Expand All @@ -161,13 +263,22 @@ 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 = 0.01;
double p_compute_rate = 1.0;

#pragma omp parallel for schedule(dynamic) num_threads(8)
#pragma omp parallel for schedule(dynamic)
for (std::size_t idx_k1 = 0; idx_k1 < m_list_vertices.size(); ++idx_k1) {
double r = dis(gen);
for (std::size_t idx_n1 = 0; idx_n1 < min_idx_conduction_band; ++idx_n1) {
if (r > p_compute_rate) {
std::array<double, 8> array = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
m_list_vertices[idx_k1].add_electron_phonon_rates(array);
continue;
}
auto hole_rate = compute_hole_phonon_rate(idx_n1, idx_k1);
auto array = hole_rate.to_array();
m_list_vertices[idx_k1].add_electron_phonon_rates(array);
}

std::cout << "Computing rates for k-point " << idx_k1 << std::endl;
for (std::size_t idx_n1 = min_idx_conduction_band; idx_n1 < max_idx_conduction_band; ++idx_n1) {
if (r > p_compute_rate) {
std::array<double, 8> array = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
Expand All @@ -179,7 +290,11 @@ void ElectronPhonon::compute_electron_phonon_rates_over_mesh() {
m_list_vertices[idx_k1].add_electron_phonon_rates(array);
}
}
if (omp_get_thread_num() == 0){
std::cout << "\rComputing rates for k-point " << idx_k1 << std::flush;}

}
std::cout << std::endl;
file.close();
// Add rate to the mesh
}
Expand All @@ -188,23 +303,28 @@ void ElectronPhonon::compute_plot_electron_phonon_rates_vs_energy_over_mesh(int
double max_energy,
double energy_step,
const std::string& filename) {
std::ofstream file(filename);
std::vector<double> energies;
std::vector<double> list_dos;
std::ofstream file(filename);
std::vector<double> energies;
std::vector<double> list_dos;
std::vector<std::array<double, 8>> list_mean_rates;
double total_dos = 0.0;
for (double energy = 0.0; energy < max_energy; energy += energy_step) {
double total_dos = 0.0;
for (double energy = 1.0; energy < max_energy; energy += energy_step) {
std::cout << "\rEnergy: " << energy << " Max energy: " << max_energy << std::flush;
double dos = 0.0;
double dos = 0.0;
std::array<double, 8> mean_rates;
std::fill(mean_rates.begin(), mean_rates.end(), 0.0);
for (int idx_band = 4; idx_band < nb_bands; ++idx_band) {
for (auto&& tetra : m_list_tetrahedra) {
for (auto&& tetra : m_list_tetrahedra) {
for (int idx_band = 4; idx_band < nb_bands; ++idx_band) {
double dos_tetra = tetra.compute_tetra_dos_energy_band(energy, idx_band);

if (std::isnan(dos_tetra) || std::isinf(dos_tetra)) {
std::cout << "Energy: " << energy << " Band: " << idx_band << " DOS tetra: " << dos_tetra << std::endl;
throw std::runtime_error("DOS tetra is NaN or Inf");
}
dos += dos_tetra;
std::array<double, 8> mean_rates = tetra.get_mean_electron_phonon_rates(idx_band);
for (std::size_t idx_rate = 0; idx_rate < mean_rates.size(); ++idx_rate) {
mean_rates[idx_rate] += mean_rates[idx_rate] * dos_tetra;
std::array<double, 8> tetra_mean_rates = tetra.get_mean_electron_phonon_rates(idx_band);
for (std::size_t idx_rate = 0; idx_rate < tetra_mean_rates.size(); ++idx_rate) {
mean_rates[idx_rate] += tetra_mean_rates[idx_rate] * dos_tetra;
}
}
}
Expand All @@ -216,6 +336,7 @@ void ElectronPhonon::compute_plot_electron_phonon_rates_vs_energy_over_mesh(int
for (std::size_t idx_energy = 0; idx_energy < energies.size(); ++idx_energy) {
file << energies[idx_energy] << " " << list_dos[idx_energy] << " ";
for (auto&& rate : list_mean_rates[idx_energy]) {
rate /= total_dos;
file << rate << " ";
}
file << std::endl;
Expand Down Expand Up @@ -300,6 +421,14 @@ void ElectronPhonon::load_phonon_parameters(const std::string& filename) {
} else {
m_optical_deformation_potential_e = deformationPotential;
}
} else if (carrierType == "hole") {
if (mode == PhononMode::acoustic) {
m_acoustic_deformation_potential_h = deformationPotential;
} else {
m_optical_deformation_potential_h = deformationPotential;
}
} else {
throw std::runtime_error("Unknown carrier type.");
}
}
}
Expand All @@ -310,7 +439,7 @@ void ElectronPhonon::export_rate_values(const std::string& filename) const {
for (auto&& vertex : m_list_vertices) {
std::vector<std::array<double, 8>> all_rates = vertex.get_electron_phonon_rates_all_bands();
for (std::size_t idx_band = 0; idx_band < all_rates.size(); ++idx_band) {
double energy = vertex.get_energy_at_band(idx_band + 4);
double energy = vertex.get_energy_at_band(idx_band);
file << energy << " ";
for (auto&& rate : all_rates[idx_band]) {
file << rate << " ";
Expand Down
3 changes: 3 additions & 0 deletions src/BZ_MESH/electron_phonon.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,6 +211,8 @@ class ElectronPhonon : public BZ_States {
HoleOverlapIntParams m_hole_overlap_int_params;
DeformationPotential m_acoustic_deformation_potential_e;
DeformationPotential m_optical_deformation_potential_e;
DeformationPotential m_acoustic_deformation_potential_h;
DeformationPotential m_optical_deformation_potential_h;

std::map<PhononModeDirection, PhononDispersion> m_phonon_dispersion;

Expand All @@ -225,6 +227,7 @@ class ElectronPhonon : public BZ_States {
double hole_overlap_integral(int n1, const Vector3D<double>& k1, int n2, const Vector3D<double>& k2);

RateValues compute_electron_phonon_rate(int idx_n1, std::size_t idx_k1);
RateValues compute_hole_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; }
Expand Down
Loading

0 comments on commit 0a9054c

Please sign in to comment.