Test if two sites are bonded, up to a certain limit.
        Args:
            site1 (Site): First site
            site2 (Site): Second site
            tol (float): Relative tolerance to test. Basically, the code
                checks if the distance between the sites is less than (1 +
                tol) * typical bond distances. Defaults to 0.2, i.e.,
                20% longer.
            bond_order: Bond order to test. If None, the code simply checks
                against all possible bond data. Defaults to None.
            default_bl: If a particular type of bond does not exist, use this
                bond length. If None, a ValueError will be thrown.
        Returns:
            Boolean indicating whether two sites are bonded.

def is_bonded(site1, site2, tol=0.2, bond_order=None, default_bl=None):
        """
        Test if two sites are bonded, up to a certain limit.
        Args:
            site1 (Site): First site
            site2 (Site): Second site
            tol (float): Relative tolerance to test. Basically, the code
                checks if the distance between the sites is less than (1 +
                tol) * typical bond distances. Defaults to 0.2, i.e.,
                20% longer.
            bond_order: Bond order to test. If None, the code simply checks
                against all possible bond data. Defaults to None.
            default_bl: If a particular type of bond does not exist, use this
                bond length. If None, a ValueError will be thrown.
        Returns:
            Boolean indicating whether two sites are bonded.
        """
        sp1 = list(site1.species.keys())[0]
        sp2 = list(site2.species.keys())[0]
        dist = site1.distance(site2)
        syms = tuple(sorted([sp1.symbol, sp2.symbol]))
        if syms in bond_lengths:
            all_lengths = bond_lengths[syms]
            if bond_order:
                return dist < (1 + tol) * all_lengths[bond_order]
            for v in all_lengths.values():
                if dist < (1 + tol) * v:
                    return True
            return False
        elif default_bl:
            return dist < (1 + tol) * default_bl
        raise ValueError("No bond data for elements {} - {}".format(*syms))