diff --git a/package/AUTHORS b/package/AUTHORS index 7d4d5e22f34..c2d41018dce 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -227,6 +227,7 @@ Chronological list of authors - Geongi Moon - Sumit Gupta - Heet Vekariya + - Johannes Stöckelmaier External code ------------- diff --git a/package/CHANGELOG b/package/CHANGELOG index ba292289607..7b6cbc709f3 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -15,7 +15,8 @@ The rules for this file: ------------------------------------------------------------------------------- ??/??/?? IAlibay, ianmkenney, PicoCentauri, pgbarletta, p-j-smith, - richardjgowers, lilyminium, ALescoulie, hmacdope, HeetVekariya + richardjgowers, lilyminium, ALescoulie, hmacdope, HeetVekariya, + JoStoe * 2.7.0 @@ -34,6 +35,7 @@ Fixes * Fix atom charge reading in PDBQT parser (Issue #4282, PR #4283) Enhancements + * Added a reader for GROMOS11 TRC trajectories (PR #4292 , Issue #4291) * Document the usage of NoDataError in its docstring (Issue #3901, PR #4359) * Refactor c_distances backend to have a cython .pxd header and expose in libmdanalysis (Issue #4315, PR #4324) diff --git a/package/MDAnalysis/coordinates/TRC.py b/package/MDAnalysis/coordinates/TRC.py new file mode 100644 index 00000000000..e1d29d2a84c --- /dev/null +++ b/package/MDAnalysis/coordinates/TRC.py @@ -0,0 +1,401 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + + +"""GROMOS11 trajectory reader --- :mod:`MDAnalysis.coordinates.TRC` +==================================================================== + +Reads coordinates, timesteps and box-sizes from GROMOS11 TRC trajectories. + +To load the trajectory into :class:`~MDAnalysis.core.universe.Universe`, +you need to provide topology information using a topology file such as +a PDB:: + + import MDAnalysis as mda + u = mda.Universe("topology.pdb", ["md_1.trc.gz","md_2.trc.gz"], + continuous=True) + +.. Note:: + The GROMOS trajectory format organizes its data in blocks. A block starts + with a blockname in capital letters (example: POSITIONRED) and ends with + a line containing only ''END''. Only the TITLE-block at the beginning of + each file is mandatory, others blocks can be chosen depending on the task. + + The trajectory format as well as all permitted blocks and their data are + documented in the GROMOS Manual Vol. 4, chapter 2 and 4. + The manual can be downloaded here: + https://gromos.net/gromos11_pdf_manuals/vol4.pdf + +This reader is designed to read the blocks "TIMESTEP", "POSITIONRED" and +"GENBOX" from the trajectory which covers the standard trajectories of most +simulations . If a trajectory includes other blocks, a warning is served +and the blocks are ignored. + +.. Note:: + MDAnalysis requires the blocks to be in the same order for each frame + and ignores non-supported blocks. + + + +Classes +------- + +.. autoclass:: TRCReader + :members: + +""" + +import pathlib +import errno +import warnings +import numpy as np + +from . import base +from .timestep import Timestep +from ..lib import util +from ..lib.util import cached, store_init_arguments + +import logging + +logger = logging.getLogger("MDAnalysis.coordinates.GROMOS11") + + +class TRCReader(base.ReaderBase): + """Coordinate reader for the GROMOS11 format""" + + format = "TRC" + units = {"time": "ps", "length": "nm"} + _Timestep = Timestep + + SUPPORTED_BLOCKS = ["TITLE", "TIMESTEP", "POSITIONRED", "GENBOX"] + NOT_SUPPORTED_BLOCKS = [ + "POSITION", + "REFPOSITION", + "VELOCITY", + "VELOCITYRED", + "FREEFORCE", + "FREEFORCERED", + "CONSFORCE", + "CONSFORCERED", + "SHAKEFAILPOSITION", + "SHAKEFAILPREVPOSITION", + "LATTICESHIFTS", + "COSDISPLACEMENTS", + "STOCHINT", + "NHCVARIABLES", + "ROTTRANSREFPOS", + "PERTDATA", + "DISRESEXPAVE", + "JVALUERESEXPAVE", + "JVALUERESEPS", + "JVALUEPERSCALE", + "ORDERPARAMRESEXPAVE", + "XRAYRESEXPAVE", + "LEUSBIAS", + "LEUSBIASBAS", + "ENERGY03", + "VOLUMEPRESSURE03", + "FREEENERGYDERIVS03", + "BFACTOR", + "AEDSS", + ] + + @store_init_arguments + def __init__(self, filename, convert_units=True, **kwargs): + super(TRCReader, self).__init__(filename, **kwargs) + + # GROMOS11 trajectories are usually either *.trc or *.trc.gz. + # (trj suffix can come up when trajectory obtained from clustering) + ext = pathlib.Path(self.filename).suffix + if (ext[1:] == "trc") or (ext[1:] == "trj"): + self.compressed = None + else: + self.compressed = ext[1:] + + self.trcfile = util.anyopen(self.filename) + self.convert_units = convert_units + + # Read and calculate some information about the trajectory + self.traj_properties = self._read_traj_properties() + + self._cache = {} + self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) + + self._reopen() + self.ts.dt = self.traj_properties["dt"] + + self._read_frame(0) + + @property + @cached("n_atoms") + def n_atoms(self): + """The number of atoms in one frame.""" + return self.traj_properties["n_atoms"] + + @property + @cached("n_frames") + def n_frames(self): + """The number of frames in the trajectory.""" + return self.traj_properties["n_frames"] + + def _frame_to_ts(self, frameDat, ts): + """Convert a frame to a :class:`TimeStep`""" + ts.frame = self._frame + ts.time = frameDat["time"] + + ts.data["time"] = frameDat["time"] + ts.data["step"] = frameDat["step"] + + ts.dimensions = frameDat["dimensions"] + ts.positions = frameDat["positions"] + + # Convert the units + if self.convert_units: + if ts.has_positions: + self.convert_pos_from_native(ts.positions) + if ts.dimensions is not None: + self.convert_pos_from_native(ts.dimensions[:3]) + + return ts + + def _read_traj_properties(self): + """ + * Reads the number of atoms per frame (n_atoms) + * Reads the number of frames (n_frames) + * Reads the startposition of the timestep block + for each frame (l_blockstart_offset) + """ + + traj_properties = {} + + # + # Check which of the supported blocks comes first in the trajectory + # + first_block = None + with util.anyopen(self.filename) as f: + for line in iter(f.readline, ""): + for blockname in TRCReader.SUPPORTED_BLOCKS: + if (blockname == line.strip()) and (blockname != "TITLE"): + # Save the name of the first non-title block + # in the trajectory file + first_block = blockname + + if first_block is not None: + break # First block found + + # + # Calculate meta-data of the trajectory + # + in_positionred_block = False + lastline_was_timestep = False + + atom_counter = 0 + n_atoms = 0 + frame_counter = 0 + + l_blockstart_offset = [] + l_timestep_timevalues = [] + + with util.anyopen(self.filename) as f: + for line in iter(f.readline, ""): + # + # First block of frame + # + if first_block == line.strip(): + l_blockstart_offset.append(f.tell() - len(line)) + frame_counter += 1 + + # + # Timestep-block + # + if "TIMESTEP" == line.strip(): + lastline_was_timestep = True + + elif lastline_was_timestep is True: + l_timestep_timevalues.append(float(line.split()[1])) + lastline_was_timestep = False + + # + # Coordinates-block + # + if "POSITIONRED" == line.strip(): + in_positionred_block = True + + elif (in_positionred_block is True) and (n_atoms == 0): + if len(line.split()) == 3: + atom_counter += 1 + + if ("END" == line.strip()) and (in_positionred_block is True): + n_atoms = atom_counter + in_positionred_block = False + + if frame_counter == 0: + raise ValueError( + "No supported blocks were found within the GROMOS trajectory!" + ) + + traj_properties["n_atoms"] = n_atoms + traj_properties["n_frames"] = frame_counter + traj_properties["l_blockstart_offset"] = l_blockstart_offset + + if len(l_timestep_timevalues) >= 2: + traj_properties["dt"] = l_timestep_timevalues[1] - l_timestep_timevalues[0] + else: + traj_properties["dt"] = 0 + warnings.warn( + "The trajectory does not contain TIMESTEP blocks!", UserWarning + ) + + return traj_properties + + def _read_GROMOS11_trajectory(self): + frameDat = {} + frameDat["step"] = int(self._frame) + frameDat["time"] = float(0.0) + frameDat["positions"] = None + frameDat["dimensions"] = None + self.periodic = False + + # Read the trajectory + f = self.trcfile + for line in iter(f.readline, ""): + if "TIMESTEP" == line.strip(): + tmp_step, tmp_time = f.readline().split() + frameDat["step"] = int(tmp_step) + frameDat["time"] = float(tmp_time) + + elif "POSITIONRED" == line.strip(): + tmp_buf = [] + while True: + coords_str = f.readline() + if "#" in coords_str: + continue + elif "END" in coords_str: + break + else: + tmp_buf.append(coords_str.split()) + + if np.array(tmp_buf).shape[0] == self.n_atoms: + frameDat["positions"] = np.asarray(tmp_buf, dtype=np.float64) + else: + raise ValueError( + "The trajectory contains the wrong number of atoms!" + ) + + elif "GENBOX" == line.strip(): + ntb_setting = int(f.readline()) + if ntb_setting == 0: + frameDat["dimensions"] = None + self.periodic = False + + elif ntb_setting in [1, 2]: + tmp_a, tmp_b, tmp_c = f.readline().split() + tmp_alpha, tmp_beta, tmp_gamma = f.readline().split() + frameDat["dimensions"] = [ + float(tmp_a), + float(tmp_b), + float(tmp_c), + float(tmp_alpha), + float(tmp_beta), + float(tmp_gamma), + ] + self.periodic = True + + gb_line3 = f.readline().split() + if np.sum(np.abs(np.array(gb_line3).astype(np.float64))) > 1e-10: + raise ValueError( + "This reader doesnt't support a shifted origin!" + ) + + gb_line4 = f.readline().split() + if np.sum(np.abs(np.array(gb_line4).astype(np.float64))) > 1e-10: + raise ValueError( + "This reader " + "doesnt't support " + "yawed, pitched or " + "rolled boxes!" + ) + + else: + raise ValueError( + "This reader doesn't support " + "truncated-octahedral " + "periodic boundary conditions" + ) + break + + elif any( + non_supp_bn in line for non_supp_bn in TRCReader.NOT_SUPPORTED_BLOCKS + ): + for non_supp_bn in TRCReader.NOT_SUPPORTED_BLOCKS: + if non_supp_bn == line.strip(): + warnings.warn( + non_supp_bn + " block is not supported!", UserWarning + ) + + return frameDat + + def _read_frame(self, i): + """read frame i""" + self._frame = i - 1 + + # Move position in file just (-2 byte) before the start of the block + self.trcfile.seek(self.traj_properties["l_blockstart_offset"][i] - 2, 0) + + return self._read_next_timestep() + + def _read_next_timestep(self): + self._frame += 1 + if self._frame >= self.n_frames: + raise IOError("Trying to go over trajectory limit") + + raw_framedata = self._read_GROMOS11_trajectory() + self._frame_to_ts(raw_framedata, self.ts) + + return self.ts + + def _reopen(self): + """Close and reopen the trajectory""" + self.close() + self.open_trajectory() + + def open_trajectory(self): + if self.trcfile is not None: + raise IOError(errno.EALREADY, "TRC file already opened", self.filename) + + # Reload trajectory file + self.trcfile = util.anyopen(self.filename) + + # Reset ts + self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) + + # Set _frame to -1, so next timestep is zero + self._frame = -1 + + return self.trcfile + + def close(self): + """Close the trc trajectory file if it was open.""" + if self.trcfile is not None: + self.trcfile.close() + self.trcfile = None diff --git a/package/MDAnalysis/coordinates/__init__.py b/package/MDAnalysis/coordinates/__init__.py index 7e1abaddcf8..9b6a7121bc9 100644 --- a/package/MDAnalysis/coordinates/__init__.py +++ b/package/MDAnalysis/coordinates/__init__.py @@ -235,6 +235,10 @@ class can choose an appropriate reader automatically. | | | | the first frame present will be read. | | [#a]_ | | | Module :mod:`MDAnalysis.coordinates.GRO` | +---------------+-----------+-------+------------------------------------------------------+ + | GROMOS11 | trc | r | Basic GROMOS11 trajectory format. | + | | | | Can read positions, box-sizes and timesteps. | + | | | | Module :mod:`MDAnalysis.coordinates.TRC` | + +---------------+-----------+-------+------------------------------------------------------+ | CHARMM | crd | r/w | "CARD" coordinate output from CHARMM; deals with | | CARD [#a]_ | | | either standard or EXTended format. Module | | | | | :mod:`MDAnalysis.coordinates.CRD` | @@ -772,6 +776,7 @@ class can choose an appropriate reader automatically. from . import PDB from . import PDBQT from . import PQR +from . import TRC from . import TRJ from . import TRR from . import H5MD diff --git a/package/MDAnalysis/coordinates/chain.py b/package/MDAnalysis/coordinates/chain.py index 33f5a80b642..bdfdeff3523 100644 --- a/package/MDAnalysis/coordinates/chain.py +++ b/package/MDAnalysis/coordinates/chain.py @@ -299,7 +299,8 @@ def __init__(self, filenames, skip=1, dt=None, continuous=False, # calculate new start_frames to have a time continuous trajectory. if continuous: - check_allowed_filetypes(self.readers, ['XTC', 'TRR', "LAMMPSDUMP"]) + check_allowed_filetypes(self.readers, ['XTC', 'TRR', 'LAMMPSDUMP', + 'TRC']) if np.any(np.array(n_frames) == 1): raise RuntimeError("ChainReader: Need at least two frames in " "every trajectory with continuous=True") diff --git a/package/doc/sphinx/source/documentation_pages/coordinates/TRC.rst b/package/doc/sphinx/source/documentation_pages/coordinates/TRC.rst new file mode 100644 index 00000000000..9312cc177df --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/coordinates/TRC.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.coordinates.TRC diff --git a/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst b/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst index 073ec44758f..af726d0f857 100644 --- a/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst @@ -36,6 +36,7 @@ provide the format in the keyword argument *format* to coordinates/PDBQT coordinates/PQR coordinates/TNG + coordinates/TRC coordinates/TRJ coordinates/TRR coordinates/TRZ diff --git a/testsuite/MDAnalysisTests/coordinates/test_trc.py b/testsuite/MDAnalysisTests/coordinates/test_trc.py new file mode 100644 index 00000000000..dd64f6d0a93 --- /dev/null +++ b/testsuite/MDAnalysisTests/coordinates/test_trc.py @@ -0,0 +1,284 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# +import pytest +import numpy as np +from numpy.testing import assert_allclose + +import MDAnalysis as mda +from MDAnalysisTests.datafiles import TRC_PDB_VAC, TRC_TRAJ1_VAC, TRC_TRAJ2_VAC +from MDAnalysisTests.datafiles import TRC_CLUSTER_VAC, TRC_EMPTY +from MDAnalysisTests.datafiles import TRC_GENBOX_ORIGIN, TRC_GENBOX_EULER +from MDAnalysisTests.datafiles import TRC_PDB_SOLV, TRC_TRAJ_SOLV +from MDAnalysisTests.datafiles import TRC_TRICLINIC_SOLV, TRC_TRUNCOCT_VAC + + +class TestTRCReaderVacuumBox: + @pytest.fixture(scope="class") + def TRC_U(self): + return mda.Universe( + TRC_PDB_VAC, [TRC_TRAJ1_VAC, TRC_TRAJ2_VAC], continuous=True + ) + + def test_initial_frame_is_0(self, TRC_U): + assert TRC_U.trajectory.ts.frame == 0 + + def test_trc_positions(self, TRC_U): + # first frame first particle + TRC_U.trajectory[0] + assert_allclose( + TRC_U.atoms.positions[0], [2.19782507, 24.65064345, 29.39783426] + ) + # fith frame first particle + TRC_U.trajectory[4] + assert_allclose(TRC_U.atoms.positions[0], [0.37026654, 22.78805010, 3.69695262]) + + def test_trc_dimensions(self, TRC_U): + assert TRC_U.trajectory[0].dimensions is None + + def test_trc_n_frames(self, TRC_U): + assert len(TRC_U.trajectory) == 6 + assert (TRC_U.trajectory.n_frames) == 6 + + def test_trc_n_atoms(self, TRC_U): + assert (TRC_U.trajectory.n_atoms) == 73 + + def test_trc_frame(self, TRC_U): + assert TRC_U.trajectory[0].frame == 0 + assert TRC_U.trajectory[4].frame == 4 + + def test_trc_time(self, TRC_U): + assert TRC_U.trajectory[0].time == 0 + assert TRC_U.trajectory[4].time == 80 + + def test_trc_dt(self, TRC_U): + time_array = np.array([ts.time for ts in TRC_U.trajectory]) + assert_allclose(time_array, np.arange(6) * 20.0) + + dt_array = np.diff(time_array) + assert_allclose(dt_array, np.full(5, 20.0)) + + def test_trc_data_step(self, TRC_U): + assert TRC_U.trajectory[0].data["step"] == 0 + assert TRC_U.trajectory[4].data["step"] == 10000 + + def test_periodic(self, TRC_U): + assert TRC_U.trajectory.periodic is False + + def test_rewind(self, TRC_U): + TRC_U.trajectory[0] + trc = TRC_U.trajectory + trc.next() + trc.next() + trc.next() + trc.next() + assert trc.ts.frame == 4, "trajectory.next() did not forward to frameindex 4" + trc.rewind() + assert trc.ts.frame == 0, "trajectory.rewind() failed to rewind to first frame" + + assert np.any(TRC_U.atoms.positions != 0), ( + "The atom positions are not populated" + ) + + def test_random_access(self, TRC_U): + TRC_U.trajectory[0] + pos0 = TRC_U.atoms.positions + TRC_U.trajectory.next() + TRC_U.trajectory.next() + pos2 = TRC_U.atoms.positions + + TRC_U.trajectory[0] + assert_allclose(TRC_U.atoms.positions, pos0) + + TRC_U.trajectory[2] + assert_allclose(TRC_U.atoms.positions, pos2) + + def test_read_frame_reopens(self, TRC_U): + TRC_U.trajectory._reopen() + TRC_U.trajectory[4] + assert TRC_U.trajectory.ts.frame == 4 + + +class TestTRCReaderSolvatedBox: + @pytest.fixture(scope="class") + def TRC_U(self): + return mda.Universe(TRC_PDB_SOLV, TRC_TRAJ_SOLV) + + def test_trc_n_atoms(self, TRC_U): + assert TRC_U.trajectory.n_atoms == 2797 + + def test_periodic(self, TRC_U): + assert TRC_U.trajectory.periodic is True + + def test_trc_dimensions(self, TRC_U): + ts = TRC_U.trajectory[1] + assert_allclose( + ts.dimensions, [30.54416298, 30.54416298, 30.54416298, 90.0, 90.0, 90.0] + ) + + def test_open_twice(self, TRC_U): + TRC_U.trajectory._reopen() + with pytest.raises(IOError): + TRC_U.trajectory.open_trajectory() + + +class TestTRCReaderSolvatedBoxNoConvert: + @pytest.fixture(scope="class") + def TRC_U(self): + return mda.Universe(TRC_PDB_SOLV, TRC_TRAJ_SOLV, convert_units=False) + + def test_trc_dimensions(self, TRC_U): + ts = TRC_U.trajectory[1] + assert_allclose( + ts.dimensions, [3.054416298, 3.054416298, 3.054416298, 90.0, 90.0, 90.0] + ) + + +class TestTRCReaderTriclinicBox: + @pytest.fixture(scope="class") + def TRC_U(self): + return mda.Universe(TRC_PDB_SOLV, TRC_TRICLINIC_SOLV) + + def test_trc_n_atoms(self, TRC_U): + assert TRC_U.trajectory.n_atoms == 2797 + + def test_periodic(self, TRC_U): + assert TRC_U.trajectory.periodic is True + + def test_trc_dimensions(self, TRC_U): + ts = TRC_U.trajectory[0] + assert_allclose( + ts.dimensions, + [ + 33.72394463, + 38.87568624, + 26.74177871, + 91.316862341, + 93.911031544, + 54.514000084, + ], + ) + + def test_trc_distances(self, TRC_U): + import MDAnalysis.analysis.distances as distances + import MDAnalysis.analysis.atomicdistances as atomicdistances + + ts = TRC_U.trajectory[0] + + atom_1a = TRC_U.select_atoms("resname LYSH and name O and resid 4") + atom_1b = TRC_U.select_atoms("resname ARG and name H and resid 3") + atom_1c = TRC_U.select_atoms("resname SOLV and name OW and resid 718") + atom_1d = TRC_U.select_atoms("resname SOLV and name OW and resid 305") + + atom_2a = TRC_U.select_atoms("resname VAL and name CG1 and resid 1") + atom_2b = TRC_U.select_atoms("resname SOLV and name OW and resid 497") + atom_2c = TRC_U.select_atoms("resname SOLV and name OW and resid 593") + atom_2d = TRC_U.select_atoms("resname SOLV and name OW and resid 497") + + ag1 = atom_1a + atom_1b + atom_1c + atom_1d + ag2 = atom_2a + atom_2b + atom_2c + atom_2d + + dist_A = distances.dist(ag1, ag2, box=ts.dimensions)[2] # return distance + dist_B = atomicdistances.AtomicDistances(ag1, ag2, pbc=True).run().results[0] + + assert_allclose( + dist_A, [5.9488481, 4.4777278, 20.8165518, 7.5727112], rtol=1e-06 + ) # Reference values calculated with gromos++ tser + assert_allclose( + dist_B, [5.9488481, 4.4777278, 20.8165518, 7.5727112], rtol=1e-06 + ) # Reference values calculated with gromos++ tser + + +class TestTRCReaderTruncOctBox: + def test_universe(self): + with pytest.raises(ValueError, match="truncated-octahedral"): + mda.Universe(TRC_PDB_VAC, TRC_TRUNCOCT_VAC) + + +class TestTRCGenboxOrigin: + def test_universe(self): + with pytest.raises( + ValueError, match="doesnt't support a shifted origin!" + ): + mda.Universe(TRC_PDB_VAC, TRC_GENBOX_ORIGIN) + + +class TestTRCGenboxEuler: + def test_universe(self): + with pytest.raises( + ValueError, + match=("doesnt't support yawed, pitched or rolled boxes!"), + ): + mda.Universe(TRC_PDB_VAC, TRC_GENBOX_EULER) + + +class TestTRCEmptyFile: + def test_universe(self): + with pytest.raises( + ValueError, + match=("No supported blocks were found within the GROMOS trajectory!"), + ): + mda.Universe(TRC_PDB_VAC, TRC_EMPTY) + + +class TestTRCReaderClusterTrajectory: + @pytest.fixture(scope="class") + def TRC_U(self): + with pytest.warns(UserWarning) as record_of_warnings: + TRC_U = mda.Universe(TRC_PDB_VAC, TRC_CLUSTER_VAC, format="TRC") + + warning_strings = [ + "The trajectory does not contain TIMESTEP blocks!", + "POSITION block is not supported!", + ] + + warning_strings_found = 0 + for w in record_of_warnings: + if str(w.message) in warning_strings: + warning_strings_found += 1 + + assert warning_strings_found == 2 + return TRC_U + + def test_trc_n_atoms(self, TRC_U): + assert TRC_U.trajectory.n_atoms == 73 + + def test_periodic(self, TRC_U): + assert TRC_U.trajectory.periodic is False + + def test_trc_dimensions(self, TRC_U): + ts = TRC_U.trajectory[1] + assert ts.dimensions is None + + def test_trc_n_frames(self, TRC_U): + assert len(TRC_U.trajectory) == 3 + assert TRC_U.trajectory.n_frames == 3 + + def test_trc_frame(self, TRC_U): + with pytest.warns(UserWarning, match="POSITION block is not supported!"): + assert TRC_U.trajectory[0].frame == 0 + assert TRC_U.trajectory[2].frame == 2 + + def test_trc_time(self, TRC_U): + with pytest.warns(UserWarning, match="POSITION block is not supported!"): + assert TRC_U.trajectory[0].time == 0 + assert TRC_U.trajectory[2].time == 0 diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_cluster_vac.trj.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_cluster_vac.trj.gz new file mode 100644 index 00000000000..46312316a49 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_cluster_vac.trj.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_empty.trc b/testsuite/MDAnalysisTests/data/gromos11/gromos11_empty.trc new file mode 100644 index 00000000000..8d1c8b69c3f --- /dev/null +++ b/testsuite/MDAnalysisTests/data/gromos11/gromos11_empty.trc @@ -0,0 +1 @@ + diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_genbox_euler.trc.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_genbox_euler.trc.gz new file mode 100644 index 00000000000..1cf1950f08c Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_genbox_euler.trc.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_genbox_origin.trc.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_genbox_origin.trc.gz new file mode 100644 index 00000000000..7f571f4a0b0 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_genbox_origin.trc.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_solv.pdb.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_solv.pdb.gz new file mode 100644 index 00000000000..2eb6122c24b Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_solv.pdb.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_solv.trc.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_solv.trc.gz new file mode 100644 index 00000000000..45b98f6a81e Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_solv.trc.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac.pdb.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac.pdb.gz new file mode 100644 index 00000000000..8eb622bdd27 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac.pdb.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac_1.trc.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac_1.trc.gz new file mode 100644 index 00000000000..84f1edc7f0a Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac_1.trc.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac_2.trc.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac_2.trc.gz new file mode 100644 index 00000000000..f9042a2868e Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_traj_vac_2.trc.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_triclinic_solv.trc.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_triclinic_solv.trc.gz new file mode 100644 index 00000000000..f5ea5064c5e Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_triclinic_solv.trc.gz differ diff --git a/testsuite/MDAnalysisTests/data/gromos11/gromos11_truncOcta_vac.trc.gz b/testsuite/MDAnalysisTests/data/gromos11/gromos11_truncOcta_vac.trc.gz new file mode 100644 index 00000000000..79c0de2fb16 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/gromos11/gromos11_truncOcta_vac.trc.gz differ diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index 90cf2d17d18..8ae48dfa41d 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -211,6 +211,12 @@ "legacy_DCD_NAMD_coords", # frame 0 read in for SiN_tric_namd.dcd using legacy DCD reader "legacy_DCD_c36_coords", # frames 1 and 4 read in for tip125_tric_C36.dcd using legacy DCD reader "GSD", "GSD_bonds", "GSD_long", + "TRC_PDB_VAC", "TRC_TRAJ1_VAC", "TRC_TRAJ2_VAC", # 2x 3 frames of vacuum trajectory from GROMOS11 tutorial + "TRC_CLUSTER_VAC", # three frames without TIMESTEP and GENBOX block but with unsupported POSITION block + "TRC_TRICLINIC_SOLV", "TRC_TRUNCOCT_VAC", + "TRC_GENBOX_ORIGIN", "TRC_GENBOX_EULER", + "TRC_EMPTY", # Empty file containing only one space + "TRC_PDB_SOLV", "TRC_TRAJ_SOLV", # 2 frames of solvated trajectory from GROMOS11 tutorial "GRO_MEMPROT", "XTC_MEMPROT", # YiiP transporter in POPE:POPG lipids with Na+, Cl-, Zn2+ dummy model without water "DihedralArray", "DihedralsArray", # time series of single dihedral "RamaArray", "GLYRamaArray", # time series of phi/psi angles @@ -599,6 +605,18 @@ GSD_bonds = (_data_ref / 'example_bonds.gsd').as_posix() GSD_long = (_data_ref / 'example_longer.gsd').as_posix() +TRC_PDB_VAC = (_data_ref / 'gromos11/gromos11_traj_vac.pdb.gz').as_posix() +TRC_TRAJ1_VAC = (_data_ref / 'gromos11/gromos11_traj_vac_1.trc.gz').as_posix() +TRC_TRAJ2_VAC = (_data_ref / 'gromos11/gromos11_traj_vac_2.trc.gz').as_posix() +TRC_PDB_SOLV = (_data_ref / 'gromos11/gromos11_traj_solv.pdb.gz').as_posix() +TRC_TRAJ_SOLV = (_data_ref / 'gromos11/gromos11_traj_solv.trc.gz').as_posix() +TRC_CLUSTER_VAC = (_data_ref / 'gromos11/gromos11_cluster_vac.trj.gz').as_posix() +TRC_TRICLINIC_SOLV = (_data_ref / 'gromos11/gromos11_triclinic_solv.trc.gz').as_posix() +TRC_TRUNCOCT_VAC = (_data_ref / 'gromos11/gromos11_truncOcta_vac.trc.gz').as_posix() +TRC_GENBOX_ORIGIN = (_data_ref / 'gromos11/gromos11_genbox_origin.trc.gz').as_posix() +TRC_GENBOX_EULER = (_data_ref / 'gromos11/gromos11_genbox_euler.trc.gz').as_posix() +TRC_EMPTY = (_data_ref / 'gromos11/gromos11_empty.trc').as_posix() + DihedralArray = (_data_ref / 'adk_oplsaa_dihedral.npy').as_posix() DihedralsArray = (_data_ref / 'adk_oplsaa_dihedral_list.npy').as_posix() RamaArray = (_data_ref / 'adk_oplsaa_rama.npy').as_posix() diff --git a/testsuite/MDAnalysisTests/parallelism/test_multiprocessing.py b/testsuite/MDAnalysisTests/parallelism/test_multiprocessing.py index 638132d457c..6b29aec6a57 100644 --- a/testsuite/MDAnalysisTests/parallelism/test_multiprocessing.py +++ b/testsuite/MDAnalysisTests/parallelism/test_multiprocessing.py @@ -56,6 +56,7 @@ PDB, PDB_small, PDB_multiframe, PDBQT_input, PQR, + TRC_PDB_VAC, TRC_TRAJ1_VAC, TRC_TRAJ2_VAC, TRR, TRJ, TRZ, @@ -83,6 +84,8 @@ (GRO, [GRO, GRO, GRO, GRO, GRO]), (PDB, [PDB, PDB, PDB, PDB, PDB]), (GRO, [XTC, XTC]), + (TRC_PDB_VAC, TRC_TRAJ1_VAC), + (TRC_PDB_VAC, [TRC_TRAJ1_VAC, TRC_TRAJ2_VAC]), ]) def u(request): if len(request.param) == 1: @@ -193,6 +196,8 @@ def test_creating_multiple_universe_without_offset(temp_xtc, ncopies=3): ('NCDF', NCDF, dict()), ('TXYZ', TXYZ, dict()), ('memory', np.arange(60).reshape(2, 10, 3).astype(np.float64), dict()), + ('TRC', TRC_TRAJ1_VAC, dict()), + ('CHAIN', [TRC_TRAJ1_VAC, TRC_TRAJ2_VAC], dict()), ('CHAIN', [GRO, GRO, GRO], dict()), ('CHAIN', [PDB, PDB, PDB], dict()), ('CHAIN', [XTC, XTC, XTC], dict()), diff --git a/testsuite/setup.py b/testsuite/setup.py index 2e026ed5b9b..154b7dba99b 100755 --- a/testsuite/setup.py +++ b/testsuite/setup.py @@ -149,6 +149,7 @@ def run(self): 'data/Amber/*.trj', 'data/Amber/*.mdcrd', 'data/Amber/*.ncdf', 'data/Amber/*.nc', 'data/Amber/*.inpcrd', + 'data/gromos11/*.gz', 'data/gromos11/*.trc', 'data/*.pqr', 'data/*.pdbqt', 'data/*.bz2', 'data/*.gz', 'data/*.ent', 'data/*.fasta',