Recursive algorithm to permute through all possible combinations of the
    initially supplied symmetry operations to arrive at a complete set of
    operations mapping a single atom to all other equivalent atoms in the
    point group.  This assumes that the initial number already uniquely
    identifies all operations.

    Args:
        symmops ([SymmOp]): Initial set of symmetry operations.

    Returns:
        Full set of symmetry operations.

def generate_full_symmops(symmops, tol):
    """
    Recursive algorithm to permute through all possible combinations of the
    initially supplied symmetry operations to arrive at a complete set of
    operations mapping a single atom to all other equivalent atoms in the
    point group.  This assumes that the initial number already uniquely
    identifies all operations.

    Args:
        symmops ([SymmOp]): Initial set of symmetry operations.

    Returns:
        Full set of symmetry operations.
    """
    # Uses an algorithm described in:
    # Gregory Butler. Fundamental Algorithms for Permutation Groups.
    # Lecture Notes in Computer Science (Book 559). Springer, 1991. page 15
    UNIT = np.eye(4)
    generators = [op.affine_matrix for op in symmops
                  if not np.allclose(op.affine_matrix, UNIT)]
    if not generators:
        # C1 symmetry breaks assumptions in the algorithm afterwards
        return symmops
    else:
        full = list(generators)

        for g in full:
            for s in generators:
                op = np.dot(g, s)
                d = np.abs(full - op) < tol
                if not np.any(np.all(np.all(d, axis=2), axis=1)):
                    full.append(op)

        d = np.abs(full - UNIT) < tol
        if not np.any(np.all(np.all(d, axis=2), axis=1)):
            full.append(UNIT)
        return [SymmOp(op) for op in full]