anchor.diff:      double precision :: rsq, ssq, rs
       double precision :: i_mass1, i_mass2, inv_mass
       double precision, allocatable :: i_mass(:)
  +    integer, allocatable :: im(:,:), im_r(:,:)
   
       integer :: rattle_i, rattle_max, i_link, n_link
       integer :: i1, i2
  
positive.code:     end do rattle_loop
      call p%time_rattle_vel%tac()
  
      deallocate(i_mass)
  
      if (rattle_i==rattle_max) write(*,*) 'rattle_max reached in rattle_body_vel'
  
  
negative.code:     rattle_max = 1000
      allocate(i_mass(size(p%mass)))
      i_mass = 1/p%mass
  
      allocate(im(3,p%Nmax))
  
      s = p%pos(:,1)
      do i1 = 2, p%Nmax
         im(:,i1) = floor(((p%pos(:,i1)-s)/edges) + 0.5d0)
         p%pos(:,i1) = p%pos(:,i1) - im(:,i1)*edges
      end do
  
      call p%time_rattle_vel%tic()
      rattle_loop: do rattle_i = 1, rattle_max
  