Continue SOC implementation

apps/fullstates.cpp

@@ -69,19 +69,36 @@ int main(int argc, char *argv[]) {
     my_bz_mesh.set_basis_vectors(basis);
 
     my_bz_mesh.read_mesh_geometry_from_msh_file(arg_mesh_file.getValue());
+    std::cout << "Mesh volume: " << my_bz_mesh.compute_mesh_volume() << std::endl;
+
+
     my_bz_mesh.compute_eigenstates(my_options.nrThreads);
 
-    Vector3D<double> q_shift = Vector3D<double>{1e-12, 0.0, 0.0};
+    Vector3D<double> q_shift = Vector3D<double>{1e-10, 0.0, 0.0};
     my_bz_mesh.compute_shifted_eigenstates(q_shift, my_options.nrThreads);
+    std::cout << "\n\n" << std::endl;
 
-    std::cout << "Mesh volume: " << my_bz_mesh.compute_mesh_volume() << std::endl;
 
     auto end = std::chrono::high_resolution_clock::now();
 
     std::chrono::duration<double> elapsed_seconds = end - start;
     std::cout << "Elapsed time: " << elapsed_seconds.count() << "s" << std::endl;
 
-    my_bz_mesh.export_full_eigenstates();
-
-    return 0;
+    std::vector<double> list_energy;
+    double min_energy   = 0.0;
+    double max_energy   = 20.0;
+    double energy_step  = 0.01;
+    for (double energy = min_energy; energy <= max_energy + energy_step; energy += energy_step) {
+        list_energy.push_back(energy);
+    }
+    double eta_smearing = 0.01;
+    std::cout << "Number of energies to compute: " << list_energy.size() << std::endl;
+    auto start2 = std::chrono::high_resolution_clock::now();
+    my_bz_mesh.compute_dielectric_function(list_energy, eta_smearing, my_options.nrThreads);
+    my_bz_mesh.export_dielectric_function("./TEST_DIELECTRIC_FUNCTION_");
+    auto end2 = std::chrono::high_resolution_clock::now();
+    std::chrono::duration<double> elapsed_seconds2 = end2 - start2;
+    std::cout << "Time Dielectric Function : " << elapsed_seconds2.count() << "s" << std::endl;
+
+        return 0;
 }

python/plots/plot_eps_vs_energy.py

@@ -20,7 +20,14 @@
 
 
 def plot_dielectric_function_vs_energy(filename: str):
-    energies, eps_r, eps_i = np.loadtxt(filename, unpack=True, skiprows=1, delimiter=',')
+    # energies, eps_r, eps_i = np.loadtxt(filename, unpack=True, skiprows=1, delimiter=',')
+    data = np.loadtxt(filename, unpack=True, skiprows=1, delimiter=',')
+    energies = data[0]
+    eps_r = data[1]
+    if len(data) == 3:
+        eps_i = data[2]
+    else:
+        eps_i = np.zeros_like(eps_r)
     fig, ax = plt.subplots()
     ax.plot(energies, eps_r, label="Real", c='b')
     ax.plot(energies, eps_i, "-", label="Imaginary", c='r')

src/BZ_MESH/bz_mesh.cpp

