diff --git a/src/helio/helio_coord.f90 b/src/helio/helio_coord.f90
deleted file mode 100644
index f40781810..000000000
--- a/src/helio/helio_coord.f90
+++ /dev/null
@@ -1,113 +0,0 @@
-submodule (helio_classes) s_helio_coord
-   use swiftest
-contains
-
-   module subroutine helio_coord_vb2vh_pl(self, cb)
-      !! author: David A. Minton
-      !!
-      !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
-      !!
-      !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh.f90 
-      !! Adapted from Hal Levison's Swift routine coord_vb2vh.f 
-      implicit none
-      ! Arguments
-      class(helio_pl),              intent(inout) :: self !! Helio massive body object
-      class(swiftest_cb),           intent(inout) :: cb   !! Swiftest central body object
-      ! Internals
-      integer(I4B)              :: i
-
-      if (self%nbody == 0) return
-
-      associate(pl => self, npl => self%nbody)
-         do i = 1, NDIM
-            cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vb(i, 1:npl)) / cb%Gmass
-            pl%vh(i, 1:npl) = pl%vb(i, 1:npl) - cb%vb(i)
-         end do
-      end associate
-
-      return
-   end subroutine helio_coord_vb2vh_pl
-
-
-   module subroutine helio_coord_vb2vh_tp(self, vbcb)
-      !! author: David A. Minton
-      !!
-      !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
-      !!
-      !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh_tp.f90 
-      !! Adapted from Hal Levison's Swift routine coord_vb2h_tp.f 
-      implicit none
-      ! Arguments
-      class(helio_tp),              intent(inout) :: self !! Helio massive body object
-      real(DP), dimension(:),       intent(in)    :: vbcb  !! Barycentric velocity of the central body
-
-      if (self%nbody == 0) return
-
-      associate(tp => self, ntp => self%nbody)
-         where (tp%lmask(1:ntp))
-            tp%vh(1, 1:ntp) = tp%vb(1, 1:ntp) - vbcb(1)
-            tp%vh(2, 1:ntp) = tp%vb(2, 1:ntp) - vbcb(2)
-            tp%vh(3, 1:ntp) = tp%vb(3, 1:ntp) - vbcb(3)
-         end where
-      end associate
-
-      return
-   end subroutine helio_coord_vb2vh_tp
-
-
-   module subroutine helio_coord_vh2vb_pl(self, cb)
-      !! author: David A. Minton
-      !!
-      !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
-      !!
-      !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb.f90 
-      !! Adapted from Hal Levison's Swift routine coord_vh2b.f 
-      implicit none
-      ! Arguments
-      class(helio_pl),        intent(inout) :: self !! Helio massive body object
-      class(swiftest_cb),     intent(inout) :: cb   !! Swiftest central body object
-      ! Internals
-      integer(I4B)  :: i
-      real(DP)      :: Gmtot
-
-      if (self%nbody == 0) return
-
-      associate(pl => self, npl => self%nbody)
-         Gmtot = cb%Gmass + sum(pl%Gmass(1:npl))
-         do i = 1, NDIM
-            cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vh(i, 1:npl)) / Gmtot
-            pl%vb(i, 1:npl) = pl%vh(i, 1:npl) + cb%vb(i)
-         end do
-      end associate
-
-      return
-   end subroutine helio_coord_vh2vb_pl
-
-
-   module subroutine helio_coord_vh2vb_tp(self, vbcb)
-      !! author: David A. Minton
-      !!
-      !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
-      !!
-      !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb_tp.f90
-      !! Adapted from Hal Levison's Swift routine coord_vh2b_tp.f 
-      implicit none
-      ! Arguments
-      class(helio_tp),        intent(inout) :: self !! Helio massive body object
-      real(DP), dimension(:), intent(in)    :: vbcb !! Barycentric velocity of the central body
-
-      if (self%nbody == 0) return
-
-      associate(tp => self, ntp => self%nbody)
-         where (tp%lmask(1:ntp))
-            tp%vb(1, 1:ntp) = tp%vh(1, 1:ntp) + vbcb(1)
-            tp%vb(2, 1:ntp) = tp%vh(2, 1:ntp) + vbcb(2)
-            tp%vb(3, 1:ntp) = tp%vh(3, 1:ntp) + vbcb(3)
-         end where
-      end associate
-
-      return
-   end subroutine helio_coord_vh2vb_tp
-   
-end submodule s_helio_coord
-
diff --git a/src/helio/helio_setup.f90 b/src/helio/helio_setup.f90
new file mode 100644
index 000000000..2a8c24019
--- /dev/null
+++ b/src/helio/helio_setup.f90
@@ -0,0 +1,22 @@
+submodule(helio_classes) s_helio_setup
+   use swiftest
+contains
+
+   module subroutine helio_setup_initialize_system(self, param)
+      !! author: David A. Minton
+      !!
+      !! Initialize a Helio nbody system from files, converting all heliocentric quantities to barycentric.
+      !!
+      implicit none
+      ! Arguments
+      class(helio_nbody_system),  intent(inout) :: self   !! Helio nbody system object
+      class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
+
+      call whm_setup_initialize_system(self, param)
+      call self%pl%h2b(self%cb)
+      call self%tp%h2b(self%cb)
+
+      return
+   end subroutine helio_setup_initialize_system
+
+end submodule s_helio_setup
\ No newline at end of file
diff --git a/src/io/io.f90 b/src/io/io.f90
index 3d1521680..73e951c76 100644
--- a/src/io/io.f90
+++ b/src/io/io.f90
@@ -40,7 +40,8 @@ module subroutine io_conservation_report(self, param, lterminal)
                open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append", err = 667, iomsg = errmsg)
             end if
          end if
