Use simpler form of harmonic torsion restraint that respects boundary conditions

openmmapi/include/openmm/CustomTorsionForce.h

@@ -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

@@ -907,9 +907,7 @@ class CharmmPsfFile(object):
         if verbose: print('Adding impropers...')
         # Ick. OpenMM does not have an improper torsion class. Need to
         # construct one from CustomTorsionForce that respects toroidal boundaries
-        energy_function  = 'k*dtheta_torus^2;'
-        energy_function += 'dtheta_torus = dtheta - floor(dtheta/(2*pi)+0.5)*(2*pi);'
-        energy_function += 'dtheta = theta - theta0;'
+        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')