diff --git a/src/specmicp/reaction_path.cpp b/src/specmicp/reaction_path.cpp
index b269b37..2c5c1f9 100644
--- a/src/specmicp/reaction_path.cpp
+++ b/src/specmicp/reaction_path.cpp
@@ -1,157 +1,168 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : reaction_path.cpp
  - 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.
 
 ---------------------------------------------------------*/
 #include "reaction_path.hpp"
 
 #include "thermodata.hpp"
 #include "micpsolver/micpsolver.hpp"
 #include "utils/log.hpp"
 
 #include <iostream>
 
 namespace specmicp {
 
 void ReactionPathDriver::read_database()
 {
     m_database.parse_database(m_model->database_path);
     m_database.make_canonical();
     m_data = m_database.get_database();
 }
 
+//! \brief initialize
+void ReactionPathDriver::init_solution(Eigen::VectorXd& x)
+{
+    x.resize(m_data->nb_component+m_data->nb_mineral); // Set the right size
+    x.setZero();
+    m_must_init = true; // the solution vector must be initiated by the solver
+    // prepare the solution
+    dissolve_to_components();
+}
+
 void ReactionPathDriver::dissolve_to_components()
 {
     m_tot_conc = Eigen::VectorXd::Zero(m_data->nb_component);
     m_tot_conc_increase = Eigen::VectorXd::Zero(m_data->nb_component);
 
     // everything is dissolved ...
     for (auto it=m_model->amount_aqueous.begin(); it!=m_model->amount_aqueous.end(); ++it)
     {
         try
         {
             int id = database::Database::safe_label_to_id(it->first, m_data->labels_basis);
             m_tot_conc(id) += it->second.first;
             m_tot_conc_increase(id) += it->second.second;
         }
         catch(std::invalid_argument& e)
         {
             int id = database::Database::safe_label_to_id(it->first, m_data->labels_aqueous);
             for (int idc=0; idc<m_data->nb_component; ++idc)
             {
                 m_tot_conc(idc) += m_data->nu_aqueous(id, idc)*it->second.first;
                 m_tot_conc_increase(id) += m_data->nu_aqueous(id, idc)*it->second.second;
             }
         }
 
     }
     for (auto it=m_model->amount_minerals.begin(); it!=m_model->amount_minerals.end(); ++it)
     {
         try
         {
             int id = database::Database::safe_label_to_id(it->first, m_data->labels_minerals);
             for (int idc=0; idc<m_data->nb_component; ++idc)
             {
                 m_tot_conc(idc) += m_data->nu_mineral(id, idc)*it->second.first;
                 m_tot_conc_increase(idc) += m_data->nu_mineral(id, idc)*it->second.second;
             }
         }
         catch(std::invalid_argument& e)
         {
             int id = database::Database::safe_label_to_id(it->first, m_data->labels_minerals_kinetic);
             for (int idc=0; idc<m_data->nb_component; ++idc)
             {
                 m_tot_conc(idc) += m_data->nu_mineral_kinetic(id, idc)*it->second.first;
                 m_tot_conc_increase(idc) += m_data->nu_mineral_kinetic(id, idc)*it->second.second;
             }
         }
     }
     // Do we need to remove some components ?
     std::vector<int> to_remove;
     to_remove.reserve(10);
     int new_i = 0;
     for (int i=0; i<m_data->nb_component; ++i)
     {
         if (m_tot_conc(i) == 0 and m_tot_conc_increase(i) == 0)
         {
             to_remove.push_back(i);
             continue;
         }
         m_tot_conc(new_i) = m_tot_conc(i);
         m_tot_conc_increase(new_i) = m_tot_conc_increase(i);
         ++new_i;
     }
     m_tot_conc.conservativeResize(new_i);
     m_tot_conc_increase.conservativeResize(new_i);
     if (to_remove.size() > 0) m_database.remove_components(to_remove);
 
     if (m_model->minerals_to_keep.size() != 0)
     {
         m_database.minerals_keep_only(m_model->minerals_to_keep);
     }
 }
 
 micpsolver::MiCPPerformance  ReactionPathDriver::one_step(Eigen::VectorXd& x, bool init)
 {
     ReducedSystemSolver solver;
     if (m_current_solution.is_valid())
     {
         solver = ReducedSystemSolver(m_data, m_tot_conc, m_current_solution);
     }
     else
     {
         solver = ReducedSystemSolver(m_data, m_tot_conc);
     }
 
-
+    if (m_must_init == true) {init = true;}
     solver.get_options() = m_solver_options;
     micpsolver::MiCPPerformance perf = solver.solve(x, init);
     if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized)
     {
         ERROR << "Failed to solve the system ! (most probably)";
     }
     m_current_solution = solver.get_solution(x);
     m_tot_conc += m_tot_conc_increase;
+    m_must_init = false;
     return perf;
 }
 
 void ReactionPathDriver::run_all()
 {
     Eigen::VectorXd x;
     // first step
     one_step(x, true);
     for (int i=1; i<m_model->nb_step; ++i)
     {
         one_step(x);
     }
 }
 
 
 } // end namespace specmicp
