diff --git a/src/specmicp/adimensional/adimensional_system_solver.cpp b/src/specmicp/adimensional/adimensional_system_solver.cpp
index cc3d463..3402e30 100644
--- a/src/specmicp/adimensional/adimensional_system_solver.cpp
+++ b/src/specmicp/adimensional/adimensional_system_solver.cpp
@@ -1,547 +1,551 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. 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:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. 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.
 
     3. Neither the name of the copyright holder 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 HOLDER 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 "adimensional_system_solver.hpp"
 #include "adimensional_system.hpp"
 #include "adimensional_system_solution.hpp"
 #include "adimensional_system_pcfm.hpp"
 
 #include "specmicp_common/micpsolver/micpsolver.hpp"
 
 #include "specmicp_common/log.hpp"
 
 #include <iostream>
 
 namespace specmicp {
 
 static constexpr int max_solver_restart = 3; // max number of time the solver is
                                              // restarted in case of failure
 
 
 // solve_equilibrium function
 // ==========================
 
 AdimensionalSystemSolution solve_equilibrium(
         std::shared_ptr<database::DataContainer> data,
         const AdimensionalSystemConstraints& constraints,
         const AdimensionalSystemSolverOptions& options
         )
 {
     AdimensionalSystemSolver solver(data, constraints, options);
     Vector variables;
 
     micpsolver::MiCPPerformance perf = solver.solve(variables, true);
 
     if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
         throw std::runtime_error("Failed to solve the problem");
 
     return solver.get_raw_solution(variables);
 }
 
 AdimensionalSystemSolution solve_equilibrium(
         RawDatabasePtr data_ptr,
         const AdimensionalSystemConstraints& constraints,
         const AdimensionalSystemSolution& previous_solution,
         const AdimensionalSystemSolverOptions& options
         )
 {
     AdimensionalSystemSolver solver(data_ptr, constraints, previous_solution, options);
     Vector variables;
 
     micpsolver::MiCPPerformance perf = solver.solve(variables, true);
 
     if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
         throw std::runtime_error("Failed to solve the problem");
 
     return solver.get_raw_solution(variables);
 }
 
 // The solver
 // ==========
 
 
 // constructor
 // -----------
 AdimensionalSystemSolver::AdimensionalSystemSolver(
         RawDatabasePtr data,
         const AdimensionalSystemConstraints& constraints
         ):
     OptionsHandler<AdimensionalSystemSolverOptions>(),
     m_data(data),
     m_system(std::make_shared<AdimensionalSystem>(
                  data,
                  constraints,
                  get_options().system_options,
                  get_options().units_set
                  )),
     m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral()))
 {}
 AdimensionalSystemSolver::AdimensionalSystemSolver(
         RawDatabasePtr data,
         const AdimensionalSystemConstraints& constraints,
         const AdimensionalSystemSolverOptions& options
         ):
     OptionsHandler<AdimensionalSystemSolverOptions>(options),
     m_data(data),
     m_system(std::make_shared<AdimensionalSystem>(
                  data,
                  constraints,
                  options.system_options,
                  options.units_set
                  )),
     m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral()))
 {}
 AdimensionalSystemSolver::AdimensionalSystemSolver(
         RawDatabasePtr data,
         const AdimensionalSystemConstraints& constraints,
         const AdimensionalSystemSolution& previous_solution
         ):
     OptionsHandler<AdimensionalSystemSolverOptions>(),
     m_data(data),
     m_system(std::make_shared<AdimensionalSystem>(
                  data,
                  constraints,
                  previous_solution,
                  get_options().system_options,
                  get_options().units_set
                  )),
     m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral()))
 {}
 AdimensionalSystemSolver::AdimensionalSystemSolver(
         RawDatabasePtr data,
         const AdimensionalSystemConstraints& constraints,
         const AdimensionalSystemSolution& previous_solution,
         const AdimensionalSystemSolverOptions& options
         ):
     OptionsHandler<AdimensionalSystemSolverOptions>(options),
     m_data(data),
     m_system(std::make_shared<AdimensionalSystem>(
                  data,
                  constraints,
                  previous_solution,
                  options.system_options,
                  options.units_set
                  )),
     m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral()))
 {}
 
 AdimensionalSystemSolution AdimensionalSystemSolver::get_raw_solution(
         Vector& x
         )
 {
     set_true_variable_vector(x);
     return m_system->get_solution(x, m_var);
 }
 
 
 AdimensionalSystemSolution AdimensionalSystemSolver::unsafe_get_raw_solution(
         Vector& x
         )
 {
     set_true_variable_vector(x);
     return m_system->unsafe_get_solution(x, m_var);
 }
 
 AdimensionalSystemSolution AdimensionalSystemSolver::merge_raw_solution_with_immobile_species(
         Vector& x,
         const AdimensionalSystemSolution& sol_immobile
         )
 {
     set_true_variable_vector(x);
     AdimensionalSystemSolution sol = m_system->unsafe_get_solution(x, m_var);
 
+    specmicp_assert(sol.main_variables.size() == sol_immobile.main_variables.size());
+    specmicp_assert(sol.sorbed_molalities.size() == sol_immobile.sorbed_molalities.size());
+    specmicp_assert(sol.log_gamma.size() == sol.log_gamma.size());
+
     sol.sorbed_molalities = sol_immobile.sorbed_molalities;
     if (m_data->nb_mineral() > 0)
     {
         auto m_start = m_system->dof_mineral(0);
         sol.main_variables.segment(m_start, m_data->nb_mineral()) =
             sol_immobile.main_variables.segment(m_start, m_data->nb_mineral());
     }
     if (m_data->nb_ssites() > 0)
     {
         auto s_start = m_system->dof_surface(0);
         sol.main_variables.segment(s_start, m_data->nb_ssites()) =
             sol_immobile.main_variables.segment(s_start, m_data->nb_ssites());
 
     }
     return sol;
 }
 
 
 // Solving the system
 // ------------------
 
 micpsolver::MiCPPerformance AdimensionalSystemSolver::solve(Vector& x, bool init)
 {
     // reset the units
     // - - - - - - - -
     m_system->set_units(get_options().units_set);
     // initialize the system
     // - - - - - - - - - - -
     if (init)
     {
         m_system->reasonable_starting_guess(x);
         if (get_options().use_pcfm)
         {
             run_pcfm(x);
         }
     }
     else if (get_options().force_pcfm)
     {
         run_pcfm(x);
     }
     // Solve the system
     // - - - - - - - -
     micpsolver::MiCPPerformance perf = solve_system(x);
     // If failed : try again
     // - - - - - - - - - - -
     int cnt = 0;
     while (perf.return_code < micpsolver::MiCPSolverReturnCode::Success
            and get_options().allow_restart
            and cnt < max_solver_restart
            )
     {
         WARNING << "Failed to solve the system ! Return code :"
                 << (int) perf.return_code
                 <<   ". We shake it up and start again";
         // Try safer options
         const scalar_t save_penalization_factor
                 = get_options().solver_options.penalization_factor;
         const bool save_scaling = get_options().solver_options.use_scaling;
         get_options().solver_options.use_scaling = true;
         if (save_penalization_factor == 1.0)
             get_options().solver_options.penalization_factor = 0.8;
         // Re-initialize the system
         set_return_vector(x);
         m_system->reasonable_restarting_guess(x);
         if (get_options().use_pcfm or get_options().force_pcfm)
             run_pcfm(x);
         // Solve the system
         micpsolver::MiCPPerformance perf2 = solve_system(x);
         // Restore options
         get_options().solver_options.penalization_factor = save_penalization_factor;
         get_options().solver_options.use_scaling = save_scaling;
         // Record the work performed
         perf += perf2;
         ++cnt;
     }
     // If solution is good
     // - - - - - - - - - -
     if (perf.return_code >= micpsolver::MiCPSolverReturnCode::Success)
     {
         set_return_vector(x);
     }
     return perf;
 }
 
 micpsolver::MiCPPerformance AdimensionalSystemSolver::solve_system(Vector& x)
 {
     set_true_variable_vector(x);
     if (not m_system->system_is_box_constrained()) {
         micpsolver::MiCPSolver<AdimensionalSystem, micpsolver::ReformulationF::CCK> solver(m_system);
         solver.set_options(get_options().solver_options);
         solver.solve(m_var);
         return solver.get_perfs();
     } else {
         micpsolver::MiCPSolver<AdimensionalSystem, micpsolver::ReformulationF::BoxVI> solver(m_system);
         solver.set_options(get_options().solver_options);
         solver.solve(m_var);
         return solver.get_perfs();
     }
 
 
 }
 
 // Variables management
 // ---------------------
 
 void AdimensionalSystemSolver::set_true_variable_vector(const Vector& x)
 {
     // This function sets the true variable for the specmicp system
     // the order is important
     // this is synchronised with the specmicp system, but needs to be set here
     const std::vector<index_t>& non_active_component = m_system->get_non_active_component();
     uindex_t new_i {0};
 
     if (m_system->is_active_component(0))
     {
         m_var(new_i) = x(m_system->dof_water());
         ++new_i;
     }
     for (index_t i: m_data->range_aqueous_component())
     {
         const auto it = std::find(non_active_component.cbegin(),
                                   non_active_component.cend(),
                                   i );
         if (it != non_active_component.cend()) continue;
         scalar_t value = x(m_system->dof_component(i));
         if (value == -HUGE_VAL) // check for previously undefined value
         {
             value = m_system->get_options().new_component_concentration;
         }
         m_var(new_i) = value;
         ++new_i;
     }
     if (m_system->surface_model() != SurfaceEquationType::NoEquation)
     {
         for (index_t q: m_data->range_ssites())
         {
             if (m_system->ideq_surf(q) != no_equation)
             {
                 m_var(new_i) = x(m_system->dof_surface(q));
                 ++new_i;
             }
         }
         if (m_system->ideq_surf_pot() != no_equation)
         {
             m_var(new_i) = x(m_system->dof_surface_potential());
             ++new_i;
         }
     }
     if (m_system->ideq_electron() != no_equation)
     {
         m_var(new_i) = x(m_system->dof_electron());
         ++new_i;
     }
     for (index_t m: m_data->range_mineral())
     {
         if (m_system->is_active_mineral(m))
         {
             m_var(new_i) = x(m_system->dof_mineral(m));
             ++new_i;
         }
     }
     m_var.conservativeResize(new_i);
     specmicp_assert(new_i == (unsigned) m_system->total_variables());
 }
 
 
 void AdimensionalSystemSolver::set_return_vector(Vector& x)
 {
     // This function sets the variable from the specmicp "raw" solution
     // the order is important
     // this is synchronised with the specmicp system, but needs to be set here
     const std::vector<index_t>& non_active_component = m_system->get_non_active_component();
 
     uindex_t new_i = 0;
     if (m_system->is_active_component(0))
     {
         x(m_system->dof_water()) = m_var(new_i);
         ++new_i;
     }
     else
     {
         x(m_system->dof_water()) = m_system->volume_fraction_water(x);
     }
     for (index_t i: m_data->range_aqueous_component())
     {
         const auto it = std::find(non_active_component.cbegin(), non_active_component.cend(),i);
         if (it != non_active_component.cend())
         {
             x(m_system->dof_component(i)) = -HUGE_VAL;
             continue;
         }
         x(m_system->dof_component(i)) = m_var(new_i) ;
         ++new_i;
     }
     for (index_t q: m_data->range_ssites())
     {
         if (m_system->ideq_surf(q) != no_equation)
         {
             x(m_system->dof_surface(q)) = m_var(new_i);
             ++new_i;
         } else {
             x(m_system->dof_surface(q)) = -HUGE_VAL;
         }
     }
     if (m_system->ideq_surf_pot() != no_equation)
     {
         x(m_system->dof_surface_potential()) = m_var(new_i);
         ++new_i;
     } else {
         x(m_system->dof_surface_potential()) = -HUGE_VAL;
     }
     if (m_system->ideq_electron() != no_equation)
     {
         x(m_system->dof_electron()) = m_var(new_i);
         ++new_i;
     }
     else
         x(m_system->dof_electron()) = -HUGE_VAL;
     for (index_t m: m_data->range_mineral())
     {
         if (m_system->is_active_mineral(m))
         {
             x(m_system->dof_mineral(m)) =m_var(new_i);
             ++new_i;
         }
         else
         {
             x(m_system->dof_mineral(m)) = 0.0;
         }
     }
 }
 
 // PCFM
 // ----
 
 void AdimensionalSystemSolver::run_pcfm(Vector &x)
 {
     DEBUG << "Start PCFM initialization.";
 
     // we set up the true variable
     set_true_variable_vector(x);
     // The residual is computed to have some point of comparison
     Vector residuals(m_system->total_variables());
     residuals.setZero();
     m_system->get_residuals(m_var, residuals);
     const scalar_t res_0 = residuals.norm();
 
     // the pcfm iterations are executed
     AdimensionalSystemPCFM pcfm_solver(m_system);
     PCFMReturnCode retcode = pcfm_solver.solve(m_var);
     // Check the answer
     if (retcode < PCFMReturnCode::Success)
     {
         // small prograss is still good enough
         m_system->get_residuals(m_var, residuals);
         const scalar_t final_residual = residuals.norm();
         DEBUG << "Final pcfm residuals : " << final_residual << " <? " << res_0;
         if (std::isfinite(final_residual) and final_residual < res_0)
             retcode = PCFMReturnCode::Success;
     }
     // if we did a good job, use the solution as initial state
     if (retcode == PCFMReturnCode::Success)
     {
         DEBUG << "Successful PCFM iterations.";
         set_return_vector(x);
     }
     else // we reset
     {
         DEBUG << "Unsuccessful PCFM iterations. The problem is reset to the initial state.";
         set_true_variable_vector(x);
         m_system->set_secondary_variables(m_var);
     }
 }
 
 // Initialisation of variables
 // ---------------------------
 
 void AdimensionalSystemSolver::initialize_variables(
         Vector& x,
         scalar_t volume_fraction_water,
         std::unordered_map<std::string, scalar_t> log_molalities,
         std::unordered_map<std::string, scalar_t> volume_fraction_minerals,
         scalar_t log_free_sorption_site_concentration
         )
 {
 
     m_system->reasonable_starting_guess(x);
 
     if (volume_fraction_water < 0 or volume_fraction_water > 1)
     {
         WARNING << "Initial guess for the volume fraction of water is not between 0 and 1";
     }
     x(m_system->dof_water()) = volume_fraction_water;
     for (auto pair: log_molalities)
     {
         index_t idc = m_data->get_id_component(pair.first);
         if (idc == no_species or idc == m_data->electron_index() or idc == m_data->water_index())
         {
             throw std::invalid_argument("This is not an aqueous component : "+pair.first);
         }
         if (pair.second > 0)
         {
             WARNING << "Initial molality for : " << pair.first << "is bigger than 1 molal.";
         }
         x(m_system->dof_component(idc)) = pair.second;
     }
     for (auto pair: volume_fraction_minerals)
     {
         index_t idm = m_data->get_id_mineral(pair.first);
         if (idm == no_species )
         {
             throw std::invalid_argument("This is not a mineral at equilibrium : "+pair.first);
         }
         if (pair.second < 0 or pair.second > 1)
         {
             WARNING << "Initial volume fraction for : " << pair.first << "is not between 0 and 1.";
         }
         x(m_system->dof_mineral(idm)) = pair.second;
     }
     if (log_free_sorption_site_concentration != 0.0) {
         for (index_t q: m_data->range_ssites()) {
             x(m_system->dof_surface(q)) = log_free_sorption_site_concentration;
         }
     }
 }
 
 void AdimensionalSystemSolver::initialize_variables(
         Vector& x,
         scalar_t volume_fraction_water,
         scalar_t log_molalities
         )
 {
     m_system->reasonable_starting_guess(x);
 
     if (volume_fraction_water < 0 or volume_fraction_water > 1)
     {
         WARNING << "Initial guess for the volume fraction of water is not between 0 and 1";
     }
     x(m_system->dof_water()) = volume_fraction_water;
 
     if (log_molalities > 0)
     {
         WARNING << "Initial molality for : " << log_molalities << "is bigger than 1 molal.";
     }
     x.segment(1, m_data->nb_component()-1).setConstant(log_molalities);
 }
 
 void AdimensionalSystemSolver::initialize_variables(
         Vector& x,
         scalar_t volume_fraction_water,
         scalar_t log_molalities,
         scalar_t log_free_sorption_site_concentration
         )
 {
     initialize_variables(x, volume_fraction_water, log_molalities);
     if (log_free_sorption_site_concentration != 0.0) {
         for (index_t q: m_data->range_ssites()) {
             x(m_system->dof_surface(q)) = log_free_sorption_site_concentration;
         }
     }
 }
 
 } // end namespace specmicp