Okkkay

README.md

@@ -3,6 +3,7 @@ Band structure and DOS computation using empirical pseudopotentials on the full
 * Standard __EPM over high symmetry k-points__ in the Brillouin zone.
 * Calculation of the __band structure and DOS on a mesh of k-points in the Brillouin zone__.
 * __Nonlocal corrections__ to the EPM.
+* __Spin-orbit coupling__ (SOC) for the EPM band structure.
 * MPI and OpenMP __parallelization__.
 
 [![Build & Unit Test](https://github.com/RemiHelleboid/EmpiricalPseudopotential/actions/workflows/build_code.yaml/badge.svg)](https://github.com/RemiHelleboid/EmpiricalPseudopotential/actions/workflows/build_code.yaml)
@@ -22,6 +23,10 @@ You can do three types of calculations:
 __Compute the electronic band structure over a path of high-symmetry points (e.g. $L \Gamma XWKULWXK \Gamma$) for a given material, and plot the results.__ 
 ![Band structure of Silicon](doc/band_structure_Si.png)
 
+The SOC can be included in the computation:
+
+![Band structure of Ge with SOC](doc/band_structure_Ge_soc.png)
+
 ---
 
 __Compute the electronic band structure over all k-points of an input mesh of the Brillouin Zone (or a fraction of it). The result can then be visualized, for example, through iso-energy surface.__  

apps/EmpiricalPseudoPotentialMain.cpp

@@ -235,12 +235,12 @@ int main(int argc, char* argv[]) {
     }
 
     EmpiricalPseudopotential::Materials materials;
-    std::string                         file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials.yaml";
+    std::string                         file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-local.yaml";
     if (!arg_enable_nonlocal_correction.isSet()) {
         file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-local.yaml";
     }
     materials.load_material_parameters(file_material_parameters);
-    // materials.print_material_parameters();
+    materials.print_material_parameters();
 
     print_arguments(path_list,
                     arg_material.getValue(),

apps/fullstates.cpp

@@ -43,7 +43,7 @@ int main(int argc, char *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";
+    std::string                         file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-cohen.yaml";
     if (nonlocal_epm) {
         file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials.yaml";
     }
@@ -74,9 +74,10 @@ int main(int argc, char *argv[]) {
 
     my_bz_mesh.compute_eigenstates(my_options.nrThreads);
 
-    Vector3D<double> q_shift = Vector3D<double>{1e-10, 0.0, 0.0};
+    Vector3D<double> q_shift = Vector3D<double>{1e-14, 0.0, 0.0};
     my_bz_mesh.compute_shifted_eigenstates(q_shift, my_options.nrThreads);
     std::cout << "\n\n" << std::endl;
+    my_bz_mesh.export_full_eigenstates();
 
     auto end = std::chrono::high_resolution_clock::now();
 

apps/mpi_BandsOnBZ.cpp

@@ -87,7 +87,7 @@ int main(int argc, char* argv[]) {
     cmd.parse(argc, argv);
 
     EmpiricalPseudopotential::Materials materials;
-    const std::string                   file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials.yaml";
+    const std::string                   file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-cohen.yaml";
     materials.load_material_parameters(file_material_parameters);
 
     Options my_options;
@@ -195,6 +195,9 @@ int main(int argc, char* argv[]) {
     std::vector<double> all_energies_all_bands;
     if (process_rank == MASTER) {
         all_energies_all_bands.resize(number_k_vectors * number_bands);
+        std::cout << "number_k_vectors: " << number_k_vectors << std::endl;
+        std::cout << "number_bands: " << number_bands << std::endl;
+        std::cout << "all_energies_all_bands.size(): " << all_energies_all_bands.size() << std::endl;
     }
 
     MPI_Gatherv(chunk_list_energies.data(),
@@ -222,7 +225,7 @@ int main(int argc, char* argv[]) {
     if (process_rank == MASTER) {
         std::filesystem::path in_path(mesh_filename);
         std::string out_file_bands = in_path.stem().replace_extension("").string() + "_MPI_" + my_bandstructure.path_band_filename();
-        // my_mesh.add_all_bands_on_mesh(out_file_bands + "_all_bands.msh", all_energies_all_bands, number_bands);
+        my_mesh.add_all_bands_on_mesh(out_file_bands + "_all_bands.msh", all_energies_all_bands, number_bands);
         my_mesh.export_bands_as_csv(all_energies_all_bands, number_bands);
     }
 

parameter_files/materials-local.yaml

@@ -17,6 +17,13 @@ materials:
       V3A: 0.0
       V4A: 0.0
       V11A: 0.0
+    spin-orbit-parameters:
+      period_anion: 3
+      period_cation: 3
+      radial_extent_anion: 4.60
+      radial_extent_cation: 4.60
+      alpha_soc: 1.00
+      mu_soc: 0.0
 
   - name: Germanium
     symbol: Ge
@@ -28,6 +35,13 @@ materials:
       V3A: 0.0
       V4A: 0.0
       V11A: 0.0
+    spin-orbit-parameters:
+      period_anion: 4
+      period_cation: 4
+      radial_extent_anion: 5.34
+      radial_extent_cation: 5.34
+      alpha_soc: 1.00
+      mu_soc: 0.00069
 
   - name: Tin
     symbol: Sn

python/plots/plot_band_structure.py

@@ -191,10 +191,11 @@ def plot_band_structure(filename: str, ax, index_plot, nb_bands=10):
 
     lines = [Line2D([0], [0], color='k', lw=0.5, ls=lstyle)
              for lstyle in ["-", "--"]]
-    labels = ["non-local", "local"]
+    # labels = ["non-local", "local"]
+    # labels = ["non-local", "local"]
 
-    ax.legend(lines, labels, loc='lower right', fancybox=True,  frameon=True,
-              edgecolor='k', facecolor='w', fontsize=6, framealpha=0.75)
+    # ax.legend(lines, labels, loc='lower right', fancybox=True,  frameon=True,
+    #           edgecolor='k', facecolor='w', fontsize=6, framealpha=0.75)
     fig.tight_layout()
     fig.savefig(f"{OUT_DIR}/band_structure_{material}.png", dpi=300)
 

src/BZ_MESH/bz_states.cpp

@@ -31,7 +31,6 @@ 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;
         auto k_point    = Vector3D<double>(m_list_vertices[idx_k].get_position().x(),
@@ -59,7 +58,6 @@ void BZ_States::compute_shifted_eigenstates(const Vector3D<double>& q_shift, int
     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+q = " << idx_k << "/" << m_list_vertices.size() << std::flush;
         auto k_point = Vector3D<double>(m_list_vertices[idx_k].get_position().x(),
@@ -93,13 +91,13 @@ 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);
-
-#pragma omp parallel for schedule(dynamic) num_threads(nb_threads)
+    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();
         const std::array<Vertex*, 4>& list_vertices     = m_list_tetrahedra[idx_tetra].get_list_vertices();
         double                        volume_tetra      = std::fabs(m_list_tetrahedra[idx_tetra].compute_signed_volume());
+        total_volume += volume_tetra;
         // std::cout << "Volume tetra = " << volume_tetra << std::endl;
         std::vector<double> sum_dielectric_function_real_tetra_at_energies(list_energies.size(), 0.0);
         // Loop over the vertices of the tetrahedron
@@ -130,30 +128,17 @@ void BZ_States::compute_dielectric_function(const std::vector<double>& list_ener
             dielectric_function_real_at_energies[index_energy] += sum_dielectric_function_real_tetra_at_energies[index_energy];
         }
     }
+    std::cout << "\n";
+    std::cout << "Total volume: " << total_volume << std::endl;
 
     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;
+    double pre_factor = 4.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] = 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;
+        m_dielectric_function_real[index_energy] = 1.0 + pre_factor * dielectric_function_real_at_energies[index_energy] / total_volume;
     }
-    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;
+    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).

src/EPP/BandStructure.cpp

@@ -71,6 +71,12 @@ void BandStructure::Initialize(const Material&                 material,
     m_kpoints.reserve(m_nb_points);
     m_results.reserve(m_nb_points);
 
+    if (m_enable_spin_orbit_coupling) {
+        m_nb_bands *= 2;
+        std::cout << "Spin-orbit coupling enabled. Number of bands doubled to " << m_nb_bands << std::endl;
+        m_material.get_spin_orbit_parameters().print_parameters();
+    }
+
     if (!GenerateBasisVectors(nearestNeighborsNumber)) {
         throw std::runtime_error("BandStructure::Initialize: GenerateBasisVectors failed");
     }
@@ -96,6 +102,13 @@ void BandStructure::Initialize(const Material&               material,
     m_enable_non_local_correction = enable_non_local_correction;
     m_enable_spin_orbit_coupling  = enable_soc;
 
+    if (m_enable_spin_orbit_coupling) {
+        m_nb_bands *= 2;
+        std::cout << "Spin-orbit coupling enabled. Number of bands doubled to " << m_nb_bands << std::endl;
+        m_material.get_spin_orbit_parameters().print_parameters();
+
+    }
+
     m_kpoints.clear();
     m_results.clear();
     m_kpoints = list_k_points;
@@ -115,7 +128,8 @@ void BandStructure::Compute() {
     for (unsigned int i = 0; i < m_nb_points; ++i) {
         // std::cout << "\rComputing band structure at point " << i + 1 << "/" << m_nb_points << std::flush;
         // std::cout << "Computing band structure at point " << m_kpoints[i] << std::endl;
-        hamiltonian.SetMatrix(m_kpoints[i], m_enable_non_local_correction);
+
+        hamiltonian.SetMatrix(m_kpoints[i], m_enable_non_local_correction, m_enable_spin_orbit_coupling);
         hamiltonian.Diagonalize();
 
         const Eigen::VectorXd& eigenvals = hamiltonian.eigenvalues();
@@ -146,7 +160,7 @@ void BandStructure::Compute_parallel(int nb_threads) {
 #pragma omp parallel for schedule(dynamic) num_threads(nb_threads)
     for (unsigned int index_k = 0; index_k < m_nb_points; ++index_k) {
         int tid = omp_get_thread_num();
-        hamiltonian_per_thread[tid].SetMatrix(m_kpoints[index_k], m_enable_non_local_correction);
+        hamiltonian_per_thread[tid].SetMatrix(m_kpoints[index_k], m_enable_non_local_correction, m_enable_spin_orbit_coupling);
         hamiltonian_per_thread[tid].Diagonalize();
 
         const Eigen::VectorXd& eigenvals = hamiltonian_per_thread[tid].eigenvalues();

src/EPP/DielectricFunction.cpp

@@ -54,11 +54,11 @@ void DielectricFunction::generate_k_points_random(std::size_t nb_points) {
 
 void DielectricFunction::generate_k_points_grid(std::size_t Nx, std::size_t Ny, std::size_t Nz, double shift, bool irreducible_wedge) {
     m_kpoints.clear();
-    double min = -1.0 - shift;
+    double min = -0.0 - shift;
     double max = 1.0 + shift;
-    for (std::size_t i = 0; i < Nx; ++i) {
-        for (std::size_t j = 0; j < Ny; ++j) {
-            for (std::size_t k = 0; k < Nz; ++k) {
+    for (std::size_t i = 0; i < Nx - 1; ++i) {
+        for (std::size_t j = 0; j < Ny - 1; ++j) {
+            for (std::size_t k = 0; k < Nz - 1; ++k) {
                 Vector3D<double> k_vect(min + (max - min) * i / static_cast<double>(Nx - 1),
                                         min + (max - min) * j / static_cast<double>(Ny - 1),
                                         min + (max - min) * k / static_cast<double>(Nz - 1));
@@ -99,6 +99,7 @@ void DielectricFunction::compute_dielectric_function(double eta_smearing) {
         for (std::size_t index_k = m_offset_k_index; index_k < m_offset_k_index + m_nb_kpoints; ++index_k) {
             if (index_q == 0) {
                 auto k_vect = m_kpoints[index_k];
+                std::cout << "Computing eigenvalues for k = " << k_vect << std::endl;
                 hamiltonian_k.SetMatrix(k_vect, m_nonlocal_epm);
                 hamiltonian_k.Diagonalize(keep_eigenvectors);
                 m_eigenvalues_k[index_k]  = hamiltonian_k.eigenvalues();
@@ -115,10 +116,12 @@ 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) {
-                    double overlap_integral =
-                        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];
+                    double overlap_integral = 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 << "Energy val: " << m_eigenvalues_k[index_k][idx_valence_band] << 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 =
@@ -181,7 +184,7 @@ DielectricFunction DielectricFunction::merge_results(DielectricFunction
         }
     }
     // Renormalization
-    double renormalization = 1.0 / static_cast<double>(total_number_kpoints);
+    double       renormalization = 1.0 / static_cast<double>(total_number_kpoints);
     const double factor_discrete = 2.0 / std::pow(2.0 * M_PI, 3.0);
     std::cout << "Renormalization: " << renormalization << std::endl;
     for (std::size_t index_q = 0; index_q < total_dielectric_function.size(); ++index_q) {

src/EPP/Hamiltonian.cpp

@@ -74,19 +74,36 @@ void Hamiltonian::SetMatrix(const Vector3D<double>& k, bool add_non_local_correc
     }
 
     if (enable_soc) {
-        std::size_t      matrix_size = matrix.rows();
-        std::size_t new_matrix_size = 2 * matrix_size;
-        Eigen::MatrixXcd ZerosMat = Eigen::MatrixXcd::Zero(matrix_size, matrix_size);
-        Eigen::MatrixXcd new_matrix(new_matrix_size, new_matrix_size);
-        new_matrix << matrix, ZerosMat, ZerosMat, matrix;
-        matrix = new_matrix;
-
+        // std::cout << "Spin-orbit correction enabled" << std::endl;
+        SpinOrbitParameters SpinParams = m_material.get_spin_orbit_parameters();
+        SpinOrbitCorrection soc_correction(m_material, SpinParams);
+        std::size_t         matrix_size    = matrix.rows();
+        Eigen::MatrixXcd    UpDownMatrix   = Eigen::MatrixXcd::Zero(matrix_size, matrix_size);
+        Eigen::MatrixXcd    DownUpMatrix   = Eigen::MatrixXcd::Zero(matrix_size, matrix_size);
+        Eigen::MatrixXcd    UpUpMatrix     = matrix;
+        Eigen::MatrixXcd    DownDownMatrix = matrix;
+        for (unsigned int i = 0; i < matrix_size; ++i) {
+            for (unsigned int j = 0; j < matrix_size; ++j) {
+                Vector3D<double>                          k_vector_i = (k + m_basisVectors[i]);
+                Vector3D<double>                          k_vector_j = (k + m_basisVectors[j]);
+                Eigen::Matrix<std::complex<double>, 2, 2> soc_contribution =
+                    soc_correction.compute_soc_contribution(k_vector_i, k_vector_j, m_basisVectors[i], m_basisVectors[j], tau);
+                // std::cout << soc_contribution << std::endl;
+                UpUpMatrix(i, j) += soc_contribution(0, 0);
+                UpDownMatrix(i, j) += soc_contribution(1, 0);
+                DownUpMatrix(i, j) += soc_contribution(0, 1);
+                DownDownMatrix(i, j) += soc_contribution(1, 1);
+            }
+        }
+        matrix.resize(2 * matrix_size, 2 * matrix_size);
+        matrix << UpUpMatrix, UpDownMatrix, DownUpMatrix, DownDownMatrix;
     }
 }
 
 void Hamiltonian::Diagonalize(bool keep_eigenvectors) {
     solver.compute(matrix, keep_eigenvectors ? Eigen::ComputeEigenvectors : Eigen::EigenvaluesOnly);
     if (solver.info() != Eigen::Success) {
+        std::cout << matrix << std::endl;
         throw std::runtime_error("Eigenvalue decomposition failed");
     }
 }

src/EPP/Material.cpp

@@ -65,6 +65,8 @@ void Materials::load_material_parameters(const std::string& filename) {
         if (node_spin_orbit_parameters) {
             materials[symbol].populate_spin_orbit_parameters(node_spin_orbit_parameters);
             materials[symbol].set_is_spin_orbit_parameters_populated(true);
+            std::cout << "Spin-orbit parameters for " << symbol << " are set." << std::endl;
+            materials[symbol].get_spin_orbit_parameters().print_parameters();
 
         }
     }