diff --git a/src/specmicp/reaction_path.hpp b/src/specmicp/reaction_path.hpp
index a1a8160..eb564c0 100644
--- a/src/specmicp/reaction_path.hpp
+++ b/src/specmicp/reaction_path.hpp
@@ -1,135 +1,156 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : reaction_path.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_REACTIONPATH_HPP
 #define SPECMICP_SPECMICP_REACTIONPATH_HPP
 
 //! \file reaction_path.hpp reaction path driver
 
 #include <memory>
 #include <map>
 
 #include "database/database.hpp"
 #include "micpsolver/micpsolver_structs.hpp"
 #include "reduced_system_solver.hpp"
 #include "equilibrium_data.hpp"
 
 namespace specmicp {
 
 //! \brief typedef to a reaction path amount
 //!
 //! The first element is the inital amount,
 //! The second element is the increase to add at every step
 using reaction_amount_t = std::pair<double, double>;
 
 //! \brief A reaction path model
 //!
 //! Contains every data we need to run a reaction path model
 struct ReactionPathModel
 {
     //! The amount of components
     std::map<std::string, reaction_amount_t> amount_aqueous;
     //! The amount of minerals
     std::map<std::string, reaction_amount_t> amount_minerals;
 
     //! List of minerals to take into account
     std::vector<std::string> minerals_to_keep;
 
     //! the number of step
     int nb_step;
     //! the path to the database - optional
     std::string database_path;
 };
 
 //! \brief Driver to solve a reaction path model
 class ReactionPathDriver
 {
 public:
     //! \brief Initialise from a model
     //!
     //! Build a database from the attribute database_path of the model
     ReactionPathDriver(std::shared_ptr<ReactionPathModel> model,
                        ReducedSystemSolverOptions options=ReducedSystemSolverOptions()):
         m_model(model),
-        m_solver_options(options)
+        m_solver_options(options),
+        m_must_init(false)
     {
         read_database();
     }
 
     //! \brief Initialise the driver from a model and a raw database
     ReactionPathDriver(std::shared_ptr<ReactionPathModel> model,
                        std::shared_ptr<database::DataContainer> data,
                        ReducedSystemSolverOptions options=ReducedSystemSolverOptions()):
         m_model(model),
         m_database(data),
         m_data(data),
-        m_solver_options(options)
+        m_solver_options(options),
+        m_must_init(false)
     {}
 
+    //! \brief Initialise the driver from a model and a raw database
+    ReactionPathDriver(std::shared_ptr<ReactionPathModel> model,
+                       std::shared_ptr<database::DataContainer> data,
+                       Eigen::VectorXd& x,
+                       ReducedSystemSolverOptions options=ReducedSystemSolverOptions()):
+        m_model(model),
+        m_database(data),
+        m_data(data),
+        m_solver_options(options),
+        m_must_init(false)
+    {
+        init_solution(x);
+    }
+
+
     //! Initialize the database
     void read_database();
 
     //! Dissolve everything into components
     void dissolve_to_components();
 
     //! \brief Perform one step
     micpsolver::MiCPPerformance  one_step(Eigen::VectorXd& x, bool init=false);
 
     //! \brief run all steps
     void run_all();
 
     EquilibriumState get_current_solution() const {return m_current_solution;}
 
+    //! \brief initialize
+    void init_solution(Eigen::VectorXd& x);
 
     const ReducedSystemSolverOptions& get_options() const {return m_solver_options;}
     ReducedSystemSolverOptions& get_options() {return m_solver_options;}
 
 private:
     std::shared_ptr<ReactionPathModel> m_model;
 
     specmicp::database::Database m_database;
     std::shared_ptr<specmicp::database::DataContainer> m_data;
 
     Eigen::VectorXd m_tot_conc;
     Eigen::VectorXd m_tot_conc_increase;
 
     ReducedSystemSolverOptions m_solver_options;
 
     EquilibriumState m_current_solution;
 
+    bool m_must_init;
+
     //std::shared_ptr<ThermoData> m_thermo;
 };
 
 } // end namespace specmicp
 
 #endif // SPECMICP_SPECMICP_REACTIONPATH_HPP
