Updated API and documentation to permit triclinic boxes

docs-source/usersguide/theory.rst

@@ -1261,6 +1261,61 @@ Other Features
 ##############
 
 
+Periodic Boundary Conditions
+****************************
+
+Many Force objects support periodic boundary conditions.  They act as if space
+were tiled with infinitely repeating copies of the system, then compute the
+forces acting on a single copy based on the infinite periodic copies.  In most
+(but not all) cases, they apply a cutoff so that each particle only interacts
+with a single copy of each other particle.
+
+OpenMM supports triclinic periodic boxes.  This means the periodicity is defined
+by three vectors, :math:`\mathbf{a}`\ , :math:`\mathbf{b}`\ , and
+:math:`\mathbf{c}`\ .  Given a particle position, the infinite periodic copies
+of that particle are generated by adding vectors of the form
+:math:`i \mathbf{a}+j \mathbf{b}+k \mathbf{c}`\ , where :math:`i`\ ,
+:math:`j`\ , and :math:`k` are arbitrary integers.
+
+The periodic box vectors must be chosen to satisfy certain requirements.
+Roughly speaking, :math:`\mathbf{a}`\ , :math:`\mathbf{b}`\ , and
+:math:`\mathbf{c}` need to "mostly" correspond to the x, y, and z axes.  They
+must have the form
+
+.. math::
+   \mathbf{a} = (a_x, 0, 0)
+
+   \mathbf{b} = (b_x, b_y, 0)
+
+   \mathbf{c} = (c_x, c_y, c_z)
+
+It is always possible to put the box vectors into this form by rotating the
+system until :math:`\mathbf{a}` is parallel to x and :math:`\mathbf{b}` lies in
+the xy plane.
+
+Furthermore, they must obey the following constraints:
+
+.. math::
+   a_x > 0, b_y > 0, c_z > 0
+
+   a_x \ge 2 |b_x|
+
+   a_x \ge 2 |c_x|
+
+   b_y \ge 2 |c_y|
+
+This effectively requires the box vectors to be specified in a particular
+reduced form.  By forming combinations of box vectors (a process known as
+"lattice reduction"), it is always possible to put them in this form without
+changing the periodic system they represent.
+
+These requirements have an important consequence: the periodic unit cell can
+always be treated as an axis-aligned rectangular box of size
+:math:`(a_x, b_y, c_z)`\ .  The remaining non-zero elements of the box vectors
+cause the repeating copies of the system to be staggered relative to each other,
+but they do not affect the shape or size of each copy.  The volume of the unit
+cell is simply given by :math:`a_x b_y c_z`\ .
+
 LocalEnergyMinimizer
 ********************
 

openmmapi/include/openmm/Context.h

@@ -182,8 +182,9 @@ public:
      * Set the vectors defining the axes of the periodic box (measured in nm).  They will affect
      * any Force that uses periodic boundary conditions.
      *
-     * Currently, only rectangular boxes are supported.  This means that a, b, and c must be aligned with the
-     * x, y, and z axes respectively.  Future releases may support arbitrary triclinic boxes.
+     * Triclinic boxes are supported, but the vectors must satisfy certain requirements.  In particular,
+     * a must point in the x direction, b must point "mostly" in the y direction, and c must point "mostly"
+     * in the z direction.  See the documentation for details.
      *
      * @param a      the vector defining the first edge of the periodic box
      * @param b      the vector defining the second edge of the periodic box

openmmapi/include/openmm/System.h

@@ -203,9 +203,6 @@ public:
      * created Context will have its box vectors set to these.  They will affect
      * any Force added to the System that uses periodic boundary conditions.
      *
-     * Currently, only rectangular boxes are supported.  This means that a, b, and c must be aligned with the
-     * x, y, and z axes respectively.  Future releases may support arbitrary triclinic boxes.
-     *
      * @param a      on exit, this contains the vector defining the first edge of the periodic box
      * @param b      on exit, this contains the vector defining the second edge of the periodic box
      * @param c      on exit, this contains the vector defining the third edge of the periodic box
@@ -216,8 +213,9 @@ public:
      * created Context will have its box vectors set to these.  They will affect
      * any Force added to the System that uses periodic boundary conditions.
      *
