Up

apps/CMakeLists.txt

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

apps/fullstates.cpp

@@ -0,0 +1,84 @@
+/**
+ * @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;
+    EmpiricalPseudopotential::Materials materials;
+    std::string                         file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-local.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);
+    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());
+    my_bz_mesh.compute_eigenstates(my_options.nrThreads);
+
+    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;
+}

src/BZ_MESH/CMakeLists.txt

@@ -3,11 +3,13 @@ set(HEADER_FILES_LIBBZMESH
   mesh_vertex.hpp
   bz_mesh.hpp
   vector.hpp
-  iso_triangle.hpp)
+  iso_triangle.hpp
+  bz_states.hpp)
 
 set(SOURCE_FILES_LIBBZMESH
   mesh_tetra.cpp
-  bz_mesh.cpp)
+  bz_mesh.cpp
+  bz_states.cpp)
 
 add_library(lib_bzmesh STATIC ${SOURCE_FILES_LIBBZMESH} ${HEADER_FILES_LIBBZMESH})
 

src/BZ_MESH/bz_mesh.hpp

@@ -21,7 +21,7 @@ namespace bz_mesh {
 enum class BandType { valence, conduction };
 
 class MeshBZ {
- private:
+ protected:
     /**
      * @brief The material of the Brillouin zone.
      *

src/BZ_MESH/bz_states.cpp

@@ -0,0 +1,67 @@
+/**
+ * @file bz_states.cpp
+ * @author remzerrr (remi.helleboid@gmail.com)
+ * @brief
+ * @version 0.1
+ * @date 2022-12-20
+ *
+ * @copyright Copyright (c) 2022
+ *
+ */
+
+#include "bz_states.hpp"
+
+#include <Eigen/Dense>
+#include <filesystem>
+#include <fstream>
+
+#include "Hamiltonian.h"
+#include "bz_mesh.hpp"
+#include "omp.h"
+
+namespace bz_mesh {
+
+void BZ_States::compute_eigenstates(int nb_threads) {
+    double     normalization_factor = 2.0 * M_PI / m_material.get_lattice_constant_meter();
+    const bool m_nonlocal_epm       = false;
+    const bool keep_eigenvectors    = true;
+    m_eigenvalues_k.resize(m_list_vertices.size());
+    m_eigenvectors_k.resize(m_list_vertices.size());
+    std::vector<EmpiricalPseudopotential::Hamiltonian> hamiltonian_per_thread;
+    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;
+        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;
+        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);
+        m_eigenvalues_k[idx_k]  = hamiltonian_per_thread[idx_thread].eigenvalues();
+        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);
+    }
+}
+
+void BZ_States::export_full_eigenstates() const {
+    std::filesystem::remove_all("eigenstates");
+    std::filesystem::create_directory("eigenstates");
+    std::filesystem::create_directory("eigenstates/eigenvectors");
+    std::filesystem::create_directory("eigenstates/eigenvalues");
+    for (std::size_t idx_k = 0; idx_k < m_list_vertices.size(); ++idx_k) {
+        std::ofstream eigenvalues_file("eigenstates/eigenvalues/eigenvalues_" + std::to_string(idx_k) + ".txt");
+        eigenvalues_file << m_eigenvalues_k[idx_k].transpose() << std::endl;
+        eigenvalues_file.close();
+
+        std::ofstream eigenvectors_file("eigenstates/eigenvectors/eigenvectors_" + std::to_string(idx_k) + ".txt");
+        eigenvectors_file << m_eigenvectors_k[idx_k] << std::endl;
+        eigenvectors_file.close();
+    }
+}
+
+}  // namespace bz_mesh

src/BZ_MESH/bz_states.hpp

@@ -0,0 +1,42 @@
+/**
+ * @file bz_dielectric.hpp
+ * @author remzerrr (remi.helleboid@gmail.com)
+ * @brief
+ * @version 0.1
+ * @date 2022-12-20
+ *
+ * @copyright Copyright (c) 2022
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+
+#include "bz_mesh.hpp"
+#include "Material.h"
+
+namespace bz_mesh {
+
+class BZ_States : public MeshBZ {
+ private:
+    int m_nb_bands = 0;
+
+    std::vector<Vector3D<int>> m_basisVectors;
+
+    std::vector<Eigen::VectorXd>  m_eigenvalues_k;
+
+    std::vector<Eigen::MatrixXcd> m_eigenvectors_k;
+
+ public:
+    BZ_States(const EmpiricalPseudopotential::Material& material) : MeshBZ(material) {}
+
+    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; }
+
+    void compute_eigenstates(int nb_threads = 1);
+
+    void export_full_eigenstates() const;
+};
+
+}  // namespace bz_mesh