diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp
index 0fc42aeb..fe69bf0a 100644
--- a/src/CabanaPD_Particles.hpp
+++ b/src/CabanaPD_Particles.hpp
@@ -112,6 +112,7 @@ class Particles<DeviceType, PMB, Dimension>
     // FIXME: enable variable aosoa.
     using aosoa_x_type = Cabana::AoSoA<vector_type, memory_space, 1>;
     using aosoa_u_type = Cabana::AoSoA<vector_type, memory_space, 1>;
+    using aosoa_y_type = Cabana::AoSoA<vector_type, memory_space, 1>;
     using aosoa_f_type = Cabana::AoSoA<vector_type, memory_space, 1>;
     using aosoa_vol_type = Cabana::AoSoA<scalar_type, memory_space, 1>;
     using aosoa_nofail_type = Cabana::AoSoA<int_type, memory_space, 1>;
@@ -211,12 +212,15 @@ class Particles<DeviceType, PMB, Dimension>
         int num_particles = particles_per_cell * owned_cells.size();
         resize( num_particles, 0 );
 
+        // It could be good to keep a consistent order of variables along the
+        // function
         auto x = slice_x();
         auto v = slice_v();
         auto f = slice_f();
         auto type = slice_type();
         auto rho = slice_rho();
         auto u = slice_u();
+        auto y = slice_y();
         auto vol = slice_vol();
         auto nofail = slice_nofail();
 
@@ -251,6 +255,7 @@ class Particles<DeviceType, PMB, Dimension>
                 {
                     x( pid, d ) = cell_coord[d];
                     u( pid, d ) = 0.0;
+                    y( pid, d ) = 0.0;
                     v( pid, d ) = 0.0;
                     f( pid, d ) = 0.0;
                 }
@@ -293,13 +298,24 @@ class Particles<DeviceType, PMB, Dimension>
             KOKKOS_LAMBDA( const int pid ) { init_functor( pid ); } );
     }
 
-    auto slice_x() { return Cabana::slice<0>( _aosoa_x, "positions" ); }
-    auto slice_x() const { return Cabana::slice<0>( _aosoa_x, "positions" ); }
+    auto slice_x()
+    {
+        return Cabana::slice<0>( _aosoa_x, "reference positions" );
+    }
+    auto slice_x() const
+    {
+        return Cabana::slice<0>( _aosoa_x, "reference positions" );
+    }
     auto slice_u() { return Cabana::slice<0>( _aosoa_u, "displacements" ); }
     auto slice_u() const
     {
         return Cabana::slice<0>( _aosoa_u, "displacements" );
     }
+    auto slice_y() { return Cabana::slice<0>( _aosoa_y, "current positions" ); }
+    auto slice_y() const
+    {
+        return Cabana::slice<0>( _aosoa_y, "current positions" );
+    }
     auto slice_f() { return Cabana::slice<0>( _aosoa_f, "forces" ); }
     auto slice_f_a()
     {
@@ -349,6 +365,8 @@ class Particles<DeviceType, PMB, Dimension>
 
         _aosoa_x.resize( new_local + new_ghost );
         _aosoa_u.resize( new_local + new_ghost );
+        _aosoa_y.resize( new_local + new_ghost );
+        // What about resizing velocities?
         _aosoa_vol.resize( new_local + new_ghost );
         _aosoa_f.resize( new_local );
         _aosoa_other.resize( new_local );
@@ -356,19 +374,28 @@ class Particles<DeviceType, PMB, Dimension>
         size = _aosoa_x.size();
     };
 
-    void output( const int output_step, const double output_time )
+    auto getPosition( const bool use_reference )
+    {
+        if ( use_reference )
+            return slice_x();
+        else
+            return slice_y();
+    }
+
+    void output( const int output_step, const double output_time,
+                 const bool use_reference = true )
     {
 #ifdef Cabana_ENABLE_HDF5
         Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
             h5_config, "particles", MPI_COMM_WORLD, output_step, output_time,
-            n_local, slice_x(), slice_W(), slice_f(), slice_u(), slice_v(),
-            slice_phi() );
+            n_local, getPosition( use_reference ), slice_W(), slice_f(),
+            slice_u(), slice_v(), slice_phi() );
 #else
 #ifdef Cabana_ENABLE_SILO
         Cajita::Experimental::SiloParticleOutput::writePartialRangeTimeStep(
             "particles", local_grid->globalGrid(), output_step, output_time, 0,
-            n_local, slice_x(), slice_W(), slice_f(), slice_u(), slice_v(),
-            slice_phi() );
+            n_local, getPosition( use_reference ), slice_W(), slice_f(),
+            slice_u(), slice_v(), slice_phi() );
 #else
         log( std::cout, "No particle output enabled for step ", output_step,
              " (", output_time, ")" );
@@ -381,6 +408,8 @@ class Particles<DeviceType, PMB, Dimension>
   protected:
     aosoa_x_type _aosoa_x;
     aosoa_u_type _aosoa_u;
+    aosoa_y_type _aosoa_y;
+    // What about v and other variables, like W, rho, ...?
     aosoa_f_type _aosoa_f;
     aosoa_vol_type _aosoa_vol;
     aosoa_nofail_type _aosoa_nofail;
@@ -484,21 +513,24 @@ class Particles<DeviceType, LPS, Dimension>
         _aosoa_m.resize( new_local + new_ghost );
     }
 
-    void output( const int output_step, const double output_time )
+    void output( const int output_step, const double output_time,
+                 const bool use_reference = true )
     {
 #ifdef Cabana_ENABLE_HDF5
         Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
             h5_config, "particles", MPI_COMM_WORLD, output_step, output_time,
-            n_local, base_type::slice_x(), base_type::slice_W(),
-            base_type::slice_f(), base_type::slice_u(), base_type::slice_v(),
-            base_type::slice_phi(), slice_m(), slice_theta() );
+            n_local, base_type::getPosition( use_reference ),
+            base_type::slice_W(), base_type::slice_f(), base_type::slice_u(),
+            base_type::slice_v(), base_type::slice_phi(), slice_m(),
+            slice_theta() );
 #else
 #ifdef Cabana_ENABLE_SILO
         Cajita::Experimental::SiloParticleOutput::writePartialRangeTimeStep(
             "particles", local_grid->globalGrid(), output_step, output_time, 0,
-            n_local, base_type::slice_x(), base_type::slice_W(),
-            base_type::slice_f(), base_type::slice_u(), base_type::slice_v(),
-            base_type::slice_phi(), slice_m(), slice_theta() );
+            n_local, base_type::getPosition( use_reference ),
+            base_type::slice_W(), base_type::slice_f(), base_type::slice_u(),
+            base_type::slice_v(), base_type::slice_phi(), slice_m(),
+            slice_theta() );
 #else
         log( std::cout, "No particle output enabled for ", output_step, "(",
              output_time, ")" );