Skip to content

Net charge abfe -fix alchemical ion selection heuristic #1992

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
aqemia-benedict-tan wants to merge 2 commits into
base: net_charge_abfe
Choose a base branch
from
Draft

Net charge abfe -fix alchemical ion selection heuristic #1992

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
18c3756
Select only counterions by specifying resname in mda selection
aqemia-benedict-tan Jun 5, 2026
6fbbef5
Increase the alchemical ion search region to PBC-safe limit
aqemia-benedict-tan Jun 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 25 additions & 2 deletions openfe/protocols/openmm_afe/equil_binding_afe_method.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
import gufe
import MDAnalysis as mda
from MDAnalysis.lib.distances import calc_bonds
from MDAnalysis.lib.mdamath import triclinic_vectors
import numpy as np
import numpy.typing as npt
from gufe import (
Expand Down Expand Up @@ -1057,6 +1058,9 @@ def _get_alchemical_ion(
if abs(total_charge) > 1:
errmsg = "Cannot handle net charge correction on charges greater than one"
raise ValueError(errmsg)

# Declare resname to only select counterions (exclude ligand)
ion_resnames = [e.upper() for e in IONS[total_charge]]

univ = _get_mda_universe(
omm_topology,
Expand All @@ -1067,7 +1071,7 @@ def _get_alchemical_ion(
counter_ions_idxs = [
at.ix
for at in univ.atoms
if at.element in IONS[total_charge]
if at.resname in ion_resnames
]
counter_ions = univ.atoms[counter_ions_idxs]

Expand All @@ -1076,6 +1080,25 @@ def _get_alchemical_ion(
alchem_idxs = _get_idxs_from_residxs(topology=omm_topology, residxs=residxs)
alchem_atoms = univ.atoms[alchem_idxs]

# Cap the search at the periodic minimum-image limit
univ.trajectory[-1] # read the box from the final frame
box = univ.dimensions
if box is None or not np.all(np.isfinite(box)) or np.any(box[:3] <= 0.0):
errmsg = (
"Could not identify a periodic box for the co-alchemical ion "
"search. The search distance is bounded by the periodic "
"minimum-image convention, which requires a valid (non-zero) "
f"box in the equilibrated frame; got dimensions: {box}."
)
raise ValueError(errmsg)
bvecs = triclinic_vectors(box) # box vectors, Angstrom
vol = abs(np.dot(bvecs[0], np.cross(bvecs[1], bvecs[2])))
perp_widths = [
vol / np.linalg.norm(np.cross(bvecs[(i + 1) % 3], bvecs[(i + 2) % 3]))
for i in range(3)
]
max_dist = (0.5 * min(perp_widths) / 10.0 - 0.1) * offunit.nanometer # Å -> nm

# Re-using a utility from the restraints utilities
# TODO: rename this class!
atom_finder = FindHostAtoms(
Expand All @@ -1084,7 +1107,7 @@ def _get_alchemical_ion(
min_search_distance=settings['alchemical_settings'].alchemical_ion_min_distance,
# set max distance to just above solvent padding to avoid picking
# an ion more than half a box distance away
max_search_distance=settings['solvation_settings'].solvent_padding + 0.1 * offunit.nanometer,
max_search_distance=max_dist,
)

# only run on the final frame
Expand Down