diff --git a/src/database/database.cpp b/src/database/database.cpp
index 3b9cf57..6d147f8 100644
--- a/src/database/database.cpp
+++ b/src/database/database.cpp
@@ -1,77 +1,84 @@
 /*-------------------------------------------------------
 
  - Module : database
  - File : database.cpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 
 #include "database.hpp"
 #include <iostream>
 
 #include "selector.hpp"
 #include "canonicalizer.hpp"
 #include "switch_basis.hpp"
 #include "mineral_selector.hpp"
 #include "reader.hpp"
 
 
 namespace specmicp {
 namespace database {
 
 void Database::make_canonical()
 {
     canonicalize_database(data);
 }
 
 
 void Database::change_basis(std::vector<int>& new_basis)
 {
     BasisSwitcher(data).change_basis(new_basis);
 }
 
 
 // this is a proxy to a DatabaseSelector
 void Database::remove_components(const std::vector<int>& id_components_to_remove)
 {
     DatabaseSelector selector(data);
     selector.remove_component(id_components_to_remove);
 }
 
+// this is a proxy to a DatabaseSelector
+void Database::keep_only_components(const std::vector<int>& id_components_to_keep)
+{
+    DatabaseSelector selector(data);
+    selector.keep_only_component(id_components_to_keep);
+}
+
 void Database::minerals_keep_only(const std::vector<std::string>& minerals_to_keep)
 {
     MineralSelector(data).keep_only(minerals_to_keep);
 }
 
 void Database::swap_components(std::map<std::string, std::string>& swap_to_make)
 {
     std::vector<int> id_swap;
     id_swap.reserve(data->nb_component);
     for (int i=0; i<data->nb_component; ++i)
     {
         id_swap.push_back(i);
     }
     for (auto it=swap_to_make.begin(); it!=swap_to_make.end(); ++it)
     {
         if (it->first == "H2O") {
             throw std::invalid_argument("Basis switch : water cannot be swapped from basis");
         }
         const int idc = safe_label_to_id(it->first, data->labels_basis);
         const int idaq = data->nb_component + label_to_id(it->second, data->labels_aqueous);
         id_swap[idc] = idaq;
     }
     change_basis(id_swap);
 }
 
 void Database::parse_database(std::string filepath)
 {
     DataReader reader(filepath);
     reader.parse();
-    data = reader.get_database();
+    set_database(reader.get_database());
 }
 
 
 } // end namespace database
 } // end namespace specmicp
diff --git a/src/database/database.hpp b/src/database/database.hpp
index 7fdb0ca..2d1fab3 100644
--- a/src/database/database.hpp
+++ b/src/database/database.hpp
@@ -1,86 +1,87 @@
 /*-------------------------------------------------------
 
  - Module : database
  - File : database.hpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_DATABASE_DATABASE_HPP
 #define SPECMICP_DATABASE_DATABASE_HPP
 
 //! \file database.hpp Database management
 
 #include "module.hpp"
 
 namespace specmicp {
 //! \namespace database Database management
 namespace database {
 
 //! \brief Database management
 //!
 //! Prepare the database for the computation.
 //! This class is a proxy for one-tasked modules.
 class Database: public DatabaseModule
 {
 public:
     //! \brief Default constructor
     //!
     //! the method parse_database must be called to parse the json database
     Database() {}
 
     //! \brief Initialise the database py parsing 'filepath'
     Database(std::string filepath) {
         parse_database(filepath);
         make_canonical();
     }
 
     //! \brief initialise the database with the raw database
     Database(std::shared_ptr<DataContainer> raw_data):
-        data(raw_data)
+        DatabaseModule(raw_data)
     {}
 
 
     //! \brief Parse the database 'filepath'
     void parse_database(std::string filepath);
 
     //! \brief Transform the database into canonical form
     //!
     //! This process is necessary before any computation
     void make_canonical();
 
     //! \brief Return the database
     //!
     //! Return a smart pointer of a DataCotnainer instance
     //! Note : this is a read write access, be careful
     std::shared_ptr<DataContainer> get_database() {return data;}
 
     //! \brief Change the basis
     //!
     //! @param new_basis list of id of the new basis
     //!
     //! The new basis is a list of id, id = id_component for no swapping
     //! or id = id_aqueous + nb_component for swapping a secondary species
     void change_basis(std::vector<int>& new_basis);
 
     //! \brief Swap some component in the basis
     void swap_components(std::map<std::string, std::string>& swap_to_make);
 
-    //! \brief Remove components non present in the system
+    //! \brief Remove components not present in the system
     void remove_components(const std::vector<int>& id_components_to_remove);
 
+    //! \brief Keep only components in the id_components to keep
+    void keep_only_components(const std::vector<int>& id_components_to_keep);
+
     //! \brief Keep only some minerals at equilibrium
     //!
     //! The effect is to flag all the other minerals as "kinetic"
     void minerals_keep_only(const std::vector<std::string>& minerals_to_keep);
 
-private:
-    std::shared_ptr<DataContainer> data; //!< The data
 };
 
 } // end namespace database
 } // end namespace specmicp
 
 #endif // SPECMICP_DATABASE_DATABASE_HPP
diff --git a/src/database/module.hpp b/src/database/module.hpp
index 43649ce..ea06fca 100644
--- a/src/database/module.hpp
+++ b/src/database/module.hpp
@@ -1,102 +1,105 @@
 /*-------------------------------------------------------
 
  - Module : database
  - File : module
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_DATAABASE_MODULE_HPP
 #define SPECMICP_DATAABASE_MODULE_HPP
 
 //! \file module.hpp base class for a database module
 
 
 #include <memory>
 #include "data_container.hpp"
 
 namespace specmicp {
 namespace database {
 
 
 //! \brief Base class for a database module
 class DatabaseModule
 {
 public:
     //! The type of a smart pointer to a raw database
     using  RawDatabasePtr = std::shared_ptr<DataContainer>;
     //! Index that indicates no known species
     static const int no_species = -1;
 
     //! Create a blank database
     DatabaseModule() {
         reset_database();
     }
     //! Use the database thedata
     DatabaseModule(RawDatabasePtr thedata): data(thedata) {}
 
+    //! \brief set the database
+    void set_database(RawDatabasePtr thedata) {data = thedata;}
+
     //! \brief (Re)set the database - create a blank database
     void reset_database() {data = std::make_shared<DataContainer>();}
 
     //! \brief Return true if the database is canonical
     bool is_database_canonical() {return data->is_canonical;}
 
     //! \brief Return the id of species "label" in list_labels
     //!
     //! return no_species if the species cannot be found in the database
     static int label_to_id(const std::string& label, const std::vector<std::string>& list_labels)
     {
         auto search = std::find(list_labels.begin(), list_labels.end(), label);
         if (search == list_labels.end()) return DatabaseModule::no_species;
         else return (search - list_labels.begin());
     }
     //! \brief Return the id of species "label" in list_labels
     //!
     //! Throw std::invalid_argument if the species cannot be found in the list
     static int safe_label_to_id(const std::string& label, const std::vector<std::string>& list_labels)
     {
         auto search = std::find(list_labels.begin(), list_labels.end(), label);
         if (search == list_labels.end()) throw std::invalid_argument("Unknown species : '"+label+"'.");
         else return (search - list_labels.begin());
     }
 
     //! \brief Return the id of the mineral "label"
     //!
     //! Return no_species if label is not a mineral (at equilibrium)
     int mineral_label_to_id(const std::string& label)
     {
         return label_to_id(label, data->labels_minerals);
     }
     //! \brief Return the id of the mineral "label"
     //!
     //! Return no_species if label is not a mineral (at equilibrium)
     int mineral_kinetic_label_to_id(const std::string& label)
     {
         return label_to_id(label, data->labels_minerals_kinetic);
     }
     //! \brief Return the id of the component "label"
     //!
     //! Return no_species if "label" is not a component
     int component_label_to_id(const std::string& label)
     {
         return label_to_id(label, data->labels_basis);
     }
     //! \brief Return the id of the secondary aqueous species "label"
     //!
     //! Return no_species if "label" is not an aqueous species.
     //! A component (of the basis) is not considered as a secondary aqueous species.
     int aqueous_label_to_it(const std::string& label)
     {
         return label_to_id(label, data->labels_aqueous);
     }
 
 protected:
     RawDatabasePtr data; //! The database
 };
 
 } // end namespace database
 } // end namespace specmicp
 
 #endif // SPECMICP_DATAABASE_MODULE_HPP
diff --git a/src/database/selector.cpp b/src/database/selector.cpp
index 8b236c7..7ac5d84 100644
--- a/src/database/selector.cpp
+++ b/src/database/selector.cpp
@@ -1,181 +1,197 @@
 /*-------------------------------------------------------
 
  - Module : database
  - File : selector.cpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 
 #include "selector.hpp"
 
 namespace specmicp {
 namespace database {
 
 
 
 void DatabaseSelector::analyse_component(const std::vector<int> &id_component_to_remove)
 {
     for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
     {
         m_is_component_to_remove[*it] = true;
     }
     m_nb_component_to_keep = data->nb_component - id_component_to_remove.size();
     int new_i = 0;
     for (int i=0; i<data->nb_component; ++i)
     {
         if (is_component_to_remove(i)) continue;
         data->param_aq.row(new_i) = data->param_aq.row(i);
         ++new_i;
     }
     assert(new_i == nb_component_to_keep());
 }
 
 void DatabaseSelector::select_aqueous(const std::vector<int>& id_component_to_remove)
 {
     // first we select which aqueous species we should remove
     std::vector<bool> aqueous_to_remove(data->nb_aqueous, false);
     for (int j=0; j<data->nb_aqueous; ++j)
     {
         for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
         {
             if (data->nu_aqueous(j, *it) != 0)
             {
                 aqueous_to_remove[j] = true;
                 break;
             }
         }
     }
     // then we remove them, in two steps
     const int nb_aq_to_keep = std::count(aqueous_to_remove.begin(), aqueous_to_remove.end(), false);
     int new_j = 0;
     //      1) first we copy data in the beginning of the arrays
     for (int j=0; j<data->nb_aqueous; ++j)
     {
         if (aqueous_to_remove[j])  continue;
         data->logk_aqueous(new_j) = data->logk_aqueous(j);
         data->param_aq.row(nb_component_to_keep()+new_j) = data->param_aq.row(data->nb_component+j);
         data->labels_aqueous[new_j] = data->labels_aqueous[j];
         int new_i = 0;
         for (int i=0; i<data->nb_component; ++i)
         {
             if (is_component_to_remove(i)) continue;
             data->nu_aqueous(new_j, new_i) = data->nu_aqueous(j, i);
             ++new_i;
         }
         assert(new_i == nb_component_to_keep());
         ++new_j;
     }
     assert(new_j == nb_aq_to_keep);
     //      2) then we resize the arrays
     data->nu_aqueous.conservativeResize(nb_aq_to_keep, nb_component_to_keep());
     data->logk_aqueous.conservativeResize(nb_aq_to_keep);
     data->labels_aqueous.resize(nb_aq_to_keep);
     data->param_aq.conservativeResize(nb_component_to_keep()+nb_aq_to_keep, Eigen::NoChange);
     data->nb_aqueous = nb_aq_to_keep;
     // By doing in this order, we avoid bulk copy of arrays
 }
 
 void DatabaseSelector::select_minerals(const std::vector<int>& id_component_to_remove)
 {
     // first we select which minerals we should remove
     std::vector<bool> minerals_to_remove(data->nb_mineral, false);
     for (int m=0; m<data->nb_mineral; ++m)
     {
         for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
         {
             if (data->nu_mineral(m, *it) != 0)
             {
                 minerals_to_remove[m] = true;
                 break;
             }
         }
     }
     // then we remove them, in two steps
     const int nb_min_to_keep = std::count(minerals_to_remove.begin(), minerals_to_remove.end(), false);
     int new_m = 0;
     //      1) first we copy data in the beginning of the arrays
     for (int m=0; m<data->nb_mineral; ++m)
     {
         if (minerals_to_remove[m])  continue;
         data->logk_mineral(new_m) = data->logk_mineral(m);
         data->labels_minerals[new_m] = data->labels_minerals[m];
         int new_i = 0;
         for (int i=0; i<data->nb_component; ++i)
         {
             if (is_component_to_remove(i)) continue;
             data->nu_mineral(new_m, new_i) = data->nu_mineral(m, i);
             ++new_i;
         }
         assert(new_i == nb_component_to_keep());
         ++new_m;
     }
     assert(new_m == nb_min_to_keep);
     //      2) then we resize the arrays
     data->nu_mineral.conservativeResize(nb_min_to_keep, nb_component_to_keep());
     data->logk_mineral.conservativeResize(nb_min_to_keep);
     data->labels_minerals.resize(nb_min_to_keep);
     data->nb_mineral = nb_min_to_keep;
     // By doing in this order, we avoid bulk copy of arrays
 }
 
 void DatabaseSelector::select_minerals_kinetic(const std::vector<int>& id_component_to_remove)
 {
     // first we select which minerals we should remove
     std::vector<bool> minerals_to_remove(data->nb_mineral_kinetic, false);
     for (int m=0; m<data->nb_mineral_kinetic; ++m)
     {
         for (auto it=id_component_to_remove.begin(); it!=id_component_to_remove.end(); ++it)
         {
             if (data->nu_mineral_kinetic(m, *it) != 0)
             {
                 minerals_to_remove[m] = true;
                 break;
             }
         }
     }
     // then we remove them, in two steps
     const int nb_min_to_keep = std::count(minerals_to_remove.begin(), minerals_to_remove.end(), false);
     int new_m = 0;
     //      1) first we copy data in the beginning of the arrays
     for (int m=0; m<data->nb_mineral_kinetic; ++m)
     {
         if (minerals_to_remove[m])  continue;
         data->logk_mineral_kinetic(new_m) = data->logk_mineral_kinetic(m);
         data->labels_minerals_kinetic[new_m] = data->labels_minerals_kinetic[m];
         int new_i = 0;
         for (int i=0; i<data->nb_component; ++i)
         {
             if (is_component_to_remove(i)) continue;
             data->nu_mineral_kinetic(new_m, new_i) = data->nu_mineral_kinetic(m, i);
             ++new_i;
         }
         assert(new_i == nb_component_to_keep());
         ++new_m;
     }
     assert(new_m == nb_min_to_keep);
     //      2) then we resize the arrays
     data->nu_mineral_kinetic.conservativeResize(nb_min_to_keep, nb_component_to_keep());
     data->logk_mineral_kinetic.conservativeResize(nb_min_to_keep);
     data->labels_minerals_kinetic.resize(nb_min_to_keep);
     data->nb_mineral_kinetic = nb_min_to_keep;
     // By doing in this order, we avoid bulk copy of arrays
 }
 
 void DatabaseSelector::select_components()
 {
     int new_i = 0;
     for (int i=0; i<data->nb_component; ++i)
     {
         if (is_component_to_remove(i)) continue;
         data->labels_basis[new_i] = data->labels_basis[i];
         ++new_i;
     }
     assert(new_i == nb_component_to_keep());
     data->labels_basis.resize(nb_component_to_keep());
     data->nb_component = nb_component_to_keep();
 }
 
+void DatabaseSelector::keep_only_component(const std::vector<int>& id_to_keep)
+{
+    // First build the list of components to remove
+    std::vector<int> id_to_remove;
+    id_to_remove.reserve(data->nb_component - id_to_keep.size());
+    for (int id=0; id<data->nb_component; ++id)
+    {
+        auto search = std::find(id_to_keep.begin(), id_to_keep.end(), id);
+        if (search == id_to_keep.end())
+            id_to_remove.push_back(id);
+    }
+    // Then remove them
+    remove_component(id_to_remove);
+}
+
+
 } // end namespace database
 } // end namespace specmicp
diff --git a/src/database/selector.hpp b/src/database/selector.hpp
index 50d3233..55c6bcd 100644
--- a/src/database/selector.hpp
+++ b/src/database/selector.hpp
@@ -1,71 +1,75 @@
 /*-------------------------------------------------------
 
  - Module : database
  - File : selector.hpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_DATABASE_SELECTOR_HPP
 #define SPECMICP_DATABASE_SELECTOR_HPP
 
 //! \file selector.hpp select species int the database
 
 #include <vector>
 #include "module.hpp"
 
 namespace specmicp {
 namespace database {
 
 //! \brief Select species in the database
 //!
 //! This is a special class with the only function is to remove some components from the database
 class DatabaseSelector: public DatabaseModule
 {
 public:
     DatabaseSelector(RawDatabasePtr thedata):
         DatabaseModule(thedata),
         m_is_component_to_remove(thedata->nb_component, false)
     {}
 
+    //! Remove compoment in the list "id_component_to_remove"
     void remove_component(const std::vector<int>& id_component_to_remove)
     {
         analyse_component(id_component_to_remove);
         select_aqueous(id_component_to_remove);
         select_minerals(id_component_to_remove);
         select_minerals_kinetic(id_component_to_remove);
         select_components();
     }
 
+    //! keep only components in the "id_to_keep" list
+    void keep_only_component(const std::vector<int>& id_to_keep);
+
 private:
     //! \brief Analyse wich component should be removed
     void analyse_component(const std::vector<int>& id_component_to_remove);
 
     //! \brief Select the correct aqueous species
     void select_aqueous(const std::vector<int>& id_component_to_remove);
 
     //! \brief Select the correct minerals
     void select_minerals(const std::vector<int>& id_component_to_remove);
     //! \brief Select the correct minerals governed by kinetics
     void select_minerals_kinetic(const std::vector<int>& id_component_to_remove);
 
     //! \brief Select components
     //!
     //! This should be the last operation
     void select_components();
 
     //! \brief Return true if the component must be removed
     bool is_component_to_remove(int id) {return m_is_component_to_remove[id];}
 
     int nb_component_to_keep() {return m_nb_component_to_keep;}
 
     std::vector<bool> m_is_component_to_remove;
     int m_nb_component_to_keep;
 };
 
 } // end namespace database
 } // end namespace specmicp
 
 #endif // SPECMICP_DATABASE_SELECTOR_HPP