def consecutiveDistanceByProbability(r1,r2,p,xcontact=2):
    """
    Upper bound distance constraints for consecutive domains
    return surface to surface distance.
    
    Parameters
    -----------
    r1,r2 : int
        Radius for beads
    p : float Probability for contact
    xcontact : int
        scaling of (r1+r2) where a contact is defined. By default, 
        center to center distance D = 2*(r1+r2) is defined as contact.
    """
    if p > 0:
        d = (r1+r2)*(1. + (xcontact**3-1)/p)**(1./3.)
    else:
        d = 100*(r1+r2) # just a big number
    return d-r1-r2