-         call pl%h2b(cb)
+         call pl%vb2vh(cb)
+         call pl%xh2xb(cb)
          call system%get_energy_and_momentum(param) 
          ke_orbit_now = system%ke_orbit
          ke_spin_now = system%ke_spin
diff --git a/src/kick/kick.f90 b/src/kick/kick.f90
index ba3226eb7..c41bf7649 100644
--- a/src/kick/kick.f90
+++ b/src/kick/kick.f90
@@ -14,7 +14,7 @@ module pure subroutine kick_getacch_int_pl(self)
       class(swiftest_pl), intent(inout) :: self
       ! Internals
       integer(I4B)                      :: k
-      real(DP)                          :: rji2, irij3, faci, facj
+      real(DP)                          :: rji2, irij3, faci, facj, rlim2
       real(DP)                          :: dx, dy, dz
 
       associate(pl => self, npl => self%nbody, nplpl => self%nplpl)
@@ -25,15 +25,18 @@ module pure subroutine kick_getacch_int_pl(self)
                   dy = pl%xh(2, j) - pl%xh(2, i)
                   dz = pl%xh(3, j) - pl%xh(3, i)
                   rji2 = dx**2 + dy**2 + dz**2
-                  irij3 = 1.0_DP / (rji2 * sqrt(rji2))
-                  faci = pl%Gmass(i) * irij3
-                  facj = pl%Gmass(j) * irij3
-                  pl%ah(1, i) = pl%ah(1, i) + facj * dx
-                  pl%ah(2, i) = pl%ah(2, i) + facj * dy
-                  pl%ah(3, i) = pl%ah(3, i) + facj * dz
-                  pl%ah(1, j) = pl%ah(1, j) - faci * dx
-                  pl%ah(2, j) = pl%ah(2, j) - faci * dy
-                  pl%ah(3, j) = pl%ah(3, j) - faci * dz
+                  rlim2 = (pl%radius(i) + pl%radius(j))**2
+                  if (rji2 > rlim2) then
+                     irij3 = 1.0_DP / (rji2 * sqrt(rji2))
+                     faci = pl%Gmass(i) * irij3
+                     facj = pl%Gmass(j) * irij3
+                     pl%ah(1, i) = pl%ah(1, i) + facj * dx
+                     pl%ah(2, i) = pl%ah(2, i) + facj * dy
+                     pl%ah(3, i) = pl%ah(3, i) + facj * dz
+                     pl%ah(1, j) = pl%ah(1, j) - faci * dx
+                     pl%ah(2, j) = pl%ah(2, j) - faci * dy
+                     pl%ah(3, j) = pl%ah(3, j) - faci * dz
+                  end if
                end if
             end associate
          end do
diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90
index 89f4aa055..efa326a0f 100644
--- a/src/modules/helio_classes.f90
+++ b/src/modules/helio_classes.f90
@@ -15,7 +15,8 @@ module helio_classes
    !********************************************************************************************************************************
    type, extends(whm_nbody_system) :: helio_nbody_system
    contains
-      procedure :: step => helio_step_system  !! Advance the Helio nbody system forward in time by one step
+      procedure :: step       => helio_step_system             !! Advance the Helio nbody system forward in time by one step
+      procedure :: initialize => helio_setup_initialize_system !! Performs Helio-specific initilization steps, including converting to DH coordinates
    end type helio_nbody_system
 
    !********************************************************************************************************************************
