def recurse_atoms(atom, _visited_indices=set()):
    """
    Recursively visit all atoms within a molecule for the purposes
    of determining the molecules (all sets of atoms connected by bonds)
    in the system.
    """
    
    _visited_indices.add(atom.idx)
    for bonded_atom in atom.bond_partners:
        if not bonded_atom.idx in _visited_indices:
            branch_indices = recurse_atoms(bonded_atom, _visited_indices)
            _visited_indices.update(branch_indices)
    return _visited_indices