Skip to content

Commit

Permalink
Up
Browse files Browse the repository at this point in the history
  • Loading branch information
@RemiHelleboid
RemiHelleboid committed May 22, 2024
1 parent 1874321 commit 88efb75
Show file tree
Hide file tree
Showing 6 changed files with 112 additions and 115 deletions.
9 changes: 5 additions & 4 deletions apps/impact_io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ int main(int argc, char *argv[]) {
TCLAP::ValueArg<std::string> arg_dielectric_file("y", "dielectric", "File (.msh) with the dielectric function.", true, "", "string");
TCLAP::ValueArg<std::string> arg_material("m", "material", "Symbol of the material to use (Si, Ge, GaAs, ...)", true, "Si", "string");
TCLAP::ValueArg<int> arg_nb_energies("e", "nenergy", "Number of energies to compute", false, 250, "int");
TCLAP::ValueArg<int> arg_nb_bands("b", "nbands", "Number of bands to consider", false, 12, "int");
TCLAP::ValueArg<int> arg_nb_bands("b", "nbands", "Number of bands to consider", false, 16, "int");
TCLAP::ValueArg<int> arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int");
TCLAP::ValueArg<double> arg_radius_BZ("R",
"radiusBZ",
Expand Down Expand Up @@ -67,16 +67,17 @@ int main(int argc, char *argv[]) {
my_options.nrThreads = arg_nb_threads.getValue();
my_options.print_options();
int nb_bands_to_use = arg_nb_bands.getValue();
int nb_threads = arg_nb_threads.getValue();
auto start = std::chrono::high_resolution_clock::now();

EmpiricalPseudopotential::Material current_material = materials.materials.at(arg_material.getValue());

bz_mesh::ImpactIonization my_impact_ionization(current_material, arg_mesh_file.getValue());
my_impact_ionization.read_dielectric_file(arg_dielectric_file.getValue());

my_impact_ionization.read_dielectric_file(arg_dielectric_file.getValue());
my_impact_ionization.interp_test_dielectric_function("test_dielectric_function.csv");

my_impact_ionization.set_max_radius_G0_BZ(arg_radius_BZ.getValue());
my_impact_ionization.compute_eigenstates();
my_impact_ionization.compute_eigenstates(nb_threads);

auto end = std::chrono::high_resolution_clock::now();

Expand Down
15 changes: 6 additions & 9 deletions src/BZ_MESH/bz_states.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,11 @@ void BZ_States::compute_eigenstates(int nb_threads) {
for (int i = 0; i < nb_threads; i++) {
hamiltonian_per_thread.push_back(EmpiricalPseudopotential::Hamiltonian(m_material, m_basisVectors));
}
#pragma omp parallel for schedule(dynamic) num_threads(nb_threads)
for (std::size_t idx_k = 0; idx_k < m_list_vertices.size(); ++idx_k) {
std::cout << "\rComputing eigenstates for k = " << idx_k << "/" << m_list_vertices.size() << std::flush;
if (omp_get_thread_num() == 0) {
std::cout << "\rComputing eigenstates for k = " << idx_k << "/" << m_list_vertices.size() << std::flush;
}
auto k_point = Vector3D<double>(m_list_vertices[idx_k].get_position().x(),
m_list_vertices[idx_k].get_position().y(),
m_list_vertices[idx_k].get_position().z());
Expand Down Expand Up @@ -93,7 +96,7 @@ void BZ_States::compute_dielectric_function(const std::vector<double>& list_ener
constexpr double one_fourth = 1.0 / 4.0;

std::vector<double> dielectric_function_real_at_energies(list_energies.size(), 0.0);
double total_volume = 0.0;
double total_volume = 0.0;
for (std::size_t idx_tetra = 0; idx_tetra < nb_tetra; ++idx_tetra) {
std::cout << "\rComputing dielectric function for tetrahedron " << idx_tetra << "/" << nb_tetra << std::flush;
std::array<std::size_t, 4> list_idx_vertices = m_list_tetrahedra[idx_tetra].get_list_indices_vertices();
Expand Down Expand Up @@ -133,14 +136,12 @@ void BZ_States::compute_dielectric_function(const std::vector<double>& list_ener
std::cout << "\n";
std::cout << "Total volume: " << total_volume << std::endl;

double q_squared = m_q_shift.Length() * m_q_shift.Length();
double q_squared = m_q_shift.Length() * m_q_shift.Length();
double pre_factor = 2.0 * M_PI / q_squared;
for (std::size_t index_energy = 0; index_energy < list_energies.size(); ++index_energy) {
m_dielectric_function_real[index_energy] = 1.0 + pre_factor * dielectric_function_real_at_energies[index_energy] / total_volume;
}
std::cout << "EPS[0] = " << m_dielectric_function_real[0] << std::endl;


}

// Export the dielectric function to a file in the format (energy, dielectric function) (csv format).
Expand Down Expand Up @@ -175,10 +176,6 @@ void BZ_States::export_full_eigenstates() const {
}
}


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




} // namespace bz_mesh
43 changes: 30 additions & 13 deletions src/BZ_MESH/dielectric_mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,14 +82,17 @@ void DielectricMesh::read_dielectric_file(const std::string& filename) {
}
gmsh::finalize();
m_energies = energies;
m_dielectric_function.resize(real_dielectric.size());
for (std::size_t idx_node = 0; idx_node < real_dielectric.size(); ++idx_node) {
m_dielectric_function.resize(m_list_vertices.size());
for (std::size_t idx_node = 0; idx_node < m_list_vertices.size(); ++idx_node) {
m_dielectric_function[idx_node].resize(energies.size());
for (std::size_t idx_energy = 0; idx_energy < energies.size(); ++idx_energy) {
m_dielectric_function[idx_node][idx_energy] =
std::complex<double>(real_dielectric[idx_node][idx_energy], imag_dielectric[idx_node][idx_energy]);
std::complex<double>(real_dielectric[idx_energy][idx_node], imag_dielectric[idx_energy][idx_node]);
}
}
std::cout << "Size of the dielectric function: " << m_dielectric_function.size() << " x " << m_dielectric_function[0].size()
<< std::endl;
std::cout << "Dielectric function read." << std::endl;
}

std::pair<std::size_t, double> DielectricMesh::find_closest_energy(double energy) const {
Expand All @@ -115,21 +118,35 @@ std::pair<std::size_t, double> DielectricMesh::find_closest_energy(double energy
}

complex_d DielectricMesh::interpolate_dielectric_function(const vector3& k, double energy) const {
auto k_positive{vector3{std::abs(k.x()), std::abs(k.y()), std::abs(k.z())}};
auto k_positive{vector3{std::abs(k.x()), std::abs(k.y()), std::abs(k.z())}};
Tetra* p_tetra = find_tetra_at_location(k_positive);
if (p_tetra == nullptr) {
std::cerr << "Tetra not found at location " << k_positive << std::endl;
return 0.0;
}
const std::array<double, 4> barycentric_coordinates = p_tetra->compute_barycentric_coordinates(k);
std::array<double, 4> barycentric_coordinates = p_tetra->compute_barycentric_coordinates(k_positive);

const std::size_t idx_node = p_tetra->get_index();
const std::pair<std::size_t, double> closest_energy = find_closest_energy(energy);
const std::size_t idx_energy = closest_energy.first;
const double t = closest_energy.second;
const complex_d dielectric_function_0 = m_dielectric_function[idx_node][idx_energy];
const complex_d dielectric_function_1 = m_dielectric_function[idx_node][idx_energy + 1];
complex_d dielectric_function = dielectric_function_0 * (1.0 - t) + dielectric_function_1 * t;
return dielectric_function;
std::array<std::size_t, 4> list_indices_vertices = p_tetra->get_list_indices_vertices();
std::pair<std::size_t, double> closest_energy = find_closest_energy(energy);
std::size_t idx_energy = closest_energy.first;
double t = closest_energy.second;
std::vector<std::complex<double>> dielectric_function_low(4);
std::vector<std::complex<double>> dielectric_function_high(4);
for (std::size_t idx_vertex = 0; idx_vertex < 4; ++idx_vertex) {
std::cout << "Vertex: " << list_indices_vertices[idx_vertex] << std::endl;
dielectric_function_low[idx_vertex] = m_dielectric_function[idx_energy][list_indices_vertices[idx_vertex]];
dielectric_function_high[idx_vertex] = m_dielectric_function[idx_energy + 1][list_indices_vertices[idx_vertex]];
}
std::complex<double> dielectric_function_interpolated_low =
p_tetra->interpolate_at_position(barycentric_coordinates, dielectric_function_low);
std::complex<double> dielectric_function_interpolated_high =
p_tetra->interpolate_at_position(barycentric_coordinates, dielectric_function_high);

std::cout << "Idx energy: " << idx_energy << std::endl;
std::cout << "Energy: " << energy << std::endl;
std::cout << "Closest energy: " << m_energies[idx_energy] << std::endl;
std::cout << "Interpolated energy: " << (1 - t) * m_energies[idx_energy] + t * m_energies[idx_energy + 1] << std::endl;
return (1 - t) * dielectric_function_interpolated_low + t * dielectric_function_interpolated_high;
}

} // namespace bz_mesh
141 changes: 55 additions & 86 deletions src/BZ_MESH/impact_ionization.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,12 +48,42 @@ void ImpactIonization::read_dielectric_file(const std::string& filename) {
m_dielectric_mesh.read_dielectric_file(filename);
}

void ImpactIonization::compute_eigenstates() {
int nb_bands_to_use = 4;
void ImpactIonization::interp_test_dielectric_function(std::string filename) {
double eps = 1e-6;
double x0 = eps;
double y0 = eps;
double z0 = eps;
double x1 = 3.0;
double y1 = 3.0;
double z1 = 3.0;
int nx = 20;
int ny = 20;
int nz = 20;
double dx = (x1 - x0) / nx;
double dy = (y1 - y0) / ny;
double dz = (z1 - z0) / nz;
std::ofstream file(filename);
for (int i = 0; i < nx; ++i) {
for (int j = 0; j < ny; ++j) {
for (int k = 0; k < nz; ++k) {
double x = x0 + i * dx;
double y = y0 + j * dy;
double z = z0 + k * dz;
vector3 position(x, y, z);
complex_d epsilon = m_dielectric_mesh.interpolate_dielectric_function(position, 0.0102);
file << x << ", " << y << ", " << z << ", " << epsilon.real() << ", " << epsilon.imag() << std::endl;
std::cout << "Position: " << position << " epsilon: " << epsilon << std::endl;
}
}
}
file.close();
}

void ImpactIonization::compute_eigenstates(int nb_threads) {
int nb_bands_to_use = 16;
bz_mesh::BZ_States my_bz_mesh(m_material);
my_bz_mesh.set_nb_bands(nb_bands_to_use);
EmpiricalPseudopotential::BandStructure band_structure{};
int nb_threads = 1;
int nb_nearest_neighbors = 10;
bool nonlocal_epm = false;
bool enable_soc = false;
Expand All @@ -68,7 +98,7 @@ void ImpactIonization::compute_eigenstates() {
const vector3 b3 = {1.0, 1.0, -1.0};

// std::cout << "k: " << k << std::endl;

double factor = 2.0 * M_PI / m_material.get_lattice_constant_meter();
// test k + G
std::vector<int> list_n_k = {0, 1, -1, 2, -2, 3, -3, 4, -4};
// std::vector<int> list_n_k = {0, 1, -1, 2, -2};
Expand All @@ -81,7 +111,9 @@ void ImpactIonization::compute_eigenstates() {
if (G_BZ.norm() > m_max_radius_G0_BZ) {
continue;
}
std::cout << "G_BZ: " << G_BZ << std::endl;
std::cout << "G_BZ: " << G_BZ << " --> " << G_BZ.norm() << std::endl;
G_BZ = G_BZ * factor;

auto ptr_BZ_states = std::make_unique<BZ_States>(m_material);
ptr_BZ_states->set_nb_bands(nb_bands_to_use);
ptr_BZ_states->set_basis_vectors(basis);
Expand All @@ -96,88 +128,25 @@ void ImpactIonization::compute_eigenstates() {
}

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
constexpr int nb_valence_bands = 3;
constexpr int nb_conduction_bands = 4;

return rate;
}
std::size_t nb_vtx = m_list_BZ_states[0]->get_list_vertices().size();

// /**
// * @brief Compute the direct impact ionization matrix element (Ma) for a given set of indices.
// * k2 is not an input parameter, it is computed from k1, k1_prime and k2_prime, using the conservation of momentum.
// * k2 = k1_prime + k2_prime - k1
// *
// *
// * @param idx_n1
// * @param idx_n1_prime
// * @param idx_n2
// * @param idx_n2_prime
// * @param idx_k1
// * @param idx_k1_prime
// * @param idx_k2_prime
// * @return double
// */
// std::array<complex_d, 2> ImpactIonization::compute_direct_indirect_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,
// int idx_k2_prime) const {
// vector3 k1 = m_list_vertices[idx_k1].get_position();
// vector3 k2 = m_list_vertices[idx_k2].get_position();
// vector3 k1_prime = m_list_vertices[idx_k1_prime].get_position();
// vector3 k2_prime = m_list_vertices[idx_k2_prime].get_position();
// complex_d Ma = 0.0;
// complex_d Mb = 0.0;

// 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;

// 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 g2_prime = {basis_vectors[idx_g2_prime].X * fourier_factor,
// basis_vectors[idx_g2_prime].Y * fourier_factor,
// basis_vectors[idx_g2_prime].Z * fourier_factor};
// vector3 q_a = k1_prime + g1_prime - k1 - g1;
// vector3 q_b = k2_prime + g2_prime - k1 - g1;
// auto state_k1 = m_eigenvectors_k[idx_k1](idx_n1, idx_g1);
// auto state_k2 = m_eigenvectors_k[idx_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);

// 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);
// complex_d dielectric_func_a = get_dielectric_function(q_a, omega_a);
// complex_d pre_factor_a = e_charge * e_charge / (eps_zero * dielectric_func_a * q_a.norm() * q_a.norm() *
// volume); complex_d dielectric_func_b = get_dielectric_function(q_b, omega_b); complex_d core_overlap =
// std::conj(state_k1_prime) * std::conj(state_k2_prime) * state_k1 * state_k2; complex_d micro_matrix_element_a =
// core_overlap * pre_factor_a; complex_d pre_factor_b = e_charge * e_charge / (eps_zero * dielectric_func_b *
// q_b.norm() * q_b.norm() * volume); complex_d micro_matrix_element_b = core_overlap * pre_factor_b;

// Ma += micro_matrix_element_a;
// Mb += micro_matrix_element_b;
// }
// }
// }
// }

// return {Ma, Mb};
// }
for (int n1_prime = 0; n1_prime < nb_valence_bands; ++n1_prime) {
for (int n2 = 0; n2 < nb_conduction_bands; ++n2) {
for (int n2_prime = 0; n2_prime < nb_conduction_bands; ++n2_prime) {
for (int k1_prime = 0; k1_prime < nb_vtx; ++k1_prime) {

}
}
}
}





}

} // namespace bz_mesh
3 changes: 2 additions & 1 deletion src/BZ_MESH/impact_ionization.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,10 +67,11 @@ class ImpactIonization {
public:
ImpactIonization(const EmpiricalPseudopotential::Material& material, const std::string& initial_mesh_path);
void read_dielectric_file(const std::string& filename);
void interp_test_dielectric_function(std::string filename);

double get_max_radius_G0_BZ() const { return m_max_radius_G0_BZ; }
void set_max_radius_G0_BZ(double max_radius_G0_BZ) { m_max_radius_G0_BZ = max_radius_G0_BZ; }
void compute_eigenstates();
void compute_eigenstates(int nb_threads = 1);

std::array<complex_d, 2> compute_direct_indirect_impact_ionization_matrix_element(int idx_n1,
int idx_n1_prime,
Expand Down
16 changes: 14 additions & 2 deletions src/BZ_MESH/mesh_tetra.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#pragma once

#include <array>
#include <complex>
#include <functional>
#include <iostream>
#include <numeric>
Expand All @@ -20,6 +21,7 @@
#include "bbox_mesh.hpp"
#include "mesh_vertex.hpp"


namespace bz_mesh {

class Tetra {
Expand Down Expand Up @@ -133,8 +135,9 @@ class Tetra {

std::array<double, 8> get_mean_electron_phonon_rates(int band_index) const;

double interpolate_scalar_at_position(const std::array<double, 4>& barycentric_coordinates,
const std::vector<double>& scalar_field) const;
double interpolate_scalar_at_position(const std::array<double, 4>& barycentric_coordinates,
const std::vector<double>& scalar_field) const;

vector3 interpolate_vector_at_position(const std::array<double, 4>& barycentric_coordinates,
const std::vector<vector3>& vector_field) const;
template <typename T>
Expand All @@ -145,6 +148,15 @@ class Tetra {
}
return interpolated_value;
}
std::complex<double> interpolate_at_position(const std::array<double, 4>& barycentric_coordinates,
const std::vector<std::complex<double>>& field) const {
std::complex<double> interpolated_value = 0.0;
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<double>(5, 0.0); }

Expand Down

0 comments on commit 88efb75

Please sign in to comment.