@@ -35,8 +36,6 @@ module helio_classes
    !! Helio massive body particle class
    type, extends(swiftest_pl) :: helio_pl
    contains
-      procedure :: vh2vb       => helio_coord_vh2vb_pl     !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
-      procedure :: vb2vh       => helio_coord_vb2vh_pl     !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
       procedure :: drift       => helio_drift_pl           !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure :: lindrift    => helio_drift_linear_pl    !! Method for linear drift of massive bodies due to barycentric momentum of Sun
       procedure :: accel_gr    => helio_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction
@@ -53,8 +52,6 @@ module helio_classes
    !! Helio test particle class
    type, extends(swiftest_tp) :: helio_tp
    contains
-      procedure :: vh2vb       => helio_coord_vh2vb_tp     !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
-      procedure :: vb2vh       => helio_coord_vb2vh_tp     !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
       procedure :: lindrift    => helio_drift_linear_tp    !! Method for linear drift of massive bodies due to barycentric momentum of Sun
       procedure :: drift       => helio_drift_tp           !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure :: accel_gr    => helio_gr_kick_getacch_tp !! Acceleration term arising from the post-Newtonian correction
@@ -65,32 +62,6 @@ module helio_classes
    end type helio_tp
 
    interface
-      module subroutine helio_coord_vb2vh_pl(self, cb)
-         use swiftest_classes, only : swiftest_cb
-         implicit none
-         class(helio_pl),              intent(inout) :: self !! Helio massive body object
-         class(swiftest_cb),           intent(inout) :: cb   !! Swiftest central body object
-      end subroutine helio_coord_vb2vh_pl
-   
-      module subroutine helio_coord_vb2vh_tp(self, vbcb)
-         implicit none
-         class(helio_tp),              intent(inout) :: self !! Helio massive body object
-         real(DP), dimension(:),       intent(in)    :: vbcb  !! Barycentric velocity of the central body
-      end subroutine helio_coord_vb2vh_tp
-   
-      module subroutine helio_coord_vh2vb_pl(self, cb)
-         use swiftest_classes, only : swiftest_cb
-         implicit none
-         class(helio_pl),        intent(inout) :: self !! Helio massive body object
-         class(swiftest_cb),     intent(inout) :: cb   !! Swiftest central body object
-      end subroutine helio_coord_vh2vb_pl
-   
-      module subroutine helio_coord_vh2vb_tp(self, vbcb)
-         implicit none
-         class(helio_tp),        intent(inout) :: self !! Helio massive body object
-         real(DP), dimension(:), intent(in)    :: vbcb !! Barycentric velocity of the central body
-      end subroutine helio_coord_vh2vb_tp
-
       module subroutine helio_drift_body(self, system, param, dt)
          use swiftest_classes, only : swiftest_body, swiftest_nbody_system, swiftest_parameters
          implicit none
@@ -206,6 +177,13 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, lbeg)
          logical,                      intent(in)    :: lbeg   !! Logical flag indicating whether this is the beginning of the half step or not. 
       end subroutine helio_kick_vb_tp
 
+      module subroutine helio_setup_initialize_system(self, param)
+         use swiftest_classes, only : swiftest_parameters
+         implicit none
+         class(helio_nbody_system),  intent(inout) :: self   !! Helio nbody system object
+         class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
+      end subroutine helio_setup_initialize_system
+
       module subroutine helio_step_pl(self, system, param, t, dt)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
          implicit none
diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90
index 4cf180e84..f255e1a38 100644
--- a/src/modules/swiftest_classes.f90
+++ b/src/modules/swiftest_classes.f90
@@ -147,8 +147,8 @@ module swiftest_classes
       logical,               dimension(:),   allocatable :: ldiscard        !! Body should be discarded
       logical,               dimension(:),   allocatable :: lmask           !! Logical mask used to select a subset of bodies when performing certain operations (drift, kick, accel, etc.)
       real(DP),              dimension(:),   allocatable :: mu              !! G * (Mcb + [m])
-      real(DP),              dimension(:,:), allocatable :: xh              !! Heliocentric position
-      real(DP),              dimension(:,:), allocatable :: vh              !! Heliocentric velocity
+      real(DP),              dimension(:,:), allocatable :: xh              !! Swiftestcentric position
+      real(DP),              dimension(:,:), allocatable :: vh              !! Swiftestcentric velocity
       real(DP),              dimension(:,:), allocatable :: xb              !! Barycentric position
       real(DP),              dimension(:,:), allocatable :: vb              !! Barycentric velocity
       real(DP),              dimension(:,:), allocatable :: ah              !! Total heliocentric acceleration
