diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index 11ed0e7fde825f4627c22fdf52a1cbaca64245a5..495b8b7425f5fbd97991b41bf44f1f9274dc64fd 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -5847,8 +5847,7 @@ void nonLinearMechSolver::computeBodyForceVector(){
     this->getMicroBC()->setTime(1.); 
     G1 = this->getMicroBC()->getSecondOrderKinematicalVariable();
     this->getMicroBC()->setTime(curtime);
-    //G1 = this->getMicroBC()->getSecondOrderKinematicalVariable();
-    
+       
     G1 -=G0;
     for (int idom=0; idom<domainVector.size(); idom++){
         partDomain* dom = domainVector[idom];
@@ -9374,7 +9373,7 @@ double nonLinearMechSolver::computeRightHandSide()
     if (_systemType == nonLinearMechSolver::MULT_ELIM){
       _pAl->assembleConstraintResidualToSystem();
       this->computeStiffMatrix();
-      if(computedUdF() and _pathFollowing and usePresentModuli()){
+      if(computedUdF() and _pathFollowing){
          this->computeBodyForceVector();
       }
     }
@@ -13057,7 +13056,7 @@ void nonLinearMechSolver::extractAveragePropertiesByInSystemCondensation(homogen
     
    if(needdUdF)
   {
-   // _ipf->setdFmdFM(homoData->getdUnknowndF(), IPStateBase::current);
+    _ipf->setdFmdFM(homoData->getdUnknowndF(), IPStateBase::current);
     if(needdUdG)
       _ipf->setdFmdGM(homoData->getdUnknowndG(), IPStateBase::current);
 
@@ -15231,7 +15230,7 @@ double nonLinearMechSolver::solveMicroSolverSNL(){
       if(dUdF_Flag){
          homogenizedData* homoData = this->getHomogenizationState(IPStateBase::current);
          this->extractdUdFByInSystemCondensation(homoData);
-       //  _ipf->setdFmdFM(homoData->getdUnknowndF(), IPStateBase::current);
+         _ipf->setdFmdFM(homoData->getdUnknowndF(), IPStateBase::current);
      }
      else{
          dUdF_Flag = true;
@@ -15242,9 +15241,7 @@ double nonLinearMechSolver::solveMicroSolverSNL(){
     int niteNR = 0;
     if (_pathFollowing)
     {
-      if(computedUdF() and !usePresentModuli()){
-         this->computeBodyForceVector();
-      }
+
       pathSys->setPathFollowingIncrement(dt);
       _pAl->updateConstraint(_ipf);
       
diff --git a/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h b/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
index a3e266cf5ed63842c3eea8610357aa8aaf5ddbb4..8cc34fb380dd635f092efd524c27ce1fb174bc21 100644
--- a/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
+++ b/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
@@ -40,7 +40,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
     int _iter; //  specify predictor and corrector
     double uScale2; // equal to _uScale^2
     bool _isAllocatedBodyForceVector;
-    double _loadReduceRatio =1.0;
+    
 
 
   public:
@@ -201,7 +201,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
       _controlParameterPrev = 0.;
       _controlStep = 0.;
       _setScale = false;
-      _loadReduceRatio = 1.0;
+  
   
       _try(VecAssemblyBegin(_stateStep));
       _try(VecAssemblyEnd(_stateStep));
@@ -215,7 +215,6 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
     
     virtual void resetScaleParameter(const double val){
        _setScale = false;
-       _loadReduceRatio = _loadReduceRatio*val;
     }   
 
     virtual int systemSolve(){
@@ -234,13 +233,17 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
       _try(MatMultAdd(this->_CT,_q,_qeff,_qeff));
       _try(VecScale(_qeff,this->_scale));
       
-      if (_isAllocatedBodyForceVector)
-          _try(MatMultAdd(this->_Q,_BF,_qeff,_qeff));
-          
-      //  try to reduce load vector
-      PetscScalar reduce = 1./_loadReduceRatio;  
-      _try(VecScale(_qeff,reduce));    
-
+      if (_isAllocatedBodyForceVector){
+       //  if(_iter==0){
+       //    _try(MatMultAdd(this->_Q,_BF,_qeff,_qeff));
+       //  }
+       //  else{
+           PetscScalar inv_dt = 1./_pseudoTimeIncrement;
+           _try(VecScale(_BF, inv_dt));
+           _try(MatMultAdd(this->_Q,_BF,_qeff,_qeff));
+      //   }
+       }
+      
       _try(VecAssemblyBegin(_qeff));
       _try(VecAssemblyEnd(_qeff));
 
@@ -260,7 +263,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
         }
         PetscScalar A0 =  _aU*vTv/uScale2+_aL;
 
-        double a1 = _loadReduceRatio*_pseudoTimeIncrement/sqrt(A0);
+        double a1 = _pseudoTimeIncrement/sqrt(A0);
         double a2 = -1.*a1;
 
         PetscScalar uprevTu;
@@ -268,14 +271,14 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
         
         if (uprevTu>= 0){
            _controlStep = a1;
-            _controlParameter +=  a1/_loadReduceRatio;
+            _controlParameter +=  a1;
             _try(VecAXPY(_stateStep,a1,this->_x));
             _try(VecAXPY(this->_xsol,a1,this->_x));
 
         }
         else{
           _controlStep = a2;
-          _controlParameter +=  a2/_loadReduceRatio;
+          _controlParameter +=  a2;
           _try(VecAXPY(_stateStep,a2,this->_x));
           _try(VecAXPY(this->_xsol,a2,this->_x));
         }
@@ -342,7 +345,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
         PetscScalar D = _aU*_aU*(vTdr/uScale2)*(vTdr/uScale2) - A*drTdr/uScale2*_aU;
         PetscScalar E = (_aL*_controlStep+ _aU*vTdu/uScale2)*(vTdr/uScale2)*_aU - A*duTdr/uScale2*_aU;
         PetscScalar F = (_aL*_controlStep+ _aU*vTdu/uScale2)*(_aL*_controlStep+ _aU*vTdu/uScale2)
-                        -A*(_aU*duTdu/uScale2+_aL*_controlStep*_controlStep - _loadReduceRatio*_loadReduceRatio*_pseudoTimeIncrement*_pseudoTimeIncrement);
+                        -A*(_aU*duTdu/uScale2+_aL*_controlStep*_controlStep - _pseudoTimeIncrement*_pseudoTimeIncrement);
 
         double delta = E*E -D*F;
         PetscScalar _beta = 1.;
@@ -373,7 +376,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
         C*= _aU;
         C /= uScale2;
         double dlamda2 = _controlStep*_controlStep*_aL;
-        dlamda2 -= _loadReduceRatio*_loadReduceRatio*_pseudoTimeIncrement*_pseudoTimeIncrement;
+        dlamda2 -= _pseudoTimeIncrement*_pseudoTimeIncrement;
 
         C += dlamda2;
 
@@ -403,7 +406,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
             s = a2;
           }
           _controlStep += s;
-          _controlParameter += s/_loadReduceRatio;
+          _controlParameter += s;
           _try(VecAXPY(_stateStep,s,this->_x));
           _try(VecAXPY(this->_xsol,_beta,dr));
           _try(VecAXPY(this->_xsol,s,this->_x));
@@ -447,8 +450,6 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
     virtual void setPathFollowingIncrement(const double dt) {};
     virtual double getStateParameter() const {return 0.;};
 		virtual void resetControlParameter() {};
-		virtual void resetScaleParameter(){};
-		virtual void resetScaleParameter() {};
 		virtual void addToPathFollowingContraint(const scalar& val) {};
 		virtual void addToDPathFollowingConstraintDUnknown(const int col, const scalar& val) {};
 		virtual void addToDPathFollowingConstraintDControlParameter(const scalar& val) {};
diff --git a/dG3D/src/dG3DHomogenizedTangentTerms.cpp b/dG3D/src/dG3DHomogenizedTangentTerms.cpp
index f08cc0a27c3dccad5728e64fc9bbc9206bca11e9..79eb95eef755a0c69b8d436c8c4a606baf0cd5e0 100644
--- a/dG3D/src/dG3DHomogenizedTangentTerms.cpp
+++ b/dG3D/src/dG3DHomogenizedTangentTerms.cpp
@@ -447,20 +447,25 @@ void dG3DBodyForceTermPathFollowing::get(MElement *ele,int npts,IntPt *GP,fullVe
   {
     const IPStateBase *ips  = (*vips)[i];
     const dG3DIPVariableBase *ipv  = static_cast<const dG3DIPVariableBase*>(ips->getState(IPStateBase::current));
+    const dG3DIPVariableBase *ipvprev  = static_cast<const dG3DIPVariableBase*>(ips->getState(IPStateBase::previous));
 		
     const std::vector<TensorialTraits<double>::ValType>&  Vals = ipv->f(sp,ele,GP[i]);		
     const double weight = GP[i].weight;
     double &detJ = ipv->getJacobianDeterminant(ele,GP[i]);
 		
     const double ratio = detJ*weight;
-    SVector3 Bm;
-    ipv->computeBodyForce(G, Bm);
+  //  SVector3 Bm;
+  //  ipv->computeBodyForce(G, Bm);
+    const SVector3 &Bm0 = ipvprev->getConstRefToBm();
+    const SVector3 &Bm = ipv->getConstRefToBm();
+  
+    
       
     if (_dom->getHomogenizationOrder() == 2 and ipv->hasBodyForceForHO())
     {
        for (int j=0; j<nbFF; j++)
           for (int k=0; k<3; k++){
-              bm(j+k*nbFF)+= Bm(k)*Vals[j+k*nbFF]*ratio;
+              bm(j+k*nbFF)+= (Bm(k)-Bm0(k))*Vals[j+k*nbFF]*ratio;
           }           
     }	
   }
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 139927232418615b9ca3da63cc56ebf217854f28..016944b880bd11b77ff359cd264bc8b04bd42e05 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -1176,7 +1176,7 @@ void dG3DIPVariable::computeBodyForce(const STensor33& G, SVector3& B) const{
 
 
 void dG3DIPVariable::setValueForBodyForce(const STensor43& K, const bool usePresentModuli){
-   STensor43& K = elasticTangentModuli;
+  // STensor43& K = elasticTangentModuli;
 
  if(hasBodyForceForHO()){
    static STensor33 tmp1;
@@ -1194,8 +1194,7 @@ void dG3DIPVariable::setValueForBodyForce(const STensor43& K, const bool usePres
         }
       }
     }
-   // K.print("Elastic tensor");
-  //  Msg::Info("Body Force is  %f,%f\n", _bodyForceForHO->Bm(0),_bodyForceForHO->Bm(1));
+   
     STensorOperation::zero(_bodyForceForHO->dBmdGM);
     for (int i=0; i<3; i++){
       for (int j=0; j<3; j++){