diff --git a/fidimag/common/chain_method_integrators.py b/fidimag/common/chain_method_integrators.py
index 738250bb..823b8400 100644
--- a/fidimag/common/chain_method_integrators.py
+++ b/fidimag/common/chain_method_integrators.py
@@ -249,7 +249,7 @@ def run_for(self, n_steps):
             # Creep stage: minimise with a fixed eta
             while creepCount < self.maxCreep:
                 # Update spin. Avoid pinned or zero-Ms sites
-                self.band[:] = self.band_old + eta * self.etaScale * self.forces
+                self.band[:] = self.band_old + eta * self.etaScale * self.forces_old
                 normalise_spins(self.band)
 
                 self.trailAction
@@ -285,6 +285,7 @@ def run_for(self, n_steps):
                       f'eta = {eta:>5.4e}  '
                       f'action = {self.action:>5.4e}  action_old = {self.action_old:>5.4e}'
                       )
+                # print(self.forces)
 
                 # 10 seems like a magic number; we set here a minimum number of evaulations
                 if (nStart > self.nTrail * 10) and (deltaAction < self.actionTol):
diff --git a/fidimag/common/nebm_FS.py b/fidimag/common/nebm_FS.py
index 3522dd40..958c4c38 100644
--- a/fidimag/common/nebm_FS.py
+++ b/fidimag/common/nebm_FS.py
@@ -253,15 +253,6 @@ def compute_effective_field_and_energy(self, y):
 
             self.energies[i] = self.sim.compute_energy()
 
-        # TODO: move this calc to the action function
-        # Compute the gradient norm per every image
-        Gnorms2 = np.sum(self.gradientE**2, axis=1) / self.n_images
-        # Compute the root mean square per image
-        self.gradientENorm[:] = np.sqrt(Gnorms2)
-
-        # DEBUG:
-        # print('gradEnorm', self.gradientENorm)
-
         y.shape = (-1)
         self.gradientE.shape = (-1)
 
@@ -328,10 +319,20 @@ def compute_action(self):
         #                                       self._material_int,
         #                                       self.n_dofs_image_material
         #                                       )
+        
+        # NOTE: Gradient here is projected in the S2^N tangent space
+        self.gradientE.shape = (self.n_images, -1)
+        Gnorms2 = np.sum(self.gradientE**2, axis=1) / self.n_images
+        # Compute the root mean square per image
+        self.gradientENorm[:] = np.sqrt(Gnorms2)
+        self.gradientE.shape = (-1)
+
+        # DEBUG:
+        # print('gradEnorm', self.gradientENorm)
 
         # TODO: we can use a better quadrature such as Gaussian
         # notice that the gradient norm here is using the RMS
-        action = spi.trapezoid(self.gradientENorm, self.path_distances)
+        action = spi.simpson(self.gradientENorm, x=self.path_distances)
 
         # DEBUG:
         # print('action from gradE', action)
@@ -348,9 +349,9 @@ def compute_action(self):
         dist_minus_norm = self.distances[:-1]
         # dY_plus_norm = distances[i];
         # dY_minus_norm = distances[i - 1];
-        springF2 = self.k[1:-1] * ((dist_plus_norm - dist_minus_norm)**2)
+        springF2 = 0.5 * self.k[1:-1] * ((dist_plus_norm - dist_minus_norm)**2)
         # CHECK: do we need to scale?
-        action += np.sum(springF2) / (self.n_images - 2)
+        # action += np.sum(springF2) / (self.n_images - 2)
 
         # DEBUG:
         # print('action spring term', np.sum(springF2) / (self.n_images - 2))
diff --git a/tests/test_two_particles_nebm-fs.py b/tests/test_two_particles_nebm-fs.py
index 7ddcb971..9c46382d 100644
--- a/tests/test_two_particles_nebm-fs.py
+++ b/tests/test_two_particles_nebm-fs.py
@@ -60,7 +60,8 @@ def relax_string(maxst, simname, init_im, interp, save_every=10000):
     # equal to 'the number of initial states specified', minus one.
     interpolations = interp
 
-    nebm = NEBM_FS(sim, init_im, interpolations=interpolations, name=simname)
+    nebm = NEBM_FS(sim, init_im, interpolations=interpolations, name=simname,
+                   interpolation_method='linear')
 
     # dt = integrator.stepsize means after every integrator step, the images
     # are rescaled. We can run more integrator steps if we decrease the
@@ -83,8 +84,8 @@ def mid_m(pos):
 
 def test_energy_barrier_2particles_string():
     # Initial images: we set here a rotation interpolating
-    init_im = [(-1, 0, 0), (1, 0, 0)]
-    interp = [13]
+    init_im = [(-1, 0, 0), (0.0, 0.0, 1.0), (1, 0, 0)]
+    interp = [6, 6]
 
     barriers = []