anchor.diff:  
            end do
   
  +         k = npoints + 1 - i
  +
  +         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
  +
  +         end do
  +
         end do
         !$OMP END DO
         !$OMP END PARALLEL
  
positive.code:       real(kind = xp), intent(in)  :: r(:,:)
        real(kind = xp), intent(out) :: v(:,:)
  
        real(kind = xp) :: omega
  
        omega = 0.5_xp / pi
        v(1, :) =   omega * r(2, :)
        v(2, :) = - omega * r(1, :)
        v(3, :) = 0._xp
  
     end subroutine initial_speeds
  
negative.code:       real(kind=xp), intent(out) :: a(:,:)
  
        real(kind=xp), dimension(3) :: vec, tmp
        integer :: i, j, k
  
        a = 0._xp
  
        !$OMP PARALLEL PRIVATE(j, k, vec, tmp) REDUCTION(+:a)
        !$OMP DO SCHEDULE(RUNTIME)
        do i = 1, npoints/2
  
           do j = 1, i-1
  
  