@@ -53,7 +53,11 @@ void MeshBZ::read_mesh_geometry_from_msh_file(const std::string& filename) {
 
     m_list_vertices.reserve(size_nodes_tags);
     double lattice_constant     = m_material.get_lattice_constant_meter();
-    double normalization_factor = 2.0 * M_PI / lattice_constant;
+    std::cout << "Lattice const: " << lattice_constant << std::endl;
+    std::cout << "V: " << std::pow(2.0 * M_PI, 3) / std::pow(lattice_constant, 3.0) << std::endl;
+
+    double normalization_factor =  2.0 * M_PI / lattice_constant;
+    // double normalization_factor =  1.0;
     for (std::size_t index_vertex = 0; index_vertex < size_nodes_tags; ++index_vertex) {
         m_list_vertices.push_back(Vertex(index_vertex,
                                          normalization_factor * nodeCoords[3 * index_vertex],
@@ -164,9 +168,10 @@ void MeshBZ::compute_energy_gradient_at_tetras() {
 double MeshBZ::compute_mesh_volume() const {
     double total_volume = 0.0;
     for (auto&& tetra : m_list_tetrahedra) {
-        total_volume += fabs(tetra.get_signed_volume());
+        total_volume += std::fabs(tetra.get_signed_volume());
+        // std::cout << "Tetra " << tetra.get_index() << " volume: " << tetra.get_signed_volume() << std::endl;
     }
-    total_volume *= (1.0 / pow(2.0 * M_PI, 3.0));
+    // total_volume *= (1.0 / pow(2.0 * M_PI, 3.0));    
     return total_volume;
 }
 

src/BZ_MESH/bz_mesh.hpp

@@ -93,7 +93,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_overlapp_integral_impact_ionization_electrons(double energy);
+    double compute_overlap_integral_impact_ionization_electrons(double energy);
 
     std::vector<std::vector<double>> compute_dos_band_at_band(int         band_index,
                                                               double      min_energy,

src/BZ_MESH/bz_states.cpp

@@ -49,6 +49,7 @@ void BZ_States::compute_eigenstates(int nb_threads) {
 }
 
 void BZ_States::compute_shifted_eigenstates(const Vector3D<double>& q_shift, int nb_threads) {
+    m_q_shift                       = q_shift;
     double     normalization_factor = 2.0 * M_PI / m_material.get_lattice_constant_meter();
     const bool m_nonlocal_epm       = false;
     const bool keep_eigenvectors    = true;
@@ -64,8 +65,8 @@ void BZ_States::compute_shifted_eigenstates(const Vector3D<double>& q_shift, int
         auto k_point = Vector3D<double>(m_list_vertices[idx_k].get_position().x(),
                                         m_list_vertices[idx_k].get_position().y(),
                                         m_list_vertices[idx_k].get_position().z());
+        k_point      = k_point * 1.0 / normalization_factor;
         k_point += q_shift;
-        k_point         = k_point * 1.0 / normalization_factor;
         auto idx_thread = omp_get_thread_num();
         hamiltonian_per_thread[idx_thread].SetMatrix(k_point, m_nonlocal_epm);
         hamiltonian_per_thread[idx_thread].Diagonalize(keep_eigenvectors);
@@ -76,6 +77,95 @@ void BZ_States::compute_shifted_eigenstates(const Vector3D<double>& q_shift, int
     }
 }
 
+/**
+ * @brief Compute the dielectric function for a given list of energies and a given smearing.
+ * The integration is performed by summing the contribution of each tetrahedron to the dielectric function.
+ *
+ * @param energies
+ * @param eta_smearing
+ * @param nb_threads
+ */
+void BZ_States::compute_dielectric_function(const std::vector<double>& list_energies, double eta_smearing, int nb_threads) {
+    m_list_energies                         = list_energies;
+    const int   index_first_conduction_band = 4;
+    std::size_t nb_tetra                    = m_list_tetrahedra.size();
+    m_dielectric_function_real.resize(list_energies.size());
+    constexpr double one_fourth = 1.0 / 4.0;
+
+    std::vector<double> dielectric_function_real_at_energies(list_energies.size(), 0.0);
+
+#pragma omp parallel for schedule(dynamic) num_threads(nb_threads)
+    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<std::size_t, 4>    list_idx_vertices = m_list_tetrahedra[idx_tetra].get_list_indices_vertices();
+        const std::array<Vertex*, 4>& list_vertices     = m_list_tetrahedra[idx_tetra].get_list_vertices();
+        double                        volume_tetra      = std::fabs(m_list_tetrahedra[idx_tetra].compute_signed_volume());
+        // std::cout << "Volume tetra = " << volume_tetra << std::endl;
+        std::vector<double> sum_dielectric_function_real_tetra_at_energies(list_energies.size(), 0.0);
+        // Loop over the vertices of the tetrahedron
+        for (std::size_t idx_vertex = 0; idx_vertex < 4; ++idx_vertex) {
+            std::size_t index_k = list_idx_vertices[idx_vertex];
+            for (int idx_conduction_band = index_first_conduction_band; idx_conduction_band < m_nb_bands; ++idx_conduction_band) {
+                for (int idx_valence_band = 0; idx_valence_band < index_first_conduction_band; ++idx_valence_band) {
+                    double overlap_integral = pow(std::fabs(m_eigenvectors_k_plus_q[index_k]
+                                                                .col(idx_conduction_band)
+                                                                .adjoint()
+                                                                .dot(m_eigenvectors_k[index_k].col(idx_valence_band))),
+                                                  2);
+                    double delta_energy = m_eigenvalues_k_plus_q[index_k][idx_conduction_band] - m_eigenvalues_k[index_k][idx_valence_band];
+                    for (std::size_t index_energy = 0; index_energy < list_energies.size(); ++index_energy) {
+                        double energy = list_energies[index_energy];
+                        double factor_1 =
+                            (delta_energy - energy) / ((delta_energy - energy) * (delta_energy - energy) + eta_smearing * eta_smearing);
+                        double factor_2 =
+                            (delta_energy + energy) / ((delta_energy + energy) * (delta_energy + energy) + eta_smearing * eta_smearing);
+                        double total_factor = factor_1 + factor_2;
+                        sum_dielectric_function_real_tetra_at_energies[index_energy] += overlap_integral * total_factor;
+                    }
+                }
+            }
+        }
+        for (std::size_t index_energy = 0; index_energy < list_energies.size(); ++index_energy) {
+            sum_dielectric_function_real_tetra_at_energies[index_energy] *= volume_tetra * one_fourth;
+            dielectric_function_real_at_energies[index_energy] += sum_dielectric_function_real_tetra_at_energies[index_energy];
+        }
+    }
+
+    double q_squared = m_q_shift.Length() * m_q_shift.Length();
+    std::cout << "\nq_squared = " << q_squared << std::endl;
+    // double normalization_1 = (EmpiricalPseudopotential::Constants::q * EmpiricalPseudopotential::Constants::q * 4.0 * M_PI);
+    double normalization_1 = (M_PI) / (2.0 / std::pow(2.0 * M_PI, 3));
+    // double normalization_2 = compute_mesh_volume();
+    double normalization_3 = (EmpiricalPseudopotential::Constants::q * EmpiricalPseudopotential::Constants::q);
+    std::cout << std::endl;
+    std::cout << "Normalization factor 1: " << normalization_1 << std::endl;
+    std::cout << "Normalization factor 2: " << normalization_3 << std::endl;
+    std::cout << "Normalization factor q: " << 1.0 / q_squared << std::endl;
+    for (std::size_t index_energy = 0; index_energy < list_energies.size(); ++index_energy) {
+        m_dielectric_function_real[index_energy] = dielectric_function_real_at_energies[index_energy];
+        // std::cout << "Energy = " << list_energies[index_energy] << " eV, dielectric function = " <<
+        // m_dielectric_function_real[index_energy] << std::endl; m_dielectric_function_real[index_energy] *= normalization_1;
+        m_dielectric_function_real[index_energy] *= normalization_3;
+        // m_dielectric_function_real[index_energy] /= normalization_2;
+        m_dielectric_function_real[index_energy] /= q_squared;
+        m_dielectric_function_real[index_energy] += 1.0;
+    }
+    std::cout << "E = 0 --> " << dielectric_function_real_at_energies[0] << std::endl;
+    std::cout << "E = 0 --> " << dielectric_function_real_at_energies[0] / q_squared << std::endl;
+    std::cout << "E = 0 --> " << dielectric_function_real_at_energies[0] / (q_squared * compute_mesh_volume()) << std::endl;
+    std::cout << "E = 0 --> " << m_dielectric_function_real[0] << std::endl;
+}
+
+// Export the dielectric function to a file in the format (energy, dielectric function) (csv format).
+void BZ_States::export_dielectric_function(const std::string& prefix) const {
+    std::ofstream dielectric_function_file(prefix + "_dielectric_function.csv");
+    dielectric_function_file << "Energy (eV), Dielectric function" << std::endl;
+    for (std::size_t index_energy = 0; index_energy < m_dielectric_function_real.size(); ++index_energy) {
+        dielectric_function_file << m_list_energies[index_energy] << ", " << m_dielectric_function_real[index_energy] << std::endl;
+    }
+    dielectric_function_file.close();
+}
+
 void BZ_States::export_full_eigenstates() const {
     std::filesystem::remove_all("eigenstates");
     std::filesystem::create_directory("eigenstates");

src/BZ_MESH/bz_states.hpp

@@ -29,6 +29,16 @@ class BZ_States : public MeshBZ {
     std::vector<Eigen::MatrixXcd> m_eigenvectors_k;
     std::vector<Eigen::MatrixXcd> m_eigenvectors_k_plus_q;
 
+    Vector3D<double>    m_q_shift;
+    std::vector<double> m_list_energies;
+
+    /**
+     * @brief Real part of the dielectric function.
+     * m_dielectric_function_real[idx_energy] is the real part of the dielectric function at the energy m_energies[idx_energy].
+     *
+     */
+    std::vector<double> m_dielectric_function_real;
+
  public:
     BZ_States(const EmpiricalPseudopotential::Material& material) : MeshBZ(material) {}
 
@@ -37,7 +47,11 @@ class BZ_States : public MeshBZ {
     void compute_eigenstates(int nb_threads = 1);
     void compute_shifted_eigenstates(const Vector3D<double>& q_shift, int nb_threads = 1);
 
+    const std::vector<double>& get_energies() const { return m_list_energies; }
+    void                       set_energies(const std::vector<double>& energies) { m_list_energies = energies; }
+
     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;
 
     double compute_direct_impact_ionization_matrix_element(int                     idx_n1,
                                                            int                     idx_n1_prime,

src/BZ_MESH/mesh_tetra.cpp

@@ -79,6 +79,8 @@ void Tetra::compute_min_max_energies_at_bands() {
  * @return double
  */
 double Tetra::compute_signed_volume() const {
+    // std::cout << m_list_edges[0] << std::endl;
+    // std::cout << m_list_edges[1] << std::endl;
     return (1.0 / 6.0) * scalar_triple_product(m_list_edges[0], m_list_edges[1], m_list_edges[2]);
 }
 

src/BZ_MESH/mesh_tetra.hpp

@@ -48,7 +48,7 @@ class Tetra {
      * The sign depends on the "orientation" of the tetrahedra.
      *
      */
-    double m_signed_volume;
+    double m_signed_volume = 0.0;
 
     /**
      * @brief Number of conduction bands.
@@ -84,6 +84,17 @@ class Tetra {
     Tetra(std::size_t index, const std::array<Vertex*, 4>& list_vertices);
     void compute_min_max_energies_at_bands();
 
+    std::size_t                   get_index() const { return m_index; }
+    const std::array<Vertex*, 4>& get_list_vertices() const { return m_list_vertices; }
+    std::array<std::size_t, 4>    get_list_indices_vertices() const {
+           return {m_list_vertices[0]->get_index(),
+                   m_list_vertices[1]->get_index(),
+                   m_list_vertices[2]->get_index(),
+                   m_list_vertices[3]->get_index()};
+    }
+    std::array<vector3, 6> get_list_edges() const { return m_list_edges; }
+    std::size_t            get_nb_bands() const { return m_nb_bands; }
+
     std::array<double, 4> get_band_energies_at_vertices(std::size_t index_band) const;
 
     double  compute_signed_volume() const;

src/EPP/Constants.hpp

@@ -21,6 +21,7 @@ constexpr double Ryd_to_eV      = 13.6;
 constexpr double h              = 6.62606957e-34;
 constexpr double h_bar          = 6.62606957e-34 / (2 * M_PI);
 constexpr double m0             = 9.10938356e-31;
+// constexpr double q              = 1.602176634e-19;
 constexpr double q              = 1.6e-19;
 constexpr double eps_zero       = 8.854187e-12;
 constexpr double k_b            = 1.38064852e-23;

src/EPP/DielectricFunction.cpp

@@ -71,7 +71,7 @@ void DielectricFunction::generate_k_points_grid(std::size_t Nx, std::size_t Ny,
 }
 
 /**
- * @brief Compute the energy and wavevector dependent dielectric function.
+ * @brief Compute the energy and wave vector dependent dielectric function.
  * The formula used is the one from the paper "
  *
  * @param eta_smearing