From 88efb759c56dbc3f70a5a2383244dad6f7771c59 Mon Sep 17 00:00:00 2001 From: RemiHelleboid Date: Wed, 22 May 2024 10:59:46 +0200 Subject: [PATCH] Up --- apps/impact_io.cpp | 9 +- src/BZ_MESH/bz_states.cpp | 15 ++-- src/BZ_MESH/dielectric_mesh.cpp | 43 ++++++--- src/BZ_MESH/impact_ionization.cpp | 141 ++++++++++++------------------ src/BZ_MESH/impact_ionization.hpp | 3 +- src/BZ_MESH/mesh_tetra.hpp | 16 +++- 6 files changed, 112 insertions(+), 115 deletions(-) diff --git a/apps/impact_io.cpp b/apps/impact_io.cpp index 01d098e..ef541cc 100644 --- a/apps/impact_io.cpp +++ b/apps/impact_io.cpp @@ -31,7 +31,7 @@ int main(int argc, char *argv[]) { TCLAP::ValueArg arg_dielectric_file("y", "dielectric", "File (.msh) with the dielectric function.", true, "", "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_bands("b", "nbands", "Number of bands to consider", false, 16, "int"); TCLAP::ValueArg arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int"); TCLAP::ValueArg arg_radius_BZ("R", "radiusBZ", @@ -67,16 +67,17 @@ int main(int argc, char *argv[]) { my_options.nrThreads = arg_nb_threads.getValue(); my_options.print_options(); int nb_bands_to_use = arg_nb_bands.getValue(); + int nb_threads = arg_nb_threads.getValue(); auto start = std::chrono::high_resolution_clock::now(); EmpiricalPseudopotential::Material current_material = materials.materials.at(arg_material.getValue()); bz_mesh::ImpactIonization my_impact_ionization(current_material, arg_mesh_file.getValue()); - my_impact_ionization.read_dielectric_file(arg_dielectric_file.getValue()); - + my_impact_ionization.read_dielectric_file(arg_dielectric_file.getValue()); + my_impact_ionization.interp_test_dielectric_function("test_dielectric_function.csv"); my_impact_ionization.set_max_radius_G0_BZ(arg_radius_BZ.getValue()); - my_impact_ionization.compute_eigenstates(); + my_impact_ionization.compute_eigenstates(nb_threads); auto end = std::chrono::high_resolution_clock::now(); diff --git a/src/BZ_MESH/bz_states.cpp b/src/BZ_MESH/bz_states.cpp index 6d41427..427c9ef 100644 --- a/src/BZ_MESH/bz_states.cpp +++ b/src/BZ_MESH/bz_states.cpp @@ -31,8 +31,11 @@ void BZ_States::compute_eigenstates(int nb_threads) { for (int i = 0; i < nb_threads; i++) { hamiltonian_per_thread.push_back(EmpiricalPseudopotential::Hamiltonian(m_material, m_basisVectors)); } +#pragma omp parallel for schedule(dynamic) num_threads(nb_threads) for (std::size_t idx_k = 0; idx_k < m_list_vertices.size(); ++idx_k) { - std::cout << "\rComputing eigenstates for k = " << idx_k << "/" << m_list_vertices.size() << std::flush; + if (omp_get_thread_num() == 0) { + std::cout << "\rComputing eigenstates for k = " << idx_k << "/" << m_list_vertices.size() << std::flush; + } auto k_point = Vector3D(m_list_vertices[idx_k].get_position().x(), m_list_vertices[idx_k].get_position().y(), m_list_vertices[idx_k].get_position().z()); @@ -93,7 +96,7 @@ void BZ_States::compute_dielectric_function(const std::vector& list_ener constexpr double one_fourth = 1.0 / 4.0; std::vector dielectric_function_real_at_energies(list_energies.size(), 0.0); - double total_volume = 0.0; + double total_volume = 0.0; for (std::size_t idx_tetra = 0; idx_tetra < nb_tetra; ++idx_tetra) { std::cout << "\rComputing dielectric function for tetrahedron " << idx_tetra << "/" << nb_tetra << std::flush; std::array list_idx_vertices = m_list_tetrahedra[idx_tetra].get_list_indices_vertices(); @@ -133,14 +136,12 @@ void BZ_States::compute_dielectric_function(const std::vector& list_ener std::cout << "\n"; std::cout << "Total volume: " << total_volume << std::endl; - double q_squared = m_q_shift.Length() * m_q_shift.Length(); + double q_squared = m_q_shift.Length() * m_q_shift.Length(); double pre_factor = 2.0 * M_PI / q_squared; for (std::size_t index_energy = 0; index_energy < list_energies.size(); ++index_energy) { m_dielectric_function_real[index_energy] = 1.0 + pre_factor * dielectric_function_real_at_energies[index_energy] / total_volume; } std::cout << "EPS[0] = " << m_dielectric_function_real[0] << std::endl; - - } // Export the dielectric function to a file in the format (energy, dielectric function) (csv format). @@ -175,10 +176,6 @@ void BZ_States::export_full_eigenstates() const { } } - // void BZ_States::populate_vtx_dielectric_function(const std::vector& energies, double eta_smearing); - - - } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/dielectric_mesh.cpp b/src/BZ_MESH/dielectric_mesh.cpp index 8a7f04e..a33f00a 100644 --- a/src/BZ_MESH/dielectric_mesh.cpp +++ b/src/BZ_MESH/dielectric_mesh.cpp @@ -82,14 +82,17 @@ void DielectricMesh::read_dielectric_file(const std::string& filename) { } gmsh::finalize(); m_energies = energies; - m_dielectric_function.resize(real_dielectric.size()); - for (std::size_t idx_node = 0; idx_node < real_dielectric.size(); ++idx_node) { + m_dielectric_function.resize(m_list_vertices.size()); + for (std::size_t idx_node = 0; idx_node < m_list_vertices.size(); ++idx_node) { m_dielectric_function[idx_node].resize(energies.size()); for (std::size_t idx_energy = 0; idx_energy < energies.size(); ++idx_energy) { m_dielectric_function[idx_node][idx_energy] = - std::complex(real_dielectric[idx_node][idx_energy], imag_dielectric[idx_node][idx_energy]); + std::complex(real_dielectric[idx_energy][idx_node], imag_dielectric[idx_energy][idx_node]); } } + std::cout << "Size of the dielectric function: " << m_dielectric_function.size() << " x " << m_dielectric_function[0].size() + << std::endl; + std::cout << "Dielectric function read." << std::endl; } std::pair DielectricMesh::find_closest_energy(double energy) const { @@ -115,21 +118,35 @@ std::pair DielectricMesh::find_closest_energy(double energy } complex_d DielectricMesh::interpolate_dielectric_function(const vector3& k, double energy) const { - auto k_positive{vector3{std::abs(k.x()), std::abs(k.y()), std::abs(k.z())}}; + auto k_positive{vector3{std::abs(k.x()), std::abs(k.y()), std::abs(k.z())}}; Tetra* p_tetra = find_tetra_at_location(k_positive); if (p_tetra == nullptr) { + std::cerr << "Tetra not found at location " << k_positive << std::endl; return 0.0; } - const std::array barycentric_coordinates = p_tetra->compute_barycentric_coordinates(k); + std::array barycentric_coordinates = p_tetra->compute_barycentric_coordinates(k_positive); - const std::size_t idx_node = p_tetra->get_index(); - const std::pair closest_energy = find_closest_energy(energy); - const std::size_t idx_energy = closest_energy.first; - const double t = closest_energy.second; - const complex_d dielectric_function_0 = m_dielectric_function[idx_node][idx_energy]; - const complex_d dielectric_function_1 = m_dielectric_function[idx_node][idx_energy + 1]; - complex_d dielectric_function = dielectric_function_0 * (1.0 - t) + dielectric_function_1 * t; - return dielectric_function; + std::array list_indices_vertices = p_tetra->get_list_indices_vertices(); + std::pair closest_energy = find_closest_energy(energy); + std::size_t idx_energy = closest_energy.first; + double t = closest_energy.second; + std::vector> dielectric_function_low(4); + std::vector> dielectric_function_high(4); + for (std::size_t idx_vertex = 0; idx_vertex < 4; ++idx_vertex) { + std::cout << "Vertex: " << list_indices_vertices[idx_vertex] << std::endl; + dielectric_function_low[idx_vertex] = m_dielectric_function[idx_energy][list_indices_vertices[idx_vertex]]; + dielectric_function_high[idx_vertex] = m_dielectric_function[idx_energy + 1][list_indices_vertices[idx_vertex]]; + } + std::complex dielectric_function_interpolated_low = + p_tetra->interpolate_at_position(barycentric_coordinates, dielectric_function_low); + std::complex dielectric_function_interpolated_high = + p_tetra->interpolate_at_position(barycentric_coordinates, dielectric_function_high); + + std::cout << "Idx energy: " << idx_energy << std::endl; + std::cout << "Energy: " << energy << std::endl; + std::cout << "Closest energy: " << m_energies[idx_energy] << std::endl; + std::cout << "Interpolated energy: " << (1 - t) * m_energies[idx_energy] + t * m_energies[idx_energy + 1] << std::endl; + return (1 - t) * dielectric_function_interpolated_low + t * dielectric_function_interpolated_high; } } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/impact_ionization.cpp b/src/BZ_MESH/impact_ionization.cpp index 95c6f12..8e3f706 100644 --- a/src/BZ_MESH/impact_ionization.cpp +++ b/src/BZ_MESH/impact_ionization.cpp @@ -48,12 +48,42 @@ void ImpactIonization::read_dielectric_file(const std::string& filename) { m_dielectric_mesh.read_dielectric_file(filename); } -void ImpactIonization::compute_eigenstates() { - int nb_bands_to_use = 4; +void ImpactIonization::interp_test_dielectric_function(std::string filename) { + double eps = 1e-6; + double x0 = eps; + double y0 = eps; + double z0 = eps; + double x1 = 3.0; + double y1 = 3.0; + double z1 = 3.0; + int nx = 20; + int ny = 20; + int nz = 20; + double dx = (x1 - x0) / nx; + double dy = (y1 - y0) / ny; + double dz = (z1 - z0) / nz; + std::ofstream file(filename); + for (int i = 0; i < nx; ++i) { + for (int j = 0; j < ny; ++j) { + for (int k = 0; k < nz; ++k) { + double x = x0 + i * dx; + double y = y0 + j * dy; + double z = z0 + k * dz; + vector3 position(x, y, z); + complex_d epsilon = m_dielectric_mesh.interpolate_dielectric_function(position, 0.0102); + file << x << ", " << y << ", " << z << ", " << epsilon.real() << ", " << epsilon.imag() << std::endl; + std::cout << "Position: " << position << " epsilon: " << epsilon << std::endl; + } + } + } + file.close(); +} + +void ImpactIonization::compute_eigenstates(int nb_threads) { + int nb_bands_to_use = 16; bz_mesh::BZ_States my_bz_mesh(m_material); my_bz_mesh.set_nb_bands(nb_bands_to_use); EmpiricalPseudopotential::BandStructure band_structure{}; - int nb_threads = 1; int nb_nearest_neighbors = 10; bool nonlocal_epm = false; bool enable_soc = false; @@ -68,7 +98,7 @@ void ImpactIonization::compute_eigenstates() { const vector3 b3 = {1.0, 1.0, -1.0}; // std::cout << "k: " << k << std::endl; - + double factor = 2.0 * M_PI / m_material.get_lattice_constant_meter(); // test k + G std::vector list_n_k = {0, 1, -1, 2, -2, 3, -3, 4, -4}; // std::vector list_n_k = {0, 1, -1, 2, -2}; @@ -81,7 +111,9 @@ void ImpactIonization::compute_eigenstates() { if (G_BZ.norm() > m_max_radius_G0_BZ) { continue; } - std::cout << "G_BZ: " << G_BZ << std::endl; + std::cout << "G_BZ: " << G_BZ << " --> " << G_BZ.norm() << std::endl; + G_BZ = G_BZ * factor; + auto ptr_BZ_states = std::make_unique(m_material); ptr_BZ_states->set_nb_bands(nb_bands_to_use); ptr_BZ_states->set_basis_vectors(basis); @@ -96,88 +128,25 @@ void ImpactIonization::compute_eigenstates() { } double ImpactIonization::compute_impact_ionization_rate(int idx_n1, const Vector3D& k1) { - double rate = 0.0; - int nb_conduction_bands = 4; - int nb_valence_bands = 4; - // Pseudo code + constexpr int nb_valence_bands = 3; + constexpr int nb_conduction_bands = 4; - return rate; -} + std::size_t nb_vtx = m_list_BZ_states[0]->get_list_vertices().size(); -// /** -// * @brief Compute the direct impact ionization matrix element (Ma) for a given set of indices. -// * k2 is not an input parameter, it is computed from k1, k1_prime and k2_prime, using the conservation of momentum. -// * k2 = k1_prime + k2_prime - k1 -// * -// * -// * @param idx_n1 -// * @param idx_n1_prime -// * @param idx_n2 -// * @param idx_n2_prime -// * @param idx_k1 -// * @param idx_k1_prime -// * @param idx_k2_prime -// * @return double -// */ -// std::array ImpactIonization::compute_direct_indirect_impact_ionization_matrix_element(int idx_n1, -// int idx_n1_prime, -// int idx_n2, -// int idx_n2_prime, -// int idx_k1, -// int idx_k1_prime, -// int idx_k2, -// int idx_k2_prime) const { -// vector3 k1 = m_list_vertices[idx_k1].get_position(); -// vector3 k2 = m_list_vertices[idx_k2].get_position(); -// vector3 k1_prime = m_list_vertices[idx_k1_prime].get_position(); -// vector3 k2_prime = m_list_vertices[idx_k2_prime].get_position(); -// complex_d Ma = 0.0; -// complex_d Mb = 0.0; - -// double e_charge = EmpiricalPseudopotential::Constants::q; -// double eps_zero = EmpiricalPseudopotential::Constants::eps_zero; -// double volume = compute_mesh_volume(); -// double fourier_factor = 2.0 * M_PI / volume; - -// auto basis_vectors = get_basis_vectors(); -// std::size_t size_basis_vectors = basis_vectors.size(); -// for (std::size_t idx_g1 = 0; idx_g1 < size_basis_vectors; ++idx_g1) { -// for (std::size_t idx_g2 = 0; idx_g2 < size_basis_vectors; ++idx_g2) { -// for (std::size_t idx_g1_prime = 0; idx_g1_prime < size_basis_vectors; ++idx_g1_prime) { -// for (std::size_t idx_g2_prime = 0; idx_g2_prime < size_basis_vectors; ++idx_g2_prime) { -// vector3 g1 = {basis_vectors[idx_g1].X * fourier_factor, -// basis_vectors[idx_g1].Y * fourier_factor, -// basis_vectors[idx_g1].Z * fourier_factor}; -// vector3 g1_prime = {basis_vectors[idx_g1_prime].X * fourier_factor, -// basis_vectors[idx_g1_prime].Y * fourier_factor, -// basis_vectors[idx_g1_prime].Z * fourier_factor}; -// vector3 g2_prime = {basis_vectors[idx_g2_prime].X * fourier_factor, -// basis_vectors[idx_g2_prime].Y * fourier_factor, -// basis_vectors[idx_g2_prime].Z * fourier_factor}; -// vector3 q_a = k1_prime + g1_prime - k1 - g1; -// vector3 q_b = k2_prime + g2_prime - k1 - g1; -// auto state_k1 = m_eigenvectors_k[idx_k1](idx_n1, idx_g1); -// auto state_k2 = m_eigenvectors_k[idx_k2](idx_n2, idx_g2); -// auto state_k1_prime = m_eigenvectors_k[idx_k1_prime](idx_n1_prime, idx_g1_prime); -// auto state_k2_prime = m_eigenvectors_k[idx_k2_prime](idx_n2_prime, idx_g2_prime); - -// double omega_a = m_eigenvalues_k[idx_k1](idx_n1) - m_eigenvalues_k[idx_k1_prime](idx_n1_prime); -// double omega_b = m_eigenvalues_k[idx_k2_prime](idx_n2_prime) - m_eigenvalues_k[idx_k1](idx_n1); -// complex_d dielectric_func_a = get_dielectric_function(q_a, omega_a); -// complex_d pre_factor_a = e_charge * e_charge / (eps_zero * dielectric_func_a * q_a.norm() * q_a.norm() * -// volume); complex_d dielectric_func_b = get_dielectric_function(q_b, omega_b); complex_d core_overlap = -// std::conj(state_k1_prime) * std::conj(state_k2_prime) * state_k1 * state_k2; complex_d micro_matrix_element_a = -// core_overlap * pre_factor_a; complex_d pre_factor_b = e_charge * e_charge / (eps_zero * dielectric_func_b * -// q_b.norm() * q_b.norm() * volume); complex_d micro_matrix_element_b = core_overlap * pre_factor_b; - -// Ma += micro_matrix_element_a; -// Mb += micro_matrix_element_b; -// } -// } -// } -// } - -// return {Ma, Mb}; -// } + for (int n1_prime = 0; n1_prime < nb_valence_bands; ++n1_prime) { + for (int n2 = 0; n2 < nb_conduction_bands; ++n2) { + for (int n2_prime = 0; n2_prime < nb_conduction_bands; ++n2_prime) { + for (int k1_prime = 0; k1_prime < nb_vtx; ++k1_prime) { + + } + } + } + } + + + + + +} } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/impact_ionization.hpp b/src/BZ_MESH/impact_ionization.hpp index 2cc9aa8..85de6cc 100644 --- a/src/BZ_MESH/impact_ionization.hpp +++ b/src/BZ_MESH/impact_ionization.hpp @@ -67,10 +67,11 @@ class ImpactIonization { public: ImpactIonization(const EmpiricalPseudopotential::Material& material, const std::string& initial_mesh_path); void read_dielectric_file(const std::string& filename); + void interp_test_dielectric_function(std::string filename); double get_max_radius_G0_BZ() const { return m_max_radius_G0_BZ; } void set_max_radius_G0_BZ(double max_radius_G0_BZ) { m_max_radius_G0_BZ = max_radius_G0_BZ; } - void compute_eigenstates(); + void compute_eigenstates(int nb_threads = 1); std::array compute_direct_indirect_impact_ionization_matrix_element(int idx_n1, int idx_n1_prime, diff --git a/src/BZ_MESH/mesh_tetra.hpp b/src/BZ_MESH/mesh_tetra.hpp index 928a5e1..bbf62c1 100644 --- a/src/BZ_MESH/mesh_tetra.hpp +++ b/src/BZ_MESH/mesh_tetra.hpp @@ -11,6 +11,7 @@ #pragma once #include +#include #include #include #include @@ -20,6 +21,7 @@ #include "bbox_mesh.hpp" #include "mesh_vertex.hpp" + namespace bz_mesh { class Tetra { @@ -133,8 +135,9 @@ class Tetra { std::array get_mean_electron_phonon_rates(int band_index) const; - double interpolate_scalar_at_position(const std::array& barycentric_coordinates, - const std::vector& scalar_field) 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 @@ -145,6 +148,15 @@ class Tetra { } return interpolated_value; } + std::complex interpolate_at_position(const std::array& barycentric_coordinates, + const std::vector>& field) const { + std::complex interpolated_value = 0.0; + 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); }