diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp
index f5dc605..01ba331 100644
--- a/c++/pomerol_ed.cpp
+++ b/c++/pomerol_ed.cpp
@@ -127,17 +127,6 @@ namespace pomerol2triqs {
       std::cout << "Number of Blocks is " << states_class->NumberOfBlocks() << std::endl;
     }
 
-    // Save quantum numbers and block size
-    if (verbose && !comm.rank()) { // TODO : enable no output
-      std::string filename("quantum_numbers.dat");  // TODO : specify externally
-      std::ofstream fout(filename);
-      fout << "# block_size quantum_numbers" << std::endl;
-      for (Pomerol::BlockNumber i=0; i<states_class->NumberOfBlocks(); i++)
-        fout << states_class->getBlockSize(i) << "  " << states_class->getQuantumNumbers(i) << std::endl;
-      fout.close();
-      std::cout << "'" << filename << "'" << std::endl;
-    }
-
     if (verbose && !comm.rank()) { std::cout << "\nPomerol: diagonalizing Hamiltonian" << std::endl; }
     // Matrix representation of the Hamiltonian
     matrix_h.reset(new Pomerol::Hamiltonian(index_info, *storage, *states_class));
@@ -147,27 +136,6 @@ namespace pomerol2triqs {
     // Get ground state energy
     if (verbose && !comm.rank()) { std::cout << "\nPomerol: ground state energy is " << matrix_h->getGroundEnergy() + gs_shift << std::endl; }
 
-    // Save all eigenvalues and corresponding quantum numbers
-    if (verbose && !comm.rank()) { // TODO : enable no output
-      // create a list of pairs of eigenvalue and quantum numers to sort
-      std::vector< std::pair<double, Pomerol::QuantumNumbers> > eigen;
-      for (Pomerol::BlockNumber i=0; i<states_class->NumberOfBlocks(); i++){
-          Pomerol::HamiltonianPart H_part = matrix_h->getPart(i);
-          for (Pomerol::InnerQuantumState j=0; j<H_part.getSize(); j++)
-              eigen.push_back( std::make_pair( H_part.getEigenValue(j), H_part.getQuantumNumbers() ) );
-      }
-      // sort eigenvalues in ascending order
-      std::sort(eigen.begin(), eigen.end());
-      // write into a file
-      std::string filename("eigenvalues.dat");  // TODO : specify externally
-      std::ofstream fout(filename);
-      fout << "# eigenvalues quantum_numbers" << std::endl;
-      for (auto e : eigen)
-        fout << e.first << "  " << e.second << std::endl;
-      fout.close();
-      std::cout << "'" << filename << "'" << std::endl;
-    }
-
     // Reset containers, we will compute them later if needed
     rho.release();
     ops_container.release();
@@ -222,6 +190,42 @@ namespace pomerol2triqs {
     diagonalize_main(gs_shift);
   }
 
+  void pomerol_ed::saveQuantumNumbers(const std::string &filename) {
+    if (!states_class) TRIQS_RUNTIME_ERROR << "saveQuantumNumbers: internal error!";
+
+    if (!comm.rank()) {
+      std::ofstream fout(filename);
+      fout << "# block_size quantum_numbers" << std::endl;
+      for (Pomerol::BlockNumber i=0; i<states_class->NumberOfBlocks(); i++)
+      fout << states_class->getBlockSize(i) << "  " << states_class->getQuantumNumbers(i) << std::endl;
+      fout.close();
+      std::cout << "'" << filename << "'" << std::endl;
+    }
+  }
+
+  void pomerol_ed::saveEigenValues(const std::string &filename) {
+    if (!states_class || !matrix_h) TRIQS_RUNTIME_ERROR << "saveEigenValues: internal error!";
+
+    if (!comm.rank()) {
+      // create a list of pairs of eigenvalue and quantum numers to sort
+      std::vector< std::pair<double, Pomerol::QuantumNumbers> > eigen;
+      for (Pomerol::BlockNumber i=0; i<states_class->NumberOfBlocks(); i++){
+          Pomerol::HamiltonianPart H_part = matrix_h->getPart(i);
+          for (Pomerol::InnerQuantumState j=0; j<H_part.getSize(); j++)
+              eigen.push_back( std::make_pair( H_part.getEigenValue(j), H_part.getQuantumNumbers() ) );
+      }
+      // sort eigenvalues in ascending order
+      std::sort(eigen.begin(), eigen.end());
+      // write into a file
+      std::ofstream fout(filename);
+      fout << "# eigenvalues quantum_numbers" << std::endl;
+      for (auto e : eigen)
+        fout << e.first << "  " << e.second << std::endl;
+      fout.close();
+      std::cout << "'" << filename << "'" << std::endl;
+    }
+  }
+
   Pomerol::ParticleIndex pomerol_ed::lookup_pomerol_index(indices_t const &i) const {
     auto it = index_converter.find(i);
     if (it == index_converter.end()) return -1;
diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp
index 25ec8c8..8f3f431 100644
--- a/c++/pomerol_ed.hpp
+++ b/c++/pomerol_ed.hpp
@@ -76,6 +76,12 @@ namespace pomerol2triqs {
     /// Diagonalize Hamiltonian using provided integrals of motion
     void diagonalize(many_body_op_t const &hamiltonian, std::vector<many_body_op_t> const& integrals_of_motion);
 
+    /// Save quantum numbers and block size
+    void saveQuantumNumbers(const std::string &filename);
+
+    /// Save all eigenvalues and corresponding quantum numbers
+    void saveEigenValues(const std::string &filename);
+
     /// Green's function in Matsubara frequencies
     block_gf<imfreq> G_iw(gf_struct_t const &gf_struct, double beta, int n_iw);
 
diff --git a/example/2band.atom.py b/example/2band.atom.py
index 4fd4ba2..7e45dde 100644
--- a/example/2band.atom.py
+++ b/example/2band.atom.py
@@ -66,6 +66,10 @@
 # Diagonalize H
 ed.diagonalize(H)
 
+# save data
+ed.saveQuantumNumbers("quantum_numbers.dat")
+ed.saveEigenValues("eigenvalues.dat")
+
 # 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 84ab3dd..1bb9e91 100644
--- a/python/pomerol2triqs_desc.py
+++ b/python/pomerol2triqs_desc.py
@@ -47,6 +47,12 @@
 c.add_method("""void diagonalize (many_body_op_t hamiltonian, std::vector<many_body_op_t> integrals_of_motion)""",
              doc = """Diagonalize Hamiltonian using provided integrals of motion """)
 
+c.add_method("""void saveQuantumNumbers (std::string filename)""",
+             doc = """Save quantum numbers and block size """)
+
+c.add_method("""void saveEigenValues (std::string filename)""",
+             doc = """Save all eigenvalues and corresponding quantum numbers """)
+
 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 """)