Added lots of checks for invalid parameter values (#3286)

* Added lots of checks for invalid parameter values

* Fixed test failures

* More checks for incorrect parameters

* Fixed test failures

openmmapi/src/HarmonicBondForceImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -29,14 +29,14 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/HarmonicBondForceImpl.h"
 #include "openmm/kernels.h"
+#include <sstream>
 
 using namespace OpenMM;
-using std::pair;
-using std::vector;
-using std::set;
+using namespace std;
 
 HarmonicBondForceImpl::HarmonicBondForceImpl(const HarmonicBondForce& owner) : owner(owner) {
 }
@@ -45,6 +45,24 @@ HarmonicBondForceImpl::~HarmonicBondForceImpl() {
 }
 
 void HarmonicBondForceImpl::initialize(ContextImpl& context) {
+    const System& system = context.getSystem();
+    for (int i = 0; i < owner.getNumBonds(); i++) {
+        int particle[2];
+        double length, k;
+        owner.getBondParameters(i, particle[0], particle[1], length, k);
+        for (int j = 0; j < 2; j++) {
+            if (particle[j] < 0 || particle[j] >= system.getNumParticles()) {
+                stringstream msg;
+                msg << "HarmonicBondForce: Illegal particle index for a bond: ";
+                msg << particle[j];
+                throw OpenMMException(msg.str());
+            }
+        }
+        if (length < 0)
+            throw OpenMMException("HarmonicBondForce: bond length cannot be negative");
+        if (k < 0)
+            throw OpenMMException("HarmonicBondForce: force constant cannot be negative");
+    }
     kernel = context.getPlatform().createKernel(CalcHarmonicBondForceKernel::Name(), context);
     kernel.getAs<CalcHarmonicBondForceKernel>().initialize(context.getSystem(), owner);
 }

openmmapi/src/Integrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2020 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -32,6 +32,7 @@
 #include "sfmt/SFMT.h"
 #include "SimTKOpenMMRealType.h"
 #include "openmm/Integrator.h"
+#include "openmm/OpenMMException.h"
 #include "openmm/System.h"
 #include "openmm/internal/ContextImpl.h"
 
@@ -58,6 +59,8 @@ double Integrator::getStepSize() const {
 }
 
 void Integrator::setStepSize(double size) {
+    if (size < 0)
+        throw OpenMMException("Step size cannot be negative");
     stepSize = size;
 }
 
@@ -66,6 +69,8 @@ double Integrator::getConstraintTolerance() const {
 }
 
 void Integrator::setConstraintTolerance(double tol) {
+    if (tol <= 0)
+        throw OpenMMException("Constraint tolerance must be positive");
     constraintTol = tol;
 }
 

openmmapi/src/LangevinIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -57,6 +57,18 @@ void LangevinIntegrator::initialize(ContextImpl& contextRef) {
     kernel.getAs<IntegrateLangevinStepKernel>().initialize(contextRef.getSystem(), *this);
 }
 
+void LangevinIntegrator::setTemperature(double temp) {
+    if (temp < 0)
+        throw OpenMMException("Temperature cannot be negative");
+    temperature = temp;
+}
+
+void LangevinIntegrator::setFriction(double coeff) {
+    if (coeff < 0)
+        throw OpenMMException("Friction cannot be negative");
+    friction = coeff;
+}
+
 void LangevinIntegrator::cleanup() {
     kernel = Kernel();
 }

openmmapi/src/LangevinMiddleIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -57,6 +57,18 @@ void LangevinMiddleIntegrator::initialize(ContextImpl& contextRef) {
     kernel.getAs<IntegrateLangevinMiddleStepKernel>().initialize(contextRef.getSystem(), *this);
 }
 
+void LangevinMiddleIntegrator::setTemperature(double temp) {
+    if (temp < 0)
+        throw OpenMMException("Temperature cannot be negative");
+    temperature = temp;
+}
+
+void LangevinMiddleIntegrator::setFriction(double coeff) {
+    if (coeff < 0)
+        throw OpenMMException("Friction cannot be negative");
+    friction = coeff;
+}
+
 void LangevinMiddleIntegrator::cleanup() {
     kernel = Kernel();
 }

openmmapi/src/MonteCarloAnisotropicBarostat.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,10 +35,31 @@
 using namespace OpenMM;
 
 MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double defaultTemperature, bool scaleX, bool scaleY, bool scaleZ, int frequency) :
