diff --git a/package/CHANGELOG b/package/CHANGELOG index 979aa607a16..91613de6c17 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -23,6 +23,10 @@ The rules for this file: * 2.0.0 Fixes + * Only GRO files with unit cells defined with 3 or 9 entries are now + supported (Issue #3305) + * Ensures that reading and writing of GRO files with missing unit cell + information is consistently treated (Issue #3305) * Fixed sometimes wrong sorting of atoms into fragments when unwrapping (Issue #3352, PR #3353) * Fixed missing array flatten preventing PCA from running when mean positions diff --git a/package/MDAnalysis/coordinates/GRO.py b/package/MDAnalysis/coordinates/GRO.py index 2d354a07832..9a4a12f95be 100644 --- a/package/MDAnalysis/coordinates/GRO.py +++ b/package/MDAnalysis/coordinates/GRO.py @@ -155,8 +155,19 @@ class GROReader(base.SingleFrameReaderBase): .. note:: This Reader will only read the first frame present in a file. + + .. note:: + GRO files with zeroed 3 entry unit cells (i.e. ``0.0 0.0 0.0``) + are read as missing unit cell information. In these cases ``dimensions`` + will be set to ``None``. + .. versionchanged:: 0.11.0 Frames now 0-based instead of 1-based + .. versionchanged:: 2.0.0 + Reader now only parses boxes defined with 3 or 9 fields. + Reader now reads a 3 entry zero unit cell (i.e. ``[0, 0, 0]``) as a + being without dimension information (i.e. will the timestep dimension + to ``None``). """ format = 'GRO' units = {'time': None, 'length': 'nm', 'velocity': 'nm/ps'} @@ -214,15 +225,19 @@ def _read_first_frame(self): if len(unitcell) == 3: # special case: a b c --> (a 0 0) (b 0 0) (c 0 0) # see docstring above for format (!) - self.ts.dimensions = np.r_[unitcell, [90., 90., 90.]] + # Treat empty 3 entry boxes as not having a unit cell + if np.allclose(unitcell, [0., 0., 0.]): + wmsg = ("Empty box [0., 0., 0.] found - treating as missing " + "unit cell. Dimensions set to `None`.") + warnings.warn(wmsg) + self.ts.dimensions = None + else: + self.ts.dimensions = np.r_[unitcell, [90., 90., 90.]] elif len(unitcell) == 9: self.ts.dimensions = _gmx_to_dimensions(unitcell) - else: # or maybe raise an error for wrong format?? - warnings.warn("GRO unitcell has neither 3 nor 9 entries --- might be wrong.") - # fill linearly ... not sure about this - unitcell2 = np.zeros(9) - unitcell2[:len(unitcell)] = unitcell - self.ts.dimensions = _gmx_to_dimensions(unitcell2) + else: # raise an error for wrong format + errmsg = "GRO unitcell has neither 3 nor 9 entries." + raise ValueError(errmsg) if self.convert_units: self.convert_pos_from_native(self.ts._pos) # in-place ! @@ -258,9 +273,14 @@ class GROWriter(base.WriterBase): - resnames (defaults to 'UNK') - resids (defaults to '1') - Note - ---- - The precision is hard coded to three decimal places + + .. note:: + The precision is hard coded to three decimal places. + + + .. note:: + When dimensions are missing (i.e. set to `None`), a zero width + unit cell box will be written (i.e. [0.0, 0.0, 0.0]). .. versionchanged:: 0.11.0 @@ -274,6 +294,9 @@ class GROWriter(base.WriterBase): .. versionchanged:: 0.18.0 Added `reindex` keyword argument to allow original atom ids to be kept. + .. versionchanged:: 2.0.0 + Raises a warning when writing timestep with missing dimension + information (i.e. set to ``None``). """ format = 'GRO' @@ -455,6 +478,9 @@ def write(self, obj): if (ag.dimensions is None or np.allclose(ag.dimensions[3:], [90., 90., 90.])): if ag.dimensions is None: + wmsg = ("missing dimension - setting unit cell to zeroed " + "box [0., 0., 0.]") + warnings.warn(wmsg) box = np.zeros(3) else: box = self.convert_pos_to_native( diff --git a/testsuite/MDAnalysisTests/coordinates/test_gro.py b/testsuite/MDAnalysisTests/coordinates/test_gro.py index 0b01daca585..e25bf969fc5 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_gro.py +++ b/testsuite/MDAnalysisTests/coordinates/test_gro.py @@ -20,6 +20,8 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # +from io import StringIO + import MDAnalysis as mda import numpy as np from MDAnalysis.coordinates.GRO import GROReader, GROWriter @@ -37,6 +39,7 @@ GRO_large, two_water_gro_multiframe, GRO_huge_box, + PDB_closed, ) from numpy.testing import ( assert_almost_equal, @@ -485,8 +488,44 @@ def test_multiframe_gro(): assert len(u.trajectory) == 1 assert_equal(u.dimensions, np.array([100, 100, 100, 90, 90, 90], dtype=np.float32)) + def test_huge_box_gro(): u = mda.Universe(GRO_huge_box) assert_equal(u.dimensions, np.array([4.e+05, 4.e+05, 4.e+05, 90, 90, 90], dtype=np.float32)) + + +gro_no_dims = """\ +Single Atom no dims GRO +1 + 1SOL OW 1 1.000 1.000 1.000 +""" + + +@pytest.mark.parametrize('dims', [1, 2, 4, 5, 6, 7, 8]) +def test_bad_box(dims): + cell = ' '.join([str(float(i)) for i in range(dims)]) + grofile = gro_no_dims + cell + + errmsg = "GRO unitcell has neither 3 nor 9 entries." + with pytest.raises(ValueError, match=errmsg): + u = mda.Universe(StringIO(grofile), format='gro') + + +def test_gro_empty_box_write_read(tmpdir): + # Issue #3305 - ensure that read/write deals with None dimensions the same + u = mda.Universe(PDB_closed) + + # Check lack of dimensions + assert u.dimensions is None + + with tmpdir.as_cwd(): + wmsg = " setting unit cell to zeroed box" + with pytest.warns(UserWarning, match=wmsg): + u.atoms.write('test.gro') + + wmsg = "treating as missing unit cell" + with pytest.warns(UserWarning, match=wmsg): + u2 = mda.Universe('test.gro') + assert u2.dimensions is None