diff --git a/src/specmicp/adimensional/adimensional_system.hpp b/src/specmicp/adimensional/adimensional_system.hpp
index 62f7351..d638412 100644
--- a/src/specmicp/adimensional/adimensional_system.hpp
+++ b/src/specmicp/adimensional/adimensional_system.hpp
@@ -1,373 +1,373 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : adimensional_system.hpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien 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:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * 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.
     * Neither the name of the Princeton University 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 OWNER 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_ADIMENSIONALSYSTEM_HPP
 #define SPECMICP_ADIMENSIONALSYSTEM_HPP
 
-//! \file reduced_system.hpp The MiCP program to solve speciation
+//! \file adimensional_system.hpp The MiCP program to solve speciation
 
 #include "common.hpp"
 #include "database.hpp"
 #include "micpsolver/micpprog.hpp"
 #include "../simulation_box.hpp"
 #include "physics/units.hpp"
 #include "utils/options_handler.hpp"
 
 #include "adimensional_system_structs.hpp"
 
 //! \namespace specmicp main namespace for the SpecMiCP solver
 namespace specmicp {
 
 class AdimensionalSystemSolution;
 
 //! \brief The equilibrium speciation problem
 //!
 //! Represent the equilibrium speciation problem as a Mixed Complementarity program
 //! It is called reduced becaused it does not solve explicitely for the secondary species.
 //!
 class AdimensionalSystem :
         public micpsolver::MiCPProg<AdimensionalSystem>,
         public OptionsHandler<AdimensionalSystemOptions>,
         public units::UnitBaseClass
 {
 public:
     //! \brief Create a reduced system
     //!
     //! \param ptrthermo Shared pointer to the thermodynamic data
     //! \param totconc Vector containing the total concentrations (in mol/volume) for each components
     AdimensionalSystem(
             RawDatabasePtr ptrdata,
             const AdimensionalSystemBC& conditions,
             const AdimensionalSystemOptions& options=AdimensionalSystemOptions()
             );
 
     AdimensionalSystem(
             RawDatabasePtr ptrdata,
             const AdimensionalSystemBC& conditions,
             const AdimensionalSystemSolution& previous_solution,
             const AdimensionalSystemOptions& options
             );
 
     // Variables
     // ==========
 
     //! \brief Return the total number of variables
     index_t total_variables() const {return m_nb_tot_vars;}
     //! \brief Return the number of 'free' variables (not subject to complementarity conditions)
     index_t nb_free_variables() const {return m_nb_free_vars;}
     //! \brief Return the number of variables subjected to the complementarity conditions
     index_t nb_complementarity_variables() const {return m_nb_compl_vars;}
 
     // The linear system
     // ==================
 
     //! \brief Return the residual for the water conservation equation
     scalar_t residual_water(const Vector& x) const;
     //! \brief Return the residual for the conservation equation of component (i)
     //!
     //! For (i==0), water use the method : 'residual_water'
     scalar_t residual_component(const Vector& x, index_t i) const;
     //! \brief Equilibrium condition for the minerals \f$ 1 - \mathrm{SI}_l \f$
     scalar_t residual_mineral(const Vector& x, index_t m) const;
     //! \brief Charge conservation
     scalar_t residual_charge_conservation(const Vector& x) const;
     scalar_t residual_fixed_fugacity(const Vector& x, index_t component) const;
     scalar_t residual_fixed_activity(const Vector& x, index_t component) const;
     //! \brief Return the residuals
     void get_residuals(const Vector& x, Vector& residual);
     //! \brief Return the jacobian
     void get_jacobian(Vector& x,
                       Matrix& jacobian);
 
 
     // Getters
     // #######
 
     // Equation id access
     // ==================
 
     //! \brief Return the equation id
     index_t ideq(index_t i) const {return m_ideq[i];}
     //! \brief Return the equation number for conservation of water
     index_t ideq_w() const {return m_ideq[0];}
     //! \brief Return the equation number of component 'i'
     index_t ideq_paq(index_t i) const {
         specmicp_assert(i < m_data->nb_component);
         return m_ideq[i];
     }
     //! \brief Return the equation number of mineral 'm'
     index_t ideq_min(index_t m) const {
         specmicp_assert(m < m_data->nb_mineral);
         return m_ideq[m+m_data->nb_component];
     }
 
     // Equation type
     // -------------
 
     WaterEquationType water_equation_type() const {
         return static_cast<WaterEquationType>(m_component_equation_type[0]);
     }
     //! \brief Return the type of equation for aqueous species 'component'
     AqueousComponentEquationType aqueous_component_equation_type(index_t component) const {
         specmicp_assert(component >= 0 and component < m_data->nb_component);
         return static_cast<AqueousComponentEquationType>(m_component_equation_type[component]);
     }
 
     //! Return true if an equation exist for this component
     bool is_active_component(index_t component) const {
         specmicp_assert(component < m_data->nb_component);
         return m_ideq[component] != no_equation;
     }
 
     // Variables
     // =========
 
     // Primary
     // -------
 
     //! \brief Return the saturation of water
     scalar_t saturation_water(const Vector& x) const;
     //! \brief Return the density of water
     scalar_t density_water() const;
     //! \brief Return the saturation (mol/volume) of a mineral
     scalar_t saturation_mineral(const Vector& x, index_t mineral) const;
     //! \brief Return the volume of minerals
     scalar_t molar_volume_mineral(index_t mineral) const;
     //! \brief Return the sum of saturation of the minerals
     scalar_t sum_saturation_minerals(const Vector& x) const;
     //! \brief Return the log_10 of component concentration
     scalar_t log_component_molality(const Vector& x, index_t component) const
     {
         specmicp_assert(ideq_paq(component) != no_equation and component < m_data->nb_component);
         return x(ideq_paq(component));
     }
     //! \brief Return  the concentration of 'component'
     scalar_t component_molality(const Vector& x, index_t component) const
     {
         return pow10(log_component_molality(x, component));
     }
 
     // Secondary
     // ---------
 
     //! \brief Return the concentration of secondary specis 'aqueous'
     scalar_t secondary_molality(index_t aqueous) const
     {
         if (not m_active_aqueous[aqueous]) return 0.0;
         return m_secondary_conc(aqueous);
     }
     //! \brief Return true if 'aqueous' is active
     //!
     //! i.e. Return true if all its component are present in the system
     bool is_aqueous_active(index_t aqueous) const
     {
         return m_active_aqueous[aqueous];
     }
     //! \brief Return log_10(γ_i) where i is a component
     scalar_t log_gamma_component(index_t component) const {
         return m_loggamma(component);
     }
     //! \brief Return log_10(γ_j) where j is a secondary aqueous species
     scalar_t log_gamma_secondary(index_t aqueous) const {
         return m_loggamma(m_data->nb_component + aqueous);
     }
     //! \brief Return the ionic strength
     scalar_t ionic_strength() const {
         return m_ionic_strength;
     }
     //! \brief Return the total concentration of 'component'
     scalar_t total_concentration_bc(index_t component) const {
         specmicp_assert(m_aqueous_component_equation_type[component]
                         == AqueousComponentEquationType::MassConservation);
         return m_fixed_values(component);
     }
     //! \brief Return the total concentration of 'component'
     scalar_t fixed_activity_bc(index_t component) const {
         specmicp_assert(m_aqueous_component_equation_type[component]
                         == AqueousComponentEquationType::FixedActivity);
         return m_fixed_values(component);
     }
     //! \brief Return the total concentration of 'component'
     scalar_t fixed_fugacity_bc(index_t component) const {
         specmicp_assert(m_aqueous_component_equation_type[component]
                         == AqueousComponentEquationType::FixedFugacity);
         return m_fixed_values(component);
     }
 
     //! \brief Return the fugacity of 'gas'
     scalar_t gas_fugacity(index_t gas) const {
         return m_gas_fugacity(gas);
     }
     //! \brief Return true if gas is active
     bool is_active_gas(index_t gas) const {
         return m_active_gas[gas];
     }
 
     //! \brief Return the saturation of the gas phase
     scalar_t saturation_gas_phase() const {
         return m_saturation_gas;
     }
 
     //! \brief Return the gas pressure
     scalar_t gas_total_pressure() const {
         return 1.01325e5;
     }
     //! \brief Return the temperature
     scalar_t temperature() const {
         return 273.16+25;
     }
 
 
     // Equilibrium State
     // =================
 
     //! \brief Return the equilibrium state of the system, the 'Solution' of the speciation problem
     AdimensionalSystemSolution get_solution(const Vector& xtot, const Vector& x);
 
 
     // Algorithm
     // #########
 
     //! \brief Compute the ionic strength
     void set_ionic_strength(const Vector& x);
     //! \brief Compute the activity coefficients
     void compute_log_gamma(const Vector& x);
 
     // For the normal algorithm
     void set_secondary_concentration(const Vector& x);
     //! \brief Compute the secondary variables
     void set_secondary_variables(const Vector& x);
     //! \brief Set the gas phase saturation
     void set_saturation_gas_phase(const Vector& x);
     //! \brief Set the fugacity pressure
     void set_pressure_fugacity(const Vector& x);
 
     //! \brief Compute the activity model, called by the solver at the beginning of an iteration
     bool hook_start_iteration(const Vector &x, scalar_t norm_residual);
 
     //! \brief Return a scale factor to avoid negative mass during Newton's iteration
     double max_lambda(const Vector &x, const Vector &update);
     //! \brief Return the component that does not exist in the solution (inactive dof)
     const std::vector<index_t>& get_non_active_component() {return m_nonactive_component;}
 
     //! \brief A reasonable (well... maybe...) starting guess
     void reasonable_starting_guess(Vector& x);
     //! \brief A reasonable (maybe...) restarting guess
     void reasonable_restarting_guess(Vector& x);
 
     //! \brief Return true if 'component' is the charge keeper
     //!
     //! The charge keeper is the component added or removed to satisfy the charge balance
     bool is_charge_keeper(index_t component)
     {
         return (aqueous_component_equation_type(component) == AqueousComponentEquationType::ChargeBalance);
     }
 
 // Private Members and attributes
 // ##############################
 private:
 
     // Jacobian
     // ########
 
     void jacobian_water(Vector& x, Matrix& jacobian);
     void jacobian_aqueous_components(Vector& x, Matrix& jacobian);
     void jacobian_minerals(Vector& x, Matrix& jacobian);
 
     // Equation numbering
     // ##################
 
     void number_eq(const AdimensionalSystemBC& condition);
     void number_eq_aqueous_component(
             const AdimensionalSystemBC& conditions,
             index_t& neq
             );
 
     // Setter
     // ######
 
     // main variables
     // --------------
 
     // they are set by the solver
 
     // Secondary variables
     // -------------------
 
     //! \brief Return a reference to log_10(γ_i) where i is a component
     scalar_t& log_gamma_component(index_t component) {
         return m_loggamma(component);
     }
     //! \brief Return a reference to log_10(γ_j) where j is a secondary aqueous species
     scalar_t& log_gamma_secondary(index_t aqueous) {
         return m_loggamma(m_data->nb_component + aqueous);
     }
     //! \brief Return a reference to the ionic strength
     scalar_t& ionic_strength() {
         return m_ionic_strength;
     }
 
     // Attributes
     // ##########
 
     RawDatabasePtr m_data;
     Vector m_fixed_values;
 
     index_t m_nb_tot_vars;
     index_t m_nb_free_vars;
     index_t m_nb_compl_vars;
 
     std::vector<index_t> m_ideq;
     std::vector<int> m_component_equation_type;
     std::vector<index_t> m_fixed_activity_species;
 
     Vector m_secondary_conc;
 
     Vector m_loggamma;
     scalar_t m_ionic_strength;
 
     scalar_t m_saturation_gas;
     Vector m_gas_fugacity;
 
     bool not_in_linesearch;
 
     std::vector<index_t> m_nonactive_component;
     std::vector<bool> m_active_aqueous;
     std::vector<bool> m_active_gas;
 
 };
 
 } // end namespace specmicp
 
 
 #endif // SPECMICP_ADIMENSIONALSYSTEM_HPP
