diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index bed9b7a..81276a0 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,25 +1,38 @@
 set(PROJECT_EXAMPLE_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 
-# ========== Reactmicp ====================================================
+
+# ========== SpecMiCP =====================================================
+
+# ---------- adimensional system ------------------------------------------
+
+#            thermocarbo
+
+add_executable(ex_adim_thermocarbo
+    ${PROJECT_EXAMPLE_DIR}/specmicp/adimensional/thermocarbo.cpp
+)
+
+target_link_libraries(ex_adim_thermocarbo specmicp specmicp_database)
+
+# ========== ReactMiCP ====================================================
 
 # ---------- leaching -----------------------------------------------------
 
 add_executable(leaching 
     ${PROJECT_EXAMPLE_DIR}/reactmicp/leaching.cpp
 )
 target_link_libraries(leaching reactmicp specmicp specmicp_database)
 
 # ---------- carbonation --------------------------------------------------
 
 add_executable(carbonation
     ${PROJECT_EXAMPLE_DIR}/reactmicp/carbonation.cpp
 )
 target_link_libraries(carbonation reactmicp specmicp specmicp_database)
 
 # ---------- diffusion ----------------------------------------------------
 
 add_executable(diffusion
     ${PROJECT_EXAMPLE_DIR}/reactmicp/diffusion.cpp
 )
 target_link_libraries(diffusion reactmicp specmicp specmicp_database)
 
diff --git a/examples/specmicp/adimensional/thermocarbo.cpp b/examples/specmicp/adimensional/thermocarbo.cpp
new file mode 100644
index 0000000..f4c161c
--- /dev/null
+++ b/examples/specmicp/adimensional/thermocarbo.cpp
@@ -0,0 +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
+    {
+        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;
+}