-        defaultPressure(defaultPressure), defaultTemperature(defaultTemperature), scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ), frequency(frequency) {
+        scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ) {
+    setDefaultPressure(defaultPressure);
+    setDefaultTemperature(defaultTemperature);
+    setFrequency(frequency);
     setRandomNumberSeed(0);
 }
 
+void MonteCarloAnisotropicBarostat::setDefaultPressure(const Vec3& pressure) {
+    if (pressure[0] < 0 || pressure[1] < 0 || pressure[2] < 0)
+        throw OpenMMException("All components of pressure must be non-negative");
+    defaultPressure = pressure;
+}
+
+void MonteCarloAnisotropicBarostat::setFrequency(int freq) {
+    if (freq <= 0)
+        throw OpenMMException("Frequency must be positive");
+    frequency = freq;
+}
+
+void MonteCarloAnisotropicBarostat::setDefaultTemperature(double temp) {
+    if (temp < 0)
+        throw OpenMMException("Temperature cannot be negative");
+    defaultTemperature = temp;
+}
+
 ForceImpl* MonteCarloAnisotropicBarostat::createImpl() const {
     return new MonteCarloAnisotropicBarostatImpl(*this);
 }

openmmapi/src/MonteCarloBarostat.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -30,15 +30,37 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/MonteCarloBarostat.h"
+#include "openmm/OpenMMException.h"
 #include "openmm/internal/MonteCarloBarostatImpl.h"
+#include "openmm/MonteCarloAnisotropicBarostat.h"
 
 using namespace OpenMM;
 
-MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double defaultTemperature, int frequency) :
-        defaultPressure(defaultPressure), defaultTemperature(defaultTemperature), frequency(frequency) {
+MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double defaultTemperature, int frequency) {
+    setDefaultPressure(defaultPressure);
+    setDefaultTemperature(defaultTemperature);
+    setFrequency(frequency);
     setRandomNumberSeed(0);
 }
 
+void MonteCarloBarostat::setDefaultPressure(double pressure) {
+    if (pressure < 0)
+        throw OpenMMException("Pressure cannot be negative");
+    defaultPressure = pressure;
+}
+
+void MonteCarloBarostat::setFrequency(int freq) {
+    if (freq <= 0)
+        throw OpenMMException("Frequency must be positive");
+    frequency = freq;
+}
+
+void MonteCarloBarostat::setDefaultTemperature(double temp) {
+    if (temp < 0)
+        throw OpenMMException("Temperature cannot be negative");
+    defaultTemperature = temp;
+}
+
 ForceImpl* MonteCarloBarostat::createImpl() const {
     return new MonteCarloBarostatImpl(*this);
 }

openmmapi/src/MonteCarloFlexibleBarostat.cpp

@@ -34,11 +34,32 @@
 
 using namespace OpenMM;
 