diff --git a/src/specmicp/adimensional/adimensional_system_solution.hpp b/src/specmicp/adimensional/adimensional_system_solution.hpp
index dc5f462..ee37a88 100644
--- a/src/specmicp/adimensional/adimensional_system_solution.hpp
+++ b/src/specmicp/adimensional/adimensional_system_solution.hpp
@@ -1,37 +1,39 @@
 #ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTION_HPP
 #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTION_HPP
 
 #include "common.hpp"
 
+//! \file adimensional_system_solution.hpp Structure to contain the solution returned by AdimensionalSystemSolution
+
 namespace specmicp {
 
 //! \brief Solution of an adimensional system
 struct AdimensionalSystemSolution
 {
     Vector main_variables;       //!< The main variables
     Vector secondary_molalities; //!< Molalities of secondary aqueous
     Vector log_gamma;            //!< Log_10 of activities coefficients
     scalar_t ionic_strength;     //!< Ionic strength
     Vector gas_fugacities;     //!< Gas fugacities
 
     AdimensionalSystemSolution() {}
 
     AdimensionalSystemSolution(
             const Vector& variables,
             const Vector& aqueous_molalities,
             const Vector& loggama,
             scalar_t the_ionic_strength,
             const Vector& fugacities
             ):
         main_variables(variables),
         secondary_molalities(aqueous_molalities),
         log_gamma(loggama),
         ionic_strength(the_ionic_strength),
         gas_fugacities(fugacities)
     {}
 
 };
 
 } // end namespace specmicp
 
 #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTION_HPP
