diff --git a/examples/specmicp/adimensional/thermocarbo.cpp b/examples/specmicp/adimensional/thermocarbo.cpp
index f4c161c..a7680ff 100644
--- a/examples/specmicp/adimensional/thermocarbo.cpp
+++ b/examples/specmicp/adimensional/thermocarbo.cpp
@@ -1,434 +1,434 @@
 
 #include <iostream>
 
 #include "utils/log.hpp"
 #include "specmicp/adimensional/adimensional_system_solver.hpp"
 #include "specmicp/adimensional/adimensional_system_solution.hpp"
 #include "specmicp/problem_solver/formulation.hpp"
 #include "specmicp/problem_solver/dissolver.hpp"
 #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
 
 #include "specmicp/adimensional/equilibrium_curve.hpp"
 
 void fixed_fugacity_thermocarbo()
 {
     specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
     std::map<std::string, std::string> swapping ({
                                               {"H[+]","HO[-]"},
                                               {"Si(OH)4", "SiO(OH)3[-]"},
                                               {"Al[3+]","Al(OH)4[-]"},
                                                 });
     thedatabase.swap_components(swapping);
     specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
 
     specmicp::Formulation formulation;
     specmicp::scalar_t mult = 5e3;
 
     specmicp::scalar_t m_c3s = mult*0.6;
     specmicp::scalar_t m_c2s = mult*0.2;
     specmicp::scalar_t m_c3a = mult*0.10;
     specmicp::scalar_t m_gypsum = mult*0.10;
     specmicp::scalar_t wc = 0.8;
     specmicp::scalar_t m_water = wc*1e-3*(
                   m_c3s*(3*56.08+60.08)
                 + m_c2s*(2*56.06+60.08)
                 + m_c3a*(3*56.08+101.96)
                 + m_gypsum*(56.08+80.06+2*18.02)
                 );
 
     formulation.mass_solution = m_water;
     formulation.amount_minerals = {
         {"C3S", m_c3s},
         {"C2S", m_c2s},
         {"C3A", m_c3a},
         {"Gypsum", m_gypsum}
     };
     formulation.extra_components_to_keep = {"HCO3[-]", };
     formulation.minerals_to_keep = {
         "Portlandite",
         "CSH,jennite",
         "CSH,tobermorite",
         "SiO2,am",
         "Calcite",
         "Al(OH)3,am",
         "Monosulfoaluminate",
         "Tricarboaluminate",
         "Monocarboaluminate",
         "Hemicarboaluminate",
         "Straetlingite",
         "Gypsum",
         "Ettringite",
         "Thaumasite"
     };
 
 
     specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
     specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
     specmicp::index_t id_ho  = thedatabase.component_label_to_id("HO[-]");
     specmicp::index_t id_hco3 = thedatabase.component_label_to_id("HCO3[-]");
     specmicp::index_t id_co2g = thedatabase.gas_label_to_id("CO2(g)");
 
     specmicp::AdimensionalSystemBC conditions(total_concentrations);
     conditions.charge_keeper = id_ho;
     conditions.add_fixed_fugacity_gas(id_co2g, id_hco3, -14);
 
     specmicp::AdimensionalSystemSolverOptions options;
     options.solver_options.maxstep = 50.0;
     options.solver_options.max_iter = 50;
     options.solver_options.maxiter_maxstep = 100;
     options.solver_options.use_crashing = false;
     options.solver_options.use_scaling = false;
     options.solver_options.factor_descent_condition = -1;
     options.solver_options.factor_gradient_search_direction = 100;
     options.solver_options.projection_min_variable = 1e-6;
     options.solver_options.fvectol = 1e-8;
     options.solver_options.steptol = 1e-10;
     options.solver_options.non_monotone_linesearch = false;
     options.system_options.non_ideality_tolerance = 1e-14;
 
     specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral);
     x.setZero();
     x(0) = 0.5;
     x.segment(1, raw_data->nb_component-1).setConstant(-6.0);
 
     for (specmicp::index_t i: raw_data->range_component())
     {
         std::cout << raw_data->labels_basis[i] << std::endl;
     }
 
     specmicp::uindex_t tot_iter = 0;
 
     specmicp::AdimensionalSystemSolver solver(raw_data, conditions, options);
     specmicp::micpsolver::MiCPPerformance perf =  solver.solve(x, false);
 
     specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
     //std::cout << perf.nb_iterations << std::endl;
     tot_iter += perf.nb_iterations;
 
     std::cout << "log fugacity" << "\t" << "pH";
     for (specmicp::index_t mineral: raw_data->range_mineral())
     {
         std::cout << "\t" << raw_data->labels_minerals[mineral];
     }
     std::cout << std::endl;
 
     specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
     specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set);
     std::cout << -14 << "\t" << sol.pH();
     for (specmicp::index_t mineral: raw_data->range_mineral())
     {
         std::cout << "\t" << sol.mole_concentration_mineral(mineral);
     }
     std::cout << std::endl;
 
     for (double logf=-13.9; logf<-4.1; logf+=0.05)
     {
         conditions.fixed_fugacity_bc[0].log_value = logf;
 
         solver = specmicp::AdimensionalSystemSolver (raw_data, conditions, solution, options);
         //std::cout << solver.get_options().solver_options.factor_descent_condition << std::endl;
         specmicp::micpsolver::MiCPPerformance perf =  solver.solve(x, false);
 
         specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
         std::cout << perf.nb_iterations << std::endl;
         tot_iter += perf.nb_iterations;
 
         solution = solver.get_raw_solution(x);
         specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set);
         std::cout << logf << "\t" << sol.pH();
         for (specmicp::index_t mineral: raw_data->range_mineral())
         {
             std::cout << "\t" << sol.mole_concentration_mineral(mineral);
         }
         std::cout << std::endl;
         //specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
         //std::cout << solution.main_variables(0) << std::endl;
         //std::cout << 14+solution.main_variables(1) << std::endl;
 
     }
     std::cout << tot_iter << std::endl;
 
 }
 
 
 
 void reaction_path_thermocarbo()
 {
     specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
     std::map<std::string, std::string> swapping ({
                                               {"H[+]","HO[-]"},
                                               {"Si(OH)4", "SiO(OH)3[-]"},
                                               {"Al[3+]","Al(OH)4[-]"},
                                                 });
     thedatabase.swap_components(swapping);
     thedatabase.remove_gas_phases();
     specmicp::RawDatabasePtr raw_data =     thedatabase.get_database();
 
     specmicp::Formulation formulation;
     specmicp::scalar_t mult = 5e3;
 
     specmicp::scalar_t m_c3s = mult*0.6;
     specmicp::scalar_t m_c2s = mult*0.2;
     specmicp::scalar_t m_c3a = mult*0.10;
     specmicp::scalar_t m_gypsum = mult*0.10;
     specmicp::scalar_t wc = 0.8;
     specmicp::scalar_t m_water = wc*1e-3*(
                   m_c3s*(3*56.08+60.08)
                 + m_c2s*(2*56.06+60.08)
                 + m_c3a*(3*56.08+101.96)
                 + m_gypsum*(56.08+80.06+2*18.02)
                 );
 
     formulation.mass_solution = m_water;
     formulation.amount_minerals = {
         {"C3S", m_c3s},
         {"C2S", m_c2s},
         {"C3A", m_c3a},
         {"Gypsum", m_gypsum}
     };
     formulation.extra_components_to_keep = {"HCO3[-]", };
     formulation.minerals_to_keep = {
         "Portlandite",
         "CSH,jennite",
         "CSH,tobermorite",
         "SiO2,am",
         "Calcite",
         "Al(OH)3,am",
         "Monosulfoaluminate",
         "Tricarboaluminate",
         "Monocarboaluminate",
         "Hemicarboaluminate",
         "Straetlingite",
         "Gypsum",
         "Ettringite",
         "Thaumasite"
     };
 
 
     specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
     specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
     specmicp::index_t id_ho  = thedatabase.component_label_to_id("HO[-]");
     specmicp::index_t id_hco3 = thedatabase.component_label_to_id("HCO3[-]");
     specmicp::index_t id_co2g = thedatabase.gas_label_to_id("CO2(g)");
     specmicp::index_t id_ca = thedatabase.component_label_to_id("Ca[2+]");
     std::cout << total_concentrations << std::endl;
     specmicp::AdimensionalSystemBC conditions(total_concentrations);
     conditions.charge_keeper = id_ho;
 
     specmicp::AdimensionalSystemSolverOptions options;
     options.solver_options.maxstep = 10.0;
     options.solver_options.max_iter = 100;
     options.solver_options.maxiter_maxstep = 100;
     options.solver_options.use_crashing = false;
     options.solver_options.use_scaling = false;
     options.solver_options.factor_descent_condition = -1;
     options.solver_options.factor_gradient_search_direction = 100;
     options.solver_options.projection_min_variable = 1e-9;
     options.solver_options.fvectol = 1e-6;
     options.solver_options.steptol = 1e-14;
     options.system_options.non_ideality_tolerance = 1e-10;
 
     specmicp::Vector x(raw_data->nb_component+raw_data->nb_mineral);
     x.setZero();
     x(0) = 0.8;
     x.segment(1, raw_data->nb_component-1).setConstant(-6.0);
 
     for (specmicp::index_t i: raw_data->range_component())
     {
         std::cout << raw_data->labels_basis[i] << std::endl;
     }
 
     specmicp::uindex_t tot_iter = 0;
 
     specmicp::AdimensionalSystemSolver solver(raw_data, conditions, options);
     specmicp::micpsolver::MiCPPerformance perf =  solver.solve(x, false);
 
     specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
     //std::cout << perf.nb_iterations << std::endl;
     tot_iter += perf.nb_iterations;
 
     std::cout << "Buffering \t" << "pH";
     for (specmicp::index_t mineral: raw_data->range_mineral())
     {
         std::cout << "\t" << raw_data->labels_minerals[mineral];
     }
     std::cout << std::endl;
 
     specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
     specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set);
     std::cout << 0 << "\t" << sol.pH();
     for (specmicp::index_t mineral: raw_data->range_mineral())
     {
         std::cout << "\t" << sol.mole_concentration_mineral(mineral);
     }
     std::cout << std::endl;
 
     for (double c_hco3=100; c_hco3<=conditions.total_concentrations(id_ca); c_hco3+=100)
     {
         conditions.total_concentrations(id_hco3) = c_hco3;
         conditions.total_concentrations(id_ho) -= 100;
         conditions.total_concentrations(id_h2o) += 100;
 
         solver = specmicp::AdimensionalSystemSolver (raw_data, conditions, solution, options);
         //std::cout << solver.get_options().solver_options.factor_descent_condition << std::endl;
         specmicp::micpsolver::MiCPPerformance perf =  solver.solve(x, false);
 
         specmicp_assert((int) perf.return_code > (int) specmicp::micpsolver::MiCPSolverReturnCode::NotConvergedYet);
         //std::cout << perf.nb_iterations << std::endl;
         tot_iter += perf.nb_iterations;
 
         solution = solver.get_raw_solution(x);
         specmicp::AdimensionalSystemSolutionExtractor sol(solution, raw_data, options.units_set);
         std::cout << c_hco3/conditions.total_concentrations(id_ca) << "\t" << sol.pH();
         for (specmicp::index_t mineral: raw_data->range_mineral())
         {
             std::cout << "\t" << sol.mole_concentration_mineral(mineral);
         }
         std::cout << std::endl;
         //specmicp::AdimensionalSystemSolution solution = solver.get_raw_solution(x);
         //std::cout << solution.main_variables(0) << std::endl;
         //std::cout << 14+solution.main_variables(1) << std::endl;
 
     }
     std::cout << tot_iter << std::endl;
 
 }
 
 class AutomaticReactionPath: public specmicp::EquilibriumCurve
 {
 public:
 
     AutomaticReactionPath()
     {
         specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
         std::map<std::string, std::string> swapping ({
                                                   {"H[+]","HO[-]"},
                                                   {"Si(OH)4", "SiO(OH)3[-]"},
                                                   {"Al[3+]","Al(OH)4[-]"},
                                                     });
         thedatabase.swap_components(swapping);
         thedatabase.remove_gas_phases();
         specmicp::RawDatabasePtr raw_data =     thedatabase.get_database();
         set_database(raw_data);
 
         specmicp::Formulation formulation;
         specmicp::scalar_t mult = 5e3;
 
         specmicp::scalar_t m_c3s = mult*0.6;
         specmicp::scalar_t m_c2s = mult*0.2;
         specmicp::scalar_t m_c3a = mult*0.10;
         specmicp::scalar_t m_gypsum = mult*0.10;
         specmicp::scalar_t wc = 0.8;
         specmicp::scalar_t m_water = wc*1e-3*(
                       m_c3s*(3*56.08+60.08)
                     + m_c2s*(2*56.06+60.08)
                     + m_c3a*(3*56.08+101.96)
                     + m_gypsum*(56.08+80.06+2*18.02)
                     );
 
         formulation.mass_solution = m_water;
         formulation.amount_minerals = {
             {"C3S", m_c3s},
             {"C2S", m_c2s},
             {"C3A", m_c3a},
             {"Gypsum", m_gypsum}
         };
         formulation.extra_components_to_keep = {"HCO3[-]", };
         formulation.minerals_to_keep = {
             "Portlandite",
             "CSH,jennite",
             "CSH,tobermorite",
             "SiO2,am",
             "Calcite",
             "Al(OH)3,am",
             "Monosulfoaluminate",
             "Tricarboaluminate",
             "Monocarboaluminate",
             "Hemicarboaluminate",
             "Straetlingite",
             "Gypsum",
             "Ettringite",
             "Thaumasite"
         };
 
 
         specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
         id_h2o = thedatabase.component_label_to_id("H2O");
         id_ho  = thedatabase.component_label_to_id("HO[-]");
         id_hco3 = thedatabase.component_label_to_id("HCO3[-]");
         id_co2g = thedatabase.gas_label_to_id("CO2(g)");
         id_ca = thedatabase.component_label_to_id("Ca[2+]");
 
         conditions() = specmicp::AdimensionalSystemBC(total_concentrations);
         conditions().charge_keeper = id_ho;
 
         specmicp::AdimensionalSystemSolverOptions& options = solver_options();
         options.solver_options.maxstep = 10.0;
         options.solver_options.max_iter = 100;
         options.solver_options.maxiter_maxstep = 100;
         options.solver_options.use_crashing = false;
         options.solver_options.use_scaling = false;
         options.solver_options.factor_descent_condition = -1;
         options.solver_options.factor_gradient_search_direction = 100;
         options.solver_options.projection_min_variable = 1e-9;
         options.solver_options.fvectol = 1e-6;
         options.solver_options.steptol = 1e-14;
         options.system_options.non_ideality_tolerance = 1e-10;
 
         solve_first_problem();
 
         std::cout << "Buffering \t" << "pH";
         for (specmicp::index_t mineral: raw_data->range_mineral())
         {
             std::cout << "\t" << raw_data->labels_minerals[mineral];
         }
         std::cout << std::endl;
 
         output();
 
     }
 
-    void output() const
+    void output()
     {
         specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set);
         std::cout <<  conditions().total_concentrations(id_hco3)/conditions().total_concentrations(id_ca) << "\t" << sol.pH();
         for (specmicp::index_t mineral: database()->range_mineral())
         {
             std::cout << "\t" << sol.mole_concentration_mineral(mineral);
         }
         std::cout << std::endl;
     }
 
     void update_problem()
     {
         conditions().total_concentrations(id_hco3) += 100;
         conditions().total_concentrations(id_ho) -= 100;
         conditions().total_concentrations(id_h2o) += 100;
     }
 
 private:
     specmicp::index_t id_h2o;
     specmicp::index_t id_ho;
     specmicp::index_t id_hco3;
     specmicp::index_t id_co2g;
     specmicp::index_t id_ca;
 };
 
 int main()
 {
     specmicp::logger::ErrFile::stream() = &std::cerr;
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
 
     //fixed_fugacity_thermocarbo();
     //reaction_path_thermocarbo();
 
     AutomaticReactionPath test_automatic;
 
     for (int i=0; i<130; ++i)
     {
         test_automatic.run_step();
     }
 
 
     return EXIT_SUCCESS;
 }
diff --git a/src/specmicp/adimensional/equilibrium_curve.hpp b/src/specmicp/adimensional/equilibrium_curve.hpp
index dd9d66a..854f039 100644
--- a/src/specmicp/adimensional/equilibrium_curve.hpp
+++ b/src/specmicp/adimensional/equilibrium_curve.hpp
@@ -1,102 +1,102 @@
 #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"
 
 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() const {}
+    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