diff --git a/tests/specmicp/alkaliactivated.cpp b/tests/specmicp/alkaliactivated.cpp
index e705b87..e23e622 100644
--- a/tests/specmicp/alkaliactivated.cpp
+++ b/tests/specmicp/alkaliactivated.cpp
@@ -1,131 +1,128 @@
 /*-------------------------------------------------------
 
  - Module : tests/
  - File : alkaliactivated.cpp
  - 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.
 
 ---------------------------------------------------------*/
 
 #include <iostream>
 #include "specmicp/reaction_path.hpp"
 
 #include "utils/log.hpp"
 
 void solve_lot_reaction_path()
 {
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
 
     specmicp::database::Database database("data/cemdata_specmicp.js");
     std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
     std::map<std::string, std::string> swapping ({
                                               {"H[+]","HO[-]"},
                                               {"Al[3+]","Al(OH)4[-]"},
                                               {"Si(OH)4", "SiO(OH)3[-]"}
                                                 });
     database.swap_components(swapping);
 
     std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
     double mult = 2;
     double m_naoh = 0.1;
     double m_sio2 = mult*0.5;
     double m_al2o3 = mult*0.5;
     double m_cao = mult*1.0;
     double m_gypsum = 0.2;
     //double wc = 0.5;
     double m_water = 1.0;
     //double delta_h2co3 = 0.1;
     double delta_cao = 0.1;
     model->amount_aqueous = {
         {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water -2*m_al2o3-m_cao-m_sio2, +delta_cao)},
         //{"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
         {"HO[-]", specmicp::reaction_amount_t(m_naoh - 2*m_al2o3 +2*m_cao-m_sio2, -2*delta_cao)},
         {"Na[+]", specmicp::reaction_amount_t(m_naoh, 0)},
         {"Al(OH)4[-]", specmicp::reaction_amount_t(2*m_al2o3, 0)},
         {"SiO(OH)3[-]", specmicp::reaction_amount_t(m_sio2, 0)},
         {"Ca[2+]", specmicp::reaction_amount_t(m_cao, -delta_cao)}
     };
     model->amount_minerals = {
         {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0)}
     };
     model->database_path = "data/cemdata_specmicp.js";
 
 
     double totiter = 0;
     double totfact = 0;
-    specmicp::ReactionPathDriver driver(model, data);
-    driver.dissolve_to_components();
+
+    Eigen::VectorXd x;
+    specmicp::ReactionPathDriver driver(model, data, x);
 
     Eigen::VectorXd totaq(data->nb_component);
