up doc

RemiHelleboid · Feb 8, 2024 · 07ca3e7 · 07ca3e7
1 parent 36e3209
commit 07ca3e7
Show file tree

Hide file tree

Showing 9 changed files with 110 additions and 56 deletions.
diff --git a/README.md b/README.md
@@ -20,12 +20,12 @@ This repository was initially a fork from : [EmpiricalPseudopotential](https://g
 
 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.__ 
+__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)
+![Band structure of Ge with SOC](doc/GeSOC.png)
 
 ---
 
@@ -40,6 +40,10 @@ __Compute the density of states over the all Brillouin Zone.__
 
 ---
 
+__Dielectric function computation (q-dependent)__
+![Band structure of Silicon](/doc/dielectric_func_Si.png)
+![Maps of the dielectric function](/doc/dielectric_map_Si.png)
+
 ## Build and Compilation
 
 ### Dependencies

diff --git a/apps/mpi_BandsOnBZ.cpp b/apps/mpi_BandsOnBZ.cpp
@@ -63,6 +63,7 @@ int main(int argc, char* argv[]) {
     TCLAP::ValueArg<std::string> arg_mesh_file("f", "meshfile", "Name to print", true, "bz.msh", "string");
     TCLAP::ValueArg<std::string> arg_material("m", "material", "Symbol of the material to use (Si, Ge, GaAs, ...)", true, "Si", "string");
     TCLAP::ValueArg<std::string> arg_outfile("o", "outfile", "Name of the output file", false, "", "string");
+    TCLAP::ValueArg<std::string> arg_material_parameters("p", "parameters", "Name of the file containing the material parameters", false, "", "string");
     TCLAP::ValueArg<int>         arg_nb_bands("b", "nbands", "Number of bands to compute", false, 12, "int");
     TCLAP::ValueArg<int>         arg_nearest_neighbors("n",
                                                "nearestNeighbors",
@@ -78,6 +79,7 @@ int main(int argc, char* argv[]) {
     cmd.add(arg_material);
     cmd.add(arg_nb_bands);
     cmd.add(arg_outfile);
+    cmd.add(arg_material_parameters);
     cmd.add(arg_nearest_neighbors);
     cmd.add(arg_nb_threads);
     cmd.add(arg_enable_nonlocal_correction);
@@ -87,7 +89,15 @@ 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-cohen.yaml";
+
+    std::string file_material_parameters = arg_material_parameters.getValue();
+    if (file_material_parameters.empty()) {
+        file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-cohen.yaml";
+        if (!arg_enable_nonlocal_correction.getValue()) {
+            file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials-local.yaml";
+        }
+    }
+
     materials.load_material_parameters(file_material_parameters);
 
     Options my_options;

diff --git a/doc/GeSOC.png b/doc/GeSOC.png
diff --git a/doc/dielectric_func_Si.png b/doc/dielectric_func_Si.png
diff --git a/doc/dielectric_map_Si.png b/doc/dielectric_map_Si.png
diff --git a/python/plots/plot_band_structure.py b/python/plots/plot_band_structure.py
@@ -202,7 +202,8 @@ def plot_band_structure(filename: str, ax, index_plot, nb_bands=10):
     if args.show_plot:
         plt.show()
 
-    ax.set_ylim(-2.0, band_gap + 1.0)
+    ax.set_ylim(-4.0, band_gap + 4.0)
+    fig.tight_layout()
     fig.savefig(f"{OUT_DIR}/band_structure_{material}_zoom.png", dpi=300)
 
     # plot_band_structure(FILE_PATH, OUT_DIR, args.nb_bands)
diff --git a/python/plots/plot_eps_vs_energy.py b/python/plots/plot_eps_vs_energy.py
@@ -8,12 +8,11 @@
 from scipy.interpolate import CubicSpline
 from pathlib import Path
 
-
-import sys
+import sys, glob
 
 try:
     import scienceplots
-    plt.style.use(['science', 'high-vis'])
+    plt.style.use(['science', 'muted'])
 except Exception:
     print('Could not load style sheet')
     None
@@ -22,34 +21,42 @@
 mpl.rcParams['figure.figsize'] = [3.5, 2.8]
 
 
-def plot_dielectric_function_vs_energy(filename: str, ax):
-    # 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)
-    label = Path(filename).stem.split('_')[2::]
-    label = ["{:.0f}".format(float(x)) for x in label]
-    s = "$\mathbf{k} = (" + " \ ".join(label) + ")$"
-    ax.plot(energies, eps_r, label=s)
-    # ax.plot(energies, eps_i, "-", label="Imaginary", c='r')
-    abs_eps = np.sqrt(eps_r**2 + eps_i**2)
-    # ax.plot(energies, abs_eps, label=s)
-
+def plot_dielectric_function_vs_energy(dirname):
+    list_files = glob.glob(f"{dirname}/*.csv")
+    list_q = []
+    for filename in list_files:
+        filestem = Path(filename).stem
+        qz = float(filestem.split('_')[-1])
+        qy = float(filestem.split('_')[-2])
+        qx = float(filestem.split('_')[-3])
+        print(qx, qy, qz)
+        list_q.append(np.sqrt(qx**2 + qy**2 + qz**2))
+
+    list_files = [x for _, x in sorted(zip(list_q, list_files))]
+    cmap = plt.get_cmap('jet')
+    norm = mpl.colors.Normalize(vmin=min(list_q), vmax=max(list_q))
 
-if __name__ == '__main__':
-    list_files = sys.argv[1:]
     fig, ax = plt.subplots()
-    for filename in list_files:
-        plot_dielectric_function_vs_energy(filename, ax)
+    for idx, filename in enumerate(list_files):
+        data = np.loadtxt(filename, unpack=True, skiprows=1, delimiter=',')
+        energies = data[0]
+        eps_r = data[1]
+
+        ax.plot(energies, eps_r, color=cmap(norm(list_q[idx])))
+
+    fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax, label='$|\mathbf{k}|$')
     ax.set_xlabel('Energy [eV]')
     ax.set_ylabel('$\epsilon$')
-    ax.legend()
+    # ax.legend(title="$|\mathbf{k}|$", fontsize=8)
     # ax.set_title('Eps vs Energy')
     fig.tight_layout()
     filename =f"{Path(list_files[0]).with_suffix('.png')}"
     fig.savefig(filename, dpi=600)
-    plt.show()
+    plt.show()
+
+
+if __name__ == '__main__':
+    dirname = Path(sys.argv[1])
+    plot_dielectric_function_vs_energy(dirname)
+
+
diff --git a/src/BZ_MESH/impact_ionization.cpp b/src/BZ_MESH/impact_ionization.cpp
@@ -32,20 +32,35 @@ double ImpactIonization::compute_impact_ionization_rate(int idx_n1, const Vector
     return rate;
 }
 
-double ImpactIonization::compute_direct_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_prime) const {
+/**
+ * @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();
-    vector3 k2       = k1_prime + k2_prime - k1;
-
-    std::size_t index_k2       = get_nearest_k_index(k2);
-    std::complex<double>      Ma = 0.0;
+    complex_d   Ma       = 0.0;
+    complex_d   Mb       = 0.0;
 
     double e_charge       = EmpiricalPseudopotential::Constants::q;
     double eps_zero       = EmpiricalPseudopotential::Constants::eps_zero;
@@ -64,23 +79,35 @@ double ImpactIonization::compute_direct_impact_ionization_matrix_element(int idx
                     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 q_a            = k1_prime + g1_prime - k1 - g1;.0000000000000000000000000000000000000000000000000000000000000000000000
+                    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[index_k2](idx_n2, idx_g2);
+                    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);
-                    std::complex<double> micro_matrix_element_a = std::conj(state_k1_prime) * std::conj(state_k2_prime) * state_k1 * state_k2;
-                    std::complex<double> dielectric_func      = get_dielectric_function(q_a, omega_a);
-                    std::complex<double> pre_factor = e_charge * e_charge / (eps_zero * dielectric_func * q_a.norm() * q_a.norm() * volume);
+                    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 += pre_factor * micro_matrix_element;
+                    Ma += micro_matrix_element_a;
+                    Mb += micro_matrix_element_b;
                 }
             }
         }
     }
+
+    return {Ma, Mb};
+
 }
 
 }  // namespace bz_mesh
diff --git a/src/BZ_MESH/impact_ionization.hpp b/src/BZ_MESH/impact_ionization.hpp
@@ -12,24 +12,29 @@
 #pragma once
 
 #include <Eigen/Dense>
+#include <array>
 
 #include "Material.h"
 #include "bz_mesh.hpp"
 #include "bz_states.hpp"
 
 namespace bz_mesh {
 
+typedef std::complex<double> complex_d;
+
 class ImpactIonization : private BZ_States {
  private:
     double m_impact_ionization_rate = 0.0;
+
  public:
-    double compute_direct_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_prime) const;
+    std::array<complex_d, 2> 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;
 
     double compute_impact_ionization_rate(int idx_n1, const Vector3D<double>& k1);
 };