diff --git a/src/basic/lattice_vectors.F90 b/src/basic/lattice_vectors.F90
index fedafbca881bcf42642de8faffaf7cf27ec3bfd2..9513afcdfd7ca5753642e780df24ac5b3c29c30f 100644
--- a/src/basic/lattice_vectors.F90
+++ b/src/basic/lattice_vectors.F90
@@ -403,12 +403,16 @@ contains
   end function lattice_vectors_red_to_cart
 
   !--------------------------------------------------------------
-  function lattice_vectors_fold_into_cell(this, xx) result(new_xx)
-    class(lattice_vectors_t), intent(in) :: this
-    FLOAT,                    intent(in) :: xx(this%space%dim)
+  function lattice_vectors_fold_into_cell(this, xx, red_translation_vector) result(new_xx)
+    class(lattice_vectors_t), intent(in)  :: this
+    FLOAT,                    intent(in)  :: xx(this%space%dim)
+    integer, optional,        intent(out) :: red_translation_vector(this%space%dim)
     FLOAT :: new_xx(this%space%dim)
 
-    integer :: idir
+    integer :: idir, red_translation_vector_(this%space%dim)
+
+    ! This vector represents the translation vector used to fold back into the unit cell
+    red_translation_vector_ = 0
 
     if (this%space%is_periodic()) then
       ! Convert the position to reduced coordinates
@@ -416,20 +420,22 @@ contains
 
       do idir = 1, this%space%periodic_dim
         ! Change of origin
-        new_xx(idir) = new_xx(idir) + M_HALF
+        new_xx(idir) = new_xx(idir) + 0.5_real64
 
         ! Fold into cell
-        new_xx(idir) = new_xx(idir) - anint(new_xx(idir))
+        red_translation_vector_(idir) = - int(anint(new_xx(idir)))
+        new_xx(idir) = new_xx(idir) + red_translation_vector_(idir)
         if (new_xx(idir) < -1.0e-6_real64) then
-          new_xx(idir) = new_xx(idir) + M_ONE
+          new_xx(idir) = new_xx(idir) + 1.0_real64
+          red_translation_vector_(idir) = red_translation_vector_(idir) + 1
         end if
 
         ! Sanity checks
         ASSERT(new_xx(idir) >= -1.0e-6_real64)
-        ASSERT(new_xx(idir) < M_ONE)
+        ASSERT(new_xx(idir) < 1.0_real64)
 
         ! Change origin back
-        new_xx(idir) = new_xx(idir) - M_HALF
+        new_xx(idir) = new_xx(idir) - 0.5_real64
       end do
 
       ! Convert back to Cartesian coordinates
@@ -438,6 +444,10 @@ contains
       new_xx = xx
     end if
 
+    if (present(red_translation_vector)) then
+      red_translation_vector = red_translation_vector_
+    end if
+
   end function lattice_vectors_fold_into_cell
 
   !--------------------------------------------------------------