-    Eigen::VectorXd x(data->nb_component+data->nb_mineral);
-    x(0) = 1.0;
-    x.block(1, 0, data->nb_component-1, 1).setConstant(-2);
-    x.block(data->nb_component, 0, data->nb_mineral-1, 1).setZero();
 
     //driver.get_options().solver_options.penalization_factor =1.0;
     driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
     //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
     driver.get_options().solver_options.use_scaling = false;
     driver.get_options().solver_options.max_factorization_step = 2;
     driver.get_options().solver_options.factor_gradient_search_direction = 50;
     driver.get_options().solver_options.maxstep = 50;
     driver.get_options().solver_options.maxiter_maxstep = 50;
     driver.get_options().solver_options.max_iter = 100;
     driver.get_options().allow_restart = true;
     std::cout << "Reaction_path  return_code  nb_iter  nb_fact   pH";
     for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     std::cout << std::endl;
     int max_step=11;
     for (int i=0; i<max_step; ++i)
     {
         specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
         specmicp::EquilibriumState solution = driver.get_current_solution();
         solution.total_aqueous_concentrations(totaq);
         std::cout << i*delta_cao << "  " <<  (int) perf.return_code << "   "
                   << perf.nb_iterations << "  " << perf.nb_factorization << "    "
                   << solution.pH() << "  "
                   << solution.mass_water() << " " << totaq.block(1,0,data->nb_component-1,1).transpose();
         for (int m=0; m < data->nb_mineral; ++m) std::cout << "   " << solution.moles_mineral(m);
         std::cout << std::endl;
         totiter += perf.nb_iterations;
         totfact += perf.nb_factorization;
     }
     std::cout << "Average iterations : " << totiter/max_step << std::endl;
     std::cout << "Average factorization : " << totfact/max_step << std::endl;
 }
 
 int main()
 {
     solve_lot_reaction_path();
 }
diff --git a/tests/specmicp/carboalu.cpp b/tests/specmicp/carboalu.cpp
index f27ceb6..ced0203 100644
--- a/tests/specmicp/carboalu.cpp
+++ b/tests/specmicp/carboalu.cpp
@@ -1,153 +1,148 @@
 /*-------------------------------------------------------
 
  - Module : tests/
  - File : carboalu.cpp
  - 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.
 
 ---------------------------------------------------------*/
 
 #include <iostream>
 #include "specmicp/reaction_path.hpp"
 #include "utils/log.hpp"
 
 double mult = 1.0;
 double m_c3s = mult*0.6;
 double m_c2s = mult*0.2;
 double m_c3a = mult*0.10;
 double m_gypsum = mult*0.10;
 double wc =0.8;
 
 void solve_carboalu()
 {
 
     specmicp::database::Database database("data/cemdata_specmicp.js");
     std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
     std::map<std::string, std::string> swapping ({
                                               {"H[+]","HO[-]"},
                                               {"Al[3+]","Al(OH)4[-]"},
                                               {"Si(OH)4", "SiO(OH)3[-]"}
                                                 });
     database.swap_components(swapping);
     std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
 
     double delta_h2co3 = 0.05;
     double delta_sio2 = 0.00;
     double 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));
     std::cout << m_water << std::endl;
     model->amount_aqueous = {
         {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water,
                 delta_h2co3+delta_sio2)},
         {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
         {"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3-delta_sio2)},
         {"SiO(OH)3[-]", specmicp::reaction_amount_t(0, delta_sio2)},
     };
     model->amount_minerals = {
         {"C3S", specmicp::reaction_amount_t(m_c3s, 0)},
         {"C2S", specmicp::reaction_amount_t(m_c2s, 0)},
         {"C3A", specmicp::reaction_amount_t(m_c3a, 0)},
         {"Gypsum", specmicp::reaction_amount_t(m_gypsum, 0)}
     };
 
     model->database_path = "data/cemdata_specmicp.js";
 
     model->minerals_to_keep = {
         "Portlandite",
         "CSH,jennite",
         "CSH,tobermorite",
         "SiO2,am",
         "Calcite",
         "Al(OH)3,am",
         "Monosulfoaluminate",
         "Tricarboaluminate",
         "Monocarboaluminate",
         "Hemicarboaluminate",
         "Straetlingite",
         "Gypsum",
         "Ettringite",
         "Thaumasite"
     };
 
-    specmicp::ReactionPathDriver driver(model, data);
-    driver.dissolve_to_components();
+    Eigen::VectorXd x;
+    specmicp::ReactionPathDriver driver(model, data, x);
 
     Eigen::VectorXd totaq(data->nb_component);
     double totiter = 0;
     double totfact = 0;
 
-    Eigen::VectorXd x(data->nb_component+data->nb_mineral);
-    x(0) = m_water;
-    x.block(1, 0, data->nb_component-1, 1).setConstant(-2);
-    x.block(data->nb_component, 0, data->nb_mineral, 1).setConstant(0.);
-
     //driver.get_options().solver_options.penalization_factor =0.8;
     driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
     //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
     driver.get_options().solver_options.use_scaling = true;
     driver.get_options().solver_options.max_factorization_step = 4;
     driver.get_options().solver_options.factor_gradient_search_direction = 100;
     driver.get_options().solver_options.max_iter = 50;
     driver.get_options().solver_options.maxstep = 50;
     driver.get_options().allow_restart = true;
     std::cout << "Reaction_path  return_code  nb_iter  nb_fact   pH";
     for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     std::cout << std::endl;
     const int nb_step = 51;
     for (int i=0; i<nb_step; ++i)
     {
         specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
         specmicp::EquilibriumState solution = driver.get_current_solution();
         solution.total_aqueous_concentrations(totaq);
         std::cout << i*(delta_sio2+delta_h2co3) << "  " <<  (int) perf.return_code << "   "
                   << perf.nb_iterations << "  " << perf.nb_factorization << "    "
                   << solution.pH() << "  "
                   << solution.mass_water() << " " << totaq.block(1,0,data->nb_component-1,1).transpose();
         for (int m=0; m < data->nb_mineral; ++m) std::cout << "   " << solution.moles_mineral(m);
         std::cout << std::endl;
 
         totiter += perf.nb_iterations;
         totfact += perf.nb_factorization;
     }
     std::cout << "Average iterations : " << totiter/nb_step << std::endl;
     std::cout << "Average factorization : " << totfact/nb_step << std::endl;
 }
 
 int main()
 {
 
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
     specmicp::logger::ErrFile::stream() = &std::cerr;
     for (int i=0; i<1; ++i)
     solve_carboalu();
 }
