diff --git a/python/PopulateMeshWithDielectric.py b/python/PopulateMeshWithDielectric.py index 8b38d10..1f2254c 100644 --- a/python/PopulateMeshWithDielectric.py +++ b/python/PopulateMeshWithDielectric.py @@ -25,6 +25,10 @@ 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) @@ -32,6 +36,16 @@ def extract_dielectric_file(filename): 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}") @@ -39,8 +53,8 @@ def extract_dielectric_file(filename): 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] diff --git a/python/sphere_mesh.py b/python/sphere_mesh.py index 2c866b9..5858847 100644 --- a/python/sphere_mesh.py +++ b/python/sphere_mesh.py @@ -9,8 +9,11 @@ # 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 @@ -18,18 +21,19 @@ 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") diff --git a/src/BZ_MESH/dielectric_mesh.hpp b/src/BZ_MESH/dielectric_mesh.hpp new file mode 100644 index 0000000..f69f8ab --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include + +#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 \ No newline at end of file diff --git a/src/EPP/DielectricFunction.cpp b/src/EPP/DielectricFunction.cpp index 909b8d8..c96239a 100644 --- a/src/EPP/DielectricFunction.cpp +++ b/src/EPP/DielectricFunction.cpp @@ -9,12 +9,12 @@ * */ +#include "DielectricFunction.hpp" + #include #include #include -#include "DielectricFunction.hpp" - #include #include #include @@ -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 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 q_vect = m_qpoints[index_q]; + if (q_vect.Length() <= 1e-13) { + q_vect = Vector3D(1e-13, 1e-13, 1e-13); + } std::vector> k_plus_q_vects(m_kpoints.size()); std::transform(m_kpoints.begin(), m_kpoints.end(), k_plus_q_vects.begin(), [&q_vect](const Vector3D& 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(end - start).count() / 1000.0; + auto end = std::chrono::high_resolution_clock::now(); + auto elapsed = std::chrono::duration_cast(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; }