anchor.diff:  
       end do rattle_loop
       call p%time_rattle_pos%tac()
  +    deallocate(i_mass)
   
       if (rattle_i==rattle_max) write(*,*) 'rattle_max reached in rattle_body_pos'
   
  
positive.code:           g = rs - sqrt(rs**2 - rsq*(ssq-d**2))
            g = g / (dt * inv_mass * rsq)
  
            p% pos(:,i1) = p% pos(:,i1) - g*dt*r/mass1
            p% pos(:,i2) = p% pos(:,i2) + g*dt*r/mass2
  
            p% vel(:,i1) = p% vel(:,i1) - g*r/mass1
            p% vel(:,i2) = p% vel(:,i2) + g*r/mass2
         end do
  
         do i_link = 1, n_link
  
negative.code:     double precision :: s(3) ! direction vector
      double precision :: r(3) ! old direction vector
      double precision :: rsq, ssq, rs
      double precision :: i_mass1, i_mass2, inv_mass
      double precision, allocatable :: i_mass(:)
  
      integer :: rattle_i, rattle_max, i_link, n_link
      integer :: i1, i2
  
      n_link = size(links, dim=2)
      rattle_max = 1000
  
      allocate(i_mass(size(p%mass)))
      i_mass = 1/p%mass
  
      call p%time_rattle_pos%tic()
      rattle_loop: do rattle_i = 1, rattle_max
  