anchor.diff:      mass2 = p%mass(p%species(2))
       inv_mass = 1/mass1 + 1/mass2
   
  -    s = p% pos(:,1) - p% pos(:,2)
  +    s = rel_pos(p% pos(:,1), p% pos(:,2), edges)
   
       k = dot_product(p%vel(:,1)-p%vel(:,2), s) / (d**2*inv_mass)
   
  
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
  