@@ -226,7 +226,9 @@ module swiftest_classes
       procedure :: append       => util_append_pl         !! Appends elements from one structure to another
       procedure :: h2b          => util_coord_h2b_pl      !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity)
       procedure :: b2h          => util_coord_b2h_pl      !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity)
-      procedure :: xh2xb        => util_coord_xh2xb_pl    !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position and velocity)
+      procedure :: vh2vb        => util_coord_vh2vb_pl    !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
+      procedure :: vb2vh        => util_coord_vb2vh_pl    !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
+      procedure :: xh2xb        => util_coord_xh2xb_pl    !! Convert massive bodies from heliocentric to barycentric coordinates (position only)
       procedure :: fill         => util_fill_pl           !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
       procedure :: resize       => util_resize_pl         !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small.
       procedure :: set_beg_end  => util_set_beg_end_pl    !! Sets the beginning and ending positions and velocities of planets.
@@ -258,6 +260,8 @@ module swiftest_classes
       procedure :: append    => util_append_tp         !! Appends elements from one structure to another
       procedure :: h2b       => util_coord_h2b_tp      !! Convert test particles from heliocentric to barycentric coordinates (position and velocity)
       procedure :: b2h       => util_coord_b2h_tp      !! Convert test particles from barycentric to heliocentric coordinates (position and velocity)
+      procedure :: vh2vb     => util_coord_vh2vb_tp    !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
+      procedure :: vb2vh     => util_coord_vb2vh_tp    !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
       procedure :: fill      => util_fill_tp           !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
       procedure :: get_peri  => util_peri_tp           !! Determine system pericenter passages for test particles 
       procedure :: resize    => util_resize_tp         !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small.
@@ -516,9 +520,9 @@ module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh)
          implicit none
          class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
          real(DP),                   intent(in)  :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
-         real(DP), dimension(:),     intent(in)  :: xh    !! Heliocentric position vector 
+         real(DP), dimension(:),     intent(in)  :: xh    !! Swiftestcentric position vector 
          real(DP), dimension(:),     intent(in)  :: pv    !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
-         real(DP), dimension(:),     intent(out) :: vh    !! Heliocentric velocity vector 
+         real(DP), dimension(:),     intent(out) :: vh    !! Swiftestcentric velocity vector 
       end subroutine gr_pseudovel2vel
 
       module pure subroutine gr_pv2vh_body(self, param)
@@ -531,8 +535,8 @@ module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
          implicit none
          class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
          real(DP),                   intent(in)  :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
-         real(DP), dimension(:),     intent(in)  :: xh    !! Heliocentric position vector 
-         real(DP), dimension(:),     intent(in)  :: vh    !! Heliocentric velocity vector 
+         real(DP), dimension(:),     intent(in)  :: xh    !! Swiftestcentric position vector 
+         real(DP), dimension(:),     intent(in)  :: vh    !! Swiftestcentric velocity vector 
          real(DP), dimension(:),     intent(out) :: pv    !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
       end subroutine gr_vel2pseudovel
 
