[ WIP ] Support for charge changes in the HybTop protocol for explicitly solvated systems#1978
[ WIP ] Support for charge changes in the HybTop protocol for explicitly solvated systems#1978hannahbaumann wants to merge 8 commits into
Conversation
…for explicitly solvated systems
| """ | ||
| Infer positive and negative ion residue names from a topology by | ||
| finding the most common monovalent ion of each charge type. | ||
| Falls back to NA/CL if none are found. |
There was a problem hiding this comment.
I'm doing this mostly because e.g. maestro doesn't always add ions (in the a2a example there were no negative ions in the system).
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
| system_mapping: dict, | ||
| charge_difference: int, | ||
| solvent_component: SolventComponent, | ||
| positive_ion_resname: str, |
There was a problem hiding this comment.
use element or formal charge instead?
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1978 +/- ##
==========================================
- Coverage 94.94% 90.54% -4.41%
==========================================
Files 216 216
Lines 20481 20540 +59
==========================================
- Hits 19445 18597 -848
- Misses 1036 1943 +907
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
Replace the monolithic `_get_ion_and_water_parameters` with three focused functions that handle ion parameter lookup in order of preference: 1. Most abundant monovalent ion of the correct charge sign already present in the topology. 2. Divalent or wrong-sign ions are skipped by the monovalent filter. 3. If no suitable ion exists (e.g. zero ionic strength or a SolvatedPDBComponent without ions), a minimal single-ion system is built from the forcefield to obtain NA/CL parameters. This approach is adapted from OpenFE PR OpenFreeEnergy#1978. The key API change is that `handle_alchemical_waters` and `_handle_net_charge` now accept a `forcefield` (from the system generator) instead of a `solvent_component`, since what we actually need is the ability to parameterize an ion — not the component metadata. The mutation only modifies NonbondedForce parameters (charge, sigma, epsilon), leaving bonded terms untouched. This is safe for rigid water models (TIP3P, SPC/E) where geometry is enforced by SETTLE/SHAKE constraints rather than harmonic forces. Added documentation warning that flexible water models would produce unphysical forces on the ghost hydrogens. Also adds SolvatedPDBComponent support to `_validate_charge_difference` so membrane systems can use explicit charge correction.
| "between the end states. This will be addressed by " | ||
| f"transforming a water into a {ion} ion" | ||
| ) | ||
| # resolve ion names from SolventComponent or topology |
There was a problem hiding this comment.
TBD: It's probably ok squash this into a single check that just says we'll convert to an ion of opposite charge.
| best_resname = ion_counts.most_common(1)[0][0] | ||
| return nbf.getParticleParameters(ion_atom_indices[best_resname]) | ||
|
|
||
| # No matching ion in topology: fall back to NA/CL from forcefield |
There was a problem hiding this comment.
TBD: how does this end up happening? I.e. it is purely a case where your user didn't add any counterions or did a strict neutralization?
There was a problem hiding this comment.
Yes, in the a2a case this was what the maestro prepped membrane from Ross et al. looked like, no counter ions, just a structural Na in the protein, but no ions in the solvent.
| continue | ||
| charge, _, _ = nbf.getParticleParameters(atom.index) | ||
| charge_val = charge.value_in_unit(omm_unit.elementary_charge) | ||
| if (abs(abs(charge_val) - 1.0) < 0.01 and np.sign(charge_val) == desired_sign): |
There was a problem hiding this comment.
Maybe np.isclose would be simpler.
…nergy/openfe into membrane_charge_change
|
No API break detected ✅ |
| ion_charge, ion_sigma, ion_epsilon = _get_ion_parameters( | ||
| topology, system, charge_difference, forcefield, | ||
| ) | ||
| o_charge, h_charge = _get_water_parameters(topology, system, water_resname) |
There was a problem hiding this comment.
Could we just use the single ion system trick to get the parameters for the water now that we have the forcefield? This would avoid the need to know the water resname and looping over the full system to find a water molecule?
| @@ -143,14 +249,10 @@ def handle_alchemical_waters( | |||
| A dictionary of system mappings between the stateA and stateB systems | |||
| charge_difference : int | |||
| The charge difference between state A and state B. | |||
There was a problem hiding this comment.
I would make it clear in all of the doc strings for charge_difference how this should be calculated, stateA - stateB
3 participants
This is just a first idea (and the
_get_ion_resnames_from_topologywas written together with Claude). This is mostly just to show the different things we have to consider, but I don't like the dictionary idea of known ions yet. Maybe it would also be ok to always use Cl and Na ions for explicitly solvated systems, no matter if other ion types were used as counterions?Checklist
newsentry, or the changes are not user-facing.pre-commit.ci autofix.Manual Tests: these are slow so don't need to be run every commit, only before merging and when relevant changes are made (generally at reviewer-discretion).
Developers certificate of origin