Up

CMakeLists.txt

@@ -221,6 +221,14 @@ if(NOT yaml-cpp_FOUND)
   FetchContent_MakeAvailable(yaml-cpp)
 endif()
 
+FetchContent_Declare(
+  fmt
+  GIT_REPOSITORY "https://github.com/fmtlib/fmt.git"
+  GIT_TAG master)
+FetchContent_MakeAvailable(fmt)
+
+
+
 # MPI SUPPORT
 if(ENABLE_MPI_BUILD)
   message("MPI support activated.")

apps/epsilon.cpp

@@ -46,6 +46,9 @@ std::vector<Vector3D<double>> read_qpoint_dat_file(const std::string& filename)
     std::vector<Vector3D<double>> kpoints;
     std::string                   line;
     while (std::getline(file, line)) {
+        if (line.empty()) {
+            continue;
+        }
         std::stringstream   ss(line);
         std::string         token;
         std::vector<double> kpoint;

python/PopulateMeshWithDielectric.py

@@ -0,0 +1,104 @@
+import gmsh
+import sys
+import numpy as np
+import os, sys
+import pathlib, shutil, glob
+from pathlib import Path
+
+gmsh.initialize(sys.argv)
+
+filename = "EightSphere_3.0.msh"
+
+gmsh.open(filename)
+model_name = gmsh.model.getCurrent()
+
+print(f"model_name = {model_name}")
+
+
+def extract_dielectric_file(filename):
+    energy, eps_r, eps_i = np.loadtxt(filename, delimiter=",", unpack=True, skiprows=1)
+    return energy, eps_r, eps_i
+
+
+
+DIR_OUT = "OUTDIR/"
+list_dielectric_files = glob.glob(DIR_OUT + "Si*.csv")
+print(f"list_dielectric_files = {list_dielectric_files[1]}")
+list_idx_dielectric = [int(Path(file).stem.split("_")[2]) for file in list_dielectric_files]
+
+NbFiles = len(list_dielectric_files)
+
+# Get list of nodes
+nodes = gmsh.model.mesh.getNodes()
+TagNodes = nodes[0]
+NbNodes = TagNodes.shape[0]
+
+if NbFiles != NbNodes:
+    raise ValueError(f"NbFiles = {NbFiles} != NbNodes = {NbNodes}")
+else:
+    print(f"NbFiles = {NbFiles} == NbNodes = {NbNodes}")
+
+# Open the first file to get the energy
+energy, eps_r, eps_i = extract_dielectric_file(list_dielectric_files[0])
+print(f"energy = {energy}")
+LIST_ENERGY = []
+LIST_EPS_R = []
+LIST_EPS_I = []
+
+for file in list_dielectric_files:
+    energy, eps_r, eps_i = extract_dielectric_file(file)
+    LIST_ENERGY.append(energy)
+    LIST_EPS_R.append(eps_r)
+    LIST_EPS_I.append(eps_i)
+
+LIST_ENERGY = np.array(LIST_ENERGY)
+LIST_EPS_R = np.array(LIST_EPS_R)
+LIST_EPS_I = np.array(LIST_EPS_I)
+
+print(f"LIST_ENERGY = {LIST_ENERGY.shape}")
+print(f"LIST_EPS_R = {LIST_EPS_R.shape}")
+print(f"LIST_EPS_I = {LIST_EPS_I.shape}")
+
+energies = LIST_ENERGY[0]
+print(f"energies = {energies}")
+NbEnergies = energies.shape[0]
+
+FileName = "EightSphere_3.0_with_dielectric.msh"
+
+for ix_e in range(NbEnergies//10):
+    print(f"idx_e = {ix_e}")
+    e = energies[ix_e]
+    eps_r = LIST_EPS_R[:, ix_e]
+    eps_i = LIST_EPS_I[:, ix_e]
+    # Add view
+    NameView = f"eps_r_{e}"
+    data_tag = gmsh.view.add(NameView)
+    gmsh.view.addHomogeneousModelData(data_tag, 0, model_name, "NodeData", TagNodes, eps_r)
+    gmsh.view.write(data_tag, FileName, True)
+
+
+
+
+# test_data = np.linalg.norm(Coordinates, axis=1)
+
+# # Add view
+# data_tag = gmsh.view.add("data_test")
+# # gmsh::view::addHomogeneousModelData(data_tag, 0, "mesh_discrete", "NodeData", index_vtx,
+# #                                                     vertex_index_values_of_function.second);
+# gmsh.view.addHomogeneousModelData(data_tag, 0, model_name, "NodeData", TagNodes, test_data)
+
+# # Export 
+# # gmsh::view::write(data_tag, output_path, true);
+# gmsh.view.write(data_tag, "test_data3.msh", True)
+
+# test_data2 = np.sin(np.linalg.norm(Coordinates, axis=1))
+
+# # Add view
+# data_tag2 = gmsh.view.add("data_test_2")
+# # gmsh::view::addHomogeneousModelData(data_tag, 0, "mesh_discrete", "NodeData", index_vtx,
+# #                                                     vertex_index_values_of_function.second);
+# gmsh.view.addHomogeneousModelData(data_tag2, 0, model_name, "NodeData", TagNodes, test_data2)
+
+# # Export 
+# # gmsh::view::write(data_tag, output_path, true);
+# gmsh.view.write(data_tag2, "test_data3.msh", True)

python/sphere_mesh.py

@@ -0,0 +1,50 @@
+import gmsh
+import sys
+import numpy as np
+
+gmsh.initialize(sys.argv)
+
+gmsh.model.add("boolean")
+
+# from http://en.wikipedia.org/wiki/Constructive_solid_geometry
+
+gmsh.option.setNumber("Mesh.Algorithm", 6)
+gmsh.option.setNumber("Mesh.MeshSizeMin", 0.1)
+gmsh.option.setNumber("Mesh.MeshSizeMax", 0.2)
+
+Rt = 1.0
+
+R = 3.0
+Rs = R * .7
+Rt = R 
+
+gmsh.model.occ.addBox(-R, -R, -R, 2 * R, 2 * R, R, 1)
+gmsh.model.occ.addSphere(0, 0, 0, Rt, 2)
+gmsh.model.occ.intersect([(3, 1)], [(3, 2)], 3)
+
+gmsh.model.occ.addBox(-R, -R, -R, 2 * R,  R, 2*R, 4)
+gmsh.model.occ.intersect([(3, 3)], [(3, 4)], 5)
+
+gmsh.model.occ.addBox(-R, -R, -R,  R, 2* R, 2*R, 6)
+gmsh.model.occ.intersect([(3, 5)], [(3, 6)], 7)
+
+
+gmsh.model.occ.synchronize()
+
+gmsh.model.mesh.generate(3)
+
+
+
+# Get list of nodes
+nodes = gmsh.model.mesh.getNodes()
+Coordinates = nodes[1].reshape((nodes[0].shape[0], 3))
+
+filename = f"nodes_eight_sphere_{Rt}.csv"
+np.savetxt(filename, Coordinates, delimiter=" ")
+
+gmsh.write(f"EightSphere_{R}.msh")
+
+if '-nopopup' not in sys.argv:
+    gmsh.fltk.run()
+
+gmsh.finalize()

src/BZ_MESH/CMakeLists.txt

@@ -19,7 +19,7 @@ target_include_directories(lib_bzmesh PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 
 # Link with reuired libraries
 target_link_libraries(lib_bzmesh PUBLIC libepp)
-target_link_libraries(lib_bzmesh PUBLIC Eigen3::Eigen rapidcsv)
+target_link_libraries(lib_bzmesh PUBLIC Eigen3::Eigen rapidcsv fmt::fmt)
 target_include_directories(lib_bzmesh PUBLIC ${GMSH_INC})
 target_link_libraries(lib_bzmesh PUBLIC ${GMSH_LIB})
 

src/BZ_MESH/bz_states.cpp

@@ -173,4 +173,10 @@ void BZ_States::export_full_eigenstates() const {
     }
 }
 
+
+// void BZ_States::populate_vtx_dielectric_function(const std::vector<double>& energies, double eta_smearing);
+
+
+
+
 }  // namespace bz_mesh

src/BZ_MESH/bz_states.hpp

@@ -39,21 +39,46 @@ class BZ_States : public MeshBZ {
      */
     std::vector<double> m_dielectric_function_real;
 
+    // m_vtx_dielectric_function_real[idx_vtx][idx_energy] is the real part of the dielectric function at the energy m_energies[idx_energy]
+    // and at the vertex m_vertices[idx_vtx].
+    std::vector<std::vector<double>> m_vtx_dielectric_function_real;
+
+    /**
+     * @brief The index of the first q-point this class is responsible for.
+     * This is used to parallelize the computation of the dielectric function where
+     * each instance of this class is responsible for a subset of the q-points of the mesh.
+     *
+     * For a calculation on a single CPU the offset is 0.
+     *
+     */
+    std::size_t m_offset_q_index = 0;
+
+    /**
+     * @brief Number of q-points this class is responsible for.
+     * This is used to parallelize the computation of the dielectric function where
+     * each instance of this class is responsible for a subset of the q-points.
+     *
+     * For a calculation on a single CPU this parameter is equal to the size of the m_vertices.size().
+     *
+     */
+    std::size_t m_nb_kpoints = 0;
+
  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                              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; }
     const std::vector<Vector3D<int>>& get_basis_vectors() const { return m_basisVectors; }
-    void compute_eigenstates(int nb_threads = 1);
-    void compute_shifted_eigenstates(const Vector3D<double>& q_shift, int nb_threads = 1);
+    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;
 
+    void populate_vtx_dielectric_function(const std::vector<double>& energies, double eta_smearing);
 
     void export_full_eigenstates() const;
 };