@@ -687,8 +691,8 @@ module subroutine io_write_frame_encounter(iu, t, id1, id2, Gmass1, Gmass2, radi
          integer(I4B),           intent(in) :: id1, id2         !! ids of the two encountering bodies
          real(DP),               intent(in) :: Gmass1, Gmass2   !! G*mass of the two encountering bodies
          real(DP),               intent(in) :: radius1, radius2 !! Radii of the two encountering bodies
-         real(DP), dimension(:), intent(in) :: xh1, xh2         !! Heliocentric position vectors of the two encountering bodies 
-         real(DP), dimension(:), intent(in) :: vh1, vh2         !! Heliocentric velocity vectors of the two encountering bodies 
+         real(DP), dimension(:), intent(in) :: xh1, xh2         !! Swiftestcentric position vectors of the two encountering bodies 
+         real(DP), dimension(:), intent(in) :: vh1, vh2         !! Swiftestcentric velocity vectors of the two encountering bodies 
       end subroutine io_write_frame_encounter
 
       module subroutine io_write_frame_system(self, iu, param)
@@ -937,12 +941,42 @@ module subroutine util_coord_h2b_tp(self, cb)
          class(swiftest_cb), intent(in)    :: cb   !! Swiftest central body object
       end subroutine util_coord_h2b_tp
 
+      module subroutine util_coord_vb2vh_pl(self, cb)
+         implicit none
+         class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+         class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
+      end subroutine util_coord_vb2vh_pl
+   
+      module subroutine util_coord_vb2vh_tp(self, vbcb)
+         implicit none
+         class(swiftest_tp),     intent(inout) :: self !! Swiftest test particle object
+         real(DP), dimension(:), intent(in)    :: vbcb !! Barycentric velocity of the central body
+      end subroutine util_coord_vb2vh_tp
+   
+      module subroutine util_coord_vh2vb_pl(self, cb)
+         implicit none
+         class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+         class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
+      end subroutine util_coord_vh2vb_pl
+   
+      module subroutine util_coord_vh2vb_tp(self, vbcb)
+         implicit none
+         class(swiftest_tp),        intent(inout) :: self !! Swiftest test particle object
+         real(DP), dimension(:), intent(in)    :: vbcb !! Barycentric velocity of the central body
+      end subroutine util_coord_vh2vb_tp
+
       module subroutine util_coord_xh2xb_pl(self, cb)
          implicit none
          class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
          class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
       end subroutine util_coord_xh2xb_pl
 
+      module subroutine util_coord_xh2xb_tp(self, cb)
+         implicit none
+         class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
+         class(swiftest_cb), intent(in) :: cb      !! Swiftest central body object
+      end subroutine util_coord_xh2xb_tp
+
       module subroutine util_copy_encounter(self, source)
          implicit none
          class(swiftest_encounter), intent(inout) :: self   !! Encounter list 
diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90
index 9ba20d243..032ccd47a 100644
--- a/src/modules/symba_classes.f90
+++ b/src/modules/symba_classes.f90
@@ -225,20 +225,24 @@ module subroutine symba_collision_resolve_mergers(self, system, param)
          class(symba_parameters),   intent(in)    :: param  !! Current run configuration parameters with SyMBA additions
       end subroutine symba_collision_resolve_mergers
 
-      module subroutine symba_collision_resolve_plplenc(self, system, param, t)
+      module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec)
          implicit none
          class(symba_plplenc),       intent(inout) :: self   !! SyMBA pl-pl encounter list
          class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
          class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
          real(DP),                   intent(in)    :: t      !! Current simulation time
+         real(DP),                   intent(in)    :: dt     !! Current simulation step size
+         integer(I4B),               intent(in)    :: irec   !! Current recursion level
       end subroutine symba_collision_resolve_plplenc
    
-      module subroutine symba_collision_resolve_pltpenc(self, system, param, t)
+      module subroutine symba_collision_resolve_pltpenc(self, system, param, t, dt, irec)
          implicit none
          class(symba_pltpenc),       intent(inout) :: self   !! SyMBA pl-tp encounter list
          class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
          class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
          real(DP),                   intent(in)    :: t      !! Current simulation time
+         real(DP),                   intent(in)    :: dt     !! Current simulation step size
+         integer(I4B),               intent(in)    :: irec   !! Current recursion level
       end subroutine symba_collision_resolve_pltpenc
 
       module subroutine symba_discard_pl(self, system, param)
diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90
index 6d5b26394..0b9caebcc 100644
--- a/src/modules/whm_classes.f90
+++ b/src/modules/whm_classes.f90
@@ -73,8 +73,8 @@ module whm_classes
    type, extends(swiftest_nbody_system) :: whm_nbody_system
    contains
       !> Replace the abstract procedures with concrete ones
-      procedure :: initialize   => whm_setup_initialize_system      !! Performs WHM-specific initilization steps, like calculating the Jacobi masses
-      procedure :: step         => whm_step_system       !! Advance the WHM nbody system forward in time by one step
+      procedure :: initialize   => whm_setup_initialize_system ! Performs WHM-specific initilization steps, like calculating the Jacobi masses
+      procedure :: step         => whm_step_system             !! Advance the WHM nbody system forward in time by one step
    end type whm_nbody_system
 
    interface
diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90
index f1f42d3cf..2f93881b1 100644
--- a/src/symba/symba_collision.f90
+++ b/src/symba/symba_collision.f90
@@ -1054,7 +1054,7 @@ module subroutine symba_collision_resolve_mergers(self, system, param)
    end subroutine symba_collision_resolve_mergers
 
 
-   module subroutine symba_collision_resolve_plplenc(self, system, param, t)
+   module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec)
       !! author: David A. Minton
       !! 
       !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the collision