-     * Currently, only rectangular boxes are supported.  This means that a, b, and c must be aligned with the
-     * x, y, and z axes respectively.  Future releases may support arbitrary triclinic boxes.
+     * Triclinic boxes are supported, but the vectors must satisfy certain requirements.  In particular,
+     * a must point in the x direction, b must point "mostly" in the y direction, and c must point "mostly"
+     * in the z direction.  See the documentation for details.
      *
      * @param a      the vector defining the first edge of the periodic box
      * @param b      the vector defining the second edge of the periodic box

openmmapi/include/openmm/internal/ContextImpl.h

@@ -142,9 +142,6 @@ public:
      * Get the vectors defining the axes of the periodic box (measured in nm).  They will affect
      * any Force that uses periodic boundary conditions.
      *
-     * Currently, only rectangular boxes are supported.  This means that a, b, and c must be aligned with the
-     * x, y, and z axes respectively.  Future releases may support arbitrary triclinic boxes.
-     *
      * @param a      the vector defining the first edge of the periodic box
      * @param b      the vector defining the second edge of the periodic box
      * @param c      the vector defining the third edge of the periodic box
@@ -154,8 +151,9 @@ public:
      * Set the vectors defining the axes of the periodic box (measured in nm).  They will affect
      * any Force that uses periodic boundary conditions.
      *
-     * Currently, only rectangular boxes are supported.  This means that a, b, and c must be aligned with the
-     * x, y, and z axes respectively.  Future releases may support arbitrary triclinic boxes.
+     * Triclinic boxes are supported, but the vectors must satisfy certain requirements.  In particular,
+     * a must point in the x direction, b must point "mostly" in the y direction, and c must point "mostly"
+     * in the z direction.  See the documentation for details.
      *
      * @param a      the vector defining the first edge of the periodic box
      * @param b      the vector defining the second edge of the periodic box

openmmapi/src/ContextImpl.cpp

@@ -40,6 +40,7 @@
 #include "openmm/VirtualSite.h"
 #include "openmm/Context.h"
 #include <algorithm>
+#include <cmath>
 #include <iostream>
 #include <map>
 #include <utility>
@@ -235,10 +236,10 @@ void ContextImpl::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) {
 void ContextImpl::setPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c) {
     if (a[1] != 0.0 || a[2] != 0.0)
         throw OpenMMException("First periodic box vector must be parallel to x.");
-    if (b[0] != 0.0 || b[2] != 0.0)
-        throw OpenMMException("Second periodic box vector must be parallel to y.");
-    if (c[0] != 0.0 || c[1] != 0.0)
-        throw OpenMMException("Third periodic box vector must be parallel to z.");
+    if (b[2] != 0.0)
+        throw OpenMMException("Second periodic box vector must be in the x-y plane.");
+    if (a[0] <= 0.0 || b[1] <= 0.0 || c[2] <= 0.0 || a[0] < 2*fabs(b[0]) || a[0] < 2*fabs(c[0]) || b[1] < 2*fabs(c[1]))
+        throw OpenMMException("Periodic box vectors must be in reduced form.");
     updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, a, b, c);
 }
 

openmmapi/src/System.cpp

@@ -34,6 +34,7 @@
 #include "openmm/System.h"
 #include "openmm/VirtualSite.h"
 #include "openmm/internal/AssertionUtilities.h"
+#include <cmath>
 
 using namespace OpenMM;
 
@@ -111,10 +112,10 @@ void System::getDefaultPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const {
 void System::setDefaultPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c) {
     if (a[1] != 0.0 || a[2] != 0.0)
         throw OpenMMException("First periodic box vector must be parallel to x.");
-    if (b[0] != 0.0 || b[2] != 0.0)
-        throw OpenMMException("Second periodic box vector must be parallel to y.");
-    if (c[0] != 0.0 || c[1] != 0.0)
-        throw OpenMMException("Third periodic box vector must be parallel to z.");
+    if (b[2] != 0.0)
+        throw OpenMMException("Second periodic box vector must be in the x-y plane.");
+    if (a[0] <= 0.0 || b[1] <= 0.0 || c[2] <= 0.0 || a[0] < 2*fabs(b[0]) || a[0] < 2*fabs(c[0]) || b[1] < 2*fabs(c[1]))
+        throw OpenMMException("Periodic box vectors must be in reduced form.");
     periodicBoxVectors[0] = a;
     periodicBoxVectors[1] = b;
     periodicBoxVectors[2] = c;