anchor.diff:  
       allocate(i_mass(size(p%mass)))
       i_mass = 1/p%mass
  +
  +    allocate(im(3,p%Nmax))
  +    allocate(im_r(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
  +
  +    s = p%pos_rattle(:,1)
  +    do i1 = 2, p%Nmax
  +       im_r(:,i1) = floor(((p%pos_rattle(:,i1)-s)/edges) + 0.5d0)
  +       p%pos_rattle(:,i1) = p%pos_rattle(:,i1) - im_r(:,i1)*edges
  +    end do
   
       call p%time_rattle_pos%tic()
       rattle_loop: do rattle_i = 1, rattle_max
  
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:           i_mass1 = i_mass(p%species(i1))
            i_mass2 = i_mass(p%species(i2))
            inv_mass = i_mass1 + i_mass2
            s = p%pos(:,i1)-p% pos(:,i2)
            k = abs(dot_product(p%vel(:,i1)-p%vel(:,i2), s) / (d**2*inv_mass))
            if (k>error) error = k
         end do
  