diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
index 2dc82d062637580a37b3fb3a38223bef9759ebed..480e9295e495a01d55ea9e034bf5cfe9abf41003 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
@@ -274,7 +274,7 @@ void mlawNonLinearTVE::setViscoElasticData_Bulk(const int i, const double Ki, co
     
     _temFunc_Ki[i-1] = new constantScalarFunction(1.);
 
-    Msg::Info("setting: element=%d Ki = %e ki = %e, ",i-1,Ki,ki);
+    // Msg::Info("setting: element=%d Ki = %e ki = %e, ",i-1,Ki,ki);
   }
 };
 
@@ -289,7 +289,7 @@ void mlawNonLinearTVE::setViscoElasticData_Shear(const int i, const double Gi, c
     _temFunc_Gi[i-1] = new constantScalarFunction(1.);
     _temFunc_Alphai[i-1] = new constantScalarFunction(1.);
 
-    Msg::Info("setting: element=%d Gi = %e gi = %e, ",i-1,Gi,gi);
+    // Msg::Info("setting: element=%d Gi = %e gi = %e, ",i-1,Gi,gi);
   }
 };
 
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index f268df7f20a0c0233715b9c53bfc4a1f26c908f1..6e9a010cc5814678f313da08dd0709e4853d9efa 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -1113,7 +1113,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
   double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
 
   // double Gamma = 0.;   // flow rule parameter
-  double& Gamma = q1->_Gamma;   // flow rule parameter
+  double& Gamma = q1->_Gamma;   // flow rule parameter // still initialised to zero 
   double PhiEq = PhiEqpr;   // current effective deviator  
 
 
@@ -1122,13 +1122,21 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
   static fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
   this->getYieldCoefficients(q1,a);
 
-  double Hb(0.), dHbdT(0.), ddHbddT(0.);
+  double Hb(0.), dHb(0.), dHbdT(0.), ddHbddT(0.);
   if (q1->_ipKinematic != NULL){
     Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+    dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter
     dHbdT = Hb * q1->_ipKinematic->getDRDT();  // make sure to normalise DRDT while defining in python file //NEW
     ddHbddT = Hb * q1->_ipKinematic->getDDRDDT();
   }
 
+  static fullVector<double> m(2);
+  getChabocheCoeffs(m,0,q1);
+  double Cx(0.), Px(0.);
+  getPlasticPotential(ptilde, PhiEq, Px);   // CHANGET
+  double Gt = Ge + pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
+  double Kt = Ke + pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
+  
   //a.print("a init");
 
   static STensor3 devN; // dev part of yield normal
@@ -1148,14 +1156,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
 
   double dDgammaDGamma = 0.;
   double Dgamma = 0.; // eqplastic strain
-  
-  static fullVector<double> m(2);
-  getChabocheCoeffs(m,0,q1);
-  double Cx(0.), Px(0.);
-  getPlasticPotential(ptilde, PhiEq, Px); 
-  Cx = m(0)*(q0->_epspbarre + Dgamma) + m(1)*Px - 1/Hb * dHbdT + 1;  
-  double Gt = Ge + pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
-  double Kt = Ke + pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
+  Cx = m(0)*Dgamma + m(1)*Px - 1/Hb * dHbdT * (T-T0)  + 1;  // CHECK -> CHANGED   m(0)*(q0->_epsbarre + Dgamma) and no del
 
   if (q1->dissipationIsBlocked()){
     q1->getRefToDissipationActive() = false;
@@ -1186,8 +1187,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
         
         double dCxdDgamma = m(0);
         double Stemp = STensorOperation::doubledot(devPhi,devXn);
-        double He = (PhiEq*6*Gamma*kk*kk*Hb - Stemp*3.*1./PhiEq)/(2*Cx*Cx*u) * dCxdDgamma; 
-        double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma;
+        double He = (PhiEq*6*Gamma*kk*kk*Hb - Stemp*3.*1./PhiEq)/(2*Cx*Cx*u) * dCxdDgamma;  // CHECK -> CHANGE
+        double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma; // CHECK
         dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
 
         if (Gamma>0 and etaOverDt>0)
@@ -1403,7 +1404,6 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
       //2. get DCeprDCepr and DCeinvprDCepr
       static STensor43 DCeprDCepr, DCeinvprDCepr;
       DCeprDCepr = _I4;
-
       for (int i=0; i<3; i++)
         for (int s=0; s<3; s++)
           for (int k=0; k<3; k++)