Merge pull request #762 from swails/is_periodic

Implements `usesPeriodicBoundaryConditions` on `Force`s and `System`

openmmapi/include/openmm/NonbondedForce.h

@@ -352,6 +352,17 @@ public:
      * to add new particles or exceptions, only to change the parameters of existing ones.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Returns whether or not this force makes use of periodic boundary
+     * conditions.
+     *
+     * @returns true if force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions() const {
+        return nonbondedMethod == NonbondedForce::CutoffPeriodic ||
+               nonbondedMethod == NonbondedForce::Ewald ||
+               nonbondedMethod == NonbondedForce::PME;
+    }
 protected:
     ForceImpl* createImpl() const;
 private:

openmmapi/include/openmm/PeriodicTorsionForce.h

@@ -108,6 +108,15 @@ public:
      * in a torsion cannot be changed, nor can new torsions be added.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Returns whether or not this force makes use of periodic boundary
+     * conditions.
+     *
+     * @returns true if force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions() const {
+        return false;
+    }
 protected:
     ForceImpl* createImpl() const;
 private:

openmmapi/include/openmm/RBTorsionForce.h

@@ -117,6 +117,15 @@ public:
      * in a torsion cannot be changed, nor can new torsions be added.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Returns whether or not this force makes use of periodic boundary
+     * conditions.
+     *
+     * @returns true if force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions() const {
+        return false;
+    }
 protected:
     ForceImpl* createImpl() const;
 private:

openmmapi/include/openmm/System.h

@@ -224,6 +224,15 @@ public:
      * @param c      the vector defining the third edge of the periodic box
      */
     void setDefaultPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c);
