anchor.diff:            i_mass1 = i_mass(p%species(i1))
             i_mass2 = i_mass(p%species(i2))
             inv_mass = i_mass1 + i_mass2
  -          s = rel_pos(p% pos(:,i1), p% pos(:,i2), edges)
  +          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
  
positive.code:           i1 = links(1,i_link)
            i2 = links(2,i_link)
            d = distances(i_link)
            r = rel_pos(p% pos(:,i1), p%pos(:,i2), edges)
            g = abs(sqrt(dot_product(r,r)) - d)
            if (g > error) error = g
         end do
  
negative.code:           i_mass2 = i_mass(p%species(i2))
            inv_mass = i_mass1 + i_mass2
  
            s = p%pos(:,i1)-p%pos(:,i2)
  
            k = dot_product(p%vel(:,i1)-p%vel(:,i2), s) / (d**2*inv_mass)
  
  