diff --git a/src/specmicp/specmicp.cpp b/src/specmicp/specmicp.cpp
index 756cf30..b327e2d 100644
--- a/src/specmicp/specmicp.cpp
+++ b/src/specmicp/specmicp.cpp
@@ -1,218 +1,257 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : specmicp.cpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 #include <iostream>
 #include "specmicp.hpp"
 namespace specmicp {
 
+const double Adebye = 0.5092;
+const double Bdebye = 0.3283;
+
 void SpecProg::number_eq(bool water_equation)
 {
     int neq = 0;
     m_ideq.reserve(m_thermo->nb_component()+m_thermo->nb_aqueous()+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)
     {
         m_ideq.push_back(neq);
         ++neq;
     }
     for (int j=0; j<m_thermo->nb_aqueous();++j)
     {
         m_ideq.push_back(neq);
         ++neq;
     }
     m_nb_free_vars = neq;
     for (int m=0; m<m_thermo->nb_mineral();++m)
     {
         m_ideq.push_back(neq);
         ++neq;
     }
     m_nb_tot_vars = neq;
     m_nb_compl_vars = m_nb_tot_vars - m_nb_free_vars;
 
 }
 
 double SpecProg::residual_water(const Eigen::VectorXd& x) const
 {
     double res = m_tot_conc(ideq_w()) - 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)*std::pow(10, x( ideq_saq(j)));
     }
     for (int m=0; m<m_thermo->nb_mineral(); ++m)
     {
         res -= m_thermo->nu_min(m, 0)*x(ideq_min(m));
     }
-    return res;
+    return res/m_tot_conc(ideq_w());
 }
 
 
 double SpecProg::residual_component(const Eigen::VectorXd& x, int i) const
 {
     const int idp = ideq_paq(i);
     double res = m_tot_conc(idp) - mass_water(x)*std::pow(10, x(idp));
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         res -= mass_water(x)*m_thermo->nu_aq(j, i)*std::pow(10, x(ideq_saq(j)));
     }
     for (int m=0; m<m_thermo->nb_mineral(); ++m)
     {
         res -= m_thermo->nu_min(m, i)*x(ideq_min(m));
     }
-    return res;
+    return res/m_tot_conc(ideq_paq(i));
 }
 
 
 double SpecProg::residual_equilibrium_aqueous(const Eigen::VectorXd& x, int j) const
 {
     const int ids = ideq_saq(j);
-    double res = m_thermo->logkr_aq(j) + x(ids);
+    double res = m_thermo->logkr_aq(j) + x(ids) + m_loggamma(j+m_thermo->nb_component());
     for (int i=1; i<m_thermo->nb_component(); ++i)
     {
-        res -= m_thermo->nu_aq(j, i)*x(ideq_paq(i));
+        res -= m_thermo->nu_aq(j, i)*(x(ideq_paq(i)) + m_loggamma(i));
     }
     return res;
 }
 
 
 double SpecProg::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)
     {
-        res -= m_thermo->nu_min(m, i)*x(ideq_paq(i));
+        res -= m_thermo->nu_min(m, i)*(x(ideq_paq(i)) + m_loggamma(i));
     }
     return res;
 }
 
 
 void SpecProg::get_residuals(const Eigen::VectorXd& x,
-                             Eigen::VectorXd& residual) const
+                             Eigen::VectorXd& residual)
 {
     if (ideq_w() != no_equation) residual(ideq_w()) = residual_water(x);
     for (int i=1; i<m_thermo->nb_component(); ++i)
         residual(ideq_paq(i)) = residual_component(x, i);
+    if (not_in_linesearch) // if not in linesearch set log_gamma
+    {
+        double normr = residual.block(0, 0, m_thermo->nb_component(), 1).norm();
+        std::cout << "Normr : " << normr << std::endl;
+        if (normr < 1) set_ionic_strength(x);
+        not_in_linesearch = false;
+    }
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
         residual(ideq_saq(j)) = residual_equilibrium_aqueous(x, j);
     for (int m=0; m<m_thermo->nb_mineral(); ++m)
         residual(ideq_min(m)) = residual_mineral(x, m);
 }
 
 // non-optimized Finite difference, for test only
 //void SpecProg::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;
 //    }
 //    return;
 //}
 
 
 void SpecProg::get_jacobian(Eigen::VectorXd& x,
                   Eigen::MatrixXd& jacobian)
 {
     const double log10 = std::log(10.0);
     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;
             const int ids = ideq_saq(j);
             tmp -= m_thermo->nu_aq(j, 0)*std::pow(10, x(ids));
-            jacobian(0, ids) =  -x(ideq_w())*m_thermo->nu_aq(j, 0)*std::pow(10, x(ids))*log10;
+            jacobian(0, ids) =  -x(ideq_w())*m_thermo->nu_aq(j, 0)*std::pow(10, x(ids))*log10/m_tot_conc(ideq_w());
         }