diff --git a/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp b/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp
index d812da3..a9c72f5 100644
--- a/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp
+++ b/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp
@@ -1,192 +1,194 @@
 #ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP
 #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP
 
 #include "common.hpp"
 #include "database.hpp"
 #include "adimensional_system_solution.hpp"
 #include "physics/units.hpp"
 
+//! \file adimensional_system_solution_extractor.hpp Obtain data from AdimensionalSystemSolution
+
 namespace specmicp {
 
+//! \brief Obtain physical variables from AdimensionalSystemSolution
 class AdimensionalSystemSolutionExtractor: public units::UnitBaseClass
 {
 public:
     // public members
     // ###############
 
     AdimensionalSystemSolutionExtractor(
             const AdimensionalSystemSolution& solution,
             RawDatabasePtr thedatabase,
             units::UnitsSet units_set):
         units::UnitBaseClass(units_set),
         m_solution(solution),
         m_data(thedatabase)
     {}
 
-
     // Primary variables
     // =================
 
     // water
     // ------
 
     //! \brief Return the total saturation of water
     scalar_t total_saturation_water() {
         return m_solution.main_variables(id_water());}
 
     // component
     // ---------
 
     //! \brief Return the log_10 of the molality of 'component'
     scalar_t log_molality_component(index_t component) {
         return m_solution.main_variables(id_component(component));}
 
     //! \brief Return the molality of 'component'
     scalar_t molality_component(index_t component) {
         return pow10(log_molality_component(component));
     }
 
     // solid phases
     // ------------
 
     //! \brief Return the total saturation of 'mineral'
     scalar_t total_saturation_mineral(index_t mineral) {
         return m_solution.main_variables(id_mineral(mineral));
     }
 
     // Secondary variables
     // ===================
 
     //! \brief Return the ionic strength of the solution
     scalar_t ionic_strength() {return m_solution.ionic_strength;}
 
     // component
     // ---------
 
     //! \brief Return the log_10 of the activity coefficient of 'component'
     scalar_t log_activity_coefficient_component(index_t component) {
         return m_solution.log_gamma(id_component_gamma(component));
     }
     //! \brief Return the activity coefficient of 'component'
     scalar_t activity_coefficient_component(index_t component) {
         return pow10(log_activity_coefficient_component(component));
     }
     //! \brief Return the log_10 of the activity of 'component
     scalar_t log_activity_component(index_t component) {
         return log_molality_component(component) + log_activity_coefficient_component(component);
     }
     //! \brief Return the activity of 'component'
     scalar_t activity_component(index_t component) {
         return pow10(log_activity_component(component));
     }
 
     // aqueous species
     // ---------------
     //! \brief Return the molality of secondary specis 'aqueous'
     scalar_t molality_aqueous(index_t aqueous) {
         return m_solution.secondary_molalities(id_aqueous(aqueous));}
     //! \brief Return the log10 of the activity ocefficient of secondary species 'aqueous'
     scalar_t log_activity_coefficient_aqueous(index_t aqueous) {
         return m_solution.log_gamma(id_aqueous_gamma(aqueous));
     }
     //! \brief Return the activity coefficient of secondary species 'aqueous'
     scalar_t activity_coefficient_aqueous(index_t aqueous) {
         return pow10(log_activity_coefficient_aqueous(aqueous));
     }
     //! \brief Return the activity of secondary species 'aqueous'
     scalar_t activity_aqueous(index_t aqueous) {
         return activity_coefficient_aqueous(aqueous)*molality_aqueous(aqueous);
     }
 
     // gas fugacity
     // -------------
     //! \brief Return fugacity for 'gas'
     scalar_t gas_fugacity(index_t gas) {
         return m_solution.gas_fugacities(gas);
     }
 
     // Tertiary variables
     // ==================
 
     // Water
     // -----
 
     //! \brief Return the saturation of water
     scalar_t saturation_water() {return total_saturation_water()*porosity();}
 
     //! \brief Return the mass concentration of water
     scalar_t mass_concentration_water();
 
     //! \brief Return the pH of the solution
     scalar_t pH();
 
     // Component
     // ---------
 
     //! \brief Return the total aqueous concentration
     scalar_t total_aqueous_concentration(index_t component);
     //! \brief Return the total solid concentration
     scalar_t total_solid_concentration(index_t component);
 
     // Gas phase
     // ---------
 
     //! \brief Return the saturation of the gas phase
     scalar_t saturation_gas_phase() {return 1-saturation_water();}
 
     //! \brief Return the total saturation of the gas phase
     scalar_t total_saturation_gas_phase() {return porosity() - saturation_water();}
 
     // Component
     // ---------
 
     // Solid phases
     // ------------
 
     //! \brief Return the concentration of 'mineral'
     scalar_t mole_concentration_mineral(index_t mineral);
 
     //! \brief Return the mass concentration of 'mineral'
     scalar_t mass_concentration_mineral(index_t mineral);
 
     //! \brief Return the total saturation of all minerals
     scalar_t total_saturation_minerals();
 
     //! \brief Return the porosity
     scalar_t porosity() {return 1 - total_saturation_minerals();}
 
 private:
     // private members
     // ###############
 
     index_t id_water() {return 0;}
     index_t id_component(index_t component) {
         specmicp_assert(component < m_data->nb_component);
         return component;
     }
     index_t id_mineral(index_t mineral) {
         specmicp_assert(mineral < m_data->nb_mineral);
         return m_data->nb_component+mineral;}
     index_t id_component_gamma(index_t component) {
         specmicp_assert(component < m_data->nb_component);
         return component;
     }
     index_t id_aqueous(index_t aqueous) {
         specmicp_assert(aqueous < m_data->nb_aqueous);
         return aqueous;
     }
     index_t id_aqueous_gamma(index_t aqueous) {
         specmicp_assert(aqueous < m_data->nb_aqueous);
         return m_data->nb_component+aqueous;
     }
 
 private:
     // private attributes
     // ##################
     AdimensionalSystemSolution m_solution;
     RawDatabasePtr m_data;
 };
 
 } // end namespace specmicp
 
 #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP
