From 36662df100f2e59ac959c96b5c881700189ab5b8 Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Fri, 18 Apr 2025 11:05:23 -0600 Subject: [PATCH] MAINT: rounding residue partial charges * Fixes gh-5032. * Given that 1e-4 is already optimistic for electron charge unit precision in MD forcefields, this patch proposes to round residue net charges to 4 decimal places to avoid the confusion described in the matching ticket. * Some obvious questions might include: why are we only doing this for AMBER topologies? why only round at the residue level and not i.e., upstream at atom level and also at "segment" levels? This is mostly for scoping, because this initial patch is targeted at the specific issue raised, and because it probably doesn't make sense to proceed much farther until things are discussed. --- package/MDAnalysis/core/topologyattrs.py | 8 ++++++ package/MDAnalysis/topology/TOPParser.py | 7 ++++- .../MDAnalysisTests/topology/test_top.py | 28 ++++++++++++++++++- 3 files changed, 41 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 359bd20c9c..1285ebb2ac 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -2146,6 +2146,11 @@ class Charges(AtomAttr): transplants = defaultdict(list) dtype = float + + def __init__(self, values, guessed=False, round_val=None): + super().__init__(values, guessed) + self.round = round_val + @staticmethod def _gen_initial_values(na, nr, ns): return np.zeros(na) @@ -2160,6 +2165,9 @@ def get_residues(self, rg): for i, row in enumerate(resatoms): charges[i] = self.values[row].sum() + if self.round is not None: + charges = np.round(charges, self.round) + return charges def get_segments(self, sg): diff --git a/package/MDAnalysis/topology/TOPParser.py b/package/MDAnalysis/topology/TOPParser.py index ba61bea0b0..603c8c70b5 100644 --- a/package/MDAnalysis/topology/TOPParser.py +++ b/package/MDAnalysis/topology/TOPParser.py @@ -453,7 +453,12 @@ def parse_charges(self, num_per_record, numlines): vals = self.parsesection_mapper(numlines, lambda x: float(x)) charges = np.array(vals, dtype=np.float32) charges /= 18.2223 # to electron charge units - attr = Charges(charges) + # in practice, the precision of partial charges + # in MD simulation forcefields is typically limited + # to a few decimal places (i.e., 1e-4 is already optimistic) + # see gh-5032 + # here, we round the residue net charges to 4 decimal places + attr = Charges(charges, round_val=4) return attr def parse_masses(self, num_per_record, numlines): diff --git a/testsuite/MDAnalysisTests/topology/test_top.py b/testsuite/MDAnalysisTests/topology/test_top.py index a67c9283b7..5e2a92dc41 100644 --- a/testsuite/MDAnalysisTests/topology/test_top.py +++ b/testsuite/MDAnalysisTests/topology/test_top.py @@ -27,7 +27,7 @@ import MDAnalysis as mda import numpy as np import pytest -from numpy.testing import assert_equal +from numpy.testing import assert_equal, assert_allclose from MDAnalysisTests.datafiles import PRM # ache.prmtop from MDAnalysisTests.datafiles import PRM7 # tz2.truncoct.parm7.bz2 @@ -837,3 +837,29 @@ def test_warning(self, parm, errmsgs): for msg in errmsgs: assert any(msg in recmsg for recmsg in messages) + + +def test_gh_5032(): + u = mda.Universe(PRM) + # check individual residue.charge + actual = [] + for res in u.residues: + actual.append(res.charge) + expected = [1.0, # N-terminal A + -1.0, # E + 0.0, # F + 0.0, # H (HIE) + 1.0, # R + 0.0, # W + 0.0, # S + 0.0, # S + 0.0, # Y + 0.0, # M + 0.0, # V + 0.0, # H (HIE) + 0.0, # W + 0.0] # C-terminal K + assert_allclose(actual, expected) + # check residues.charges, which are set + # separately + assert_allclose(u.residues.charges, expected)