Periodic box size is specified with box vectors instead of widths

openmmapi/include/StandardMMForceField.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "Force.h"
+#include "Vec3.h"
 #include <map>
 #include <vector>
 #include "internal/windowsExport.h"
@@ -138,23 +139,29 @@ public:
      */
     void setCutoffDistance(double distance);
     /**
-     * Get the dimensions of the periodic box (in nm).  If the NonbondedMethod in use does not use periodic
-     * boundary conditions, these values will have no effect.
+     * Get the vectors which define the axes of the periodic box (measured in nm).  If the NonbondedMethod
+     * in use does not use periodic boundary conditions, these values will have no effect.
      *
-     * @param x      on exit, this contains the width of the periodic box along the x axis
-     * @param y      on exit, this contains the width of the periodic box along the y axis
-     * @param z      on exit, this contains the width of the periodic box along the z axis
+     * 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
      */
-    void getPeriodicBoxSize(double& x, double& y, double& z) const;
+    void getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const;
     /**
-     * Set the dimensions of the periodic box (in nm).  If the NonbondedMethod in use does not use periodic
-     * boundary conditions, these values will have no effect.
+     * Set the vectors which define the axes of the periodic box (measured in nm).  If the NonbondedMethod
+     * in use does not use periodic boundary conditions, these values will have no effect.
+     *
+     * 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 x      the width of the periodic box along the x axis
-     * @param y      the width of the periodic box along the y axis
-     * @param z      the width of the periodic box along the z axis
+     * @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
      */
-    void setPeriodicBoxSize(double x, double y, double z);
+    void setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c);
     /**
      * Get the nonbonded force parameters for an atom.
      * 
@@ -306,7 +313,7 @@ private:
     class NB14Info;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance;
-    double periodicBoxSize[3];
+    Vec3 periodicBoxVectors[3];
 
 // Retarded visual studio compiler complains about being unable to 
 // export private stl class members.

openmmapi/src/StandardMMForceField.cpp

@@ -39,7 +39,9 @@ using namespace OpenMM;
 StandardMMForceField::StandardMMForceField(int numAtoms, int numBonds, int numAngles, int numPeriodicTorsions, int numRBTorsions, int numNonbonded14) :
         atoms(numAtoms), bonds(numBonds), angles(numAngles), periodicTorsions(numPeriodicTorsions), rbTorsions(numRBTorsions), nb14s(numNonbonded14),
         nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
-    periodicBoxSize[0] = periodicBoxSize[1] = periodicBoxSize[2] = 2.0;
+    periodicBoxVectors[0] = Vec3(2, 0, 0);
+    periodicBoxVectors[1] = Vec3(0, 2, 0);
+    periodicBoxVectors[2] = Vec3(0, 0, 2);
 }
 
 StandardMMForceField::NonbondedMethod StandardMMForceField::getNonbondedMethod() const {
@@ -58,16 +60,22 @@ void StandardMMForceField::setCutoffDistance(double distance) {
     cutoffDistance = distance;
 }
 
-void StandardMMForceField::getPeriodicBoxSize(double& x, double& y, double& z) const {
-    x = periodicBoxSize[0];
-    y = periodicBoxSize[1];
-    z = periodicBoxSize[2];
+void StandardMMForceField::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const {
+    a = periodicBoxVectors[0];
+    b = periodicBoxVectors[1];
+    c = periodicBoxVectors[2];
 }
 
-void StandardMMForceField::setPeriodicBoxSize(double x, double y, double z) {
-    periodicBoxSize[0] = x;
-    periodicBoxSize[1] = y;
-    periodicBoxSize[2] = z;
+void StandardMMForceField::setPeriodicBoxVectors(Vec3 a, Vec3 b, 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.");
+    periodicBoxVectors[0] = a;
+    periodicBoxVectors[1] = b;
+    periodicBoxVectors[2] = c;
 }
 
 void StandardMMForceField::getAtomParameters(int index, double& charge, double& radius, double& depth) const {

platforms/cuda/tests/TestCudaStandardMMForceField.cpp

@@ -392,7 +392,7 @@ void testPeriodic() {
     forceField->setNonbondedMethod(StandardMMForceField::CutoffPeriodic);
     const double cutoff = 2.0;
     forceField->setCutoffDistance(cutoff);
-    forceField->setPeriodicBoxSize(4.0, 4.0, 4.0);
+    forceField->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
     system.addForce(forceField);
     OpenMMContext context(system, integrator, platform);
     vector<Vec3> positions(3);

platforms/reference/src/ReferenceKernels.cpp

@@ -217,11 +217,11 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const System& system, c
     }
     nonbondedMethod = CalcStandardMMForceFieldKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
-    double boxSize[3];
-    force.getPeriodicBoxSize(boxSize[0], boxSize[1], boxSize[2]);
-    periodicBoxSize[0] = (RealOpenMM) boxSize[0];
-    periodicBoxSize[1] = (RealOpenMM) boxSize[1];
-    periodicBoxSize[2] = (RealOpenMM) boxSize[2];
+    Vec3 boxVectors[3];
+    force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
+    periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
+    periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
     if (nonbondedMethod == NoCutoff)
         neighborList = NULL;
     else

platforms/reference/tests/TestReferenceStandardMMForceField.cpp

@@ -393,7 +393,7 @@ void testPeriodic() {
     forceField->setNonbondedMethod(StandardMMForceField::CutoffPeriodic);
     const double cutoff = 2.0;
     forceField->setCutoffDistance(cutoff);
-    forceField->setPeriodicBoxSize(4.0, 4.0, 4.0);
+    forceField->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
     system.addForce(forceField);
     OpenMMContext context(system, integrator, platform);
     vector<Vec3> positions(3);