diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9d36cf0..ea2ca69 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,192 +1,179 @@
 ##################### External packages #########################
 
 find_package(CxxTest)          # Necessary if you want the tests
 
 # catch for unit testing
 include(ExternalProject)
 find_package(Git REQUIRED)
 
 ExternalProject_Add(
     catch
     PREFIX ${CMAKE_BINARY_DIR}/catch
     GIT_REPOSITORY https://github.com/philsquared/Catch.git
     TIMEOUT 10
     UPDATE_COMMAND ${GIT_EXECUTABLE} pull
     CONFIGURE_COMMAND ""
     BUILD_COMMAND ""
     INSTALL_COMMAND ""
     LOG_DOWNLOAD ON
    )
 
 ################## Configuration ################################
 
 set(PROJECT_TEST_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 ExternalProject_Get_Property(catch source_dir)
 set(CATCH_INCLUDE_DIR ${source_dir}/include CACHE INTERNAL "Path to include folder for Catch")
 
 
 include_directories(${CATCH_INCLUDE_DIR})
 include_directories(${PROJECT_TEST_DIR})
 ################## Tests ########################################
 
 enable_testing(true)
 
 # CXX Test
 # ========
 
 if(CXXTEST_FOUND)
     set(CXXTEST_USE_PYTHON TRUE)
     include_directories(${CXXTEST_INCLUDE_DIR})
     enable_testing()
 
     # MiCP Solver
     # ------------
     CXXTEST_ADD_TEST(test_cond_number
                      test_cond_number.cpp
                      ${PROJECT_TEST_DIR}/micpsolver/test_cond_number.hpp)
 
     CXXTEST_ADD_TEST(test_ncp_funtion
                      test_ncp_function.cpp
                      ${PROJECT_TEST_DIR}/micpsolver/test_ncp_function.hpp)
 
     CXXTEST_ADD_TEST(test_micpsolver
                      test_micpsolver.cpp
                      ${PROJECT_TEST_DIR}/micpsolver/test_micpsolver.hpp)
 
     # Database
     # --------
 
     CXXTEST_ADD_TEST(test_data_reader
                     test_data_reader.cpp
                     ${PROJECT_TEST_DIR}/database/test_data_reader.hpp)
     target_link_libraries(test_data_reader specmicp_database jsoncpp)
 
 
     CXXTEST_ADD_TEST(test_data_selector
                     test_data_selector.cpp
                     ${PROJECT_TEST_DIR}/database/test_data_selector.hpp)
     target_link_libraries(test_data_selector specmicp_database jsoncpp)
 
     # Odeint
     # ------
 
     CXXTEST_ADD_TEST(test_butchertableau
                      test_butchertableau.cpp
                      ${PROJECT_TEST_DIR}/odeint/test_butchertableau.hpp)
 
 
     CXXTEST_ADD_TEST(test_embeddedrungekutta
                      test_embeddedrungekutta.cpp
                      ${PROJECT_TEST_DIR}/odeint/test_embeddedrungekutta.hpp)
 
 endif()
 
 # Test 'Perso'
 # ============
 
 add_executable(thermocarbo ${PROJECT_TEST_DIR}/specmicp/thermocarbo.cpp)
 target_link_libraries(thermocarbo specmicp specmicp_database)
 
 add_executable(carboalu ${PROJECT_TEST_DIR}/specmicp/carboalu.cpp)
 target_link_libraries(carboalu specmicp specmicp_database)
 
 add_executable(alkaliactivated ${PROJECT_TEST_DIR}/specmicp/alkaliactivated.cpp)
 target_link_libraries(alkaliactivated specmicp specmicp_database)
 
 add_executable(rate_nicoleau ${PROJECT_TEST_DIR}/specmicp/rate_nicoleau.cpp)
 target_link_libraries(rate_nicoleau specmicp specmicp_database)
 
 add_executable(kinetics ${PROJECT_TEST_DIR}/specmicp/kinetics.cpp)
 target_link_libraries(kinetics specmicp specmicp_database)
 
 
 set(DATA_TEST ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3sm.csv
               ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st1.csv
               ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st2.csv
 )
 
 add_custom_target(data_nicoleau SOURCES ${DATA_TEST})
 
 file(INSTALL ${DATA_TEST}
      DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data_test/)
 
 # Catch test
 #===========
 
-# Global step
-# -----------
-
-add_custom_target(catch_header SOURCES test_database.hpp)
-
-add_executable(test_reactmicp_diffusion
-    ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/test_reactmicp_diffusion.cpp
-    ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/diffusion.cpp
-    ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/numbering.cpp
-    ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/secondary.cpp
-    ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/solving.cpp
-    ${PROJECT_TEST_DIR}/reactmicp/systems/diffusion/utils.cpp
 )
-target_link_libraries(test_reactmicp_diffusion reactmicp specmicp specmicp_database)
 
 # Mesh
 # ----
 
 add_executable(test_meshes
     ${PROJECT_TEST_DIR}/reactmicp/meshes/test_meshes.cpp
     ${PROJECT_TEST_DIR}/reactmicp/meshes/test_uniform_mesh_1d.cpp
     ${PROJECT_TEST_DIR}/reactmicp/meshes/test_axisymmetric_uniform_mesh_1d.cpp
 )
 
 # Saturated Diffusion
 # -------------------
 
 add_executable(test_reactmicp_saturated_diffusion
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/test_reactmicp_saturated_diffusion.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/variables.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/solver.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/neutrality_solver.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/transport_program.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/reactive_transport.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_diffusion/utils.cpp
 )
 target_link_libraries(test_reactmicp_saturated_diffusion reactmicp specmicp specmicp_database)
 
 # Specmicp : Adim system
 # ----------------------
 
 add_executable(test_specmicp_adim
     ${PROJECT_TEST_DIR}/specmicp/adim/test_specmicp_adim.cpp
 
     ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system.cpp
     ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system_solver.cpp
     ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system_conditions.cpp
     ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system_problem_solver.cpp
     ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system_thermocarbo.cpp
     ${PROJECT_TEST_DIR}/specmicp/adim/adimensional_system_carboalu.cpp
 )
 
 target_link_libraries(test_specmicp_adim specmicp specmicp_database)
 
 # Reactmicp : Reactive transport solver
 # -------------------------------------
 
 add_executable(test_reactive_transport_solver
 
     ${PROJECT_TEST_DIR}/reactmicp/solver/test_reactive_transport_solver.cpp
     ${PROJECT_TEST_DIR}/reactmicp/solver/test_coupling.cpp
 )
 
 target_link_libraries(test_reactive_transport_solver reactmicp specmicp specmicp_database)
 
 # Saturated diffusion using new reactive transport solver
 # -------------------------------------------------------
 
 add_executable(test_reactmicp_saturated_react
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_react/test_reactmicp_saturated_react.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_react/speciation_system.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_react/variables.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_react/transport_program.cpp
     ${PROJECT_TEST_DIR}/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
 )
 
 target_link_libraries(test_reactmicp_saturated_react reactmicp specmicp specmicp_database)
diff --git a/tests/reactmicp/systems/diffusion/diffusion.cpp b/tests/reactmicp/systems/diffusion/diffusion.cpp
deleted file mode 100644
index 5205cd4..0000000
--- a/tests/reactmicp/systems/diffusion/diffusion.cpp
+++ /dev/null
@@ -1,63 +0,0 @@
-#include "catch.hpp"
-
-#include "utils.hpp"
-#include "reactmicp/systems/diffusion/diffusion.hpp"
-#include "reactmicp/meshes/mesh1d.hpp"
-
-#include "database/database.hpp"
-
-using namespace specmicp;
-using namespace specmicp::reactmicp;
-using namespace specmicp::reactmicp::systems;
-using namespace specmicp::reactmicp::systems::diffusion;
-
-TEST_CASE("Diffusion system", "[Diffusion]") {
-
-    DiffusionSecondaryVariables::Variables var;
-    std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
-    EquilibriumState comp_bc;
-    comp_bc = blank_composition(database);
-    BoundaryCondition bc_node0;
-    bc_node0.node = 0;
-    bc_node0.composition = comp_bc;
-    bc_node0.fix_all();
-    ListBoundaryCondition bcs = {bc_node0,};
-
-    EquilibriumState initial_state = sample_carbo_composition(database);
-
-    std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(4, 1e-3, 0.001);
-    std::shared_ptr<DiffusionParameter> parameters = std::make_shared<DiffusionParameter>(1e-8, 0.2);
-
-    SECTION("Creating a diffusion program") {
-        DiffusionProgram theprog(themesh, database, parameters, bcs, var);
-
-        REQUIRE(theprog.get_neq() == 4*(2*(database->nb_component-1)+database->nb_mineral));
-    }
-
-    SECTION("Initializing with initial state") {
-        DiffusionProgram theprog(themesh, database, parameters, bcs, initial_state, var);
-
-        REQUIRE(theprog.get_neq() == 4*(2*(database->nb_component-1)+database->nb_mineral));
-        REQUIRE(std::abs(std::pow(10.0, var.displacement((database->nb_component-1))) - comp_bc.molality_component(1)) < 1e-7);
-        REQUIRE(std::abs(std::pow(10.0,
-                                  var.displacement(3*(database->nb_component-1)+database->nb_mineral))
-                         - initial_state.molality_component(1)) < 1e-7);
-
-        Eigen::VectorXd residual(Eigen::VectorXd::Zero(theprog.get_neq()));
-        theprog.get_residuals(var, residual);
-
-        REQUIRE(residual.cols() == 1);
-        REQUIRE(residual.rows() == theprog.get_neq());
-        for (int it=0; it<theprog.get_neq(); ++it)
-        {
-            REQUIRE(std::isfinite(residual(it)));
-        }
-
-        SparseMatrix jacob;
-        theprog.get_jacobian(var, residual, jacob, 1);
-
-        REQUIRE(jacob.cols() == theprog.get_neq());
-        REQUIRE(jacob.rows() == theprog.get_neq());
-        REQUIRE(std::isfinite(jacob.norm()));
-    }
-}
diff --git a/tests/reactmicp/systems/diffusion/numbering.cpp b/tests/reactmicp/systems/diffusion/numbering.cpp
deleted file mode 100644
index d03a8d4..0000000
--- a/tests/reactmicp/systems/diffusion/numbering.cpp
+++ /dev/null
@@ -1,61 +0,0 @@
-#include "catch.hpp"
-#include "test_database.hpp"
-#include "reactmicp/systems/diffusion/diffusion_numbering.hpp"
-#include "utils.hpp"
-
-using namespace specmicp;
-using namespace specmicp::reactmicp::systems::diffusion;
-
-TEST_CASE("Correct numbering for the reactive diffusion 1d", "[Diffusion,Numbering]") {
-
-    ListBoundaryCondition void_bcs;
-    DiffusionNumbering::Variables var;
-    std::shared_ptr<specmicp::database::DataContainer> database = get_test_carbo_database();
-
-    SECTION("3 component and 5 mineral, 5 nodes, no bc") {
-        DiffusionNumbering numberer(5, 3, 5, void_bcs, var);
-
-        CHECK( numberer.get_ndf() == 11 );
-        CHECK( numberer.get_ndof_diffusion(1) == 0 );
-        CHECK( numberer.get_ndof_diffusion(3) == 2 );
-        CHECK( numberer.get_ndof_massbalance(1) == 3 );
-        CHECK( numberer.get_ndof_massbalance(2) == 4 );
-        CHECK( numberer.get_ndof_mineral(0) == 6 );
-        CHECK( numberer.get_ndof_mineral(4) == 10 );
-
-        CHECK( numberer.get_dof_diffusion(0, 3) == numberer.get_ndof_diffusion(3) );
-        CHECK( numberer.get_dof(3, 2) == 11*3+2 );
-
-        CHECK( numberer.get_neq() == 5*11 );
-        CHECK( numberer.id_equation_diffusion(3, 3)   == numberer.get_dof_diffusion(3, 3) );
-        CHECK( numberer.id_equation_massbalance(3, 3) == numberer.get_dof_massbalance(3, 3) );
-        CHECK( numberer.id_equation_mineral(3, 3)     == numberer.get_dof_mineral(3, 3) );
-    }
-
-    SECTION("Initialization from database") {
-
-        DiffusionNumbering numberer(5, database, void_bcs, var);
-
-        REQUIRE( numberer.get_ndf() == 4+4+database->nb_mineral );
-    }
-
-    SECTION("Boundary conditions") {
-        EquilibriumState comp_bc;
-        REQUIRE_NOTHROW(comp_bc = blank_composition(database));
-        BoundaryCondition bc_node0;
-        bc_node0.node = 0;
-        bc_node0.composition = comp_bc;
-        bc_node0.fix_all();
-        CHECK(bc_node0.mineral.size() == database->nb_mineral);
-        REQUIRE(bc_node0.component.size() == database->nb_component -1);
-
-        ListBoundaryCondition bcs = {bc_node0,};
-        DiffusionNumbering numberer(5, database, bcs, var);
-        REQUIRE(numberer.get_neq() == numberer.get_ndf()*4);
-        REQUIRE(numberer.node_has_bc(0) == true);
-        REQUIRE(numberer.node_has_bc(1) == false);
-        REQUIRE(std::abs(std::pow(10.0, numberer.component_concentration(0, 2, var)) - comp_bc.molality_component(2)) < 1e-7);
-
-
-    }
-}
diff --git a/tests/reactmicp/systems/diffusion/secondary.cpp b/tests/reactmicp/systems/diffusion/secondary.cpp
deleted file mode 100644
index 1497734..0000000
--- a/tests/reactmicp/systems/diffusion/secondary.cpp
+++ /dev/null
@@ -1,73 +0,0 @@
-#include "catch.hpp"
-
-#include "utils.hpp"
-#include "test_database.hpp"
-#include "reactmicp/systems/diffusion/diffusion_secondary.hpp"
-
-#include "database/database.hpp"
-
-using namespace specmicp;
-using namespace specmicp::reactmicp::systems::diffusion;
-
-TEST_CASE("Secondary variables", "[Secondary, Diffusion]") {
-
-    DiffusionSecondaryVariables::Variables var;
-    std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
-    EquilibriumState comp_bc;
-    comp_bc = blank_composition(database);
-    BoundaryCondition bc_node0;
-    bc_node0.node = 0;
-    bc_node0.composition = comp_bc;
-    bc_node0.fix_all();
-    ListBoundaryCondition bcs = {bc_node0,};
-    DiffusionNumbering numberer(5, database, bcs, var);
-
-    SECTION("Initialization") {
-        DiffusionSecondaryVariables secondary(5, database, numberer);
-
-        CHECK( secondary.secondary_concentration(2, 2) == 0 );
-        CHECK( secondary.secondary_concentration(2, 3) == 0 );
-        CHECK( secondary.secondary_concentration(2, 4) == 0 );
-        CHECK( secondary.loggamma_component(2, 3) == 0 );
-        CHECK( secondary.loggamma_secondary(2, 3) == 0 );
-    }
-
-    SECTION("Ionic strength computation") {
-        DiffusionSecondaryVariables secondary(5, database, numberer);
-
-        database::Database database_handler(database);
-
-        var.displacement(numberer.get_dof_massbalance(1, database_handler.component_label_to_id("Ca[2+]"))) = -3;
-        var.displacement(numberer.get_dof_massbalance(1, database_handler.component_label_to_id("HO[-]"))) = -3;
-        var.displacement(numberer.get_dof_massbalance(1, database_handler.component_label_to_id("HCO3[-]"))) = -3;
-        var.displacement(numberer.get_dof_massbalance(1, database_handler.component_label_to_id("SiO(OH)3[-]"))) = -3;
-
-        secondary.nodal_ionic_strength(1, var);
-        REQUIRE(secondary.ionic_strength(1) == (4*1e-3 +3*1e-3)*0.5);
-
-        secondary.secondary_concentration(1, database_handler.aqueous_label_to_id("CO3[2-]")) = 1e-2;
-        secondary.nodal_ionic_strength(1, var);
-        REQUIRE(secondary.ionic_strength(1) == (4*1e-3 +3*1e-3 + 4*1e-2)*0.5);
-
-        secondary.secondary_concentration(1, database_handler.aqueous_label_to_id("CO2")) = 1e-1;
-        secondary.nodal_ionic_strength(1, var);
-        REQUIRE(secondary.ionic_strength(1) == (4*1e-3 +3*1e-3 + 4*1e-2)*0.5);
-    }
-
-    SECTION("Aqueous concentration") {
-
-        DiffusionSecondaryVariables secondary(5, database, numberer);
-
-        secondary.ionic_strength(0) = comp_bc.ionic_strength();
-
-        secondary.nodal_loggamma(0);
-        secondary.nodal_secondary_concentrations(0, var);
-
-        for (int aq: database->range_aqueous())
-        {
-            REQUIRE(std::pow(10.0, secondary.loggamma_secondary(0, aq)) == comp_bc.activity_coefficient_secondary(aq));
-            REQUIRE(std::abs(secondary.secondary_concentration(0, aq) - comp_bc.molality_secondary()[aq]) < 1e-7);
-        }
-    }
-
-}
diff --git a/tests/reactmicp/systems/diffusion/solving.cpp b/tests/reactmicp/systems/diffusion/solving.cpp
deleted file mode 100644
index ba44475..0000000
--- a/tests/reactmicp/systems/diffusion/solving.cpp
+++ /dev/null
@@ -1,76 +0,0 @@
-#include "catch.hpp"
-
-#include "utils.hpp"
-#include "reactmicp/systems/diffusion/diffusion.hpp"
-#include "reactmicp/meshes/mesh1d.hpp"
-
-#include "database/database.hpp"
-
-#include "reactmicp/micpsolvers/micpsolver.hpp"
-#include "utils/log.hpp"
-#include <iostream>
-using namespace specmicp;
-using namespace specmicp::reactmicp;
-using namespace specmicp::reactmicp::systems;
-using namespace specmicp::reactmicp::systems::diffusion;
-
-TEST_CASE("Solving the diffusion system", "[Diffusion,MiCP]") {
-    specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
-    specmicp::logger::ErrFile::stream() = &std::cerr;
-
-    DiffusionSecondaryVariables::Variables var;
-    std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
-    EquilibriumState comp_bc;
-    comp_bc = blank_composition(database);
-    BoundaryCondition bc_node0;
-    bc_node0.node = 0;
-    bc_node0.composition = comp_bc;
-    bc_node0.fix_all();
-    ListBoundaryCondition bcs = {bc_node0,};
-
-    EquilibriumState initial_state = sample_carbo_composition(database);
-
-    std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(4, 1e-3, 0.001);
-    std::shared_ptr<DiffusionParameter> parameters = std::make_shared<DiffusionParameter>(1e-6, 0.2);
-
-    SECTION("Try solving") {
-        auto program = std::make_shared<DiffusionProgram>(themesh, database, parameters, bcs, initial_state, var);
-
-        specmicp::reactmicp::micpsolver::MiCPParabolicSolver<DiffusionProgram> solver(program);
-
-        //var.displacement(8+13) = 0.0005;
-        const int ndf = program->get_ndf();
-        for (int i=0; i<ndf; ++i)
-        {
-            for (int node: themesh->range_nodes())
-            {
-                std::cout << var.displacement(i+ndf*node) << "\t";
-            }
-            std::cout << std::endl;
-        }
-        //std::cout << var.displacement << std::endl;
-
-        for (int k=0; k<10; ++k)
-        {
-            auto retcode = solver.solve_one_timestep(1000, var);
-            std::cout << var.displacement(8+13) << std::endl;
-            for (int i=0; i<program->get_ndf(); ++i)
-            {
-                for (int node: themesh->range_nodes())
-                {
-                    std::cout << var.displacement(i+ndf*node) << "\t";
-                }
-                std::cout << std::endl;
-            }
-            REQUIRE((int) retcode == 1);
-        }
-        for (int i=0; i<ndf; ++i)
-        {
-            for (int node: themesh->range_nodes())
-            {
-                std::cout << var.displacement(i+ndf*node) << "\t";
-            }
-            std::cout << std::endl;
-        }
-    }
-}
diff --git a/tests/reactmicp/systems/diffusion/test_reactmicp_diffusion.cpp b/tests/reactmicp/systems/diffusion/test_reactmicp_diffusion.cpp
deleted file mode 100644
index 178916e..0000000
--- a/tests/reactmicp/systems/diffusion/test_reactmicp_diffusion.cpp
+++ /dev/null
@@ -1,3 +0,0 @@
-#define CATCH_CONFIG_MAIN
-#include "catch.hpp"
-
diff --git a/tests/reactmicp/systems/diffusion/utils.cpp b/tests/reactmicp/systems/diffusion/utils.cpp
deleted file mode 100644
index 3a1bb51..0000000
--- a/tests/reactmicp/systems/diffusion/utils.cpp
+++ /dev/null
@@ -1,104 +0,0 @@
-#include "utils.hpp"
-
-EquilibriumState blank_composition(std::shared_ptr<specmicp::database::DataContainer> database)
-{
-    std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
-    double m_c3s = 1e-7;
-    double m_c2s = 1e-7;
-    double m_water = 1;
-    double c_hco3 = 1e-4;
-    model->amount_aqueous = {
-        {"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0)+c_hco3, 0)},
-        {"HCO3[-]", specmicp::reaction_amount_t(c_hco3, 0)},
-        {"HO[-]", specmicp::reaction_amount_t(-c_hco3, 0)},
-    };
-    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 x;
-    specmicp::ReactionPathDriver driver(model, database, 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;
-
-    specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
-    if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
-        throw std::runtime_error("Error - system cannot be solved");
-    return driver.get_current_solution();
-
-}
-
-
-EquilibriumState sample_carbo_composition(std::shared_ptr<specmicp::database::DataContainer> database)
-{
-    std::shared_ptr<specmicp::ReactionPathModel> model = std::make_shared<specmicp::ReactionPathModel>();
-
-    const double mult = 0.005;
-    const double m_c3s = mult*0.7;
-    const double m_c2s = mult*0.3;
-    const double wc = 1.0;
-    const double m_water = wc*((3*56.08+60.08)*m_c3s + (2*56.08+60.08)*m_c2s)*1e-3;
-    const double delta_h2co3 = 0.005;
-    model->amount_aqueous = {
-        {"H2O", specmicp::reaction_amount_t(m_water/database->molar_mass_basis_si(0), delta_h2co3)},
-        {"HCO3[-]", specmicp::reaction_amount_t(delta_h2co3, 0.0)},
-        {"HO[-]", specmicp::reaction_amount_t(-delta_h2co3, 0.0)},
-    };
-    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 x;
-    specmicp::ReactionPathDriver driver(model, database, 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;
-
-    specmicp::micpsolver::MiCPPerformance perf =  driver.one_step(x);
-    if (perf.return_code != specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized)
-        throw std::runtime_error("Error - system cannot be solved");
-    return driver.get_current_solution();
-
-}
-
-//! \brief Return a simple fully processed database for testing purposes
-//!
-//! Correspond to the carbonation problem
-std::shared_ptr<specmicp::database::DataContainer> get_test_carbo_database()
-{
-    specmicp::database::Database database("data/cemdata_specmicp.js");
-    std::map<std::string, std::string> swapping ({
-                                              {"H[+]","HO[-]"},
-                                              {"Si(OH)4", "SiO(OH)3[-]"}
-                                                });
-    database.swap_components(swapping);
-
-    std::vector<std::string> to_keep = {"H2O", "HO[-]", "Ca[2+]", "SiO(OH)3[-]", "HCO3[-]"};
-    database.keep_only_components(to_keep);
-
-    return database.get_database();
-
-}
diff --git a/tests/reactmicp/systems/diffusion/utils.hpp b/tests/reactmicp/systems/diffusion/utils.hpp
deleted file mode 100644
index 8fde552..0000000
--- a/tests/reactmicp/systems/diffusion/utils.hpp
+++ /dev/null
@@ -1,9 +0,0 @@
-#include "specmicp/reaction_path.hpp"
-
-using namespace specmicp;
-
-EquilibriumState blank_composition(std::shared_ptr<specmicp::database::DataContainer> database);
-
-EquilibriumState sample_carbo_composition(std::shared_ptr<specmicp::database::DataContainer> database);
-
-std::shared_ptr<specmicp::database::DataContainer> get_test_carbo_database();