-        jacobian(0, 0) = tmp;
+        jacobian(0, 0) = tmp/m_tot_conc(ideq_w());
         for (int m=0; m<m_thermo->nb_mineral(); ++m)
         {
             const int idm = m + nb_free_variables();
-            jacobian(0, idm) = -m_thermo->nu_min(m, 0);
+            jacobian(0, idm) = -m_thermo->nu_min(m, 0)/m_tot_conc(ideq_w());
         }
     }
     // aqueous component
     for (int i=1; i<m_thermo->nb_component(); ++i)
     {
         double tmp = -std::pow(10.0, x(i));
-        jacobian(i, i) = -x(0)*std::pow(10, x(ideq_paq(i)))*log10;
+        jacobian(i, i) = -x(0)*std::pow(10, x(ideq_paq(i)))*log10/m_tot_conc(ideq_paq(i));
         for (int j=0; j<m_thermo->nb_aqueous(); ++j)
         {
             if ( m_thermo->nu_aq(j, i) == 0 ) continue;
             const int ids = ideq_saq(j);
             tmp -= m_thermo->nu_aq(j, i)*std::pow(10, x(ids));
-            jacobian(i, ids) = -x(ideq_w())*m_thermo->nu_aq(j, i)*std::pow(10, x(ids))*log10;
+            jacobian(i, ids) = -x(ideq_w())*m_thermo->nu_aq(j, i)*std::pow(10, x(ids))*log10/m_tot_conc(ideq_paq(i));
         }
         for (int m=0; m<m_thermo->nb_mineral(); ++m)
         {
-            jacobian(i, ideq_min(m)) = - m_thermo->nu_min(m, i);
+            jacobian(i, ideq_min(m)) = - m_thermo->nu_min(m, i)/m_tot_conc(ideq_paq(i));
         }
-        jacobian(i, 0) = tmp;
+        jacobian(i, 0) = tmp/m_tot_conc(ideq_paq(i));
     }
     // Aqueous equilibrium
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         const int ids = ideq_saq(j);
         jacobian(ids, ids) = 1.0;
         for (int i=1; i<m_thermo->nb_component(); ++i)
         {
             jacobian(ids, i) = -m_thermo->nu_aq(j, i);
         }
     }
     // mineral equilibrium
     for (int m=0; m<m_thermo->nb_mineral(); ++m)
     {
         const int idm = ideq_min(m);
         for (int i=1; i<m_thermo->nb_component(); ++i)
         {
             jacobian(idm, i) = -m_thermo->nu_min(m, i);
         }
     }
 }
 
 
 void SpecProg::init_secondary_conc(Eigen::VectorXd &x)
 {
     for (int j=0; j<m_thermo->nb_aqueous(); ++j)
     {
         const int ids = ideq_saq(j);
         double sum = -m_thermo->logkr_aq(j);
         std::cout << j << " : " << m_thermo->logkr_aq(j) << std::endl;
         for (int i=1; i<m_thermo->nb_component(); ++i)
         {
             sum += m_thermo->nu_aq(j, i)*x(ideq_paq(i));
         }
         x(ids) = sum;
     }
 }
 
+void SpecProg::set_ionic_strength(const Eigen::VectorXd& x)
+{
+    double ionic = 0;
+    for (int i=1; i<m_thermo->nb_component(); ++i)
+    {
+        ionic += std::pow(10.0,x(ideq_paq(i)))*std::pow(2,m_thermo->charge_paq(i));
+    }
+    for (int j=0; j<m_thermo->nb_aqueous(); ++j)
+    {
+        ionic += std::pow(10.0, x(ideq_saq(j)))*std::pow(2,m_thermo->charge_saq(j));
+    }
+    if (ionic > 5) ionic = 5;
+    m_ionic_strength = 0; //0.5*ionic;
+    // std::cout << m_ionic_strength; << std::endl;
+}
+
+void SpecProg::compute_log_gamma(const Eigen::VectorXd& x)
+{
+    set_ionic_strength(x);
+    const double sqrti = std::sqrt(m_ionic_strength);
+    for (int j=0; j<m_thermo->nb_component()+m_thermo->nb_aqueous(); ++j)
+    {
+        double tmp = - Adebye*std::pow(2, m_thermo->charge(j))*sqrti;
+        tmp /= (1 + m_thermo->ao_debye(j)*Bdebye*sqrti);
+        tmp += m_thermo->bdot_debye(j)*m_ionic_strength;
+        m_loggamma(j) = tmp;
+    }
+}
+
 } // end namespace specmicp