src/BZ_MESH/mesh_tetra.hpp

@@ -113,6 +113,20 @@ class Tetra {
     double               compute_tetra_iso_surface_energy_band(double energy, std::size_t band_index) const;
     double               compute_tetra_dos_energy_band(double energy, std::size_t band_index) const;
 
+    double  interpolate_scalar_at_position(const std::array<double, 4>& barycentric_coordinates,
+                                           const std::vector<double>&   scalar_field) const;
+    vector3 interpolate_vector_at_position(const std::array<double, 4>& barycentric_coordinates,
+                                           const std::vector<vector3>&  vector_field) const;
+template <typename T>
+      T interpolate_at_position(const std::array<double, 4>& barycentric_coordinates,
+                                 const std::vector<T>&         field) const {
+         T interpolated_value = T::Zero();
+         for (std::size_t idx_vtx = 0; idx_vtx < 4; ++idx_vtx) {
+               interpolated_value += barycentric_coordinates[idx_vtx] * field[idx_vtx];
+         }
+         return interpolated_value;
+      }
+
     static void reset_stat_iso_computing() { ms_case_stats = std::vector<double>(5, 0.0); }
 
     static void print_stat_iso_computing() {

src/EPP/CMakeLists.txt

@@ -32,7 +32,7 @@ add_library(libepp STATIC ${SOURCE_FILES_LIBEPP} ${HEADER_FILES_LIBEPP})
 target_include_directories(libepp PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 
 # Link with reuired libraries
-target_link_libraries(libepp PUBLIC Eigen3::Eigen yaml-cpp rapidcsv)
+target_link_libraries(libepp PUBLIC Eigen3::Eigen yaml-cpp rapidcsv fmt::fmt)
 target_include_directories(libepp PUBLIC ${GMSH_INC})
 target_link_libraries(libepp PUBLIC ${GMSH_LIB})
 

src/EPP/DielectricFunction.cpp

@@ -9,6 +9,10 @@
  *
  */
 
+#include <fmt/core.h>
+#include <fmt/format.h>
+#include <fmt/ostream.h>
+
 #include "DielectricFunction.hpp"
 
 #include <algorithm>
@@ -190,6 +194,10 @@ DielectricFunction DielectricFunction::merge_results(DielectricFunction
     for (std::size_t index_q = 0; index_q < total_dielectric_function.size(); ++index_q) {
         Vector3D<double> q_vect    = RootDielectricFunction.m_qpoints[index_q];
         double           q_squared = pow(q_vect.Length(), 2);
+        // TO IMPROVE: This is a hack to avoid division by zero.
+        if (q_squared == 0.0) {
+            q_squared = 1e-14;
+        }
         for (std::size_t index_energy = 0; index_energy < total_dielectric_function[index_q].size(); ++index_energy) {
             total_dielectric_function[index_q][index_energy] *= renormalization;
             total_dielectric_function[index_q][index_energy] =
@@ -245,8 +253,9 @@ void DielectricFunction::apply_kramers_kronig() {
 void DielectricFunction::export_dielectric_function_at_q(const std::string& filename, std::size_t idx_q, bool name_auto) const {
     std::string outname;
     if (name_auto) {
-        outname = m_export_prefix + '_' + 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 = 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);
     } else {
         outname = filename;
     }