anchor.diff:      double precision :: rsq, ssq, rs
       double precision :: mass1, mass2, inv_mass
   
  -    r = p% pos_rattle(:,1) - p% pos_rattle(:,2)
  -    s = p% pos(:,1) - p% pos(:,2)
  +    r = rel_pos(p% pos_rattle(:,1),p% pos_rattle(:,2), edges)
  +    s = rel_pos(p% pos(:,1),p% pos(:,2), edges)
       mass1 = p%mass(p%species(1))
       mass2 = p%mass(p%species(2))
       inv_mass = 1/mass1 + 1/mass2
  
positive.code: 
    end subroutine md_vel
  
    subroutine rattle_dimer_pos(p, d, dt)
      type(particle_system_t), intent(inout) :: p
      double precision, intent(in) :: d
      double precision, intent(in) :: dt
  
      double precision :: g
      double precision :: s(3) ! direction vector
  
negative.code: module md
    use particle_system
    use common
    implicit none
  
    private
  