diff --git a/.disabled-travis.yml b/.disabled-travis.yml index feb2d9f8775..cde07f01320 100644 --- a/.disabled-travis.yml +++ b/.disabled-travis.yml @@ -30,7 +30,7 @@ env: - SETUP_CMD="${PYTEST_FLAGS}" - BUILD_CMD="pip install -e package/ && (cd testsuite/ && python setup.py build)" - CONDA_MIN_DEPENDENCIES="mmtf-python 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 chemfiles tqdm>=4.43.0 tidynamics>=1.0.0 rdkit>=2020.03.1 h5py" + - CONDA_DEPENDENCIES="${CONDA_MIN_DEPENDENCIES} seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12 chemfiles tqdm>=4.43.0 tidynamics>=1.0.0 rdkit>=2020.03.1 h5py openmm" - CONDA_CHANNELS='biobuilds conda-forge' - CONDA_CHANNEL_PRIORITY=True - PIP_DEPENDENCIES="duecredit parmed" @@ -43,7 +43,7 @@ env: matrix: # Run a coverage test for most versions - CODECOV="true" SETUP_CMD="${PYTEST_FLAGS} --cov=MDAnalysis" - - PYTHON_VERSION=3.8 CODECOV="true" SETUP_CMD="${PYTEST_FLAGS} --cov=MDAnalysis" + - PYTHON_VERSION=3.8 CODECOV="true" SETUP_CMD="${PYTEST_FLAGS} --cov=MDAnalysis" - PYTHON_VERSION=3.7 CODECOV="true" SETUP_CMD="${PYTEST_FLAGS} --cov=MDAnalysis" - NUMPY_VERSION=1.16.0 - NUMPY_VERSION=dev EVENT_TYPE="cron" diff --git a/.github/workflows/gh-ci.yaml b/.github/workflows/gh-ci.yaml index 70419bf5ebf..84c8c3947fc 100644 --- a/.github/workflows/gh-ci.yaml +++ b/.github/workflows/gh-ci.yaml @@ -14,7 +14,7 @@ defaults: shell: bash -l {0} env: - MDA_CONDA_MIN_DEPS: "pip pytest==6.1.2 mmtf-python biopython networkx cython matplotlib-base scipy griddataformats hypothesis gsd codecov threadpoolctl" + MDA_CONDA_MIN_DEPS: "pip pytest==6.1.2 mmtf-python biopython networkx cython matplotlib-base scipy griddataformats hypothesis gsd codecov threadpoolctl openmm" MDA_CONDA_EXTRA_DEPS: "seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12 chemfiles tqdm>=4.43.0 tidynamics>=1.0.0 rdkit>=2020.03.1 h5py" MDA_PIP_MIN_DEPS: 'coveralls coverage<5 pytest-cov pytest-xdist' MDA_PIP_EXTRA_DEPS: 'duecredit parmed' diff --git a/package/AUTHORS b/package/AUTHORS index b01c64ac42f..27e7a155623 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -165,6 +165,7 @@ Chronological list of authors - Henry Kobin - Kosuke Kudo - Sulay Shah + - Alexander Yang External code ------------- diff --git a/package/CHANGELOG b/package/CHANGELOG index f0f791f0da5..153319ddc7f 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -18,7 +18,7 @@ The rules for this file: calcraven，xiki-tempula, mieczyslaw, manuel.nuno.melo, PicoCentauri, hanatok, rmeli, aditya-kamath, tirkarthi, LeonardoBarneschi, hejamu, biogen98, orioncohen, z3y50n, hp115, ojeda-e, thadanipaarth, HenryKobin, - 1ut, sulays, PicoCentauri + 1ut, sulays, PicoCentauri, ahy3nz * 2.0.0 @@ -214,6 +214,7 @@ Changes ``analysis.helix_analysis`` (PR #2929) * Move water bridge analysis from hbonds to hydrogenbonds (Issue #2739 PR #2913) * `threadpoolctl` is required for installation (Issue #2860) + * Added OpenMM coordinate and topology converters (Issue #2863, PR #2917) Deprecations * In 2.1.0 the TRZReader will default to a dt of 1.0 ps when failing to diff --git a/package/MDAnalysis/coordinates/OpenMM.py b/package/MDAnalysis/coordinates/OpenMM.py new file mode 100644 index 00000000000..389576a7084 --- /dev/null +++ b/package/MDAnalysis/coordinates/OpenMM.py @@ -0,0 +1,198 @@ +# -*- 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 +# + +"""OpenMM structure I/O --- :mod:`MDAnalysis.coordinates.OpenMM` +================================================================ + + +Read coordinates data from a +`OpenMM `_ +:class:`simtk.openmm.app.simulation.Simulation` with :class:`OpenMMReader` +into a MDAnalysis Universe. + +Also converts other objects within the +`OpenMM Application Layer `_: + + - `simtk.openmm.app.pdbfile.PDBFile `_ + - `simtk.openmm.app.modeller.Modeller `_ + - `simtk.openmm.app.pdbxfile.PDBxFile `_ + +Example +------- +OpenMM can read various file formats into OpenMM objects. +MDAnalysis can then convert some of these OpenMM objects into MDAnalysis Universe objects. + + >>> import simtk.openmm.app as app + >>> import MDAnalysis as mda + >>> from MDAnalysis.tests.datafiles import PDBX + >>> pdbxfile = app.PDBxFile(PDBX) + >>> mda.Universe(pdbxfile) + + + + +Classes +------- + +.. autoclass:: OpenMMSimulationReader + :members: + +.. autoclass:: OpenMMAppReader + :members: + + +""" + + +import numpy as np + +from . import base +from .. import units + + +class OpenMMSimulationReader(base.SingleFrameReaderBase): + """Reader for OpenMM Simulation objects + + + .. versionadded:: 2.0.0 + """ + + format = "OPENMMSIMULATION" + units = {"time": "ps", "length": "nm", "velocity": "nm/ps", + "force": "kJ/(mol*nm)", "energy": "kJ/mol"} + + @staticmethod + def _format_hint(thing): + """Can this reader read *thing*? + """ + try: + from simtk.openmm.app import Simulation + except ImportError: + return False + else: + return isinstance(thing, Simulation) + + def _read_first_frame(self): + self.n_atoms = self.filename.topology.getNumAtoms() + + self.ts = self._mda_timestep_from_omm_context() + + if self.convert_units: + self.convert_pos_from_native(self.ts._pos) + self.ts.triclinic_dimensions = self.convert_pos_from_native( + self.ts.triclinic_dimensions, inplace=False + ) + self.ts.dimensions[3:] = _sanitize_box_angles(self.ts.dimensions[3:]) + self.convert_velocities_from_native(self.ts._velocities) + self.convert_forces_from_native(self.ts._forces) + self.convert_time_from_native(self.ts.dt) + + def _mda_timestep_from_omm_context(self): + """ Construct Timestep object from OpenMM context """ + import simtk.unit as u + + state = self.filename.context.getState(-1, getVelocities=True, + getForces=True, getEnergy=True) + + n_atoms = self.filename.context.getSystem().getNumParticles() + + ts = self._Timestep(n_atoms, **self._ts_kwargs) + ts.frame = 0 + ts.data["time"] = state.getTime()._value + ts.data["potential_energy"] = ( + state.getPotentialEnergy().in_units_of(u.kilojoule/u.mole) + ) + ts.data["kinetic_energy"] = ( + state.getKineticEnergy().in_units_of(u.kilojoule/u.mole) + ) + ts.triclinic_dimensions = state.getPeriodicBoxVectors( + asNumpy=True)._value + ts.dimensions[3:] = _sanitize_box_angles(ts.dimensions[3:]) + ts.positions = state.getPositions(asNumpy=True)._value + ts.velocities = state.getVelocities(asNumpy=True)._value + ts.forces = state.getForces(asNumpy=True)._value + + return ts + + +class OpenMMAppReader(base.SingleFrameReaderBase): + """Reader for OpenMM Application layer objects + + See also `the object definition in the OpenMM Application layer `_ + + .. versionadded:: 2.0.0 + """ + + format = "OPENMMAPP" + units = {"time": "ps", "length": "nm"} + + @staticmethod + def _format_hint(thing): + """Can this reader read *thing*? + """ + try: + from simtk.openmm import app + except ImportError: + return False + else: + return isinstance(thing, (app.PDBFile, app.Modeller, + app.PDBxFile)) + + def _read_first_frame(self): + self.n_atoms = self.filename.topology.getNumAtoms() + + self.ts = self._mda_timestep_from_omm_app() + + if self.convert_units: + self.convert_pos_from_native(self.ts._pos) + self.ts.triclinic_dimensions = self.convert_pos_from_native( + self.ts.triclinic_dimensions, inplace=False + ) + self.ts.dimensions[3:] = _sanitize_box_angles(self.ts.dimensions[3:]) + + def _mda_timestep_from_omm_app(self): + """ Construct Timestep object from OpenMM Application object """ + + omm_object = self.filename + n_atoms = omm_object.topology.getNumAtoms() + + ts = self._Timestep(n_atoms, **self._ts_kwargs) + ts.frame = 0 + if omm_object.topology.getPeriodicBoxVectors() is not None: + ts.triclinic_dimensions = np.array( + omm_object.topology.getPeriodicBoxVectors()._value + ) + ts.dimensions[3:] = _sanitize_box_angles(ts.dimensions[3:]) + ts.positions = np.array(omm_object.getPositions()._value) + + return ts + + +def _sanitize_box_angles(angles): + """ Ensure box angles correspond to first quadrant + + See `discussion on unitcell angles `_ + """ + inverted = 180 - angles + + return np.min(np.array([angles, inverted]), axis=0) diff --git a/package/MDAnalysis/coordinates/__init__.py b/package/MDAnalysis/coordinates/__init__.py index 6fb2ac4a99a..4855f5cdecb 100644 --- a/package/MDAnalysis/coordinates/__init__.py +++ b/package/MDAnalysis/coordinates/__init__.py @@ -732,6 +732,7 @@ class can choose an appropriate reader automatically. from . import INPCRD from . import LAMMPS from . import MOL2 +from . import OpenMM from . import ParmEd from . import PDB from . import PDBQT diff --git a/package/MDAnalysis/topology/OpenMMParser.py b/package/MDAnalysis/topology/OpenMMParser.py new file mode 100644 index 00000000000..01a6eb50183 --- /dev/null +++ b/package/MDAnalysis/topology/OpenMMParser.py @@ -0,0 +1,189 @@ +# -*- 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 +# + +"""OpenMM topology parser +========================= + +.. versionadded:: 2.0.0 + + +Converts an +`OpenMM `_ +:class:`simtk.openmm.app.topology.Topology` into a :class:`MDAnalysis.core.Topology`. + +Also converts some objects within the +`OpenMM Application layer `_ + + - `simtk.openmm.app.pdbfile.PDBFile `_ + - `simtk.openmm.app.simulation.Simulation `_ + - `simtk.openmm.app.modeller.Modeller `_ + - `simtk.openmm.app.pdbxfile.PDBxFile `_ + +The :class:`OpenMMTopologyParser` generates a topology from an OpenMM Topology object. + + +Classes +------- + +.. autoclass:: OpenMMTopologyParser + :members: + :inherited-members: + +.. autoclass:: OpenMMAppTopologyParser + :members: + :inherited-members: + +""" + +import numpy as np + +from .base import TopologyReaderBase +from .guessers import guess_types +from ..core.topology import Topology +from ..core.topologyattrs import ( + Atomids, + Atomnames, + Atomtypes, + Bonds, + ChainIDs, + Elements, + Masses, + Resids, + Resnums, + Resnames, + Segids, +) + + +class OpenMMTopologyParser(TopologyReaderBase): + format = "OPENMMTOPOLOGY" + + @staticmethod + def _format_hint(thing): + """Can this Parser read object *thing*? + + """ + try: + from simtk.openmm import app + except ImportError: + return False + else: + return isinstance(thing, app.Topology) + + def _mda_topology_from_omm_topology(self, omm_topology): + """ Construct mda topology from omm topology + + Can be used for any openmm object that contains a topology object + + Parameters + ---------- + omm_topology: simtk.openmm.Topology + + Returns + ------- + top : MDAnalysis.core.topology.Topology + + """ + + atom_resindex = [a.residue.index for a in omm_topology.atoms()] + residue_segindex = [r.chain.index for r in omm_topology.residues()] + atomids = [a.id for a in omm_topology.atoms()] + atomnames = [a.name for a in omm_topology.atoms()] + chainids = [a.residue.chain.id for a in omm_topology.atoms()] + elements = [a.element.symbol for a in omm_topology.atoms()] + atomtypes = [a.element.symbol for a in omm_topology.atoms()] + masses = [a.element.mass._value for a in omm_topology.atoms()] + resnames = [r.name for r in omm_topology.residues()] + resids = [r.index + 1 for r in omm_topology.residues()] + resnums = resids.copy() + segids = [c.index for c in omm_topology.chains()] + bonds = [(b.atom1.index, b.atom2.index) for b in omm_topology.bonds()] + bond_orders = [b.order for b in omm_topology.bonds()] + bond_types = [b.type for b in omm_topology.bonds()] + + n_atoms = len(atomids) + n_residues = len(resids) + n_segments = len(segids) + + attrs = [ + Atomids(np.array(atomids, dtype=np.int32)), + Atomnames(np.array(atomnames, dtype=object)), + Atomtypes(np.array(atomtypes, dtype=object)), + Bonds(bonds, types=bond_types, order=bond_orders, guessed=False), + ChainIDs(np.array(chainids, dtype=object)), + Elements(np.array(elements, dtype=object)), + Masses(np.array(masses, dtype=np.float32)), + Resids(resids), + Resnums(resnums), + Resnames(resnames), + Segids(segids), + ] + + top = Topology( + n_atoms, + n_residues, + n_segments, + attrs=attrs, + atom_resindex=atom_resindex, + residue_segindex=residue_segindex, + ) + + return top + + def parse(self, **kwargs): + omm_topology = self.filename + top = self._mda_topology_from_omm_topology(omm_topology) + + return top + + +class OpenMMAppTopologyParser(OpenMMTopologyParser): + format = "OPENMMAPP" + + @staticmethod + def _format_hint(thing): + """Can this Parser read object *thing*? + + """ + try: + from simtk.openmm import app + except ImportError: + return False + else: + return isinstance( + thing, + ( + app.PDBFile, app.Modeller, + app.Simulation, app.PDBxFile + ) + ) + + def parse(self, **kwargs): + try: + omm_topology = self.filename.getTopology() + except AttributeError: + omm_topology = self.filename.topology + top = self._mda_topology_from_omm_topology(omm_topology) + + return top + diff --git a/package/MDAnalysis/topology/__init__.py b/package/MDAnalysis/topology/__init__.py index 8d7628693e9..a8cd14a0bbd 100644 --- a/package/MDAnalysis/topology/__init__.py +++ b/package/MDAnalysis/topology/__init__.py @@ -307,8 +307,7 @@ 'CRDParser', 'TOPParser', 'PDBQTParser', 'TPRParser', 'LAMMPSParser', 'XYZParser', 'GMSParser', 'DLPolyParser', 'HoomdXMLParser','GSDParser', 'ITPParser', 'ParmEdParser', - 'RDKitParser'] - + 'RDKitParser', 'OpenMMParser'] from . import core from . import PSFParser from . import TOPParser @@ -331,6 +330,7 @@ from . import GSDParser from . import MinimalParser from . import ITPParser +from . import OpenMMParser from . import ParmEdParser from . import RDKitParser from . import FHIAIMSParser diff --git a/package/doc/sphinx/source/documentation_pages/coordinates/OpenMM.rst b/package/doc/sphinx/source/documentation_pages/coordinates/OpenMM.rst new file mode 100644 index 00000000000..ece8b877519 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/coordinates/OpenMM.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.coordinates.OpenMM diff --git a/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst b/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst index 3666d4928b3..d9fa0362666 100644 --- a/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/coordinates_modules.rst @@ -32,6 +32,7 @@ provide the format in the keyword argument *format* to coordinates/MMTF coordinates/MOL2 coordinates/NAMDBIN + coordinates/OpenMM coordinates/PDB coordinates/PDBQT coordinates/PQR diff --git a/package/doc/sphinx/source/documentation_pages/topology/OpenMM.rst b/package/doc/sphinx/source/documentation_pages/topology/OpenMM.rst new file mode 100644 index 00000000000..5c312eaea86 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/topology/OpenMM.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.topology.OpenMMParser diff --git a/package/doc/sphinx/source/documentation_pages/topology_modules.rst b/package/doc/sphinx/source/documentation_pages/topology_modules.rst index ed8caba8ce6..9b62fbf938a 100644 --- a/package/doc/sphinx/source/documentation_pages/topology_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/topology_modules.rst @@ -38,6 +38,7 @@ topology file format in the *topology_format* keyword argument to topology/MinimalParser topology/MMTFParser topology/MOL2Parser + topology/OpenMM topology/PDBParser topology/ExtendedPDBParser topology/PDBQTParser diff --git a/testsuite/MDAnalysisTests/coordinates/test_openmm.py b/testsuite/MDAnalysisTests/coordinates/test_openmm.py new file mode 100644 index 00000000000..7025658dc3c --- /dev/null +++ b/testsuite/MDAnalysisTests/coordinates/test_openmm.py @@ -0,0 +1,242 @@ +# -*- 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_almost_equal + +from MDAnalysisTests.datafiles import PDBX +from MDAnalysisTests.coordinates.reference import RefAdKSmall + +from MDAnalysisTests.coordinates.base import _SingleFrameReader + +import MDAnalysis as mda + +mm = pytest.importorskip("simtk.openmm") +unit = pytest.importorskip("simtk.unit") +app = pytest.importorskip("simtk.openmm.app") + + +class TestOpenMMBasicSimulationReader(): + @pytest.fixture + def omm_sim_uni(self): + system = mm.System() + topology = app.Topology() + chain = topology.addChain("CHAIN") + hydrogen = app.element.Element.getByAtomicNumber(1) + residue = topology.addResidue("RES", chain) + for i in range(5): + system.addParticle(1.0) + topology.addAtom(hydrogen.symbol, hydrogen, residue) + positions = np.ones((5, 3)) * unit.angstrom + integrator = mm.LangevinIntegrator( + 273 * unit.kelvin, + 1.0 / unit.picoseconds, 2.0 * unit.femtoseconds, + ) + simulation = app.Simulation(topology, system, integrator) + simulation.context.setPositions(positions) + + return mda.Universe(simulation) + + def test_dimensions(self, omm_sim_uni): + assert_almost_equal( + omm_sim_uni.trajectory.ts.dimensions, + np.array([20., 20., 20., 90., 90., 90.]), + 3, + "OpenMMBasicSimulationReader failed to get unitcell dimensions " + + "from OpenMM Simulation Object", + ) + + def test_coordinates(self, omm_sim_uni): + up = omm_sim_uni.atoms.positions + reference = np.ones((5, 3)) + assert_almost_equal(up, reference, decimal=3) + + def test_basic_topology(self, omm_sim_uni): + assert omm_sim_uni.atoms.n_atoms == 5 + assert omm_sim_uni.residues.n_residues == 1 + assert omm_sim_uni.residues.resnames[0] == "RES" + assert omm_sim_uni.segments.n_segments == 1 + assert omm_sim_uni.segments.segids[0] == 0 + assert len(omm_sim_uni.bonds.indices) == 0 + + +class TestOpenMMPDBFileReader(_SingleFrameReader): + __test__ = True + + def setUp(self): + self.universe = mda.Universe(app.PDBFile(RefAdKSmall.filename)) + self.ref = mda.Universe(RefAdKSmall.filename) + self.prec = 3 + + def test_dimensions(self): + assert_almost_equal( + self.universe.trajectory.ts.dimensions, + self.ref.trajectory.ts.dimensions, + self.prec, + "OpenMMPDBFileReader failed to get unitcell dimensions " + + "from OpenMMPDBFile", + ) + + def test_coordinates(self): + up = self.universe.atoms.positions + rp = self.ref.atoms.positions + assert_almost_equal(up, rp, decimal=3) + + +class TestOpenMMModellerReader(_SingleFrameReader): + __test__ = True + + def setUp(self): + pdb_obj = app.PDBFile(RefAdKSmall.filename) + modeller = app.Modeller(pdb_obj.topology, pdb_obj.positions) + self.universe = mda.Universe(modeller) + self.ref = mda.Universe(RefAdKSmall.filename) + self.prec = 3 + + def test_dimensions(self): + assert_almost_equal( + self.universe.trajectory.ts.dimensions, + self.ref.trajectory.ts.dimensions, + self.prec, + "OpenMMModellerReader failed to get unitcell dimensions " + + "from OpenMMModeller", + ) + + def test_coordinates(self): + up = self.universe.atoms.positions + rp = self.ref.atoms.positions + assert_almost_equal(up, rp, decimal=3) + + +class TestOpenMMSimulationReader(_SingleFrameReader): + __test__ = True + + def setUp(self): + pdb = app.PDBFile(RefAdKSmall.filename) + forcefield = app.ForceField("amber99sbildn.xml") + system = forcefield.createSystem( + pdb.topology, nonbondedMethod=app.NoCutoff, constraints=app.HBonds + ) + integrator = mm.LangevinIntegrator( + 300 * unit.kelvin, 1 / unit.picoseconds, 2.0 * unit.femtoseconds + ) + sim = app.Simulation(pdb.topology, system, integrator) + sim.context.setPositions(pdb.positions) + self.universe = mda.Universe(sim) + self.ref = mda.Universe(RefAdKSmall.filename) + self.prec = 3 + + def test_dimensions(self): + assert_almost_equal( + self.universe.trajectory.ts.dimensions, + self.ref.trajectory.ts.dimensions, + self.prec, + "OpenMMSimulationReader failed to get unitcell dimensions " + + "from OpenMMSimulation", + ) + + def test_coordinates(self): + up = self.universe.atoms.positions + rp = self.ref.atoms.positions + assert_almost_equal(up, rp, decimal=3) + + @pytest.mark.xfail(reason='OpenMM pickling not supported yet') + def test_pickle_singleframe_reader(self): + """ + See `OpenMM SwigPyObject serialisation discussion `_ + """ + super().test_pickle_singleframe_reader() + + +@pytest.fixture +def PDBX_U(): + return mda.Universe(app.PDBxFile(PDBX)) + + +def test_pdbx_coordinates(PDBX_U): + ref_pos = 10 * np.array( + [ + [0.6537, 3.7168, 1.4751], + [0.6141, 3.7813, 1.3496], + [0.4851, 3.8633, 1.3566], + [0.4428, 3.914, 1.2532], + [0.6074, 3.6749, 1.2371], + [0.7366, 3.5948, 1.2334], + [0.8485, 3.6434, 1.1667], + [0.7494, 3.477, 1.3064], + [0.9699, 3.5741, 1.1709], + [0.8706, 3.4079, 1.3107], + [0.9803, 3.4575, 1.2437], + [0.4265, 3.8799, 1.476], + [0.3019, 3.9561, 1.4964], + [0.3076, 4.1011, 1.4488], + [0.2061, 4.1541, 1.4068], + [0.2633, 3.9555, 1.6453], + [0.1999, 3.8275, 1.6944], + [0.1911, 3.7306, 1.6147], + [0.1619, 3.8224, 1.8133], + [0.424, 4.1656, 1.4575], + [0.4388, 4.3071, 1.4222], + [0.5294, 4.3324, 1.3023], + [0.5726, 4.4457, 1.2828], + [0.4853, 4.387, 1.5449], + [0.3731, 4.4096, 1.6449], + [0.4257, 4.4404, 1.7833], + [0.3124, 4.4673, 1.8792], + [0.5552, 4.2303, 1.2196], + [0.6382, 4.2504, 1.1005], + [0.555, 4.3159, 0.9885], + [0.4442, 4.2696, 0.9581], + [0.7097, 4.119, 1.0495], + [0.7974, 4.0487, 1.1586], + [0.7916, 4.1436, 0.9212], + [0.9078, 4.133, 1.2278], + [0.6112, 4.4208, 0.9256], + [0.5533, 4.4859, 0.8069], + [0.6548, 4.4763, 0.6937], + [0.6177, 4.4543, 0.5787], + [0.5139, 4.6313, 0.8343], + [0.3789, 4.6433, 0.904], + [0.7844, 4.4865, 0.7289], + [0.8992, 4.4777, 0.6384], + [0.9414, 4.3301, 0.6208], + [1.0515, 4.2927, 0.6608], + [1.0151, 4.5643, 0.6923], + [0.8528, 4.2473, 0.5599], + [0.8763, 4.1044, 0.5361], + [1.0028, 4.0746, 0.4537], + [1.0652, 3.9719, 0.478], + [0.7534, 4.0355, 0.4762], + [0.6391, 4.0245, 0.5727], + [0.5337, 4.1106, 0.5851], + [0.6204, 3.9232, 0.6731], + [0.4486, 4.0675, 0.6847], + [0.5002, 3.9534, 0.7413], + [0.695, 3.8109, 0.714], + [0.4534, 3.8758, 0.8486], + [0.6468, 3.7325, 0.8175], + [0.5277, 3.7648, 0.8837], + ] + ) + rp = PDBX_U.atoms.positions + assert_almost_equal(ref_pos, rp, decimal=3) diff --git a/testsuite/MDAnalysisTests/data/4x8u.pdbx b/testsuite/MDAnalysisTests/data/4x8u.pdbx new file mode 100644 index 00000000000..5ef564d4929 --- /dev/null +++ b/testsuite/MDAnalysisTests/data/4x8u.pdbx @@ -0,0 +1,88 @@ +# extract from 4X8U involving insertion codes +# +data_4X8U +# +_entry.id 4X8U +# +loop_ +_atom_site.group_PDB +_atom_site.id +_atom_site.type_symbol +_atom_site.label_atom_id +_atom_site.label_alt_id +_atom_site.label_comp_id +_atom_site.label_asym_id +_atom_site.label_entity_id +_atom_site.label_seq_id +_atom_site.pdbx_PDB_ins_code +_atom_site.Cartn_x +_atom_site.Cartn_y +_atom_site.Cartn_z +_atom_site.occupancy +_atom_site.B_iso_or_equiv +_atom_site.pdbx_formal_charge +_atom_site.auth_seq_id +_atom_site.auth_comp_id +_atom_site.auth_asym_id +_atom_site.auth_atom_id +_atom_site.pdbx_PDB_model_num +ATOM 315 N N . PHE A 1 43 ? 6.537 37.168 14.751 1.00 16.93 ? 59 PHE H N 1 +ATOM 316 C CA . PHE A 1 43 ? 6.141 37.813 13.496 1.00 16.64 ? 59 PHE H CA 1 +ATOM 317 C C . PHE A 1 43 ? 4.851 38.633 13.566 1.00 21.61 ? 59 PHE H C 1 +ATOM 318 O O . PHE A 1 43 ? 4.428 39.140 12.532 1.00 22.06 ? 59 PHE H O 1 +ATOM 319 C CB . PHE A 1 43 ? 6.074 36.749 12.371 1.00 17.35 ? 59 PHE H CB 1 +ATOM 320 C CG . PHE A 1 43 ? 7.366 35.948 12.334 1.00 17.99 ? 59 PHE H CG 1 +ATOM 321 C CD1 . PHE A 1 43 ? 8.485 36.434 11.667 1.00 21.29 ? 59 PHE H CD1 1 +ATOM 322 C CD2 . PHE A 1 43 ? 7.494 34.770 13.064 1.00 18.58 ? 59 PHE H CD2 1 +ATOM 323 C CE1 . PHE A 1 43 ? 9.699 35.741 11.709 1.00 22.04 ? 59 PHE H CE1 1 +ATOM 324 C CE2 . PHE A 1 43 ? 8.706 34.079 13.107 1.00 20.43 ? 59 PHE H CE2 1 +ATOM 325 C CZ . PHE A 1 43 ? 9.803 34.575 12.437 1.00 18.31 ? 59 PHE H CZ 1 +ATOM 326 N N . ASP A 1 44 ? 4.265 38.799 14.760 1.00 19.41 ? 60 ASP H N 1 +ATOM 327 C CA . ASP A 1 44 ? 3.019 39.561 14.964 1.00 21.98 ? 60 ASP H CA 1 +ATOM 328 C C . ASP A 1 44 ? 3.076 41.011 14.488 1.00 31.37 ? 60 ASP H C 1 +ATOM 329 O O . ASP A 1 44 ? 2.061 41.541 14.068 1.00 32.54 ? 60 ASP H O 1 +ATOM 330 C CB . ASP A 1 44 ? 2.633 39.555 16.453 1.00 23.09 ? 60 ASP H CB 1 +ATOM 331 C CG . ASP A 1 44 ? 1.999 38.275 16.944 1.00 28.33 ? 60 ASP H CG 1 +ATOM 332 O OD1 . ASP A 1 44 ? 1.911 37.306 16.147 1.00 29.31 ? 60 ASP H OD1 1 +ATOM 333 O OD2 . ASP A 1 44 ? 1.619 38.224 18.133 1.00 29.53 ? 60 ASP H OD2 1 +ATOM 334 N N . LYS A 1 45 A 4.240 41.656 14.575 1.00 30.66 ? 60 LYS H N 1 +ATOM 335 C CA . LYS A 1 45 A 4.388 43.071 14.222 1.00 31.68 ? 60 LYS H CA 1 +ATOM 336 C C . LYS A 1 45 A 5.294 43.324 13.023 1.00 38.11 ? 60 LYS H C 1 +ATOM 337 O O . LYS A 1 45 A 5.726 44.457 12.828 1.00 39.40 ? 60 LYS H O 1 +ATOM 338 C CB . LYS A 1 45 A 4.853 43.870 15.449 1.00 33.51 ? 60 LYS H CB 1 +ATOM 339 C CG . LYS A 1 45 A 3.731 44.096 16.449 1.00 48.38 ? 60 LYS H CG 1 +ATOM 340 C CD . LYS A 1 45 A 4.257 44.404 17.833 1.00 64.54 ? 60 LYS H CD 1 +ATOM 341 C CE . LYS A 1 45 A 3.124 44.673 18.792 1.00 76.49 ? 60 LYS H CE 1 +ATOM 342 N N . ILE A 1 46 B 5.552 42.303 12.196 1.00 34.95 ? 60 ILE H N 1 +ATOM 343 C CA . ILE A 1 46 B 6.382 42.504 11.005 1.00 35.67 ? 60 ILE H CA 1 +ATOM 344 C C . ILE A 1 46 B 5.550 43.159 9.885 1.00 40.88 ? 60 ILE H C 1 +ATOM 345 O O . ILE A 1 46 B 4.442 42.696 9.581 1.00 39.95 ? 60 ILE H O 1 +ATOM 346 C CB . ILE A 1 46 B 7.097 41.190 10.495 1.00 38.06 ? 60 ILE H CB 1 +ATOM 347 C CG1 . ILE A 1 46 B 7.974 40.487 11.586 1.00 37.99 ? 60 ILE H CG1 1 +ATOM 348 C CG2 . ILE A 1 46 B 7.916 41.436 9.212 1.00 37.97 ? 60 ILE H CG2 1 +ATOM 349 C CD1 . ILE A 1 46 B 9.078 41.330 12.278 1.00 43.63 ? 60 ILE H CD1 1 +ATOM 350 N N . LYS A 1 47 C 6.112 44.208 9.256 1.00 38.12 ? 60 LYS H N 1 +ATOM 351 C CA . LYS A 1 47 C 5.533 44.859 8.069 1.00 39.21 ? 60 LYS H CA 1 +ATOM 352 C C . LYS A 1 47 C 6.548 44.763 6.937 1.00 45.85 ? 60 LYS H C 1 +ATOM 353 O O . LYS A 1 47 C 6.177 44.543 5.787 1.00 46.49 ? 60 LYS H O 1 +ATOM 354 C CB . LYS A 1 47 C 5.139 46.313 8.343 1.00 41.67 ? 60 LYS H CB 1 +ATOM 355 C CG . LYS A 1 47 C 3.789 46.433 9.040 1.00 58.21 ? 60 LYS H CG 1 +ATOM 356 N N . ASN A 1 48 D 7.844 44.865 7.289 1.00 43.92 ? 60 ASN H N 1 +ATOM 357 C CA . ASN A 1 48 D 8.992 44.777 6.384 1.00 43.65 ? 60 ASN H CA 1 +ATOM 358 C C . ASN A 1 48 D 9.414 43.301 6.208 1.00 44.11 ? 60 ASN H C 1 +ATOM 359 O O . ASN A 1 48 D 10.515 42.927 6.608 1.00 42.37 ? 60 ASN H O 1 +ATOM 360 C CB . ASN A 1 48 D 10.151 45.643 6.923 1.00 45.00 ? 60 ASN H CB 1 +ATOM 361 N N . TRP A 1 49 ? 8.528 42.473 5.599 1.00 39.73 ? 61 TRP H N 1 +ATOM 362 C CA . TRP A 1 49 ? 8.763 41.044 5.361 1.00 39.08 ? 61 TRP H CA 1 +ATOM 363 C C . TRP A 1 49 ? 10.028 40.746 4.537 1.00 42.50 ? 61 TRP H C 1 +ATOM 364 O O . TRP A 1 49 ? 10.652 39.719 4.780 1.00 44.39 ? 61 TRP H O 1 +ATOM 365 C CB . TRP A 1 49 ? 7.534 40.355 4.762 1.00 37.55 ? 61 TRP H CB 1 +ATOM 366 C CG . TRP A 1 49 ? 6.391 40.245 5.727 1.00 38.53 ? 61 TRP H CG 1 +ATOM 367 C CD1 . TRP A 1 49 ? 5.337 41.106 5.851 1.00 41.76 ? 61 TRP H CD1 1 +ATOM 368 C CD2 . TRP A 1 49 ? 6.204 39.232 6.731 1.00 37.95 ? 61 TRP H CD2 1 +ATOM 369 N NE1 . TRP A 1 49 ? 4.486 40.675 6.847 1.00 41.50 ? 61 TRP H NE1 1 +ATOM 370 C CE2 . TRP A 1 49 ? 5.002 39.534 7.413 1.00 42.54 ? 61 TRP H CE2 1 +ATOM 371 C CE3 . TRP A 1 49 ? 6.950 38.109 7.140 1.00 38.20 ? 61 TRP H CE3 1 +ATOM 372 C CZ2 . TRP A 1 49 ? 4.534 38.758 8.486 1.00 41.36 ? 61 TRP H CZ2 1 +ATOM 373 C CZ3 . TRP A 1 49 ? 6.468 37.325 8.175 1.00 38.85 ? 61 TRP H CZ3 1 +ATOM 374 C CH2 . TRP A 1 49 ? 5.277 37.648 8.837 1.00 39.68 ? 61 TRP H CH2 1 diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index d647f435169..879c9ee664d 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -201,6 +201,7 @@ "MMTF_NOCRYST", # File with meaningless CRYST1 record (Issue #2679, PR #2685) "FHIAIMS", # to test FHIAIMS coordinate files "SDF_molecule", # MDL SDFile for rdkit + "PDBX", # PDBxfile "PDB_elements", # PDB file with elements ] @@ -559,5 +560,7 @@ PDB_elements = resource_filename(__name__, 'data/elements.pdb') +PDBX = resource_filename(__name__, "data/4x8u.pdbx") + # This should be the last line: clean up namespace del resource_filename diff --git a/testsuite/MDAnalysisTests/topology/test_openmm.py b/testsuite/MDAnalysisTests/topology/test_openmm.py new file mode 100644 index 00000000000..d4208a0c4d4 --- /dev/null +++ b/testsuite/MDAnalysisTests/topology/test_openmm.py @@ -0,0 +1,140 @@ +# -*- 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 + +import MDAnalysis as mda + +from MDAnalysisTests.topology.base import ParserBase +from MDAnalysisTests.datafiles import CONECT, PDBX + + +app = pytest.importorskip('simtk.openmm.app') + + +class OpenMMTopologyBase(ParserBase): + parser = mda.topology.OpenMMParser.OpenMMTopologyParser + expected_attrs = [ + "ids", + "names", + "resids", + "resnames", + "masses", + "bonds", + "chainIDs", + "elements", + ] + expected_n_bonds = 0 + + def test_creates_universe(self, filename): + """Check that Universe works with this Parser""" + u = mda.Universe(filename, topology_format="OPENMMTOPOLOGY") + assert isinstance(u, mda.Universe) + + def test_attr_size(self, top): + assert len(top.ids) == top.n_atoms + assert len(top.names) == top.n_atoms + assert len(top.resids) == top.n_residues + assert len(top.resnames) == top.n_residues + + def test_atoms(self, top): + assert top.n_atoms == self.expected_n_atoms + + def test_bonds(self, top): + assert len(top.bonds.values) == self.expected_n_bonds + if self.expected_n_bonds: + assert isinstance(top.bonds.values[0], tuple) + else: + assert top.bonds.values == [] + + def test_resids(self, top): + assert len(top.resids.values) == self.expected_n_residues + if self.expected_n_residues: + assert isinstance(top.resids.values, np.ndarray) + else: + assert top.resids.values == [] + + def test_resnames(self, top): + assert len(top.resnames.values) == self.expected_n_residues + if self.expected_n_residues: + assert isinstance(top.resnames.values, np.ndarray) + else: + assert top.resnames.values == [] + + def test_resnums(self, top): + assert len(top.resnums.values) == self.expected_n_residues + if self.expected_n_residues: + assert isinstance(top.resnums.values, np.ndarray) + else: + assert top.resnums.values == [] + + def test_segids(self, top): + assert len(top.segids.values) == self.expected_n_segments + if self.expected_n_segments: + assert isinstance(top.segids.values, np.ndarray) + else: + assert top.segids.values == [] + + +class OpenMMAppTopologyBase(OpenMMTopologyBase): + parser = mda.topology.OpenMMParser.OpenMMAppTopologyParser + expected_attrs = [ + "ids", + "names", + "resids", + "resnames", + "masses", + "bonds", + "chainIDs", + "elements", + ] + expected_n_bonds = 0 + + def test_creates_universe(self, filename): + """Check that Universe works with this Parser""" + u = mda.Universe(filename, topology_format="OPENMMAPP") + assert isinstance(u, mda.Universe) + + +class TestOpenMMTopologyParser(OpenMMTopologyBase): + ref_filename = app.PDBFile(CONECT).topology + expected_n_atoms = 1890 + expected_n_residues = 199 + expected_n_segments = 3 + expected_n_bonds = 1922 + + +class TestOpenMMPDBFileParser(OpenMMAppTopologyBase): + ref_filename = app.PDBFile(CONECT) + expected_n_atoms = 1890 + expected_n_residues = 199 + expected_n_segments = 3 + expected_n_bonds = 1922 + + +class TestOpenMMPDBxFileParser(OpenMMAppTopologyBase): + ref_filename = app.PDBxFile(PDBX) + expected_n_atoms = 60 + expected_n_residues = 7 + expected_n_segments = 1 + expected_n_bonds = 62