diff --git a/src/specmicp/adimensional/adimensional_system_solver.hpp b/src/specmicp/adimensional/adimensional_system_solver.hpp
index 87fc255..7b15daa 100644
--- a/src/specmicp/adimensional/adimensional_system_solver.hpp
+++ b/src/specmicp/adimensional/adimensional_system_solver.hpp
@@ -1,111 +1,114 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : reduced_system_solver.hpp
  - Author : Fabien Georget
 
 Copyright (c) 2014, Fabien 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:
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * 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.
     * Neither the name of the Princeton University 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 OWNER 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_ADIMENSIONALSYSTEMSOLVER_HPP
 #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP
 
 #include "database.hpp"
 #include "adimensional_system.hpp"
 #include "adimensional_system_solver_structs.hpp"
 #include "utils/options_handler.hpp"
 
-//! \file reduced_system_solver.hpp Solve a reduced system
+//! \file adimensional_system_solver.hpp Solve a reduced system
 
 namespace specmicp {
 
 // forward declaration
 class AdimensionalSystemSolution;
 
-//! \brief Solver a reduced system
+//! \brief Solve an adimensional system
 //!
-//! Take care of possibly non-existing component in the system
+//! Take care of non-existing component in the system
+//! and restart the computation if necessary
 class AdimensionalSystemSolver: public OptionsHandler<AdimensionalSystemSolverOptions>
 {
 public:
     using SystemPtr = std::shared_ptr<AdimensionalSystem>;
 
-    //! Default constructor - you need to use another method to initialize it correctly
+    //! Default constructor
+    //!
+    //! Another constructor is necessary
     AdimensionalSystemSolver() {}
 
     AdimensionalSystemSolver(RawDatabasePtr data,
                         const AdimensionalSystemBC& conditons,
                         AdimensionalSystemSolverOptions options=AdimensionalSystemSolverOptions()
                         );
 
     AdimensionalSystemSolver(RawDatabasePtr data,
                         const AdimensionalSystemBC& conditions,
                         const AdimensionalSystemSolution& previous_solution,
                         AdimensionalSystemSolverOptions options=AdimensionalSystemSolverOptions()
                         );
 
     //! \brief solve the problem using initial guess x
     //!
     //! \param[in,out] x   in -> initial guess, out -> solution
     //! \param init if true, the algorithm guess a starting point
     micpsolver::MiCPPerformance solve(Vector& x, bool init=false);
 
     //! \brief Return the system used for the computation
     SystemPtr get_system() {return m_system;}
 
     //! \brief Return the solution in a manageable form
     AdimensionalSystemSolution get_raw_solution(const Eigen::VectorXd& x);
 
 private:
     //! \brief set up the true variable vector
     void set_true_variable_vector(const Vector& x);
     //! \brief set up the true solution vector
     //!
     //! add zero components
     void set_return_vector(Vector& x);
     //! \brief solve the problem
     micpsolver::MiCPPerformance solve();
 
     RawDatabasePtr m_data; //! The raw database
     SystemPtr m_system; //! The system to solve
     Vector m_var; //! Copy of the solution vector (necessary in case of failing)
 };
 
 //! \brief Solve a reduced system, function provided for convenience
 inline micpsolver::MiCPPerformance solve_reduced_system(
         std::shared_ptr<database::DataContainer> data,
         Vector& totconc,
         Vector& x
         )
 {
     return AdimensionalSystemSolver(data, totconc).solve(x);
 }
 
 
 } // end namespace specmicp
 
 #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP
diff --git a/src/specmicp/adimensional/adimensional_system_solver_structs.hpp b/src/specmicp/adimensional/adimensional_system_solver_structs.hpp
index d6e6b77..a43b193 100644
--- a/src/specmicp/adimensional/adimensional_system_solver_structs.hpp
+++ b/src/specmicp/adimensional/adimensional_system_solver_structs.hpp
@@ -1,32 +1,34 @@
 #ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVERSTRUCTS_HPP
 #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVERSTRUCTS_HPP
 
 #include "common.hpp"
 #include "micpsolver/micpsolver_structs.hpp"
 #include "adimensional_system_structs.hpp"
 #include "physics/units.hpp"
 
+//! \file adimensional_system_solver_structs.hpp Options and conditions for AdimensionalSystemSolver
+
 namespace specmicp {
 
 //! \brief Options for the Equilibrium solver
 //!
 //! Most of the options are contained in the MiCP solver options
 struct AdimensionalSystemSolverOptions
 {
     micpsolver::MiCPSolverOptions solver_options; //!< Options of the MiCP solver
     AdimensionalSystemOptions     system_options; //!< Options of the system
     bool allow_restart;      //!< Allow the restarting if the problem failed the first part
     units::UnitsSet units_set;         //!< Set of units
 
     AdimensionalSystemSolverOptions():
         solver_options(), // default options of the solver
         system_options(), // default options of the system
         allow_restart(true)
     {
     }
 };
 
 
 } // end namespace specmicp
 
 #endif //SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVERSTRUCTS_HPP
