Skip to content

Commit cec335f

Browse files
authored
ENH: support GMX linear angle potential (#5372)
* ENH: support GMX linear angle potential * Fixes gh-5361. * Add support for parsing linear angle potentials in GMX TPR files. As described in the GMX docs linked in the matching issue, these are useful for planar molecules like CO2. * TPR position and angle data-related regression tests have been added. I did check that the angle data test fails before the patch applied here. * I believe there may be a serious unit/scaling issue in the TPR coordinate parser, but I'll open a separate ticket about that since it is not directly related. * Update CHANGELOG
1 parent 941d679 commit cec335f

6 files changed

Lines changed: 27 additions & 2 deletions

File tree

package/CHANGELOG

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,8 @@ Fixes
6060
DSSP by porting upstream PyDSSP 0.9.1 fix (Issue #4913)
6161

6262
Enhancements
63+
* Read Linear Angle bonded interaction entries from GROMACS TPR files and
64+
add the angles to the topology (Issue #5361, PR #5372)
6365
* Enables parallelization for analysis.atomicdistances.AtomicDistances
6466
(Issue #4662, PR #4822)
6567
* DSSP uses ``capped_distance`` and ``calc_bonds`` for faster distance

package/MDAnalysis/topology/tpr/utils.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -826,7 +826,7 @@ def do_moltype(data, symtab, fver):
826826
bonds += list(ik_obj.process(ias))
827827
elif ik_obj.name in ['ANGLES', 'G96ANGLES', 'CROSS_BOND_BOND',
828828
'CROSS_BOND_ANGLE', 'UREY_BRADLEY', 'QANGLES',
829-
'RESTRANGLES', 'TABANGLES']:
829+
'RESTRANGLES', 'TABANGLES', 'LINEAR_ANGLES']:
830830
angles += list(ik_obj.process(ias))
831831
elif ik_obj.name in ['PDIHS', 'RBDIHS', 'RESTRDIHS', 'CBTDIHS',
832832
'FOURDIHS', 'TABDIHS']:

testsuite/MDAnalysisTests/coordinates/test_tpr.py

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,7 @@
6868
TPR2020B2,
6969
INPCRD,
7070
TPR_gh_5145,
71+
TPR_linear_angle,
7172
)
7273
import MDAnalysis as mda
7374
from MDAnalysis.coordinates.TPR import TPRReader
@@ -84,6 +85,16 @@
8485
# NOTE: expected values are expressed in nm
8586
# units and converted to Angstrom in the body
8687
# of the test below, before assertions
88+
89+
# see gh-5361 for CO2/linear angle:
90+
(
91+
TPR_linear_angle,
92+
[1.250, 1.250, 1.250],
93+
[1.250, 1.260, 1.136],
94+
(3, 3),
95+
[0, 0, 0],
96+
[0, 0, 0],
97+
),
8798
# this case is an alanine dipeptide
8899
# with neural network potential active
89100
# and nonzero velocities
3.21 KB
Binary file not shown.

testsuite/MDAnalysisTests/datafiles.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -137,6 +137,7 @@
137137
"TPR2024_4",
138138
"TPR2025_0",
139139
"TPR2026_0",
140+
"TPR_linear_angle",
140141
"TPR510_bonded",
141142
"TPR2016_bonded",
142143
"TPR2018_bonded",
@@ -574,6 +575,7 @@
574575
TPR2024_4 = (_data_ref / "tprs/2lyz_gmx_2024_4.tpr").as_posix()
575576
TPR2025_0 = (_data_ref / "tprs/2lyz_gmx_2025_0.tpr").as_posix()
576577
TPR2026_0 = (_data_ref / "tprs/2lyz_gmx_2026_0.tpr").as_posix()
578+
TPR_linear_angle = (_data_ref / "tprs/Linear_Angle.tpr").as_posix()
577579
# double precision
578580
TPR455Double = (_data_ref / "tprs/drew_gmx_4.5.5.double.tpr").as_posix()
579581
TPR460 = (_data_ref / "tprs/ab42_gmx_4.6.tpr").as_posix()

testsuite/MDAnalysisTests/topology/test_tprparser.py

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,8 @@
5151
TPR2023_bonded, TPR2024_4_bonded,
5252
TPR2025_0_bonded, TPR2024_bonded,
5353
TPR2026_0_bonded,
54-
TPR_NNPOT_2025_0, TPR_NNPOT_2026_0)
54+
TPR_NNPOT_2025_0, TPR_NNPOT_2026_0,
55+
TPR_linear_angle)
5556
from numpy.testing import assert_equal
5657

5758
# fmt: on
@@ -438,3 +439,12 @@ def test_resids(resid_from_one, resid_addition):
438439
resids,
439440
err_msg="tpr_resid_from_one kwarg not switching resids",
440441
)
442+
443+
444+
def test_gh_5361():
445+
# GROMACS linear angle potential handling with CO2 example
446+
u = mda.Universe(TPR_linear_angle)
447+
actual = u.angles.to_indices()
448+
expected = [[1, 0, 2]]
449+
assert_equal(u.atoms.names, ["C", "O1", "O2"])
450+
assert_equal(actual, expected)

0 commit comments

Comments
 (0)