diff --git a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h
index 84ac25c1..12e65d24 100644
--- a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h
@@ -603,7 +603,37 @@ void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
   conformable(_grid, q_in_1.Grid());
   conformable(_grid, q_in_2.Grid());
   conformable(_grid, q_out.Grid());
-  assert(0);
+
+  PropagatorField tmp_shifted(UGrid); 
+  PropagatorField g5Lg5(UGrid); 
+  PropagatorField R(UGrid); 
+  PropagatorField gmuR(UGrid); 
+
+    Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT,
+  };
+  Gamma gmu=Gamma(Gmu[mu]);
+
+  tmp_shifted=Cshift(q_in_1,mu,1);
+  g5Lg5=g5*tmp_shifted*g5;
+  R=q_in_2;
+  gmuR=gmu*R;
+
+  qout=adj(g5Lg5)*R;
+  qout+=adj(g5Lg5)*gmuR;
+
+  g5Lg5=g5*q_in_1*g5;
+  tmp_shifted=Cshift(q_in_2,mu,1);
+  Impl::multLinkField(R,this->Umu,tmp_shifted,mu);
+  gmuR=gmu*R;
+
+  qout-=adj(g5Lg5)*R;
+  qout+=adj(g5Lg5)*gmuR;
+
+  qout/=2;
 }