diff --git a/src/core.cpp b/src/core.cpp
index b2a764058d9e6ed5237047b384260ae8ddae5df1..dc79a75e01e7fc701656053c24a44ead28c88183 100644
--- a/src/core.cpp
+++ b/src/core.cpp
@@ -367,11 +367,22 @@ public:
             Array<ADComplex, Shape<Dynamic,Dynamic,4,4,4> >
                 U({nx_, ny_, 4, 4, 4});
 
+            const size_t invdir[4] = {2, 3, 0, 1};
+            const int neigh[4][2] = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}};
+
             for (size_t i = 0; i < nx_; ++i) {
                 for (size_t j = 0; j < ny_; ++j) {
                     for (size_t k = 0; k < 4; ++k) {
-                        U.part(i,j,k).matrix() =
-                            U_.part(i,j,k).matrix() + J.part(i,j,k).matrix() * U_.part(i,j,k).matrix();
+                        U.part(i,j,k).matrix() = U_.part(i,j,k).matrix();
+
+                        const int i2 = i + neigh[k][0];
+                        const int j2 = j + neigh[k][1];
+                        if (i2 >= (int)nx_ || j2 >= (int)ny_)
+                            continue;
+
+                        U.part(i,j,k).matrix() = U_.part(i,j,k).matrix()
+                            + J.part(i,j,k).matrix() * U_.part(i,j,k).matrix()
+                            - U_.part(i,j,k).matrix() * J.part(i2,j2,invdir[k]).matrix();
                     }
                 }
             }
diff --git a/usadelndsoc/solver.py b/usadelndsoc/solver.py
index 7d27eddbcbf6ca32863a8f9b7077765eab39bb90..3d5cc7521576d71e6639910fd2db85c8854e0e1f 100644
--- a/usadelndsoc/solver.py
+++ b/usadelndsoc/solver.py
@@ -596,7 +596,7 @@ class Result:
 
     @property
     def J(self):
-        return self._core.grad_J(self.Phi).transpose(0, 1, 2, 4, 3) * (-1 / 8)
+        return self._core.grad_J(self.Phi).transpose(0, 1, 2, 4, 3) * (-1 / 16)
 
     @property
     def S(self):