From 1e0a7ad673909c5371aa3c85dd8cdecc9a61159b Mon Sep 17 00:00:00 2001 From: RemiHelleboid Date: Fri, 1 Dec 2023 01:05:54 +0100 Subject: [PATCH] upp --- python/sphere_mesh.py | 4 +- src/BZ_MESH/bz_mesh.cpp | 35 ++++++++++++ src/BZ_MESH/bz_mesh.hpp | 7 ++- src/BZ_MESH/bz_states.hpp | 8 ++- src/BZ_MESH/impact_ionization.cpp | 89 ++++++++++++++++--------------- src/BZ_MESH/impact_ionization.hpp | 21 ++++---- src/BZ_MESH/vector.hpp | 6 +++ 7 files changed, 113 insertions(+), 57 deletions(-) diff --git a/python/sphere_mesh.py b/python/sphere_mesh.py index 8334e89..2c866b9 100644 --- a/python/sphere_mesh.py +++ b/python/sphere_mesh.py @@ -9,8 +9,8 @@ # 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) +gmsh.option.setNumber("Mesh.MeshSizeMin", 0.05) +gmsh.option.setNumber("Mesh.MeshSizeMax", 0.10) Rt = 1.0 diff --git a/src/BZ_MESH/bz_mesh.cpp b/src/BZ_MESH/bz_mesh.cpp index 3e1b841..1d312e1 100644 --- a/src/BZ_MESH/bz_mesh.cpp +++ b/src/BZ_MESH/bz_mesh.cpp @@ -96,6 +96,41 @@ void MeshBZ::read_mesh_geometry_from_msh_file(const std::string& filename) { m_total_volume = compute_mesh_volume(); } +/** + * @brief Get the nearest k index object. + * Brute force search of the nearest k-point index. :( + * + * @param k + * @return std::size_t + */ +std::size_t MeshBZ::get_nearest_k_index(const Vector3D& k) const { + vector3 K(k.X, k.Y, k.Z); + std::size_t index_nearest_k = 0; + double min_distance = std::numeric_limits::max(); + for (std::size_t index_k = 0; index_k < m_list_vertices.size(); ++index_k) { + double distance = (K - m_list_vertices[index_k].get_position()).norm(); + if (distance < min_distance) { + min_distance = distance; + index_nearest_k = index_k; + } + } + return index_nearest_k; +} +std::size_t MeshBZ::get_nearest_k_index(const vector3& k) const { + std::size_t index_nearest_k = 0; + double min_distance = std::numeric_limits::max(); + for (std::size_t index_k = 0; index_k < m_list_vertices.size(); ++index_k) { + double distance = (k - m_list_vertices[index_k].get_position()).norm(); + if (distance < min_distance) { + min_distance = distance; + index_nearest_k = index_k; + } + } + return index_nearest_k; +} + + + /** * @brief Read the energy values for each band at every k-points (vertices) of the mesh. * diff --git a/src/BZ_MESH/bz_mesh.hpp b/src/BZ_MESH/bz_mesh.hpp index 3ae7e67..47abcdc 100644 --- a/src/BZ_MESH/bz_mesh.hpp +++ b/src/BZ_MESH/bz_mesh.hpp @@ -81,8 +81,13 @@ class MeshBZ { std::size_t get_number_vertices() const { return m_list_vertices.size(); } std::size_t get_number_elements() const { return m_list_tetrahedra.size(); } - std::size_t get_number_bands() const { return m_min_band.size(); } + double get_volume() const { return m_total_volume; } + + vector3 get_k_at_index(std::size_t index) const { return m_list_vertices[index].get_position(); } + std::size_t get_nearest_k_index(const Vector3D& k) const; + std::size_t get_nearest_k_index(const vector3& k) const; + std::size_t get_number_bands() const { return m_min_band.size(); } std::pair get_min_max_energy_at_band(const int& band_index) const { return std::make_pair(m_min_band[band_index], m_max_band[band_index]); } diff --git a/src/BZ_MESH/bz_states.hpp b/src/BZ_MESH/bz_states.hpp index c1dd05a..1302868 100644 --- a/src/BZ_MESH/bz_states.hpp +++ b/src/BZ_MESH/bz_states.hpp @@ -19,13 +19,14 @@ namespace bz_mesh { class BZ_States : public MeshBZ { - private: + protected: int m_nb_bands = 0; std::vector> m_basisVectors; std::vector m_eigenvalues_k; std::vector m_eigenvalues_k_plus_q; + std::vector m_eigenvectors_k; std::vector m_eigenvectors_k_plus_q; @@ -80,6 +81,11 @@ class BZ_States : public MeshBZ { void populate_vtx_dielectric_function(const std::vector& energies, double eta_smearing); + std::complex get_dielectric_function(const vector3& q, double energy) const { + // TODO: implement this function + return 1.0; + } + void export_full_eigenstates() const; }; diff --git a/src/BZ_MESH/impact_ionization.cpp b/src/BZ_MESH/impact_ionization.cpp index 996729a..bbd8ec3 100644 --- a/src/BZ_MESH/impact_ionization.cpp +++ b/src/BZ_MESH/impact_ionization.cpp @@ -15,6 +15,7 @@ #include #include +#include "Constants.hpp" #include "Hamiltonian.h" #include "bz_mesh.hpp" #include "bz_states.hpp" @@ -22,60 +23,64 @@ namespace bz_mesh { -void ImpactIonization::compute_impact_ionization_rate() { - this->compute_eigenstates(); - auto vertices = this->get_list_vertices(); +double ImpactIonization::compute_impact_ionization_rate(int idx_n1, const Vector3D& k1) { + double rate = 0.0; + int nb_conduction_bands = 4; + int nb_valence_bands = 4; + // Pseudo code + return rate; } -double ImpactIonization::compute_direct_impact_ionization_matrix_element(int idx_n1, - int idx_n1_prime, - int idx_n2, - int idx_n2_prime, - const Vector3D& k1, - const Vector3D& k1_prime, - const Vector3D& k2, - const Vector3D& k2_prime) const { - // Momentum conservation - // Vector3D q = k1 + k1_prime - k2 - k2_prime; - // if (q.norm() > 1e-12) { - // std::cout << "Momentum conservation is not satisfied" << std::endl; - // return 0.0; - // } +double ImpactIonization::compute_direct_impact_ionization_matrix_element(int idx_n1, + int idx_n1_prime, + int idx_n2, + int idx_n2_prime, + int idx_k1, + int idx_k1_prime, + int idx_k2_prime) const { + vector3 k1 = m_list_vertices[idx_k1].get_position(); + vector3 k1_prime = m_list_vertices[idx_k1_prime].get_position(); + vector3 k2_prime = m_list_vertices[idx_k2_prime].get_position(); + vector3 k2 = k1_prime + k2_prime - k1; - const auto& basis_vectors = get_basis_vectors(); - std::size_t nb_basis_vectors = basis_vectors.size(); - double matrix_element = 0.0; + std::size_t index_k2 = get_nearest_k_index(k2); + std::complex Ma = 0.0; - Eigen::VectorXd m_eigenvalues_k1; - Eigen::VectorXd m_eigenvalues_k1_prime; - Eigen::VectorXd m_eigenvalues_k2; - Eigen::VectorXd m_eigenvalues_k2_prime; - Eigen::MatrixXcd m_eigenvectors_k1; - Eigen::MatrixXcd m_eigenvectors_k1_prime; - Eigen::MatrixXcd m_eigenvectors_k2; - Eigen::MatrixXcd m_eigenvectors_k2_prime; - // Get the eigenvalues and eigenvectors for the four k-points (TODO) + double e_charge = EmpiricalPseudopotential::Constants::q; + double eps_zero = EmpiricalPseudopotential::Constants::eps_zero; + double volume = compute_mesh_volume(); + double fourier_factor = 2.0 * M_PI / volume; - double w_a = (m_eigenvalues_k1_prime[idx_n1_prime] - m_eigenvalues_k1[idx_n1]) / EmpiricalPseudopotential::Constants::h_bar; + auto basis_vectors = get_basis_vectors(); + std::size_t size_basis_vectors = basis_vectors.size(); + for (std::size_t idx_g1 = 0; idx_g1 < size_basis_vectors; ++idx_g1) { + for (std::size_t idx_g2 = 0; idx_g2 < size_basis_vectors; ++idx_g2) { + for (std::size_t idx_g1_prime = 0; idx_g1_prime < size_basis_vectors; ++idx_g1_prime) { + for (std::size_t idx_g2_prime = 0; idx_g2_prime < size_basis_vectors; ++idx_g2_prime) { + vector3 g1 = {basis_vectors[idx_g1].X * fourier_factor, + basis_vectors[idx_g1].Y * fourier_factor, + basis_vectors[idx_g1].Z * fourier_factor}; + vector3 g1_prime = {basis_vectors[idx_g1_prime].X * fourier_factor, + basis_vectors[idx_g1_prime].Y * fourier_factor, + basis_vectors[idx_g1_prime].Z * fourier_factor}; + vector3 q_a = k1_prime + g1_prime - k1 - g1;.0000000000000000000000000000000000000000000000000000000000000000000000 + auto state_k1 = m_eigenvectors_k[idx_k1](idx_n1, idx_g1); + auto state_k2 = m_eigenvectors_k[index_k2](idx_n2, idx_g2); + auto state_k1_prime = m_eigenvectors_k[idx_k1_prime](idx_n1_prime, idx_g1_prime); + auto state_k2_prime = m_eigenvectors_k[idx_k2_prime](idx_n2_prime, idx_g2_prime); - // Quadruple sum over the basis vectors - for (std::size_t idx_g1 = 0; idx_g1 < nb_basis_vectors; ++idx_g1) { - for (std::size_t idx_g2 = 0; idx_g2 < nb_basis_vectors; ++idx_g2) { - for (std::size_t idx_g3 = 0; idx_g3 < nb_basis_vectors; ++idx_g3) { - for (std::size_t idx_g4 = 0; idx_g4 < nb_basis_vectors; ++idx_g4) { - auto g1 = basis_vectors[idx_g1]; - auto g2 = basis_vectors[idx_g2]; - auto g3 = basis_vectors[idx_g3]; - auto g4 = basis_vectors[idx_g4]; + double omega_a = m_eigenvalues_k[idx_k1](idx_n1) - m_eigenvalues_k[idx_k1_prime](idx_n1_prime); + double omega_b = m_eigenvalues_k[idx_k2_prime](idx_n2_prime) - m_eigenvalues_k[idx_k1](idx_n1); + std::complex micro_matrix_element_a = std::conj(state_k1_prime) * std::conj(state_k2_prime) * state_k1 * state_k2; + std::complex dielectric_func = get_dielectric_function(q_a, omega_a); + std::complex pre_factor = e_charge * e_charge / (eps_zero * dielectric_func * q_a.norm() * q_a.norm() * volume); - auto A1_G1 = m_eigenvectors_k1(idx_n1, idx_g1); + Ma += pre_factor * micro_matrix_element; } } } } - - return matrix_element; } } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/impact_ionization.hpp b/src/BZ_MESH/impact_ionization.hpp index d804456..7d6f662 100644 --- a/src/BZ_MESH/impact_ionization.hpp +++ b/src/BZ_MESH/impact_ionization.hpp @@ -21,18 +21,17 @@ namespace bz_mesh { class ImpactIonization : private BZ_States { private: + double m_impact_ionization_rate = 0.0; public: - double compute_direct_impact_ionization_matrix_element(int idx_n1, - int idx_n1_prime, - int idx_n2, - int idx_n2_prime, - const Vector3D& k1, - const Vector3D& k1_prime, - const Vector3D& k2, - const Vector3D& k2_prime) const; - - void compute_impact_ionization_rate(); + double compute_direct_impact_ionization_matrix_element(int idx_n1, + int idx_n1_prime, + int idx_n2, + int idx_n2_prime, + int idx_k1, + int idx_k1_prime, + int idx_k2_prime) const; + + double compute_impact_ionization_rate(int idx_n1, const Vector3D& k1); }; - } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/vector.hpp b/src/BZ_MESH/vector.hpp index d710fb4..09f079d 100644 --- a/src/BZ_MESH/vector.hpp +++ b/src/BZ_MESH/vector.hpp @@ -15,6 +15,9 @@ #include #include #include + +#include "Vector3D.h" + namespace bz_mesh { enum class permutaion_type { XY, XZ, YZ, XYZ, YZX, ZXY }; @@ -29,6 +32,9 @@ class vector3 { vector3() : m_x(0u), m_y(0u), m_z(0u) {} vector3(double x, double y) : m_x(x), m_y(y), m_z(0u) {} vector3(double x, double y, double z) : m_x(x), m_y(y), m_z(z) {} + vector3(const vector3 &vector) : m_x(vector.m_x), m_y(vector.m_y), m_z(vector.m_z) {} + vector3(const Vector3D &vector) : m_x(vector.X), m_y(vector.Y), m_z(vector.Z) {} + vector3(const Vector3D &vector) : m_x(vector.X), m_y(vector.Y), m_z(vector.Z) {} double x() const { return m_x; } double y() const { return m_y; }