diff --git a/src/specmicp/reduced_system.cpp b/src/specmicp/reduced_system.cpp
index a1bcf23..78e06e5 100644
--- a/src/specmicp/reduced_system.cpp
+++ b/src/specmicp/reduced_system.cpp
@@ -1,405 +1,426 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : specmicp.cpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 #include <iostream>
 #include "reduced_system.hpp"
 #include "utils/log.hpp"
 
 namespace specmicp {
 
 inline double pow10(double x)
 {
     return std::pow(10.0, x);
 }
 
 const double Adebye = 0.5092;
 const double Bdebye = 0.3283;
 
 void ReducedSystem::number_eq(bool water_equation)
 {
     int neq = 0;
     m_ideq.reserve(m_thermo->nb_component()+m_thermo->nb_mineral());
     if (water_equation)
     {
         m_ideq.push_back(neq);
         ++neq;
     } else {
         m_ideq.push_back(no_equation);
     }
     for (int i=1;i<m_thermo->nb_component();++i)
     {
         // Remve components that does not exist
         if (m_tot_conc(i) == 0 and i!= 1)
         {
             INFO << "Component " << m_thermo->label_aqueous(i)
                  << "has total concentration equal to zero, we desactivate it" << std::endl;
             m_nonactive_component.push_back(i);
             m_ideq.push_back(no_equation);
             continue;
         }
         else
         {
             m_ideq.push_back(neq);
             ++neq;
         }
     }
     m_nb_free_vars = neq;
     for (int m=0; m<m_thermo->nb_mineral();++m)
     {
         bool can_precipitate = true;
         // Remove minerals that cannot precipitate
         for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it)
         {
             if (m_thermo->nu_min(m, *it) != 0)
             {
                 can_precipitate = false;
                 break; // this is not a mineral that can precipitate
             }
         }
         if (can_precipitate)
         {
             m_ideq.push_back(neq);
             ++neq;
         }
         else
         {
             m_ideq.push_back(no_equation);
         }
     }
     m_nb_tot_vars = neq;
     m_nb_compl_vars = m_nb_tot_vars - m_nb_free_vars;
     // aqueous species
     m_active_aqueous.reserve(m_thermo->nb_aqueous());
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         bool can_exist = true;
         for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it)
         {
             if (m_thermo->nu_aq(j,*it) != 0)
             {
                 can_exist = false;
             }
         }
         m_active_aqueous.push_back(can_exist);
     }
 
 }
 
 double ReducedSystem::residual_water(const Eigen::VectorXd& x) const
 {
     double res = m_tot_conc(0) - mass_water(x)/molar_mass_water;
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         res -= mass_water(x)*m_thermo->nu_aq(j, 0)*m_secondary_conc(j);
     }
     for (int m=0; m<m_thermo->nb_mineral(); ++m)
     {
         if (ideq_min(m) == no_equation) continue;
         res -= m_thermo->nu_min(m, 0)*x(ideq_min(m));
     }
     return res;
 }
 
 
 double ReducedSystem::residual_component(const Eigen::VectorXd& x, int i) const
 {
     const int idp = ideq_paq(i);
     double res = m_tot_conc(i) - mass_water(x)*pow10(x(idp));
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         if (not m_active_aqueous[j]) continue;
         res -= mass_water(x)*m_thermo->nu_aq(j, i)*m_secondary_conc(j);
     }
     for (int m=0; m<m_thermo->nb_mineral(); ++m)
     {
         if (ideq_min(m) == no_equation) continue;
         res -= m_thermo->nu_min(m, i)*x(ideq_min(m));
     }
     return res;
 }
 
 
 double ReducedSystem::residual_mineral(const Eigen::VectorXd& x, int m) const
 {
     double res = m_thermo->logkr_min(m);
     for (int i=1; i<m_thermo->nb_component(); ++i)
     {
         if (m_thermo->nu_min(m, i) != 0)
             res -= m_thermo->nu_min(m, i)*(x(ideq_paq(i)) + m_loggamma(i));
     }
     return res;
 }
 
 
 void ReducedSystem::get_residuals(const Eigen::VectorXd& x,
                              Eigen::VectorXd& residual)
 {
     set_secondary_concentration(x);
     if (ideq_w() != no_equation) residual(ideq_w()) = residual_water(x);
     for (int i=1; i<m_thermo->nb_component(); ++i)
     {
         if (ideq_paq(i) != no_equation) residual(ideq_paq(i)) = residual_component(x, i);
     }
     for (int m=0; m<m_thermo->nb_mineral(); ++m)
     {
         if (ideq_min(m) != no_equation) residual(ideq_min(m)) = residual_mineral(x, m);
     }
 }
 
  //non-optimized Finite difference, for test only
 //void ReducedSystem::get_jacobian(Eigen::VectorXd& x,
 //                  Eigen::MatrixXd& jacobian)
 //{
 //    const int neq = total_variables();
 //    Eigen::VectorXd res(total_variables());
 //    Eigen::VectorXd perturbed_res(total_variables());
 
 //    get_residuals(x, res);
 //    for (int j=0; j<neq; ++j)
 //    {
 //        double h = 1e-8*std::abs(x(j));
 //        //h = std::copysign(h, x(j));
 //        if (h==0) h = 1e-8;
 //        double tmp = x(j);
 //        x(j) += h;
 //        h = x(j) - tmp;
 //        get_residuals(x, perturbed_res);
 //        for (int i=0; i<neq; ++i)
 //        {
 //            jacobian(i, j) = (perturbed_res(i) - res(i))/h;
 //        }
 //        x(j) = tmp;
 //    }
 // std::cout << jacobian << std::endl;
 //    return;
 //}
 
 
 void ReducedSystem::get_jacobian(Eigen::VectorXd& x,
                   Eigen::MatrixXd& jacobian)
 {
     const double log10 = std::log(10.0);
     jacobian.resize(total_variables(), total_variables());
     jacobian.setZero();
     // water
     if (ideq_w() != no_equation)
     {
         double tmp = -1/molar_mass_water;
         for (int j=0; j<m_thermo->nb_aqueous(); ++j)
         {
             if ( m_thermo->nu_aq(j, 0) == 0 ) continue;
             tmp -= m_thermo->nu_aq(j, 0)*m_secondary_conc(j);
         }
         jacobian(0, 0) = tmp;
         for (int k=1; k<m_thermo->nb_component(); ++k)
         {
             if (ideq_paq(k) == no_equation) continue;
             double tmp = 0;
             for (int j=0; j<m_thermo->nb_aqueous(); ++j)
             {
                 tmp -= m_thermo->nu_aq(j,0)*m_thermo->nu_aq(j,k)*m_secondary_conc(j);
             }
             jacobian(0, ideq_paq(k)) = x(ideq_w())*log10 * tmp;
         }
         for (int m=0; m<m_thermo->nb_mineral(); ++m)
         {
             if (ideq_min(m) == no_equation) continue;
             jacobian(0, ideq_min(m)) = -m_thermo->nu_min(m, 0);
         }
     }
     // aqueous component
     for (int i=1; i<m_thermo->nb_component(); ++i)
     {
         if (ideq_paq(i) == no_equation) continue;
         const int idp = ideq_paq(i);
         // aqueous components
         for (int k=1; k<m_thermo->nb_component(); ++k)
         {
             if (ideq_paq(k) == no_equation) continue;
             double tmp_iip = 0;
             if (k == i) tmp_iip -= pow10(x(ideq_paq(i)))*log10;
             for (int j=0; j<m_thermo->nb_aqueous(); ++j)
             {
                 tmp_iip -= log10*m_thermo->nu_aq(j, i)*m_thermo->nu_aq(j, k)*m_secondary_conc(j);
             }
             jacobian(idp, ideq_paq(k)) = mass_water(x)*tmp_iip;
         }
         // minerals
         for (int m=0; m<m_thermo->nb_mineral(); ++m)
         {
             if (ideq_min(m) == no_equation) continue;
             jacobian(idp, ideq_min(m)) = - m_thermo->nu_min(m, i);
         }
         if (ideq_w() != no_equation) // Water
         {
             double tmp_iw = -std::pow(10.0, x(idp));
             for (int j=0; j<m_thermo->nb_aqueous(); ++j)
             {
                 if ( m_thermo->nu_aq(j, i) == 0 ) continue;
                 tmp_iw -= m_thermo->nu_aq(j, i)*m_secondary_conc(j);
             }
             jacobian(idp, ideq_w()) = tmp_iw;
         }
     }
     // mineral equilibrium
     for (int m=0; m<m_thermo->nb_mineral(); ++m)
     {
         const int idm = ideq_min(m);
         if (idm == no_equation) continue;
         for (int i=1; i<m_thermo->nb_component(); ++i)
         {
             if (ideq_paq(i) == no_equation) continue;
             jacobian(idm, ideq_paq(i)) = -m_thermo->nu_min(m, i);
         }
     }
     //std::cout << jacobian << std::endl;
 }
 
 
 void ReducedSystem::set_secondary_concentration(const Eigen::VectorXd& x)
 {
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         if (m_active_aqueous[j] == false)
         {
             m_secondary_conc(j) = 0;
             continue;
         }
         double conc = -m_thermo->logkr_aq(j) - m_loggamma(j+m_thermo->nb_component());
         for (int k=1; k<m_thermo->nb_component(); ++k)
         {
             if (m_thermo->nu_aq(j, k) == 0) continue;
             conc +=  m_thermo->nu_aq(j, k)*(x(ideq_paq(k)) + m_loggamma(k));
         }
         m_secondary_conc(j) = pow10(conc);
     }
 }
 
