Skip to content

Commit

Permalink
Up
Browse files Browse the repository at this point in the history
  • Loading branch information
@RemiHelleboid
RemiHelleboid committed Dec 26, 2022
1 parent e861bc4 commit 08e61cf
Show file tree
Hide file tree
Showing 3 changed files with 34 additions and 22 deletions.
24 changes: 19 additions & 5 deletions src/BZ_MESH/bz_states.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,20 +34,34 @@ void BZ_States::compute_eigenstates(int nb_threads) {
#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;
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());
k_point = k_point * 1.0 / normalization_factor;
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());
k_point = k_point * 1.0 / normalization_factor;
auto idx_thread = omp_get_thread_num();
hamiltonian_per_thread[idx_thread].SetMatrix(k_point, m_nonlocal_epm);
hamiltonian_per_thread[idx_thread].Diagonalize(keep_eigenvectors);
m_eigenvalues_k[idx_k] = hamiltonian_per_thread[idx_thread].eigenvalues();
m_eigenvectors_k[idx_k] = hamiltonian_per_thread[idx_thread].get_eigenvectors();
auto nb_rows = m_eigenvectors_k[idx_k].rows();
auto nb_rows = m_eigenvectors_k[idx_k].rows();
m_eigenvectors_k[idx_k].conservativeResize(nb_rows, m_nb_bands);
}
}

// double BZ_States::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>& k2,
// const Vector3D<double>& k1_prime,
// const Vector3D<double>& k2_prime) const {
// double matrix_element = 0.0;
// double normalization_factor = 2.0 * M_PI / m_material.get_lattice_constant_meter();
// const bool m_nonlocal_epm = false;
// const bool keep_eigenvectors = true;
// }

void BZ_States::export_full_eigenstates() const {
std::filesystem::remove_all("eigenstates");
std::filesystem::create_directory("eigenstates");
Expand Down
16 changes: 13 additions & 3 deletions src/BZ_MESH/bz_states.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@

#include <Eigen/Dense>

#include "bz_mesh.hpp"
#include "Material.h"
#include "bz_mesh.hpp"

namespace bz_mesh {

Expand All @@ -25,16 +25,26 @@ class BZ_States : public MeshBZ {
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;

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<Vector3D<int>>& basis_vectors) { m_basisVectors = basis_vectors; }

void compute_eigenstates(int nb_threads = 1);
void compute_dielectric_function(const std::vector<double>& energies, double eta_smearing, int nb_threads = 1);

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>& k2,
const Vector3D<double>& k1_prime,
const Vector3D<double>& k2_prime) const;

void export_full_eigenstates() const;
};
Expand Down
16 changes: 2 additions & 14 deletions src/EPP/DielectricFunction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,10 +101,8 @@ void DielectricFunction::compute_dielectric_function(double eta_smearing) {
auto k_vect = m_kpoints[index_k];
hamiltonian_k.SetMatrix(k_vect, m_nonlocal_epm);
hamiltonian_k.Diagonalize(keep_eigenvectors);

m_eigenvalues_k[index_k] = hamiltonian_k.eigenvalues();
m_eigenvectors_k[index_k] = hamiltonian_k.get_eigenvectors();
// std::cout << "k = " << k_vect << " -> " << m_eigenvectors_k[index_k] << std::endl;
// Keep only firsts columns
auto nb_rows = m_eigenvectors_k[index_k].rows();
m_eigenvectors_k[index_k].conservativeResize(nb_rows, m_nb_bands);
Expand All @@ -117,22 +115,17 @@ void DielectricFunction::compute_dielectric_function(double eta_smearing) {
std::vector<double> list_k_sum(m_energies.size());
for (int idx_conduction_band = index_first_conduction_band; idx_conduction_band < m_nb_bands; ++idx_conduction_band) {
for (int idx_valence_band = 0; idx_valence_band < index_first_conduction_band; ++idx_valence_band) {
// const auto& eigen_vect_k = m_eigenvectors_k[index_k].col(idx_valence_band);
// std::cout << "eigen_vect_k norm = " << eigenvectors_k_plus_q.mean() << std::endl;
double overlap_integral =
pow(std::abs(eigenvectors_k_plus_q.col(idx_conduction_band).dot(m_eigenvectors_k[index_k].col(idx_valence_band))),
pow(std::abs(eigenvectors_k_plus_q.col(idx_conduction_band).adjoint().dot(m_eigenvectors_k[index_k].col(idx_valence_band))),
2);
double delta_energy = eigenvalues_k_plus_q[idx_conduction_band] - m_eigenvalues_k[index_k][idx_valence_band];
// std::cout << "Overlap integral = " << overlap_integral << std::endl;
// std::cout << "delta_energy = " << delta_energy << std::endl;
for (std::size_t index_energy = 0; index_energy < m_energies.size(); ++index_energy) {
double energy = m_energies[index_energy];
double factor_1 =
(delta_energy - energy) / ((delta_energy - energy) * (delta_energy - energy) + eta_smearing * eta_smearing);
double factor_2 =
(delta_energy + energy) / ((delta_energy + energy) * (delta_energy + energy) + eta_smearing * eta_smearing);
double total_factor = factor_1 + factor_2;
// std::cout << "overlap_integral = " << overlap_integral << std::endl;
list_k_sum[index_energy] += overlap_integral * total_factor;
}
}
Expand All @@ -144,17 +137,14 @@ void DielectricFunction::compute_dielectric_function(double eta_smearing) {
}
for (std::size_t index_energy = 0; index_energy < m_energies.size(); ++index_energy) {
list_total_sum[index_energy] += list_k_sum[index_energy];
// std::cout << "Energy = " << m_energies[index_energy] << " -> epsilon = " << list_total_sum[index_energy] << std::endl;
}
}
std::vector<double> list_epsilon(m_energies.size());
for (std::size_t index_energy = 0; index_energy < m_energies.size(); ++index_energy) {
list_epsilon[index_energy] = list_total_sum[index_energy];
// std::cout << "Energy = " << m_energies[index_energy] << " -> epsilon = " << list_epsilon[index_energy] << std::endl;
}
auto end = std::chrono::high_resolution_clock::now();
// std::cout << std::endl;
std::cout << "Time elapsed: " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / 1000.0 << " s"
std::cout << "\nTime elapsed: " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / 1000.0 << " s"
<< std::endl;
start = std::chrono::high_resolution_clock::now();
m_dielectric_function_real.push_back(list_epsilon);
Expand All @@ -178,8 +168,6 @@ DielectricFunction DielectricFunction::merge_results(DielectricFunction
std::cout << "Number total kpoint to the merge : " << total_number_kpoints << std::endl;
// Add-up the k-contributions.
for (std::size_t index_instance = 0; index_instance < dielectric_function_results.size(); ++index_instance) {
// const double ratio_number_k_points =
// static_cast<double>(nb_kpoints_per_instance[index_instance]) / static_cast<double>(total_number_kpoints);
// Add the contributions to epsilon for each q-point and energy.
for (std::size_t index_q = 0; index_q < dielectric_function_results[index_instance].size(); ++index_q) {
if (index_instance == 0) {
Expand Down

0 comments on commit 08e61cf

Please sign in to comment.