@@ -1065,8 +1065,11 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t)
       class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
       class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
       real(DP),                   intent(in)    :: t      !! Current simulation time
+      real(DP),                   intent(in)    :: dt     !! Current simulation step size
+      integer(I4B),               intent(in)    :: irec   !! Current recursion level
       ! Internals
       real(DP) :: Eorbit_before, Eorbit_after
+      logical :: lplpl_collision
    
       associate(plplenc_list => self, plplcollision_list => system%plplcollision_list)
          select type(pl => system%pl)
@@ -1074,26 +1077,41 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t)
             select type(param)
             class is (symba_parameters)
                if (plplcollision_list%nenc == 0) return ! No collisions to resolve
-
-               write(*, *) "Collision between massive bodies detected at time t = ", t
-               call pl%vb2vh(system%cb)
+               ! Make sure that the heliocentric and barycentric coordinates are consistent with each other
+               call pl%vb2vh(system%cb) 
                call pl%xh2xb(system%cb)
-               if (param%lfragmentation) then
-                  call plplcollision_list%resolve_fragmentations(system, param)
-               else
-                  call plplcollision_list%resolve_mergers(system, param)
-               end if
-
-               ! Destroy the collision list now that the collisions are resolved
-               call plplcollision_list%setup(0)
-
+   
                ! Get the energy before the collision is resolved
                if (param%lenergy) then
                   call system%get_energy_and_momentum(param)
                   Eorbit_before = system%te
                end if
 
-               call pl%rearray(system, param)
+               do
+                  write(*, *) "Collision between massive bodies detected at time t = ", t
+                  if (param%lfragmentation) then
+                     call plplcollision_list%resolve_fragmentations(system, param)
+                  else
+                     call plplcollision_list%resolve_mergers(system, param)
+                  end if
+   
+                  ! Destroy the collision list now that the collisions are resolved
+                  call plplcollision_list%setup(0)
+
+                  ! Save the add/discard information to file
+                  call system%write_discard(param)
+
+                  ! Rearrange the arrays: Remove discarded bodies, add any new bodies, resort, and recompute all indices and encounter lists
+                  call pl%rearray(system, param)
+
+                  ! Destroy the add/discard list so that we don't append the same body multiple times if another collision is detected
+                  call system%pl_discards%setup(0, param)
+                  call system%pl_adds%setup(0, param)
+
+                  ! Check whether or not any of the particles that were just added are themselves in a collision state. 
+                  lplpl_collision = plplenc_list%collision_check(system, param, t, dt, irec)
+                  if (.not.lplpl_collision) exit
+               end do
 
                if (param%lenergy) then
                   call system%get_energy_and_momentum(param)
@@ -1109,7 +1127,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t)
    end subroutine symba_collision_resolve_plplenc
 
 
-   module subroutine symba_collision_resolve_pltpenc(self, system, param, t)
+   module subroutine symba_collision_resolve_pltpenc(self, system, param, t, dt, irec)
       !! author: David A. Minton
       !! 
       !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the collision
@@ -1120,6 +1138,8 @@ module subroutine symba_collision_resolve_pltpenc(self, system, param, t)
       class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
       class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
       real(DP),                   intent(in)    :: t      !! Current simulation tim
+      real(DP),                   intent(in)    :: dt     !! Current simulation step size
+      integer(I4B),               intent(in)    :: irec   !! Current recursion level
   
       call system%tp%discard(system, param)
 
diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90
index 307c37675..67f8063b5 100644
--- a/src/symba/symba_setup.f90
+++ b/src/symba/symba_setup.f90
@@ -73,7 +73,7 @@ module subroutine symba_setup_initialize_system(self, param)
 
       ! Call parent method
       associate(system => self)
-         call whm_setup_initialize_system(system, param)
+         call helio_setup_initialize_system(system, param)
          call system%pltpenc_list%setup(0)
          call system%plplenc_list%setup(0)
          call system%plplcollision_list%setup(0)
diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90
index ca5aa43cf..e07d16cc2 100644
--- a/src/symba/symba_step.f90
+++ b/src/symba/symba_step.f90
@@ -197,8 +197,8 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci)
                      lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) 
                      lpltp_collision = pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci) 
 
-                     if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl)
-                     if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl)
+                     if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci)
+                     if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci)
                   end if
 
                   call self%set_recur_levels(ireci)
diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90
index 4daad0e53..b4a77231c 100644
--- a/src/util/util_coord.f90
+++ b/src/util/util_coord.f90
@@ -42,46 +42,10 @@ module subroutine util_coord_h2b_pl(self, cb)
    end subroutine util_coord_h2b_pl
 
 
-   module subroutine util_coord_xh2xb_pl(self, cb)
-      !! author: David A. Minton
-      !!
-      !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position and velocity)
-      !!
-      !! Adapted from David E. Kaufmann's Swifter routine coord_h2b.f90 
-      !! Adapted from Hal Levison's Swift routine coord_h2b.f 
-      implicit none
-      ! Arguments
-      class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
-      class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
-      ! Internals
-      integer(I4B)  :: i
-      real(DP)      :: Gmtot
-      real(DP), dimension(NDIM) :: xtmp
-
-      if (self%nbody == 0) return
-      associate(pl => self, npl => self%nbody)
-         Gmtot = cb%Gmass
-         xtmp(:) = 0.0_DP
-         do i = 1, npl
-            if (pl%status(i) == INACTIVE) cycle
-            Gmtot = Gmtot + pl%Gmass(i)
-            xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i)
-         end do
-         cb%xb(:) = -xtmp(:) / Gmtot
-         do i = 1, npl
-            if (pl%status(i) == INACTIVE) cycle
-            pl%xb(:,i) = pl%xh(:,i) + cb%xb(:)
-         end do
-      end associate
-
-      return
-   end subroutine util_coord_xh2xb_pl
-
-
    module subroutine util_coord_h2b_tp(self, cb)
       !! author: David A. Minton
       !!
-      !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity)
+      !! Convert test particles from heliocentric to barycentric coordinates (position and velocity)
       !!
       !! Adapted from David E. Kaufmann's Swifter routine coord_h2b_tp.f90 
       !! Adapted from Hal Levison's Swift routine coord_h2b_tp.f 
@@ -156,5 +120,174 @@ module subroutine util_coord_b2h_tp(self, cb)
 
       return
    end subroutine util_coord_b2h_tp
+
+
+   module subroutine util_coord_vb2vh_pl(self, cb)
+      !! author: David A. Minton
+      !!
+      !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh.f90 
+      !! Adapted from Hal Levison's Swift routine coord_vb2vh.f 
+      implicit none
+      ! Arguments
+      class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+      class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
+      ! Internals
+      integer(I4B)              :: i
+
+      if (self%nbody == 0) return
+
+      associate(pl => self, npl => self%nbody)
+         do i = 1, NDIM
+            cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vb(i, 1:npl)) / cb%Gmass
+            pl%vh(i, 1:npl) = pl%vb(i, 1:npl) - cb%vb(i)
+         end do
+      end associate
+
+      return
+   end subroutine util_coord_vb2vh_pl
+
+
+   module subroutine util_coord_vb2vh_tp(self, vbcb)
+      !! author: David A. Minton
+      !!
+      !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh_tp.f90 
+      !! Adapted from Hal Levison's Swift routine coord_vb2h_tp.f 
+      implicit none
+      ! Arguments
+      class(swiftest_tp),     intent(inout) :: self !! Swiftest test particle object
+      real(DP), dimension(:), intent(in)    :: vbcb !! Barycentric velocity of the central body
+
+      if (self%nbody == 0) return
+
+      associate(tp => self, ntp => self%nbody)
+         where (tp%lmask(1:ntp))
+            tp%vh(1, 1:ntp) = tp%vb(1, 1:ntp) - vbcb(1)
+            tp%vh(2, 1:ntp) = tp%vb(2, 1:ntp) - vbcb(2)
+            tp%vh(3, 1:ntp) = tp%vb(3, 1:ntp) - vbcb(3)
+         end where
+      end associate
+
+      return
+   end subroutine util_coord_vb2vh_tp
+
+
+   module subroutine util_coord_vh2vb_pl(self, cb)
+      !! author: David A. Minton
+      !!
+      !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb.f90 
+      !! Adapted from Hal Levison's Swift routine coord_vh2b.f 
+      implicit none
+      ! Arguments
+      class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+      class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
+      ! Internals
+      integer(I4B)  :: i
+      real(DP)      :: Gmtot
+
+      if (self%nbody == 0) return
+
+      associate(pl => self, npl => self%nbody)
+         Gmtot = cb%Gmass + sum(pl%Gmass(1:npl))
+         do i = 1, NDIM
+            cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vh(i, 1:npl)) / Gmtot
+            pl%vb(i, 1:npl) = pl%vh(i, 1:npl) + cb%vb(i)
+         end do
+      end associate
+
+      return
+   end subroutine util_coord_vh2vb_pl
+
+
+   module subroutine util_coord_vh2vb_tp(self, vbcb)
+      !! author: David A. Minton
+      !!
+      !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb_tp.f90
+      !! Adapted from Hal Levison's Swift routine coord_vh2b_tp.f 
+      implicit none
+      ! Arguments
+      class(swiftest_tp),     intent(inout) :: self !! Swiftest test particle object
+      real(DP), dimension(:), intent(in)    :: vbcb !! Barycentric velocity of the central body
+
+      if (self%nbody == 0) return
+
+      associate(tp => self, ntp => self%nbody)
+         where (tp%lmask(1:ntp))
+            tp%vb(1, 1:ntp) = tp%vh(1, 1:ntp) + vbcb(1)
+            tp%vb(2, 1:ntp) = tp%vh(2, 1:ntp) + vbcb(2)
+            tp%vb(3, 1:ntp) = tp%vh(3, 1:ntp) + vbcb(3)
+         end where
+      end associate
+
+      return
+   end subroutine util_coord_vh2vb_tp
    
