From 4b585363f432954a003c7a891c1944a213dc5a80 Mon Sep 17 00:00:00 2001 From: RemiHelleboid Date: Fri, 15 Sep 2023 22:40:45 +0200 Subject: [PATCH] Change all files to add the SOC option --- apps/BandsOnBZ.cpp | 9 +- apps/EmpiricalPseudoPotentialMain.cpp | 30 +++++- apps/epsilon.cpp | 3 +- apps/fullstates.cpp | 15 ++- apps/mpi_BandsOnBZ.cpp | 6 +- src/EPP/BandStructure.cpp | 8 +- src/EPP/BandStructure.h | 11 ++- src/EPP/SpinOrbitFunctional.cpp | 130 ++++++++++++-------------- 8 files changed, 117 insertions(+), 95 deletions(-) diff --git a/apps/BandsOnBZ.cpp b/apps/BandsOnBZ.cpp index a3d0cca..dccaa1f 100644 --- a/apps/BandsOnBZ.cpp +++ b/apps/BandsOnBZ.cpp @@ -33,6 +33,7 @@ int main(int argc, char* argv[]) { 10, "int"); TCLAP::SwitchArg arg_enable_nonlocal_correction("C", "nonlocal-correction", "Enable the non-local-correction for the EPM model", false); + TCLAP::SwitchArg arg_enable_soc("s", "soc", "Enable the spin-orbit coupling for the EPM model", false); TCLAP::SwitchArg arg_cond_band_zero("z", "MinCondZero", "Shift the conduction band minimum to 0 eV", false); TCLAP::ValueArg arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int"); cmd.add(arg_mesh_file); @@ -42,9 +43,8 @@ int main(int argc, char* argv[]) { cmd.add(arg_nearest_neighbors); cmd.add(arg_nb_threads); cmd.add(arg_enable_nonlocal_correction); + cmd.add(arg_enable_soc); cmd.add(arg_cond_band_zero); - - cmd.parse(argc, argv); @@ -52,6 +52,7 @@ int main(int argc, char* argv[]) { const std::string file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials.yaml"; materials.load_material_parameters(file_material_parameters); bool enable_nonlocal_correction = arg_enable_nonlocal_correction.isSet(); + bool enable_soc = arg_enable_soc.isSet(); Options my_options; my_options.materialName = arg_material.getValue(); @@ -71,7 +72,8 @@ int main(int argc, char* argv[]) { auto start = std::chrono::high_resolution_clock::now(); EmpiricalPseudopotential::BandStructure my_bandstructure; - my_bandstructure.Initialize(mat, my_options.nrLevels, mesh_kpoints, my_options.nearestNeighbors, enable_nonlocal_correction); + my_bandstructure + .Initialize(mat, my_options.nrLevels, mesh_kpoints, my_options.nearestNeighbors, enable_nonlocal_correction, arg_enable_soc); if (my_options.nrThreads > 1) { my_bandstructure.Compute_parallel(my_options.nrThreads); @@ -80,7 +82,6 @@ int main(int argc, char* argv[]) { } my_bandstructure.AdjustValues(arg_cond_band_zero.getValue()); - auto end = std::chrono::high_resolution_clock::now(); auto total_time_count = std::chrono::duration_cast(end - start).count(); diff --git a/apps/EmpiricalPseudoPotentialMain.cpp b/apps/EmpiricalPseudoPotentialMain.cpp index 8915945..dde82a6 100644 --- a/apps/EmpiricalPseudoPotentialMain.cpp +++ b/apps/EmpiricalPseudoPotentialMain.cpp @@ -50,6 +50,7 @@ int compute_path_mat(const EmpiricalPseudopotential::Material& material, unsigned int nb_bands, unsigned int nearestNeighbors, bool enable_non_local_correction, + bool enable_soc, const std::string& result_dir, bool call_python_plot) { Options my_options; @@ -57,8 +58,13 @@ int compute_path_mat(const EmpiricalPseudopotential::Material& material, my_options.nrPoints = nb_points; my_options.nrLevels = nb_bands; EmpiricalPseudopotential::BandStructure my_bandstructure; - my_bandstructure - .Initialize(material, my_options.nrLevels, path, my_options.nrPoints, my_options.nearestNeighbors, enable_non_local_correction); + my_bandstructure.Initialize(material, + my_options.nrLevels, + path, + my_options.nrPoints, + my_options.nearestNeighbors, + enable_non_local_correction, + enable_soc); my_bandstructure.Compute(); my_bandstructure.AdjustValues(); std::cout << "Time to compute: " << my_bandstructure.get_computation_time_s() << " s" << std::endl; @@ -80,6 +86,7 @@ int compute_all_mat(EmpiricalPseudopotential::Materials list_materials, int nearestNeighbors, int nrPoints, bool enable_non_local_correction, + bool enable_soc, int nb_threads, const std::string& result_dir, bool call_python_plot) { @@ -97,8 +104,13 @@ int compute_all_mat(EmpiricalPseudopotential::Materials list_materials, } std::cout << std::endl; EmpiricalPseudopotential::BandStructure my_bandstructure; - my_bandstructure - .Initialize(mat, my_options.nrLevels, path, my_options.nrPoints, my_options.nearestNeighbors, enable_non_local_correction); + my_bandstructure.Initialize(mat, + my_options.nrLevels, + path, + my_options.nrPoints, + my_options.nearestNeighbors, + enable_non_local_correction, + enable_soc); my_bandstructure.Compute_parallel(my_options.nrThreads); my_bandstructure.AdjustValues(); @@ -124,6 +136,7 @@ int compute_all_path_all_mat(EmpiricalPseudopotential::Materials list_materials, int nearestNeighbors, int nrPoints, bool enable_non_local_correction, + bool enable_soc, int nb_threads, const std::string& result_dir, bool call_python_plot) { @@ -148,7 +161,8 @@ int compute_all_path_all_mat(EmpiricalPseudopotential::Materials list_materials, my_options.paths[path_index], my_options.nrPoints, my_options.nearestNeighbors, - enable_non_local_correction); + enable_non_local_correction, + enable_soc); my_bandstructure.Compute_parallel(my_options.nrThreads); my_bandstructure.AdjustValues(); @@ -190,6 +204,7 @@ int main(int argc, char* argv[]) { TCLAP::ValueArg arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int"); TCLAP::ValueArg arg_res_dir("r", "resultdir", "directory to store the results.", false, "./", "str"); TCLAP::SwitchArg arg_enable_nonlocal_correction("C", "nonlocal-correction", "Enable the non-local-correction for the EPM model", false); + TCLAP::SwitchArg arg_enable_soc("S", "soc", "Enable the spin-orbit coupling for the EPM model", false); TCLAP::SwitchArg all_path_mat("A", "all", "Compute the band structure on all the paths for all the materials", false); TCLAP::SwitchArg plot_with_python("P", "plot", "Call a python script after the computation to plot the band structure.", false); cmd.add(arg_path_sym_points); @@ -200,6 +215,7 @@ int main(int argc, char* argv[]) { cmd.add(arg_nb_threads); cmd.add(arg_res_dir); cmd.add(arg_enable_nonlocal_correction); + cmd.add(arg_enable_soc); cmd.add(all_path_mat); cmd.add(plot_with_python); @@ -237,6 +253,7 @@ int main(int argc, char* argv[]) { bool call_python_plot = plot_with_python.isSet(); bool enable_nonlocal_correction = arg_enable_nonlocal_correction.isSet(); + bool enable_soc = arg_enable_soc.isSet(); if (all_path_mat.getValue()) { std::cout << "Compute the band structure on all the paths for all the materials" << std::endl; @@ -245,6 +262,7 @@ int main(int argc, char* argv[]) { arg_nearest_neighbors.getValue(), arg_nb_points.getValue(), enable_nonlocal_correction, + enable_soc, arg_nb_threads.getValue(), arg_res_dir.getValue(), call_python_plot); @@ -257,6 +275,7 @@ int main(int argc, char* argv[]) { arg_nb_bands.getValue(), arg_nearest_neighbors.getValue(), enable_nonlocal_correction, + enable_soc, arg_res_dir.getValue(), call_python_plot); } else if (!arg_material.isSet() && arg_path_sym_points.isSet()) { @@ -267,6 +286,7 @@ int main(int argc, char* argv[]) { arg_nearest_neighbors.getValue(), arg_nb_points.getValue(), enable_nonlocal_correction, + enable_soc, arg_nb_threads.getValue(), arg_res_dir.getValue(), call_python_plot); diff --git a/apps/epsilon.cpp b/apps/epsilon.cpp index eac9c78..2b671b4 100644 --- a/apps/epsilon.cpp +++ b/apps/epsilon.cpp @@ -151,6 +151,7 @@ int main(int argc, char** argv) { int Nkz = config["Nkz"].as(); bool nonlocal_epm = false; + bool enable_soc = false; if (config["nonlocal"]) { nonlocal_epm = config["nonlocal"].as(); } @@ -184,7 +185,7 @@ int main(int argc, char** argv) { EmpiricalPseudopotential::Material current_material = materials.materials.at("Si"); EmpiricalPseudopotential::BandStructure band_structure{}; - band_structure.Initialize(current_material, nb_bands, {}, nb_nearest_neighbors, nonlocal_epm); + band_structure.Initialize(current_material, nb_bands, {}, nb_nearest_neighbors, nonlocal_epm, enable_soc); EmpiricalPseudopotential::DielectricFunction MyDielectricFunc(current_material, band_structure.get_basis_vectors(), nb_bands); double shift = 0.0; diff --git a/apps/fullstates.cpp b/apps/fullstates.cpp index 146bdc5..b56a8c8 100644 --- a/apps/fullstates.cpp +++ b/apps/fullstates.cpp @@ -41,6 +41,7 @@ int main(int argc, char *argv[]) { cmd.parse(argc, argv); bool nonlocal_epm = false; + bool enable_soc = false; EmpiricalPseudopotential::Materials materials; std::string file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-local.yaml"; if (nonlocal_epm) { @@ -64,30 +65,28 @@ int main(int argc, char *argv[]) { EmpiricalPseudopotential::BandStructure band_structure{}; int nb_nearest_neighbors = 10; - band_structure.Initialize(current_material, nb_bands_to_use, {}, nb_nearest_neighbors, nonlocal_epm); + band_structure.Initialize(current_material, nb_bands_to_use, {}, nb_nearest_neighbors, nonlocal_epm, enable_soc); auto basis = band_structure.get_basis_vectors(); 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-10, 0.0, 0.0}; my_bz_mesh.compute_shifted_eigenstates(q_shift, my_options.nrThreads); std::cout << "\n\n" << 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; std::vector list_energy; - double min_energy = 0.0; - double max_energy = 20.0; - double energy_step = 0.01; + 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); } @@ -96,9 +95,9 @@ int main(int argc, char *argv[]) { 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(); + 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; + return 0; } \ No newline at end of file diff --git a/apps/mpi_BandsOnBZ.cpp b/apps/mpi_BandsOnBZ.cpp index 8153e8b..064fbe5 100644 --- a/apps/mpi_BandsOnBZ.cpp +++ b/apps/mpi_BandsOnBZ.cpp @@ -71,6 +71,7 @@ int main(int argc, char* argv[]) { 10, "int"); TCLAP::SwitchArg arg_enable_nonlocal_correction("C", "nonlocal-correction", "Enable the non-local-correction for the EPM model", false); + TCLAP::SwitchArg arg_enable_soc("S", "soc", "Enable the spin-orbit coupling for the EPM model", false); TCLAP::SwitchArg arg_cond_band_zero("z", "MinCondZero", "Shift the conduction band minimum to 0 eV", false); TCLAP::ValueArg arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int"); cmd.add(arg_mesh_file); @@ -80,6 +81,7 @@ int main(int argc, char* argv[]) { cmd.add(arg_nearest_neighbors); cmd.add(arg_nb_threads); cmd.add(arg_enable_nonlocal_correction); + cmd.add(arg_enable_soc); cmd.add(arg_cond_band_zero); cmd.parse(argc, argv); @@ -168,8 +170,10 @@ int main(int argc, char* argv[]) { } bool enable_nonlocal_correction = arg_enable_nonlocal_correction.getValue(); + bool enable_soc = arg_enable_soc.getValue(); EmpiricalPseudopotential::BandStructure my_bandstructure; - my_bandstructure.Initialize(mat, my_options.nrLevels, Chunk_list_k_points, my_options.nearestNeighbors, enable_nonlocal_correction); + my_bandstructure + .Initialize(mat, my_options.nrLevels, Chunk_list_k_points, my_options.nearestNeighbors, enable_nonlocal_correction, enable_soc); my_bandstructure.Compute(); my_bandstructure.AdjustValues(arg_cond_band_zero.getValue()); diff --git a/src/EPP/BandStructure.cpp b/src/EPP/BandStructure.cpp index 080a793..db03b0f 100644 --- a/src/EPP/BandStructure.cpp +++ b/src/EPP/BandStructure.cpp @@ -57,13 +57,15 @@ void BandStructure::Initialize(const Material& material, const std::vector& path, unsigned int nbPoints, unsigned int nearestNeighborsNumber, - bool enable_non_local_correction) { + bool enable_non_local_correction, + bool enable_soc) { m_material = material; m_nb_bands = nb_bands; m_path = path; m_nb_points = nbPoints; m_nearestNeighborsNumber = nearestNeighborsNumber; m_enable_non_local_correction = enable_non_local_correction; + m_enable_spin_orbit_coupling = enable_soc; m_kpoints.clear(); m_results.clear(); m_kpoints.reserve(m_nb_points); @@ -85,12 +87,14 @@ void BandStructure::Initialize(const Material& material, std::size_t nb_bands, std::vector> list_k_points, unsigned int nearestNeighborsNumber, - bool enable_non_local_correction) { + bool enable_non_local_correction, + bool enable_soc) { m_material = material; m_nb_bands = nb_bands; m_nb_points = list_k_points.size(); m_nearestNeighborsNumber = nearestNeighborsNumber; m_enable_non_local_correction = enable_non_local_correction; + m_enable_spin_orbit_coupling = enable_soc; m_kpoints.clear(); m_results.clear(); diff --git a/src/EPP/BandStructure.h b/src/EPP/BandStructure.h index b809ba1..6aebe55 100644 --- a/src/EPP/BandStructure.h +++ b/src/EPP/BandStructure.h @@ -23,13 +23,15 @@ class BandStructure { * @param nrPoints * @param nearestNeighborsNumber * @param enable_non_local_correction + * @param enable_soc */ void Initialize(const Material& material, std::size_t nb_bands, const std::vector& path, unsigned int nrPoints, unsigned int nearestNeighborsNumber, - bool enable_non_local_correction); + bool enable_non_local_correction, + bool enable_soc); /** * @brief Initialize the band structure with the given material and the k-points on which we want @@ -40,19 +42,21 @@ class BandStructure { * @param list_k_points * @param nearestNeighborsNumber * @param enable_non_local_correction + * @param enable_soc */ void Initialize(const Material& material, std::size_t nb_bands, std::vector> list_k_points, unsigned int nearestNeighborsNumber, - bool enable_non_local_correction); + bool enable_non_local_correction, + bool enable_soc); const std::vector& GetPath() const { return m_path; } std::string get_path_as_string() const; unsigned int GetPointsNumber() const { return static_cast(m_kpoints.size()); } void Compute(); void Compute_parallel(int nb_threads); - double AdjustValues(bool minConductionBandToZero=false); + double AdjustValues(bool minConductionBandToZero = false); void print_results() const; std::string path_band_filename() const; @@ -82,6 +86,7 @@ class BandStructure { unsigned int m_nb_bands; unsigned int m_nearestNeighborsNumber; bool m_enable_non_local_correction; + bool m_enable_spin_orbit_coupling = false; std::vector> basisVectors; std::vector> m_kpoints; diff --git a/src/EPP/SpinOrbitFunctional.cpp b/src/EPP/SpinOrbitFunctional.cpp index ba6e700..75578cb 100644 --- a/src/EPP/SpinOrbitFunctional.cpp +++ b/src/EPP/SpinOrbitFunctional.cpp @@ -17,76 +17,64 @@ 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; - - } - - +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