diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index 440e196bc..3183b540b 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -353,6 +353,13 @@ end subroutine mechanical_results
 !--------------------------------------------------------------------------------------------------
 !> @brief calculation of stress (P) with time integration based on a residuum in Lp and
 !> intermediate acceleration of the Newton-Raphson correction
+!> @modified by Satya and achal
+!> check for detour i.e. if twinning is possible (we are not going ahead in Lp loop consistency) 
+!> checking by calling the deltaFp subroutine that should return 4 things
+!1) deltaFp= Cij (correspondance matrix) representing twinning shear and reorientation
+!2) -(twin volume fraction) for each twin system to make it harder for twinned material point to twin again by any twin system
+!3) jump in the last sampled volume fraction
+!4) logical true if twinning possible false if not occurring
 !--------------------------------------------------------------------------------------------------
 function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
 
@@ -450,7 +457,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
       Fe = matmul(matmul(A,B), invFi_new)
       call phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
                                         Fe, Fi_new, ph, en)
+      ! achal call Kinematic DeltaFp here
 
+
+
+
+
+                                        
       call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
                                          S, Fi_new, ph,en)
 
diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90
index bdee12f3e..4a9406b5d 100644
--- a/src/phase_mechanical_plastic_phenopowerlaw.f90
+++ b/src/phase_mechanical_plastic_phenopowerlaw.f90
@@ -365,7 +365,7 @@ associate(prm => param(ph), stt => state(ph), &
   call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
   dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
   call kinetics_tw(Mp,ph,en,dot_gamma_tw)
-  call plastic_kinematic_deltaFp(Mp,ph,en,twinJump,deltaFp)                                      ! delete this
+  if(en==1) call plastic_kinematic_deltaFp(Mp,ph,en,twinJump,deltaFp)                                      ! delete this
   !write(6,*)'deltaFp', deltaFp                                                                  ! delete this
   !write(6,*)'characteristicShearTwin', prm%gamma_char
   !write(6,*)'Schmid_twin',prm%P_sl
@@ -422,7 +422,12 @@ fdot_twin = (0.05_pReal*(abs(tau_tw)/state(ph)%xi_tw(:,en))**param(ph)%n_tw)/par
 
 !write(6,*) 'twin_var', twin_var                                                                      !delete this
 
-!write(6,*)'correspondanceMatrix', param(ph)%CorrespondanceMatrix(:,:,1)                              !delete this
+if (en==1) write(6,*)'correspondanceMatrix1', param(ph)%CorrespondanceMatrix(:,:,1)                              !delete this
+if (en==1) write(6,*)'correspondanceMatrix2', param(ph)%CorrespondanceMatrix(:,:,2)                              !delete this
+if (en==1) write(6,*)'correspondanceMatrix3', param(ph)%CorrespondanceMatrix(:,:,3)                              !delete this
+if (en==1) write(6,*)'correspondanceMatrix4', param(ph)%CorrespondanceMatrix(:,:,4)                              !delete this
+if (en==1) write(6,*)'correspondanceMatrix5', param(ph)%CorrespondanceMatrix(:,:,5)                              !delete this
+if (en==1) write(6,*)'correspondanceMatrix6', param(ph)%CorrespondanceMatrix(:,:,6)                              !delete this
 
 call RANDOM_NUMBER(random)