diff --git a/src/dfpmsolver/driver_structs.hpp b/src/dfpmsolver/driver_structs.hpp
index 60df068..0c369f3 100644
--- a/src/dfpmsolver/driver_structs.hpp
+++ b/src/dfpmsolver/driver_structs.hpp
@@ -1,63 +1,67 @@
 #ifndef SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
 #define SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
 
 #include "common.hpp"
 #include "utils/sparse_solvers/sparse_solver_structs.hpp"
 
 namespace specmicp {
 namespace dfpmsolver {
 
 
 //! \brief Options of a driver
 //!
 struct DriverOptions
 {
     scalar_t residuals_tolerance;    //!< Tolerance for the residual
+    scalar_t absolute_tolerance;    //!< Absolute tolerance for the residual
     scalar_t step_tolerance;        //!< Tolerance for the minimum step length
     scalar_t threshold_stationary_point; //!< if ||R||>threshold, the point is classified as stationary
     int maximum_iterations;       //!< Maximum iterations allowed
     scalar_t maximum_step_length;   //!< Maximum step length allowed
     int max_iterations_at_max_length;   //!< Maximum number of iterations at maximum step length
     scalar_t coeff_accept_newton_step; //!< Accept Newton step if enough progress is made
     sparse_solvers::SparseSolver sparse_solver;   //!< The sparse solver to use
 
     DriverOptions():
         residuals_tolerance(5e-5),
+        absolute_tolerance(1e-12),
         step_tolerance(1e-10),
         threshold_stationary_point(1e-4),
         maximum_iterations(200),
         maximum_step_length(1e3),
         max_iterations_at_max_length(50),
         coeff_accept_newton_step(0.9),
         sparse_solver(sparse_solvers::SparseSolver::SparseQR)
     {}
 };
 
 //! \brief Performance of a driver
 struct DriverPerformance
 {
     int nb_call_residuals;
     int nb_call_jacobian;
     int nb_iterations;
     int nb_consecutive_max_step_taken;
     int nb_max_step_taken;
     bool maximum_step_taken;
     scalar_t current_residual;
+    scalar_t absolute_residual;
     scalar_t current_update;
 
     DriverPerformance():
         nb_call_residuals(0),
         nb_call_jacobian(0),
         nb_iterations(0),
         nb_consecutive_max_step_taken(0),
         nb_max_step_taken(0),
         maximum_step_taken(0),
         current_residual(0.0),
+        absolute_residual(0.0),
         current_update(0.0)
     {}
 };
 
 } // end namespace dfpmsolver
 } // end namespace specmicp
 
 #endif // SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
