diff --git a/.appveyor.yml b/.appveyor.yml index 76339b3fb2f..6f07a88fc36 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -10,7 +10,7 @@ cache: environment: global: CONDA_CHANNELS: conda-forge - CONDA_DEPENDENCIES: pip setuptools wheel cython mock six biopython networkx joblib matplotlib scipy vs2015_runtime pytest mmtf-python GridDataFormats hypothesis pytest-cov codecov + CONDA_DEPENDENCIES: pip setuptools wheel cython mock six biopython networkx joblib matplotlib scipy vs2015_runtime pytest mmtf-python GridDataFormats hypothesis pytest-cov codecov chemfiles PIP_DEPENDENCIES: gsd==1.9.3 duecredit parmed DEBUG: "False" MINGW_64: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin diff --git a/.travis.yml b/.travis.yml index 2818163d465..1249fc9a70d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -28,7 +28,7 @@ env: - SETUP_CMD="${PYTEST_FLAGS}" - BUILD_CMD="pip install -e package/ && (cd testsuite/ && python setup.py build)" - CONDA_MIN_DEPENDENCIES="mmtf-python mock six biopython networkx cython matplotlib scipy griddataformats hypothesis gsd codecov" - - CONDA_DEPENDENCIES="${CONDA_MIN_DEPENDENCIES} seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12" + - CONDA_DEPENDENCIES="${CONDA_MIN_DEPENDENCIES} seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12 chemfiles" - CONDA_CHANNELS='biobuilds conda-forge' - CONDA_CHANNEL_PRIORITY=True - PIP_DEPENDENCIES="duecredit parmed" @@ -66,7 +66,7 @@ matrix: INSTALL_HOLE="false" - env: NAME="python 2.7" - PYTHON_VERSION=2.7 + PYTHON_VERSION=2.7 CODECOV="true" NUMPY_VERSION=1.16 SETUP_CMD="${PYTEST_FLAGS} --cov=MDAnalysis" diff --git a/package/AUTHORS b/package/AUTHORS index 6168416987b..05f38057159 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -128,6 +128,7 @@ Chronological list of authors 2020 - Charlie Cook - Yuanyu Chang + - Guillaume Fraux - Ivan Hristov - Michael Quevillon - Hao Tian @@ -136,7 +137,6 @@ Chronological list of authors - Shubham Sharma - External code ------------- diff --git a/package/CHANGELOG b/package/CHANGELOG index ccf66efba15..17ef7119be3 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -16,7 +16,7 @@ The rules for this file: mm/dd/yy richardjgowers, kain88-de, lilyminium, p-j-smith, bdice, joaomcteixeira, PicoCentauri, davidercruz, jbarnoud, RMeli, IAlibay, mtiberti, CCook96, Yuan-Yu, xiki-tempula, HTian1997, Iv-Hristov, hmacdope, AnshulAngaria, - ss62171 + ss62171, Luthaf * 0.21.0 @@ -66,6 +66,7 @@ Fixes * Fixed example docs for polymer persistence length (#2582) Enhancements + * Added ability to use Chemfiles as a trajectory reader backend (PR #1862) * New analysis.hole2 module for HOLE interfacing. Contains a function (hole) for running HOLE on singular PDB files and class (HoleAnalysis) for trajectories (PR #2523) diff --git a/package/MDAnalysis/coordinates/__init__.py b/package/MDAnalysis/coordinates/__init__.py index f724e257a33..5bd23b82bf8 100644 --- a/package/MDAnalysis/coordinates/__init__.py +++ b/package/MDAnalysis/coordinates/__init__.py @@ -716,6 +716,7 @@ class can choose an appropriate reader automatically. from . import base from .core import reader, writer from . import chain +from . import chemfiles from . import CRD from . import DCD from . import DLPoly diff --git a/package/MDAnalysis/coordinates/chemfiles.py b/package/MDAnalysis/coordinates/chemfiles.py new file mode 100644 index 00000000000..c7fbe842c2c --- /dev/null +++ b/package/MDAnalysis/coordinates/chemfiles.py @@ -0,0 +1,393 @@ +# -*- 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. +# +# 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 +# +""" +Reading trajectory with Chemfiles --- :mod:`MDAnalysis.coordinates.chemfiles` +============================================================================= + +Classes to read and write files using the `chemfiles`_ library. This library +provides C++ implementation of multiple formats readers and writers, the full +list if available `here `_. + +.. _chemfiles: https://chemfiles.org +.. _formats: http://chemfiles.org/chemfiles/latest/formats.html + +Classes +------- + +.. autoclass:: ChemfilesReader + +.. autoclass:: ChemfilesWriter + +""" +from __future__ import absolute_import, division +from __future__ import print_function, unicode_literals + +import numpy as np +from six import string_types +from distutils.version import LooseVersion +import warnings + +from . import base, core + +try: + import chemfiles +except ImportError: + HAS_CHEMFILES = False +else: + HAS_CHEMFILES = True + + +MIN_CHEMFILES_VERSION = LooseVersion("0.9") +MAX_CHEMFILES_VERSION = LooseVersion("0.10") + + +def check_chemfiles_version(): + """Check an appropriate Chemfiles is available + + Returns True if a usable chemfiles version is available + + .. versionadded:: 1.0.0 + """ + if not HAS_CHEMFILES: + warnings.warn( + "No Chemfiles package found. " + "Try installing with 'pip install chemfiles'" + ) + return False + version = LooseVersion(chemfiles.__version__) + wrong = version < MIN_CHEMFILES_VERSION or version >= MAX_CHEMFILES_VERSION + if wrong: + warnings.warn( + "unsupported Chemfiles version {}, we need a version >{} and <{}" + .format(version, MIN_CHEMFILES_VERSION, MAX_CHEMFILES_VERSION) + ) + return not wrong + + +class ChemfilesReader(base.ReaderBase): + """Coordinate reader using chemfiles. + + The file format to used is guessed based on the file extension. If no + matching format is found, a ``ChemfilesError`` is raised. It is also + possible to manually specify the format to use for a given file. + + .. versionadded:: 1.0.0 + """ + format = 'chemfiles' + units = {'time': 'fs', 'length': 'Angstrom'} + + def __init__(self, filename, chemfiles_format="", **kwargs): + """ + Parameters + ---------- + filename : chemfiles.Trajectory or str + the chemfiles object to read or filename to read + chemfiles_format : str (optional) + if *filename* was a string, use the given format name instead of + guessing from the extension. The `list of supported formats + `_ and the associated names is available in chemfiles + documentation. + **kwargs : dict + General reader arguments. + + + .. _formats: http://chemfiles.org/chemfiles/latest/formats.html + """ + if not check_chemfiles_version(): + raise RuntimeError("Please install Chemfiles > {}" + "".format(MIN_CHEMFILES_VERSION)) + super(ChemfilesReader, self).__init__(filename, **kwargs) + self._format = chemfiles_format + self._cached_n_atoms = None + self._open() + self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) + self.next() + + @staticmethod + def _format_hint(thing): + """Can this Reader read *thing*""" + # nb, filename strings can still get passed through if + # format='CHEMFILES' is used + return HAS_CHEMFILES and isinstance(thing, chemfiles.Trajectory) + + def _open(self): + self._closed = False + self._step = 0 + self._frame = -1 + if isinstance(self.filename, chemfiles.Trajectory): + self._file = self.filename + else: + self._file = chemfiles.Trajectory(self.filename, 'r', self._format) + + def close(self): + """close reader""" + if not self._closed: + self._file.close() + self._closed = True + + @property + def n_frames(self): + """number of frames in trajectory""" + return self._file.nsteps + + @property + def n_atoms(self): + """number of atoms in the first frame of the trajectory""" + if self._cached_n_atoms is None: + self._cached_n_atoms = len(self._file.read_step(0).atoms) + return self._cached_n_atoms + + def _reopen(self): + """reopen trajectory""" + self.close() + self._open() + + def _read_frame(self, i): + """read frame i""" + self._step = i + return self._read_next_timestep() + + def _read_next_timestep(self, ts=None): + """copy next frame into timestep""" + if self._step >= self.n_frames: + raise IOError('trying to go over trajectory limit') + if ts is None: + ts = self.ts + self.ts = ts + frame = self._file.read_step(self._step) + self._frame_to_ts(frame, ts) + ts.frame = frame.step + self._step += 1 + return ts + + def _frame_to_ts(self, frame, ts): + """convert a chemfiles frame to a :class:`TimeStep`""" + if len(frame.atoms) != self.n_atoms: + raise IOError( + "The number of atom changed in the trajectory. " + "This is not supported by MDAnalysis." + ) + + ts.dimensions[:] = frame.cell.lengths + frame.cell.angles + ts.positions[:] = frame.positions[:] + if frame.has_velocities(): + ts.has_velocities = True + ts.velocities[:] = frame.velocities[:] + + def Writer(self, filename, n_atoms=None, **kwargs): + """Return writer for trajectory format""" + if n_atoms is None: + n_atoms = self.n_atoms + return ChemfilesWriter(filename, n_atoms, **kwargs) + + +class ChemfilesWriter(base.WriterBase): + """ + Coordinate writer using chemfiles. + + The file format to used is guessed based on the file extension. If no + matching format is found, a ``ChemfilesError`` is raised. It is also + possible to manually specify the format to use for a given file. + + Chemfiles support writting to files with varying number of atoms if the + underlying format support it. This is the case of most of text-based + formats. + + .. versionadded:: 1.0.0 + """ + format = 'chemfiles' + multiframe = True + + # chemfiles mostly[1] uses these units for the in-memory representation, + # and convert into the file format units when writing. + # + # [1] mostly since some format don't have a specified unit + # (XYZ for example), so then chemfiles just assume they are in A and fs. + units = {'time': 'fs', 'length': 'Angstrom'} + + def __init__(self, filename, n_atoms=0, mode="w", chemfiles_format="", topology=None, **kwargs): + """ + Parameters + ---------- + filename: str + filename of trajectory file. + n_atoms: int + number of atoms in the trajectory to write. This value is not + used and can vary during trajectory, if the underlying format + support it + mode : str (optional) + file opening mode: use 'a' to append to an existing file or 'w' to + create a new file + chemfiles_format : str (optional) + use the given format name instead of guessing from the extension. + The `list of supported formats `_ and the associated names + is available in chemfiles documentation. + topology : Universe or AtomGroup (optional) + use the topology from this :class:`~MDAnalysis.core.groups.AtomGroup` + or :class:`~MDAnalysis.core.universe.Universe` to write all the + timesteps to the file + **kwargs : dict + General writer arguments. + + + .. _formats: http://chemfiles.org/chemfiles/latest/formats.html + """ + if not check_chemfiles_version(): + raise RuntimeError("Please install Chemfiles > {}" + "".format(MIN_CHEMFILES_VERSION)) + self.filename = filename + self.n_atoms = n_atoms + if mode != "a" and mode != "w": + raise IOError("Expected 'a' or 'w' as mode in ChemfilesWriter") + self._file = chemfiles.Trajectory(self.filename, mode, chemfiles_format) + self._closed = False + if topology is not None: + if isinstance(topology, string_types): + self._file.set_topology(topology) + else: + topology = self._topology_to_chemfiles(topology, n_atoms) + self._file.set_topology(topology) + + def close(self): + """Close the trajectory file and finalize the writing""" + if hasattr(self, "_closed") and not self._closed: + self._file.close() + self._closed = True + + def write(self, obj): + """Write object `obj` at current trajectory frame to file. + + Topology for the output is taken from the `obj` or default to the value + of the `topology` keyword supplied to the :class:`ChemfilesWriter` + constructor. + + If `obj` contains velocities, and the underlying format supports it, the + velocities are writen to the file. Writing forces is unsupported at the + moment. + + Parameters + ---------- + obj : Universe or AtomGroup + The :class:`~MDAnalysis.core.groups.AtomGroup` or + :class:`~MDAnalysis.core.universe.Universe` to write. + """ + if hasattr(obj, "atoms"): + if hasattr(obj, 'universe'): + # For AtomGroup and children (Residue, ResidueGroup, Segment) + ts_full = obj.universe.trajectory.ts + if ts_full.n_atoms == obj.atoms.n_atoms: + ts = ts_full + else: + # Only populate a time step with the selected atoms. + ts = ts_full.copy_slice(atoms.indices) + elif hasattr(obj, 'trajectory'): + # For Universe only --- get everything + ts = obj.trajectory.ts + else: + if isinstance(obj, base.Timestep): + ts = obj + topology = None + else: + raise TypeError("No Timestep found in obj argument") + + frame = self._timestep_to_chemfiles(ts) + frame.topology = self._topology_to_chemfiles(obj, len(frame.atoms)) + self._file.write(frame) + + def write_next_timestep(self, ts): + """Write timestep object into trajectory. + + Parameters + ---------- + ts: TimeStep + """ + frame = self._timestep_to_chemfiles(ts) + self._file.write(frame) + + def _timestep_to_chemfiles(self, ts): + """ + Convert a Timestep to a chemfiles Frame + """ + # TODO: CONVERTERS? + frame = chemfiles.Frame() + frame.resize(ts.n_atoms) + if ts.has_positions: + frame.positions[:] = ts.positions[:] + if ts.has_velocities: + frame.add_velocities() + frame.velocities[:] = ts.velocities[:] + frame.cell = chemfiles.UnitCell(*ts.dimensions) + return frame + + def _topology_to_chemfiles(self, obj, n_atoms): + """ + Convert an AtomGroup or an Universe to a chemfiles Topology + """ + topology = chemfiles.Topology() + if not hasattr(obj, "atoms"): + # use an empty topology + topology.resize(n_atoms) + return topology + + # (1) add all atoms to the topology + residues = {} + for atom in obj.atoms: + name = getattr(atom, 'name', "") + type = getattr(atom, 'type', name) + chemfiles_atom = chemfiles.Atom(name, type) + + if hasattr(atom, 'altLoc'): + chemfiles_atom["altloc"] = str(atom.altLoc) + + if hasattr(atom, 'segid'): + chemfiles_atom["segid"] = str(atom.segid) + + if hasattr(atom, 'segindex'): + chemfiles_atom["segindex"] = int(atom.segindex) + + if hasattr(atom, 'resid'): + resname = getattr(atom, 'resname', "") + if atom.resid not in residues.keys(): + residues[atom.resid] = chemfiles.Residue(resname, atom.resid) + residue = residues[atom.resid] + + atom_idx = len(topology.atoms) + residue.atoms.append(atom_idx) + + if hasattr(atom, "record_type"): + # set corresponding chemfiles residue property + if atom.record_type == "ATOM": + residue["is_standard_pdb"] = True + else: + residue["is_standard_pdb"] = False + + topology.atoms.append(chemfiles_atom) + + # (2) add residues to the topology + for residue in residues.values(): + topology.residues.append(residue) + + # (3) add bonds to the topology + for bond in getattr(obj, "bonds", []): + topology.add_bond(bond.atoms[0].ix, bond.atoms[1].ix) + + return topology diff --git a/package/doc/sphinx/source/documentation_pages/coordinates/chemfiles.rst b/package/doc/sphinx/source/documentation_pages/coordinates/chemfiles.rst new file mode 100644 index 00000000000..13016a9a621 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/coordinates/chemfiles.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.coordinates.chemfiles diff --git a/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst b/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst index 21f1ab32cb8..e74387a4bd6 100644 --- a/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst @@ -41,6 +41,7 @@ provide the format in the keyword argument *format* to coordinates/XTC coordinates/XYZ coordinates/memory + coordinates/chemfiles coordinates/null .. rubric:: Coordinate core modules diff --git a/testsuite/MDAnalysisTests/coordinates/base.py b/testsuite/MDAnalysisTests/coordinates/base.py index eb509c82193..ba62bd0bd36 100644 --- a/testsuite/MDAnalysisTests/coordinates/base.py +++ b/testsuite/MDAnalysisTests/coordinates/base.py @@ -236,14 +236,14 @@ def test_double_close(self, reader): def test_get_writer_1(self, ref, reader): with tempdir.in_tempdir(): - outfile = 'test-writer' + ref.ext + outfile = 'test-writer.' + ref.ext with reader.Writer(outfile) as W: assert_equal(isinstance(W, ref.writer), True) assert_equal(W.n_atoms, reader.n_atoms) def test_get_writer_2(self, ref, reader): with tempdir.in_tempdir(): - outfile = 'test-writer' + ref.ext + outfile = 'test-writer.' + ref.ext with reader.Writer(outfile, n_atoms=100) as W: assert_equal(isinstance(W, ref.writer), True) assert_equal(W.n_atoms, 100) diff --git a/testsuite/MDAnalysisTests/coordinates/test_chemfiles.py b/testsuite/MDAnalysisTests/coordinates/test_chemfiles.py new file mode 100644 index 00000000000..fa9a7ef5f4c --- /dev/null +++ b/testsuite/MDAnalysisTests/coordinates/test_chemfiles.py @@ -0,0 +1,256 @@ +# -*- 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. +# +# 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 +# +from __future__ import absolute_import + +import numpy as np +import pytest + +import MDAnalysis as mda +from MDAnalysis.coordinates.chemfiles import ( + ChemfilesReader, ChemfilesWriter, check_chemfiles_version, +) + +from MDAnalysisTests import datafiles, tempdir +from MDAnalysisTests.coordinates.base import ( + MultiframeReaderTest, BaseWriterTest, BaseReference +) +from MDAnalysisTests.coordinates.test_xyz import XYZReference + + +# skip entire test module if no appropriate chemfiles +chemfiles = pytest.importorskip('chemfiles') +@pytest.mark.skipif(not check_chemfiles_version(), + reason="Wrong version of chemfiles") +class TestChemFileXYZ(MultiframeReaderTest): + @staticmethod + @pytest.fixture + def ref(): + base = XYZReference() + base.writer = ChemfilesWriter + base.dimensions = np.array([0, 0, 0, 90, 90, 90], dtype=np.float32) + + return base + + @pytest.fixture + def reader(self, ref): + reader = ChemfilesReader(ref.trajectory) + reader.add_auxiliary('lowf', ref.aux_lowf, dt=ref.aux_lowf_dt, initial_time=0, time_selector=None) + reader.add_auxiliary('highf', ref.aux_highf, dt=ref.aux_highf_dt, initial_time=0, time_selector=None) + return reader + + + +class ChemfilesXYZReference(BaseReference): + def __init__(self): + super(ChemfilesXYZReference, self).__init__() + self.trajectory = datafiles.COORDINATES_XYZ + self.topology = datafiles.COORDINATES_XYZ + self.reader = ChemfilesReader + self.writer = ChemfilesWriter + self.ext = 'xyz' + self.volume = 0 + self.dimensions = np.zeros(6) + self.dimensions[3:] = 90.0 + + +@pytest.mark.skipif(not check_chemfiles_version(), + reason="Wrong version of chemfiles") +class TestChemfilesReader(MultiframeReaderTest): + @staticmethod + @pytest.fixture() + def ref(): + return ChemfilesXYZReference() + + +@pytest.mark.skipif(not check_chemfiles_version(), + reason="Wrong version of chemfiles") +class TestChemfilesWriter(BaseWriterTest): + @staticmethod + @pytest.fixture() + def ref(): + return ChemfilesXYZReference() + + # Disable 'test_no_container' as it try to open a file for writing without + # extension. + def test_no_container(self, ref): + pass + + def test_no_extension_raises(self, ref): + with pytest.raises(chemfiles.ChemfilesError): + ref.writer('foo') + + +@pytest.mark.skipif(not check_chemfiles_version(), + reason="Wrong version of chemfiles") +class TestChemfiles(object): + def test_read_chemfiles_format(self): + u = mda.Universe( + datafiles.LAMMPSdata, + format="chemfiles", + topology_format="data", + chemfiles_format="LAMMPS Data" + ) + + for ts in u.trajectory: + assert ts.n_atoms == 18364 + + def test_changing_system_size(self): + outfile = tempdir.TempDir().name + "chemfiles-changing-size.xyz" + with open(outfile, "w") as fd: + fd.write(VARYING_XYZ) + + u = mda.Universe(outfile, format="chemfiles", topology_format="XYZ") + + with pytest.raises(IOError): + u.trajectory._read_next_timestep() + + def test_wrong_open_mode(self): + with pytest.raises(IOError): + _ = ChemfilesWriter("", mode="r") + + def check_topology(self, reference, file): + u = mda.Universe(reference) + atoms = set([ + (atom.name, atom.type, atom.record_type) + for atom in u.atoms + ]) + bonds = set([ + (bond.atoms[0].ix, bond.atoms[1].ix) + for bond in u.bonds + ]) + + check = mda.Universe(file) + np.testing.assert_equal( + u.trajectory.ts.positions, + check.trajectory.ts.positions + ) + + for atom in check.atoms: + assert (atom.name, atom.type, atom.record_type) in atoms + + for bond in check.bonds: + assert (bond.atoms[0].ix, bond.atoms[1].ix) in bonds + + def test_write_topology(self): + u = mda.Universe(datafiles.CONECT) + outfile = tempdir.TempDir().name + "chemfiles-write-topology.pdb" + with ChemfilesWriter(outfile) as writer: + writer.write(u) + self.check_topology(datafiles.CONECT, outfile) + + # Manually setting the topology when creating the ChemfilesWriter + # (1) from an object + with ChemfilesWriter(outfile, topology=u) as writer: + writer.write_next_timestep(u.trajectory.ts) + self.check_topology(datafiles.CONECT, outfile) + + # (2) from a file + with ChemfilesWriter(outfile, topology=datafiles.CONECT) as writer: + writer.write_next_timestep(u.trajectory.ts) + # FIXME: this does not work, since chemfiles also insert the bonds + # which are implicit in PDB format (between standard residues), while + # MDAnalysis only read the explicit CONNECT records. + + # self.check_topology(datafiles.CONECT, outfile) + + def test_write_velocities(self): + ts = mda.coordinates.base.Timestep(4, velocities=True) + ts.dimensions = [20, 30, 41, 90, 90, 90] + + ts.positions = [ + [1, 1, 1], + [2, 2, 2], + [3, 3, 3], + [4, 4, 4], + ] + ts.velocities = [ + [10, 10, 10], + [20, 20, 20], + [30, 30, 30], + [40, 40, 40], + ] + + u = mda.Universe.empty(4) + u.add_TopologyAttr('type') + u.atoms[0].type = "H" + u.atoms[1].type = "H" + u.atoms[2].type = "H" + u.atoms[3].type = "H" + + outfile = tempdir.TempDir().name + "chemfiles-write-velocities.lmp" + with ChemfilesWriter(outfile, topology=u, chemfiles_format='LAMMPS Data') as writer: + writer.write_next_timestep(ts) + + with open(outfile) as file: + content = file.read() + assert content == EXPECTED_LAMMPS_DATA + + +VARYING_XYZ = """2 + +A 0 0 0 +A 0 0 0 +4 + +A 0 0 0 +A 0 0 0 +A 0 0 0 +A 0 0 0 +""" + + +EXPECTED_LAMMPS_DATA = """LAMMPS data file -- atom_style full -- generated by chemfiles +4 atoms +0 bonds +0 angles +0 dihedrals +0 impropers +1 atom types +0 bond types +0 angle types +0 dihedral types +0 improper types +0 20.0 xlo xhi +0 30.0 ylo yhi +0 41.0 zlo zhi + +# Pair Coeffs +# 1 H + +Masses + +1 0.0 # H + +Atoms # full + +1 1 1 0.0 1.0 1.0 1.0 # H +2 2 1 0.0 2.0 2.0 2.0 # H +3 3 1 0.0 3.0 3.0 3.0 # H +4 4 1 0.0 4.0 4.0 4.0 # H + +Velocities + +1 10.0 10.0 10.0 +2 20.0 20.0 20.0 +3 30.0 30.0 30.0 +4 40.0 40.0 40.0 +""" diff --git a/testsuite/MDAnalysisTests/coordinates/test_xyz.py b/testsuite/MDAnalysisTests/coordinates/test_xyz.py index 6436233fa1d..89e3ac975d9 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_xyz.py +++ b/testsuite/MDAnalysisTests/coordinates/test_xyz.py @@ -175,4 +175,4 @@ def test_elements_and_names(self, outfile, attr): with open(outfile, "r") as r: names = ''.join(l.split()[0].strip() for l in r.readlines()[2:-1]) - assert names[:-1].lower() == 'testing' \ No newline at end of file + assert names[:-1].lower() == 'testing'