mda: extend MDAnalysis plugin by bond information.

samples/MDAnalysisIntegration.py

@@ -36,6 +36,7 @@
 for i in range(10):
     system.part.add(id=i, pos=np.random.random(3) * system.box_l,
                     v=np.random.random(3))
+
 for i in range(5, 10):
     system.part[i].q = 1.0
     system.part[i].type = 1
@@ -75,13 +76,13 @@
 
 #
 # ========================================================="
-# Example #2: Write the configuration on a PDB file        "
+# Example #2: Write the configuration to a PDB file        "
 # ========================================================="
 #
 
 
 u.atoms.write("system.pdb")
-print("===> The initial configuration has been written on system.pdb ")
+print("===> The initial configuration has been written to system.pdb ")
 
 
 #

src/python/espressomd/MDA_ESP/__init__.py

@@ -24,21 +24,19 @@
 
 A minimal working example is the following:
 
->>> # imports
 >>> import espressomd
 >>> from espressomd import MDA_ESP
 >>> import MDAnalysis as mda
 >>> # system setup
->>> system = espressomd.System()
+>>> system = espressomd.System(box_l=[10., 10., 10.])
 >>> system.time_step = 1.
 >>> system.cell_system.skin = 1.
->>> system.box_l = [10.,10.,10.]
->>> system.part.add(id=0,pos=[1.,2.,3.])
+>>> system.part.add(id=0, pos=[1., 2., 3.])
 >>> # set up the stream
 >>> eos = MDA_ESP.Stream(system)