diff --git a/tests/specmicp/thermocarbo.cpp b/tests/specmicp/thermocarbo.cpp
index 7a5a784..416f57f 100644
--- a/tests/specmicp/thermocarbo.cpp
+++ b/tests/specmicp/thermocarbo.cpp
@@ -1,211 +1,206 @@
 /*-------------------------------------------------------
 
  - Module : tests/
  - File : thermocarbo.cpp
  - 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.
 
 ---------------------------------------------------------*/
 
 
 #include <iostream>
 #include "specmicp/reduced_system.hpp"
 #include "micpsolver/micpsolver.hpp"
 #include "micpsolver/micpsolver_min.hpp"
 #include "specmicp/reduced_system_solver.hpp"
 #include "database/database.hpp"
 
 #include "specmicp/reaction_path.hpp"
 
 void solve_lot_reaction_path()
 {
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
     specmicp::database::Database database("data/cemdata_specmicp.js");
     std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
     std::map<std::string, std::string> swapping ({
                                               {"H[+]","HO[-]"},
                                               {"Si(OH)4", "SiO(OH)3[-]"}
                                                 });
     database.swap_components(swapping);
 
     std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
     double m_c3s = 0.7;
     double m_c2s = 0.3;
     double wc = 0.5;
     double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3;
     double delta_h2co3 = 0.1;
     model->amount_aqueous = {
         {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, delta_h2co3)},
         {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
         {"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3)},
     };
     model->amount_minerals = {
         {"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
         {"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}
     };
 
     model->database_path = "data/cemdata_specmicp.js";
 
     Eigen::VectorXd totaq(5);
     double totiter = 0;
     double totfact = 0;
-    specmicp::ReactionPathDriver driver(model, data);
-    driver.dissolve_to_components();
-
 
-    Eigen::VectorXd x(data->nb_component+data->nb_mineral);
-    x(0) = 1.0;
-    x.block(1, 0, data->nb_component-1, 1).setConstant(-2.0);
-    x.block(data->nb_component, 0, data->nb_mineral, 1).setZero();
+    Eigen::VectorXd x;
+    specmicp::ReactionPathDriver driver(model, data, x);
 
     //driver.get_options().solver_options.penalization_factor =1.0;
     driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
     //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
     driver.get_options().solver_options.use_scaling = false;
     driver.get_options().solver_options.max_factorization_step = 1;
     driver.get_options().solver_options.factor_gradient_search_direction = 50;
     driver.get_options().solver_options.maxstep = 50;
     driver.get_options().solver_options.maxiter_maxstep = 50;
     driver.get_options().solver_options.max_iter = 100;
     driver.get_options().allow_restart = true;
     std::cout << "Reaction_path  return_code  nb_iter  nb_fact   pH";
     for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     std::cout << std::endl;
    int max_step=41;
     for (int i=0; i<max_step; ++i)
     {
         specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
         specmicp::EquilibriumState solution = driver.get_current_solution();
         solution.total_aqueous_concentrations(totaq);
         std::cout << i*(0.1) << "  " <<  (int) perf.return_code << "   "
                   << perf.nb_iterations << "  " << perf.nb_factorization << "    "
                   << solution.pH() << "  "
                   << solution.mass_water() << " " << totaq.block(1,0,data->nb_component-1,1).transpose();
                   //<< solution.molality_secondary(database.aqueous_label_to_id("CO2")) << "\t" ;
         for (int m=0; m < data->nb_mineral; ++m) std::cout << "   " << solution.moles_mineral(m);
         std::cout << std::endl;
         totiter += perf.nb_iterations;
         totfact += perf.nb_factorization;
     }
     std::cout << "Average iterations : " << totiter/max_step << std::endl;
     std::cout << "Average factorization : " << totfact/max_step << std::endl;
 
 }
 
 
 void solve()
 {
     specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
 
     specmicp::database::Database database("data/cemdata_specmicp.js");
     std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
     std::map<std::string, std::string> swapping ({
                                               {"H[+]","HO[-]"},
                                               {"Si(OH)4", "SiO(OH)3[-]"}
                                                 });
     database.swap_components(swapping);
 
     std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
     double m_c3s = 0.7;
     double m_c2s = 0.3;
     double wc = 1.0;
     double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3;
     double delta_h2co3 = 0.1;
     model->amount_aqueous = {
         {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, delta_h2co3)},
         {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
         {"HO[-]", specmicp::reaction_amount_t(0, -delta_h2co3)},
     };
     model->amount_minerals = {
         {"C3S", specmicp::reaction_amount_t(m_c3s, 0.0)},
         {"C2S", specmicp::reaction_amount_t(m_c2s, 0.0)}
     };
 
     //x(0) = m_water;
 
     model->database_path = "data/cemdata_specmicp.js";
 
     specmicp::ReactionPathDriver driver(model, data);
     driver.dissolve_to_components();
 
     Eigen::VectorXd totaq(data->nb_component);
     Eigen::VectorXd x(data->nb_component+data->nb_mineral);
     x(0) = 1.0;
     x.block(1, 0, data->nb_component-1, 1).setConstant(-2);
     x.block(data->nb_component, 0, data->nb_mineral, 1).setZero();
 
     driver.get_options().solver_options.penalization_factor = 0.8;
     driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
     //driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::min;
     driver.get_options().solver_options.use_scaling = false;
     driver.get_options().solver_options.max_factorization_step = 1;
     driver.get_options().solver_options.factor_gradient_search_direction = 50;
     driver.get_options().solver_options.maxstep = 50;
     driver.get_options().solver_options.maxiter_maxstep = 50;
     driver.get_options().solver_options.max_iter = 100;
     driver.get_options().allow_restart = false;
     driver.get_options().solver_options.power_descent_condition = 2.0;
     std::cout << "Reaction_path  return_code  nb_iter  nb_fact   pH";
     for (auto it = data->labels_basis.begin(); it!= data->labels_basis.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     std::cout << std::endl;
 
     specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
     specmicp::EquilibriumState solution = driver.get_current_solution();
     solution.total_aqueous_concentrations(totaq);
     std::cout << 0.0  << "  " <<  (int) perf.return_code << "   "
               << perf.nb_iterations << "  " << perf.nb_factorization << "    "
               << solution.pH() << "  "
               << solution.mass_water() << " " << totaq.block(1,0,data->nb_component-1,1).transpose();
     for (int m=0; m < data->nb_mineral; ++m) std::cout << "   " << solution.moles_mineral(m);
     std::cout << std::endl;
         std::cout << "Solution \n" << x << std::endl;
 
 }
 
 
 
 int main()
 {
     specmicp::logger::ErrFile::stream() = &std::cerr;
     //solve();
     solve_lot_reaction_path();
     //for (int i=0; i<5000; ++i) solve_lot_reaction_path();
 }