Continue implementation of ii

RemiHelleboid · May 25, 2024 · 0de85b3 · 0de85b3
1 parent 88efb75
commit 0de85b3
Show file tree

Hide file tree

Showing 12 changed files with 139 additions and 47 deletions.
diff --git a/apps/impact_io.cpp b/apps/impact_io.cpp
@@ -78,6 +78,9 @@ int main(int argc, char *argv[]) {
 
     my_impact_ionization.set_max_radius_G0_BZ(arg_radius_BZ.getValue());
     my_impact_ionization.compute_eigenstates(nb_threads);
+    int idxn1 = 0;
+    int idxk1 = 0;
+    my_impact_ionization.compute_impact_ionization_rate(idxn1, idxk1);
 
     auto end = std::chrono::high_resolution_clock::now();
 

diff --git a/apps/mpi_DOS_MeshBZ.cpp b/apps/mpi_DOS_MeshBZ.cpp
@@ -239,7 +239,6 @@ int main(int argc, char *argv[]) {
         std::cout << "------------------------------------------------------------------\n";
 
         std::cout << "Exporting DOS values to file." << std::endl;
-        // dos_at_bands[band_index][energy_index] = dos_value;
         std::vector<std::vector<double>> dos_at_bands(number_bands * 2);
         for (std::size_t index_value_dos = 0; index_value_dos < total_number_dos; index_value_dos++) {
             int    band_index = list_bands[index_value_dos];

diff --git a/src/BZ_MESH/bz_mesh.hpp b/src/BZ_MESH/bz_mesh.hpp
@@ -89,6 +89,9 @@ class MeshBZ {
     MeshBZ(const EmpiricalPseudopotential::Material& material) : m_material(material){};
     MeshBZ(const MeshBZ&) = default;
 
+
+    vector3 get_vertex_position(std::size_t idx_vtx) const { return m_list_vertices[idx_vtx].get_position(); }
+
     vector3             get_center() const { return m_center; }
     void                shift_bz_center(const vector3& shift);
 
@@ -132,7 +135,7 @@ class MeshBZ {
     double compute_mesh_volume() const;
     double compute_iso_surface(double iso_energy, int band_index) const;
     double compute_dos_at_energy_and_band(double iso_energy, int band_index) const;
-    double compute_overlap_integral_impact_ionization_electrons(double energy);
+    double compute_dos_like_integral(double iso_energy, const std::vector<double>& f_values) const;
 
     std::vector<std::vector<double>> compute_dos_band_at_band(int         band_index,
                                                               double      min_energy,

diff --git a/src/BZ_MESH/bz_states.cpp b/src/BZ_MESH/bz_states.cpp
@@ -47,6 +47,7 @@ void BZ_States::compute_eigenstates(int nb_threads) {
         m_eigenvectors_k[idx_k] = hamiltonian_per_thread[idx_thread].get_eigenvectors();
         auto nb_rows            = m_eigenvectors_k[idx_k].rows();
         m_eigenvectors_k[idx_k].conservativeResize(nb_rows, m_nb_bands);
+        // std::cout << "\r Nb rows = " << m_eigenvectors_k[idx_k].rows() << " Nb cols = " << m_eigenvectors_k[idx_k].cols() << std::endl;
     }
     std::cout << std::endl;
 }

diff --git a/src/BZ_MESH/bz_states.hpp b/src/BZ_MESH/bz_states.hpp
@@ -22,11 +22,10 @@ class BZ_States : public MeshBZ {
  protected:
     int m_nb_bands = 0;
 
-
     std::vector<Vector3D<int>> m_basisVectors;
 
-    std::vector<Eigen::VectorXd>  m_eigenvalues_k;
-    std::vector<Eigen::VectorXd>  m_eigenvalues_k_plus_q;
+    std::vector<Eigen::VectorXd> m_eigenvalues_k;
+    std::vector<Eigen::VectorXd> m_eigenvalues_k_plus_q;
 
     std::vector<Eigen::MatrixXcd> m_eigenvectors_k;
     std::vector<Eigen::MatrixXcd> m_eigenvectors_k_plus_q;
@@ -78,14 +77,16 @@ class BZ_States : public MeshBZ {
     const std::vector<double>& get_energies() const { return m_list_energies; }
     void                       set_energies(const std::vector<double>& energies) { m_list_energies = energies; }
 
+    const std::vector<Eigen::MatrixXcd>& get_eigen_states() const { return m_eigenvectors_k; }
+
     void compute_dielectric_function(const std::vector<double>& energies, double eta_smearing, int nb_threads = 1);
     void export_dielectric_function(const std::string& prefix) const;
 
     void populate_vtx_dielectric_function(const std::vector<double>& energies, double eta_smearing);
 
     std::complex<double> get_dielectric_function(const vector3& q, double energy) const {
-      // TODO: implement this function
-      return 1.0;
+        // TODO: implement this function
+        return 1.0;
     }
 
     void export_full_eigenstates() const;

diff --git a/src/BZ_MESH/electron_phonon.cpp b/src/BZ_MESH/electron_phonon.cpp
@@ -128,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_e;
 
                     if (rate_value < 0.0 || rate_value > 1e50 || std::isnan(rate_value) || std::isinf(rate_value)) {
                         std::cout << "Rate value: " << rate_value << std::endl;
@@ -228,7 +228,7 @@ RateValues ElectronPhonon::compute_hole_phonon_rate(int idx_n1, std::size_t idx_
                                         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_e;
 
                     if (rate_value < 0.0 || rate_value > 1e50 || std::isnan(rate_value) || std::isinf(rate_value)) {
                         std::cout << "Rate value: " << rate_value << std::endl;

diff --git a/src/BZ_MESH/impact_ionization.cpp b/src/BZ_MESH/impact_ionization.cpp
@@ -49,27 +49,27 @@ void ImpactIonization::read_dielectric_file(const std::string& filename) {
 }
 
 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;
+    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);
+                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;
@@ -127,26 +127,111 @@ void ImpactIonization::compute_eigenstates(int nb_threads) {
     std::cout << "Number of BZ states: " << m_list_BZ_states.size() << std::endl;
 }
 
-double ImpactIonization::compute_impact_ionization_rate(int idx_n1, const Vector3D<double>& k1) {
-    constexpr int nb_valence_bands = 3;
-    constexpr int nb_conduction_bands = 4;
-
-    std::size_t nb_vtx = m_list_BZ_states[0]->get_list_vertices().size();
+double ImpactIonization::compute_impact_ionization_rate(int idx_n1, std::size_t idx_k1) {
+    constexpr int                     nb_valence_bands    = 3;
+    constexpr int                     nb_conduction_bands = 4;
+    constexpr int                     min_conduction_band = 4;
+    const std::size_t                 nb_vtx              = m_list_BZ_states[0]->get_list_vertices().size();
+    const std::vector<Vector3D<int>>& basis_vector_G      = m_list_BZ_states[0]->get_basis_vectors();
+
+    std::cout << "Start computing impact ionization rate for band " << idx_n1 << " and k-point " << idx_k1 << std::endl;
+    auto start_precompute = std::chrono::high_resolution_clock::now();
+    std::cout << "Nb of Bz states: " << m_list_BZ_states.size() << std::endl;
+    std::cout << "Nb cols: " << m_list_BZ_states[0]->get_eigen_states()[0].cols() << std::endl;
+
+    // Sum_2_prime[idx_band][idx_node]
+    std::vector<std::vector<complex_d>> Sum_2_prime(nb_conduction_bands);
+    for (int idx_n2_prime = 0; idx_n2_prime < nb_conduction_bands; ++idx_n2_prime) {
+        int n2_prime = idx_n2_prime + min_conduction_band;
+        Sum_2_prime[idx_n2_prime].resize(nb_vtx);
+        for (std::size_t idx_node = 0; idx_node < nb_vtx; ++idx_node) {
+            const Eigen::MatrixXcd& A_2_prime   = m_list_BZ_states[0]->get_eigen_states()[idx_node];
+            Sum_2_prime[idx_n2_prime][idx_node] = A_2_prime.col(n2_prime).sum();
+        }
+    }
+    std::cout << "Sum_2_prime done" << std::endl;
+    // Sum_2[idx_band][idx_node]
+    std::vector<std::vector<complex_d>> Sum_2(nb_valence_bands);
+    for (int idx_n2 = 0; idx_n2 < nb_valence_bands; ++idx_n2) {
+        Sum_2_prime[idx_n2].resize(nb_vtx);
+        for (std::size_t idx_node = 0; idx_node < nb_vtx; ++idx_node) {
+            const Eigen::MatrixXcd& A_2 = m_list_BZ_states[0]->get_eigen_states()[idx_node];
+            Sum_2_prime[idx_n2][idx_node] = A_2.col(idx_n2).sum();
+        }
+    }
 
-    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) {
-
+    std::cout << "Sum_2 done" << std::endl;
+
+    // Sum_2_prime[idx_band][idx_node] (n1, k1 are in an outter loop)
+    std::vector<std::vector<complex_d>> Sum_1_prime_1(nb_conduction_bands);
+    for (int idx_n1_prime = 0; idx_n1_prime < nb_conduction_bands; ++idx_n1_prime) {
+        std::size_t n1_prime = idx_n1_prime + min_conduction_band;
+        Sum_1_prime_1[idx_n1_prime].resize(nb_vtx);
+        for (std::size_t idx_k1_prime = 0; idx_k1_prime < nb_vtx; ++idx_k1_prime) {
+            const Eigen::MatrixXcd& A_1_prime = m_list_BZ_states[0]->get_eigen_states()[idx_k1_prime];
+            const Eigen::MatrixXcd& A_1       = m_list_BZ_states[0]->get_eigen_states()[idx_k1];
+            int                     nb_Gvect  = A_1_prime.rows();
+            complex_d               sum       = 0.0;
+            for (int idx_G1_prime = 0; idx_G1_prime < nb_Gvect; ++idx_G1_prime) {
+                for (int idx_G1 = 0; idx_G1 < nb_Gvect; ++idx_G1) {
+                    auto   G1_prime = basis_vector_G[idx_G1_prime];
+                    auto   G1       = basis_vector_G[idx_G1];
+                    auto   k1       = m_list_BZ_states[0]->get_vertex_position(idx_k1);
+                    auto   k1_prime = m_list_BZ_states[0]->get_vertex_position(idx_k1_prime);
+                    auto   q_a      = k1 - k1_prime + G1 - G1_prime;
+                    double energy_w = m_list_BZ_states[0]->get_energies()[idx_n1] - m_list_BZ_states[0]->get_energies()[n1_prime];
+                    // complex_d epsilon    = m_dielectric_mesh.interpolate_dielectric_function(q_a, energy_w);
+                    complex_d epsilon    = 1.0;
+                    complex_d factor_eps = EmpiricalPseudopotential::Constants::q_e * EmpiricalPseudopotential::Constants::q_e /
+                                           (EmpiricalPseudopotential::Constants::eps_zero * epsilon * q_a.norm() * q_a.norm());
+                    sum += std::conj(A_1_prime(idx_G1_prime, n1_prime)) * A_1(idx_G1, idx_n1) * factor_eps;
                 }
-            }   
+            }
+            Sum_1_prime_1[idx_n1_prime][idx_k1_prime] = sum;
         }
     }
-
-
-
-
 
+    auto end_precompute = std::chrono::high_resolution_clock::now();
+    std::cout << "Precompute time: " << std::chrono::duration_cast<std::chrono::milliseconds>(end_precompute - start_precompute).count()
+              << " ms" << std::endl;
+    std::cout << "DONE PRECOMPUTE\n Start computing matrix element" << std::endl;
+
+    auto start_compute = std::chrono::high_resolution_clock::now();
+
+    const std::vector<Vertex>& list_vertices = m_list_BZ_states[0]->get_list_vertices();
+
+    for (int idx_n2 = 0; idx_n2 < nb_valence_bands; ++idx_n2) {
+        for (int idx_n1_prime = 0; idx_n1_prime < nb_conduction_bands; ++idx_n1_prime) {
+            std::size_t n1_prime = idx_n1_prime + min_conduction_band;
+            for (int idx_n2_prime = 0; idx_n2_prime < nb_conduction_bands; ++idx_n2_prime) {
+                std::size_t n2_prime = idx_n2_prime + min_conduction_band;
+                for (std::size_t idx_k1_prime = 0; idx_k1_prime < nb_vtx; ++idx_k1_prime) {
+                    std::cout << "idx_n1: " << idx_n1 << " idx_n1_prime: " << idx_n1_prime << " idx_n2: " << idx_n2 << " idx_n2_prime: " << idx_n2_prime
+                              << " idx_k1: " << idx_k1 << " idx_k1_prime: " << idx_k1_prime << std::endl;
+                    std::vector<double> FullMatrixElement(nb_vtx);
+                    for (std::size_t idx_k2_prime = 0; idx_k2_prime < nb_vtx; ++idx_k2_prime) {
+                        // std::cout << "\r" << idx_k1_prime << " / " << nb_vtx << " --> " << idx_k2_prime << " / " << nb_vtx << std::flush;
+                        vector3 k_2_momentum = list_vertices[idx_k1].get_position() - list_vertices[idx_k1_prime].get_position() -
+                                               list_vertices[idx_k2_prime].get_position();
+                        if (k_2_momentum.norm() > m_max_radius_G0_BZ || k_2_momentum.norm() < 1e-12) {
+                            continue;
+                        }
+                        std::size_t idx_k2 = 18;
+                        complex_d   Ma =
+                            Sum_2_prime[idx_n2_prime][idx_k2_prime] * Sum_1_prime_1[idx_n1_prime][idx_k1_prime] * Sum_2[idx_n2][idx_k2];
+                        complex_d Mb =
+                            Sum_2_prime[idx_n1_prime][idx_k1_prime] * Sum_1_prime_1[idx_n2_prime][idx_k2_prime] * Sum_2[idx_n2][idx_k2];
+                        FullMatrixElement[idx_k2_prime] =
+                            std::abs(Ma) * std::abs(Ma) + std::abs(Mb) * std::abs(Mb) + std::abs(Ma - Mb) * std::abs(Ma - Mb);
+                    }
+                }
+            }
+        }
+    }
+    std::cout << std::endl;
+    auto end_compute = std::chrono::high_resolution_clock::now();
+    std::cout << "Compute time: " << std::chrono::duration_cast<std::chrono::milliseconds>(end_compute - start_compute).count() << " ms" << std::endl;
+    return 0.0;
 }
 
 }  // namespace bz_mesh
diff --git a/src/BZ_MESH/impact_ionization.hpp b/src/BZ_MESH/impact_ionization.hpp
@@ -82,7 +82,7 @@ class ImpactIonization {
                                                                                       int idx_k2,
                                                                                       int idx_k2_prime) const;
 
-    double compute_impact_ionization_rate(int idx_n1, const Vector3D<double>& k1);
+    double compute_impact_ionization_rate(int idx_n1, std::size_t idx_k1);
 };
 
 }  // namespace bz_mesh
diff --git a/src/EPP/Constants.hpp b/src/EPP/Constants.hpp
@@ -22,7 +22,7 @@ constexpr double h              = 6.62606957e-34;
 constexpr double h_bar          = 6.62606957e-34 / (2 * M_PI);
 constexpr double h_bar_eV       = 6.582119514e-16;
 constexpr double m0             = 9.10938356e-31;
-constexpr double q              = 1.602176634e-19;
+constexpr double q_e              = 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;

diff --git a/src/EPP/Hamiltonian.cpp b/src/EPP/Hamiltonian.cpp
@@ -50,7 +50,7 @@ void Hamiltonian::SetMatrix(const Vector3D<double>& k, bool add_non_local_correc
     const double       fourier_factor = 2.0 * M_PI / latticeConstant;
     matrix                            = m_constant_non_diagonal_matrix;
     // diagonal elements
-    const double diag_factor = pow(Constants::h_bar, 2) / (2.0 * Constants::m0 * Constants::q);
+    const double diag_factor = pow(Constants::h_bar, 2) / (2.0 * Constants::m0 * Constants::q_e);
     for (unsigned int i = 0; i < basisSize; ++i) {
         Vector3D<double> real_k_vector = (k + m_basisVectors[i]);
         const double     KG2           = fourier_factor * fourier_factor * diag_factor * (real_k_vector * real_k_vector);

diff --git a/src/EPP/Material.cpp b/src/EPP/Material.cpp
@@ -139,7 +139,7 @@ double F_2_function_gaussian(const Vector3D<double>& K1, const Vector3D<double>&
 std::complex<double> Material::compute_pseudopotential_non_local_correction(const Vector3D<double>& K1_normalized,
                                                                             const Vector3D<double>& K2_normalized,
                                                                             const Vector3D<double>& tau) const {
-    const double           diag_factor       = pow(Constants::h_bar, 2) / (2.0 * Constants::m0 * Constants::q);
+    const double           diag_factor       = pow(Constants::h_bar, 2) / (2.0 * Constants::m0 * Constants::q_e);
     const double           fourier_factor    = 2.0 * M_PI / get_lattice_constant_meter();
     const Vector3D<double> G_diff_normalized = (K1_normalized - K2_normalized);
     const Vector3D<double> K1                = K1_normalized * fourier_factor;

diff --git a/src/EPP/NonLocalFunctional.cpp b/src/EPP/NonLocalFunctional.cpp
@@ -27,7 +27,7 @@
 //     : m_non_local_parameters(non_local_parameters),
 //       m_material(material),
 //       m_tau(tau),
-//       m_cinetic_factor{(Constants::h_bar * Constants::h_bar) / (2.0 * Constants::m0 * Constants::q)},
+//       m_cinetic_factor{(Constants::h_bar * Constants::h_bar) / (2.0 * Constants::m0 * Constants::q_e)},
 //       m_fourrier_factor{2.0 * M_PI / m_material.get_lattice_constant_meter()},
 //       m_V0_pref_factor(2.0 * M_PI / m_material.get_lattice_constant_meter()),
 //       m_V2_square_well_pref_factor{m_V0_pref_factor},
@@ -101,7 +101,7 @@
 // std::complex<double> Material::compute_pseudopotential_non_local_correction(const Vector3D<double>& K1_normalized,
 //                                                                             const Vector3D<double>& K2_normalized,
 //                                                                             const Vector3D<double>& tau) const {
-//     const double           diag_factor       = pow(Constants::h_bar, 2) / (2.0 * Constants::m0 * Constants::q);
+//     const double           diag_factor       = pow(Constants::h_bar, 2) / (2.0 * Constants::m0 * Constants::q_e);
 //     const double           fourier_factor    = 2.0 * M_PI / get_lattice_constant_meter();
 //     const Vector3D<double> G_diff_normalized = (K1_normalized - K2_normalized);
 //     const Vector3D<double> K1                = K1_normalized * fourier_factor;