diff --git a/apps/elph.cpp b/apps/elph.cpp
index 4e46ba3..2894544 100644
--- a/apps/elph.cpp
+++ b/apps/elph.cpp
@@ -44,6 +44,8 @@ int main(int argc, char const *argv[])
 
     cmd.parse(argc, argv);
 
+    auto start = std::chrono::high_resolution_clock::now();
+
     EmpiricalPseudopotential::Materials materials;
     const std::string                   file_material_parameters = std::string(CMAKE_SOURCE_DIR) + "/parameter_files/materials.yaml";
     materials.load_material_parameters(file_material_parameters);
@@ -65,6 +67,16 @@ int main(int argc, char const *argv[])
     ElectronPhonon.read_mesh_geometry_from_msh_file(mesh_band_input_file);
     ElectronPhonon.read_mesh_bands_from_msh_file(mesh_band_input_file);
 
-    // ElectronPhonon.compute_electron_phonon_rates_over_mesh();
+    ElectronPhonon.compute_electron_phonon_rates_over_mesh();
+
+    ElectronPhonon.export_rate_values("rates_all.csv");
+
+    ElectronPhonon.compute_plot_electron_phonon_rates_vs_energy_over_mesh(my_options.nrLevels, 10.0, 0.01, "rates_vs_energy.csv");
+
+
+    auto stop = std::chrono::high_resolution_clock::now();
+    std::cout << "Time taken: " << std::chrono::duration_cast<std::chrono::seconds>(stop - start).count() << " seconds" << std::endl;
+
+    return 0;
 
 }
diff --git a/parameter_files/phonon_michaillat.yaml b/parameter_files/phonon_michaillat.yaml
index 22d9d56..ce61e3f 100644
--- a/parameter_files/phonon_michaillat.yaml
+++ b/parameter_files/phonon_michaillat.yaml
@@ -3,7 +3,7 @@
 
 materials:
   - name: Si
-    Radius-WS: 2.122
+    Radius-WS: 2.122e-10
     dispersion:
       longitudinal:
         acoustic:
@@ -29,8 +29,8 @@ materials:
           A: 6.0
           B: -2.0
         optic:
-          A: 2.5e17
-          B: -8.0e16
+          B: -8.0e20
+          A: 2.5e21
         energy-threshold: 2.0
       hole:
         acoustic:
diff --git a/python/plots/plot_density_of_states.py b/python/plots/plot_density_of_states.py
index cd27ede..359698f 100644
--- a/python/plots/plot_density_of_states.py
+++ b/python/plots/plot_density_of_states.py
@@ -13,6 +13,8 @@
 import matplotlib.style
 import matplotlib as mpl
 
+import scienceplots
+
 plt.style.use(['science', 'muted'])
 mpl.rcParams['figure.figsize'] = [3.5, 2.8]
 mpl.rcParams['figure.dpi'] = 300
diff --git a/python/plots/plot_phonon_rate.py b/python/plots/plot_phonon_rate.py
new file mode 100644
index 0000000..cf02808
--- /dev/null
+++ b/python/plots/plot_phonon_rate.py
@@ -0,0 +1,57 @@
+import numpy as np
+import scipy.stats as st
+import pandas as pd
+from pathlib import Path
+import matplotlib.pyplot as plt
+from matplotlib.pyplot import cm, title
+import scipy.stats as st
+from scipy.stats import skewnorm
+import glob
+import os
+from argparse import ArgumentParser
+
+
+import matplotlib.style
+import matplotlib as mpl
+
+import scienceplots
+
+plt.style.use(['science', 'muted', 'scatter', 'grid'])
+
+
+def scatter_plot_rates(filename):
+    fig, ax = plt.subplots()
+    data = np.loadtxt(filename)
+    energy = data[:,0]
+    energy -= np.min(energy)
+    for i in range(1, data.shape[1]):
+        ax.scatter(energy, data[:,i], label=f"Mode {i}", s=1)   
+    ax.set_xlabel("Energy (eV)")
+    ax.set_ylabel("Rate (s$^-1$)")
+    
+    fig.tight_layout()
+    
+    ax.legend()
+    
+def plot_rates(filename):
+    fig, ax = plt.subplots()
+    data = np.loadtxt(filename)
+    energy = data[:,0]
+    energy -= np.min(energy)
+    dos = data[:,1]
+    ax.plot(energy, dos, label="DOS")
+    ax.set_xlabel("Energy (eV)")
+    ax.set_ylabel("DOS")
+    fig.tight_layout()
+    ax.legend()
+    
+    
+    
+    
+if __name__ == "__main__":
+    parser = ArgumentParser()
+    parser.add_argument("-f", "--filename", type=str, required=True, help="Filename to plot")
+    args = parser.parse_args()
+    # scatter_plot_rates(args.filename)
+    plot_rates(args.filename)
+    plt.show()
\ No newline at end of file
diff --git a/src/BZ_MESH/bz_mesh.cpp b/src/BZ_MESH/bz_mesh.cpp
index 1d312e1..14ef31d 100644
--- a/src/BZ_MESH/bz_mesh.cpp
+++ b/src/BZ_MESH/bz_mesh.cpp
@@ -16,10 +16,9 @@
 #include <fstream>
 #include <vector>
 
