diff --git a/src/specmicp/problem_solver/reactant_box.cpp b/src/specmicp/problem_solver/reactant_box.cpp
index 77037e8..54ce33f 100644
--- a/src/specmicp/problem_solver/reactant_box.cpp
+++ b/src/specmicp/problem_solver/reactant_box.cpp
@@ -1,754 +1,771 @@
 /* =============================================================================
 
  Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
  Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
  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 SPECMICP_DLL_LOCAL si_container {
     std::string mineral;
     std::string component;
     scalar_t saturation_index;
 
     si_container(
             const std::string& amineral,
             const std::string& acomponent,
             scalar_t t_si
             ):
         mineral(amineral), component(acomponent), saturation_index(t_si)
     {}
 };
 
 struct ReactantBox::ReactantBoxImpl {
     bool immobile_disabled {false};
     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};
     scalar_t water_param {-1.0};
     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;
     std::vector<si_container> fix_saturation_index;
     std::vector<std::pair<std::string, scalar_t>> mineral_upper_bound;
+    SurfaceConstraint surface_model {};
 
     // 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) const {
         throw std::invalid_argument(error_msg);
     }
     void assert_mass_solution() const {
         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;
 
     scalar_t parse_unit_aqueous_species(
             std::string                name,
             const units::AmountUnit&   unit
             ) const;
     scalar_t parse_unit_solid_phase(
             std::string                name,
             const units::AmountUnit&   unit
             ) 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
         )
 {
 
     const scalar_t factor = m_impl->parse_unit_aqueous_species(name, unit);
     const scalar_t molality = factor*value;
     m_impl->aqueous_species[name] += molality;
 }
 
 void ReactantBox::set_aqueous_species(
         std::string name,
         scalar_t    value,
         std::string str_unit
         )
 {
     set_aqueous_species(name, value, units::parse_amount_unit(str_unit));
 }
 
 void ReactantBox::set_aqueous_species(
         std::string                name,
         scalar_t                   value,
         const units::AmountUnit&   unit
         )
 {
 
     const scalar_t factor = m_impl->parse_unit_aqueous_species(name, unit);
     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
         )
 {
 
     const scalar_t factor = m_impl->parse_unit_solid_phase(name, unit);
     m_impl->solid_phase[name] += factor*value;
 }
 
 void ReactantBox::set_solid_phase(
         std::string name,
         scalar_t    value,
         std::string str_unit
         )
 {
     return set_solid_phase(name, value, units::parse_amount_unit(str_unit));
 }
 
 void ReactantBox::set_solid_phase(
         std::string                name,
         scalar_t                   value,
         const units::AmountUnit&   unit
         )
 {
 
     const scalar_t factor = m_impl->parse_unit_solid_phase(name, unit);
     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(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::set_fixed_saturation(scalar_t saturation)
 {
     if (saturation <= 0.0 or saturation >= 1.0) {
         throw std::invalid_argument("The saturation must be between 0 and 1"
                                 "(Value : " + std::to_string(saturation) + ")");
     }
     m_impl->water_eq_type = WaterEquationType::FixedSaturation;
     m_impl->water_param = saturation;
 }
 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);
 }
 void ReactantBox::add_fixed_saturation_index(
         std::string mineral,
         std::string component,
         scalar_t saturation_index
         )
 {
     m_impl->fix_saturation_index.emplace_back(mineral, component, saturation_index);
     add_component(component);
 }
 
 void ReactantBox::set_mineral_upper_bound(
         std::string mineral,
         scalar_t max_volume_fraction
         )
 {
     if (max_volume_fraction < 0 or (not std::isfinite(max_volume_fraction))) {
         throw std::invalid_argument("The saturation must a finite positive number"
                        "(Value : " + std::to_string(max_volume_fraction) + ")");
     }
     m_impl->mineral_upper_bound.emplace_back(mineral, max_volume_fraction);
 }
 
 void ReactantBox::disable_immobile_species()
 {
     m_impl->immobile_disabled = true;
 }
 
+
+void ReactantBox::disable_surface_model()
+{
+    m_impl->surface_model = SurfaceConstraint();
+}
+
+void ReactantBox::set_equilibrium_surface_model(Vector ssites_concentrations)
+{
+    m_impl->surface_model = make_equilibrium_surface_model(ssites_concentrations);
+}
+
+void ReactantBox::set_EDL_surface_model(Vector ssites_concentrations, scalar_t sorption_surface)
+{
+    m_impl->surface_model = make_EDL_surface_model(ssites_concentrations, sorption_surface);
+}
+
 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 ( m_impl->immobile_disabled ) {
         constraints.disable_immobile_species();
     }
     // Water
     switch (m_impl->water_eq_type) {
     case WaterEquationType::MassConservation:
             constraints.enable_conservation_water();
         break;
     case WaterEquationType::FixedSaturation:
         specmicp_assert(m_impl->water_param > 0.0 and m_impl->water_param < 1.0);
         constraints.set_fixed_saturation(m_impl->water_param);
         break;
     case WaterEquationType::SaturatedSystem:
         constraints.set_saturated_system();
         break;
     case WaterEquationType::NoEquation:
         constraints.disable_conservation_water();
     }
     if (m_impl->h2o_partial_pressure_model != nullptr)
     {
         constraints.set_water_partial_pressure_model(m_impl->h2o_partial_pressure_model);
     }
 
     // volume fraction
     constraints.set_inert_volume_fraction(
                 m_impl->inert_vol_frac
                 );
 
 
     // Components
 
     if (m_impl->charge_keeper != "")
     {
         constraints.set_charge_keeper(
                     raw_data->get_id_component(m_impl->charge_keeper)
                     );
     }
 
     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)
                     );
     }
     for (auto& si: m_impl->fix_saturation_index)
     {
         constraints.add_fixed_saturation_index(
                     raw_data->get_id_mineral(si.mineral),
                     raw_data->get_id_component(si.component),
                     si.saturation_index
                     );
     }
 
     for (auto& mup: m_impl->mineral_upper_bound)
     {
         constraints.set_mineral_upper_bound(
                     raw_data->get_id_mineral(mup.first),
                     mup.second
                     );
     }
-
+    constraints.surface_model = m_impl->surface_model;
     return constraints;
 }
 
 
 
 scalar_t ReactantBox::ReactantBoxImpl::parse_unit_aqueous_species(
         std::string                name,
         const units::AmountUnit&   unit
         ) const
 {
     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 = 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)
     {
         assert_mass_solution();
         const scalar_t molar_mass =  db_search.molar_mass_aqueous(name);
         factor = unit.factor_si / molar_mass / mass_solution;
     }
     else if (unit.type == units::AmountUnitType::NumberOfMoles)
     {
         assert_mass_solution();
         factor = unit.factor_si / mass_solution;
     }
     else
     {
         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.");
     }
     return factor;
 }
 
 scalar_t ReactantBox::ReactantBoxImpl::parse_unit_solid_phase(
         std::string                name,
         const units::AmountUnit&   unit
         ) const
 {
     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 = 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 = 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 = db_search.molar_volume_solid_phase(name);
         factor = unit.factor_si / molar_volume;
 
     }
     else
     {
         raise_unknown_unit("Unit for solid phase " + name + " is invalid.");
     }
     return factor;
 }
 
 } // end namespace specmicp
diff --git a/src/specmicp/problem_solver/reactant_box.hpp b/src/specmicp/problem_solver/reactant_box.hpp
index c4e448e..f9061c7 100644
--- a/src/specmicp/problem_solver/reactant_box.hpp
+++ b/src/specmicp/problem_solver/reactant_box.hpp
@@ -1,217 +1,238 @@
 /* =============================================================================
 
  Copyright (c) 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
  Copyright (c) 2017-2018 F. Georget <fabien.georget@epfl.ch> EPFL
  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:
     //! \brief Constructor
     ReactantBox(
         std::shared_ptr<database::DataContainer> raw_data,
         const units::UnitsSet& units_set
     );
     //! \brief Destructor
     ~ReactantBox();
     //! \brief Move constructor
     ReactantBox(ReactantBox&& other);
     //! \brief Move constructor
     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
             );
     //! \brief Set the amount of an aqueous species
     void set_aqueous_species(
             std::string name,
             scalar_t    value,
             std::string str_unit
             );
     //! \brief Set the amount an aqueous species
     void set_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
             );
     //! \brief Set the amount of a solid phase
     void set_solid_phase(
             std::string name,
             scalar_t    value,
             std::string str_unit
             );
     //! \brief Set the amount of a solid phase
     void set_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(std::string charge_keeper);
     //! \brief Set the inert volume fraction (0 by default)
+
+    // Water
+    // -----
+
     void set_inert_volume_fraction(scalar_t inert_volume_fraction);
     //! \brief The system is saturated
     void set_saturated_system();
     //! \brief The saturation of the system is fixed
     void set_fixed_saturation(scalar_t saturation);
     //! \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);
 
+    // Fixed activity/fugacity/molality
+    // --------------------------------
+
     //! \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 Add a fixed saturation index (SI) solid phase
     //!
     //! \param mineral the solid phase with fixed saturation index
     //! \param component the adjustable component to obtain the given SI
     //! \param saturation_index the fixed SI
     void add_fixed_saturation_index(
             std::string mineral,
             std::string component,
             scalar_t saturation_index
             );
 
+    // Surface
+    // -------
+
+    void disable_surface_model();
+    //! \brief Create a simple equilibrium surface model
+    void set_equilibrium_surface_model(Vector ssites_concentrations);
+    //! \brief Create an EDL surface model
+    void set_EDL_surface_model(Vector ssites_concentrations, scalar_t sorption_surface);
+
+    // Mineral
+    // -------
 
     //! \brief Set an upper bound for a mineral
     //!
     //! It's an upper bound on the volume fraction
     void set_mineral_upper_bound(
             std::string mineral,
             scalar_t max_volume_fraction
             );
 
     //! \brief Disable the immobile species while solving the system
     void disable_immobile_species();
 
+    // Access to low level API
+    // =======================
+
     //! \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