+    /**
+     * Returns whether or not any forces in this System use periodic boundaries.
+     *
+     * If a force in this System does not implement usesPeriodicBoundaryConditions
+     * a OpenMM::OpenMMException is thrown
+     *
+     * @return true if at least one force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions();
 private:
     class ConstraintInfo;
     Vec3 periodicBoxVectors[3];

openmmapi/src/Force.cpp

@@ -49,6 +49,10 @@ void Force::setForceGroup(int group) {
     forceGroup = group;
 }
 
+bool Force::usesPeriodicBoundaryConditions() const {
+    throw OpenMMException("usesPeriodicBoundaryConditions is not implemented");
+}
+
 ForceImpl& Force::getImplInContext(Context& context) {
     const vector<ForceImpl*>& impls = context.getImpl().getForceImpls();
     for (int i = 0; i < (int) impls.size(); i++)

openmmapi/src/System.cpp

@@ -119,3 +119,25 @@ void System::setDefaultPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Ve
     periodicBoxVectors[1] = b;
     periodicBoxVectors[2] = c;
 }
+
+bool System::usesPeriodicBoundaryConditions() {
+
+    bool uses_pbc = false;
+    bool all_forces_implement = true;
+    for (std::vector<Force*>::const_iterator it = forces.begin();
+            it != forces.end(); it++) {
+        try {
+            if ((*it)->usesPeriodicBoundaryConditions())
+                uses_pbc = true;
+        }
+        catch (OpenMMException &e) {
+            all_forces_implement = false;
+        }
+    }
+
+    if (!all_forces_implement && !uses_pbc) {
+        throw OpenMMException("not all forces implement usesPeriodicBoundaryConditions");
+    }
+
+    return uses_pbc;
+}

platforms/reference/tests/TestReferenceAndersenThermostat.cpp

@@ -62,8 +62,9 @@ void testTemperature() {
         forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
     }
     system.addForce(forceField);
-    AndersenThermostat* thermstat = new AndersenThermostat(temp, collisionFreq);
-    system.addForce(thermstat);
+    AndersenThermostat* thermostat = new AndersenThermostat(temp, collisionFreq);
+    system.addForce(thermostat);
+    ASSERT(!thermostat->usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     vector<Vec3> positions(numParticles);
     for (int i = 0; i < numParticles; ++i)
@@ -110,8 +111,8 @@ void testConstraints() {
     system.addConstraint(5, 6, 1);
     system.addConstraint(6, 7, 1);
     system.addConstraint(7, 4, 1);
-    AndersenThermostat* thermstat = new AndersenThermostat(temp, collisionFreq);
-    system.addForce(thermstat);
+    AndersenThermostat* thermostat = new AndersenThermostat(temp, collisionFreq);
+    system.addForce(thermostat);
     Context context(system, integrator, platform);
     vector<Vec3> positions(numParticles);
     positions[0] = Vec3(0, 0, 0);

platforms/reference/tests/TestReferenceCMAPTorsionForce.cpp

@@ -64,6 +64,8 @@ void testCMAPTorsions() {
     periodic->addTorsion(0, 1, 2, 3, 2, M_PI/4, 1.5);
     periodic->addTorsion(1, 2, 3, 4, 3, M_PI/3, 2.0);
     system1.addForce(periodic);
+    ASSERT(!periodic->usesPeriodicBoundaryConditions());
+    ASSERT(!system1.usesPeriodicBoundaryConditions());
     System system2;
     for (int i = 0; i < 5; i++)
         system2.addParticle(1.0);
@@ -81,6 +83,8 @@ void testCMAPTorsions() {
     cmap->addMap(mapSize, mapEnergy);
     cmap->addTorsion(0, 0, 1, 2, 3, 1, 2, 3, 4);
     system2.addForce(cmap);
+    ASSERT(!cmap->usesPeriodicBoundaryConditions());
+    ASSERT(!system2.usesPeriodicBoundaryConditions());
 
     // Set the atoms in various positions, and verify that both systems give equal forces and energy.
 

platforms/reference/tests/TestReferenceCustomAngleForce.cpp

@@ -71,6 +71,8 @@ void testAngles() {
     parameters[1] = 0.5;
     custom->addAngle(1, 2, 3, parameters);
     customSystem.addForce(custom);
+    ASSERT(!custom->usesPeriodicBoundaryConditions());
+    ASSERT(!customSystem.usesPeriodicBoundaryConditions());
 
     // Create an identical system using a HarmonicAngleForce.
 

platforms/reference/tests/TestReferenceCustomBondForce.cpp

@@ -67,6 +67,8 @@ void testBonds() {
     parameters[1] = 0.7;
     forceField->addBond(1, 2, parameters);
     system.addForce(forceField);
+    ASSERT(!forceField->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     vector<Vec3> positions(3);
     positions[0] = Vec3(0, 2, 0);

platforms/reference/tests/TestReferenceCustomCompoundBondForce.cpp

@@ -85,6 +85,8 @@ void testBond() {
     parameters[5] = 1.3;
     custom->addBond(particles, parameters);
     customSystem.addForce(custom);
+    ASSERT(!custom->usesPeriodicBoundaryConditions());
+    ASSERT(!customSystem.usesPeriodicBoundaryConditions());
 
     // Create an identical system using standard forces.
 

platforms/reference/tests/TestReferenceCustomExternalForce.cpp

@@ -67,6 +67,8 @@ void testForce() {
     parameters[1] = 3.0;
     forceField->addParticle(2, parameters);
     system.addForce(forceField);
+    ASSERT(!forceField->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     vector<Vec3> positions(3);
     positions[0] = Vec3(0, 2, 0);

platforms/reference/tests/TestReferenceCustomGBForce.cpp

@@ -135,6 +135,18 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
     custom->setNonbondedMethod(customMethod);
     standardSystem.addForce(obc);
     customSystem.addForce(custom);
+    if (customMethod == CustomGBForce::CutoffPeriodic) {
+        ASSERT(custom->usesPeriodicBoundaryConditions());
+        ASSERT(obc->usesPeriodicBoundaryConditions());
+        ASSERT(standardSystem.usesPeriodicBoundaryConditions());
+        ASSERT(customSystem.usesPeriodicBoundaryConditions());
+    }
+    else {
+        ASSERT(!custom->usesPeriodicBoundaryConditions());
+        ASSERT(!obc->usesPeriodicBoundaryConditions());
+        ASSERT(!standardSystem.usesPeriodicBoundaryConditions());
+        ASSERT(!customSystem.usesPeriodicBoundaryConditions());
+    }
     VerletIntegrator integrator1(0.01);
     VerletIntegrator integrator2(0.01);
     Context context1(standardSystem, integrator1, platform);
@@ -528,7 +540,8 @@ static void buildEthane( GBVIForce* gbviForce, std::vector<Vec3>& positions ) {
        C_gamma  = -0.2863;
        H_radius =  0.125;
        H_gamma  =  0.2437;
-    } else {
+    }
+    else {
        C_radius =  0.215;
        C_gamma  = -1.1087;
        H_radius =  0.150;
@@ -606,7 +619,8 @@ static void buildDimer( GBVIForce* gbviForce, std::vector<Vec3>& positions ) {
        C_gamma  = -0.2863;
        H_radius =  0.125;
        H_gamma  =  0.2437;
-    } else {
+    }
+    else {
        C_radius =  0.215;
        C_gamma  = -1.1087;
        H_radius =  0.150;
@@ -733,7 +747,8 @@ static void findScaledRadii( GBVIForce& gbviForce, std::vector<double> & scaledR
             }
             scaledRadiusJ = radiusJ;
 //             errors++;
-        } else {
+        }
+        else {
 
             double rJ2    = radiusJ*radiusJ;
     
@@ -763,7 +778,8 @@ static void findScaledRadii( GBVIForce& gbviForce, std::vector<double> & scaledR
             scaledRadiusJ  = (radiusJ*radiusJ*radiusJ) - 0.125*scaledRadiusJ; 
             if( scaledRadiusJ > 0.0 ){
                 scaledRadiusJ  = 0.95*pow( scaledRadiusJ, (1.0/3.0) );
-            } else {
+            }
+            else {
                 scaledRadiusJ  = 0.0;
             }
         }
@@ -838,9 +854,11 @@ void testGBVI(GBVIForce::NonbondedMethod gbviMethod, CustomGBForce::NonbondedMet
 
     if( molecule == "Monomer" ){
         buildMonomer( gbvi, positions );
-    } else if( molecule == "Dimer" ){
+    }
+    else if( molecule == "Dimer" ){
         buildDimer( gbvi, positions );
-    } else {
+    }
+    else {
         buildEthane( gbvi, positions );
     }
 

platforms/reference/tests/TestReferenceCustomHbondForce.cpp

@@ -83,6 +83,8 @@ void testHbond() {
     custom->addAcceptor(2, 3, 4, parameters);
     custom->setCutoffDistance(10.0);
     customSystem.addForce(custom);
+    ASSERT(!custom->usesPeriodicBoundaryConditions());
+    ASSERT(!customSystem.usesPeriodicBoundaryConditions());
 
     // Create an identical system using HarmonicBondForce, HarmonicAngleForce, and PeriodicTorsionForce.
 

platforms/reference/tests/TestReferenceCustomManyParticleForce.cpp

@@ -63,6 +63,14 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
         system.addParticle(1.0);
     system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(force);
+    if (force->getNonbondedMethod() == CustomManyParticleForce::CutoffPeriodic) {
+        ASSERT(force->usesPeriodicBoundaryConditions());
+        ASSERT(system.usesPeriodicBoundaryConditions());
+    }
+    else {
+        ASSERT(!force->usesPeriodicBoundaryConditions());
+        ASSERT(!system.usesPeriodicBoundaryConditions());
+    }
     VerletIntegrator integrator(0.001);
     ReferencePlatform platform;
     Context context(system, integrator, platform);

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -184,6 +184,8 @@ void testCutoff() {
     forceField->setNonbondedMethod(CustomNonbondedForce::CutoffNonPeriodic);
     forceField->setCutoffDistance(2.5);
     system.addForce(forceField);
+    ASSERT(!forceField->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     vector<Vec3> positions(3);
     positions[0] = Vec3(0, 0, 0);
@@ -213,6 +215,8 @@ void testPeriodic() {
     forceField->setCutoffDistance(2.0);
     system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
     system.addForce(forceField);
+    ASSERT(forceField->usesPeriodicBoundaryConditions());
+    ASSERT(system.usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     vector<Vec3> positions(3);
     positions[0] = Vec3(0, 0, 0);
@@ -531,6 +535,8 @@ void testCoulombLennardJones() {
     customNonbonded->setNonbondedMethod(CustomNonbondedForce::NoCutoff);
     standardSystem.addForce(standardNonbonded);
     customSystem.addForce(customNonbonded);
+    ASSERT(!customNonbonded->usesPeriodicBoundaryConditions());
+    ASSERT(!customSystem.usesPeriodicBoundaryConditions());
     VerletIntegrator integrator1(0.01);
     VerletIntegrator integrator2(0.01);
     Context context1(standardSystem, integrator1, platform);

platforms/reference/tests/TestReferenceCustomTorsionForce.cpp

@@ -75,6 +75,8 @@ void testTorsions() {
     parameters[1] = 2;
     custom->addTorsion(1, 2, 3, 4, parameters);
     customSystem.addForce(custom);
+    ASSERT(!custom->usesPeriodicBoundaryConditions());
+    ASSERT(!customSystem.usesPeriodicBoundaryConditions());
 
     // Create an identical system using a PeriodicTorsionForce.
 

platforms/reference/tests/TestReferenceGBSAOBCForce.cpp

@@ -58,6 +58,8 @@ void testSingleParticle() {
     GBSAOBCForce* forceField = new GBSAOBCForce();
     forceField->addParticle(0.5, 0.15, 1);
     system.addForce(forceField);
+    ASSERT(!forceField->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     vector<Vec3> positions(1);
     positions[0] = Vec3(0, 0, 0);
@@ -96,6 +98,8 @@ void testGlobalSettings() {
     forceField->setSolventDielectric(solventDielectric);
     forceField->setSurfaceAreaEnergy(surfaceAreaEnergy);
     system.addForce(forceField);
+    ASSERT(!forceField->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     vector<Vec3> positions(1);
     positions[0] = Vec3(0, 0, 0);
@@ -136,22 +140,34 @@ void testCutoffAndPeriodic() {
 
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
     gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
+    ASSERT(!nonbonded->usesPeriodicBoundaryConditions());
+    ASSERT(!gbsa->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     context.setPositions(positions);
     State state1 = context.getState(State::Forces);
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
+    ASSERT(nonbonded->usesPeriodicBoundaryConditions());
+    ASSERT(gbsa->usesPeriodicBoundaryConditions());
+    ASSERT(system.usesPeriodicBoundaryConditions());
     context.reinitialize();
     context.setPositions(positions);
     State state2 = context.getState(State::Forces);
     positions[1][0]+= boxSize;
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
     gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
+    ASSERT(!nonbonded->usesPeriodicBoundaryConditions());
+    ASSERT(!gbsa->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
     context.reinitialize();
     context.setPositions(positions);
     State state3 = context.getState(State::Forces);
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
+    ASSERT(nonbonded->usesPeriodicBoundaryConditions());
+    ASSERT(gbsa->usesPeriodicBoundaryConditions());
+    ASSERT(system.usesPeriodicBoundaryConditions());
     context.reinitialize();
     context.setPositions(positions);
     State state4 = context.getState(State::Forces);

platforms/reference/tests/TestReferenceGBVIForce.cpp

@@ -66,6 +66,8 @@ void testSingleParticle() {
     double gamma          =  1.0;
     forceField->addParticle(charge, radius, gamma);
     system.addForce(forceField);
+    ASSERT(!forceField->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
 
     Context context(system, integrator, platform);
     vector<Vec3> positions(1);
@@ -126,7 +128,8 @@ void testEnergyEthane( int applyBornRadiiScaling ) {
        C_gamma  = -0.2863;
        H_radius =  0.125;
        H_gamma  =  0.2437;
-    } else {
+    }
+    else {
        C_radius =  0.215;
        C_gamma  = -1.1087;
        H_radius =  0.150;
@@ -142,7 +145,8 @@ void testEnergyEthane( int applyBornRadiiScaling ) {
     GBVIForce* forceField = new GBVIForce();
     if( applyBornRadiiScaling ){
         forceField->setBornRadiusScalingMethod( GBVIForce::QuinticSpline );
-    } else {
+    }
+    else {
         forceField->setBornRadiusScalingMethod( GBVIForce::NoScaling );
     }
     for( int i = 0; i < numParticles; i++ ){

platforms/reference/tests/TestReferenceHarmonicAngleForce.cpp

@@ -60,6 +60,8 @@ void testAngles() {
     forceField->addAngle(0, 1, 2, PI_M/3, 1.1);
     forceField->addAngle(1, 2, 3, PI_M/2, 1.2);
     system.addForce(forceField);
+    ASSERT(!forceField->usesPeriodicBoundaryConditions());
+    ASSERT(!system.usesPeriodicBoundaryConditions());
     Context context(system, integrator, platform);
     vector<Vec3> positions(4);
     positions[0] = Vec3(0, 1, 0);