diff --git a/src/specmicp/problem_solver/reactant_box.cpp b/src/specmicp/problem_solver/reactant_box.cpp
index 7d39052..adf9655 100644
--- a/src/specmicp/problem_solver/reactant_box.cpp
+++ b/src/specmicp/problem_solver/reactant_box.cpp
@@ -1,576 +1,595 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 
 #include "reactant_box.hpp"
 
 #include "specmicp/adimensional/adimensional_system_structs.hpp"
 
 #include "specmicp_database/data_container.hpp"
 #include "specmicp_database/database.hpp"
 #include "specmicp_database/unknown_class.hpp"
 
 #include "specmicp_common/physics/laws.hpp"
 
 #include <vector>
 #include <unordered_map>
 #include <stdexcept>
 
 #include <iostream>
 
 namespace specmicp {
 
 
 struct SPECMICP_DLL_LOCAL fug_container {
     std::string gas;
     std::string component;
     scalar_t fugacity;
 
     fug_container(
             const std::string& agas,
             const std::string& acomponent,
             scalar_t t_fugacity
             ):
         gas(agas), component(acomponent), fugacity(t_fugacity)
     {}
 };
 
 struct ReactantBox::ReactantBoxImpl {
     scalar_t mass_solution {-1};
     std::vector<std::string> keep_components {};
     std::unordered_map<std::string, scalar_t> aqueous_species {};
     std::unordered_map<std::string, scalar_t> solid_phase {};
 
     // constraints
     std::string charge_keeper {""};
     scalar_t inert_vol_frac {0.0};
     WaterEquationType water_eq_type {WaterEquationType::MassConservation};
     water_partial_pressure_f h2o_partial_pressure_model;
     std::vector<fug_container> fix_fugacity;
     std::vector<std::pair<std::string, scalar_t>> fix_activity;
     std::vector<std::pair<std::string, scalar_t>> fix_molality;
 
     // database lookup
     database::SpeciesTypeFinder db_search;
 
     ReactantBoxImpl(RawDatabasePtr thedb):
         db_search(thedb)
     {}
 
     std::vector<index_t> get_components_to_keep() const;
 
     void dissolve_aqueous_species(
             Vector& total_concentration
             ) const;
     void dissolve_solid_phases(
             Vector& total_concentration
             ) const;
 
     void raise_unknown_unit(const std::string& error_msg) {
         throw std::invalid_argument(error_msg);
     }
     void assert_mass_solution() {
         if (mass_solution == -1) {
             throw std::logic_error("Mass of solution is not set !");
         }
     }
 
     scalar_t factor_total_concentration(const units::UnitsSet& units_set) const;
 
 };
 
 ReactantBox::ReactantBox(
     std::shared_ptr<database::DataContainer> raw_data,
     const units::UnitsSet& units_set
 ):
     units::UnitBaseClass(units_set),
     database::DatabaseHolder(raw_data),
     m_impl(utils::make_pimpl<ReactantBoxImpl>(raw_data))
 {
 
 }
 
 
 ReactantBox::ReactantBox(ReactantBox&& other):
     units::UnitBaseClass(other.get_units()),
     database::DatabaseHolder(other.m_data),
     m_impl(std::move(other.m_impl))
 {
 
 }
 
 
 ReactantBox& ReactantBox::operator=(ReactantBox&& other)
 {
     m_impl = std::move(other.m_impl);
     return *this;
 }
 
 ReactantBox::~ReactantBox() = default;
 
 void ReactantBox::set_solution(scalar_t value, std::string str_unit)
 {
     set_solution(value, units::parse_amount_unit(str_unit));
 }
 void ReactantBox::set_solution(
         scalar_t                   value,
         const units::AmountUnit&   unit
         )
 {
     scalar_t factor {1.0};
     if (unit.type == units::AmountUnitType::Mass) {
         factor = unit.factor_si;
     } else if (unit.type == units::AmountUnitType::Volume) {
         factor = unit.factor_si*laws::density_water(units::SI_units);
     } else if (unit.type == units::AmountUnitType::VolumeFraction) {
         factor = unit.factor_si*laws::density_water(units::SI_units);
     } else if (unit.type == units::AmountUnitType::MassConcentration) {
         factor = unit.factor_si;
     } else {
         m_impl->raise_unknown_unit(
                     "Only mass. volume, volume fraction,"
                     " or mass concentration units "
                     "are accepted for the amount of solution"
                     );
     }
     m_impl->mass_solution = factor*value;
 
 }
 
 void ReactantBox::add_aqueous_species(
         std::string name,
         scalar_t    value,
         std::string str_unit
         )
 {
     add_aqueous_species(name, value, units::parse_amount_unit(str_unit));
 }
 
 void ReactantBox::add_aqueous_species(
         std::string                name,
         scalar_t                   value,
         const units::AmountUnit&   unit
         )
 {
     scalar_t factor = 1.0;
     if (unit.type == units::AmountUnitType::Molality)
     {
         factor = unit.factor_si;
     }
     else if (unit.type == units::AmountUnitType::MoleConcentration)
     {
         factor = unit.factor_si/laws::density_water(units::SI_units);
     }
     else if (unit.type == units::AmountUnitType::MassConcentration)
     {
         const scalar_t molar_mass = m_impl->db_search.molar_mass_aqueous(name);
         factor = unit.factor_si / molar_mass
                 / laws::density_water(units::SI_units);
     }
     else if (unit.type == units::AmountUnitType::Mass)
     {
         m_impl->assert_mass_solution();
         const scalar_t molar_mass =  m_impl->db_search.molar_mass_aqueous(name);
         factor = unit.factor_si / molar_mass / m_impl->mass_solution;
     }
     else if (unit.type == units::AmountUnitType::NumberOfMoles)
     {
         m_impl->assert_mass_solution();
         factor = unit.factor_si / m_impl->mass_solution;
     }
     else
     {
         m_impl->raise_unknown_unit("Unit for aqueous species " + name + "is invalid."
                                    " Only accepted units for aqueous species are"
                                    " molality, mole and mass concettration,"
                                    " mass and number of moles.");
     }
     const scalar_t molality = factor*value;
     m_impl->aqueous_species[name] = molality;
 }
 
 void ReactantBox::add_solid_phase(
         std::string name,
         scalar_t    value,
         std::string str_unit
         )
 {
     return add_solid_phase(name, value, units::parse_amount_unit(str_unit));
 }
 
 void ReactantBox::add_solid_phase(
         std::string                name,
         scalar_t                   value,
         const units::AmountUnit&   unit
         )
 {
 
     scalar_t factor = 1.0;
 
     // convert to mole concentration is SI
     if (unit.type == units::AmountUnitType::MoleConcentration or
             unit.type == units::AmountUnitType::NumberOfMoles)
     {
         factor = unit.factor_si;
     }
     else if (unit.type == units::AmountUnitType::MassConcentration or
              unit.type == units::AmountUnitType::Mass
              )
     {
         const scalar_t molar_mass = m_impl->db_search.molar_mass_solid_phase(name);
         factor = unit.factor_si / molar_mass;
     }
     else if (unit.type == units::AmountUnitType::Volume)
     {
         const scalar_t molar_volume = m_impl->db_search.molar_volume_solid_phase(name);
         factor = unit.factor_si / molar_volume;
     }
     else if (unit.type == units::AmountUnitType::VolumeFraction)
     {
         const scalar_t molar_volume = m_impl->db_search.molar_volume_solid_phase(name);
         factor = unit.factor_si / molar_volume;
 
     }
     else
     {
         m_impl->raise_unknown_unit("Unit for solid phase " + name + " is invalid.");
     }
     m_impl->solid_phase[name] = factor*value;
 }
 
 void ReactantBox::add_component(const std::string& name)
 {
     m_impl->keep_components.push_back(name);
 }
 
 Vector ReactantBox::get_total_concentration(bool modify_db) const
 {
     std::shared_ptr<database::DataContainer> raw_data = get_database();
 
     if (modify_db) {
         // remove extra components
         const auto id_comp_to_keep = m_impl->get_components_to_keep();
         database::Database db_manager(raw_data);
         db_manager.keep_only_components(id_comp_to_keep);
     }
 
 
     Vector total_concentration = Vector::Zero(raw_data->nb_component());
     total_concentration[database::water_id] =
             m_impl->mass_solution/raw_data->molar_mass_basis(database::water_id, units::SI_units);
     m_impl->dissolve_aqueous_species(total_concentration);
     m_impl->dissolve_solid_phases(total_concentration);
 
     total_concentration *= m_impl->factor_total_concentration(get_units());
 
     return total_concentration;
 }
 
 void ReactantBox::ReactantBoxImpl::dissolve_aqueous_species(
         Vector& total_concentration
         ) const
 {
     auto raw_data = db_search.get_database();
     for (auto& it: aqueous_species) {
         database::GenericAqueousSpecies aq =  db_search.find_aqueous_species(it.first);
         const scalar_t conc = mass_solution*it.second;
         if (aq.type == database::AqueousSpeciesClass::Component)
         {
             total_concentration[aq.id] += conc;
         }
         else if (aq.type == database::AqueousSpeciesClass::Aqueous)
         {
             for (auto idc: raw_data->range_component())
             {
                 const auto nu = raw_data->nu_aqueous(aq.id, idc);
                 if (nu != 0)
                 {
                     total_concentration[idc] += nu*conc;
                 }
             }
         }
         else if (aq.type == database::AqueousSpeciesClass::Compound)
         {
             for (auto idc: raw_data->range_component())
             {
                 const auto nu = raw_data->nu_compound(aq.id, idc);
                 if (nu != 0)
                 {
                     total_concentration[idc] += nu*conc;
                 }
             }
         }
         else
         {
             throw std::invalid_argument("'"+it.first+"' is not a valid aqueous species.");
         }
     }
 }
 
 
 void ReactantBox::ReactantBoxImpl::dissolve_solid_phases(
         Vector& total_concentration
         ) const
 {
     auto raw_data = db_search.get_database();
     for (auto& it: solid_phase) {
         database::GenericSolidPhase sp = db_search.find_solid_phase(it.first);
         if (sp.type == database::SolidPhaseClass::EquilibriumMineral)
         {
             for (auto idc: raw_data->range_component()) {
                 const auto nu = raw_data->nu_mineral(sp.id, idc);
                 if (nu != 0) {
                     total_concentration[idc] += nu*it.second;
                 }
             }
         }
         else if (sp.type == database::SolidPhaseClass::MineralKinetics)
         {
             for (auto idc: raw_data->range_component()) {
                 const auto nu = raw_data->nu_mineral_kinetic(sp.id, idc);
                 if (nu != 0) {
                     total_concentration[idc] += nu*it.second;
                 }
             }
         }
         else
         {
             throw std::invalid_argument("'"+it.first+"' is not a valid solid phase.");
         }
     }
 }
 
 std::vector<index_t> ReactantBox::ReactantBoxImpl::get_components_to_keep(
         ) const
 {
     auto raw_data = db_search.get_database();
 
     std::vector<bool> to_keep(raw_data->nb_component(), false);
     to_keep[database::water_id] = true;
     to_keep[database::electron_id] = true;
 
     // aqueous species
     // ---------------
     for (auto& it: aqueous_species) {
         database::GenericAqueousSpecies aq =  db_search.find_aqueous_species(it.first);
         if (aq.type == database::AqueousSpeciesClass::Component)
         {
             to_keep[aq.id] = true;
         }
         else if (aq.type == database::AqueousSpeciesClass::Aqueous)
         {
             for (auto idc: raw_data->range_aqueous_component())
             {
                 if (raw_data->nu_aqueous(aq.id, idc) != 0)
                 {
                     to_keep[idc] = true;
                 }
             }
         }
         else if (aq.type == database::AqueousSpeciesClass::Compound)
         {
             for (auto idc: raw_data->range_aqueous_component())
             {
                 if (raw_data->nu_compound(aq.id, idc) != 0)
                 {
                     to_keep[idc] = true;
                 }
             }
         }
         else
         {
             throw std::invalid_argument("'"+it.first+"' is not a valid aqueous species.");
         }
     }
     // solid phases
     // ------------
     for (auto& it: solid_phase) {
         database::GenericSolidPhase sp =  db_search.find_solid_phase(it.first);
         if (sp.type == database::SolidPhaseClass::EquilibriumMineral)
         {
             for (auto idc: raw_data->range_aqueous_component())
             {
                 if (raw_data->nu_mineral(sp.id, idc) != 0)
                 {
                     to_keep[idc] = true;
                 }
             }
         }
         else if (sp.type == database::SolidPhaseClass::MineralKinetics)
         {
             for (auto idc: raw_data->range_aqueous_component())
             {
                 if (raw_data->nu_mineral_kinetic(sp.id, idc) != 0)
                 {
                     to_keep[idc] = true;
                 }
             }
         }
         else
         {
             throw std::invalid_argument("'"+it.first+"' is not a valid solid phase.");
         }
     }
     // extra components
     // -----------------
     for (auto& comp_label: keep_components) {
         const index_t idc = raw_data->get_id_component(comp_label);
+        specmicp_assert(idc != no_species);
         to_keep[idc] = true;
     }
 
     // get ids of components to keep
     std::vector<index_t> id_comp_to_keep;
     id_comp_to_keep.reserve(raw_data->nb_component());
 
     for (auto id: raw_data->range_component()) {
         if (to_keep[id]) {
             id_comp_to_keep.push_back(id);
         }
     }
     return id_comp_to_keep;
 
 }
 
 scalar_t ReactantBox::ReactantBoxImpl::factor_total_concentration(
         const units::UnitsSet& units_set
         ) const
 {
     scalar_t mol_fac    = units::scaling_factor(units_set.quantity);
     scalar_t length_fac = units::scaling_factor(units_set.length);
 
     return std::pow(length_fac, 3)/mol_fac;
 }
 
 
-void ReactantBox::set_charge_keeper(const std::string& charge_keeper)
+void ReactantBox::set_charge_keeper(std::string charge_keeper)
 {
+    auto id = m_data->get_id_component(charge_keeper);
+    if (id == no_species) // check that the component exist
+    {
+        // if not check if it is an element and find the corresponding component
+        id = m_data->get_id_component_from_element(charge_keeper);
+        if (id == no_species)
+        {
+            throw std::invalid_argument(
+                        "ReactantBox::set_charge_keeper : "
+                        "the argument (" + charge_keeper
+                        + ") is not a component nor an element");
+        }
+        charge_keeper = m_data->get_label_component(id);
+    }
+
     m_impl->charge_keeper = charge_keeper;
     add_component(charge_keeper); // so it is always in the system
 }
 void ReactantBox::set_inert_volume_fraction(scalar_t inert_volume_fraction)
 {
     m_impl->inert_vol_frac = inert_volume_fraction;
 }
 void ReactantBox::set_saturated_system()
 {
     m_impl->water_eq_type = WaterEquationType::SaturatedSystem;
 }
 void ReactantBox::disable_conservation_water()
 {
     m_impl->water_eq_type = WaterEquationType::NoEquation;
 }
 
 void ReactantBox::set_water_partial_pressure_model(
         water_partial_pressure_f h2o_pressure_model
         )
 {
     m_impl->h2o_partial_pressure_model = h2o_pressure_model;
 }
 void ReactantBox::add_fixed_fugacity_gas(
         std::string gas,
         std::string component,
         scalar_t fugacity
         )
 {
+    specmicp_assert(fugacity > 0);
     m_impl->fix_fugacity.emplace_back(gas, component, fugacity);
     add_component(component);
 }
 void ReactantBox::add_fixed_activity_component(
         std::string component,
         scalar_t activity
         )
 {
+    specmicp_assert(activity > 0);
     m_impl->fix_activity.emplace_back(component, activity);
     add_component(component);
 }
 void ReactantBox::add_fixed_molality_component(
         std::string component,
         scalar_t molality
         )
 {
+    specmicp_assert(molality > 0);
     m_impl->fix_molality.emplace_back(component, molality);
     add_component(component);
 }
 
 
 AdimensionalSystemConstraints ReactantBox::get_constraints(bool modify_db) const
 {
     database::DataContainer* raw_data = get_database().get();
 
     AdimensionalSystemConstraints constraints;
     constraints.total_concentrations = get_total_concentration(modify_db);
 
-    if (not m_impl->charge_keeper.empty())
+    if (m_impl->charge_keeper != "")
     {
         constraints.set_charge_keeper(
                     raw_data->get_id_component(m_impl->charge_keeper)
                     );
     }
     constraints.set_inert_volume_fraction(
                 m_impl->inert_vol_frac
                 );
     constraints.water_equation = m_impl->water_eq_type;
     if (m_impl->h2o_partial_pressure_model != nullptr)
     {
         constraints.set_water_partial_pressure_model(m_impl->h2o_partial_pressure_model);
     }
 
     for (auto& fug: m_impl->fix_fugacity)
     {
         constraints.add_fixed_fugacity_gas(
                     raw_data->get_id_gas(fug.gas),
                     raw_data->get_id_component(fug.component),
                     std::log10(fug.fugacity)
                     );
     }
     for (auto& act: m_impl->fix_activity)
     {
         constraints.add_fixed_activity_component(
                     raw_data->get_id_component(act.first),
                     std::log10(act.second)
                     );
     }
     for (auto& mol: m_impl->fix_molality)
     {
         constraints.add_fixed_molality_component(
                     raw_data->get_id_component(mol.first),
                     std::log10(mol.second)
                     );
     }
 
     return constraints;
 }
 
 
 } // end namespace specmicp
diff --git a/src/specmicp/problem_solver/reactant_box.hpp b/src/specmicp/problem_solver/reactant_box.hpp
index 853b333..2c19005 100644
--- a/src/specmicp/problem_solver/reactant_box.hpp
+++ b/src/specmicp/problem_solver/reactant_box.hpp
@@ -1,167 +1,167 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 #ifndef SPECMICP_SPECMICP_PBSOLVER_REACTANTBOX_HPP
 #define SPECMICP_SPECMICP_PBSOLVER_REACTANTBOX_HPP
 
 //! \file problem_solver/reactant_box.hpp
 //! \brief Handle initial composition, with units
 
 #include "specmicp_common/types.hpp"
 #include "specmicp_common/pimpl_ptr.hpp"
 
 #include "specmicp_common/physics/units.hpp"
 #include "specmicp_database/database_holder.hpp"
 
 namespace specmicp {
 
 struct AdimensionalSystemConstraints;
 
 using water_partial_pressure_f = std::function<scalar_t (scalar_t)>;
 
 //! \brief A set of reactants
 //!
 //! Take care of the units
 class SPECMICP_DLL_PUBLIC ReactantBox :
     public units::UnitBaseClass,
     public database::DatabaseHolder
 
 {
 public:
     ReactantBox(
         std::shared_ptr<database::DataContainer> raw_data,
         const units::UnitsSet& units_set
     );
     ~ReactantBox();
     ReactantBox(ReactantBox&& other);
 
     ReactantBox& operator=(ReactantBox&& other);
 
     // Solution
     // =========
     //! \brief Set the solution amount
     void set_solution(scalar_t value, std::string str_unit);
     //! \brief Set the solution amount
     void set_solution(
             scalar_t                   value,
             const units::AmountUnit&   unit
             );
 
     // Aqueous
     // =======
     //! \brief Add an aqueous species
     void add_aqueous_species(
             std::string name,
             scalar_t    value,
             std::string str_unit
             );
     //! \brief Add an aqueous species
     void add_aqueous_species(
             std::string                name,
             scalar_t                   value,
             const units::AmountUnit&   unit
             );
 
     // Solid
     // =====
     //! \brief Add a solid phase
     void add_solid_phase(
             std::string name,
             scalar_t    value,
             std::string str_unit
             );
     //! \brief Add a solid phase
     void add_solid_phase(
             std::string                name,
             scalar_t                   value,
             const units::AmountUnit&   unit
             );
 
     // Total concentration
     // ===================
     //! \brief add a component without corresponding species
     //!
     //! It is useful if the database must be modified
     void add_component(const std::string& name);
     //! \brief Return the total concentration vector
     Vector get_total_concentration(bool modify_db=false)  const;
 
     // Constraints
     // ===========
 
     //! \brief Set the charge keeper
-    void set_charge_keeper(const std::string& charge_keeper);
+    void set_charge_keeper(std::string charge_keeper);
     //! \brief Set the inert volume fraction (0 by default)
     void set_inert_volume_fraction(scalar_t inert_volume_fraction);
     //! \brief The system is saturated
     void set_saturated_system();
     //! \brief Disable the water conservation equation
     //!
     //! \warning probably not a good idead
     void disable_conservation_water();
     //! \brief Set the water partial pressure model
     void set_water_partial_pressure_model(water_partial_pressure_f h2o_pressure_model);
 
     //! \brief Add a fixed fugacity gas
     void add_fixed_fugacity_gas(
             std::string gas,
             std::string component,
             scalar_t    fugacity
             );
     //! \brief Add a fixed activity component
     void add_fixed_activity_component(
             std::string component,
             scalar_t    activity
             );
     //! \brief Add a fixed molality component
     //!
     //! Molality is in mol/kg !
     void add_fixed_molality_component(
             std::string component,
             scalar_t    molality
             );
 
 
     //! \brief Return a set of constraints to be passed to the specmicp solver
     AdimensionalSystemConstraints get_constraints(bool modify_db=true) const;
 
 private:
     struct SPECMICP_DLL_LOCAL ReactantBoxImpl;
     utils::pimpl_ptr<ReactantBoxImpl> m_impl;
 
 };
 
 } // end namespace specmicp
 
 #endif // SPECMICP_SPECMICP_PBSOLVER_REACTANTBOX_HPP
