anchor.diff:           do j = 1, k-1
   
               vec = r(:, k) - r(:, j)
  -            tmp = G / (norm2(vec)**2 + epsilon2)**1.5_xp
  -
  -            a(:, k) = a(:, k) - tmp*m(j)*vec
  -            a(:, j) = a(:, j) + tmp*m(k)*vec
  +            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
   
  
positive.code: 
           do j = jstart, jend
  
              ! When i = j, this is zero if r1=r2, but important otherwise
              vec = r1(:, i) - r2(:, j)
              tmp = G / (norm2(vec)**2 + epsilon2)**1.5_xp
  
              a(:, i) = a(:, i) - tmp*m(j)*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
  