Skip to content

BUG: use base units for TPR x/v parsing #5374

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

Merged
IAlibay merged 2 commits into from
May 16, 2026
+44 −6
Merged

BUG: use base units for TPR x/v parsing #5374

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c620345
BUG: use base units for TPR x/v parsing
tylerjereddy May 3, 2026
902d22f
BUG: PR 5374 revisions
tylerjereddy May 14, 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
4 changes: 4 additions & 0 deletions package/MDAnalysis/coordinates/TPR.py
View file
Open in desktop
Copy link
Copy Markdown
Member

@IAlibay IAlibay May 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we know if dimensions setting also needs fixing?

Copy link
Copy Markdown
Member Author

@tylerjereddy tylerjereddy May 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On the latest version of this feature branch,

import MDAnalysis as mda
from MDAnalysisTests.datafiles import TPR, XTC

u = mda.Universe(TPR, XTC)
print(u.dimensions)
u = mda.Universe(TPR)
print(u.dimensions)
u = mda.Universe(TPR, TPR)
print(u.dimensions)

produces:

[80.017006 80.017006 80.017006 60.       60.       90.      ]
None
None

That looks like more like a feature request than the same category of bug, so scope might be argued. I'd probably suggest a separate issue for that one.

Copy link
Copy Markdown
Member

@IAlibay IAlibay May 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That sounds good - I was mostly just looking at other single file readers and checking what they were doing. If we don't have to, that works.

Original file line number Diff line number Diff line change
Expand Up @@ -127,8 +127,12 @@ def _read_first_frame(self):
self.ts._pos = np.asarray(
tpr_utils.ndo_rvec(data, th.natoms), dtype=np.float32
)
if self.convert_units:
self.convert_pos_from_native(self.ts._pos)
if th.bV:
self.ts.velocities = np.asarray(
tpr_utils.ndo_rvec(data, th.natoms), dtype=np.float32
)
if self.convert_units:
self.convert_velocities_from_native(self.ts.velocities)
self.ts.has_velocities = True
46 changes: 40 additions & 6 deletions testsuite/MDAnalysisTests/coordinates/test_tpr.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@
TPR_gh_5145,
)
import MDAnalysis as mda
from MDAnalysis.coordinates.TPR import TPRReader


import pytest
Expand All @@ -80,6 +81,9 @@
@pytest.mark.parametrize(
"tpr_file, exp_first_atom, exp_last_atom, exp_shape, exp_vel_first_atom, exp_vel_last_atom",
[
# NOTE: expected values are expressed in nm
# units and converted to Angstrom in the body
# of the test below, before assertions
# this case is an alanine dipeptide
# with neural network potential active
# and nonzero velocities
Expand Down Expand Up @@ -460,11 +464,15 @@ def test_basic_read_tpr(
u = mda.Universe(tpr_file, tpr_file)
else:
u = mda.Universe(tpr_file)
assert_allclose(u.atoms.positions[0, ...], exp_first_atom)
assert_allclose(u.atoms.positions[-1, ...], exp_last_atom)
assert_allclose(u.atoms.positions[0, ...], np.asarray(exp_first_atom) * 10)
assert_allclose(u.atoms.positions[-1, ...], np.asarray(exp_last_atom) * 10)
assert_equal(u.atoms.positions.shape, exp_shape)
assert_allclose(u.atoms.velocities[0, ...], exp_vel_first_atom)
assert_allclose(u.atoms.velocities[-1, ...], exp_vel_last_atom)
assert_allclose(
u.atoms.velocities[0, ...], np.asarray(exp_vel_first_atom) * 10
)
assert_allclose(
u.atoms.velocities[-1, ...], np.asarray(exp_vel_last_atom) * 10
)
assert_equal(u.atoms.velocities.shape, exp_shape)


Expand All @@ -483,8 +491,8 @@ def test_different_versions():
# reading topology and positions/velocities
# for the same system with different TPR
# file versions
exp_first_atom = [3.25000e-01, 1.00400e00, 1.03800e00]
exp_last_atom = [-2.56000e-01, 1.37300e00, 3.59800e00]
exp_first_atom = [3.25000, 10.0400, 10.3800]
exp_last_atom = [-2.56000, 13.7300, 35.9800]
exp_shape = (2263, 3)
exp_vel = np.zeros(3)
u = mda.Universe(TPR2020, TPR2024_4)
Expand All @@ -494,3 +502,29 @@ def test_different_versions():
assert_allclose(u.atoms.velocities[-1, ...], exp_vel)
assert_equal(u.atoms.positions.shape, exp_shape)
assert_equal(u.atoms.velocities.shape, exp_shape)


@pytest.mark.parametrize("convert_units", [True, False])
def test_unit_swapping(convert_units):
reader = TPRReader(TPR_xvf_2024_4, convert_units=convert_units)
ts = reader.ts
actual_positions = ts.positions
actual_velocities = ts.velocities
# nm units:
expected_first_pos = np.asarray([3.19900e00, 1.62970e00, 1.54480e00])
expected_last_pos = np.asarray([3.39350e00, 3.49420e00, 3.02400e00])
expected_first_vel = np.asarray([-0.20668714, 0.26678202, -0.10564042])
expected_last_vel = np.asarray(
[-3.38010e-02, -3.22064e-01, -1.9863836e-01]
)
if convert_units:
# convert_units=True is the MDA default for readers, so we
# expect the nm to get converted to A
factor = 10
else:
# GROMACS "native" nm units are read in
factor = 1
assert_allclose(actual_positions[0, ...], expected_first_pos * factor)
assert_allclose(actual_positions[-1, ...], expected_last_pos * factor)
assert_allclose(actual_velocities[0, ...], expected_first_vel * factor)
assert_allclose(actual_velocities[-1, ...], expected_last_vel * factor)
Loading