-MonteCarloFlexibleBarostat::MonteCarloFlexibleBarostat(double defaultPressure, double defaultTemperature, int frequency, bool scaleMoleculesAsRigid) :
-        defaultPressure(defaultPressure), defaultTemperature(defaultTemperature), frequency(frequency), scaleMoleculesAsRigid(scaleMoleculesAsRigid) {
+MonteCarloFlexibleBarostat::MonteCarloFlexibleBarostat(double defaultPressure, double defaultTemperature,
+                                                       int frequency, bool scaleMoleculesAsRigid) : scaleMoleculesAsRigid(scaleMoleculesAsRigid) {
+    setDefaultPressure(defaultPressure);
+    setDefaultTemperature(defaultTemperature);
+    setFrequency(frequency);
     setRandomNumberSeed(0);
 }
 
+void MonteCarloFlexibleBarostat::setDefaultPressure(double pressure) {
+    if (pressure < 0)
+        throw OpenMMException("Pressure cannot be negative");
+    defaultPressure = pressure;
+}
+
+void MonteCarloFlexibleBarostat::setFrequency(int freq) {
+    if (freq <= 0)
+        throw OpenMMException("Frequency must be positive");
+    frequency = freq;
+}
+
+void MonteCarloFlexibleBarostat::setDefaultTemperature(double temp) {
+    if (temp < 0)
+        throw OpenMMException("Temperature cannot be negative");
+    defaultTemperature = temp;
+}
+
 ForceImpl* MonteCarloFlexibleBarostat::createImpl() const {
     return new MonteCarloFlexibleBarostatImpl(*this);
 }

openmmapi/src/MonteCarloMembraneBarostat.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,11 +35,38 @@
 using namespace OpenMM;
 
 MonteCarloMembraneBarostat::MonteCarloMembraneBarostat(double defaultPressure, double defaultSurfaceTension, double defaultTemperature, XYMode xymode, ZMode zmode, int frequency) :
-        defaultPressure(defaultPressure), defaultSurfaceTension(defaultSurfaceTension), defaultTemperature(defaultTemperature),
-        xymode(xymode), zmode(zmode), frequency(frequency) {
+        xymode(xymode), zmode(zmode) {
+    setDefaultPressure(defaultPressure);
+    setDefaultSurfaceTension(defaultSurfaceTension);
+    setDefaultTemperature(defaultTemperature);
+    setFrequency(frequency);
     setRandomNumberSeed(0);
 }
 
+void MonteCarloMembraneBarostat::setDefaultPressure(double pressure) {
+    if (pressure < 0)
+        throw OpenMMException("Pressure cannot be negative");
+    defaultPressure = pressure;
+}
+
+void MonteCarloMembraneBarostat::setDefaultSurfaceTension(double surfaceTension) {
+    if (surfaceTension < 0)
+        throw OpenMMException("Surface tension cannot be negative");
+    defaultSurfaceTension = surfaceTension;
+}
+
+void MonteCarloMembraneBarostat::setFrequency(int freq) {
+    if (freq <= 0)
+        throw OpenMMException("Frequency must be positive");
+    frequency = freq;
+}
+
+void MonteCarloMembraneBarostat::setDefaultTemperature(double temp) {
+    if (temp < 0)
+        throw OpenMMException("Temperature cannot be negative");
+    defaultTemperature = temp;
+}
+
 ForceImpl* MonteCarloMembraneBarostat::createImpl() const {
     return new MonteCarloMembraneBarostatImpl(*this);
 }

openmmapi/src/NonbondedForceImpl.cpp

@@ -62,33 +62,73 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
         if (owner.getSwitchingDistance() < 0 || owner.getSwitchingDistance() >= owner.getCutoffDistance())
             throw OpenMMException("NonbondedForce: Switching distance must satisfy 0 <= r_switch < r_cutoff");
     }
+    for (int i = 0; i < owner.getNumParticles(); i++) {
+        double charge, sigma, epsilon;
+        owner.getParticleParameters(i, charge, sigma, epsilon);
+        if (sigma < 0)
+            throw OpenMMException("NonbondedForce: sigma for a particle cannot be negative");
+        if (epsilon < 0)
+            throw OpenMMException("NonbondedForce: epsilon for a particle cannot be negative");
+    }
     vector<set<int> > exceptions(owner.getNumParticles());
     for (int i = 0; i < owner.getNumExceptions(); i++) {
-        int particle1, particle2;
+        int particle[2];
         double chargeProd, sigma, epsilon;
-        owner.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
-        if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
+        owner.getExceptionParameters(i, particle[0], particle[1], chargeProd, sigma, epsilon);
+        for (int j = 0; j < 2; j++) {
+            if (particle[j] < 0 || particle[j] >= owner.getNumParticles()) {
                 stringstream msg;
                 msg << "NonbondedForce: Illegal particle index for an exception: ";
-            msg << particle1;
+                msg << particle[j];
                 throw OpenMMException(msg.str());
             }
-        if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
-            stringstream msg;
-            msg << "NonbondedForce: Illegal particle index for an exception: ";
-            msg << particle2;
-            throw OpenMMException(msg.str());
         }
-        if (exceptions[particle1].count(particle2) > 0 || exceptions[particle2].count(particle1) > 0) {
+        if (exceptions[particle[0]].count(particle[1]) > 0 || exceptions[particle[1]].count(particle[0]) > 0) {
             stringstream msg;
             msg << "NonbondedForce: Multiple exceptions are specified for particles ";
-            msg << particle1;
+            msg << particle[0];
             msg << " and ";
-            msg << particle2;
+            msg << particle[1];
+            throw OpenMMException(msg.str());
+        }
+        exceptions[particle[0]].insert(particle[1]);
+        exceptions[particle[1]].insert(particle[0]);
+        if (sigma <= 0)
+            throw OpenMMException("NonbondedForce: sigma for an exception must be positive");
+        if (epsilon < 0)
+            throw OpenMMException("NonbondedForce: epsilon for an exception cannot be negative");
+    }
+    for (int i = 0; i < owner.getNumParticleParameterOffsets(); i++) {
+        string parameter;
+        int particleIndex;
+        double chargeScale, sigmaScale, epsilonScale;
+        owner.getParticleParameterOffset(i, parameter, particleIndex, chargeScale, sigmaScale, epsilonScale);
+        if (particleIndex < 0 || particleIndex >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "NonbondedForce: Illegal particle index for a particle parameter offset: ";
+            msg << particleIndex;
+            throw OpenMMException(msg.str());
+        }
+        if (sigmaScale < 0)
+            throw OpenMMException("NonbondedForce: sigma scale for a particle parameter offset cannot be negative");
+        if (epsilonScale < 0)
+            throw OpenMMException("NonbondedForce: epsilon scale for a particle parameter offset cannot be negative");
+    }
+    for (int i = 0; i < owner.getNumExceptionParameterOffsets(); i++) {
+        string parameter;
+        int exceptionIndex;
+        double chargeScale, sigmaScale, epsilonScale;
+        owner.getExceptionParameterOffset(i, parameter, exceptionIndex, chargeScale, sigmaScale, epsilonScale);
+        if (exceptionIndex < 0 || exceptionIndex >= owner.getNumExceptions()) {
+            stringstream msg;
+            msg << "NonbondedForce: Illegal exception index for an exception parameter offset: ";
+            msg << exceptionIndex;
             throw OpenMMException(msg.str());
         }
-        exceptions[particle1].insert(particle2);
-        exceptions[particle2].insert(particle1);
+        if (sigmaScale < 0)
+            throw OpenMMException("NonbondedForce: sigma scale for a particle parameter offset cannot be negative");
+        if (epsilonScale < 0)
+            throw OpenMMException("NonbondedForce: epsilon scale for a particle parameter offset cannot be negative");
     }
     if (owner.getNonbondedMethod() != NonbondedForce::NoCutoff && owner.getNonbondedMethod() != NonbondedForce::CutoffNonPeriodic) {
         Vec3 boxVectors[3];

openmmapi/src/NoseHooverIntegrator.cpp

@@ -70,6 +70,16 @@ NoseHooverIntegrator::~NoseHooverIntegrator() {}
 
 int NoseHooverIntegrator::addThermostat(double temperature, double collisionFrequency,
                                         int chainLength, int numMTS, int numYoshidaSuzuki) {
+    if (temperature < 0)
+        throw OpenMMException("NoseHooverIntegrator: temperature cannot be negative");
+    if (collisionFrequency <= 0)
+        throw OpenMMException("NoseHooverIntegrator: collisionFrequency must be positive");
+    if (collisionFrequency < 0)
+        throw OpenMMException("NoseHooverIntegrator: chainLength must be positive");
+    if (numMTS < 0)
+        throw OpenMMException("NoseHooverIntegrator: numMTS must be positive");
+    if (numYoshidaSuzuki != 1 && numYoshidaSuzuki != 3 && numYoshidaSuzuki != 5 && numYoshidaSuzuki != 7)
+        throw OpenMMException("NoseHooverIntegrator: numYoshidaSuzuki must be 1, 3, 5, or 7");
     hasSubsystemThermostats_ = false;
     return addSubsystemThermostat(std::vector<int>(), std::vector<std::pair<int, int>>(), temperature,
                                   collisionFrequency, temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
@@ -80,6 +90,20 @@ int NoseHooverIntegrator::addSubsystemThermostat(const std::vector<int>& thermos
                                                  double temperature, double collisionFrequency,
                                                  double relativeTemperature, double relativeCollisionFrequency,
                                                  int chainLength, int numMTS, int numYoshidaSuzuki) {
+    if (temperature < 0)
+        throw OpenMMException("NoseHooverIntegrator: temperature cannot be negative");
+    if (relativeTemperature < 0)
+        throw OpenMMException("NoseHooverIntegrator: relativeTemperature cannot be negative");
+    if (collisionFrequency <= 0)
+        throw OpenMMException("NoseHooverIntegrator: collisionFrequency must be positive");
+    if (relativeCollisionFrequency <= 0)
+        throw OpenMMException("NoseHooverIntegrator: relativeCollisionFrequency must be positive");
+    if (collisionFrequency < 0)
+        throw OpenMMException("NoseHooverIntegrator: chainLength must be positive");
+    if (numMTS < 0)
+        throw OpenMMException("NoseHooverIntegrator: numMTS must be positive");
+    if (numYoshidaSuzuki != 1 && numYoshidaSuzuki != 3 && numYoshidaSuzuki != 5 && numYoshidaSuzuki != 7)
+        throw OpenMMException("NoseHooverIntegrator: numYoshidaSuzuki must be 1, 3, 5, or 7");
     int chainID = noseHooverChains.size();
     // check if one thermostat already applies to all atoms or pairs
     if ( (chainID > 0) && (noseHooverChains[0].getThermostatedAtoms().size()+noseHooverChains[0].getThermostatedPairs().size() == 0) ) {

openmmapi/src/PeriodicTorsionForceImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -29,14 +29,14 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/PeriodicTorsionForceImpl.h"
 #include "openmm/kernels.h"
+#include <sstream>
 
 using namespace OpenMM;
-using std::pair;
-using std::vector;
-using std::set;
+using namespace std;
 
 PeriodicTorsionForceImpl::PeriodicTorsionForceImpl(const PeriodicTorsionForce& owner) : owner(owner) {
 }
@@ -45,6 +45,22 @@ PeriodicTorsionForceImpl::~PeriodicTorsionForceImpl() {
 }
 
 void PeriodicTorsionForceImpl::initialize(ContextImpl& context) {
+    const System& system = context.getSystem();
+    for (int i = 0; i < owner.getNumTorsions(); i++) {
+        int particle[4], periodicity;
+        double phase, k;
+        owner.getTorsionParameters(i, particle[0], particle[1], particle[2], particle[3], periodicity, phase, k);
+        for (int j = 0; j < 4; j++) {
+            if (particle[j] < 0 || particle[j] >= system.getNumParticles()) {
+                stringstream msg;
+                msg << "PeriodicTorsionForce: Illegal particle index for a torsion: ";
+                msg << particle[j];
+                throw OpenMMException(msg.str());
+            }
+        }
+        if (periodicity < 1)
+            throw OpenMMException("PeriodicTorsionForce: periodicity must be positive");
+    }
     kernel = context.getPlatform().createKernel(CalcPeriodicTorsionForceKernel::Name(), context);
     kernel.getAs<CalcPeriodicTorsionForceKernel>().initialize(context.getSystem(), owner);
 }

openmmapi/src/RBTorsionForceImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -29,14 +29,14 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/RBTorsionForceImpl.h"
 #include "openmm/kernels.h"
+#include <sstream>
 
 using namespace OpenMM;
-using std::pair;
-using std::vector;
-using std::set;
+using namespace std;
 
 RBTorsionForceImpl::RBTorsionForceImpl(const RBTorsionForce& owner) : owner(owner) {
 }
@@ -45,6 +45,20 @@ RBTorsionForceImpl::~RBTorsionForceImpl() {
 }
 
 void RBTorsionForceImpl::initialize(ContextImpl& context) {
+    const System& system = context.getSystem();
+    for (int i = 0; i < owner.getNumTorsions(); i++) {
+        int particle[4];
+        double c0, c1, c2, c3, c4, c5;
+        owner.getTorsionParameters(i, particle[0], particle[1], particle[2], particle[3], c0, c1, c2, c3, c4, c5);
+        for (int j = 0; j < 4; j++) {
+            if (particle[j] < 0 || particle[j] >= system.getNumParticles()) {
+                stringstream msg;
+                msg << "RBTorsionForce: Illegal particle index for a torsion: ";
+                msg << particle[j];
+                throw OpenMMException(msg.str());
+            }
+        }
+    }
     kernel = context.getPlatform().createKernel(CalcRBTorsionForceKernel::Name(), context);
     kernel.getAs<CalcRBTorsionForceKernel>().initialize(context.getSystem(), owner);
 }

openmmapi/src/VariableLangevinIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -61,6 +61,30 @@ void VariableLangevinIntegrator::initialize(ContextImpl& contextRef) {
     kernel.getAs<IntegrateVariableLangevinStepKernel>().initialize(contextRef.getSystem(), *this);
 }
 
+void VariableLangevinIntegrator::setTemperature(double temp) {
+    if (temp < 0)
+        throw OpenMMException("Temperature cannot be negative");
+    temperature = temp;
+}
+
+void VariableLangevinIntegrator::setFriction(double coeff) {
+    if (coeff < 0)
+        throw OpenMMException("Friction cannot be negative");
+    friction = coeff;
+}
+
+void VariableLangevinIntegrator::setErrorTolerance(double tol) {
+    if (tol <= 0)
+        throw OpenMMException("Error tolerance must be positive");
+    errorTol = tol;
+}
+
+void VariableLangevinIntegrator::setMaximumStepSize(double size) {
+    if (size < 0)
+        throw OpenMMException("Maximum step size cannot be negative");
+    maxStepSize = size;
+}
+
 void VariableLangevinIntegrator::cleanup() {
     kernel = Kernel();
 }

openmmapi/src/VariableVerletIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -55,6 +55,18 @@ void VariableVerletIntegrator::initialize(ContextImpl& contextRef) {
     kernel.getAs<IntegrateVariableVerletStepKernel>().initialize(contextRef.getSystem(), *this);
 }
 
+void VariableVerletIntegrator::setErrorTolerance(double tol) {
+    if (tol <= 0)
+        throw OpenMMException("Error tolerance must be positive");
+    errorTol = tol;
+}
+
+void VariableVerletIntegrator::setMaximumStepSize(double size) {
+    if (size < 0)
+        throw OpenMMException("Maximum step size cannot be negative");
+    maxStepSize = size;
+}
+
 void VariableVerletIntegrator::cleanup() {
     kernel = Kernel();
 }

plugins/drude/openmmapi/include/openmm/DrudeIntegrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -69,9 +69,7 @@ public:
      *
      * @param temp    the temperature of the heat bath, measured in Kelvin
      */
-    void setDrudeTemperature(double temp) {
-        drudeTemperature = temp;
-    }
+    void setDrudeTemperature(double temp);
     /**
      * Get the maximum distance a Drude particle can ever move from its parent particle, measured in nm.  This is implemented
      * with a hard wall constraint.  If this distance is set to 0 (the default), the hard wall constraint is omitted.

plugins/drude/openmmapi/include/openmm/DrudeLangevinIntegrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -78,9 +78,7 @@ public:
      *
      * @param temp    the temperature of the heat bath, measured in Kelvin
      */
-    void setTemperature(double temp) {
-        temperature = temp;
-    }
+    void setTemperature(double temp);
     /**
      * Get the friction coefficient which determines how strongly the system is coupled to
      * the main heat bath (in inverse ps).
@@ -96,9 +94,7 @@ public:
      *
      * @param coeff    the friction coefficient, measured in 1/ps
      */
-    void setFriction(double coeff) {
-        friction = coeff;
-    }
+    void setFriction(double coeff);
     /**
      * Get the friction coefficient which determines how strongly the internal coordinates of Drude particles
      * are coupled to the heat bath (in inverse ps).
@@ -114,9 +110,7 @@ public:
      *
      * @param coeff    the friction coefficient, measured in 1/ps
      */
-    void setDrudeFriction(double coeff) {
-        drudeFriction = coeff;
-    }
+    void setDrudeFriction(double coeff);
     /**
      * Advance a simulation through time by taking a series of time steps.
      *

plugins/drude/openmmapi/include/openmm/DrudeSCFIntegrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -70,9 +70,7 @@ public:
      *
      * @param tol    the error tolerance to use, measured in kJ/mol/nm
      */
-    void setMinimizationErrorTolerance(double tol) {
-        tolerance = tol;
-    }
+    void setMinimizationErrorTolerance(double tol);
     /**
      * Advance a simulation through time by taking a series of time steps.
      *

plugins/drude/openmmapi/src/DrudeForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Portions copyright (c) 2013-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -42,6 +42,10 @@ DrudeForce::DrudeForce() {
 }
 
 int DrudeForce::addParticle(int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34) {
+    if (polarizability <= 0)
+        throw OpenMMException("Polarizability must be positive");
+    if ((aniso12 <= 0 && particle2 != -1) || (aniso34 <= 0 && particle3 != -1 && particle4 != -1))
+        throw OpenMMException("Anisotropy factors must be positive");
     particles.push_back(ParticleInfo(particle, particle1, particle2, particle3, particle4, charge, polarizability, aniso12, aniso34));
     return particles.size()-1;
 }
@@ -61,6 +65,10 @@ void DrudeForce::getParticleParameters(int index, int& particle, int& particle1,
 
 void DrudeForce::setParticleParameters(int index, int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34) {
     ASSERT_VALID_INDEX(index, particles);
+    if (polarizability <= 0)
+        throw OpenMMException("Polarizability must be positive");
+    if ((aniso12 <= 0 && particle2 != -1) || (aniso34 <= 0 && particle3 != -1 && particle4 != -1))
+        throw OpenMMException("Anisotropy factors must be positive");
     particles[index].particle = particle;
     particles[index].particle1 = particle1;
     particles[index].particle2 = particle2;

plugins/drude/openmmapi/src/DrudeForceImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Portions copyright (c) 2013-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -58,84 +58,59 @@ void DrudeForceImpl::initialize(ContextImpl& context) {
 
     set<int> usedParticles;
     for (int i = 0; i < owner.getNumParticles(); i++) {
-        int particle, particle1, particle2, particle3, particle4;
+        int particle[5];
         double charge, k, k2, k3;
-        owner.getParticleParameters(i, particle, particle1, particle2, particle3, particle4, charge, k, k2, k3);
-        if (particle < 0 || particle >= system.getNumParticles()) {
+        owner.getParticleParameters(i, particle[0], particle[1], particle[2], particle[3], particle[4], charge, k, k2, k3);
+        for (int i = 0; i < 2; i++) {
+            if (particle[i] < 0 || particle[i] >= system.getNumParticles()) {
                 stringstream msg;
                 msg << "DrudeForce: Illegal particle index: ";
-            msg << particle;
+                msg << particle[i];
                 throw OpenMMException(msg.str());
             }
-        if (particle1 < 0 || particle1 >= system.getNumParticles()) {
+            if (usedParticles.find(particle[i]) != usedParticles.end()) {
                 stringstream msg;
-            msg << "DrudeForce: Illegal particle index: ";
-            msg << particle1;
+                msg << "DrudeForce: Particle index is used by two different Drude particles: ";
+                msg << particle[i];
                 throw OpenMMException(msg.str());
             }
-        if (particle2 < -1 || particle2 >= system.getNumParticles()) {
-            stringstream msg;
-            msg << "DrudeForce: Illegal particle index: ";
-            msg << particle2;
-            throw OpenMMException(msg.str());
+            usedParticles.insert(particle[i]);
         }
-        if (particle3 < -1 || particle3 >= system.getNumParticles()) {
+        for (int i = 2; i < 5; i++) {
+            if (particle[i] < -1 || particle[i] >= system.getNumParticles()) {
                 stringstream msg;
                 msg << "DrudeForce: Illegal particle index: ";
-            msg << particle3;
+                msg << particle[i];
                 throw OpenMMException(msg.str());
             }
-        if (particle4 < -1 || particle4 >= system.getNumParticles()) {
-            stringstream msg;
-            msg << "DrudeForce: Illegal particle index: ";
-            msg << particle4;
-            throw OpenMMException(msg.str());
         }
-        if (usedParticles.find(particle) != usedParticles.end()) {
-            stringstream msg;
-            msg << "DrudeForce: Particle index is used by two different Drude particles: ";
-            msg << particle;
-            throw OpenMMException(msg.str());
-        }
-        usedParticles.insert(particle);
-        if (usedParticles.find(particle1) != usedParticles.end()) {
-            stringstream msg;
-            msg << "DrudeForce: Particle index is used by two different Drude particles: ";
-            msg << particle1;
-            throw OpenMMException(msg.str());
-        }
-        usedParticles.insert(particle1);
     }
 
     // Check for errors in the specification of screened pairs.
 
     vector<set<int> > screenedPairs(owner.getNumParticles());
     for (int i = 0; i < owner.getNumScreenedPairs(); i++) {
-        int particle1, particle2;
+        int particle[2];
         double thole;
-        owner.getScreenedPairParameters(i, particle1, particle2, thole);
-        if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
+        owner.getScreenedPairParameters(i, particle[0], particle[1], thole);
+        for (int i = 0; i < 2; i++) {
+            if (particle[i] < 0 || particle[i] >= owner.getNumParticles()) {
                 stringstream msg;
                 msg << "DrudeForce: Illegal particle index for a screened pair: ";
-            msg << particle1;
+                msg << particle[i];
                 throw OpenMMException(msg.str());
             }
-        if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
-            stringstream msg;
-            msg << "DrudeForce: Illegal particle index for a screened pair: ";
-            msg << particle2;
-            throw OpenMMException(msg.str());
         }
-        if (screenedPairs[particle1].count(particle2) > 0 || screenedPairs[particle2].count(particle1) > 0) {
+        if (screenedPairs[particle[0]].count(particle[1]) > 0 || screenedPairs[particle[1]].count(particle[0]) > 0) {
             stringstream msg;
             msg << "DrudeForce: Multiple screened pairs are specified for particles ";
-            msg << particle1;
+            msg << particle[0];
             msg << " and ";
-            msg << particle2;
+            msg << particle[1];
             throw OpenMMException(msg.str());
         }
-        screenedPairs[particle1].insert(particle2);
-        screenedPairs[particle2].insert(particle1);
+        screenedPairs[particle[0]].insert(particle[1]);
+        screenedPairs[particle[1]].insert(particle[0]);
     }
     kernel.getAs<CalcDrudeForceKernel>().initialize(context.getSystem(), owner);
 }

plugins/drude/openmmapi/src/DrudeIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -47,6 +47,12 @@ std::vector<Vec3> DrudeIntegrator::getVelocitiesForTemperature(const System &sys
     return assignDrudeVelocities(system, temperature, drudeTemperature, randomSeedIn);
 }
 
+void DrudeIntegrator::setDrudeTemperature(double temp) {
+    if (temp < 0)
+        throw OpenMMException("Drude temperature cannot be negative");
+    drudeTemperature = temp;
+}
+
 double DrudeIntegrator::getMaxDrudeDistance() const {
     return maxDrudeDistance;
 }