diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90
index 270f48a7a..688ac1a6e 100644
--- a/src/fragmentation/fragmentation.f90
+++ b/src/fragmentation/fragmentation.f90
@@ -167,7 +167,7 @@ subroutine set_scale_factors()
          ! Set scale factors
          Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2)))
          rscale = sum(radius(:))
-         mscale = sqrt(Escale * rscale) / param%GU
+         mscale = sqrt(Escale * rscale / param%GU) 
          vscale = sqrt(Escale / mscale) 
          tscale = rscale / vscale 
          Lscale = mscale * rscale * vscale
@@ -341,11 +341,25 @@ subroutine calculate_system_energy(linclude_fragments)
             ltmp(1:npl) = .true.
             call tmpsys%pl%setup(npl_new, param)
             call tmpsys%pl%fill(pl, ltmp)
+            tmpsys%pl%mass(1:npl) = tmpsys%pl%mass(1:npl) / mscale
+            tmpsys%pl%Gmass(1:npl) = tmpsys%pl%Gmass(1:npl) * tscale**2 / rscale**3
+            tmpsys%cb%mass = tmpsys%cb%mass / mscale
+            tmpsys%cb%Gmass = tmpsys%cb%Gmass * tscale**2 / rscale**3
+            tmpsys%cb%radius = tmpsys%cb%radius / rscale
+            tmpsys%pl%radius(1:npl) = tmpsys%pl%radius(1:npl) / rscale
+            tmpsys%pl%xh(:,1:npl) = tmpsys%pl%xh(:,1:npl) / rscale
+            tmpsys%pl%vh(:,1:npl) = tmpsys%pl%vh(:,1:npl) / vscale
+            tmpsys%pl%xb(:,1:npl) = tmpsys%pl%xb(:,1:npl) / rscale
+            tmpsys%pl%vb(:,1:npl) = tmpsys%pl%vb(:,1:npl) / vscale
+            tmpsys%pl%rot(:,1:npl) = tmpsys%pl%rot(:,1:npl) * tscale
+            tmpsys%cb%xb(:) = tmpsys%cb%xb(:) / rscale
+            tmpsys%cb%vb(:) = tmpsys%cb%vb(:) / vscale
+
       
             if (linclude_fragments) then ! Append the fragments if they are included
                ! Energy calculation requires the fragments to be in the system barcyentric frame, s
                tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag)
-               tmpsys%pl%Gmass(npl+1:npl_new) = param%GU * m_frag(1:nfrag)
+               tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * rscale**3 * mscale / tscale**2 
                tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag)
                tmpsys%pl%xb(:,npl+1:npl_new) =  xb_frag(:,1:nfrag)
                tmpsys%pl%vb(:,npl+1:npl_new) =  vb_frag(:,1:nfrag)
diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90
index 4ec9c6c36..05c3e2436 100644
--- a/src/util/util_get_energy_momentum.f90
+++ b/src/util/util_get_energy_momentum.f90
@@ -76,14 +76,14 @@ module subroutine util_get_energy_momentum_system(self, param)
          !$omp simd 
          do i = 1, npl
             associate(px => pl%xh(1,i), py => pl%xh(2,i), pz => pl%xh(3,i))
-               pecb(i) = -param%GU * cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2)
+               pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2)
             end associate
          end do
    
          ! Do the potential energy between pairs of massive bodies
          do k = 1, pl%nplpl
             associate(ik => pl%k_plpl(1, k), jk => pl%k_plpl(2, k))
-               pepl(k) = -param%GU * pl%mass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) 
+               pepl(k) = -pl%Gmass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) 
                lstatpl(k) = (lstatus(ik) .and. lstatus(jk))
             end associate
          end do
@@ -99,8 +99,8 @@ module subroutine util_get_energy_momentum_system(self, param)
             do i = 1, npl
                irh(i) = 1.0_DP / norm2(pl%xh(:,i))
             end do
-            call obl_pot(npl, cb%mass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot)
-            system%pe = system%pe + param%GU * oblpot
+            call obl_pot(npl, cb%Gmass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot)
+            system%pe = system%pe + oblpot
          end if
    
          system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl))