From 00273793ef823ba9d134c9e392d92e638cd76e83 Mon Sep 17 00:00:00 2001 From: RemiHelleboid Date: Fri, 15 Sep 2023 22:23:04 +0200 Subject: [PATCH] Continue SOC implementation --- apps/fullstates.cpp | 27 +++++- python/plots/plot_eps_vs_energy.py | 9 +- src/BZ_MESH/bz_mesh.cpp | 11 ++- src/BZ_MESH/bz_mesh.hpp | 2 +- src/BZ_MESH/bz_states.cpp | 92 +++++++++++++++++- src/BZ_MESH/bz_states.hpp | 14 +++ src/BZ_MESH/mesh_tetra.cpp | 2 + src/BZ_MESH/mesh_tetra.hpp | 13 ++- src/EPP/Constants.hpp | 1 + src/EPP/DielectricFunction.cpp | 2 +- src/EPP/SpinOrbitFunctional.cpp | 149 ++++++++++++++++++----------- src/EPP/SpinOrbitFunctional.hpp | 115 +++++++++++----------- 12 files changed, 311 insertions(+), 126 deletions(-) diff --git a/apps/fullstates.cpp b/apps/fullstates.cpp index 6a2c020..146bdc5 100644 --- a/apps/fullstates.cpp +++ b/apps/fullstates.cpp @@ -69,19 +69,36 @@ int main(int argc, char *argv[]) { my_bz_mesh.set_basis_vectors(basis); my_bz_mesh.read_mesh_geometry_from_msh_file(arg_mesh_file.getValue()); + std::cout << "Mesh volume: " << my_bz_mesh.compute_mesh_volume() << std::endl; + + my_bz_mesh.compute_eigenstates(my_options.nrThreads); - Vector3D q_shift = Vector3D{1e-12, 0.0, 0.0}; + Vector3D q_shift = Vector3D{1e-10, 0.0, 0.0}; my_bz_mesh.compute_shifted_eigenstates(q_shift, my_options.nrThreads); + std::cout << "\n\n" << std::endl; - std::cout << "Mesh volume: " << my_bz_mesh.compute_mesh_volume() << std::endl; auto end = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed_seconds = end - start; std::cout << "Elapsed time: " << elapsed_seconds.count() << "s" << std::endl; - my_bz_mesh.export_full_eigenstates(); - - return 0; + std::vector list_energy; + double min_energy = 0.0; + double max_energy = 20.0; + double energy_step = 0.01; + for (double energy = min_energy; energy <= max_energy + energy_step; energy += energy_step) { + list_energy.push_back(energy); + } + double eta_smearing = 0.01; + std::cout << "Number of energies to compute: " << list_energy.size() << std::endl; + auto start2 = std::chrono::high_resolution_clock::now(); + my_bz_mesh.compute_dielectric_function(list_energy, eta_smearing, my_options.nrThreads); + my_bz_mesh.export_dielectric_function("./TEST_DIELECTRIC_FUNCTION_"); + auto end2 = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed_seconds2 = end2 - start2; + std::cout << "Time Dielectric Function : " << elapsed_seconds2.count() << "s" << std::endl; + + return 0; } \ No newline at end of file diff --git a/python/plots/plot_eps_vs_energy.py b/python/plots/plot_eps_vs_energy.py index 1879944..f6ced20 100644 --- a/python/plots/plot_eps_vs_energy.py +++ b/python/plots/plot_eps_vs_energy.py @@ -20,7 +20,14 @@ def plot_dielectric_function_vs_energy(filename: str): - energies, eps_r, eps_i = np.loadtxt(filename, unpack=True, skiprows=1, delimiter=',') + # energies, eps_r, eps_i = np.loadtxt(filename, unpack=True, skiprows=1, delimiter=',') + data = np.loadtxt(filename, unpack=True, skiprows=1, delimiter=',') + energies = data[0] + eps_r = data[1] + if len(data) == 3: + eps_i = data[2] + else: + eps_i = np.zeros_like(eps_r) fig, ax = plt.subplots() ax.plot(energies, eps_r, label="Real", c='b') ax.plot(energies, eps_i, "-", label="Imaginary", c='r') diff --git a/src/BZ_MESH/bz_mesh.cpp b/src/BZ_MESH/bz_mesh.cpp index cb8e25f..3e1b841 100644 --- a/src/BZ_MESH/bz_mesh.cpp +++ b/src/BZ_MESH/bz_mesh.cpp @@ -53,7 +53,11 @@ void MeshBZ::read_mesh_geometry_from_msh_file(const std::string& filename) { m_list_vertices.reserve(size_nodes_tags); double lattice_constant = m_material.get_lattice_constant_meter(); - double normalization_factor = 2.0 * M_PI / lattice_constant; + std::cout << "Lattice const: " << lattice_constant << std::endl; + std::cout << "V: " << std::pow(2.0 * M_PI, 3) / std::pow(lattice_constant, 3.0) << std::endl; + + double normalization_factor = 2.0 * M_PI / lattice_constant; + // double normalization_factor = 1.0; for (std::size_t index_vertex = 0; index_vertex < size_nodes_tags; ++index_vertex) { m_list_vertices.push_back(Vertex(index_vertex, normalization_factor * nodeCoords[3 * index_vertex], @@ -164,9 +168,10 @@ void MeshBZ::compute_energy_gradient_at_tetras() { double MeshBZ::compute_mesh_volume() const { double total_volume = 0.0; for (auto&& tetra : m_list_tetrahedra) { - total_volume += fabs(tetra.get_signed_volume()); + total_volume += std::fabs(tetra.get_signed_volume()); + // std::cout << "Tetra " << tetra.get_index() << " volume: " << tetra.get_signed_volume() << std::endl; } - total_volume *= (1.0 / pow(2.0 * M_PI, 3.0)); + // total_volume *= (1.0 / pow(2.0 * M_PI, 3.0)); return total_volume; } diff --git a/src/BZ_MESH/bz_mesh.hpp b/src/BZ_MESH/bz_mesh.hpp index f9424b9..3ae7e67 100644 --- a/src/BZ_MESH/bz_mesh.hpp +++ b/src/BZ_MESH/bz_mesh.hpp @@ -93,7 +93,7 @@ class MeshBZ { double compute_mesh_volume() const; double compute_iso_surface(double iso_energy, int band_index) const; double compute_dos_at_energy_and_band(double iso_energy, int band_index) const; - double compute_overlapp_integral_impact_ionization_electrons(double energy); + double compute_overlap_integral_impact_ionization_electrons(double energy); std::vector> compute_dos_band_at_band(int band_index, double min_energy, diff --git a/src/BZ_MESH/bz_states.cpp b/src/BZ_MESH/bz_states.cpp index 22c3a3a..85dc9a1 100644 --- a/src/BZ_MESH/bz_states.cpp +++ b/src/BZ_MESH/bz_states.cpp @@ -49,6 +49,7 @@ void BZ_States::compute_eigenstates(int nb_threads) { } void BZ_States::compute_shifted_eigenstates(const Vector3D& q_shift, int nb_threads) { + m_q_shift = q_shift; double normalization_factor = 2.0 * M_PI / m_material.get_lattice_constant_meter(); const bool m_nonlocal_epm = false; const bool keep_eigenvectors = true; @@ -64,8 +65,8 @@ void BZ_States::compute_shifted_eigenstates(const Vector3D& q_shift, int auto k_point = Vector3D(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; k_point += q_shift; - 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); @@ -76,6 +77,95 @@ void BZ_States::compute_shifted_eigenstates(const Vector3D& q_shift, int } } +/** + * @brief Compute the dielectric function for a given list of energies and a given smearing. + * The integration is performed by summing the contribution of each tetrahedron to the dielectric function. + * + * @param energies + * @param eta_smearing + * @param nb_threads + */ +void BZ_States::compute_dielectric_function(const std::vector& list_energies, double eta_smearing, int nb_threads) { + m_list_energies = list_energies; + const int index_first_conduction_band = 4; + std::size_t nb_tetra = m_list_tetrahedra.size(); + m_dielectric_function_real.resize(list_energies.size()); + constexpr double one_fourth = 1.0 / 4.0; + + std::vector dielectric_function_real_at_energies(list_energies.size(), 0.0); + +#pragma omp parallel for schedule(dynamic) num_threads(nb_threads) + 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 list_idx_vertices = m_list_tetrahedra[idx_tetra].get_list_indices_vertices(); + const std::array& list_vertices = m_list_tetrahedra[idx_tetra].get_list_vertices(); + double volume_tetra = std::fabs(m_list_tetrahedra[idx_tetra].compute_signed_volume()); + // std::cout << "Volume tetra = " << volume_tetra << std::endl; + std::vector sum_dielectric_function_real_tetra_at_energies(list_energies.size(), 0.0); + // Loop over the vertices of the tetrahedron + for (std::size_t idx_vertex = 0; idx_vertex < 4; ++idx_vertex) { + std::size_t index_k = list_idx_vertices[idx_vertex]; + 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) { + double overlap_integral = pow(std::fabs(m_eigenvectors_k_plus_q[index_k] + .col(idx_conduction_band) + .adjoint() + .dot(m_eigenvectors_k[index_k].col(idx_valence_band))), + 2); + double delta_energy = m_eigenvalues_k_plus_q[index_k][idx_conduction_band] - m_eigenvalues_k[index_k][idx_valence_band]; + for (std::size_t index_energy = 0; index_energy < list_energies.size(); ++index_energy) { + double energy = list_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; + sum_dielectric_function_real_tetra_at_energies[index_energy] += overlap_integral * total_factor; + } + } + } + } + for (std::size_t index_energy = 0; index_energy < list_energies.size(); ++index_energy) { + sum_dielectric_function_real_tetra_at_energies[index_energy] *= volume_tetra * one_fourth; + dielectric_function_real_at_energies[index_energy] += sum_dielectric_function_real_tetra_at_energies[index_energy]; + } + } + + double q_squared = m_q_shift.Length() * m_q_shift.Length(); + std::cout << "\nq_squared = " << q_squared << std::endl; + // double normalization_1 = (EmpiricalPseudopotential::Constants::q * EmpiricalPseudopotential::Constants::q * 4.0 * M_PI); + double normalization_1 = (M_PI) / (2.0 / std::pow(2.0 * M_PI, 3)); + // double normalization_2 = compute_mesh_volume(); + double normalization_3 = (EmpiricalPseudopotential::Constants::q * EmpiricalPseudopotential::Constants::q); + std::cout << std::endl; + std::cout << "Normalization factor 1: " << normalization_1 << std::endl; + std::cout << "Normalization factor 2: " << normalization_3 << std::endl; + std::cout << "Normalization factor q: " << 1.0 / q_squared << std::endl; + for (std::size_t index_energy = 0; index_energy < list_energies.size(); ++index_energy) { + m_dielectric_function_real[index_energy] = dielectric_function_real_at_energies[index_energy]; + // std::cout << "Energy = " << list_energies[index_energy] << " eV, dielectric function = " << + // m_dielectric_function_real[index_energy] << std::endl; m_dielectric_function_real[index_energy] *= normalization_1; + m_dielectric_function_real[index_energy] *= normalization_3; + // m_dielectric_function_real[index_energy] /= normalization_2; + m_dielectric_function_real[index_energy] /= q_squared; + m_dielectric_function_real[index_energy] += 1.0; + } + std::cout << "E = 0 --> " << dielectric_function_real_at_energies[0] << std::endl; + std::cout << "E = 0 --> " << dielectric_function_real_at_energies[0] / q_squared << std::endl; + std::cout << "E = 0 --> " << dielectric_function_real_at_energies[0] / (q_squared * compute_mesh_volume()) << std::endl; + std::cout << "E = 0 --> " << m_dielectric_function_real[0] << std::endl; +} + +// Export the dielectric function to a file in the format (energy, dielectric function) (csv format). +void BZ_States::export_dielectric_function(const std::string& prefix) const { + std::ofstream dielectric_function_file(prefix + "_dielectric_function.csv"); + dielectric_function_file << "Energy (eV), Dielectric function" << std::endl; + for (std::size_t index_energy = 0; index_energy < m_dielectric_function_real.size(); ++index_energy) { + dielectric_function_file << m_list_energies[index_energy] << ", " << m_dielectric_function_real[index_energy] << std::endl; + } + dielectric_function_file.close(); +} + void BZ_States::export_full_eigenstates() const { std::filesystem::remove_all("eigenstates"); std::filesystem::create_directory("eigenstates"); diff --git a/src/BZ_MESH/bz_states.hpp b/src/BZ_MESH/bz_states.hpp index 38a749b..cee17c3 100644 --- a/src/BZ_MESH/bz_states.hpp +++ b/src/BZ_MESH/bz_states.hpp @@ -29,6 +29,16 @@ class BZ_States : public MeshBZ { std::vector m_eigenvectors_k; std::vector m_eigenvectors_k_plus_q; + Vector3D m_q_shift; + std::vector m_list_energies; + + /** + * @brief Real part of the dielectric function. + * m_dielectric_function_real[idx_energy] is the real part of the dielectric function at the energy m_energies[idx_energy]. + * + */ + std::vector m_dielectric_function_real; + public: BZ_States(const EmpiricalPseudopotential::Material& material) : MeshBZ(material) {} @@ -37,7 +47,11 @@ class BZ_States : public MeshBZ { void compute_eigenstates(int nb_threads = 1); void compute_shifted_eigenstates(const Vector3D& q_shift, int nb_threads = 1); + const std::vector& get_energies() const { return m_list_energies; } + void set_energies(const std::vector& energies) { m_list_energies = energies; } + void compute_dielectric_function(const std::vector& energies, double eta_smearing, int nb_threads = 1); + void export_dielectric_function(const std::string& prefix) const; double compute_direct_impact_ionization_matrix_element(int idx_n1, int idx_n1_prime, diff --git a/src/BZ_MESH/mesh_tetra.cpp b/src/BZ_MESH/mesh_tetra.cpp index 9e8230c..45f5c24 100644 --- a/src/BZ_MESH/mesh_tetra.cpp +++ b/src/BZ_MESH/mesh_tetra.cpp @@ -79,6 +79,8 @@ void Tetra::compute_min_max_energies_at_bands() { * @return double */ double Tetra::compute_signed_volume() const { + // std::cout << m_list_edges[0] << std::endl; + // std::cout << m_list_edges[1] << std::endl; return (1.0 / 6.0) * scalar_triple_product(m_list_edges[0], m_list_edges[1], m_list_edges[2]); } diff --git a/src/BZ_MESH/mesh_tetra.hpp b/src/BZ_MESH/mesh_tetra.hpp index 0666722..8f5b9aa 100644 --- a/src/BZ_MESH/mesh_tetra.hpp +++ b/src/BZ_MESH/mesh_tetra.hpp @@ -48,7 +48,7 @@ class Tetra { * The sign depends on the "orientation" of the tetrahedra. * */ - double m_signed_volume; + double m_signed_volume = 0.0; /** * @brief Number of conduction bands. @@ -84,6 +84,17 @@ class Tetra { Tetra(std::size_t index, const std::array& list_vertices); void compute_min_max_energies_at_bands(); + std::size_t get_index() const { return m_index; } + const std::array& get_list_vertices() const { return m_list_vertices; } + std::array get_list_indices_vertices() const { + return {m_list_vertices[0]->get_index(), + m_list_vertices[1]->get_index(), + m_list_vertices[2]->get_index(), + m_list_vertices[3]->get_index()}; + } + std::array get_list_edges() const { return m_list_edges; } + std::size_t get_nb_bands() const { return m_nb_bands; } + std::array get_band_energies_at_vertices(std::size_t index_band) const; double compute_signed_volume() const; diff --git a/src/EPP/Constants.hpp b/src/EPP/Constants.hpp index 60f97e0..b822653 100644 --- a/src/EPP/Constants.hpp +++ b/src/EPP/Constants.hpp @@ -21,6 +21,7 @@ constexpr double Ryd_to_eV = 13.6; constexpr double h = 6.62606957e-34; constexpr double h_bar = 6.62606957e-34 / (2 * M_PI); constexpr double m0 = 9.10938356e-31; +// constexpr double q = 1.602176634e-19; constexpr double q = 1.6e-19; constexpr double eps_zero = 8.854187e-12; constexpr double k_b = 1.38064852e-23; diff --git a/src/EPP/DielectricFunction.cpp b/src/EPP/DielectricFunction.cpp index b3dacb7..af3bc11 100644 --- a/src/EPP/DielectricFunction.cpp +++ b/src/EPP/DielectricFunction.cpp @@ -71,7 +71,7 @@ void DielectricFunction::generate_k_points_grid(std::size_t Nx, std::size_t Ny, } /** - * @brief Compute the energy and wavevector dependent dielectric function. + * @brief Compute the energy and wave vector dependent dielectric function. * The formula used is the one from the paper " * * @param eta_smearing diff --git a/src/EPP/SpinOrbitFunctional.cpp b/src/EPP/SpinOrbitFunctional.cpp index 637adba..ba6e700 100644 --- a/src/EPP/SpinOrbitFunctional.cpp +++ b/src/EPP/SpinOrbitFunctional.cpp @@ -1,59 +1,92 @@ -/** - * @file SpinOrbitDunctional.cpp - * @author remzerrr (remi.helleboid@gmail.com) - * @brief - * @version 0.1 - * @date 2022-09-30 - * - * @copyright Copyright (c) 2022 - * - */ - -#include "SpinOrbitFunctional.hpp" - -#include -#include -#include - -namespace EmpiricalPseudopotential { - -// double SpinOrbitCorrection::compute_B2_cation(const Vector3D& K) const { -// const Vector3D kappa = K * (Constants::bohr_radius / m_soc_parameters.m_zeta_cation); -// double B2 = 1.0 / std::pow((1.0 + kappa * kappa), 3.0); -// return B2; -// } - -// double SpinOrbitCorrection::compute_B2_anion(const Vector3D& K) const { -// const Vector3D kappa = K * (Constants::bohr_radius / m_soc_parameters.m_zeta_anion); -// double B2 = 1.0 / std::pow((1.0 + kappa * kappa), 3.0); -// return B2; -// } - -// double SpinOrbitCorrection::compute_B3_cation(const Vector3D& K) const { -// const Vector3D kappa = K * (Constants::bohr_radius / m_soc_parameters.m_zeta_cation); -// double B3 = (5 - kappa * kappa) / (5.0 * pow((1.0 * kappa * kappa), 4.0)); -// return B3; -// } - -// double SpinOrbitCorrection::compute_B3_anion(const Vector3D& K) const { -// const Vector3D kappa = K * (Constants::bohr_radius / m_soc_parameters.m_zeta_anion); -// double B3 = (5 - kappa * kappa) / (5.0 * pow((1.0 * kappa * kappa), 4.0)); -// return B3; -// } - -// double SpinOrbitCorrection::compute_B4_cation(const Vector3D& K) const { -// const Vector3D kappa = K * (Constants::bohr_radius / m_soc_parameters.m_zeta_cation); -// double B4 = (5.0 - 3.0 * kappa * kappa) / (5.0 * pow((1.0 + kappa * kappa), 5.0)); - -// return B4; -// } - -// double SpinOrbitCorrection::compute_B4_anion(const Vector3D& K) const { -// const Vector3D kappa = K * (Constants::bohr_radius / m_soc_parameters.m_zeta_anion); -// double B4 = (5.0 - 3.0 * kappa * kappa) / (5.0 * pow((1.0 + kappa * kappa), 5.0)); -// return B4; -// } - - - +/** + * @file SpinOrbitDunctional.cpp + * @author remzerrr (remi.helleboid@gmail.com) + * @brief + * @version 0.1 + * @date 2022-09-30 + * + * @copyright Copyright (c) 2022 + * + */ + +#include "SpinOrbitFunctional.hpp" + +#include +#include +#include + +namespace EmpiricalPseudopotential { + + double SpinOrbitCorrection::compute_B2_cation(const Vector3D& K) const { + double kappa = K.Length() * (Constants::bohr_radius / m_soc_parameters.m_radial_extent_cation); + double B2 = 1.0 / std::pow((1.0 + kappa * kappa), 3.0); + return B2; + } + + double SpinOrbitCorrection::compute_B2_anion(const Vector3D& K) const { + double kappa = K.Length() * (Constants::bohr_radius / m_soc_parameters.m_radial_extent_anion); + double B2 = 1.0 / std::pow((1.0 + kappa * kappa), 3.0); + return B2; + } + + double SpinOrbitCorrection::compute_B3_cation(const Vector3D& K) const { + double kappa = K.Length() * (Constants::bohr_radius / m_soc_parameters.m_radial_extent_cation); + double B3 = (5 - kappa * kappa) / (5.0 * pow((1.0 + kappa * kappa), 4.0)); + return B3; + } + + double SpinOrbitCorrection::compute_B3_anion(const Vector3D& K) const { + double kappa = K.Length() * (Constants::bohr_radius / m_soc_parameters.m_radial_extent_anion); + double B3 = (5 - kappa * kappa) / (5.0 * pow((1.0 + kappa * kappa), 4.0)); + return B3; + } + + double SpinOrbitCorrection::compute_B4_cation(const Vector3D& K) const { + double kappa = K.Length() * (Constants::bohr_radius / m_soc_parameters.m_radial_extent_cation); + double B4 = (5.0 - 3.0 * kappa * kappa) / (5.0 * pow((1.0 + kappa * kappa), 5.0)); + return B4; + } + + double SpinOrbitCorrection::compute_B4_anion(const Vector3D& K) const { + double kappa = K.Length() * (Constants::bohr_radius / m_soc_parameters.m_radial_extent_anion); + double B4 = (5.0 - 3.0 * kappa * kappa) / (5.0 * pow((1.0 + kappa * kappa), 5.0)); + return B4; + } + + + double SpinOrbitCorrection::compute_lambda_1(const Vector3D& K, const Vector3D& Kp) const { + + double lambda_1 = m_soc_parameters.m_mu * compute_B2_cation(K) * compute_B2_cation(Kp); + return lambda_1; + + } + + double SpinOrbitCorrection::compute_lambda_2(const Vector3D& K, const Vector3D& Kp) const { + + double lambda_2 = m_soc_parameters.m_alpha * m_soc_parameters.m_mu * compute_B2_anion(K) * compute_B2_anion(Kp); + return lambda_2; + + } + + + double SpinOrbitCorrection::compute_lambda_sym(const Vector3D& K, const Vector3D& Kp) const { + + double lambda_1 = compute_lambda_1(K, Kp); + double lambda_2 = compute_lambda_2(K, Kp); + double lambda_sym = (lambda_1 + lambda_2)/2; + return lambda_sym; + + } + + double SpinOrbitCorrection::compute_lambda_antisym(const Vector3D& K, const Vector3D& Kp) const { + + double lambda_1 = compute_lambda_1(K, Kp); + double lambda_2 = compute_lambda_2(K, Kp); + double lambda_antisym = (lambda_1 - lambda_2) / 2; + return lambda_antisym; + + } + + + } // namespace EmpiricalPseudopotential \ No newline at end of file diff --git a/src/EPP/SpinOrbitFunctional.hpp b/src/EPP/SpinOrbitFunctional.hpp index e4694f5..605bcfb 100644 --- a/src/EPP/SpinOrbitFunctional.hpp +++ b/src/EPP/SpinOrbitFunctional.hpp @@ -1,55 +1,60 @@ -/** - * @file SpinOrbitFunctions.hpp - * @author remzerrr (remi.helleboid@gmail.com) - * @brief - * @version 0.1 - * @date 2022-09-30 - * - * @copyright Copyright (c) 2022 - * - */ - -#pragma once - -#include -#include -#include -#include - -#include "Constants.hpp" -#include "Material.h" -#include "SpinOrbitParameters.hpp" -#include "Vector3D.h" - -namespace EmpiricalPseudopotential { - -/** - * @brief Functor class to compute the spin-orbit Hamiltonian. - * The name of the functions follow the notation of the paper: - * Pötz, W. & Vogl, P. Theory of optical-phonon deformation - * potentials in tetrahedral semiconductors. Phys. Rev. B 24, 2025–2037 (1981). - * - */ -class SpinOrbitCorrection { - protected: - SpinOrbitParameters m_soc_parameters; - Material m_material; - - public: - SpinOrbitCorrection() = delete; - SpinOrbitCorrection(const SpinOrbitParameters& parameters, const Material& material) - : m_soc_parameters(parameters), - m_material(material) {} - - double compute_B2_cation(const Vector3D& K) const; - double compute_B2_anion(const Vector3D& K) const; - double compute_B3_cation(const Vector3D& K) const; - double compute_B3_anion(const Vector3D& K) const; - double compute_B4_cation(const Vector3D& K) const; - double compute_B4_anion(const Vector3D& K) const; - - double compute_lambda_1(const Vector3D& K) const; - double compute_lambda_2(const Vector3D& K) const; -}; - -} // namespace EmpiricalPseudopotential +/** + * @file SpinOrbitFunctions.hpp + * @author remzerrr (remi.helleboid@gmail.com) + * @brief + * @version 0.1 + * @date 2022-09-30 + * + * @copyright Copyright (c) 2022 + * + */ + +#pragma once + +#include +#include +#include +#include + +#include "Constants.hpp" +#include "Material.h" +#include "SpinOrbitParameters.hpp" +#include "Vector3D.h" + +namespace EmpiricalPseudopotential { + +/** + * @brief Functor class to compute the spin-orbit Hamiltonian. + * The name of the functions follow the notation of the paper: + * Pötz, W. & Vogl, P. Theory of optical-phonon deformation + * potentials in tetrahedral semiconductors. Phys. Rev. B 24, 2025–2037 (1981). + * + */ +class SpinOrbitCorrection { + protected: + SpinOrbitParameters m_soc_parameters; + Material m_material; + + public: + SpinOrbitCorrection() = delete; + SpinOrbitCorrection(const SpinOrbitParameters& parameters, const Material& material) + : m_soc_parameters(parameters), + m_material(material) {} + + double compute_B2_cation(const Vector3D& K) const; + double compute_B2_anion(const Vector3D& K) const; + double compute_B3_cation(const Vector3D& K) const; + double compute_B3_anion(const Vector3D& K) const; + double compute_B4_cation(const Vector3D& K) const; + double compute_B4_anion(const Vector3D& K) const; + + double compute_lambda_1(const Vector3D& K, const Vector3D& Kp) const; + double compute_lambda_2(const Vector3D& K, const Vector3D& Kp) const; + + + double compute_lambda_sym(const Vector3D& K, const Vector3D& Kp) const; + double compute_lambda_antisym(const Vector3D& K, const Vector3D& Kp) const; + +}; + +} // namespace EmpiricalPseudopotential