Skip to content

Commit

Permalink
upp
Browse files Browse the repository at this point in the history
  • Loading branch information
@RemiHelleboid
RemiHelleboid committed Dec 1, 2023
1 parent 23d873f commit 1e0a7ad
Show file tree
Hide file tree
Showing 7 changed files with 113 additions and 57 deletions.
4 changes: 2 additions & 2 deletions python/sphere_mesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
35 changes: 35 additions & 0 deletions src/BZ_MESH/bz_mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>& k) const {
vector3 K(k.X, k.Y, k.Z);
std::size_t index_nearest_k = 0;
double min_distance = std::numeric_limits<double>::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<double>::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.
*
Expand Down
7 changes: 6 additions & 1 deletion src/BZ_MESH/bz_mesh.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>& 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<double, double> 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]);
}
Expand Down
8 changes: 7 additions & 1 deletion src/BZ_MESH/bz_states.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,14 @@
namespace bz_mesh {

class BZ_States : public MeshBZ {
private:
protected:
int m_nb_bands = 0;

std::vector<Vector3D<int>> m_basisVectors;

std::vector<Eigen::VectorXd> m_eigenvalues_k;
std::vector<Eigen::VectorXd> m_eigenvalues_k_plus_q;

std::vector<Eigen::MatrixXcd> m_eigenvectors_k;
std::vector<Eigen::MatrixXcd> m_eigenvectors_k_plus_q;

Expand Down Expand Up @@ -80,6 +81,11 @@ class BZ_States : public MeshBZ {

void populate_vtx_dielectric_function(const std::vector<double>& energies, double eta_smearing);

std::complex<double> get_dielectric_function(const vector3& q, double energy) const {
// TODO: implement this function
return 1.0;
}

void export_full_eigenstates() const;
};

Expand Down
89 changes: 47 additions & 42 deletions src/BZ_MESH/impact_ionization.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,67 +15,72 @@
#include <filesystem>
#include <fstream>

#include "Constants.hpp"
#include "Hamiltonian.h"
#include "bz_mesh.hpp"
#include "bz_states.hpp"
#include "omp.h"

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<double>& 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<double>& k1,
const Vector3D<double>& k1_prime,
const Vector3D<double>& k2,
const Vector3D<double>& k2_prime) const {
// Momentum conservation
// Vector3D<double> 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<double> 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<double> micro_matrix_element_a = std::conj(state_k1_prime) * std::conj(state_k2_prime) * state_k1 * state_k2;
std::complex<double> dielectric_func = get_dielectric_function(q_a, omega_a);
std::complex<double> 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
21 changes: 10 additions & 11 deletions src/BZ_MESH/impact_ionization.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>& k1,
const Vector3D<double>& k1_prime,
const Vector3D<double>& k2,
const Vector3D<double>& 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<double>& k1);
};


} // namespace bz_mesh
6 changes: 6 additions & 0 deletions src/BZ_MESH/vector.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,9 @@
#include <iostream>
#include <optional>
#include <vector>

#include "Vector3D.h"

namespace bz_mesh {

enum class permutaion_type { XY, XZ, YZ, XYZ, YZX, ZXY };
Expand All @@ -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<double> &vector) : m_x(vector.X), m_y(vector.Y), m_z(vector.Z) {}
vector3(const Vector3D<int> &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; }
Expand Down

0 comments on commit 1e0a7ad

Please sign in to comment.