diff --git a/tests/specmicp/carboalu.cpp b/tests/specmicp/carboalu.cpp
index d5eb678..2978f93 100644
--- a/tests/specmicp/carboalu.cpp
+++ b/tests/specmicp/carboalu.cpp
@@ -1,90 +1,93 @@
 #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::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 delta_h2co3 = 0.1;
+    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_components = {
-        {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water- 4*m_c3s - 3*m_c2s -6*m_c3a+2*m_gypsum,
+    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)},
-        {"Al(OH)4[-]", specmicp::reaction_amount_t(2*m_c3a, 0)},
-        {"Ca[2+]", specmicp::reaction_amount_t(3*m_c3s +2*m_c2s+3*m_c3a+m_gypsum, 0)},
-        {"HO[-]", specmicp::reaction_amount_t(5*m_c3s + 3*m_c2s + 4* m_c3a, -delta_h2co3-delta_sio2)},
-        {"SiO(OH)3[-]", specmicp::reaction_amount_t(m_c3s + m_c2s, delta_sio2)},
-        {"SO4[2-]", specmicp::reaction_amount_t(m_gypsum, 0)}
+        {"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";
 
     specmicp::ReactionPathDriver driver(model, data);
     driver.dissolve_to_components();
 
     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 =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 = true;
+    //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 = 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 = false;
+    driver.get_options().allow_restart = true;
     std::cout << "Reaction_path  return_code  nb_iter  nb_fact   pH  H2O  HO  Al C Ca Si SO4";
     for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     std::cout << std::endl;
-    const int nb_step = 26;
+    const int nb_step = 51;
     for (int i=0; i<nb_step; ++i)
     {
         specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
         driver.total_aqueous_concentration(x, totaq);
         std::cout << i*(delta_sio2+delta_h2co3) << "  " <<  (int) perf.return_code << "   "
                   << perf.nb_iterations << "  " << perf.nb_factorization << "    "
                   << 14+x(1) << "  "
                   << x(0) << " " << totaq.block(1,0,data->nb_component-1,1).transpose() << "   "
             << x.block(data->nb_component, 0, data->nb_mineral, 1).transpose() << 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()
 {
-   for (int i=0; i<100; ++i)
+   for (int i=0; i<1; ++i)
     solve_carboalu();
 }
diff --git a/tests/specmicp/thermocarbo.cpp b/tests/specmicp/thermocarbo.cpp
index 2e8081b..5818209 100644
--- a/tests/specmicp/thermocarbo.cpp
+++ b/tests/specmicp/thermocarbo.cpp
@@ -1,588 +1,604 @@
 #include <iostream>
 #include "specmicp/extended_system.hpp"
 #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"
 
 std::shared_ptr<specmicp::ThermoData>  set_thermodata_from_database()
 {
     specmicp::database::Database database("data/cemdata_specmicp.js");
     std::shared_ptr<specmicp::database::DataContainer> data = database.get_database();
     //std::cout << data->nu_aqueous << std::endl;
     //std::cout << data->nu_mineral << std::endl;
 
     std::vector<int> new_basis({0, 5, 2, 3, 8});
     database.change_basis(new_basis);
 
     //std::cout << data->nu_aqueous << std::endl;
     //std::cout << data->nu_mineral << std::endl;
 
     std::cout << "Basis : " << std::endl;
     for (auto it = data->labels_basis.begin(); it!=data->labels_basis.end(); ++it)
     {
         std::cout << *it << "\t" << data->param_aq.row(it - data->labels_basis.begin()) << std::endl;
     }
     std::cout << std::endl;
     std::cout << "Aqueous : \n";
     for (auto it = data->labels_aqueous.begin(); it!=data->labels_aqueous.end(); ++it)
     {
         const int j = it - data->labels_aqueous.begin();
         std::cout << *it << "\t\t" << data->nu_aqueous.row(j) << "\t"
                   << data->logk_aqueous(j) << "\t" << data->param_aq.row(j+data->nb_component) << std::endl;
     }
     std::cout << "Minerals : \n";
     for (auto it = data->labels_minerals.begin(); it!=data->labels_minerals.end(); ++it)
     {
         const int j = it - data->labels_minerals.begin();
         std::cout << *it << "\t\t" << data->nu_mineral.row(j) << "\t" << data->logk_mineral(j) << std::endl;
     }
     std::cout << std::endl;
     std::shared_ptr<specmicp::ThermoData> ptr = std::make_shared<specmicp::ThermoData>(data);
     return ptr;
 
 }
 
 specmicp::micpsolver::MiCPPerformance solve_one(
         std::shared_ptr<specmicp::ThermoData> thermoptr,
         const Eigen::VectorXd& totconc,
                Eigen::VectorXd& x,
               Eigen::VectorXd& totaq)
 {
     x(15) =  std::max(0.0, x(15)-1.0);
     x(16) =  std::max(0.0, x(16)-1.0);
     x(17) =  std::max(0.0, x(17)-1.0);
     x(18) =  std::max(0.0, x(18)-1.0);
     x(19) =  std::max(0.0, x(19)-1.0);
     std::shared_ptr<specmicp::SpecProg> prog = std::make_shared<specmicp::SpecProg>(thermoptr, totconc);
     prog->init_secondary_conc(x);
     specmicp::micpsolver::MiCPSolver<specmicp::SpecProg, specmicp::micpsolver::NCPfunction::penalizedFB> solver(prog);
     //solver.get_options().fvectol = 1e-6;
     solver.get_options().maxstep  =  200.0;
     solver.get_options().maxiter_maxstep = 100;
     solver.get_options().use_crashing = false;
     solver.get_options().use_scaling = false;
     solver.get_options().projection_min_variable = 1e-6;
     solver.solve(x);
     prog->aqueous_tot_conc(x, totaq);
     return solver.get_performance();
 }
 
 void solve_lot()
 {
 
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
     Eigen::VectorXd x(20);
     x(0)  = 10.0;
     x(1)  = -2.0;
     x(2)  = -2.0;
     x(3)  = -2.0;
     x(4)  = -2.0;
     x(15) =  0.0;
     x(16) =  0.0;
     x(17) =  0.0;
     x(18) =  0.0;
     x(19) =  0.0;
     //x.setConstant(-1);
     //prog->init_secondary_conc(x);
     Eigen::VectorXd totaq(5);
 
     double m_c3s = 0.70;
     double m_c2s = 0.3;
     double m_h2co3 = 0.0;
 
     double totiter = 0;
     double totfact = 0;
     std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
     std::cout << "Reaction_path  return_code  nb_iter  pH  H2O  HO  C Ca Si "
               << " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
     for (int i=0; i<40; ++i)
     {
         m_h2co3 = 0.1*(1+i);
         Eigen::VectorXd totconc(5);
         totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
                 5*m_c3s + 3*m_c2s-m_h2co3,
                 m_h2co3,
                 3*m_c3s +2*m_c2s,
                 m_c3s + m_c2s;
         specmicp::micpsolver::MiCPPerformance perf = solve_one(thermoptr, totconc, x, totaq);
         std::cout << m_h2co3 << "  " <<  (int) perf.return_code << "   " << perf.nb_iterations << "  " << 14+x(1) << "  "
                   << x(0) << " " << totaq.block(1,0,4,1).transpose() << "   " << x.block(15, 0, 5, 1).transpose() << std::endl;
         totiter += perf.nb_iterations;
         totfact += perf.nb_factorization;
     }
     std::cout << "Average iterations : " << totiter/40 << std::endl;
     std::cout << "Average factorization : " << totfact/40 << std::endl;
 }
 
 specmicp::micpsolver::MiCPPerformance solve_one_reduced(
         std::shared_ptr<specmicp::database::DataContainer> data,
         Eigen::VectorXd& totconc,
         Eigen::VectorXd& x,
         Eigen::VectorXd& totaq)
 {
     x(5) =  std::max(0.0, x(5)-1.0);
     x(6) =  std::max(0.0, x(6)-1.0);
     x(7) =  std::max(0.0, x(7)-1.0);
     x(8) =  std::max(0.0, x(8)-1.0);
     x(9) =  std::max(0.0, x(9)-1.0);
     specmicp::micpsolver::MiCPSolverOptions options;
     options.maxstep  =  100.0;
     options.maxiter_maxstep = 100;
     options.use_crashing = false;
     options.use_scaling = true;
     options.projection_min_variable = 1e-6;
 
     specmicp::ReducedSystemSolver solver(data, totconc);
     solver.get_options().solver_options = options;
     specmicp::micpsolver::MiCPPerformance perf = solver.solve(x);
     solver.total_aqueous_concentration(x, totaq);
     return perf;
 }
 
 
 void solve_lot_reduced()
 {
 
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
     Eigen::VectorXd x(10);
     x(0)  = 1.0;
     x(1)  = -2.0;
     x(2)  = -2.0;
     x(3)  = -2.0;
     x(4)  = -2.0;
     x(5) =  0.0;
     x(6) =  0.0;
     x(7) =  0.0;
     x(8) =  0.0;
     x(9) =  0.0;
     //x.setConstant(-1);
     //prog->init_secondary_conc(x);
     Eigen::VectorXd totaq(5);
 
     double m_c3s = 0.70;
     double m_c2s = 0.3;
     double m_h2co3 = 0.0;
     double totiter = 0;
     double totfact = 0;
     std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
     Eigen::VectorXd totconc(5);
     totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
             5*m_c3s + 3*m_c2s-m_h2co3,
             m_h2co3,
             3*m_c3s +2*m_c2s,
             m_c3s + m_c2s;
     std::cout << "Reaction_path  return_code  nb_iter  pH  H2O  HO  C Ca Si "
               << " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
     for (int i=0; i<41; ++i)
     {
         m_h2co3 = 0.1*(i);
         totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
                 5*m_c3s + 3*m_c2s-m_h2co3,
                 m_h2co3,
                 3*m_c3s +2*m_c2s,
                 m_c3s + m_c2s;
         specmicp::micpsolver::MiCPPerformance perf = solve_one_reduced(thermoptr->get_raw_database(), totconc, x, totaq);
         //std::cout << x << std::endl;
         std::cout << m_h2co3 << "  " <<  (int) perf.return_code << "   " << perf.nb_iterations << "  " << 14+x(1) << "  "
                   << x(0) << " " << totaq.block(1,0,4,1).transpose() << "   " << x.block(5, 0, 5, 1).transpose() << std::endl;
         totiter += perf.nb_iterations;
         totfact += perf.nb_factorization;
     }
     std::cout << "Average iterations : " << totiter/41 << std::endl;
     std::cout << "Average factorization : " << totfact/41 << std::endl;
 }
 
 void solve_lot_reaction_path()
 {
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
     Eigen::VectorXd x(10);
     x(0) =  1.0;
-    x(1) = -2.0;
-    x(2) = -2.0;
-    x(3) = -2.0;
-    x(4) = -2.0;
-    x(5) =  0.0;
-    x(6) =  0.0;
-    x(7) =  0.0;
-    x(8) =  0.0;
-    x(9) =  0.0;
+    x(1) = -6.0;
+    x(2) = -6.0;
+    x(3) = -6.0;
+    x(4) = -6.0;
+    x(5) =  1.0;
+    x(6) =  1.0;
+    x(7) =  1.0;
+    x(8) =  1.0;
+    x(9) =  1.0;
     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 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_components = {
-        {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water- 4*m_c3s - 3*m_c2s, delta_h2co3)},
+    model->amount_aqueous = {
+        {"H2O", specmicp::reaction_amount_t(m_water/specmicp::molar_mass_water, delta_h2co3)},
         {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
-        {"Ca[2+]", specmicp::reaction_amount_t(3*m_c3s +2*m_c2s, 0)},
-        {"HO[-]", specmicp::reaction_amount_t(5*m_c3s + 3*m_c2s, -delta_h2co3)},
-        {"SiO(OH)3[-]", specmicp::reaction_amount_t(m_c3s + m_c2s, 0)}
+        {"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();
     //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 = 4;
+    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 = 100;
-    driver.get_options().solver_options.maxiter_maxstep = 100;
+    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().allow_restart = true;
     std::cout << "Reaction_path  return_code  nb_iter  nb_fact   pH  H2O  HO  C Ca Si "
                << " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
     int max_step=41;
     for (int i=0; i<max_step; ++i)
     {
-        x(0)  = 1.0;
-        x(1)  = -6.0;
-        x(2)  = -6.0;
-        x(3)  = -6.0;
-        x(4)  = -6.0;
-        x(5) =  1.0;
-        x(6) =  1.0;
-        x(7) =  1.0;
-        x(8) =  1.0;
-        x(9) =  1.0;
         specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
         driver.total_aqueous_concentration(x, totaq);
         std::cout << i*(0.1) << "  " <<  (int) perf.return_code << "   "
                   << perf.nb_iterations << "  " << perf.nb_factorization << "    "
                   << 14+x(1) << "  "
                   << x(0) << " " << totaq.block(1,0,4,1).transpose() << "   " << x.block(5, 0, 5, 1).transpose() << 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_lot_reaction_path_with_aluminum()
 {
     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)2[+]"},
                                               {"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.70;
     double m_c2s = 0.3;
     double m_aloh3 = 0.2;
     double m_caso4 = 0.0;
     double delta_h2co3 = 0.05;
     double delta_sio2 = 0.00;
-    model->amount_components = {
+    model->amount_aqueous = {
         {"H2O", specmicp::reaction_amount_t(1/specmicp::molar_mass_water- 4*m_c3s - 3*m_c2s, delta_h2co3+delta_sio2)},
         {"HCO3[-]", specmicp::reaction_amount_t(0, delta_h2co3)},
         {"Al(OH)2[+]", specmicp::reaction_amount_t(m_aloh3, 0)},
         {"Ca[2+]", specmicp::reaction_amount_t(3*m_c3s +2*m_c2s+m_caso4, 0)},
         {"HO[-]", specmicp::reaction_amount_t(5*m_c3s + 3*m_c2s + m_aloh3, -delta_h2co3-delta_sio2)},
         {"SiO(OH)3[-]", specmicp::reaction_amount_t(m_c3s + m_c2s, delta_sio2)},
         {"SO4[2-]", specmicp::reaction_amount_t(m_caso4, 0)}
     };
 
     model->database_path = "data/cemdata_specmicp.js";
 
     specmicp::ReactionPathDriver driver(model, data);
     driver.dissolve_to_components();
     Eigen::VectorXd totaq(data->nb_component);
     double totiter = 0;
     double totfact = 0;
 
     Eigen::VectorXd x(data->nb_component+data->nb_mineral);
     x(0) = 1.0;
     x.block(1, 0, data->nb_component-1, 1).setConstant(-6);
     x.block(data->nb_component, 0, data->nb_mineral, 1).setZero();
     //driver.get_options().solver_options.penalization_factor =0.5;
     driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
     driver.get_options().solver_options.use_scaling = true;
     driver.get_options().solver_options.max_factorization_step = 1;
     driver.get_options().solver_options.factor_gradient_search_direction = 100;
     driver.get_options().solver_options.max_iter = 45;
     std::cout << "Reaction_path  return_code  nb_iter  nb_fact   pH  H2O  HO  Al C Ca Si SO4";
     for (auto it = data->labels_minerals.begin(); it!= data->labels_minerals.end(); ++it)
     {
             std::cout << "   " << *it;
     }
     std::cout << std::endl;
     const int nb_step = 81;
     for (int i=0; i<nb_step; ++i)
     {
         specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
         driver.total_aqueous_concentration(x, totaq);
         std::cout << i*(delta_sio2) << "  " <<  (int) perf.return_code << "   "
                   << perf.nb_iterations << "  " << perf.nb_factorization << "    "
                   << 14+x(1) << "  "
                   << x(0) << " " << totaq.block(1,0,4,1).transpose() << "   "
             << x.block(data->nb_component, 0, data->nb_mineral, 1).transpose() << 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;
 }
 
 
 void solve_lot_reaction_path_scarce_water()
 {
     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)2[+]"},
                                               {"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.70;
     double m_c2s = 0.3;
     double m_aloh3 = 0.2;
     double m_caso4 = 0.1;
     double delta_h2o = -0.02/specmicp::molar_mass_water;
-    model->amount_components = {
+    model->amount_aqueous = {
         {"H2O", specmicp::reaction_amount_t(1/specmicp::molar_mass_water- 4*m_c3s - 3*m_c2s, delta_h2o)},
         {"HCO3[-]", specmicp::reaction_amount_t(0, 0)},
         {"Al(OH)2[+]", specmicp::reaction_amount_t(m_aloh3, 0)},
         {"Ca[2+]", specmicp::reaction_amount_t(3*m_c3s +2*m_c2s+m_caso4, 0)},
         {"HO[-]", specmicp::reaction_amount_t(5*m_c3s + 3*m_c2s + m_aloh3, 0)},
         {"SiO(OH)3[-]", specmicp::reaction_amount_t(m_c3s + m_c2s, 0)},
         {"SO4[2-]", specmicp::reaction_amount_t(m_caso4, 0)}
     };
 
     model->database_path = "data/cemdata_specmicp.js";
 
     specmicp::ReactionPathDriver driver(model, data);
     driver.dissolve_to_components();
     Eigen::VectorXd totaq(data->nb_component);
     double totiter = 0;
     double totfact = 0;
 
     Eigen::VectorXd x(data->nb_component+data->nb_mineral);
     x(0) = 1.0;
     x.block(1, 0, data->nb_component-1, 1).setConstant(-6);
     x.block(data->nb_component, 0, data->nb_mineral, 1).setZero();
     //driver.get_options().solver_options.penalization_factor =0.5;
     driver.get_options().ncp_function = specmicp::micpsolver::NCPfunction::penalizedFB;
     driver.get_options().solver_options.use_scaling = true;
     driver.get_options().solver_options.max_factorization_step = 1;
     driver.get_options().solver_options.factor_gradient_search_direction = 100;
     driver.get_options().solver_options.max_iter = 100;
     std::cout << "Reaction_path  return_code  nb_iter  nb_fact   pH  H2O  HO  Al Ca Si SO4";
     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);
         driver.total_aqueous_concentration(x, totaq);
         std::cout << (1-i*0.02)/specmicp::molar_mass_water << "  " <<  (int) perf.return_code << "   "
                   << perf.nb_iterations << "  " << perf.nb_factorization << "    "
                   << 14+x(1) << "  "
                   << x(0) << " " << totaq.block(1,0,5,1).transpose() << "   "
             << x.block(data->nb_component, 0, data->nb_mineral, 1).transpose() << 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;
 }
 
 
 void solve()
 {
     specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
-    std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
-    Eigen::VectorXd totconc(5);
-    double m_c3s = 0.70;
+    Eigen::VectorXd x(10);
+    x(0) =  1;
+    x(1) =  -2;
+    x(2) = -2;
+    x(3) = -2;
+    x(4) = -2;
+    x(5) =  0.0;
+    x(6) =  0.0;
+    x(7) =  0.0;
+    x(8) =  0.0;
+    x(9) =  0.0;
+    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 m_h2co3 = 1.9;
-    totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
-            5*m_c3s + 3*m_c2s-m_h2co3,
-            m_h2co3,
-            3*m_c3s +2*m_c2s,
-            m_c3s + m_c2s;
-    std::shared_ptr<specmicp::SpecProg> prog = std::make_shared<specmicp::SpecProg>(thermoptr, totconc);
-    specmicp::micpsolver::MiCPSolver<specmicp::SpecProg, specmicp::micpsolver::NCPfunction::penalizedFB> solver(prog);
-    //solver.get_options().fvectol = 1e-6;
-    solver.get_options().max_iter = 100;
-    solver.get_options().maxstep  =  200.0;
-    solver.get_options().maxiter_maxstep = 100;
-    solver.get_options().factor_descent_condition = 1e-4;
-    solver.get_options().use_crashing = 0;
-    Eigen::VectorXd x(prog->total_variables());
-    //x.setConstant(-1);
-    x(0)  = 1.0;
-    x(1)  = -2.0;
-    x(2)  = -2.0;
-    x(3)  = -2.0;
-    x(4)  = -2.0;
-    x(15) =  0.0;
-    x(16) =  0.0;
-    x(17) =  0.0;
-    x(18) =  0.0;
-    x(19) =  0.0;
-    prog->init_secondary_conc(x);
+    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)}
+    };
 
-    std::cout << x << std::endl;
+    //x(0) = m_water;
 
-//    Eigen::VectorXd residual(20);
-//    prog->get_residuals(x, residual);
-//    std::cout << "Residual \n ------ \n" << residual << std::endl;
-//    Eigen::MatrixXd jacob(20, 20);
-//    prog->get_jacobian(x, jacob);
-//    std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
+    model->database_path = "data/cemdata_specmicp.js";
 
-    std::cout << "Return Code : " <<  (int) solver.solve(x) << std::endl;
+    Eigen::VectorXd totaq(5);
+    specmicp::ReactionPathDriver driver(model, data);
+    driver.dissolve_to_components();
+    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  H2O  HO  C Ca Si "
+               << " Portlandite SiO2(am) Jennite Tobermorite Calcite" << std::endl;
 
+        specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
+        driver.total_aqueous_concentration(x, totaq);
+        std::cout << 0*(0.1) << "  " <<  (int) perf.return_code << "   "
+                  << perf.nb_iterations << "  " << perf.nb_factorization << "    "
+                  << 14+x(1) << "  "
+                  << x(0) << " " << totaq.block(1,0,4,1).transpose() << "   " << x.block(5, 0, 5, 1).transpose() << std::endl;
+
+        std::cout << "Solution \n" << x << std::endl;
 
-    std::cout << "Solution \n -------- \n" << x << std::endl;
 }
 
 void solve_reduced()
 {
 
     specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
     std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
     Eigen::VectorXd totconc(5);
     double m_c3s = 0.70;
     double m_c2s = 0.3;
     double m_h2co3 = 2.8;
     totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
             5*m_c3s + 3*m_c2s-m_h2co3,
             m_h2co3,
             3*m_c3s +2*m_c2s,
             m_c3s + m_c2s;
     std::shared_ptr<specmicp::ReducedSystem> prog = std::make_shared<specmicp::ReducedSystem>(thermoptr, totconc);
     specmicp::micpsolver::MiCPSolver<specmicp::ReducedSystem, specmicp::micpsolver::NCPfunction::penalizedFB> solver(prog);
     //solver.get_options().fvectol = 1e-6;
     solver.get_options().max_iter = 100;
     solver.get_options().maxstep  =  200;
     solver.get_options().maxiter_maxstep = 100;
     solver.get_options().factor_descent_condition = 1e-4;
     solver.get_options().use_crashing = false;
     solver.get_options().use_scaling = true;
     Eigen::VectorXd x(prog->total_variables());
     //x.setConstant(-1);
     x(0)  = 1.0485;
     x(1)  = -5.01517;
     x(2)  = -3.52536;
     x(3)  = -3.50741;
     x(4)  = -3.5383;
     x(5) =  0.0;
     x(6) =  0.997242;
     x(7) =  0.0;
     x(8) =  0.0;
     x(9) =  2.69967;
 
     std::cout << x << std::endl;
 
 //    Eigen::VectorXd residual(20);
 //    prog->get_residuals(x, residual);
 //    std::cout << "Residual \n ------ \n" << residual << std::endl;
 //    Eigen::MatrixXd jacob(20, 20);
 //    prog->get_jacobian(x, jacob);
 //    std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
 
     std::cout << "Return Code : " <<  (int) solver.solve(x) << std::endl;
 
 
     std::cout << "Solution \n -------- \n" << x << std::endl;
 }
 
 
 void solve_reduced_min()
 {
 
     specmicp::stdlog::ReportLevel() = specmicp::logger::Debug;
     std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_from_database();
     Eigen::VectorXd totconc(5);
     double m_c3s = 0.70;
     double m_c2s = 0.3;
     double m_h2co3 = 2.8;
     totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
             5*m_c3s + 3*m_c2s-m_h2co3,
             m_h2co3,
             3*m_c3s +2*m_c2s,
             m_c3s + m_c2s;
     std::shared_ptr<specmicp::ReducedSystem> prog = std::make_shared<specmicp::ReducedSystem>(thermoptr, totconc);
     specmicp::micpsolver::MiCPSolver<specmicp::ReducedSystem, specmicp::micpsolver::NCPfunction::min> solver(prog);
     //solver.get_options().fvectol = 1e-6;
     solver.get_options().max_iter = 100;
     solver.get_options().maxstep  =  200;
     solver.get_options().maxiter_maxstep = 100;
     solver.get_options().factor_descent_condition = 1e-4;
     solver.get_options().use_crashing = false;
     solver.get_options().use_scaling = true;
     Eigen::VectorXd x(prog->total_variables());
     //x.setConstant(-1);
     x(0)  = 1.0;
     x(1)  = -2;
     x(2)  = -2;
     x(3)  = -2;
     x(4)  = -2;
     x(5) =  0.0;
     x(6) =  0.0;
     x(7) =  0.0;
     x(8) =  0.0;
     x(9) =  0.0;
 
     std::cout << x << std::endl;
 
 //    Eigen::VectorXd residual(20);
 //    prog->get_residuals(x, residual);
 //    std::cout << "Residual \n ------ \n" << residual << std::endl;
 //    Eigen::MatrixXd jacob(20, 20);
 //    prog->get_jacobian(x, jacob);
 //    std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
 
     std::cout << "Return Code : " <<  (int) solver.solve(x) << std::endl;
 
 
     std::cout << "Solution \n -------- \n" << x << std::endl;
 }
 
+
+
+
 int main()
 {
     specmicp::logger::ErrFile::stream() = &std::cerr;
     //set_thermodata_from_database();
     //solve();
     //solve_lot();
     //solve_reduced();
     //solve_reduced_min();
     //solve_lot_reduced();
     //solve_aluminium();
     solve_lot_reaction_path();
-    //for (int i=0; i<1000; ++i) solve_lot_reaction_path();
+    //for (int i=0; i<5000; ++i) solve_lot_reaction_path();
     //solve_lot_reaction_path_with_aluminum();
     //solve_lot_reaction_path_scarce_water();
 }