Merge pull request #1761 from jchodera/fix-charmm-impropers

Fix bug where CHARMM impropers did not respect toroidal boundary conditions

devtools/packaging/README.md

 # Packaging OpenMM into ZIP installers
 
+Set your environment variable `TAG` to the git tag for the release:
+```bash
+# OpenMM 7.1.0rc1
+export TAG="9567ddb"
+```
+
 ## Source
 
 Start the docker container:
 ```bash
-docker run -i -t --rm -v `pwd`:/io jchodera/omnia-build-box:cuda80-amd30-clang38 bash
+docker run -i -t --rm -e TAG -v `pwd`:/io jchodera/omnia-build-box:cuda80-amd30-clang38 bash
 ```
 Inside the docker container:
 ```bash
@@ -23,7 +29,7 @@ cp packaging/compressed/* /io
 
 Start the docker container:
 ```bash
-docker run -i -t --rm -v `pwd`:/io jchodera/omnia-build-box:cuda80-amd30-clang38 bash
+docker run -i -t --rm -e TAG -v `pwd`:/io jchodera/omnia-build-box:cuda80-amd30-clang38 bash
 ```
 Inside the docker container:
 ```bash
@@ -50,4 +56,3 @@ source openmm/devtools/packaging/scripts/osx/prepare.sh
 source openmm/devtools/packaging/scripts/osx/build.sh
 source openmm/devtools/packaging/scripts/osx/package.sh
 ```
-

openmmapi/include/openmm/CustomTorsionForce.h

@@ -52,11 +52,11 @@ namespace OpenMM {
  * part of the system definition, while values of global parameters may be modified during a simulation by calling Context::setParameter().
  * Finally, call addTorsion() once for each torsion.  After an torsion has been added, you can modify its parameters by calling setTorsionParameters().
  * This will have no effect on Contexts that already exist unless you call updateParametersInContext().
- * theta is guaranteed to be in the range [-pi,+pi]
+ * Note that theta is guaranteed to be in the range [-pi,+pi], which may cause issues with force discontinuities if the energy function does not respect this domain.
  *
- * As an example, the following code creates a CustomTorsionForce that implements a harmonic potential:
+ * As an example, the following code creates a CustomTorsionForce that implements a periodic potential:
  *
- * <tt>CustomTorsionForce* force = new CustomTorsionForce("0.5*k*(theta-theta0)^2");</tt>
+ * <tt>CustomTorsionForce* force = new CustomTorsionForce("0.5*k*(1-cos(theta-theta0))");</tt>
  *
  * This force depends on two parameters: the spring constant k and equilibrium angle theta0.  The following code defines these parameters:
  *
@@ -65,6 +65,10 @@ namespace OpenMM {
  * force->addPerTorsionParameter("theta0");
  * </pre></tt>
  *
+ * If a harmonic restraint is desired, it is important to be careful of the domain for theta, using an idiom like this:
+ *
+ * <tt>CustomTorsionForce* force = new CustomTorsionForce("0.5*k*min(dtheta, 2*pi-dtheta)^2; dtheta = abs(theta-theta0); pi = 3.1415926535");</tt>
+ *
  * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.
  * Call addEnergyParameterDerivative() to request that the derivative with respect to a particular parameter be
  * computed.  You can then query its value in a Context by calling getState() on it.

wrappers/python/simtk/openmm/app/charmmpsffile.py

@@ -906,8 +906,10 @@ class CharmmPsfFile(object):
 
         if verbose: print('Adding impropers...')
         # Ick. OpenMM does not have an improper torsion class. Need to
-        # construct one from CustomTorsionForce
-        force = mm.CustomTorsionForce('k*(theta-theta0)^2')
+        # construct one from CustomTorsionForce that respects toroidal boundaries
+        energy_function = 'k*min(dtheta, 2*pi-dtheta)^2; dtheta = abs(theta-theta0);'
+        energy_function += 'pi = %f;' % pi
+        force = mm.CustomTorsionForce(energy_function)
         force.addPerTorsionParameter('k')
         force.addPerTorsionParameter('theta0')
         force.setForceGroup(self.IMPROPER_FORCE_GROUP)

wrappers/python/tests/TestCharmmFiles.py

@@ -4,6 +4,7 @@ from simtk.openmm.app import *
 from simtk.openmm import *
 from simtk.unit import *
 import simtk.openmm.app.element as elem
+import math
 import warnings
 
 class TestCharmmFiles(unittest.TestCase):
@@ -206,6 +207,30 @@ class TestCharmmFiles(unittest.TestCase):
         ene_strict     = state_strict.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
         ene_permissive = state_permissive.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
         self.assertAlmostEqual(ene_strict, ene_permissive, delta=0.00001)
+    
+    def test_Impropers(self):
+        """Test CHARMM improper torsions."""
+        psf = CharmmPsfFile('systems/improper.psf')
+        system = psf.createSystem(self.params)
+        force = [f for f in system.getForces() if isinstance(f, CustomTorsionForce)][0]
+        group = force.getForceGroup()
+        integrator = VerletIntegrator(0.001)
+        context = Context(system, integrator, Platform.getPlatformByName("Reference"))
+        angle = 0.1
+        pos1 = [Vec3(0,0,0), Vec3(1,0,0), Vec3(1,1,0), Vec3(0,1,math.tan(angle))] # theta = angle
+        pos2 = [Vec3(0,0,0), Vec3(1,0,0), Vec3(1,1,0), Vec3(2,1,math.tan(angle))] # theta = pi-angle
+        pos3 = [Vec3(0,0,0), Vec3(1,0,0), Vec3(1,1,0), Vec3(2,1,-math.tan(angle))] # theta = -pi+angle
+        for theta0 in (0, math.pi):
+            force.setTorsionParameters(0, 0, 1, 2, 3, [1.0, theta0])
+            force.updateParametersInContext(context)
+            for pos in (pos1, pos2, pos3):
+                context.setPositions(pos)
+                energy = context.getState(getEnergy=True, groups={group}).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
+                if (theta0 == 0 and pos == pos1) or (theta0 == math.pi and pos in (pos2, pos3)):
+                    dtheta = angle
+                else:
+                    dtheta = math.pi-angle
+                self.assertAlmostEqual(energy, dtheta**2, delta=1e-5)
 
 
 if __name__ == '__main__':

wrappers/python/tests/systems/improper.psf

+PSF CMAP CHEQ
+
+       2 !NTITLE
+**                                                                              
+*  Minimal file defining a single improper torsion                     
+
+       4 !NATOM
+       1 AAL  3    ALA  C      32   0.340000       12.0110           0   0.00000     -0.301140E-02
+       2 AAL  3    ALA  CA     22   0.700000E-01   12.0110           0   0.00000     -0.301140E-02
+       3 AAL  3    ALA  OT2    72  -0.670000       15.9990           0   0.00000     -0.301140E-02
+       4 AAL  3    ALA  OT1    72  -0.670000       15.9990           0   0.00000     -0.301140E-02
+
+       0 !NBOND: bonds
+
+       0 !NTHETA: angles
+
+       0 !NPHI: dihedrals
+
+       1 !NIMPHI: impropers
+       1       2       3       4
+
+       0 !NDON: donors
+
+       0 !NACC: acceptors
+
+       0 !NNB
+
+       0       0       0       0
+
+       1       0 !NGRP NST2
+       0       0       0
+
+       1 !MOLNT
+       1       1       1       1
+
+       0       0 !NUMLP NUMLPH
+
+       0 !NCRTERM: cross-terms