diff --git a/src/specmicp/adimensional/adimensional_system_structs.hpp b/src/specmicp/adimensional/adimensional_system_structs.hpp
index e009693..4e57b66 100644
--- a/src/specmicp/adimensional/adimensional_system_structs.hpp
+++ b/src/specmicp/adimensional/adimensional_system_structs.hpp
@@ -1,123 +1,125 @@
 #ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP
 #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP
 
 #include "common.hpp"
 
+//! \file adimensional_system_structs.hpp Options and BCs for the AdimensionalSystem
+
 namespace specmicp {
 
 //! \brief Options for the Adimensional Systems
 struct AdimensionalSystemOptions
 {
-    bool non_ideality;       //!< Solve for non ideality
-    scalar_t non_ideality_tolerance; //!< Tolerance for non ideality
-    index_t non_ideality_max_iter;  //!< Max iterations fornon ideality
+    bool non_ideality;                //!< Solve for non ideality
+    scalar_t non_ideality_tolerance;  //!< Tolerance for non ideality
+    index_t non_ideality_max_iter;    //!< Max iterations fornon ideality
     scalar_t under_relaxation_factor; //!< Under relaxation factor for the conservation of water
 
 
     AdimensionalSystemOptions():
         non_ideality(true),
         non_ideality_tolerance(1e-8),
         non_ideality_max_iter(10),
         under_relaxation_factor(0.9)
     {}
 };
 
 //! \brief Type of an aqueous component equation
 enum class AqueousComponentEquationType
 {
-    NoEquation = no_equation,         //!< Not an equation, component is not present in the system
+    NoEquation = no_equation, //!< Not an equation, component is not present in the system
     MassConservation,   //!< Mass balance
     ChargeBalance,      //!< M.B. replaced by charge balance
     FixedFugacity,      //!< M.B. replaced by a fixed fugacity equation
     FixedActivity       //!< M.B. replaced by a fixed activity equation
 };
 
 enum class WaterEquationType
 {
-    NoEquation = no_equation,   //!< Amount of water is not solved
+    NoEquation = no_equation,  //!< Amount of water is not solved
     MassConservation,   //!< Water is conserved
     SaturatedSystem     //!< System is saturated
 };
 
 //! \brief Struct to contain information needed to solve a fix fugacity problem
 struct FixedFugacityBC
 {
-    index_t id_gas;
-    index_t id_component;
-    scalar_t log_value;
+    index_t id_gas;       //!< Id of the fixed-fugacity gas
+    index_t id_component; //!< Id of the corresponding component
+    scalar_t log_value;   //!< Logarithm of the fugacity
 
     FixedFugacityBC(index_t gas, index_t component, scalar_t logvalue):
         id_gas(gas),
         id_component(component),
         log_value(logvalue)
     {}
 };
 
 //! \brief Struct to contain information needed to solve a fix activity problem.
 struct FixedActivityBC
 {
-    index_t id_component;
-    scalar_t log_value;
+    index_t id_component; //!< Id of the fixed-activity component
+    scalar_t log_value;   //!< Log of the activity
 
     FixedActivityBC(index_t component, scalar_t logvalue):
         id_component(component),
         log_value(logvalue)
     {}
 
 };
 
 //! \brief Struct to contains the "Boundary conditions" for
 struct AdimensionalSystemBC
 {
-    Vector total_concentrations; //!< Total concentrations
+    Vector total_concentrations;      //!< Total concentrations
     WaterEquationType water_equation; //!< Water equation
     index_t charge_keeper;       //!< The equation for this component is replace by the charge balance
     bool saturated_system;       //!> System is saturated - no gas phase
     std::vector<FixedFugacityBC> fixed_fugacity_bc; //!< Contains information about fixed fugacity gas
     std::vector<FixedActivityBC> fixed_activity_bc; //!< Contains information about fixed activity component
 
 
     AdimensionalSystemBC():
         water_equation(WaterEquationType::MassConservation),
         charge_keeper(no_species)
     {}
 
    AdimensionalSystemBC(const Vector& total_concs):
        total_concentrations(total_concs),
        water_equation(WaterEquationType::MassConservation),
        charge_keeper(no_species)
      {}
 
     //! \brief Enable the conservation of water
     void enable_conservation_water() {water_equation = WaterEquationType::MassConservation;}
     //! \brief Disable the conservation of water
     void disable_conservation_water() {water_equation = WaterEquationType::NoEquation;}
     //! \brief The system is saturated
     void set_saturated_system() {water_equation = WaterEquationType::SaturatedSystem;}
 
     //! \brief Set the charge keeper to 'component'
     void set_charge_keeper(index_t component) {
         charge_keeper = component;
     }
 
     //! \brief Add a fixed fugacity gas condition
     void add_fixed_fugacity_gas(const FixedFugacityBC& bc) {
         fixed_fugacity_bc.push_back(bc);
     }
     //! \brief Add a fixed fugacity gas condition
     void add_fixed_fugacity_gas(index_t gas, index_t component, scalar_t logvalue) {
         fixed_fugacity_bc.push_back(FixedFugacityBC(gas, component, logvalue));
     }
     //! \brief Add a fixed activity component condition
     void add_fixed_activity_component(const FixedActivityBC& bc) {
         fixed_activity_bc.push_back(bc);
     }
     //! \brief Add a fixed activity component condition
     void add_fixed_activity_component(index_t component, scalar_t log_value) {
         fixed_activity_bc.push_back(FixedActivityBC(component, log_value));
     }
 };
 
 } // end namespace specmicp
 
 #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP
diff --git a/src/specmicp/adimensional/equilibrium_curve.hpp b/src/specmicp/adimensional/equilibrium_curve.hpp
index 854f039..7167db6 100644
--- a/src/specmicp/adimensional/equilibrium_curve.hpp
+++ b/src/specmicp/adimensional/equilibrium_curve.hpp
@@ -1,102 +1,104 @@
 #ifndef SPECMICP_SPECMICP_EQUILIBRIUMCURVE_HPP
 #define SPECMICP_SPECMICP_EQUILIBRIUMCURVE_HPP
 
 #include "common.hpp"
 #include "database.hpp"
 #include "adimensional_system_solver_structs.hpp"
 #include "adimensional_system_solution.hpp"
 
+//! \file equilibrium_curve.hpp ABC to compute a reaction path
+
 namespace specmicp {
 
 //! \brief Abstract Base Class to compute an equilibrium curve
 //!
 //!
 class EquilibriumCurve
 {
 public:
     EquilibriumCurve() {}
 
     EquilibriumCurve(RawDatabasePtr database,
                      AdimensionalSystemBC conditions
                      ):
         m_data(database),
         m_conditions(conditions)
 
     {}
 
     //! \brief Run one step of the equilibrium curve
     void run_step();
     //! \brief Return the raw database
     RawDatabasePtr database() const {return m_data;}
     //! \brief Return the current solution
     const AdimensionalSystemSolution& current_solution() const {
         return m_current_solution;
     }
     //! \brief Return the current solution vector
     Vector& solution_vector() {
         return m_solution_vector;
     }
     //! \brief Return the current MiCPsolver performance
     const micpsolver::MiCPPerformance& solver_performance() const {
         return m_current_perf;
     }
     //! \brief Return a const reference to the solver options
     const AdimensionalSystemSolverOptions& solver_options() const {
         return m_solver_options;
     }
     //! \brief Return a const reference to the BC of AdimensionalSystem
     const AdimensionalSystemBC& conditions() const {
         return m_conditions;
     }
     //! \brief update the problem;
     virtual void update_problem() = 0;
 
     //! \brief Output - do nothing by default
     virtual void output() {}
 
     //! \brief How to handle an error
     virtual void error_handling(std::string msg) const;
 
     //! \brief Solve the first problem
     void solve_first_problem();
 
 protected:
     //! \brief Set the database
     void set_database(RawDatabasePtr the_database) {
         m_data =the_database;
     }
     //! \brief Read-write reference to the conditions
     AdimensionalSystemBC& conditions() {
         return m_conditions;
     }
     //! \brief Read-write reference to the options of AdimensionalSystemSolver
     AdimensionalSystemSolverOptions& solver_options() {
         return m_solver_options;
     }
 private:
     //! \brief Set the problem, according user input
     void set_problem();
     //! \brief Solve the problem
     //!
     //! May fail, call error_handling
     void solve_problem();
     //! \brief Post processing of the data
     void post_processing();
 
     RawDatabasePtr m_data;
 
     AdimensionalSystemBC m_conditions;
 
     AdimensionalSystemSolverOptions m_solver_options;
 
     AdimensionalSystemSolution m_current_solution;
 
     micpsolver::MiCPPerformance m_current_perf;
 
     Vector m_solution_vector;
 
 };
 
 } // end namespace specmicp
 
 #endif // SPECMICP_SPECMICP_EQUILIBRIUMCURVE_HPP