diff --git a/src/specmicp/specmicp.hpp b/src/specmicp/specmicp.hpp
index ae0da9f..12c6c75 100644
--- a/src/specmicp/specmicp.hpp
+++ b/src/specmicp/specmicp.hpp
@@ -1,87 +1,103 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : specmicp.hpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_SPECMICP_HPP
 #define SPECMICP_SPECMICP_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 SpecProg : public micpsolver::MiCPProg<SpecProg>
 {
 public:
     const int no_equation  = -1;
 
     SpecProg(std::shared_ptr<ThermoData> ptrthermo,
              Eigen::VectorXd totconc,
              bool water_eq=true):
         m_thermo(ptrthermo),
-        m_tot_conc(totconc)
+        m_tot_conc(totconc),
+        m_loggamma(ptrthermo->nb_component()+ptrthermo->nb_aqueous())
     {
         number_eq(water_eq);
+        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) const;
+                            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 {return m_ideq[i];}
     int ideq_saq(int j) const {return m_ideq[j+m_thermo->nb_component()];}
     int ideq_min(int m) const {return m_ideq[m+m_thermo->nb_component()+m_thermo->nb_aqueous()];}
 
     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 types::vector_t& x) {
+        //compute_log_gamma(x);
+        not_in_linesearch = true;
+    }
+
 private:
     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;
+
+    double m_ionic_strength;
+    Eigen::VectorXd m_loggamma;
+
+    bool not_in_linesearch;
 };
 
 } // end namespace specmicp
 
 
 #endif // SPECMICP_SPECMICP_HPP
diff --git a/src/specmicp/thermodata.hpp b/src/specmicp/thermodata.hpp
index 1a077f1..bf79a79 100644
--- a/src/specmicp/thermodata.hpp
+++ b/src/specmicp/thermodata.hpp
@@ -1,76 +1,88 @@
 /*-------------------------------------------------------
 
  - Module : specmicp
  - File : thermodata.hpp
  - Author : Fabien Georget
 
  Copyright (c) 2014, Fabien Georget, Princeton University
 
 ---------------------------------------------------------*/
 
 #ifndef SPECMICP_SPECMICP_THERMODATA_HPP
 #define SPECMICP_SPECMICP_THERMODATA_HPP
 
 //! \file thermodata.hpp Thermodynamic data
 
 #include <Eigen/Dense>
 
 namespace specmicp {
 
 const double molar_mass_water = 18e-3;
 
 //! Contains the thermodynamic information
 //!
 //! Huge storage class
 class ThermoData
 {
 public:
     using id_t = int;
 
     ThermoData() {}
     ThermoData(int number_aqueous, int number_min,
-               Eigen::MatrixXd nu, Eigen::VectorXd logkr):
+               Eigen::MatrixXd& nu, Eigen::VectorXd& logkr,
+               Eigen::MatrixXd& param_aq):
         m_nb_component(nu.cols()),
         m_nb_aqueous(number_aqueous),
         m_nb_mineral(number_min),
         m_nu(nu),
-        m_logkr(logkr)
+        m_logkr(logkr),
+        m_aqparam(param_aq)
     {}
 
     int nb_component() {return m_nb_component;}
     int nb_aqueous() {return m_nb_aqueous;}
     int nb_mineral() {return m_nb_mineral;}
 
     double logkr(id_t id) {return m_logkr(id);}
     double logkr_aq(id_t id) {
         assert(id < nb_aqueous());
         return m_logkr(id);
     }
     double logkr_min(id_t id) {
         assert(id < nb_mineral());
         return m_logkr(id + nb_aqueous());
     }
     double nu(id_t ids, id_t idp) {return m_nu(ids, idp);}
     double nu_aq(id_t ids, id_t idp) {
         assert(ids < nb_aqueous());
         return m_nu(ids, idp);
     }
     double nu_min(id_t idm, id_t idp) {
         assert(idm < nb_mineral());
         return m_nu(idm + nb_aqueous(), idp);
     }
 
+    double charge(id_t id) {return m_aqparam(id, 0);}
+    double charge_paq(id_t idp) {return m_aqparam(idp, 0);}
+    double charge_saq(id_t ids) {return m_aqparam(ids+m_nb_component, 0);}
+    double ao_debye(id_t id) {return m_aqparam(id, 1);}
+    double ao_debye_paq(id_t idp) {return m_aqparam(idp, 1);}
+    double ao_debye_saq(id_t ids) {return m_aqparam(ids+m_nb_component, 1);}
+    double bdot_debye(id_t id) {return m_aqparam(id, 2);}
+    double bdot_debye_paq(id_t idp) {return m_aqparam(idp, 2);}
+    double bdot_debye_saq(id_t ids) {return m_aqparam(ids+m_nb_component, 2);}
 
 private:
     id_t m_nb_component;
     id_t m_nb_aqueous;
     id_t m_nb_mineral;
     std::vector<id_t> m_basis;
     Eigen::MatrixXd m_nu;
     Eigen::VectorXd m_logkr;
+    Eigen::MatrixXd m_aqparam;
 };
 
 } // end namespace specmicp
 
 
 #endif // SPECMICP_SPECMICP_THERMODATA_HPP