-#include "rapidcsv.h"
-
 #include "gmsh.h"
 #include "omp.h"
+#include "rapidcsv.h"
 
 #pragma omp declare reduction(merge : std::vector <double> : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()))
 
@@ -52,11 +51,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 lattice_constant = m_material.get_lattice_constant_meter();
     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 = 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,
@@ -99,12 +98,12 @@ void MeshBZ::read_mesh_geometry_from_msh_file(const std::string& filename) {
 /**
  * @brief Get the nearest k index object.
  * Brute force search of the nearest k-point index. :(
- * 
- * @param k 
- * @return std::size_t 
+ *
+ * @param k
+ * @return std::size_t
  */
 std::size_t MeshBZ::get_nearest_k_index(const Vector3D<double>& k) const {
-    vector3 K(k.X, k.Y, k.Z);
+    vector3     K(k.X, k.Y, k.Z);
     std::size_t index_nearest_k = 0;
     double      min_distance    = std::numeric_limits<double>::max();
     for (std::size_t index_k = 0; index_k < m_list_vertices.size(); ++index_k) {
@@ -129,8 +128,6 @@ std::size_t MeshBZ::get_nearest_k_index(const vector3& k) const {
     return index_nearest_k;
 }
 
-
-
 /**
  * @brief Read the energy values for each band at every k-points (vertices) of the mesh.
  *
@@ -206,7 +203,7 @@ double MeshBZ::compute_mesh_volume() const {
         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;
 }
 
@@ -287,4 +284,52 @@ void MeshBZ::export_k_points_to_file(const std::string& filename) const {
     file.close();
 }
 
+bool MeshBZ::is_inside_mesh_geometry(const vector3& k) const {
+    double kx     = k.x();
+    double ky     = k.y();
+    double kz     = k.z();
+    bool   cond_1 = std::abs(kx) <= 1.0 and std::abs(ky) <= 1.0 and std::abs(kz) <= 1.0;
+    bool   cond_2 = std::abs(kx) + std::abs(ky) + std::abs(kz) <= 3.0 / 2.0;
+    return cond_1 and cond_2;
+}
+
+bool MeshBZ::is_inside_mesh_geometry(const Vector3D<double>& k) const {
+    double kx     = k.X;
+    double ky     = k.Y;
+    double kz     = k.Z;
+    bool   cond_1 = std::abs(kx) <= 1.0 and std::abs(ky) <= 1.0 and std::abs(kz) <= 1.0;
+    bool   cond_2 = std::abs(kx) + std::abs(ky) + std::abs(kz) <= 3.0 / 2.0;
+    return cond_1 and cond_2;
+}
+
+vector3 MeshBZ::retrieve_k_inside_mesh_geometry(const vector3& k) const {
+    const vector3 b1 = {-1.0, 1.0, 1.0};
+    const vector3 b2 = {1.0, -1.0, 1.0};
+    const vector3 b3 = {1.0, 1.0, -1.0};
+
+    // std::cout << "k: " << k << std::endl;
+
+    // test k + G
+    std::vector<int> list_n_k = {0, 1, -1, 2, -2, 3, -3, 4, -4};
+    // std::vector<int> list_n_k = {0, 1, -1, 2, -2};
+    for (auto&& n_k_x : list_n_k) {
+        for (auto&& n_k_y : list_n_k) {
+            for (auto&& n_k_z : list_n_k) {
+                vector3 k_plus_G = k + n_k_x * b1 + n_k_y * b2 + n_k_z * b3;
+                if (is_inside_mesh_geometry(k_plus_G)) {
+                    if (k_plus_G.norm() <= 1e-6) {
+                        // std::cout << "k: " << k << std::endl;
+                        // std::cout << " G: " << n_k_x * b1 + n_k_y * b2 + n_k_z * b3 << std::endl;
+                        // std::cout << "k + G: " << k_plus_G << std::endl;
+                    }
+                    return k_plus_G;
+                }
+            }
+        }
+    }
+    std::cout << "No k-point inside the mesh geometry found." << std::endl;
+    throw std::runtime_error("No k-point inside the mesh geometry found.");
+
+}
+
 }  // namespace bz_mesh
\ No newline at end of file
diff --git a/src/BZ_MESH/bz_mesh.hpp b/src/BZ_MESH/bz_mesh.hpp
index 47abcdc..99e0440 100644
--- a/src/BZ_MESH/bz_mesh.hpp
+++ b/src/BZ_MESH/bz_mesh.hpp
@@ -83,6 +83,10 @@ class MeshBZ {
     std::size_t get_number_elements() const { return m_list_tetrahedra.size(); }
     double     get_volume() const { return m_total_volume; }
 
+    bool is_inside_mesh_geometry(const vector3& k) const;
+    bool is_inside_mesh_geometry(const Vector3D<double>& k) const;
+    vector3 retrieve_k_inside_mesh_geometry(const vector3& k) const;
+
     vector3 get_k_at_index(std::size_t index) const { return m_list_vertices[index].get_position(); }
     std::size_t get_nearest_k_index(const Vector3D<double>& k) const;
     std::size_t get_nearest_k_index(const vector3& k) const;
diff --git a/src/BZ_MESH/electron_phonon.cpp b/src/BZ_MESH/electron_phonon.cpp
index 81b6f05..c6cdfc2 100644
--- a/src/BZ_MESH/electron_phonon.cpp
+++ b/src/BZ_MESH/electron_phonon.cpp
@@ -14,7 +14,9 @@
 #include <Eigen/Dense>
 #include <array>
 #include <cmath>
+#include <fstream>
 #include <iostream>
+#include <random>
 #include <tuple>
 
 #include "Constants.hpp"
@@ -24,7 +26,14 @@
 namespace bz_mesh {
 
 double ElectronPhonon::bose_einstein_distribution(double energy, double temperature) {
+    if (energy < 1e-9) {
+        return 0.0;
+    }
     double f = 1.0 / (std::expm1(energy / (EmpiricalPseudopotential::Constants::k_b_eV * temperature)));
+    if (f < 0 || std::isnan(f) || std::isinf(f)) {
+        std::cout << "Energy: " << energy << " Temperature: " << temperature << " f: " << f << std::endl;
+        throw std::runtime_error("Bose-Einstein distribution is negative or NaN");
+    }
     return f;
 }
 
@@ -56,62 +65,122 @@ RateValues ElectronPhonon::compute_electron_phonon_rate(int idx_n1, std::size_t
             Vector3D<double> k2{k_2.x(), k_2.y(), k_2.z()};
 
             double           overlap_integral = electron_overlap_integral(k1, k2);
-            auto             q                = k1 - k2;
+            auto             q                = k2 - k1;
+            auto             initial_q        = q;
             Vector3D<double> q_ph{q.X, q.X, q.X};
-            // for (auto&& mode_phonon : m_phonon_dispersion) {
-            //     double e_ph = mode_phonon.get_phonon_dispersion(q_ph.Length());
-            //     for (auto sign_phonon : {-1.0, 1.0}) {
-            //         double               bose_part    = (bose_einstein_distribution(e_ph, m_temperature) + sign_phonon);
-            //         double               energy_final = energy_n1_k1 + e_ph * sign_phonon;
-            //         double               dos_tetra    = tetra.compute_tetra_dos_energy_band(energy_final, idx_n2);
-            //         DeformationPotential deformation_potential =
-            //             (mode_phonon.m_mode == PhononMode::acoustic) ? m_acoustic_deformation_potential :
-            //             m_optical_deformation_potential;
-            //         double          deformation_potential_value = deformation_potential.get_deformation_potential(q_ph, energy_final);
-            //         PhononMode      phonon_mode                 = mode_phonon.m_mode;
-            //         PhononDirection phonon_direction            = mode_phonon.m_direction;
-            //         PhononEvent     phonon_event                = (sign_phonon == 1.0) ? PhononEvent::absorption : PhononEvent::emission;
-            //         double          rate_value = (EmpiricalPseudopotential::Constants::pi / (m_rho * e_ph)) * deformation_potential_value
-            //         *
-            //                             deformation_potential_value * overlap_integral * overlap_integral * bose_part * dos_tetra;
-            //         RateValue rate(phonon_mode, phonon_direction, phonon_event, rate_value);
-            //         rates_k1_n1.add_rate(rate);
-            //     }
-            // }
+            bool             is_in_bz = is_inside_mesh_geometry(q_ph / m_material.get_fourier_factor());
+
+            // No Umklapp process for now.
+            if (!is_in_bz) {
+                auto new_q = retrieve_k_inside_mesh_geometry(q_ph / m_material.get_fourier_factor()) * m_material.get_fourier_factor();
+                q_ph       = Vector3D<double>{new_q.x(), new_q.y(), new_q.z()};
+                // std::cout << "q: " << q_ph / m_material.get_fourier_factor() << " new_q: " << new_q / m_material.get_fourier_factor() <<
+                // std::endl;
+            }
+            is_in_bz = is_inside_mesh_geometry(q_ph / m_material.get_fourier_factor());
+            if (!is_in_bz) {
+                throw std::runtime_error("Q is not inside the BZ");
+            }
+
             for (auto&& mode_phonon : m_phonon_dispersion) {
                 PhononModeDirection mode_direction = mode_phonon.first;
-                double              e_ph           = mode_phonon.second.get_phonon_dispersion(q_ph.Length());
+                double e_ph = mode_phonon.second.get_phonon_dispersion(q_ph.Length()) * EmpiricalPseudopotential::Constants::h_bar_eV;
+                if (e_ph > 100e-3) {
+                    std::cout << "Energy phonon: " << e_ph << std::endl << std::endl;
+                    throw std::runtime_error("Energy phonon too high");
+                }
+                if (e_ph <= 0.0) {
+                    std::cout << "Energy phonon: " << e_ph << std::endl << std::endl;
+                    std::cout << "Q: " << q_ph << std::endl << std::endl;
+                    std::cout << "Q: " << q_ph.Length() / m_material.get_fourier_factor() << std::endl << std::endl;
+                    std::cout << "Energy (n1, k1): " << energy_n1_k1 << std::endl << std::endl;
+                    bool is_in_bz = is_inside_mesh_geometry(q_ph / m_material.get_fourier_factor());
+                    std::cout << "Is in BZ: " << is_in_bz << std::endl << std::endl;
+                    // throw std::runtime_error("Energy phonon negative");
+                    continue;
+                }
+                // std::cout << "Energy phonon: " << e_ph << std::endl;
                 for (auto sign_phonon : {-1.0, 1.0}) {
-                    double               bose_part    = (bose_einstein_distribution(e_ph, m_temperature) + sign_phonon);
-                    double               energy_final = energy_n1_k1 + e_ph * sign_phonon;
-                    double               dos_tetra    = tetra.compute_tetra_dos_energy_band(energy_final, idx_n2);
-                    DeformationPotential deformation_potential =
-                        (mode_direction.first == PhononMode::acoustic) ? m_acoustic_deformation_potential : m_optical_deformation_potential;
-                    double          deformation_potential_value = deformation_potential.get_deformation_potential(q_ph, energy_final);
-                    PhononMode      phonon_mode                 = mode_direction.first;
-                    PhononDirection phonon_direction            = mode_direction.second;
-                    PhononEvent     phonon_event                = (sign_phonon == 1.0) ? PhononEvent::absorption : PhononEvent::emission;
-                    double          rate_value = (EmpiricalPseudopotential::Constants::pi / (m_rho * e_ph)) * deformation_potential_value *
+                    double bose_part    = (bose_einstein_distribution(e_ph, m_temperature) + 0.5 + 0.5 * std::pow(-1.0, sign_phonon));
+                    double energy_final = energy_n1_k1 + e_ph * sign_phonon;
+                    // std::cout << "Energy (n1, k1): " << energy_n1_k1 << " Energy phonon: " << e_ph << " Energy final: " << energy_final
+                    //           << std::endl;
+                    if (!tetra.is_energy_inside_band(energy_final, idx_n2)) {
+                        continue;
+                    }
+                    // std::cout << "Volume: " << this->get_volume() << std::endl;
+                    double dos_tetra = tetra.compute_tetra_dos_energy_band(energy_final, idx_n2);
+                    // std::cout << "Energy final: " << energy_final << " dos tetra: " << dos_tetra
+                    //           << std::endl;
+                    DeformationPotential deformation_potential       = (mode_direction.first == PhononMode::acoustic)
+                                                                           ? m_acoustic_deformation_potential_e
+                                                                           : m_optical_deformation_potential_e;
+                    double               deformation_potential_value = deformation_potential.get_deformation_potential(q_ph, energy_final);
+                    PhononMode           phonon_mode                 = mode_direction.first;
+                    PhononDirection      phonon_direction            = mode_direction.second;
+                    PhononEvent          phonon_event = (sign_phonon == 1.0) ? PhononEvent::absorption : PhononEvent::emission;
+                    double rate_value = (EmpiricalPseudopotential::Constants::pi / (m_rho * e_ph)) * deformation_potential_value *
                                         deformation_potential_value * overlap_integral * overlap_integral * bose_part * dos_tetra;
+                    // std::cout << this->get_volume() << std::endl;
+                    rate_value /= this->get_volume();
+                    rate_value *= EmpiricalPseudopotential::Constants::q;
+
+                    if (rate_value < 0.0 || rate_value > 1e50 || std::isnan(rate_value) || std::isinf(rate_value)) {
+                        std::cout << "Rate value: " << rate_value << std::endl;
+                        std::cout << "Overlap integral: " << overlap_integral << std::endl;
+                        std::cout << "Deformation potential: " << deformation_potential_value << std::endl;
+                        std::cout << "Bose part: " << bose_part << std::endl;
+                        std::cout << "DOS tetra: " << dos_tetra << std::endl;
+                        std::cout << "Energy final: " << energy_final << std::endl;
+                        std::cout << "q: " << q_ph << std::endl;
+                        std::cout << "initial q: " << initial_q / m_material.get_fourier_factor() << std::endl;
+                        std::cout << "Energy phonon: " << e_ph << std::endl;
+                        std::cout << "Energy (n1, k1): " << energy_n1_k1 << std::endl;
+                        std::cout << "--------------------------------------------------------------------------------" << std::endl;
+                    }
+                    // std::cout << "Rate value: " << (deformation_potential_value) << std::endl;
                     RateValue rate(phonon_mode, phonon_direction, phonon_event, rate_value);
                     rates_k1_n1.add_rate(rate);
                 }
             }
         }
     }
-    rates_k1_n1.print_rates();
+    // rates_k1_n1.print_rates();
+    // std::cout << "*************************************************************************************************************" <<
+    // std::endl;
     return rates_k1_n1;
 }
 
 void ElectronPhonon::compute_electron_phonon_rates_over_mesh() {
+    auto indices_conduction_bands = m_indices_conduction_bands;
+    auto min_idx_conduction_band  = *std::min_element(indices_conduction_bands.begin(), indices_conduction_bands.end());
+    auto max_idx_conduction_band  = *std::max_element(indices_conduction_bands.begin(), indices_conduction_bands.end());
+    std::cout << "Min index conduction band: " << min_idx_conduction_band << std::endl;
+    std::ofstream file("rates.txt");
+
+    std::random_device                     rd;
+    std::mt19937                           gen(rd());
+    std::uniform_real_distribution<double> dis(0.0, 1.0);
+    double                                 p_compute_rate = 0.01;
+
+#pragma omp parallel for schedule(dynamic) num_threads(8)
     for (std::size_t idx_k1 = 0; idx_k1 < m_list_vertices.size(); ++idx_k1) {
+        double r = dis(gen);
+
         std::cout << "Computing rates for k-point " << idx_k1 << std::endl;
-        for (std::size_t idx_n1 = 0; idx_n1 < m_indices_valence_bands.size(); ++idx_n1) {
-            auto rate = compute_electron_phonon_rate(idx_n1, idx_k1);
-            m_list_vertices[idx_k1].add_electron_phonon_rates(rate.to_array());
+        for (std::size_t idx_n1 = min_idx_conduction_band; idx_n1 < max_idx_conduction_band; ++idx_n1) {
+            if (r > p_compute_rate) {
+                std::array<double, 8> array = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+                m_list_vertices[idx_k1].add_electron_phonon_rates(array);
+                continue;
+            } else {
+                auto rate  = compute_electron_phonon_rate(idx_n1, idx_k1);
+                auto array = rate.to_array();
+                m_list_vertices[idx_k1].add_electron_phonon_rates(array);
+            }
         }
     }
-
+    file.close();
     // Add rate to the mesh
 }
 
@@ -119,19 +188,53 @@ void ElectronPhonon::compute_plot_electron_phonon_rates_vs_energy_over_mesh(int
                                                                             double             max_energy,
                                                                             double             energy_step,
                                                                             const std::string& filename) {
-    std::vector<double>              energies;
-    std::vector<std::vector<double>> rates;
+    std::ofstream file(filename);
+    std::vector<double> energies;
+    std::vector<double> list_dos;
+    std::vector<std::array<double, 8>> list_mean_rates;
+    double total_dos = 0.0;
     for (double energy = 0.0; energy < max_energy; energy += energy_step) {
-        energies.push_back(energy);
-        for (int idx_band = 0; idx_band < nb_bands; ++idx_band) {
-            double rate = 0.0;
-            for (auto&& vertex : m_list_vertices) {
-                for (std::size_t idx_band = 0; idx_band < vertex.get_number_bands(); ++idx_band) {
-                    rate += vertex.get_energy_at_band(idx_band);
+        std::cout << "\rEnergy: " << energy << " Max energy: " << max_energy << std::flush;
+        double dos = 0.0;
+        std::array<double, 8> mean_rates;
+        std::fill(mean_rates.begin(), mean_rates.end(), 0.0);
+        for (int idx_band = 4; idx_band < nb_bands; ++idx_band) {
+            for (auto&& tetra : m_list_tetrahedra) {
+                double dos_tetra = tetra.compute_tetra_dos_energy_band(energy, idx_band);
+                dos += dos_tetra;
+                std::array<double, 8> mean_rates = tetra.get_mean_electron_phonon_rates(idx_band);
+                for (std::size_t idx_rate = 0; idx_rate < mean_rates.size(); ++idx_rate) {
+                    mean_rates[idx_rate] += mean_rates[idx_rate] * dos_tetra;
                 }
             }
-            // rates[idx_band].push_back(rate);
         }
+        energies.push_back(energy);
+        list_mean_rates.push_back(mean_rates);
+        list_dos.push_back(dos);
+        total_dos += dos;
+    }
+    for (std::size_t idx_energy = 0; idx_energy < energies.size(); ++idx_energy) {
+        file << energies[idx_energy] << " " << list_dos[idx_energy] << " ";
+        for (auto&& rate : list_mean_rates[idx_energy]) {
+            file << rate << " ";
+        }
+        file << std::endl;
+    }
+    file.close();
+}
+
+void ElectronPhonon::plot_phonon_dispersion(const std::string& filename) const {
+    std::ofstream file(filename);
+    for (auto&& vtx : m_list_vertices) {
+        auto k = vtx.get_position();
+        file << k.x() << " " << k.y() << " " << k.z() << " ";
+        for (auto&& mode_direction : m_phonon_dispersion) {
+            auto q = k;
+            q /= m_material.get_fourier_factor();
+            auto e_ph = mode_direction.second.get_phonon_dispersion(q.norm()) * EmpiricalPseudopotential::Constants::h_bar_eV;
+            file << e_ph << " ";
+        }
+        file << std::endl;
     }
 }
 
@@ -152,8 +255,6 @@ void ElectronPhonon::load_phonon_parameters(const std::string& filename) {
 
     auto same_material = [my_material](const YAML::Node& node) { return node["name"].as<std::string>() == my_material; };
     auto it_material   = std::find_if(list_materials.begin(), list_materials.end(), same_material);
-    
-
 
     if (it_material == list_materials.end()) {
         throw std::runtime_error("Material " + my_material + " not found in file " + filename);
@@ -193,13 +294,31 @@ void ElectronPhonon::load_phonon_parameters(const std::string& filename) {
 
             PhononMode           mode = (wave == "acoustic") ? PhononMode::acoustic : PhononMode::optical;
             DeformationPotential deformationPotential(mode, A, B);
-            if (mode == PhononMode::acoustic) {
-                m_acoustic_deformation_potential = deformationPotential;
-            } else {
-                m_optical_deformation_potential = deformationPotential;
+            if (carrierType == "electron") {
+                if (mode == PhononMode::acoustic) {
+                    m_acoustic_deformation_potential_e = deformationPotential;
+                } else {
+                    m_optical_deformation_potential_e = deformationPotential;
+                }
+            }
+        }
+    }
+}
+
+void ElectronPhonon::export_rate_values(const std::string& filename) const {
+    std::ofstream file(filename);
+    for (auto&& vertex : m_list_vertices) {
+        std::vector<std::array<double, 8>> all_rates = vertex.get_electron_phonon_rates_all_bands();
+        for (std::size_t idx_band = 0; idx_band < all_rates.size(); ++idx_band) {
+            double energy = vertex.get_energy_at_band(idx_band + 4);
+            file << energy << " ";
+            for (auto&& rate : all_rates[idx_band]) {
+                file << rate << " ";
             }
+            file << std::endl;
         }
     }
+    file.close();
 }
 
 }  // namespace bz_mesh
\ No newline at end of file
diff --git a/src/BZ_MESH/electron_phonon.hpp b/src/BZ_MESH/electron_phonon.hpp
index a3a5d64..9d199c7 100644
--- a/src/BZ_MESH/electron_phonon.hpp
+++ b/src/BZ_MESH/electron_phonon.hpp
@@ -18,11 +18,10 @@
 #pragma once
 
 #include <Eigen/Dense>
-
 #include <array>
 #include <cmath>
-#include <map>
 #include <iostream>
+#include <map>
 #include <vector>
 
 #include "bz_states.hpp"
@@ -58,6 +57,13 @@ struct PhononDispersion {
 
     double get_phonon_dispersion(double q) const {
         double e_ph = m_omega + m_vs * q + m_c * q * q;
+        if (e_ph < 0.0) {
+            std::cout << "ERROR: negative phonon energy: " << e_ph << std::endl;
+            std::cout << "Q: " << q << std::endl;
+            std::cout << "OMEGA: " << m_omega << std::endl;
+            std::cout << "VS: " << m_vs << std::endl;
+            std::cout << "C: " << m_c << std::endl;
+        }
         return e_ph;
     }
 };
@@ -72,7 +78,7 @@ struct DeformationPotential {
     DeformationPotential(PhononMode type, double A, double B) : m_mode(type), m_A(A), m_B(B) {}
 
     double get_deformation_potential(Vector3D<double> q, double energy) const {
-        double clamp_energy = std::max(energy, m_energy_threshold);
+        double clamp_energy = std::min(energy, m_energy_threshold);
         if (m_mode == PhononMode::acoustic) {
             return std::sqrt(m_A + clamp_energy * m_B) * q.Length();
         } else {
@@ -199,12 +205,12 @@ typedef std::pair<PhononMode, PhononDirection> PhononModeDirection;
 
 class ElectronPhonon : public BZ_States {
  private:
-    double               m_temperature        = 0.0;
-    double               m_rho                = 0.0;
+    double               m_temperature        = 300.0;
+    double               m_rho                = 2.329e3;
     double               m_radii_wigner_seitz = 0.0;
     HoleOverlapIntParams m_hole_overlap_int_params;
-    DeformationPotential m_acoustic_deformation_potential;
-    DeformationPotential m_optical_deformation_potential;
+    DeformationPotential m_acoustic_deformation_potential_e;
+    DeformationPotential m_optical_deformation_potential_e;
 
     std::map<PhononModeDirection, PhononDispersion> m_phonon_dispersion;
 
@@ -212,6 +218,7 @@ class ElectronPhonon : public BZ_States {
     explicit ElectronPhonon(const EmpiricalPseudopotential::Material& material) : BZ_States(material) {}
 
     void load_phonon_parameters(const std::string& filename);
+    void plot_phonon_dispersion(const std::string& filename) const;
 
     double bose_einstein_distribution(double energy, double temperature);
     double electron_overlap_integral(const Vector3D<double>& k1, const Vector3D<double>& k2);
@@ -224,6 +231,8 @@ class ElectronPhonon : public BZ_States {
 
     void compute_electron_phonon_rates_over_mesh();
 
+    void export_rate_values(const std::string& filename) const;
+
     void compute_plot_electron_phonon_rates_vs_energy_over_mesh(int                nb_bands,
                                                                 double             max_energy,
                                                                 double             energy_step,
diff --git a/src/BZ_MESH/mesh_tetra.cpp b/src/BZ_MESH/mesh_tetra.cpp
index eecd908..089e972 100644
--- a/src/BZ_MESH/mesh_tetra.cpp
+++ b/src/BZ_MESH/mesh_tetra.cpp
@@ -59,7 +59,13 @@ void Tetra::compute_gradient_energy_at_bands() {
         const double                eps_14               = (e_0 - energies_at_vertices[indices_sort[3]]);
         const double                gradient_energy      = sqrt((eps_12 * eps_12 + eps_13 * eps_13 + eps_14 * eps_14));
         if (gradient_energy == 0) {
-            std::cout << "Warning: gradient energy is zero for tetra " << m_index << " at band " << band_index << std::endl;
+            std::cout << "Gradient energy is zero for tetra " << m_index << " at band " << band_index << std::endl;
+            std::cout << "Energies: " << energies_at_vertices[0] << " " << energies_at_vertices[1] << " " << energies_at_vertices[2] << " "
+                      << energies_at_vertices[3] << std::endl;
+            std::cout << m_list_vertices[0]->get_position() << "\n"
+                      << m_list_vertices[1]->get_position() << "\n"
+                      << m_list_vertices[2]->get_position() << "\n"
+                      << m_list_vertices[3]->get_position() << std::endl;
         }
         m_gradient_energy_per_band.push_back(gradient_energy);
     }
@@ -337,4 +343,18 @@ double Tetra::compute_tetra_dos_energy_band(double energy, std::size_t band_inde
            this->compute_tetra_iso_surface_energy_band(energy, band_index);
 }
 
+bool Tetra::is_energy_inside_band(double energy, std::size_t index_band) const {
+    return (energy >= m_min_energy_per_band[index_band] && energy <= m_max_energy_per_band[index_band]);
+}
+
+std::array<double, 8> Tetra::get_mean_electron_phonon_rates(int band_index) const {
+    std::array<double, 8> mean_rates;
+    std::fill(mean_rates.begin(), mean_rates.end(), 0.0);
+    for (std::size_t i = 0; i < 4; i++) {
+        const std::array<double, 8>& rates = m_list_vertices[i]->get_electron_phonon_rates(band_index);
+        std::transform(mean_rates.begin(), mean_rates.end(), rates.begin(), mean_rates.begin(), std::plus<double>());
+    }
+    return mean_rates;
+}
+
 }  // namespace bz_mesh
\ No newline at end of file
diff --git a/src/BZ_MESH/mesh_tetra.hpp b/src/BZ_MESH/mesh_tetra.hpp
index b65854c..914a97a 100644
--- a/src/BZ_MESH/mesh_tetra.hpp
+++ b/src/BZ_MESH/mesh_tetra.hpp
@@ -109,10 +109,13 @@ class Tetra {
     vector3               compute_euclidean_coordinates_with_indices(const std::array<double, 4>& barycentric_coordinates,
                                                                      const std::array<int, 4>&    indices_vertex) const;
 
+    bool is_energy_inside_band(double energy, std::size_t index_band) const;
     std::array<int, 4>   get_index_vertices_with_sorted_energy_at_band(std::size_t index_band) const;
     std::vector<vector3> compute_band_iso_energy_surface(double iso_energy, std::size_t band_index) const;
     double               compute_tetra_iso_surface_energy_band(double energy, std::size_t band_index) const;
     double               compute_tetra_dos_energy_band(double energy, std::size_t band_index) const;
+    
+    std::array<double, 8> get_mean_electron_phonon_rates(int band_index) const;
 
     
 
diff --git a/src/BZ_MESH/mesh_vertex.hpp b/src/BZ_MESH/mesh_vertex.hpp
index 8681015..8cbcc86 100644
--- a/src/BZ_MESH/mesh_vertex.hpp
+++ b/src/BZ_MESH/mesh_vertex.hpp
@@ -49,7 +49,7 @@ class Vertex {
      * [ALO, ALA, ATO, ATA, ELO, ELA, ETO, ETA] (ALO: absorption longitudinal optical, ALA: absorption longitudinal acoustic, ...)
      * 
      */
-    std::vector<std::array<double, 8>> m_electron_phonon_rates;
+    std::vector<std::array<double, 8>> m_electron_phonon_rates = {};
 
  public:
     /**
@@ -150,6 +150,8 @@ class Vertex {
      */
     const std::array<double, 8>& get_electron_phonon_rates(std::size_t band_index) const { return m_electron_phonon_rates[band_index]; }
 
+    const std::vector<std::array<double, 8>>& get_electron_phonon_rates_all_bands() const { return m_electron_phonon_rates; }
+
 };
 
 }  // namespace bz_mesh
\ No newline at end of file
diff --git a/src/EPP/Constants.hpp b/src/EPP/Constants.hpp
index 0cdce8b..b319f2c 100644
--- a/src/EPP/Constants.hpp
+++ b/src/EPP/Constants.hpp
@@ -20,6 +20,7 @@ namespace Constants {
 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 h_bar_eV       = 6.582119514e-16;
 constexpr double m0             = 9.10938356e-31;
 constexpr double q              = 1.602176634e-19;
 constexpr double eps_zero       = 8.854187e-12;
diff --git a/src/EPP/Material.h b/src/EPP/Material.h
index e9f1cdf..3870851 100644
--- a/src/EPP/Material.h
+++ b/src/EPP/Material.h
@@ -110,6 +110,11 @@ class Material {
         return angstrom_to_m * m_lattice_constant;
     }
 
+    double get_fourier_factor() const { 
+        constexpr double angstrom_to_m = 1.0e-10;
+        return 2.0 * M_PI / (angstrom_to_m * m_lattice_constant);
+    }
+
     /**
      * @brief Get the atomic volume of the material in angstrom^3.
      *