+
+void ReducedSystem::unsafe_set_secondary_concentration(const Eigen::VectorXd& x)
+{
+    for (int j=0; j<m_thermo->nb_aqueous(); ++j)
+    {
+        if (m_active_aqueous[j] == false)
+        {
+            m_secondary_conc(j) = 0;
+            continue;
+        }
+        double conc = -m_thermo->logkr_aq(j) - m_loggamma(j+m_thermo->nb_component());
+        for (int k=1; k<m_thermo->nb_component(); ++k)
+        {
+            if (m_thermo->nu_aq(j, k) == 0) continue;
+            conc +=  m_thermo->nu_aq(j, k)*(x(k) + m_loggamma(k)); // we don't use ideq_paq here !
+        }
+        m_secondary_conc(j) = pow10(conc);
+    }
+}
+
 void ReducedSystem::set_ionic_strength(const Eigen::VectorXd& x)
 {
     double ionic = 0;
     for (int i=1; i<m_thermo->nb_component(); ++i)
     {
         if (ideq_paq(i) == no_equation or m_thermo->charge_paq(i) == 0) continue;
         ionic += pow10(x(ideq_paq(i)))*std::pow(m_thermo->charge_paq(i),2);
     }
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         if (not m_active_aqueous[j] or m_thermo->charge_saq(j) == 0) continue;
         ionic += m_secondary_conc(j)*std::pow(m_thermo->charge_saq(j),2);
     }
     m_ionic_strength = ionic;
 }
 
 void ReducedSystem::compute_log_gamma(const Eigen::VectorXd& x)
 {
     set_ionic_strength(x);
     const double sqrti = std::sqrt(m_ionic_strength);
     for (int i=0; i<m_thermo->nb_component(); ++i)
     {
         if (ideq_paq(i) == no_equation)
         {
             m_loggamma(i) = 0;
             continue;
         }
         double tmp = 0;
         if (m_thermo->charge(i) != 0)
         {
             tmp = - Adebye*std::pow(2, m_thermo->charge(i))*sqrti;
             tmp /= (1 + m_thermo->ao_debye(i)*Bdebye*sqrti);
         }
         tmp += m_thermo->bdot_debye(i)*m_ionic_strength;
         m_loggamma(i) = tmp;
     }
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         if (not m_active_aqueous[j])
         {
             m_loggamma(m_thermo->nb_component() + j) = 0;
             continue;
         }
         double tmp = 0;
         if (m_thermo->charge_saq(j) != 0)
         {
             tmp = - Adebye*std::pow(2, m_thermo->charge_saq(j))*sqrti;
             tmp /= (1 + m_thermo->ao_debye_saq(j)*Bdebye*sqrti);
         }
         tmp += m_thermo->bdot_debye_saq(j)*m_ionic_strength;
         m_loggamma(m_thermo->nb_component() + j) = tmp;
     }
 }
 
 
 void ReducedSystem::aqueous_tot_conc(const Eigen::VectorXd& x,
                       Eigen::VectorXd& aq_totconc)
 {
     assert(aq_totconc.rows() == m_thermo->nb_component());
-    set_secondary_concentration(x); // just to be sure
+    unsafe_set_secondary_concentration(x); // just to be sure
+    aq_totconc(0) = x(ideq_w());
     for (int i=1; i<m_thermo->nb_component();++i)
     {
         if (ideq_paq(i) != no_equation)
         {
-            double tmp = pow10(x(ideq_paq(i)));
+            double tmp = pow10(x(i));
             for (int j=0; j<m_thermo->nb_aqueous();++j)
             {
                 if (m_thermo->nu_aq(j, i) == 0) continue;
                 tmp += m_thermo->nu_aq(j, i)*m_secondary_conc(j);
             }
             aq_totconc(i) = tmp;
         }
         else
         {
             aq_totconc(i) = 0;
         }
     }
 }
 
 
 double ReducedSystem::max_lambda(const Eigen::VectorXd& x, const Eigen::VectorXd& update)
 {
     //return 1.0/std::max(1.0,(-update.block(0,0,5,1).cwiseQuotient(0.95*x.block(0,0,5,1))).maxCoeff());
     if (ideq_w() != no_equation)
     {
      //   std::cout << "max lambda : " << 1.0/std::max(1.0, -update(0)/(0.95*x(0))) << std::endl;
         return 1.0/std::max(1.0, -update(0)/(0.5*x(0)));
     }
     else
     {
         return 1.0;
     }
 }
 
 void ReducedSystem::reasonable_starting_guess(Eigen::VectorXd& x)
 {
     x.resize(total_variables());
     x(0) = m_tot_conc(0)*specmicp::molar_mass_water;
     for (int i=1; i<m_thermo->nb_component(); ++i)
     {
         x(i) = std::log10(0.1*m_tot_conc(i)/x(0));
     }
     x.block(m_thermo->nb_component(), 0, m_thermo->nb_mineral(), 1).setZero();
 }
 
 void ReducedSystem::reasonable_restarting_guess(Eigen::VectorXd& x)
 {
-    x.resize(total_variables());
+    //x.conservativeResize(total_variables());
     x(0) = m_tot_conc(0)*specmicp::molar_mass_water;
     //for (int i=1; i<m_thermo->nb_component(); ++i)
     //{
     //    x(i) = std::log10(0.1*m_tot_conc(i)/x(0));
     //}
     //x.block(m_thermo->nb_component(), 0, m_thermo->nb_mineral(), 1).setZero();
     for (int i=m_thermo->nb_component(); i<m_thermo->nb_component()+m_thermo->nb_mineral(); ++i)
     {
         x(i) = 0.0;
     }
 }
 
 
 } // end namespace specmicp
diff --git a/src/specmicp/reduced_system.hpp b/src/specmicp/reduced_system.hpp
index 008f8a6..6a2ba12 100644
--- a/src/specmicp/reduced_system.hpp
+++ b/src/specmicp/reduced_system.hpp
@@ -1,136 +1,139 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : specmicp.hpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_REDUCEDSYSTEM_HPP
 #define SPECMICP_REDUCEDSYSTEM_HPP
 
 //! \file specmicp.hpp The MiCP program to solve speciation
 
 #include <memory>
 #include "thermodata.hpp"
 #include "micpsolver/micpprog.hpp"
 
 //! \namespace specmicp main namespace for the SpecMiCP solver
 namespace specmicp {
 
 class ReducedSystem : public micpsolver::MiCPProg<ReducedSystem>
 {
 public:
     const int no_equation  = -1;
 
     ReducedSystem(std::shared_ptr<ThermoData> ptrthermo,
              Eigen::VectorXd& totconc,
              bool water_eq=true):
         m_thermo(ptrthermo),
         m_tot_conc(totconc),
         m_secondary_conc(ptrthermo->nb_aqueous()),
         m_loggamma(ptrthermo->nb_component()+ptrthermo->nb_aqueous())
     {
         number_eq(water_eq);
         m_secondary_conc.setZero();
         m_loggamma.setZero();
 
     }
 
     // Variables
 
     int total_variables() const {return m_nb_tot_vars;}
     int nb_free_variables() const {return m_nb_free_vars;}
     int nb_complementarity_variables() const {return m_nb_compl_vars;}
 
     // Residual
 
     double residual_water(const Eigen::VectorXd& x) const;
     double residual_component(const Eigen::VectorXd& x, int i) const;
     double residual_equilibrium_aqueous(const Eigen::VectorXd& x, int j) const;
     double residual_mineral(const Eigen::VectorXd& x, int m) const;
 
     void get_residuals(const Eigen::VectorXd& x,
                             Eigen::VectorXd& residual);
 
     void get_jacobian(Eigen::VectorXd &x,
                       Eigen::MatrixXd& jacobian);
 
     void init_secondary_conc(Eigen::VectorXd& x);
 
     int ideq(int i) const {return m_ideq[i];}
     int ideq_w() const {return m_ideq[0];}
     int ideq_paq(int i) const {
         assert(i < m_thermo->nb_component());
         return m_ideq[i];
     }
     int ideq_min(int m) const {
         assert(m < m_thermo->nb_mineral());
         return m_ideq[m+m_thermo->nb_component()];
     }
 
     double mass_water(const Eigen::VectorXd& x) const
     {
         if (ideq_w() != no_equation) return x(ideq_w());
         else return 1.0;
     }
 
     void set_ionic_strength(const Eigen::VectorXd& x);
     void compute_log_gamma(const Eigen::VectorXd& x);
 
 
     void hook_start_iteration(const Eigen::VectorXd& x, double norm_residual) {
         //compute_log_gamma(x);
         not_in_linesearch = true;
         if (norm_residual < nb_free_variables())
         {
             // Use fix point iterations for non-ideality
             for (int i=0; i<2; ++i) // 2 seems to reduce the number of iterations, more does not help
             {
                 set_secondary_concentration(x);
                 compute_log_gamma(x);
             }
         }
     }
 
     void aqueous_tot_conc(const Eigen::VectorXd& x,
                           Eigen::VectorXd& aq_totconc);
 
     double max_lambda(const Eigen::VectorXd& x, const Eigen::VectorXd& update);
 
     const std::vector<int>& get_non_active_component() {return m_nonactive_component;}
 
     void reasonable_starting_guess(Eigen::VectorXd& x);
     void reasonable_restarting_guess(Eigen::VectorXd& x);
 
 private:
 
+    // For the normal algorithm
     void set_secondary_concentration(const Eigen::VectorXd& x);
+    // For the tot aq_conc pb
+    void unsafe_set_secondary_concentration(const Eigen::VectorXd& x);
     void number_eq(bool water_equation);
 
     std::shared_ptr<ThermoData> m_thermo;
     Eigen::VectorXd m_tot_conc;
 
 
     int m_nb_tot_vars;
     int m_nb_free_vars;
     int m_nb_compl_vars;
 
     std::vector<int> m_ideq;
 
     Eigen::VectorXd m_secondary_conc;
     double m_ionic_strength;
     Eigen::VectorXd m_loggamma;
 
     bool not_in_linesearch;
 
     std::vector<int> m_nonactive_component;
     std::vector<bool> m_active_aqueous;
 };
 
 } // end namespace specmicp
 
 
 #endif // SPECMICP_REDUCEDSYSTEM_HPP