Added large membrane protien test system.

README.md

@@ -20,4 +20,14 @@ It has not yet been tested for:
 
 There has not yet been any effort to support:
 - Martini v3
-- Polarizable Martini
+- Polarizable Martini
+
+Limitations:
+- Math in .itp files:
+    - Gromacs allows for mathmatical expressions to be used in .itp files
+        - `    1     6      1   0.98112 RUBBER_FC*1.000000`
+    - This is not allowed in OpenMM and must be edited
+        - `    1     6      1   0.98112 RUBBER_FC`
+- The only supported electrostatic option is `coulombtype = reaction-field`
+- The only supported option for Van der Waals is `vdw_type = cutoff`
+- The cutoff must be the same for `rcoulomb` and `rvdw`

martini_openmm/martini.py

@@ -276,7 +276,21 @@ def _processLine(self, line, file):
                         "Found %s before [ moleculetype ]" % self._currentCategory
                     )
                 raise ValueError(
-                    f"[ {self._currentCategory} ] is currently not supported."
+                    f"[ {self._currentCategory} ] is currently not supported.\n"
+                    "Please file an issue at github.com/maccallumlab/martini_openmm"
+                )
+            elif self._currentCategory == "system":  # ignore the system description
+                pass
+            elif self._currentCategory == "position_restraints":
+                raise ValueError(
+                    "[ position_restraints ] is not currently supported.\n"
+                    "The same effect can be achieved using a CustomExternalForce in OpenMM.\n"
+                    "Please see github.com/maccallumlab/martini_openmm for details."
+                )
+            else:
+                raise ValueError(
+                    f"[ {self._currentCategory} ] is not currently supported.\n"
+                    "Please file an issue at github.com/maccallumlab/martini_openmm"
                 )
 
     def _processDefaults(self, line):
@@ -522,7 +536,6 @@ def _addAtomsToSystem(self, sys, moleculeType):
             sys.addParticle(mass)
 
     def _addBondsToSystem(self, sys, moleculeType, bondedTypes, baseAtomIndex):
-        bonds = None
         for fields in moleculeType.bonds:
             atoms = [int(x) - 1 for x in fields[:2]]
             types = tuple(bondedTypes[i] for i in atoms)
@@ -539,19 +552,17 @@ def _addBondsToSystem(self, sys, moleculeType, bondedTypes, baseAtomIndex):
 
             length = float(params[0])
 
-            if bonds is None:
-                bonds = mm.HarmonicBondForce()
-                sys.addForce(bonds)
-            bonds.addBond(
+            if self.bonds is None:
+                self.bonds = mm.HarmonicBondForce()
+                sys.addForce(self.bonds)
+            self.bonds.addBond(
                 baseAtomIndex + atoms[0],
                 baseAtomIndex + atoms[1],
                 length,
                 float(params[1]),
             )
 
     def _addAnglesToSystem(self, sys, moleculeType, bondedTypes, baseAtomIndex):
-        angles_harmonic = None
-        angles_gmx = None
         degToRad = math.pi / 180
 
         for fields in moleculeType.angles:
@@ -577,23 +588,25 @@ def _addAnglesToSystem(self, sys, moleculeType, bondedTypes, baseAtomIndex):
             # Need to implment gmx angle types
             theta = float(params[0]) * degToRad
             if int(fields[3]) == 2:
-                if angles_gmx is None:
-                    angles_gmx = mm.CustomAngleForce("0.5*k*(cos(theta)-cos(theta0))^2")
-                    angles_gmx.addPerAngleParameter("theta0")
-                    angles_gmx.addPerAngleParameter("k")
-                    sys.addForce(angles_gmx)
-                angles_gmx.addAngle(
+                if self.angles_gmx is None:
+                    self.angles_gmx = mm.CustomAngleForce(
+                        "0.5*k*(cos(theta)-cos(theta0))^2"
+                    )
+                    self.angles_gmx.addPerAngleParameter("theta0")
+                    self.angles_gmx.addPerAngleParameter("k")
+                    sys.addForce(self.angles_gmx)
+                self.angles_gmx.addAngle(
                     baseAtomIndex + atoms[0],
                     baseAtomIndex + atoms[1],
                     baseAtomIndex + atoms[2],
                     [theta, float(params[1])],
                 )
 
             elif int(fields[3]) == 1:
-                if angles_harmonic is None:
-                    angles_harmonic = mm.HarmonicAngleForce()
-                    sys.addForce(angles_harmonic)
-                angles_harmonic.addAngle(
+                if self.angles_harmonic is None:
+                    self.angles_harmonic = mm.HarmonicAngleForce()
+                    sys.addForce(self.angles_harmonic)
+                self.angles_harmonic.addAngle(
                     baseAtomIndex + atoms[0],
                     baseAtomIndex + atoms[1],
                     baseAtomIndex + atoms[2],
@@ -613,9 +626,6 @@ def _addTorsionToSystem(
         wildcardDihedralTypes,
         baseAtomIndex,
     ):
-        periodic = None
-        harmonicTorsion = None
-        rb = None
         degToRad = math.pi / 180
 
         for fields in moleculeType.dihedrals:
@@ -660,10 +670,10 @@ def _addTorsionToSystem(
                     # Periodic torsion
                     k = float(params[6])
                     if k != 0:
-                        if periodic is None:
-                            periodic = mm.PeriodicTorsionForce()
-                            sys.addForce(periodic)
-                        periodic.addTorsion(
+                        if self.periodic is None:
+                            self.periodic = mm.PeriodicTorsionForce()
+                            sys.addForce(self.periodic)
+                        self.periodic.addTorsion(
                             baseAtomIndex + atoms[0],
                             baseAtomIndex + atoms[1],
                             baseAtomIndex + atoms[2],
@@ -677,17 +687,17 @@ def _addTorsionToSystem(
                     k = float(params[6])
                     phi0 = float(params[5])
                     if k != 0:
-                        if harmonicTorsion is None:
-                            harmonicTorsion = mm.CustomTorsionForce(
+                        if self.harmonicTorsion is None:
+                            self.harmonicTorsion = mm.CustomTorsionForce(
                                 "0.5*k*(thetap-theta0)^2; thetap = step(-(theta-theta0+pi))*2*pi+theta+step(theta-theta0-pi)*(-2*pi); pi = %.15g"
                                 % math.pi
                             )
-                            harmonicTorsion.addPerTorsionParameter("theta0")
-                            harmonicTorsion.addPerTorsionParameter("k")
-                            sys.addForce(harmonicTorsion)
+                            self.harmonicTorsion.addPerTorsionParameter("theta0")
+                            self.harmonicTorsion.addPerTorsionParameter("k")
+                            sys.addForce(self.harmonicTorsion)
                         # map phi0 into correct space
                         phi0 = phi0 - 360 if phi0 > 180 else phi0
-                        harmonicTorsion.addTorsion(
+                        self.harmonicTorsion.addTorsion(
                             baseAtomIndex + atoms[0],
                             baseAtomIndex + atoms[1],
                             baseAtomIndex + atoms[2],
@@ -698,9 +708,9 @@ def _addTorsionToSystem(
                     # RB Torsion
                     c = [float(x) for x in params[5:11]]
                     if any(x != 0 for x in c):
-                        if rb is None:
-                            rb = mm.RBTorsionForce()
-                            sys.addForce(rb)
+                        if self.rb is None:
+                            self.rb = mm.RBTorsionForce()
+                            sys.addForce(self.rb)
                         if dihedralType == "5":
                             # Convert Fourier coefficients to RB coefficients.
                             c = [
@@ -711,7 +721,7 @@ def _addTorsionToSystem(
                                 -4 * c[3],
                                 0,
                             ]
-                        rb.addTorsion(
+                        self.rb.addTorsion(
                             baseAtomIndex + atoms[0],
                             baseAtomIndex + atoms[1],
                             baseAtomIndex + atoms[2],
@@ -962,7 +972,7 @@ def __init__(
         self._defines = OrderedDict()
         self._genpairs = True
         if defines is not None:
-            for define, value in defines.iteritems():
+            for define, value in defines.items():
                 self._defines[define] = value
 
         # Parse the file.
@@ -982,6 +992,12 @@ def __init__(
         self._nonbondTypes = {}
         self._processFile(file)
         self._all_vsites = []
+        self.angles_harmonic = None
+        self.angles_gmx = None
+        self.bonds = None
+        self.periodic = None
+        self.harmonicTorsion = None
+        self.rb = None
 
         top = Topology()
         self.topology = top

tests/Scripts/compare.py

@@ -5,6 +5,13 @@
 import sys
 
 
+# tolerances for energy and force comparisons
+energy_atol = 1e-4
+energy_rtol = 1e-5
+f_atol = 1e-4
+f_rtol = 1e-5
+
+
 def get_gmx_energy():
     f = open("gmx/energy.xvg")
     line = f.readlines()[-1]
@@ -39,42 +46,47 @@ def get_omm_forces():
 omm_energy = get_omm_energy()
 gmx_forces = get_gmx_forces()
 omm_forces = get_omm_forces()
-failed = False
 
-if not math.isclose(gmx_energy, omm_energy, rel_tol=1e-6):
+failed = False
+relative_energy_error = abs(gmx_energy - omm_energy) / abs(gmx_energy)
+abs_energy_error = abs(gmx_energy - omm_energy)
+if relative_energy_error > 1e-5 and abs_energy_error > 1e-3:
     print("Gromacs and OpenMM energies do not match!")
-    print(f"    Gromacs: {gmx_energy}")
-    print(f"     OpenMM: {omm_energy}")
-    print(f"      Delta: {gmx_energy - omm_energy}")
+    print(f"    Gromacs: {gmx_energy:16.3f}")
+    print(f"     OpenMM: {omm_energy:16.3f}")
+    print(f"      Delta: {(gmx_energy - omm_energy):16.3f}")
+    print(f"   Relative: {relative_energy_error:16.3f}")
     print()
     failed = True
 
-if not np.allclose(gmx_forces, omm_forces, atol=2e-2):
-    diffs = np.where(np.abs(gmx_forces - omm_forces) > 1e-2)
-    n_diff = len(diffs[0]) // 3
-    i = diffs[0][0]
-    g_f = gmx_forces[i, :]
-    o_f = omm_forces[i, :]
-    max_diff = np.max(np.abs(gmx_forces - omm_forces))
-    max_diff_ind = np.unravel_index(
-        np.argmax(np.abs(gmx_forces - omm_forces)), gmx_forces.shape
+if not np.allclose(omm_forces, gmx_forces, rtol=f_rtol, atol=f_atol):
+    errors = np.logical_not(
+        np.isclose(omm_forces, gmx_forces, rtol=f_rtol, atol=f_atol)
     )
-    ind = max_diff_ind[0]
+    n_diff = np.sum(errors)
+    bad_ind = np.where(errors)
+    first_row = bad_ind[0][0]
+    first_col = bad_ind[1][0]
+    g_f = gmx_forces[first_row, :]
+    o_f = omm_forces[first_row, :]
+    d_f = np.abs(g_f - o_f)
+    errors = np.abs(omm_forces - gmx_forces) - f_atol - f_rtol * np.abs(gmx_forces)
+    max_ind = np.unravel_index(np.argmax(errors, axis=None), errors.shape)
+    max_g = gmx_forces[max_ind[0], :]
+    max_o = omm_forces[max_ind[0], :]
+    max_d = np.abs(max_g - max_o)
     print("Gromacs and OpenMM forces do not match!")
     print(f"    Values differ at {n_diff} of {gmx_forces.shape[0]} positions.")
     print()
-    print(f"    First differing position: {i}")
+    print(f"    First differing position: ({first_row}, {first_col})")
     print(f"    Gromacs: {g_f[0]:20f} {g_f[1]:20f} {g_f[2]:20f}")
     print(f"     OpenMM: {o_f[0]:20f} {o_f[1]:20f} {o_f[2]:20f}")
+    print(f" Difference: {d_f[0]:20f} {d_f[1]:20f} {d_f[2]:20f}")
     print()
-    print(f"    Largest difference is: {max_diff}.")
-    print(f"    Largest difference is at: {ind}.")
-    print(
-        f"    Gromacs: {gmx_forces[ind, 0]:20f} {gmx_forces[ind, 1]:20f} {gmx_forces[ind, 2]:20f}"
-    )
-    print(
-        f"     OpenMM: {omm_forces[ind, 0]:20f} {omm_forces[ind, 1]:20f} {omm_forces[ind, 2]:20f}"
-    )
+    print(f"    Largest violation is at: ({max_ind[0]}, {max_ind[1]}).")
+    print(f"    Gromacs: {max_g[0]:20f} {max_g[1]:20f} {max_g[2]:20f}")
+    print(f"     OpenMM: {max_o[0]:20f} {max_o[1]:20f} {max_o[2]:20f}")
+    print(f" Difference: {max_d[0]:20f} {max_d[1]:20f} {max_d[2]:20f}")
     print()
     failed = True
 

tests/Scripts/gen_openmm.py

@@ -4,23 +4,38 @@
 import martini_openmm as martini
 import numpy as np
 
+
 # do all tests with reference platform to avoid messing
 # with cuda / drivers / etc
+# this also uses full double precision
 platform = mm.Platform.getPlatformByName("Reference")
 
 conf = app.GromacsGroFile("minimized.gro")
 box_vectors = conf.getPeriodicBoxVectors()
-top = martini.GromacsMartiniV2TopFile("system.top", periodicBoxVectors=box_vectors)
+
+# get any defines
+defines = {}
+try:
+    def_file = open("defines.txt")
+    for line in def_file:
+        line = line.strip()
+        defines[line] = True
+except FileNotFoundError:
+    pass
+
+top = martini.GromacsMartiniV2TopFile(
+    "system.top", periodicBoxVectors=box_vectors, defines=defines
+)
 
 system = top.createSystem(nonbondedCutoff=1.1 * u.nanometer)
 integrator = mm.LangevinIntegrator(
     300 * u.kelvin, 1.0 / u.picosecond, 2 * u.femtosecond
 )
+
 simulation = mm.app.Simulation(top.topology, system, integrator, platform)
 
 simulation.context.setPositions(conf.getPositions())
 state = simulation.context.getState(getEnergy=True, getForces=True)
-
 energy = state.getPotentialEnergy().value_in_unit(u.kilojoule_per_mole)
 forces = np.array(state.getForces().value_in_unit(u.kilojoule / u.nanometer / u.mole))
 
@@ -33,6 +48,6 @@
     site = atoms[0]
     forces[site, :] = 0.0
 
-with open('energy.txt', 'w') as efile:
+with open("energy.txt", "w") as efile:
     print(energy, file=efile)
-np.savetxt('forces.txt', forces)
+np.savetxt("forces.txt", forces)

tests/Scripts/run_gmx.sh

@@ -1,9 +1,10 @@
 set -e
 cd gmx
-gmx grompp -f martini_md.mdp -c minimized.gro -p system.top -n system.ndx -o run.tpr -maxwarn 1
-gmx mdrun -deffnm run -rerun minimized.gro -nt 1
-echo pot | gmx energy -f run > energy.txt
-echo 0 | gmx traj -f run -s run -of forces
+# use double precision for everything
+gmx_d grompp -f martini_md.mdp -c minimized.gro -p system.top -n system.ndx -o run.tpr -r minimized.gro -maxwarn 1 -v
+gmx_d mdrun -deffnm run -rerun minimized.gro -nt 1 -v
+echo pot | gmx_d energy -f run > energy.txt
+echo 0 | gmx_d traj -f run -s run -of forces
 rm mdout.mdp
 rm run.edr run.log run.tpr run.trr
 cd ..