Update epsilon

python/PopulateMeshWithDielectric.py

@@ -25,22 +25,36 @@ def extract_dielectric_file(filename):
 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]
+print(f"list_idx_dielectric = {list_idx_dielectric}")
+arg_sorted = np.argsort(list_idx_dielectric)
+list_dielectric_files = [list_dielectric_files[i] for i in arg_sorted]
+
 
 NbFiles = len(list_dielectric_files)
 
 # Get list of nodes
 nodes = gmsh.model.mesh.getNodes()
 TagNodes = nodes[0]
 NbNodes = TagNodes.shape[0]
+Coordinates = nodes[1].reshape((NbNodes, 3))
+
+# Debug
+for k in range(10):
+    print(f"Coordinates[{k}] = {Coordinates[k]}")
+    print(f"Filename = {list_dielectric_files[k]}")
+    print(f"\n")
+
+
+
 
 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}")
+# energy, eps_r, eps_i = extract_dielectric_file(list_dielectric_files[0])
+# print(f"energy = {energy}")
 LIST_ENERGY = []
 LIST_EPS_R = []
 LIST_EPS_I = []
@@ -65,7 +79,7 @@ def extract_dielectric_file(filename):
 
 FileName = "EightSphere_3.0_with_dielectric.msh"
 
-for ix_e in range(NbEnergies//10):
+for ix_e in range(NbEnergies):
     print(f"idx_e = {ix_e}")
     e = energies[ix_e]
     eps_r = LIST_EPS_R[:, ix_e]

python/sphere_mesh.py

@@ -9,27 +9,31 @@
 # from http://en.wikipedia.org/wiki/Constructive_solid_geometry
 
 gmsh.option.setNumber("Mesh.Algorithm", 6)
-gmsh.option.setNumber("Mesh.MeshSizeMin", 0.05)
-gmsh.option.setNumber("Mesh.MeshSizeMax", 0.10)
+
+minSize = 0.2
+maxSize = 0.25
+gmsh.option.setNumber("Mesh.MeshSizeMin", minSize)
+gmsh.option.setNumber("Mesh.MeshSizeMax", maxSize)
 
 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.addBox(0, 0, 0, 2 * R, 2 * R, 2 * 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, 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.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.occ.removeAllDuplicates()
 
 gmsh.model.mesh.generate(3)
 
@@ -39,7 +43,7 @@
 nodes = gmsh.model.mesh.getNodes()
 Coordinates = nodes[1].reshape((nodes[0].shape[0], 3))
 
-filename = f"nodes_eight_sphere_{Rt}.csv"
+filename = f"nodes_eight_sphere_{R}.csv"
 np.savetxt(filename, Coordinates, delimiter=" ")
 
 gmsh.write(f"EightSphere_{R}.msh")

src/BZ_MESH/dielectric_mesh.hpp

@@ -0,0 +1,36 @@
+/**
+ * @file dielectric_mesh.hpp
+ * @author your name (you@domain.com)
+ * @brief
+ * @version 0.1
+ * @date 2024-05-17
+ *
+ * @copyright Copyright (c) 2024
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include "Material.h"
+#include "bz_mesh.hpp"
+
+namespace bz_mesh {
+
+class DielectricMesh : public MeshBZ {
+ protected:
+    int m_nb_bands = 0;
+
+ public:
+    DielectricMesh(const EmpiricalPseudopotential::Material& material) : MeshBZ(material) {}
+    void read_dielectric_file(const std::string& filename);
+};
+
+}  // namespace bz_mesh

src/EPP/DielectricFunction.cpp

@@ -9,12 +9,12 @@
  *
  */
 
+#include "DielectricFunction.hpp"
+
 #include <fmt/core.h>
 #include <fmt/format.h>
 #include <fmt/ostream.h>
 
-#include "DielectricFunction.hpp"
-
 #include <algorithm>
 #include <chrono>
 #include <cmath>
@@ -47,7 +47,7 @@ DielectricFunction::DielectricFunction(const Material& material, const std::vect
 void DielectricFunction::generate_k_points_random(std::size_t nb_points) {
     std::random_device               rd;
     std::mt19937                     gen(rd());
-    std::uniform_real_distribution<> dis(-1 - .0, 1.0);
+    std::uniform_real_distribution<> dis(-1, 1.0);
     while (m_kpoints.size() < nb_points) {
         Vector3D<double> k(dis(gen), dis(gen), dis(gen));
         if (is_in_first_BZ(k)) {
@@ -93,6 +93,9 @@ void DielectricFunction::compute_dielectric_function(double eta_smearing, int mp
     Hamiltonian hamiltonian_k_plus_q(m_material, m_basisVectors);
     for (std::size_t index_q = 0; index_q < m_qpoints.size(); ++index_q) {
         Vector3D<double>              q_vect = m_qpoints[index_q];
+        if (q_vect.Length() <= 1e-13) {
+            q_vect = Vector3D<double>(1e-13, 1e-13, 1e-13);
+        }
         std::vector<Vector3D<double>> k_plus_q_vects(m_kpoints.size());
         std::transform(m_kpoints.begin(), m_kpoints.end(), k_plus_q_vects.begin(), [&q_vect](const Vector3D<double>& k) {
             return k + q_vect;
@@ -146,12 +149,11 @@ void DielectricFunction::compute_dielectric_function(double eta_smearing, int mp
         for (std::size_t index_energy = 0; index_energy < m_energies.size(); ++index_energy) {
             list_epsilon[index_energy] = list_total_sum[index_energy];
         }
-        auto end = std::chrono::high_resolution_clock::now();
-        auto elapsed =
-            std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / 1000.0;
+        auto end     = std::chrono::high_resolution_clock::now();
+        auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / 1000.0;
         if (mpi_rank == 0) {
-        std::cout << "Computed dielectric function for q = " << m_qpoints[index_q] << " -> " << index_q + 1 << "/" << m_qpoints.size() << " in "
-                  << elapsed << " s" << std::endl;
+            std::cout << "Computed dielectric function for q = " << m_qpoints[index_q] << " -> " << index_q + 1 << "/" << m_qpoints.size()
+                      << " in " << elapsed << " s" << std::endl;
         }
         start = std::chrono::high_resolution_clock::now();
         m_dielectric_function_real.push_back(list_epsilon);
@@ -253,9 +255,15 @@ 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(idx_q) + '_' + std::to_string(m_qpoints[idx_q].X) + "_" + std::to_string(m_qpoints[idx_q].Y) + "_" +
+        // 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}_{:.6f}_{:.6f}_{:.6f}.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;
     }