anchor.diff:              vec = r1(:, i) - r2(:, j)
               tmp = vec * G / sqrt(sum(vec**2) + epsilon2)**3
   
  -            a1(:, i) = a1(:, i) - m(j)*tmp
  -            a2(:, j) = a2(:, j) + m(i)*tmp
  +            a1(:, i) = a1(:, i) - m2(j)*tmp
  +            a2(:, j) = a2(:, j) + m1(i)*tmp
   
            end do
   
  
positive.code:             allocate(r_i(3, N))
              allocate(r_np_i(3, N))
              allocate(r_right(3, N))
  
              a_right = 0._xp
  
  
negative.code: 
     end subroutine compute_force
  
     subroutine compute_force_mpi (m1, r1, istart, iend, m2, r2, jstart, jend, a1, a2)
        implicit none
  
        real(xp), intent(in) :: m1(:), m2(:)
        real(xp), intent(in) :: r1(:, :), r2(:, :)
        integer,  intent(in) :: istart, iend, jstart, jend
  
  