Fix export of epsilon

apps/CMakeLists.txt

@@ -26,6 +26,11 @@ target_link_libraries(elph.epm PUBLIC libepp lib_bzmesh)
 target_compile_features(elph.epm PRIVATE cxx_std_20)
 target_link_libraries(elph.epm PUBLIC MPI::MPI_CXX)
 
+add_executable(impact_io.epm impact_io.cpp)
+target_link_libraries(impact_io.epm PUBLIC libepp lib_bzmesh)
+target_compile_features(impact_io.epm PRIVATE cxx_std_20)
+target_link_libraries(impact_io.epm PUBLIC MPI::MPI_CXX)
+
 if(USE_MPI_ACCELERATION)
   add_executable(mpiBandsOnBZ mpi_BandsOnBZ.cpp)
   target_link_libraries(mpiBandsOnBZ PUBLIC libepp Eigen3::Eigen)

apps/impact_io.cpp

@@ -0,0 +1,85 @@
+/**
+ * @file fullstates.cpp
+ * @author remzerrr (remi.helleboid@gmail.com)
+ * @brief
+ * @version 0.1
+ * @date 2022-12-20
+ *
+ * @copyright Copyright (c) 2022
+ *
+ */
+
+#include <tclap/CmdLine.h>
+
+#include <chrono>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+
+#include "BandStructure.h"
+#include "Material.h"
+#include "Options.h"
+#include "bz_mesh.hpp"
+#include "bz_meshfile.hpp"
+#include "bz_states.hpp"
+
+int main(int argc, char *argv[]) {
+    TCLAP::CmdLine               cmd("EPP PROGRAM. COMPUTE BAND STRUCTURE ON A BZ MESH.", ' ', "1.0");
+    TCLAP::ValueArg<std::string> arg_mesh_file("f", "meshbandfile", "File with BZ mesh and bands energy.", true, "bz.msh", "string");
+    TCLAP::ValueArg<std::string> arg_material("m", "material", "Symbol of the material to use (Si, Ge, GaAs, ...)", true, "Si", "string");
+    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::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.parse(argc, argv);
+
+    bool                                nonlocal_epm = false;
+    bool                                enable_soc   = false;
+    EmpiricalPseudopotential::Materials materials;
+    std::string                         file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-cohen.yaml";
+    if (nonlocal_epm) {
+        file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials.yaml";
+    }
+    std::cout << "Loading material parameters from " << file_material_parameters << std::endl;
+    materials.load_material_parameters(file_material_parameters);
+
+    Options my_options;
+    my_options.materialName = arg_material.getValue();
+    my_options.nrLevels     = arg_nb_bands.getValue();
+    my_options.nrThreads    = arg_nb_threads.getValue();
+    my_options.print_options();
+    int  nb_bands_to_use = arg_nb_bands.getValue();
+    auto start           = std::chrono::high_resolution_clock::now();
+
+    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);
+
+
+    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;
+}

python/plots/plot_eps_vs_energy.py

@@ -26,11 +26,14 @@ def plot_dielectric_function_vs_energy(dirname):
     list_q = []
     for filename in list_files:
         filestem = Path(filename).stem
-        qz = float(filestem.split('_')[-1])
-        qy = float(filestem.split('_')[-2])
-        qx = float(filestem.split('_')[-3])
-        print(qx, qy, qz)
-        list_q.append(np.sqrt(qx**2 + qy**2 + qz**2))
+        try:
+            qz = float(filestem.split('_')[-1])
+            qy = float(filestem.split('_')[-2])
+            qx = float(filestem.split('_')[-3])
+            print(qx, qy, qz)
+            list_q.append(np.sqrt(qx**2 + qy**2 + qz**2))
+        except:
+            continue
     print(list_q)
 
     list_files = [x for _, x in sorted(zip(list_q, list_files))]

src/BZ_MESH/bz_mesh.cpp

@@ -115,6 +115,7 @@ std::size_t MeshBZ::get_nearest_k_index(const Vector3D<double>& k) const {
     }
     return index_nearest_k;
 }
+
 std::size_t MeshBZ::get_nearest_k_index(const vector3& k) const {
     std::size_t index_nearest_k = 0;
     double      min_distance    = std::numeric_limits<double>::max();

src/BZ_MESH/bz_mesh.hpp

@@ -72,6 +72,7 @@ class MeshBZ {
 
  public:
     MeshBZ(const EmpiricalPseudopotential::Material& material) : m_material(material){};
+    MeshBZ(const MeshBZ&) = default;
 
     void read_mesh_geometry_from_msh_file(const std::string& filename);
     void read_mesh_bands_from_msh_file(const std::string& filename);

src/BZ_MESH/bz_states.hpp

@@ -22,6 +22,8 @@ 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;
 
     std::vector<Eigen::VectorXd>  m_eigenvalues_k;

src/BZ_MESH/impact_ionization.hpp

@@ -13,6 +13,7 @@
 
 #include <Eigen/Dense>
 #include <array>
+#include <memory>
 
 #include "Material.h"
 #include "bz_mesh.hpp"
@@ -21,10 +22,20 @@
 namespace bz_mesh {
 
 typedef std::complex<double> complex_d;
+typedef std::unique_ptr<BZ_States> uptr_BZstates;
 
-class ImpactIonization : private BZ_States {
+class ImpactIonization {
  private:
-    double m_impact_ionization_rate = 0.0;
+    /**
+     * @brief List of Brillouin Zone with their precomputed electronic states.
+     * Each of them correspond to a 1BZ shifted by a vector G of the reciprocal lattice.
+     * 
+     */
+    std::vector<uptr_BZstates> m_list_BZ_states;
+
+    std::vector<double> m_impact_ionization_results;
+
+
 
  public:
     std::array<complex_d, 2> compute_direct_indirect_impact_ionization_matrix_element(int idx_n1,

src/EPP/DielectricFunction.cpp

@@ -255,7 +255,7 @@ void DielectricFunction::export_dielectric_function_at_q(const std::string& file
     if (name_auto) {
         // outname = m_export_prefix + '_' + std::to_string(idx_q) + '_' + std::to_string(m_qpoints[idx_q].X) + "_" + std::to_string(m_qpoints[idx_q].Y) + "_" +
         //           std::to_string(m_qpoints[idx_q].Z) + ".csv";
-        outname = fmt::format("{}_{:05}_{:.1f}_{:.1f}_{:.1f}.csv", m_export_prefix, idx_q, m_qpoints[idx_q].X, m_qpoints[idx_q].Y, m_qpoints[idx_q].Z);
+        outname = fmt::format("{}_{:05}_{:.6f}_{:.6f}_{:.6f}.csv", m_export_prefix, idx_q, m_qpoints[idx_q].X, m_qpoints[idx_q].Y, m_qpoints[idx_q].Z);
     } else {
         outname = filename;
     }