diff --git a/python/kpoints.py b/python/kpoints.py
index 9fc2677..fa68bd7 100644
--- a/python/kpoints.py
+++ b/python/kpoints.py
@@ -85,6 +85,7 @@ def main():
         
     directions = [int(x) for x in args.directions]
     shift = [float(x) for x in args.shift.split(',')]
+    print(shift)
     irreducible_wedge = args.irreducible_wedge
     output = args.output
     
@@ -103,6 +104,7 @@ def main():
         kx = np.linspace(min, max, Nx)
         ky = np.linspace(min, max, Ny)
         list_kpoints = [[x, y, 0.0] for x in kx for y in ky if np.linalg.norm(np.array([x, y, 0.0]) - center) <= radius]
+        
     elif args.grid == '3D':
         kx = np.linspace(min, max, Nx)
         ky = np.linspace(min, max, Ny)
@@ -113,6 +115,9 @@ def main():
     else:
         raise ValueError('Unknown grid type: {}'.format(args.grid))
     
+    # Apply shift
+    list_kpoints = [np.array(k) + np.array(shift) for k in list_kpoints]
+    
     print(f"Generated {len(list_kpoints)} k-points.")
     print(f"Writing k-points to {output}.")
     np.savetxt(output, list_kpoints, fmt='%.8e')
diff --git a/python/plots/plot_epsilon_map.py b/python/plots/plot_epsilon_map.py
index a102ef9..d2eaeff 100644
--- a/python/plots/plot_epsilon_map.py
+++ b/python/plots/plot_epsilon_map.py
@@ -28,7 +28,7 @@ def get_crystal_dir(list_q):
 
 
 def parse_epsilon_files(dirname):
-    files = glob.glob(dirname + "/*.csv")
+    files = glob.glob(dirname + "/Si*.csv")
     list_qxyz = []
     list_norm_q = []
     data = []
@@ -37,13 +37,15 @@ def parse_epsilon_files(dirname):
     for idx, filename in enumerate(files):
         print(f"\r File n° {idx} / {len(files)}", end="", flush=True)
         stem_name = Path(filename).stem
+        split_name = stem_name.split("_")
+        print(split_name)
         qx, qy, qz = stem_name.split(
             "_")[-3], stem_name.split("_")[-2], stem_name.split("_")[-1]
         norm_q = np.sqrt(float(qx)**2 + float(qy)**2 + float(qz)**2)
         list_norm_q.append(norm_q)
         # print(f"q = (qx, qy, qz) = ({qx}, {qy}, {qz})")
         list_qxyz.append((qx, qy, qz))
-        energies, eps_r = np.loadtxt(
+        energies, eps_r, eps_i = np.loadtxt(
             filename, unpack=True, skiprows=1, delimiter=',')
         # print("Type and shape eps_r: ", type(eps_r), eps_r.shape)
         list_energies = energies
diff --git a/src/BZ_MESH/bz_states.cpp b/src/BZ_MESH/bz_states.cpp
index c232cb4..a59b130 100644
--- a/src/BZ_MESH/bz_states.cpp
+++ b/src/BZ_MESH/bz_states.cpp
@@ -132,7 +132,7 @@ void BZ_States::compute_dielectric_function(const std::vector<double>& list_ener
     std::cout << "Total volume: " << total_volume << std::endl;
 
     double q_squared = m_q_shift.Length() * m_q_shift.Length();
-    double pre_factor = 4.0 * M_PI / q_squared;
+    double pre_factor = 2.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] = 1.0 + pre_factor * dielectric_function_real_at_energies[index_energy] / total_volume;
     }
diff --git a/src/EPP/DielectricFunction.cpp b/src/EPP/DielectricFunction.cpp
index bea5983..d846382 100644
--- a/src/EPP/DielectricFunction.cpp
+++ b/src/EPP/DielectricFunction.cpp
@@ -99,7 +99,6 @@ 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();
@@ -121,7 +120,6 @@ void DielectricFunction::compute_dielectric_function(double eta_smearing) {
                             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 =
@@ -193,7 +191,7 @@ DielectricFunction DielectricFunction::merge_results(DielectricFunction
         for (std::size_t index_energy = 0; index_energy < total_dielectric_function[index_q].size(); ++index_energy) {
             total_dielectric_function[index_q][index_energy] *= renormalization;
             total_dielectric_function[index_q][index_energy] =
-                1.0 + (4.0 * M_PI / q_squared) * total_dielectric_function[index_q][index_energy];
+                1.0 + (2.0 * M_PI / q_squared) * total_dielectric_function[index_q][index_energy];
         }
     }
 
@@ -251,6 +249,7 @@ void DielectricFunction::export_dielectric_function_at_q(const std::string& file
         outname = filename;
     }
     std::ofstream outfile(outname);
+    std::cout << m_energies[0] << " " << m_dielectric_function_real[idx_q][0] << std::endl;
     outfile << "Energy (eV),EpsilonReal,EpsilonImaginary" << std::endl;
     for (std::size_t index_energy = 0; index_energy < m_energies.size(); ++index_energy) {
         outfile << m_energies[index_energy] << "," << m_dielectric_function_real[idx_q][index_energy] << ","