->>> # feed Universe with a topology and with coordinates
->>> u = mda.Universe(eos.topology,eos.trajectory)
->>> print u
+>>> # feed Universe with a topology and coordinates
+>>> u = mda.Universe(eos.topology, eos.trajectory)
+>>> print(u)
 <Universe with 1 atoms>
 
 """
@@ -66,7 +64,7 @@
 from MDAnalysis.core.topologyattrs import (
     Atomnames, Atomids, Atomtypes, Masses,
     Resids, Resnums, Segids, Resnames, AltLocs,
-    ICodes, Occupancies, Tempfactors, Charges
+    ICodes, Occupancies, Tempfactors, Charges, Bonds, Angles, Dihedrals
 )
 
 
@@ -76,7 +74,7 @@ class Stream:
     Create an object that provides a MDAnalysis topology and a coordinate reader
 
     >>> eos = MDA_ESP.Stream(system)
-    >>> u = mda.Universe(eos.topology,eos.trajectory)
+    >>> u = mda.Universe(eos.topology, eos.trajectory)
 
     Parameters
     ----------
@@ -144,12 +142,27 @@ def parse(self):
         atomtypes = []
         masses = []
         charges = []
+        bonds = []
+        angles = []
+        dihedrals = []
 
         for p in espresso.part:
             names.append("A" + repr(p.type))
             atomtypes.append("T" + repr(p.type))
             masses.append(p.mass)
             charges.append(p.q)
+            for bond in p.bonds:
+                partner_ids = bond[1:]
+                n_partner = len(partner_ids)
+                if n_partner == 1:
+                    bonds.append((p.id, partner_ids[0]))
+                elif n_partner == 2:
+                    angles.append((partner_ids[0], p.id, partner_ids[1]))
+                elif n_partner == 3:
+                    dihedrals.append(
+                        (partner_ids[0], p.id, partner_ids[1], partner_ids[2]))
+                else:
+                    continue
         natoms = len(espresso.part)
         attrs = [Atomnames(np.array(names, dtype=object)),
                  Atomids(np.arange(natoms) + 1),
@@ -164,6 +177,9 @@ def parse(self):
                  Tempfactors(np.zeros(natoms)),
                  ICodes(np.array([' '], dtype=object)),
                  Charges(np.array(charges)),
+                 Bonds(bonds),
+                 Angles(angles),
+                 Dihedrals(dihedrals)
                  ]
 
         top = Topology(natoms, 1, 1, attrs=attrs)

testsuite/python/CMakeLists.txt

@@ -208,6 +208,7 @@ python_test(FILE elc_vs_analytic.py MAX_NUM_PROC 2)
 python_test(FILE rotation.py MAX_NUM_PROC 1)
 python_test(FILE shapes.py MAX_NUM_PROC 1)
 python_test(FILE h5md.py MAX_NUM_PROC 2)
+python_test(FILE mdanalysis.py MAX_NUM_PROC 2)
 
 add_custom_target(
   python_test_data

testsuite/python/mdanalysis.py

@@ -0,0 +1,85 @@
+#
+# Copyright (C) 2019-2020 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+"""
+Testmodule for the MDAnalysis interface.
+"""
+import espressomd
+import espressomd.interactions
+import numpy as np
+import unittest as ut
+try:
+    import MDAnalysis as mda
+    from espressomd import MDA_ESP
+    skipIfMissingPythonPackage = ut.case._id
+except ImportError:
+    skipIfMissingPythonPackage = ut.skip(
+        "Python module MDAnalysis not available, skipping test!")
+
+
+@skipIfMissingPythonPackage
+class TestMDAnalysis(ut.TestCase):
+    system = espressomd.System(box_l=[10.0, 10.0, 10.0])
+    system.time_step = 0.001
+    system.cell_system.skin = 0.1
+
+    for i in range(10):
+        system.part.add(id=i, pos=[i, i % 2, 0], v=[0, i, -i], f=[1, 2 * i, 0],
+                        type=i % 2, q=i % 3 - 1)
+
+    bond = espressomd.interactions.HarmonicBond(k=1., r_0=1.5, r_cut=2.5)
+    angle = espressomd.interactions.AngleCosine(bend=1., phi0=2 * np.pi / 3)
+    dihe = espressomd.interactions.Dihedral(bend=1., mult=2, phase=np.pi / 3)
+    system.bonded_inter.add(bond)
+    system.bonded_inter.add(angle)
+    system.bonded_inter.add(dihe)
+    system.part[1].add_bond((bond, 0))
+    system.part[3].add_bond((angle, 2, 4))
+    system.part[6].add_bond((dihe, 5, 7, 8))
+
+    def test_universe(self):
+        system = self.system
+        eos = MDA_ESP.Stream(system)
+        u = mda.Universe(eos.topology, eos.trajectory)
+        # check atoms
+        self.assertEqual(len(u.atoms), 10)
+        np.testing.assert_equal(u.atoms.ids, np.arange(10) + 1)
+        np.testing.assert_equal(u.atoms.types, 5 * ['T0', 'T1'])
+        np.testing.assert_almost_equal(
+            u.atoms.charges, system.part[:].q, decimal=6)
+        np.testing.assert_almost_equal(
+            u.atoms.positions, system.part[:].pos, decimal=6)
+        np.testing.assert_almost_equal(
+            u.atoms.velocities, system.part[:].v, decimal=6)
+        np.testing.assert_almost_equal(
+            u.atoms.forces, system.part[:].f, decimal=6)
+        # check bonds
+        self.assertEqual(len(u.bonds), 1)
+        self.assertEqual(len(u.angles), 1)
+        self.assertEqual(len(u.dihedrals), 1)
+        np.testing.assert_equal(u.bonds[0].atoms.ids, [1, 2])
+        np.testing.assert_equal(u.angles[0].atoms.ids, [3, 4, 5])
+        np.testing.assert_equal(u.dihedrals[0].atoms.ids, [6, 7, 8, 9])
+        self.assertEqual(len(u.atoms[9].bonds), 0)
+        self.assertEqual(len(u.atoms[9].angles), 0)
+        self.assertEqual(len(u.atoms[9].dihedrals), 0)
+
+
+if __name__ == "__main__":
+    ut.main()

testsuite/scripts/samples/test_MDAnalysisIntegration.py

@@ -17,6 +17,7 @@
 
 import unittest as ut
 import importlib_wrapper
+import numpy as np
 
 try:
     import MDAnalysis  # pylint: disable=unused-import
@@ -33,6 +34,21 @@
 class Sample(ut.TestCase):
     system = sample.system
 
+    def test_universe(self):
+        system = self.system
+        u = sample.u
+        self.assertEqual(len(u.atoms), 10)
+        np.testing.assert_equal(u.atoms.ids, np.arange(10) + 1)
+        np.testing.assert_equal(u.atoms.types, sorted(5 * ['T0', 'T1']))
+        np.testing.assert_almost_equal(
+            u.atoms.charges, system.part[:].q, decimal=6)
+        np.testing.assert_almost_equal(
+            u.atoms.positions, system.part[:].pos, decimal=6)
+        np.testing.assert_almost_equal(
+            u.atoms.velocities, system.part[:].v, decimal=6)
+        np.testing.assert_almost_equal(
+            u.atoms.forces, system.part[:].f, decimal=6)
+
 
 if __name__ == "__main__":
     ut.main()