diff --git a/src/dfpmsolver/parabolic_driver.inl b/src/dfpmsolver/parabolic_driver.inl
index 28a5e52..6440f61 100644
--- a/src/dfpmsolver/parabolic_driver.inl
+++ b/src/dfpmsolver/parabolic_driver.inl
@@ -1,476 +1,481 @@
 #include "parabolic_driver.hpp" // for syntaxic coloration only
 #include "utils/log.hpp"
 //#include "utils/sparse_solver.hpp"
 
 namespace specmicp {
 namespace dfpmsolver {
 
 template <class Program>
 void ParabolicDriver<Program>::compute_jacobian(Vector& displacement,
                       Vector& velocity,
                       Eigen::SparseMatrix<scalar_t>& jacobian
                       )
 {
     program().compute_jacobian(displacement, velocity, jacobian, get_options().alpha*m_current_dt);
     jacobian = jacobian*scaling().asDiagonal();
     jacobian.makeCompressed();
 
     get_perfs().nb_call_jacobian += 1;
 }
 
 template <class Program>
 ParabolicDriverReturnCode  ParabolicDriver<Program>::solve_timestep(scalar_t dt, Eigen::VectorXd& displacement)
 {
     initialize_timestep(dt, displacement);
     return restart_timestep(displacement);
 }
 
 template <class Program>
 void  ParabolicDriver<Program>::initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement)
 {
     initialize_scaling();
     m_residuals = Eigen::VectorXd::Zero(get_neq());
     m_current_dt = dt;
     set_predictor(displacement);
     reset_velocity();
     program().apply_bc(dt, displacement, m_velocity);
 
     compute_residuals(displacement, m_velocity, m_residuals);
     m_norm_0 = norm_residuals();
     get_perfs().nb_iterations = 0;
 }
 
 template <class Program>
 void  ParabolicDriver<Program>::set_predictor(Vector& displacement)
 {
     if (get_options().alpha < 1)
         m_predictor = displacement + (1-get_options().alpha)*m_current_dt*m_velocity;
     else
         m_predictor = displacement;
 }
 
 template <class Program>
 scalar_t ParabolicDriver<Program>::update_norm(const Vector& update)
 {
     // l-∞ scaled norm
     scalar_t norm = 0.0;
     for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
     {
         const index_t id_eq = program().id_equation(dof);
         if (id_eq == no_species or update(id_eq) == 0.0) continue;
         norm = std::max(norm, std::abs(update(id_eq))/(std::max(std::abs(m_velocity(dof)), 1.0/scaling(id_eq))));
     }
     return norm;
 }
 
 template <class Program>
 ParabolicDriverReturnCode  ParabolicDriver<Program>::restart_timestep(Vector& displacement)
 {
     ParabolicDriverReturnCode return_code = ParabolicDriverReturnCode::NotConvergedYet;
 
     //m_solver.reset(nullptr);
 
     Eigen::VectorXd update(get_neq());
     update.setZero();
     get_perfs().current_update = 0;
     bool force_recompute_jacobian = true;
     while (return_code == ParabolicDriverReturnCode::NotConvergedYet)
     {
         compute_residuals(displacement, m_velocity, m_residuals);
+        get_perfs().absolute_residual = m_residuals.norm();
         get_perfs().current_residual = m_residuals.norm()/m_norm_0;
         get_perfs().current_update = update_norm(update);
         DEBUG << " NB iterations : " << get_perfs().nb_iterations
                   << " - res : " << get_perfs().current_residual
                   << " - update : " << get_perfs().current_update;
         return_code = check_convergence();
 
         if (return_code != ParabolicDriverReturnCode::NotConvergedYet) break;
         if (m_solver == nullptr)
         {
             m_solver = sparse_solvers::get_sparse_solver<
                     Eigen::SparseMatrix<scalar_t>, Vector, Vector>(get_options().sparse_solver);
             compute_jacobian(displacement, m_velocity, m_jacobian);
             m_solver->analyse_pattern(m_jacobian);
             m_solver->decompose(m_jacobian);
             force_recompute_jacobian =false;
         }
         else if (force_recompute_jacobian
                 or get_perfs().nb_iterations % get_options().quasi_newton == 0)
         {
             compute_jacobian(displacement, m_velocity, m_jacobian);
             m_solver->decompose(m_jacobian);
         }
         get_perfs().nb_iterations += 1;
         m_gradient = m_jacobian.transpose()*m_residuals;
         sparse_solvers::SparseSolverReturnCode retcode = m_solver->solve_scaling(m_residuals, scaling(), update);
         if (retcode != sparse_solvers::SparseSolverReturnCode::Success)
         {
             ERROR << "Error when solving linear system : " << (int) retcode << std::endl;
             return ParabolicDriverReturnCode::ErrorLinearSystem;
         }
         //if (update.norm() < get_options().step_tolerance) return_code = ParabolicDriverReturnCode::ErrorMinimized;
         //else
             return_code = linesearch(update, displacement);
     }
     return return_code;
 
 }
 
 template <class Program>
 ParabolicDriverReturnCode  ParabolicDriver<Program>::check_convergence()
 {
     ParabolicDriverReturnCode termcode = ParabolicDriverReturnCode::NotConvergedYet;
     const int nb_iterations = get_perfs().nb_iterations;
     const scalar_t norm_residuals = get_perfs().current_residual;
     const scalar_t norm_update = get_perfs().current_update;
     //std::cout << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0 << std::endl;
     DEBUG << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0;
     if (norm_residuals < get_options().residuals_tolerance)
     {
         termcode = ParabolicDriverReturnCode::ResidualMinimized;
     }
+    else if (get_perfs().absolute_residual < get_options().absolute_tolerance)
+    {
+        termcode = ParabolicDriverReturnCode::ResidualMinimized;
+    }
     else if (nb_iterations > 0 and norm_update > 0.0 and norm_update < 1.01*get_options().step_tolerance)
     {
         if (norm_residuals  > get_options().threshold_stationary_point)
         {
             ERROR << "Stationary point detected !";
             termcode = ParabolicDriverReturnCode::StationaryPoint;
         }
         WARNING << "Error is minimized - may indicate a stationnary point";
         termcode = ParabolicDriverReturnCode::ErrorMinimized;
 
     }
     else if (nb_iterations >  get_options().maximum_iterations)
     {
         ERROR << "Maximum number of iteration reached (" << get_options().maximum_iterations << ")";
         termcode = ParabolicDriverReturnCode::MaxIterations;
     }
     else if (get_perfs().maximum_step_taken)
     {
         get_perfs().nb_consecutive_max_step_taken += 1;
         get_perfs().nb_max_step_taken += 1;
         if (get_perfs().nb_consecutive_max_step_taken == get_options().max_iterations_at_max_length) {
             ERROR << "Divergence detected - Maximum step length taken two many times";
             termcode = ParabolicDriverReturnCode::MaxStepTakenTooManyTimes;
         }
     }
     else
     {
         get_perfs().nb_consecutive_max_step_taken = 0;
     }
     get_perfs().return_code = termcode;
     return termcode;
 }
 
 template <class Program>
 double ParabolicDriver<Program>::compute_residuals_linesearch(
         Vector &update,
         scalar_t lambda,
         Vector &displacement
         )
 {
     Eigen::VectorXd velocity(m_velocity);
     Eigen::VectorXd residual = Eigen::VectorXd::Zero(get_neq());
     program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, velocity);
     compute_residuals(displacement, velocity, residual);
     return 0.5*residual.squaredNorm();
 }
 
 template <class Program>
 double ParabolicDriver<Program>::compute_residuals_strang_linesearch(
         Vector &update,
         scalar_t lambda,
         Vector &displacement
         )
 {
     Eigen::VectorXd velocity(m_velocity);
     Eigen::VectorXd residual = Eigen::VectorXd::Zero(get_neq());
     program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, velocity);
     compute_residuals(displacement, velocity, residual);
     return update.dot(residual);
 }
 
 
 template <class Program>
 void ParabolicDriver<Program>::update_variable(
         const Vector& update,
         scalar_t lambda,
         Vector& displacement
         )
 {
     program().update_solution(update, lambda, get_options().alpha*m_current_dt,
                               m_predictor, displacement, m_velocity);
 }
 
 template <class Program>
 ParabolicDriverReturnCode ParabolicDriver<Program>::linesearch(
         Vector &update,
         Vector &displacements
         )
 {
     base::is_step_too_long(update);
     get_perfs().maximum_step_taken = false;
 
     scalar_t lambda;
     ParabolicLinesearchReturnCode retcode;
     switch (get_options().linesearch) {
     case ParabolicLinesearch::Bactracking:
         retcode = backtracking_linesearch(update, displacements, lambda);
         break;
     case ParabolicLinesearch::Strang:
         retcode = strang_linesearch(update, displacements, lambda);
         break;
     default:
         throw std::runtime_error("Linesearch type for Parabolic driver is not recognized");
         break;
     }
     if (retcode != ParabolicLinesearchReturnCode::Success)
     {
         return ParabolicDriverReturnCode::LinesearchFailed;
     }
 
     update_variable(update, lambda, displacements);
     update *= lambda;
 
     return ParabolicDriverReturnCode::NotConvergedYet;
 }
 
 namespace internal
 {
 //! \brief Return true if a and b have the same sign
 inline bool have_same_sign(double a, double b)
 {
     return (std::copysign(1., a) * std::copysign(1., b) > 0);
 }
 } // end namespace internal
 
 template <class Program>
 ParabolicLinesearchReturnCode ParabolicDriver<Program>::strang_linesearch(
         Vector& update,
         Vector& displacements,
         scalar_t& lambda
         )
 {
     DEBUG << "Strang linesearch";
     Eigen::VectorXd xp(displacements);
 
     const scalar_t s_tol = 0.5;
     const scalar_t s_max = 2.0;
     const int lin_max = 10;
 
     scalar_t s_b = 0.0;
     scalar_t g_b = 0.5*m_residuals.squaredNorm();
     scalar_t s_a = 1.0;
     scalar_t g_a = compute_residuals_strang_linesearch(update, 1.0, xp);
     scalar_t newtlen = (scaling().asDiagonal()*update).norm();
 
    // const scalar_t s_r = s_a;
     const scalar_t g_r = g_a;
 
     if (std::abs(g_a) <= s_tol*std::abs(g_b))
     {
         DEBUG << "Skip linesearch ";
         lambda = 1.0;
         if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
             get_perfs().maximum_step_taken = true;
         }
         return ParabolicLinesearchReturnCode::Success;
     }
 
     while (internal::have_same_sign(g_a, g_b) and s_a < s_max)
     {
         s_b = s_a;
         s_a = 2*s_a;
         g_b = g_a;
         g_a = compute_residuals_strang_linesearch(update, s_a, xp);
     }
 
     scalar_t g_ = g_a;
     scalar_t g_0 = g_a;
     scalar_t s = s_a;
 
     int l;
     for (l=0; l<lin_max; ++l)
     {
         if (internal::have_same_sign(g_a, g_b) or
                     std::abs(g_) < s_tol*std::abs(g_0)
                     or std::abs(s_a - s_b) < s_tol*(s_a + s_b)/2 )
             break;
 
         s = s_a - g_a * ( s_a - s_b)/(g_a - g_b);
         g_ = compute_residuals_strang_linesearch(update, s, xp);
         if (internal::have_same_sign(g_, g_a))
             g_b = g_b/2.;
         else
         {
             s_b = s_a;
             g_b = g_a;
         }
         s_a = s;
         g_a = g_;
 
     }
     if (l == lin_max) return ParabolicLinesearchReturnCode::MaximumIterations;
     if (g_a >= g_r) {
         WARNING << "Failed to find better update in Strang linesearch";
         //lambda = 0.1*s_r;
         return backtracking_linesearch(update, displacements, lambda);
     }
     lambda = s_a;
     if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
         get_perfs().maximum_step_taken = true;
     }
     return ParabolicLinesearchReturnCode::Success;
 }
 
 template <class Program>
 ParabolicLinesearchReturnCode ParabolicDriver<Program>::backtracking_linesearch(
         Vector& update,
         Vector& displacements,
         scalar_t& lambda_out
         )
 {
     // References
     // ----------
     //    - Algo A6.3.1 : Dennis and Schnabel (1983)
     //    - Nocedal & Wrigth (2006)
     DEBUG << "Linesearch";
     Eigen::VectorXd xp(displacements);
     double fcp;
 
     int retcode = 2; // 2 not converged, 1 problem, 0 success
     const scalar_t alpha = 1e-4;
     scalar_t newtlen = (scaling().asDiagonal()*update).norm();
     scalar_t init_slope = m_gradient.dot(update);
     const scalar_t rellength = update_norm(update);
     const scalar_t minlambda = get_options().step_tolerance / rellength;
     scalar_t lambda = 1.0;
     scalar_t lambda_prev = lambda;
 
     scalar_t merit_value = 0.5*m_residuals.squaredNorm();
     // new residual
     fcp = compute_residuals_linesearch(update, lambda, xp);
     // Skip linesearch if enough progress is done
     // ------------------------------------------
     if (fcp < get_options().coeff_accept_newton_step *merit_value)
     {
         DEBUG << "Skip linesearch ";
         lambda_out = 1.0;
         return ParabolicLinesearchReturnCode::Success;
     }
     // The linesearch
     // --------------
     scalar_t fc = merit_value;
     scalar_t fcp_prev;
     int cnt = 0;
     do
     {
         SPAM << "cnt : " <<cnt << " - lambda : " << lambda << " # " << fc << " # " << fcp << " # " << init_slope;
         if (fcp <= fc + alpha*lambda*init_slope) //pg760 Fachinei2003
         {
             retcode = 0;
             if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
                 get_perfs().maximum_step_taken = true;
             }
             break;
         }
         else if (lambda < minlambda)
         {
             retcode = 0;
             lambda = minlambda;
             //retcode = 1;
             break;
         }
         else
         {
             //WARNING << "denom : " << fcp - fc -init_slope;
             // Select a new step length
             // - - - - - - - - - - - -
             double lambdatmp;
             if (cnt == 0) { // only a quadratic at the first
                 lambdatmp = - init_slope / (2*(fcp - fc -init_slope));
             }
             else
             {
                 const scalar_t factor = 1.0 /(lambda - lambda_prev);
                 const scalar_t x1 = fcp - fc - lambda*init_slope;
                 const scalar_t x2 = fcp_prev - fc - lambda_prev*init_slope;
 
                 const scalar_t a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev));
                 const scalar_t b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev));
 
                 if (a == 0)
                 { // cubic interpolation is in fact a quadratic interpolation
                     lambdatmp = - init_slope/(2*b);
                 }
                 else
                 {
                     const scalar_t disc = b*b-3*a*init_slope;
                     lambdatmp = (-b+std::sqrt(disc))/(3*a);
                 }
                 if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda;
             }
             //WARNING << "lambdatmp : " << lambdatmp;
             lambda_prev = lambda;
             fcp_prev = fcp;
             if (lambdatmp < 0.1*lambda) {
                 lambda = 0.1 * lambda;
             } else {
                 lambda = lambdatmp;
             }
         }
         if (not std::isfinite(lambda))
         {
             ERROR << "Lambda is non finite - we stop";
             return ParabolicLinesearchReturnCode::Divergence;
         }
         fcp = compute_residuals_linesearch(update, lambda, xp);
         //xp.velocity = x.velocity;
         ++cnt;
     } while(retcode == 2 and cnt < 50);
     //WARNING << "Lambda : " << lambda << " - iterations : " << cnt;
     if (cnt == 50)
     {
         ERROR << "Too much linesearch iterations ! We stop";
         return ParabolicLinesearchReturnCode::MaximumIterations;
     }
     lambda_out = lambda;
     switch (retcode) {
     case 0:
         return ParabolicLinesearchReturnCode::Success;
         break;
     case 1:
         return ParabolicLinesearchReturnCode::LambdaTooSmall;
     default:
         return ParabolicLinesearchReturnCode::NotSupposedToHappen;
         break;
     }
 }
 
 template <class Program>
 void ParabolicDriver<Program>::set_velocity(Vector& velocity_vector)
 {
     m_velocity.resize(program().get_tot_ndf());
     m_velocity.setZero();
     for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
     {
         if (program().id_equation(dof) == no_equation)
         {
             m_velocity(dof) = velocity_vector(dof);
         }
     }
 }
 
 template <class Program>
 void ParabolicDriver<Program>::reset_velocity()
 {
     for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
     {
         if (program().id_equation(dof) == no_equation) continue;
         m_velocity(dof) = 0.0;
 
     }
 }
 
 } // end namespace dfpmsolver
 } // end namespace specmicp