Description

I am encountering a systematic failure during ABFE setup when working with a GROMACS system solvated with KCl, using a ligand with net charge -1. The issue is fully reproducible, and the setup hangs at the find_alchemical_counterions stage without any error. The log shows:

ligand_atoms net charge: -1
Ions net charges: [('CL', -1), ...]

No cations are detected, and the process does not progress.

Environment

• YANK version: 0.25.2
• Python: 3.9.19
• CUDA: available
• Execution: via SLURM
• System: GROMACS .top / .gro with KCl solvent

What fails

• The pipeline hangs during counterion assignment
• The failure happens before any MCMC or sampling steps start
• The log stops immediately after reading the system and computing ligand_atoms net charge

Relevant log excerpt:

2026-01-29 19:32:42,879: Single node: executing <function find_alchemical_counterions at 0x2b659d3e0dc0>
2026-01-29 19:32:42,883: ligand_atoms net charge: -1
2026-01-29 19:32:43,343: Ions net charges: [('CL', -1), ('CL', -1), ...]

Root Cause

In yank/yank.py, the Topography class defines the ions_atoms property as:

ION_RESIDUE_NAMES = {'NA', 'CL'}

It checks for ions via:

• residue.name in ION_RESIDUE_NAMES, or
• Presence of + / - in residue.name

GROMACS K⁺ ions are labeled 'K' (no +) and are not in ION_RESIDUE_NAMES. Therefore:

• K⁺ ions are not detected
• Assignment of alchemical counterions fails
• Initialization hangs without raising an exception

Workaround / Temporary Fix

Manually modify the code:

ION_RESIDUE_NAMES = {'NA', 'CL', 'K'}

After this change, K⁺ ions are detected, counterions are assigned correctly, and ABFE setup completes.

Expected Behavior

• All cations/anions defined in the YAML (positive_ion, negative_ion) should be automatically detected regardless of residue name
• Initialization should either proceed or raise a clear exception if ions cannot be detected

Suggested PR / Fix

A robust solution could involve:

1. Automatically including all ions specified in the YAML (positive_ion / negative_ion) in ions_atoms
2. Or detecting ions based on formal charge instead of hardcoded residue names

This would remove the need for hardcoding residue names like 'K', 'NA', 'CL' and make the code compatible with a wider variety of GROMACS systems.

Proposed Code Change

Option 1: Extend ION_RESIDUE_NAMES dynamically

In yank/yank.py, modify the Topography class to include ions from YAML configuration:

@property
def ions_atoms(self):
    """Indices of ions in the system, including user-specified ions."""
    
    # Base set of common ions
    ION_RESIDUE_NAMES = {'NA', 'CL', 'K', 'CA', 'MG', 'ZN'}
    
    # Optionally extend with ions from configuration
    # (requires access to YAML config, may need refactoring)
    
    ions_atoms = set()
    for residue in self._topology.residues():
        residue_name = residue.name.upper()
        if (residue_name in ION_RESIDUE_NAMES or 
            '+' in residue.name or 
            '-' in residue.name):
            ions_atoms.update(atom.index for atom in residue.atoms())
    
    return ions_atoms

Option 2: Detect ions by formal charge

@property
def ions_atoms(self):
    """Indices of ions in the system, detected by single-atom residues with charge."""
    
    ions_atoms = set()
    for residue in self._topology.residues():
        # Check if single-atom residue (typical for ions)
        atoms_in_residue = list(residue.atoms())
        if len(atoms_in_residue) == 1:
            atom = atoms_in_residue[0]
            # Check if charged (requires access to partial charges)
            # This may require extracting charge from force field parameters
            ions_atoms.add(atom.index)
    
    return ions_atoms

Additional Context

This issue affects any GROMACS user working with:

• K⁺ / Cl⁻ systems (biochemistry standard)
• Other ions like Ca²⁺, Mg²⁺, Zn²⁺
• Custom ion residue names

A more generic solution would benefit the broader YANK community.