diff --git a/package/CHANGELOG b/package/CHANGELOG index 7fa27bb5178..3d2a64b41ee 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -133,6 +133,8 @@ Enhancements explicit in the topology (Issue #2468, PR #2775) Changes + * TPRParser now loads TPR files with `tpr_resid_from_one=True` by default, + which starts TPR resid indexing from 1 (instead of 0 as in 1.x) (Issue #2364, PR #3152) * Introduces encore specific C compiler arguments to allow for lowering of optimisations on non-x86 platforms (Issue #1389, PR #3149) * Continuous integration uses mamba rather than conda to install the diff --git a/package/MDAnalysis/topology/TPRParser.py b/package/MDAnalysis/topology/TPRParser.py index 786f27c94d6..6ff6f185847 100644 --- a/package/MDAnalysis/topology/TPRParser.py +++ b/package/MDAnalysis/topology/TPRParser.py @@ -176,12 +176,27 @@ class TPRParser(TopologyReaderBase): """ format = 'TPR' - def parse(self, **kwargs): + def parse(self, tpr_resid_from_one=True, **kwargs): """Parse a Gromacs TPR file into a MDAnalysis internal topology structure. + Parameters + ---------- + tpr_resid_from_one: bool (optional) + Toggle whether to index resids from 1 or 0 from TPR files. + TPR files index resids from 0 by default, even though GRO and ITP + files index from 1. + Returns ------- structure : dict + + + .. versionchanged:: 1.1.0 + Added the ``tpr_resid_from_one`` keyword to control if + resids are indexed from 0 or 1. Default ``False``. + + .. versionchanged:: 2.0.0 + Changed to ``tpr_resid_from_one=True`` by default. """ with openany(self.filename, mode='rb') as infile: tprf = infile.read() @@ -219,7 +234,8 @@ def parse(self, **kwargs): tpr_utils.ndo_real(data, state_ngtc) # relevant to Berendsen tcoupl_lambda if th.bTop: - tpr_top = tpr_utils.do_mtop(data, th.fver) + tpr_top = tpr_utils.do_mtop(data, th.fver, + tpr_resid_from_one=tpr_resid_from_one) else: msg = f"{self.filename}: No topology found in tpr file" logger.critical(msg) diff --git a/package/MDAnalysis/topology/tpr/utils.py b/package/MDAnalysis/topology/tpr/utils.py index 5aab588ea1e..b764e7ec59e 100644 --- a/package/MDAnalysis/topology/tpr/utils.py +++ b/package/MDAnalysis/topology/tpr/utils.py @@ -46,6 +46,7 @@ The module also contains the :func:`do_inputrec` to read the TPR header with. """ + import numpy as np import xdrlib import struct @@ -284,7 +285,7 @@ def extract_box_info(data, fver): return obj.Box(box, box_rel, box_v) -def do_mtop(data, fver): +def do_mtop(data, fver, tpr_resid_from_one=False): # mtop: the topology of the whole system symtab = do_symtab(data) do_symstr(data, symtab) # system_name @@ -371,6 +372,9 @@ def do_mtop(data, fver): molnums = np.array(molnums, dtype=np.int32) segids = np.array(segids, dtype=object) resids = np.array(resids, dtype=np.int32) + if tpr_resid_from_one: + resids += 1 + resnames = np.array(resnames, dtype=object) (residx, new_resids, (new_resnames, diff --git a/testsuite/MDAnalysisTests/analysis/test_density.py b/testsuite/MDAnalysisTests/analysis/test_density.py index 33e2c7c9e73..937f25f03c4 100644 --- a/testsuite/MDAnalysisTests/analysis/test_density.py +++ b/testsuite/MDAnalysisTests/analysis/test_density.py @@ -148,7 +148,7 @@ class DensityParameters(object): @pytest.fixture() def universe(self): - return mda.Universe(self.topology, self.trajectory) + return mda.Universe(self.topology, self.trajectory, tpr_resid_from_one=False) class TestDensityAnalysis(DensityParameters): def check_DensityAnalysis(self, ag, ref_meandensity, diff --git a/testsuite/MDAnalysisTests/core/test_atomselections.py b/testsuite/MDAnalysisTests/core/test_atomselections.py index 831428da852..6d52ccf0a8f 100644 --- a/testsuite/MDAnalysisTests/core/test_atomselections.py +++ b/testsuite/MDAnalysisTests/core/test_atomselections.py @@ -485,7 +485,7 @@ class TestSelectionsTPR(object): @staticmethod @pytest.fixture(scope='class') def universe(): - return MDAnalysis.Universe(TPR,XTC) + return MDAnalysis.Universe(TPR, XTC, tpr_resid_from_one=False) @pytest.mark.parametrize('selstr', [ 'same fragment as bynum 1', @@ -1299,7 +1299,7 @@ def test_mass_sel_warning(u_fake_masses): @pytest.mark.parametrize("selstr,n_res", [ - ("resnum -10 to 3", 14), + ("resnum -10 to 3", 13), ("resnum -5--3", 3), # select -5 to -3 ("resnum -3 : -5", 0), # wrong way around ]) diff --git a/testsuite/MDAnalysisTests/topology/test_tprparser.py b/testsuite/MDAnalysisTests/topology/test_tprparser.py index 7826a4d6613..f6435ccf8a4 100644 --- a/testsuite/MDAnalysisTests/topology/test_tprparser.py +++ b/testsuite/MDAnalysisTests/topology/test_tprparser.py @@ -39,6 +39,7 @@ ) from MDAnalysisTests.topology.base import ParserBase import MDAnalysis.topology.TPRParser +import MDAnalysis as mda BONDED_TPRS = ( TPR510_bonded, @@ -287,3 +288,14 @@ def test_elements(): 'H', '', 'Na', 'Na', 'Na', 'Na', ], dtype=object) assert_equal(topology.elements.values[-20:], reference) + + +@pytest.mark.parametrize("resid_from_one,resid_addition", [ + (False, 0), + (True, 1), # status quo for 2.x + ]) +def test_resids(resid_from_one, resid_addition): + u = mda.Universe(TPR, tpr_resid_from_one=resid_from_one) + resids = np.arange(len(u.residues)) + resid_addition + assert_equal(u.residues.resids, resids, + err_msg="tpr_resid_from_one kwarg not switching resids")