Start construction ImpactIonization class

apps/impact_io.cpp

@@ -22,6 +22,7 @@
 #include "bz_mesh.hpp"
 #include "bz_meshfile.hpp"
 #include "bz_states.hpp"
+#include "impact_ionization.hpp"
 
 int main(int argc, char *argv[]) {
     TCLAP::CmdLine               cmd("EPP PROGRAM. COMPUTE BAND STRUCTURE ON A BZ MESH.", ' ', "1.0");
@@ -30,13 +31,20 @@ int main(int argc, char *argv[]) {
     TCLAP::ValueArg<int>         arg_nb_energies("e", "nenergy", "Number of energies to compute", false, 250, "int");
     TCLAP::ValueArg<int>         arg_nb_bands("b", "nbands", "Number of bands to consider", false, 12, "int");
     TCLAP::ValueArg<int>         arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int");
+    TCLAP::ValueArg<double>      arg_radius_BZ("R",
+                                           "radiusBZ",
+                                           "Max norm of G vector for which the BZ center in G is taken into account",
+                                           false,
+                                           0,
+                                           "double");
     TCLAP::SwitchArg plot_with_python("P", "plot", "Call a python script after the computation to plot the band structure.", false);
     cmd.add(plot_with_python);
     cmd.add(arg_mesh_file);
     cmd.add(arg_material);
     cmd.add(arg_nb_bands);
     cmd.add(arg_nb_energies);
     cmd.add(arg_nb_threads);
+    cmd.add(arg_radius_BZ);
 
     cmd.parse(argc, argv);
 
@@ -60,26 +68,14 @@ int main(int argc, char *argv[]) {
 
     EmpiricalPseudopotential::Material current_material = materials.materials.at(arg_material.getValue());
 
-    bz_mesh::BZ_States my_bz_mesh(current_material);
-    my_bz_mesh.set_nb_bands(nb_bands_to_use);
-    EmpiricalPseudopotential::BandStructure band_structure{};
-
-    int nb_nearest_neighbors = 10;
-    band_structure.Initialize(current_material, nb_bands_to_use, {}, nb_nearest_neighbors, nonlocal_epm, enable_soc);
-    auto basis = band_structure.get_basis_vectors();
-    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);
-
+    bz_mesh::ImpactIonization my_impact_ionization(current_material, arg_mesh_file.getValue());
+    my_impact_ionization.set_max_radius_G0_BZ(arg_radius_BZ.getValue());
+    my_impact_ionization.compute_eigenstates();
 
     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;
 
-
     return 0;
 }

src/BZ_MESH/CMakeLists.txt

@@ -6,6 +6,7 @@ set(HEADER_FILES_LIBBZMESH
   iso_triangle.hpp
   bz_states.hpp
   electron_phonon.hpp
+  impact_ionization.hpp
   )
 
 
@@ -14,6 +15,7 @@ set(SOURCE_FILES_LIBBZMESH
   bz_mesh.cpp
   bz_states.cpp
   electron_phonon.cpp
+  impact_ionization.cpp
 )
 
 add_library(lib_bzmesh STATIC ${SOURCE_FILES_LIBBZMESH} ${HEADER_FILES_LIBBZMESH})

src/BZ_MESH/bz_mesh.cpp

@@ -24,6 +24,14 @@
 
 namespace bz_mesh {
 
+
+void MeshBZ::shift_bz_center(const vector3& center) {
+    m_center = center;
+    for (auto&& vtx : m_list_vertices) {
+        vtx.shift_position(center);
+    }
+}
+
 /**
  * @brief Read the geometry of the mesh from the .msh file: the vertices and the elements are added to
  * the m_list_vertices and m_list_elements lists.

src/BZ_MESH/bz_mesh.hpp

@@ -28,6 +28,12 @@ class MeshBZ {
      */
     EmpiricalPseudopotential::Material m_material;
 
+    /**
+     * @brief G vector on which the BZ is centered.
+     *
+     */
+    vector3 m_center = {0.0, 0.0, 0.0};
+
     /**
      * @brief List of the vertices of the BZ mesh. Each vertices represent a vector k within the Brillouin Zone.
      *
@@ -74,6 +80,9 @@ class MeshBZ {
     MeshBZ(const EmpiricalPseudopotential::Material& material) : m_material(material){};
     MeshBZ(const MeshBZ&) = default;
 
+    vector3 get_center() const { return m_center; }
+    void    shift_bz_center(const vector3& shift);
+
     void read_mesh_geometry_from_msh_file(const std::string& filename);
     void read_mesh_bands_from_msh_file(const std::string& filename);
     void add_new_band_energies_to_vertices(const std::vector<double>& energies_at_vertices);
@@ -82,17 +91,17 @@ class MeshBZ {
 
     std::size_t get_number_vertices() const { return m_list_vertices.size(); }
     std::size_t get_number_elements() const { return m_list_tetrahedra.size(); }
-    double     get_volume() const { return m_total_volume; }
+    double      get_volume() const { return m_total_volume; }
 
-    bool is_inside_mesh_geometry(const vector3& k) const;
-    bool is_inside_mesh_geometry(const Vector3D<double>& k) const;
+    bool    is_inside_mesh_geometry(const vector3& k) const;
+    bool    is_inside_mesh_geometry(const Vector3D<double>& k) const;
     vector3 retrieve_k_inside_mesh_geometry(const vector3& k) const;
 
-    vector3 get_k_at_index(std::size_t index) const { return m_list_vertices[index].get_position(); }
+    vector3     get_k_at_index(std::size_t index) const { return m_list_vertices[index].get_position(); }
     std::size_t get_nearest_k_index(const Vector3D<double>& k) const;
     std::size_t get_nearest_k_index(const vector3& k) const;
 
-    std::size_t get_number_bands() const { return m_min_band.size(); }
+    std::size_t               get_number_bands() const { return m_min_band.size(); }
     std::pair<double, double> get_min_max_energy_at_band(const int& band_index) const {
         return std::make_pair(m_min_band[band_index], m_max_band[band_index]);
     }

src/BZ_MESH/bz_states.cpp

@@ -45,6 +45,7 @@ void BZ_States::compute_eigenstates(int nb_threads) {
         auto nb_rows            = m_eigenvectors_k[idx_k].rows();
         m_eigenvectors_k[idx_k].conservativeResize(nb_rows, m_nb_bands);
     }
+    std::cout << std::endl;
 }
 
 void BZ_States::compute_shifted_eigenstates(const Vector3D<double>& q_shift, int nb_threads) {
@@ -73,6 +74,7 @@ void BZ_States::compute_shifted_eigenstates(const Vector3D<double>& q_shift, int
         auto nb_rows                   = m_eigenvectors_k[idx_k].rows();
         m_eigenvectors_k_plus_q[idx_k].conservativeResize(nb_rows, m_nb_bands);
     }
+    std::cout << std::endl;
 }
 
 /**

src/BZ_MESH/bz_states.hpp

@@ -22,7 +22,6 @@ class BZ_States : public MeshBZ {
  protected:
     int m_nb_bands = 0;
 
-    vector3 m_center = {0.0,0.0,0.0};
 
     std::vector<Vector3D<int>> m_basisVectors;
 
@@ -68,6 +67,7 @@ class BZ_States : public MeshBZ {
 
  public:
     BZ_States(const EmpiricalPseudopotential::Material& material) : MeshBZ(material) {}
+    BZ_States(const BZ_States& other) = default;
 
     void                              set_nb_bands(int nb_bands) { m_nb_bands = nb_bands; }
     void                              set_basis_vectors(const std::vector<Vector3D<int>>& basis_vectors) { m_basisVectors = basis_vectors; }