+
+   module subroutine util_coord_xh2xb_pl(self, cb)
+      !! author: David A. Minton
+      !!
+      !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position only)
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine coord_h2b.f90 
+      !! Adapted from Hal Levison's Swift routine coord_h2b.f 
+      implicit none
+      ! Arguments
+      class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+      class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
+      ! Internals
+      integer(I4B)  :: i
+      real(DP)      :: Gmtot
+      real(DP), dimension(NDIM) :: xtmp
+
+      if (self%nbody == 0) return
+      associate(pl => self, npl => self%nbody)
+         Gmtot = cb%Gmass
+         xtmp(:) = 0.0_DP
+         do i = 1, npl
+            if (pl%status(i) == INACTIVE) cycle
+            Gmtot = Gmtot + pl%Gmass(i)
+            xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i)
+         end do
+         cb%xb(:) = -xtmp(:) / Gmtot
+         do i = 1, npl
+            if (pl%status(i) == INACTIVE) cycle
+            pl%xb(:,i) = pl%xh(:,i) + cb%xb(:)
+         end do
+      end associate
+
+      return
+   end subroutine util_coord_xh2xb_pl
+
+
+   module subroutine util_coord_xh2xb_tp(self, cb)
+      !! author: David A. Minton
+      !!
+      !! Convert test particles from heliocentric to barycentric coordinates (position only)
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine coord_h2b_tp.f90 
+      !! Adapted from Hal Levison's Swift routine coord_h2b_tp.f 
+      implicit none
+      ! Arguments
+      class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
+      class(swiftest_cb), intent(in) :: cb      !! Swiftest central body object
+      ! Internals
+      integer(I4B) :: i
+
+      if (self%nbody == 0) return
+      associate(tp => self, ntp => self%nbody)
+         do concurrent (i = 1:ntp, tp%status(i) /= INACTIVE)
+            tp%xb(:, i) = tp%xh(:, i) + cb%xb(:)
+         end do
+      end associate
+
+      return
+   end subroutine util_coord_xh2xb_tp
+
 end submodule s_util_coord
\ No newline at end of file
diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90
index 7535ea67b..9a4060dc1 100644
--- a/src/util/util_get_energy_momentum.f90
+++ b/src/util/util_get_energy_momentum.f90
@@ -93,7 +93,7 @@ module subroutine util_get_energy_momentum_system(self, param)
             end do
 
             do concurrent (k = 1:pl%nplpl, lstatpl(k))
-               pepl(k) = -pl%Gmass(indi(k)) * pl%mass(indj(k)) / norm2(pl%xb(:, indi(k)) - pl%xb(:, indj(k))) 
+               pepl(k) = -(pl%Gmass(indi(k)) * pl%mass(indj(k))) / norm2(pl%xb(:, indi(k)) - pl%xb(:, indj(k)))
             end do
          end associate
    
diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90
index eaed16c14..ad44fbefb 100644
--- a/src/whm/whm_setup.f90
+++ b/src/whm/whm_setup.f90
@@ -73,7 +73,7 @@ module subroutine whm_setup_initialize_system(self, param)
       !!
       implicit none
       ! Arguments
-      class(whm_nbody_system),    intent(inout) :: self    !! Swiftest system object
+      class(whm_nbody_system),    intent(inout) :: self   !! WHM nbody system object
       class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
 
       call setup_initialize_system(self, param)