diff --git a/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp b/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp
index 464e992..7117fcb 100644
--- a/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp
+++ b/tests/reactmicp/systems/unsaturated/hdf5_unsaturated.cpp
@@ -1,155 +1,155 @@
 #include "catch.hpp"
 
 #include "reactmicp/io/hdf5_unsaturated.hpp"
 
 #include "specmicp/adimensional/adimensional_system_solution.hpp"
 #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
 
 #include "specmicp_database/database.hpp"
 
 #include "specmicp/problem_solver/reactant_box.hpp"
 #include "specmicp/adimensional/adimensional_system_solver.hpp"
 
 
 #include "reactmicp/systems/unsaturated/variables_interface.hpp"
 
 #include "dfpm/mesh.hpp"
 
 #include <string>
 
 #define PATH_TEST_FILES "test_hdf5_unsaturated"
 
 using namespace specmicp;
 using namespace specmicp::reactmicp::systems::unsaturated;
 
 static specmicp::database::RawDatabasePtr get_database()
 {
     static database::RawDatabasePtr raw_data {nullptr};
     if (raw_data == nullptr)
     {
         specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
         thedatabase.keep_only_components(
             {"H2O", "H[+]", "Ca[2+]", "Si(OH)4", "HCO3[-]"}
             );
         std::map<std::string, std::string> swap = {{"HCO3[-]", "CO2"},};
         thedatabase.swap_components(swap);
         raw_data = thedatabase.get_database();
         raw_data->freeze_db();
     }
     return raw_data;
 }
 
 static AdimensionalSystemSolution
 get_init_solution(
         database::RawDatabasePtr the_database,
         units::UnitsSet the_units
         )
 {
 
     specmicp::ReactantBox reactant_box(the_database, the_units);
     reactant_box.set_solution(0.5, "L/L");
     reactant_box.add_solid_phase("C3S", 0.350, "L/L");
     reactant_box.add_solid_phase("C2S",  0.50, "L/L");
     reactant_box.add_aqueous_species("CO2", 32.72, "mol/kg");
-    reactant_box.set_charge_keeper("HO[-]");
+    reactant_box.set_charge_keeper("H");
 
 
     auto constraints = reactant_box.get_constraints(true);
 
     AdimensionalSystemSolver adim_solver(the_database, constraints);
     adim_solver.get_options().solver_options.set_tolerance(1e-14);
     adim_solver.get_options().units_set = the_units;
 
     Vector x;
     adim_solver.initialize_variables(x, 0.8, -3);
     micpsolver::MiCPPerformance perf = adim_solver.solve(x);
     if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
     {
         throw std::runtime_error("Error solving the initial problem");
     }
 
     return adim_solver.get_raw_solution(x);
 }
 
 
 
 TEST_CASE("hdf5 output", "[hdf5]") {
 
     std::string filename = "test_hdf5_unsaturated.hdf5";
 
     database::RawDatabasePtr raw_data = get_database();
 
     mesh::Uniform1DMeshGeometry mesh_geom;
     mesh_geom.nb_nodes = 10;
     mesh_geom.dx = 1;
     mesh_geom.section = 1;
     mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d(mesh_geom);
 
 
     units::UnitsSet the_units;
     the_units.length = units::LengthUnit::centimeter;
 
     specmicp::AdimensionalSystemSolution init_sol
             = get_init_solution(raw_data, the_units);
 
     specmicp::AdimensionalSystemSolutionExtractor extr(init_sol, raw_data, the_units);
 
     SECTION("Writer") {
 
 
         VariablesInterface vars_init(the_mesh, raw_data,
             {raw_data->get_id_component_from_element("C")}, the_units);
 
         for (auto node: the_mesh->range_nodes())
         {
             vars_init.initialize_variables(node, extr);
         }
 
         auto vars = vars_init.get_variables();
         io::UnsaturatedHDF5Saver saver(filename, vars, the_units);
         saver.save_timestep(10.0);
     }
 
 
 
     SECTION("Reader") {
 
         specmicp::io::UnsaturatedHDF5Reader reader(filename);
 
         // units
         units::UnitsSet read_unit = reader.get_units();
         CHECK(read_unit.length   == the_units.length);
         CHECK(read_unit.mass     == the_units.mass);
         CHECK(read_unit.quantity == the_units.quantity);
 
         // adimensional solution
         specmicp::AdimensionalSystemSolution sol = reader.get_adim_solution(10.0, 1);
 
         REQUIRE(sol.is_valid);
         for (auto r: range(sol.main_variables.rows()))
         {
             if  (std::isfinite(sol.main_variables(r)))
             {
                 CHECK(sol.main_variables(r)       == Approx(init_sol.main_variables(r)).epsilon(1e-10));
             }
         }
         CHECK(sol.log_gamma.norm()            == Approx(init_sol.log_gamma.norm()).epsilon(1e-10));
         CHECK(sol.gas_fugacities.norm()       == Approx(init_sol.gas_fugacities.norm()).epsilon(1e-10));
         CHECK(sol.secondary_molalities.norm() == Approx(init_sol.secondary_molalities.norm()).epsilon(1e-10));
 
 
         // main variables
         Vector liq_sat = reader.main_variable_vs_x(10.0, "H2O", "liquid_saturation");
         CHECK(liq_sat.rows() == the_mesh->nb_nodes());
         CHECK(liq_sat(1) == Approx(extr.saturation_water()));
 
         Vector ca_s = reader.main_variable_vs_x(10.0, "Ca[2+]", "solid_concentration");
         CHECK(ca_s.rows() == the_mesh->nb_nodes());
         CHECK(ca_s(1) == Approx(extr.total_solid_concentration(raw_data->get_id_component("Ca[2+]"))));
 
 
 
 
 
     }
 }