anchor.diff:        !$OMP DO SCHEDULE(RUNTIME)
         do i = istart, iend
   
  -         Ec = Ec + 0.5_xp * m(i) * norm2(v(:,i))**2
  +         Ec = Ec + 0.5_xp * m(i) * sum(v(:,i)**2)
   
            do j = 1, npoints
   
               ! TODO: Check this condition when using mpi_flag=1 and distributed energy computation  
               if (i /= j) then
  -               Ep = Ep - 0.5_xp * G * m(j) * m(i) / sqrt(norm2(r(:, i) - r(:, j))**2 + epsilon2)
  +               Ep = Ep - 0.5_xp * G * m(j) * m(i) / sqrt(sum((r(:, i) - r(:, j))**2) + epsilon2)
               end if
   
            end do
  
positive.code:          do j = 1, i-1
  
              vec = r(:, i) - r(:, j)
              tmp = G / (norm2(vec)**2 + epsilon2)**1.5_xp
  
              a(:, i) = a(:, i) - tmp*m(j)*vec
              a(:, j) = a(:, j) + tmp*m(i)*vec
  
           end do
  
  
negative.code:          do j = 1, k-1
  
              vec = r(:, k) - r(:, j)
              tmp = vec * G / sqrt(sum(vec**2) + epsilon2)**3
  
              a(:, k) = a(:, k) - m(j)*tmp
              a(:, j) = a(:, j) + m(k)*tmp
  
           end do
  
  