def find_groups(sl1, sl2):
    """
    Given two lists of atom lists, find the groups of sets in each
    list that only contain each others' elements (i.e. if somehow we
    have two parallel reactions in one.)

    Parameters
    ----------
    sl1, sl2 : list of lists
        List of lists of atom indices corresponding to molecules.

    Returns
    -------
    list of lists
        Atom indices corresponding to separately reacting groups of molecules.
    """
    # Convert to sets
    sl1c = [set(s) for s in sl1]
    sl2c = [set(s) for s in sl2]
    # Iterate this while loop until we find a single set of atom groups
    while set([tuple(sorted(list(s))) for s in sl1c]) != set([tuple(sorted(list(s))) for s in sl2c]):
        # Double loop over molecule atom indices
        for s1 in sl1c:
            for s2 in sl2c:
                # For any pair of molecules that have any overlapping atoms,
                # add all atoms in each atom set to the other.
                if len(s1.intersection(s2)) > 0:
                    s1.update(s2)
                    s2.update(s1)
    result = sorted([list(t) for t in list(set([tuple(sorted(list(s))) for s in sl1c]))])
    return result