From 23d873f20a586206a23d7e5b248122687238a641 Mon Sep 17 00:00:00 2001 From: RemiHelleboid Date: Thu, 30 Nov 2023 16:39:43 +0100 Subject: [PATCH] Up --- CMakeLists.txt | 8 +++ apps/epsilon.cpp | 3 + python/PopulateMeshWithDielectric.py | 104 +++++++++++++++++++++++++++ python/sphere_mesh.py | 50 +++++++++++++ src/BZ_MESH/CMakeLists.txt | 2 +- src/BZ_MESH/bz_states.cpp | 6 ++ src/BZ_MESH/bz_states.hpp | 33 +++++++-- src/BZ_MESH/mesh_tetra.hpp | 14 ++++ src/EPP/CMakeLists.txt | 2 +- src/EPP/DielectricFunction.cpp | 13 +++- 10 files changed, 227 insertions(+), 8 deletions(-) create mode 100644 python/PopulateMeshWithDielectric.py create mode 100644 python/sphere_mesh.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 14d805f..0d69f4c 100644 --- a/CMakeLists.txt +++ b/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.") diff --git a/apps/epsilon.cpp b/apps/epsilon.cpp index 2f7d5da..67e6607 100644 --- a/apps/epsilon.cpp +++ b/apps/epsilon.cpp @@ -46,6 +46,9 @@ std::vector> read_qpoint_dat_file(const std::string& filename) std::vector> kpoints; std::string line; while (std::getline(file, line)) { + if (line.empty()) { + continue; + } std::stringstream ss(line); std::string token; std::vector kpoint; diff --git a/python/PopulateMeshWithDielectric.py b/python/PopulateMeshWithDielectric.py new file mode 100644 index 0000000..8b38d10 --- /dev/null +++ b/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) diff --git a/python/sphere_mesh.py b/python/sphere_mesh.py new file mode 100644 index 0000000..8334e89 --- /dev/null +++ b/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() diff --git a/src/BZ_MESH/CMakeLists.txt b/src/BZ_MESH/CMakeLists.txt index c62e1a7..e2e5e5d 100644 --- a/src/BZ_MESH/CMakeLists.txt +++ b/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}) diff --git a/src/BZ_MESH/bz_states.cpp b/src/BZ_MESH/bz_states.cpp index a59b130..9a5bf7b 100644 --- a/src/BZ_MESH/bz_states.cpp +++ b/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& energies, double eta_smearing); + + + + } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/bz_states.hpp b/src/BZ_MESH/bz_states.hpp index 5551903..c1dd05a 100644 --- a/src/BZ_MESH/bz_states.hpp +++ b/src/BZ_MESH/bz_states.hpp @@ -39,14 +39,38 @@ class BZ_States : public MeshBZ { */ std::vector 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> 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>& 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>& basis_vectors) { m_basisVectors = basis_vectors; } const std::vector>& get_basis_vectors() const { return m_basisVectors; } - void compute_eigenstates(int nb_threads = 1); - void compute_shifted_eigenstates(const Vector3D& q_shift, int nb_threads = 1); + void compute_eigenstates(int nb_threads = 1); + void compute_shifted_eigenstates(const Vector3D& q_shift, int nb_threads = 1); const std::vector& get_energies() const { return m_list_energies; } void set_energies(const std::vector& energies) { m_list_energies = energies; } @@ -54,6 +78,7 @@ class BZ_States : public MeshBZ { void compute_dielectric_function(const std::vector& 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& energies, double eta_smearing); void export_full_eigenstates() const; }; diff --git a/src/BZ_MESH/mesh_tetra.hpp b/src/BZ_MESH/mesh_tetra.hpp index 8f5b9aa..38b13dd 100644 --- a/src/BZ_MESH/mesh_tetra.hpp +++ b/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& barycentric_coordinates, + const std::vector& scalar_field) const; + vector3 interpolate_vector_at_position(const std::array& barycentric_coordinates, + const std::vector& vector_field) const; +template + T interpolate_at_position(const std::array& barycentric_coordinates, + const std::vector& 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(5, 0.0); } static void print_stat_iso_computing() { diff --git a/src/EPP/CMakeLists.txt b/src/EPP/CMakeLists.txt index ffcccfe..207c0ae 100644 --- a/src/EPP/CMakeLists.txt +++ b/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}) diff --git a/src/EPP/DielectricFunction.cpp b/src/EPP/DielectricFunction.cpp index 8755a42..a66af4f 100644 --- a/src/EPP/DielectricFunction.cpp +++ b/src/EPP/DielectricFunction.cpp @@ -9,6 +9,10 @@ * */ +#include +#include +#include + #include "DielectricFunction.hpp" #include @@ -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 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; }