Periodic box size is specified by System rather than NonbondedForce

examples/HelloWaterBox.cpp

@@ -248,7 +248,7 @@ myInitializeOpenMM( int                 numWatersAlongEdge,
     // Create periodic box
     nonbond.setNonbondedMethod(OpenMM::NonbondedForce::CutoffPeriodic);
     nonbond.setCutoffDistance(CutoffDistanceInAng * OpenMM::NmPerAngstrom);
-    nonbond.setPeriodicBoxVectors(Vec3(boxEdgeLengthInNm,0,0), 
+    system.setPeriodicBoxVectors(Vec3(boxEdgeLengthInNm,0,0),
                                   Vec3(0,boxEdgeLengthInNm,0), 
                                   Vec3(0,0,boxEdgeLengthInNm));
 

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -170,30 +170,6 @@ public:
      * is NoCutoff, this value will have no effect.
      */
     void setCutoffDistance(double distance);
-    /**
-     * 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.
-     *
-     * 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 getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const;
-    /**
-     * 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 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 setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c);
     /**
      * Add a new per-particle parmeter that the interaction may depend on.
      *
@@ -369,7 +345,6 @@ private:
     class FunctionInfo;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance;
-    Vec3 periodicBoxVectors[3];
     std::string energyExpression;
     std::vector<ParameterInfo> parameters;
     std::vector<GlobalParameterInfo> globalParameters;

openmmapi/include/openmm/GBSAOBCForce.h

@@ -46,13 +46,29 @@ namespace OpenMM {
  * be exactly equal to the number of particles in the System, or else an exception will be thrown when you
  * try to create a Context.  After a particle has been added, you can modify its force field parameters
  * by calling setParticleParameters().
- * <p>
- * If the System also contains a NonbondedForce, this force will use the cutoffs
- * and periodic boundary conditions specified in it.
  */
 
 class OPENMM_EXPORT GBSAOBCForce : public Force {
 public:
+    /**
+     * This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
+     */
+    enum NonbondedMethod {
+        /**
+         * No cutoff is applied to nonbonded interactions.  The full set of N^2 interactions is computed exactly.
+         * This necessarily means that periodic boundary conditions cannot be used.  This is the default.
+         */
+        NoCutoff = 0,
+        /**
+         * Interactions beyond the cutoff distance are ignored.
+         */
+        CutoffNonPeriodic = 1,
+        /**
+         * Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
+         * each other particle.  Interactions beyond the cutoff distance are ignored.
+         */
+        CutoffPeriodic = 2,
+    };
     /*
      * Create a GBSAOBCForce.
      */
@@ -115,26 +131,31 @@ public:
     void setSoluteDielectric(double dielectric) {
         soluteDielectric = dielectric;
     }
+    /**
+     * Get the method used for handling long range nonbonded interactions.
+     */
+    NonbondedMethod getNonbondedMethod() const;
+    /**
+     * Set the method used for handling long range nonbonded interactions.
+     */
+    void setNonbondedMethod(NonbondedMethod method);
+    /**
+     * Get the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     */
+    double getCutoffDistance() const;
+    /**
+     * Set the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     */
+    void setCutoffDistance(double distance);
 protected:
     ForceImpl* createImpl();
 private:
     class ParticleInfo;
-    double solventDielectric, soluteDielectric;
-
-// Retarded visual studio compiler complains about being unable to 
-// export private stl class members.
-// This stanza explains that it should temporarily shut up.
-#if defined(_MSC_VER)
-#pragma warning(push)
-#pragma warning(disable:4251)
-#endif
-
+    NonbondedMethod nonbondedMethod;
+    double cutoffDistance, solventDielectric, soluteDielectric;
     std::vector<ParticleInfo> particles;
-
-#if defined(_MSC_VER)
-#pragma warning(pop)
-#endif
-
 };
 
 class GBSAOBCForce::ParticleInfo {

openmmapi/include/openmm/GBVIForce.h

@@ -46,13 +46,29 @@ namespace OpenMM {
  * be exactly equal to the number of particles in the System, or else an exception will be thrown when you
  * try to create a Context.  After a particle has been added, you can modify its force field parameters
  * by calling setParticleParameters().
- * <p>
- * If the System also contains a NonbondedForce, this force will use the cutoffs
- * and periodic boundary conditions specified in it.
  */
 
 class OPENMM_EXPORT GBVIForce : public Force {
 public:
+    /**
+     * This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
+     */
+    enum NonbondedMethod {
+        /**
+         * No cutoff is applied to nonbonded interactions.  The full set of N^2 interactions is computed exactly.
+         * This necessarily means that periodic boundary conditions cannot be used.  This is the default.
+         */
+        NoCutoff = 0,
+        /**
+         * Interactions beyond the cutoff distance are ignored.
+         */
+        CutoffNonPeriodic = 1,
+        /**
+         * Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
+         * each other particle.  Interactions beyond the cutoff distance are ignored.
+         */
+        CutoffPeriodic = 2,
+    };
     /*
      * Create a GBVIForce.
      */
@@ -115,26 +131,31 @@ public:
     void setSoluteDielectric(double dielectric) {
         soluteDielectric = dielectric;
     }
+    /**
+     * Get the method used for handling long range nonbonded interactions.
+     */
+    NonbondedMethod getNonbondedMethod() const;
+    /**
+     * Set the method used for handling long range nonbonded interactions.
+     */
+    void setNonbondedMethod(NonbondedMethod method);
+    /**
+     * Get the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     */
+    double getCutoffDistance() const;
+    /**
+     * Set the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     */
+    void setCutoffDistance(double distance);
 protected:
     ForceImpl* createImpl();
 private:
     class ParticleInfo;
-    double solventDielectric, soluteDielectric;
-
-// Retarded visual studio compiler complains about being unable to 
-// export private stl class members.
-// This stanza explains that it should temporarily shut up.
-#if defined(_MSC_VER)
-#pragma warning(push)
-#pragma warning(disable:4251)
-#endif
-
+    NonbondedMethod nonbondedMethod;
+    double cutoffDistance, solventDielectric, soluteDielectric;
     std::vector<ParticleInfo> particles;
-
-#if defined(_MSC_VER)
-#pragma warning(pop)
-#endif
-
 };
 
 class GBVIForce::ParticleInfo {

openmmapi/include/openmm/NonbondedForce.h

@@ -33,7 +33,6 @@
  * -------------------------------------------------------------------------- */
 
 #include "Force.h"
-#include "Vec3.h"
 #include <map>
 #include <set>
 #include <utility>
@@ -155,30 +154,6 @@ public:
      * rigorous guarantee that all forces on all atoms will be less than the tolerance, however.
      */
     void setEwaldErrorTolerance(double tol);
-    /**
-     * 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.
-     *
-     * 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 getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const;
-    /**
-     * 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 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 setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c);
     /**
      * Add the nonbonded force parameters for a particle.  This should be called once for each particle
      * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.
@@ -273,24 +248,10 @@ private:
     class ExceptionInfo;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance, rfDielectric, ewaldErrorTol;
-    Vec3 periodicBoxVectors[3];
     void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
-
-
-// Retarded visual studio compiler complains about being unable to 
-// export private stl class members.
-// This stanza explains that it should temporarily shut up.
-#if defined(_MSC_VER)
-#pragma warning(push)
-#pragma warning(disable:4251)
-#endif
     std::vector<ParticleInfo> particles;
     std::vector<ExceptionInfo> exceptions;
     std::map<std::pair<int, int>, int> exceptionMap;
-#if defined(_MSC_VER)
-#pragma warning(pop)
-#endif
-
 };
 
 class NonbondedForce::ParticleInfo {

openmmapi/include/openmm/System.h

@@ -32,6 +32,7 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include "Vec3.h"
 #include <vector>
 #include "internal/windowsExport.h"
 
@@ -41,12 +42,13 @@ class OPENMM_EXPORT Force;
 
 /**
  * This class represents a molecular system.  The definition of a System involves
- * three elements:
+ * four elements:
  * 
  * <ol>
  * <li>The set of particles in the system</li>
  * <li>The forces acting on them</li>
  * <li>Pairs of particles whose separation should be constrained to a fixed value</li>
+ * <li>For periodic systems, the dimensions of the periodic box</li>
  * </ol>
  * 
  * The particles and constraints are defined directly by the System object, while
@@ -162,25 +164,36 @@ public:
     Force& getForce(int index) {
         return *forces[index];
     }
+    /**
+     * Get the vectors which define the axes of the periodic box (measured in nm).  These 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
+     */
+    void getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const;
+    /**
+     * Set the vectors which define the axes of the periodic box (measured in nm).  These 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      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 setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c);
 private:
     class ConstraintInfo;
-
-// Retarded visual studio compiler complains about being unable to 
-// export private stl class members.
-// This stanza explains that it should temporarily shut up.
-#if defined(_MSC_VER)
-#pragma warning(push)
-#pragma warning(disable:4251)
-#endif
-
+    Vec3 periodicBoxVectors[3];
     std::vector<double> masses;
     std::vector<ConstraintInfo> constraints;
     std::vector<Force*> forces;
-
-#if defined(_MSC_VER)
-#pragma warning(pop)
-#endif
-
 };
 
 class System::ConstraintInfo {

openmmapi/src/CustomNonbondedForce.cpp

@@ -47,9 +47,6 @@ using std::stringstream;
 using std::vector;
 
 CustomNonbondedForce::CustomNonbondedForce(const string& energy) : energyExpression(energy), nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
-    periodicBoxVectors[0] = Vec3(2, 0, 0);
-    periodicBoxVectors[1] = Vec3(0, 2, 0);
-    periodicBoxVectors[2] = Vec3(0, 0, 2);
 }
 
 const string& CustomNonbondedForce::getEnergyFunction() const {
@@ -76,24 +73,6 @@ void CustomNonbondedForce::setCutoffDistance(double distance) {
     cutoffDistance = distance;
 }
 
-void CustomNonbondedForce::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const {
-    a = periodicBoxVectors[0];
-    b = periodicBoxVectors[1];
-    c = periodicBoxVectors[2];
-}
-
-void CustomNonbondedForce::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;
-}
-
 int CustomNonbondedForce::addParameter(const string& name, const string& combiningRule) {
     parameters.push_back(ParameterInfo(name, combiningRule));
     return parameters.size()-1;

openmmapi/src/GBSAOBCForce.cpp

@@ -36,7 +36,7 @@
 
 using namespace OpenMM;
 
-GBSAOBCForce::GBSAOBCForce() : solventDielectric(78.3), soluteDielectric(1.0) {
+GBSAOBCForce::GBSAOBCForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), solventDielectric(78.3), soluteDielectric(1.0) {
 }
 
 int GBSAOBCForce::addParticle(double charge, double radius, double scalingFactor) {
@@ -56,6 +56,22 @@ void GBSAOBCForce::setParticleParameters(int index, double charge, double radius
     particles[index].scalingFactor = scalingFactor;
 }
 
+GBSAOBCForce::NonbondedMethod GBSAOBCForce::getNonbondedMethod() const {
+    return nonbondedMethod;
+}
+
+void GBSAOBCForce::setNonbondedMethod(NonbondedMethod method) {
+    nonbondedMethod = method;
+}
+
+double GBSAOBCForce::getCutoffDistance() const {
+    return cutoffDistance;
+}
+
+void GBSAOBCForce::setCutoffDistance(double distance) {
+    cutoffDistance = distance;
+}
+
 ForceImpl* GBSAOBCForce::createImpl() {
     return new GBSAOBCForceImpl(*this);
 }

openmmapi/src/GBVIForce.cpp

@@ -36,7 +36,7 @@
 
 using namespace OpenMM;
 
-GBVIForce::GBVIForce() : solventDielectric(78.3), soluteDielectric(1.0) {
+GBVIForce::GBVIForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), solventDielectric(78.3), soluteDielectric(1.0) {
 }
 
 int GBVIForce::addParticle(double charge, double radius, double gamma) {
@@ -56,6 +56,22 @@ void GBVIForce::setParticleParameters(int index, double charge, double radius, d
     particles[index].gamma  = gamma;
 }
 
+GBVIForce::NonbondedMethod GBVIForce::getNonbondedMethod() const {
+    return nonbondedMethod;
+}
+
+void GBVIForce::setNonbondedMethod(NonbondedMethod method) {
+    nonbondedMethod = method;
+}
+
+double GBVIForce::getCutoffDistance() const {
+    return cutoffDistance;
+}
+
+void GBVIForce::setCutoffDistance(double distance) {
+    cutoffDistance = distance;
+}
+
 ForceImpl* GBVIForce::createImpl() {
     return new GBVIForceImpl(*this);
 }

openmmapi/src/NonbondedForce.cpp

@@ -47,9 +47,6 @@ using std::stringstream;
 using std::vector;
 
 NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), rfDielectric(78.3), ewaldErrorTol(1e-4) {
-    periodicBoxVectors[0] = Vec3(2, 0, 0);
-    periodicBoxVectors[1] = Vec3(0, 2, 0);
-    periodicBoxVectors[2] = Vec3(0, 0, 2);
 }
 
 NonbondedForce::NonbondedMethod NonbondedForce::getNonbondedMethod() const {
@@ -85,24 +82,6 @@ void NonbondedForce::setEwaldErrorTolerance(double tol)
     ewaldErrorTol = tol;
 }
 
-void NonbondedForce::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const {
-    a = periodicBoxVectors[0];
-    b = periodicBoxVectors[1];
-    c = periodicBoxVectors[2];
-}
-
-void NonbondedForce::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;
-}
-
 int NonbondedForce::addParticle(double charge, double sigma, double epsilon) {
     particles.push_back(ParticleInfo(charge, sigma, epsilon));
     return particles.size()-1;

openmmapi/src/System.cpp

@@ -30,11 +30,15 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/Force.h"
+#include "openmm/OpenMMException.h"
 #include "openmm/System.h"
 
 using namespace OpenMM;
 
 System::System() {
+    periodicBoxVectors[0] = Vec3(2, 0, 0);
+    periodicBoxVectors[1] = Vec3(0, 2, 0);
+    periodicBoxVectors[2] = Vec3(0, 0, 2);
 }
 
 System::~System() {
@@ -58,3 +62,21 @@ void System::setConstraintParameters(int index, int particle1, int particle2, do
     constraints[index].particle2 = particle2;
     constraints[index].distance = distance;
 }
+
+void System::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const {
+    a = periodicBoxVectors[0];
+    b = periodicBoxVectors[1];
+    c = periodicBoxVectors[2];
+}
+
+void System::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;
+}

platforms/cuda/src/CudaKernels.cpp

@@ -288,7 +288,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             exclusionList[exclusions[i].second].push_back(exclusions[i].first);
         }
         Vec3 boxVectors[3];
-        force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
         gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
         CudaNonbondedMethod method = NO_CUTOFF;
         if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
@@ -410,13 +410,13 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
         exclusionList[exclusions[i].first].push_back(exclusions[i].second);
         exclusionList[exclusions[i].second].push_back(exclusions[i].first);
     }
+    Vec3 boxVectors[3];
+    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
     CudaNonbondedMethod method = NO_CUTOFF;
     if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
         method = CUTOFF;
     if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
-        Vec3 boxVectors[3];
-        force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-        gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
         method = PERIODIC;
     }
     data.customNonbondedMethod = method;

platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp

@@ -186,7 +186,7 @@ void testPeriodic() {
     forceField->addParticle(vector<double>());
     forceField->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
     forceField->setCutoffDistance(2.0);
-    forceField->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    system.setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
     system.addForce(forceField);
     Context context(system, integrator, platform);
     vector<Vec3> positions(3);
@@ -251,4 +251,3 @@ int main() {
     cout << "Done" << endl;
     return 0;
 }
-

platforms/cuda/tests/TestCudaGBSAOBCForce.cpp

@@ -90,7 +90,8 @@ void testCutoffAndPeriodic() {
     const double cutoffDistance = 3.0;
     const double boxSize = 10.0;
     nonbonded->setCutoffDistance(cutoffDistance);
-    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    gbsa->setCutoffDistance(cutoffDistance);
+    system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(gbsa);
     system.addForce(nonbonded);
     vector<Vec3> positions(2);
@@ -100,19 +101,23 @@ void testCutoffAndPeriodic() {
     // Calculate the forces for both cutoff and periodic with two different atom positions.
 
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
     Context context(system, integrator, cuda);
     context.setPositions(positions);
     State state1 = context.getState(State::Forces);
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
     context.reinitialize();
     context.setPositions(positions);
     State state2 = context.getState(State::Forces);
     positions[1][0]+= boxSize;
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
     context.reinitialize();
     context.setPositions(positions);
     State state3 = context.getState(State::Forces);
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
     context.reinitialize();
     context.setPositions(positions);
     State state4 = context.getState(State::Forces);
@@ -127,7 +132,7 @@ void testCutoffAndPeriodic() {
     ASSERT_EQUAL_VEC(state3.getForces()[1], Vec3(0, 0, 0), 0.01);
 }
 
-void testForce(int numParticles, NonbondedForce::NonbondedMethod method) {
+void testForce(int numParticles, NonbondedForce::NonbondedMethod method, GBSAOBCForce::NonbondedMethod method2) {
     CudaPlatform cuda;
     ReferencePlatform reference;
     System system;
@@ -141,11 +146,13 @@ void testForce(int numParticles, NonbondedForce::NonbondedMethod method) {
         nonbonded->addParticle(charge, 1, 0);
     }
     nonbonded->setNonbondedMethod(method);
+    gbsa->setNonbondedMethod(method2);
     nonbonded->setCutoffDistance(3.0);
+    gbsa->setCutoffDistance(3.0);
     int grid = (int) floor(0.5+pow(numParticles, 1.0/3.0));
     if (method == NonbondedForce::CutoffPeriodic) {
         double boxSize = (grid+1)*1.1;
-        nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+        system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     }
     system.addForce(gbsa);
     system.addForce(nonbonded);
@@ -208,9 +215,9 @@ int main() {
         testSingleParticle();
         testCutoffAndPeriodic();
         for (int i = 5; i < 11; i++) {
-            testForce(i*i*i, NonbondedForce::NoCutoff);
-            testForce(i*i*i, NonbondedForce::CutoffNonPeriodic);
-            testForce(i*i*i, NonbondedForce::CutoffPeriodic);
+            testForce(i*i*i, NonbondedForce::NoCutoff, GBSAOBCForce::NoCutoff);
+            testForce(i*i*i, NonbondedForce::CutoffNonPeriodic, GBSAOBCForce::CutoffNonPeriodic);
+            testForce(i*i*i, NonbondedForce::CutoffPeriodic, GBSAOBCForce::CutoffPeriodic);
         }
     }
     catch(const exception& e) {

platforms/cuda/tests/TestCudaNonbondedForce.cpp

@@ -339,7 +339,7 @@ void testPeriodic() {
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     const double cutoff = 2.0;
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    system.setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
     system.addForce(nonbonded);
     Context context(system, integrator, platform);
     vector<Vec3> positions(3);
@@ -372,7 +372,7 @@ void testEwald() {
     const double cutoff = 2.0;
     nonbonded->setCutoffDistance(cutoff);
     nonbonded->setEwaldErrorTolerance(TOL);
-    nonbonded->setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
     system.addForce(nonbonded);
     Context context(system, integrator, platform);
     vector<Vec3> positions(2);
@@ -447,7 +447,7 @@ void testLargeSystem() {
     // Now do the same thing with periodic boundary conditions.
 
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
-    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     cudaContext.reinitialize();
     referenceContext.reinitialize();
     cudaContext.setPositions(positions);
@@ -485,7 +485,7 @@ void testBlockInteractions(bool periodic) {
     }
     nonbonded->setNonbondedMethod(periodic ? NonbondedForce::CutoffPeriodic : NonbondedForce::CutoffNonPeriodic);
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(nonbonded);
     Context context(system, integrator, cuda);
     context.setPositions(positions);

platforms/cuda/tests/TstCudaUsingParameterFile.cpp

@@ -1529,7 +1529,7 @@ if( log ){
                   boxIndex++;
                }
             }
-            nonbondedForce->setPeriodicBoxVectors( box[0], box[1], box[2] );
+            system.setPeriodicBoxVectors( box[0], box[1], box[2] );
 
          } else if( field.compare( "NonbondedForceMethod" ) == 0 ){
             std::string nonbondedForceMethod = tokens[1];
@@ -1569,7 +1569,7 @@ if( log ){
 
       (void) fprintf( log, "CutoffDistance %14.7e\n", nonbondedForce->getCutoffDistance() );
       Vec3 a, b, c;
-      nonbondedForce->getPeriodicBoxVectors( a, b, c);
+      system.getPeriodicBoxVectors( a, b, c);
       (void) fprintf( log, "Box [%14.7f %14.7f %14.7f]\n    [%14.7f %14.7f %14.7f]\n    [%14.7f %14.7f %14.7f]\n",
                       a[0], a[1], a[2], b[0], b[1], b[2], c[0], c[1], c[2] );
   

platforms/reference/src/ReferenceKernels.cpp

@@ -389,7 +389,7 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
     nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
     Vec3 boxVectors[3];
-    force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    system.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];
@@ -598,7 +598,7 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
     nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
     Vec3 boxVectors[3];
-    force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    system.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];
@@ -700,25 +700,16 @@ void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBS
     obcParameters->setScaledRadiusFactors(scaleFactors);
     obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
     obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
-
-    // If there is a NonbondedForce in this system, use it to initialize cutoffs and periodic boundary conditions.
-
-    for (int i = 0; i < system.getNumForces(); i++) {
-        const NonbondedForce* nonbonded = dynamic_cast<const NonbondedForce*>(&system.getForce(i));
-        if (nonbonded != NULL) {
-            if (nonbonded->getNonbondedMethod() != NonbondedForce::NoCutoff)
-                obcParameters->setUseCutoff(static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
-            if (nonbonded->getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
-                Vec3 boxVectors[3];
-                nonbonded->getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-                RealOpenMM periodicBoxSize[3];
-                periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
-                periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
-                periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
-                obcParameters->setPeriodic(periodicBoxSize);
-            }
-            break;
-        }
+    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
+        obcParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
+    if (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic) {
+        Vec3 boxVectors[3];
+        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+        RealOpenMM periodicBoxSize[3];
+        periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
+        periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
+        periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
+        obcParameters->setPeriodic(periodicBoxSize);
     }
     obc = new CpuObc(obcParameters);
     obc->setIncludeAceApproximation(true);
@@ -745,14 +736,11 @@ ReferenceCalcGBVIForceKernel::~ReferenceCalcGBVIForceKernel() {
 }
 
 void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIForce& force, const std::vector<double> & inputScaledRadii ) {
-
     int numParticles = system.getNumParticles();
-
     charges.resize(numParticles);
     vector<RealOpenMM> atomicRadii(numParticles);
     vector<RealOpenMM> scaledRadii(numParticles);
     vector<RealOpenMM> gammas(numParticles);
-
     for (int i = 0; i < numParticles; ++i) {
         double charge, radius, gamma;
         force.getParticleParameters(i, charge, radius, gamma);
@@ -761,35 +749,24 @@ void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIFo
         gammas[i]        = static_cast<RealOpenMM>(gamma);
         scaledRadii[i]   = static_cast<RealOpenMM>(inputScaledRadii[i]);
     }
-
     GBVIParameters * gBVIParameters = new GBVIParameters(numParticles);
     gBVIParameters->setAtomicRadii(atomicRadii);
     gBVIParameters->setGammaParameters(gammas);
     gBVIParameters->setScaledRadii(scaledRadii);
-    gBVIParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
-    gBVIParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
-
-    // If there is a NonbondedForce in this system, use it to initialize cutoffs and periodic boundary conditions.
-
-    for (int i = 0; i < system.getNumForces(); i++) {
-        const NonbondedForce* nonbonded = dynamic_cast<const NonbondedForce*>(&system.getForce(i));
-        if (nonbonded != NULL) {
-            if (nonbonded->getNonbondedMethod() != NonbondedForce::NoCutoff)
-                gBVIParameters->setUseCutoff( static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
-            if (nonbonded->getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
-                Vec3 boxVectors[3];
-                nonbonded->getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-                RealOpenMM periodicBoxSize[3];
-                periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
-                periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
-                periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
-                gBVIParameters->setPeriodic(periodicBoxSize);
-            }
-            break;
-        }
+    gBVIParameters->setSolventDielectric(static_cast<RealOpenMM>(force.getSolventDielectric()));
+    gBVIParameters->setSoluteDielectric(static_cast<RealOpenMM>(force.getSoluteDielectric()));
+    if (force.getNonbondedMethod() != GBVIForce::NoCutoff)
+        gBVIParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
+    if (force.getNonbondedMethod() == GBVIForce::CutoffPeriodic) {
+        Vec3 boxVectors[3];
+        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+        RealOpenMM periodicBoxSize[3];
+        periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
+        periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
+        periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
+        gBVIParameters->setPeriodic(periodicBoxSize);
     }
     gbvi = new CpuGBVI(gBVIParameters);
-
 }
 
 void ReferenceCalcGBVIForceKernel::executeForces(ContextImpl& context) {

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -186,7 +186,7 @@ void testPeriodic() {
     forceField->addParticle(vector<double>());
     forceField->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
     forceField->setCutoffDistance(2.0);
-    forceField->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    system.setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
     system.addForce(forceField);
     Context context(system, integrator, platform);
     vector<Vec3> positions(3);

platforms/reference/tests/TestReferenceEwald.cpp

@@ -67,7 +67,7 @@ void testLargeSystem() {
     nonbonded->setNonbondedMethod(NonbondedForce::PME);
     const double cutoff = 1.0;
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setPeriodicBoxVectors(Vec3(2.82, 0, 0), Vec3(0, 2.82, 0), Vec3(0, 0, 2.82));
+    system.setPeriodicBoxVectors(Vec3(2.82, 0, 0), Vec3(0, 2.82, 0), Vec3(0, 0, 2.82));
     nonbonded->setEwaldErrorTolerance(TOL);
     system.addForce(nonbonded);
     Context context(system, integrator, platform);
@@ -99,7 +99,7 @@ void testEwald() {
     nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
     const double cutoff = 2.0;
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
     nonbonded->setEwaldErrorTolerance(TOL);
     system.addForce(nonbonded);
     Context context(system, integrator, platform);
@@ -134,7 +134,7 @@ void testPME() {
     nonbonded->setNonbondedMethod(NonbondedForce::PME);
     const double cutoff = 0.8;
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setPeriodicBoxVectors(Vec3(1.86206, 0, 0), Vec3(0, 1.86206, 0), Vec3(0, 0, 1.86206));
+    system.setPeriodicBoxVectors(Vec3(1.86206, 0, 0), Vec3(0, 1.86206, 0), Vec3(0, 0, 1.86206));
     nonbonded->setEwaldErrorTolerance(TOL);
     system.addForce(nonbonded);
     Context context(system, integrator, platform);

platforms/reference/tests/TestReferenceGBSAOBCForce.cpp

@@ -86,7 +86,8 @@ void testCutoffAndPeriodic() {
     const double cutoffDistance = 3.0;
     const double boxSize = 10.0;
     nonbonded->setCutoffDistance(cutoffDistance);
-    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    gbsa->setCutoffDistance(cutoffDistance);
+    system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(gbsa);
     system.addForce(nonbonded);
     vector<Vec3> positions(2);
@@ -96,19 +97,23 @@ void testCutoffAndPeriodic() {
     // Calculate the forces for both cutoff and periodic with two different atom positions.
 
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
     Context context(system, integrator, platform);
     context.setPositions(positions);
     State state1 = context.getState(State::Forces);
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
     context.reinitialize();
     context.setPositions(positions);
     State state2 = context.getState(State::Forces);
     positions[1][0]+= boxSize;
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
     context.reinitialize();
     context.setPositions(positions);
     State state3 = context.getState(State::Forces);
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
     context.reinitialize();
     context.setPositions(positions);
     State state4 = context.getState(State::Forces);