diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp
index 01ba331..67464ca 100644
--- a/c++/pomerol_ed.cpp
+++ b/c++/pomerol_ed.cpp
@@ -258,7 +258,6 @@ namespace pomerol2triqs {
       rho.reset(new Pomerol::DensityMatrix(*states_class, *matrix_h, beta));
       rho->prepare();
       rho->compute();
-      double DensityMatrixCutoff = 1e-10;  // TODO : specify externally
       rho->truncateBlocks(DensityMatrixCutoff, verbose);
     }
   }
diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp
index 8f3f431..60e37a3 100644
--- a/c++/pomerol_ed.hpp
+++ b/c++/pomerol_ed.hpp
@@ -66,6 +66,8 @@ namespace pomerol2triqs {
     // block2_gf<Mesh, tensor_valued<4>> compute_g2(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, block_order_t block_order,
     //                                              g2_blocks_t const &g2_blocks, Filler filler) const;
 
+    double DensityMatrixCutoff = 1e-15;
+
     public:
     /// Create a new solver object
     pomerol_ed(index_converter_t const &index_converter, bool verbose = false);
@@ -82,6 +84,9 @@ namespace pomerol2triqs {
     /// Save all eigenvalues and corresponding quantum numbers
     void saveEigenValues(const std::string &filename);
 
+    /// Set DensityMatrixCutoff (default: 1e-15, 0 = No cutoff). DensityMatrix will be recomputed.
+    void setDensityMatrixCutoff(const double DensityMatrixCutoff);
+
     /// Green's function in Matsubara frequencies
     block_gf<imfreq> G_iw(gf_struct_t const &gf_struct, double beta, int n_iw);
 
@@ -99,4 +104,9 @@ namespace pomerol2triqs {
     // TRIQS_WRAP_ARG_AS_DICT
     // block2_gf<w_l_lp_t, tensor_valued<4>> G2_iw_l_lp(g2_iw_l_lp_params_t const &p);
   };
+
+  inline void pomerol_ed::setDensityMatrixCutoff(const double DensityMatrixCutoff){
+    this->DensityMatrixCutoff = DensityMatrixCutoff;
+    rho.release();
+  }
 }
diff --git a/example/2band.atom.py b/example/2band.atom.py
index 7e45dde..0439dda 100644
--- a/example/2band.atom.py
+++ b/example/2band.atom.py
@@ -70,6 +70,9 @@
 ed.saveQuantumNumbers("quantum_numbers.dat")
 ed.saveEigenValues("eigenvalues.dat")
 
+# set DensityMatrixCutoff
+ed.setDensityMatrixCutoff(1e-10)
+
 # Compute G(i\omega)
 G_iw = ed.G_iw(gf_struct, beta, n_iw)
 
diff --git a/python/pomerol2triqs_desc.py b/python/pomerol2triqs_desc.py
index 1bb9e91..f71b47d 100644
--- a/python/pomerol2triqs_desc.py
+++ b/python/pomerol2triqs_desc.py
@@ -53,6 +53,9 @@
 c.add_method("""void saveEigenValues (std::string filename)""",
              doc = """Save all eigenvalues and corresponding quantum numbers """)
 
+c.add_method("""void setDensityMatrixCutoff (double DensityMatrixCutoff)""",
+             doc = """Set DensityMatrixCutoff (default: 1e-15, 0 = No cutoff). DensityMatrix will be recomputed. """)
+
 c.add_method("""block_gf<imfreq> G_iw (gf_struct_t gf_struct, double beta, int n_iw)""",
              doc = """Green\'s function in Matsubara frequencies """)