diff --git a/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp b/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
index 56d9318..16c3e82 100644
--- a/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
+++ b/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
@@ -1,143 +1,142 @@
 #include "equilibrium_stagger.hpp"
 #include "variables.hpp"
 
 #include "specmicp/adimensional/adimensional_system_solver.hpp"
 #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
 #include "utils/log.hpp"
 #include "reactmicp/solver/staggers_base/stagger_structs.hpp"
 
 #ifdef SPECMICP_USE_OPENMP
 #include <omp.h>
 #endif // SPECMICP_USE_OPENMP
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 // EquilibriumStagger::
 
 //! \brief Initialize the stagger at the beginning of an iteration
 void EquilibriumStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var)
 {
     m_dt = dt;
 
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
 
     // Initialize velocity using values from previous timestep
     for (index_t node=0; node<true_var->nb_nodes(); ++node)
     {
         if (true_var->is_fixed_composition(node)) continue;
         scalar_t alpha = 1.0;
         for (index_t component=1; component<true_var->nb_component(); ++component)
         {
             alpha = std::max(alpha,
                              dt*true_var->solid_concentration(node, component, true_var->chemistry_rate())
                                /true_var->solid_concentration(node, component, true_var->displacement())
                              );
         }
         auto solid_velocity = true_var->velocity().segment(
                     true_var->offset_node(node)+true_var->offset_solid_concentration(),
                     true_var->nb_component());
         auto solid_chemistry_rate = true_var->chemistry_rate().segment(
                     true_var->offset_node(node)+true_var->offset_solid_concentration(),
                     true_var->nb_component());
         solid_velocity = 1/alpha * solid_chemistry_rate;
     }
 }
 
 //! \brief Solve the equation for the timestep
 solver::StaggerReturnCode EquilibriumStagger::restart_timestep(VariablesBasePtr var)
 {
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
     int failed_chemistry = 0;
 #ifdef SPECMICP_USE_OPENMP
-#pragma omp parallel
+#pragma omp parallel  default(none) shared(true_var, failed_chemistry)
     {
-#pragma omp for schedule(runtime)
+#pragma omp for schedule(dynamic, 5)
         for (index_t node=0; node<true_var->nb_nodes(); ++node)
         {
-            if (true_var->is_fixed_composition(node)) continue;
-            int retcode = solve_one_node(node, true_var);
+            // only solve if necessary
+            if (true_var->is_fixed_composition(node) or failed_chemistry > 0) continue;
+            const auto retcode = solve_one_node(node, true_var);
             if (retcode > 0)
             {
-#pragma omp atomic
                 ++failed_chemistry;
-
             }
         }
     }
 #else
     {
         for (index_t node=0; node<true_var->nb_nodes(); ++node)
         {
             if (true_var->is_fixed_composition(node)) continue;
-            int retcode = solve_one_node(node, true_var);
+            const auto retcode = solve_one_node(node, true_var);
             if (retcode > 0)
             {
                 ++failed_chemistry;
                 break;
             }
         }
     }
 #endif // SPECMICP_USE_OPENMP
     if (failed_chemistry > 0)
         return solver::StaggerReturnCode::UnknownError;
     return solver::StaggerReturnCode::ResidualMinimized;
 }
 
 //!
 int EquilibriumStagger::solve_one_node(index_t node, SaturatedVariablesPtr true_var)
 {
     AdimensionalSystemConstraints constraints(get_constraints(node));
 
     constraints.total_concentrations = true_var->total_concentrations(node);
 
     AdimensionalSystemSolver adim_solver(true_var->get_database(),
                                          constraints,
                                          true_var->equilibrium_solution(node),
                                          m_options);
     Vector variables(true_var->equilibrium_solution(node).main_variables);
     micpsolver::MiCPPerformance perf =  adim_solver.solve(variables);
 
     micpsolver::MiCPSolverReturnCode retcode = perf.return_code;
     if (retcode <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
     {
 
         ERROR << "Failed to solve chemistry problem at node " << node
               << ", return code = " << static_cast<int>(retcode)
               << ", residual = " << perf.current_residual;
         ERROR << "Total concentration : \n" << constraints.total_concentrations;
         return  1;
     }
 
     true_var->equilibrium_solution(node) = adim_solver.get_raw_solution(variables);
     AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
                                                   true_var->get_database(),
                                                   m_options.units_set);
 
     for (index_t component=0; component<true_var->nb_component(); ++component)
     {
         const scalar_t c_aq = extractor.density_water()*extractor.total_aqueous_concentration(component);
         true_var->aqueous_concentration(node, component, true_var->displacement()) = c_aq;
 
         const scalar_t c_aq_0 = true_var->aqueous_concentration(node, component, true_var->predictor());
         const scalar_t vel_aq = (c_aq - c_aq_0)/m_dt;
         true_var->aqueous_concentration(node, component, true_var->velocity()) = vel_aq;
 
         const scalar_t c_sol = extractor.total_immobile_concentration(component);
         true_var->solid_concentration(node, component, true_var->displacement()) = c_sol;
 
         const scalar_t c_sol_0 = true_var->solid_concentration(node, component, true_var->predictor());
         const scalar_t vel_sol = (c_sol - c_sol_0)/m_dt;
         true_var->solid_concentration(node, component, true_var->velocity()) = vel_sol;
         true_var->solid_concentration(node, component, true_var->chemistry_rate()) =  vel_sol;
     }
 
     return 0;
 }
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp