anchor.diff:            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
  +          p% pos(:,i1) = p% pos(:,i1) - g*dt*r*i_mass1
  +          p% pos(:,i2) = p% pos(:,i2) + g*dt*r*i_mass2
  +
  +          p% vel(:,i1) = p% vel(:,i1) - g*r*i_mass1
  +          p% vel(:,i2) = p% vel(:,i2) + g*r*i_mass2
          end do
   
          do i_link = 1, n_link
  
positive.code: 
      end do rattle_loop
      call p%time_rattle_pos%tac()
  
      if (rattle_i==rattle_max) write(*,*) 'rattle_max reached in rattle_body_pos'
  
  
negative.code: 
            r = rel_pos(p% pos_rattle(:,i1),p% pos_rattle(:,i2), edges)
            s = rel_pos(p% pos(:,i1),p% pos(:,i2), edges)
            i_mass1 = i_mass(p%species(i1))
            i_mass2 = i_mass(p%species(i2))
            inv_mass = i_mass1 + i_mass2
  
            rsq = dot_product(r,r)
            ssq = dot_product(s,s)
  