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

Open
tylerjereddy wants to merge 1 commit into
base: develop
Choose a base branch
from
+15 −6
Open

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
2 changes: 2 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,10 @@ def _read_first_frame(self):
self.ts._pos = np.asarray(
tpr_utils.ndo_rvec(data, th.natoms), dtype=np.float32
)
self.convert_pos_from_native(self.ts._pos)
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.

I think you need to wrap it around self.convert_units like we do here https://github.com/MDAnalysis/mdanalysis/blob/develop/package%2FMDAnalysis%2Fcoordinates%2FGRO.py#L255

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.

I noticed that. However, its use does need motivating with a regression test that fails without it, which wasn't immediately obvious to me since we don't write TPR files.

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.

Also, GROWriter __init__ has:

        convert_units : bool (optional)
            units are converted to the MDAnalysis base format; [``True``]

but farther down it has:

        if self.convert_units:
            # Convert back to nm from Angstroms,
            # Not inplace because AtomGroup is not a copy
            positions = self.convert_pos_to_native(positions, inplace=False)

so while the default behavior makes sense, the docs are garbled--basically wrong, right? Should it even be "allowed" to write with A units for GRO?

Copy link
Copy Markdown
Member

@IAlibay IAlibay May 6, 2026

Choose a reason for hiding this comment

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

I'm so sorry, I've been staring at this for a while and I just can't seem to think this one through. What do you recommend doing here? I am under the impression that the self.convert_units bits are part of the base reader class API - is this correct?

Copy link
Copy Markdown
Member Author

@tylerjereddy tylerjereddy May 7, 2026

Choose a reason for hiding this comment

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

I think my point was just that convert_units is a bit weird in that it can allow you to do things that are not file format compliant, and that our docs are sometimes wrong on the matter -- I've opened gh-5378 on that topic. Personally, I'm not a huge fan of some of the things it allows, but I suppose that's a separate discussion.

For this PR, I believe you have a point that I should conform to the option of reading in either A or nm units, and add a test for that for TPRReader, probably parametrized over the two options.

if th.bV:
self.ts.velocities = np.asarray(
tpr_utils.ndo_rvec(data, th.natoms), dtype=np.float32
)
self.convert_velocities_from_native(self.ts.velocities)
self.ts.has_velocities = True
19 changes: 13 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 @@ -80,6 +80,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 +463,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 +490,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 Down
Loading