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, i-1
   
  -            Ep = Ep -  G * m(j) * m(i) / sqrt(norm2(r(:, i) - r(:, j))**2 + epsilon2)
  +            Ep = Ep - G * m(j) * m(i) / sqrt(sum((r(:, i) - r(:, j))**2) + epsilon2)
   
            end do
   
            k = length + 1 - i
  -         Ec = Ec + 0.5_xp * m(k) * norm2(v(:,k))**2
  +         Ec = Ec + 0.5_xp * m(k) * sum(v(:,k)**2)
   
            do j = 1, k-1
   
  -            Ep = Ep -  G * m(j) * m(k) / sqrt(norm2(r(:, k) - r(:, j))**2 + epsilon2)
  +            Ep = Ep - G * m(j) * m(k) / sqrt(sum((r(:, k) - r(:, j))**2) + epsilon2)
   
            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:       !$OMP DO SCHEDULE(RUNTIME)
        do i = istart, iend
  
           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(sum((r(:, i) - r(:, j))**2) + epsilon2)
              end if
  
           end do
  