diff --git a/tests/specmicp/thermocarbo.cpp b/tests/specmicp/thermocarbo.cpp
index 9f62e87..eccb967 100644
--- a/tests/specmicp/thermocarbo.cpp
+++ b/tests/specmicp/thermocarbo.cpp
@@ -1,138 +1,166 @@
 #include <iostream>
 #include "specmicp/specmicp.hpp"
 #include "micpsolver/micpsolver.hpp"
 
 std::shared_ptr<specmicp::ThermoData>  set_thermodata()
 {
     int nmin = 5;
     Eigen::MatrixXd nu(10+nmin, 5);
     nu << 1, 0, -1, 0, 0,
           1, 0, -1, 1, 0,
          -1, 0,  1, 1, 0,
           0, 1, 1, 0, 0,
           -1, 0, 1, 0, 1,
           0, 0, -1, 0, 1,
           -1, 1, 1, 0, 1,
           0, 1, 0, 0, 1,
           0, 1, 0, 1, 0,
          -1, 1, 1, 1, 0,
        // minerals
           0, 1, 2, 0, 0,
           -1, 0, -1, 1, 0,
           -0.5567, 1.6667, 2.3333, 1, 0,
           -0.5, 0.8333, 0.6667, 1, 0,
            -1, 1, 1, 0, 1
             ;
 
     Eigen::VectorXd logkr(10+nmin);
     logkr << 13.9995, 4.8595, 0.1005, -1.2195,
             -3.6706, 7.6476, -6.8947, -1.1057,
             -1.2, -4.4995,
             // minerals
             -5.1995,
             1.476,
             -13.17,
             -8.0,
             -12.1505;
 
-    std::shared_ptr<specmicp::ThermoData> ptr = std::make_shared<specmicp::ThermoData>(10, nmin, nu, logkr);
+    Eigen::MatrixXd param_aqueous(Eigen::MatrixXd::Zero(5+10,3));
+    param_aqueous << 0,  0.0 , 0.0 ,  // H2O
+            2,  4.86, 0.15,  // Ca+2
+            -1, 10.65, 0.0 ,  // HO-
+            -1,  4.0 , 0.0 ,  // SiO(HO)3-
+            -1,  5.40, 0.0 ,  // HCO3-
+            1,  9.0 , 0.0 ,  // H+
+            0,  0.0 , 0.0 ,  // Si(OH)4
+            -2,  4.0 , 0.0 ,  // SiO2(OH)-2
+            1,  4.0 , 0.0 ,  // CaOH+
+            -2,  5.40, 0.0 ,  // CO3-2
+            0,  0.0 , 0.0 ,  // CO2(aq)
+            0,  0.0 , 0.0 ,  // CaCO3(aq)
+            1,  0.0 , 0.0 ,  // CaHCO3+
+            1,  4.0 , 0.0 ,  // CaSiO(HO)3+
+            0,  0.0 , 0.0    // CaSiO2(HO)2(aq)
+            ;
+
+
+    std::shared_ptr<specmicp::ThermoData> ptr = std::make_shared<specmicp::ThermoData>(10, nmin, nu, logkr, param_aqueous);
     return ptr;
 }
 
 std::shared_ptr<specmicp::ThermoData>  set_thermodata_aluminum()
 {
     Eigen::MatrixXd nu(5, 3);
     nu << 1, 1, -1,
           2, 1, -2,
           3, 1, -3,
           4, 1, -4,
           1, 0, -1
             ;
 
     Eigen::VectorXd logkr(5);
     logkr << 5.0, 10.1, 16.9, 22.7, 14.0;
 
-    std::shared_ptr<specmicp::ThermoData> ptr = std::make_shared<specmicp::ThermoData>(5, 0, nu, logkr);
+    Eigen::MatrixXd paramaq(Eigen::MatrixXd::Zero(5+3,3)); // ideal case
+
+    std::shared_ptr<specmicp::ThermoData> ptr = std::make_shared<specmicp::ThermoData>(5, 0, nu, logkr, paramaq);
     return ptr;
 }
 
 void solve_aluminium()
 {
     std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata_aluminum();
     Eigen::VectorXd totconc(3);
     totconc << 1.0/specmicp::molar_mass_water, 1e-6, -3e-6;
     //totconc << 1e-6, -3e-6;
 
     std::shared_ptr<specmicp::SpecProg> prog = std::make_shared<specmicp::SpecProg>(thermoptr, totconc, true);
     specmicp::micpsolver::MiCPSolver<specmicp::SpecProg> solver(prog);
     //solver.get_options().fvectol = 1e-6;
-    solver.get_options().maxstep  =  10.0;
+    solver.get_options().maxstep  =  20.0;
     solver.get_options().maxiter_maxstep = 100;
     //solver.get_options().factor_descent_condition = -1e-6;
     Eigen::VectorXd x(prog->total_variables());
     //x.setConstant(-1);
     x(0) = 1.0;
     x(1) = -9;
-    x(2) = -2;
+    x(2) = -13;
     prog->init_secondary_conc(x);
 
     std::cout << x << std::endl;
 
 //    Eigen::VectorXd residual(20);
 //    prog->get_residuals(x, residual);
 //    std::cout << "Residual \n ------ \n" << residual << std::endl;
 //    Eigen::MatrixXd jacob(20, 20);
 //    prog->get_jacobian(x, jacob);
 //    std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
 
     std::cout << "Return Code : " <<  (int) solver.solve(x) << std::endl;
 
 
     std::cout << "Solution \n -------- \n" << x << std::endl;
 }
 
 void solve()
 {
     std::shared_ptr<specmicp::ThermoData> thermoptr = set_thermodata();
     Eigen::VectorXd totconc(5);
-    totconc << 1.0/specmicp::molar_mass_water, 1.3, 1.9, 0.5, 0.2;
-
+    double m_c3s = 0.70;
+    double m_c2s = 0.3;
+    double m_h2co3 = 0.2;
+    totconc << 1/specmicp::molar_mass_water - 4*m_c3s - 3*m_c2s + m_h2co3,
+                        3*m_c3s +2*m_c2s,
+                        5*m_c3s + 3*m_c2s-m_h2co3,
+                        m_c3s + m_c2s,
+                        m_h2co3;
     std::shared_ptr<specmicp::SpecProg> prog = std::make_shared<specmicp::SpecProg>(thermoptr, totconc);
     specmicp::micpsolver::MiCPSolver<specmicp::SpecProg> solver(prog);
     //solver.get_options().fvectol = 1e-6;
     solver.get_options().maxstep  =  200.0;
     solver.get_options().maxiter_maxstep = 100;
+    solver.get_options().factor_descent_condition = 1e-4;
     Eigen::VectorXd x(prog->total_variables());
     //x.setConstant(-1);
-    x(0) = 1.0;
-    x(1) = -3;
-    x(2) = -3;
-    x(3) = -6;
-    x(4) = -6;
-    x(15) = 1.0;
-    x(16) = 0.0;
-    x(17) = 1.0;
-    x(18) = 0.0;
-    x(19) = 1.0;
+    x(0)  = 1.0;
+    x(1)  = 0.0;
+    x(2)  = 0.0;
+    x(3)  = 0.0;
+    x(4)  = 0.0;
+    x(15) =  0.0;
+    x(16) =  0.0;
+    x(17) =  0.0;
+    x(18) =  0.0;
+    x(19) =  0.0;
     prog->init_secondary_conc(x);
 
     std::cout << x << std::endl;
 
 //    Eigen::VectorXd residual(20);
 //    prog->get_residuals(x, residual);
 //    std::cout << "Residual \n ------ \n" << residual << std::endl;
 //    Eigen::MatrixXd jacob(20, 20);
 //    prog->get_jacobian(x, jacob);
 //    std::cout << "Jacobian \n ------ \n" << jacob << std::endl;
 
     std::cout << "Return Code : " <<  (int) solver.solve(x) << std::endl;
 
 
     std::cout << "Solution \n -------- \n" << x << std::endl;
 }
 
 int main()
 {
     solve();
     //solve_aluminium();
 }