Added index bounds checking to lots of methods

openmmapi/include/openmm/System.h

@@ -97,9 +97,7 @@ public:
      * 
      * @param index the index of the particle for which to get the mass
      */
-    double getParticleMass(int index) const {
-        return masses[index];
-    }
+    double getParticleMass(int index) const;
     /**
      * Set the mass (in atomic mass units) of a particle.  If the mass is 0, Integrators will ignore
      * the particle and not modify its position or velocity.  This is most often
@@ -109,9 +107,7 @@ public:
      * @param index  the index of the particle for which to set the mass
      * @param mass   the mass of the particle
      */
-    void setParticleMass(int index, double mass) {
-        masses[index] = mass;
-    }
+    void setParticleMass(int index, double mass);
     /**
      * Set a particle to be a virtual site.  The VirtualSite object should have
      * been created on the heap with the "new" operator.  The System takes over
@@ -195,17 +191,13 @@ public:
      * 
      * @param index  the index of the Force to get
      */
-    const Force& getForce(int index) const {
-        return *forces[index];
-    }
+    const Force& getForce(int index) const;
     /**
      * Get a writable reference to one of the Forces in this System.
      *
      * @param index  the index of the Force to get
      */
-    Force& getForce(int index) {
-        return *forces[index];
-    }
+    Force& getForce(int index);
     /**
      * Get the default values of the vectors defining the axes of the periodic box (measured in nm).  Any newly
      * created Context will have its box vectors set to these.  They will affect

tests/AssertionUtilities.h → openmmapi/include/openmm/internal/AssertionUtilities.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -36,23 +36,15 @@
  * This file provides a variety of macros useful in test cases.
  */
 
-#include "openmm/OpenMMException.h"
 #include <cmath>
 #include <string>
 #include <sstream>
 
-void throwException(const char* file, int line, const std::string& details) {
-    std::string fn(file);
-    std::string::size_type pos = fn.find_last_of("/\\");
-    if (pos+1>=fn.size())
-        pos=0;
-    std::string filename(fn,(int)(pos+1),(int)(fn.size()-(pos+1)));
-    std::stringstream message;
-    message << "Assertion failure at "<<filename<<":"<<line;
-    if (details.size() > 0)
-        message << ".  "<<details;
-    throw OpenMM::OpenMMException(message.str());
-}
+namespace OpenMM {
+
+void throwException(const char* file, int line, const std::string& details);
+
+} // namespace OpenMM
 
 #define ASSERT(cond) {if (!(cond)) throwException(__FILE__, __LINE__, "");};
 
@@ -64,4 +56,6 @@ void throwException(const char* file, int line, const std::string& details) {
 
 #define ASSERT_USUALLY_EQUAL_TOL(expected, found, tol) {double _scale_ = std::abs(expected) > 1.0 ? std::abs(expected) : 1.0; if (!(std::abs((expected)-(found))/_scale_ <= (tol))) {std::stringstream details; details << "Expected "<<(expected)<<", found "<<(found)<<" (This test is stochastic and may occasionally fail)"; throwException(__FILE__, __LINE__, details.str());}};
 
+#define ASSERT_VALID_INDEX(index, vector) {if (index < 0 || index >= vector.size()) throwException(__FILE__, __LINE__, "Index out of range");};
+
 #endif /*OPENMM_ASSERTIONUTILITIES_H_*/

openmmapi/src/AssertionUtilities.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This file provides a variety of macros useful in test cases.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/OpenMMException.h"
+#include <string>
+#include <sstream>
+
+namespace OpenMM {
+
+void throwException(const char* file, int line, const std::string& details) {
+    std::string fn(file);
+    std::string::size_type pos = fn.find_last_of("/\\");
+    if (pos+1>=fn.size())
+        pos=0;
+    std::string filename(fn,(int)(pos+1),(int)(fn.size()-(pos+1)));
+    std::stringstream message;
+    message << "Assertion failure at "<<filename<<":"<<line;
+    if (details.size() > 0)
+        message << ".  "<<details;
+    throw OpenMMException(message.str());
+}
+
+} // namespace OpenMM

openmmapi/src/CMAPTorsionForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/CMAPTorsionForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/CMAPTorsionForceImpl.h"
 
 using namespace OpenMM;
@@ -47,11 +48,13 @@ int CMAPTorsionForce::addMap(int size, const std::vector<double>& energy) {
 }
 
 void CMAPTorsionForce::getMapParameters(int index, int& size, std::vector<double>& energy) const {
+    ASSERT_VALID_INDEX(index, maps);
     size = maps[index].size;
     energy = maps[index].energy;
 }
 
 void CMAPTorsionForce::setMapParameters(int index, int size, const std::vector<double>& energy) {
+    ASSERT_VALID_INDEX(index, maps);
     if (energy.size() != size*size)
         throw OpenMMException("CMAPTorsionForce: incorrect number of energy values");
     maps[index].size = size;
@@ -64,6 +67,7 @@ int CMAPTorsionForce::addTorsion(int map, int a1, int a2, int a3, int a4, int b1
 }
 
 void CMAPTorsionForce::getTorsionParameters(int index, int& map, int& a1, int& a2, int& a3, int& a4, int& b1, int& b2, int& b3, int& b4) const {
+    ASSERT_VALID_INDEX(index, torsions);
     map = torsions[index].map;
     a1 = torsions[index].a1;
     a2 = torsions[index].a2;
@@ -76,6 +80,7 @@ void CMAPTorsionForce::getTorsionParameters(int index, int& map, int& a1, int& a
 }
 
 void CMAPTorsionForce::setTorsionParameters(int index, int map, int a1, int a2, int a3, int a4, int b1, int b2, int b3, int b4) {
+    ASSERT_VALID_INDEX(index, torsions);
     torsions[index].map = map;
     torsions[index].a1 = a1;
     torsions[index].a2 = a2;

openmmapi/src/CustomAngleForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/CustomAngleForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/CustomAngleForceImpl.h"
 #include <cmath>
 #include <map>
@@ -60,10 +61,12 @@ int CustomAngleForce::addPerAngleParameter(const string& name) {
 }
 
 const string& CustomAngleForce::getPerAngleParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, parameters);
     return parameters[index].name;
 }
 
 void CustomAngleForce::setPerAngleParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, parameters);
     parameters[index].name = name;
 }
 
@@ -73,18 +76,22 @@ int CustomAngleForce::addGlobalParameter(const string& name, double defaultValue
 }
 
 const string& CustomAngleForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].name;
 }
 
 void CustomAngleForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].name = name;
 }
 
 double CustomAngleForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].defaultValue;
 }
 
 void CustomAngleForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].defaultValue = defaultValue;
 }
 
@@ -94,6 +101,7 @@ int CustomAngleForce::addAngle(int particle1, int particle2, int particle3, cons
 }
 
 void CustomAngleForce::getAngleParameters(int index, int& particle1, int& particle2, int& particle3, std::vector<double>& parameters) const {
+    ASSERT_VALID_INDEX(index, angles);
     particle1 = angles[index].particle1;
     particle2 = angles[index].particle2;
     particle3 = angles[index].particle3;
@@ -101,6 +109,7 @@ void CustomAngleForce::getAngleParameters(int index, int& particle1, int& partic
 }
 
 void CustomAngleForce::setAngleParameters(int index, int particle1, int particle2, int particle3, const vector<double>& parameters) {
+    ASSERT_VALID_INDEX(index, angles);
     angles[index].parameters = parameters;
     angles[index].particle1 = particle1;
     angles[index].particle2 = particle2;

openmmapi/src/CustomBondForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/CustomBondForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/CustomBondForceImpl.h"
 #include <cmath>
 #include <map>
@@ -60,10 +61,12 @@ int CustomBondForce::addPerBondParameter(const string& name) {
 }
 
 const string& CustomBondForce::getPerBondParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, parameters);
     return parameters[index].name;
 }
 
 void CustomBondForce::setPerBondParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, parameters);
     parameters[index].name = name;
 }
 
@@ -73,18 +76,22 @@ int CustomBondForce::addGlobalParameter(const string& name, double defaultValue)
 }
 
 const string& CustomBondForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].name;
 }
 
 void CustomBondForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].name = name;
 }
 
 double CustomBondForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].defaultValue;
 }
 
 void CustomBondForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].defaultValue = defaultValue;
 }
 
@@ -94,12 +101,14 @@ int CustomBondForce::addBond(int particle1, int particle2, const vector<double>&
 }
 
 void CustomBondForce::getBondParameters(int index, int& particle1, int& particle2, std::vector<double>& parameters) const {
+    ASSERT_VALID_INDEX(index, bonds);
     particle1 = bonds[index].particle1;
     particle2 = bonds[index].particle2;
     parameters = bonds[index].parameters;
 }
 
 void CustomBondForce::setBondParameters(int index, int particle1, int particle2, const vector<double>& parameters) {
+    ASSERT_VALID_INDEX(index, bonds);
     bonds[index].parameters = parameters;
     bonds[index].particle1 = particle1;
     bonds[index].particle2 = particle2;

openmmapi/src/CustomExternalForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/CustomExternalForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/CustomExternalForceImpl.h"
 #include <cmath>
 #include <map>
@@ -60,10 +61,12 @@ int CustomExternalForce::addPerParticleParameter(const string& name) {
 }
 
 const string& CustomExternalForce::getPerParticleParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, parameters);
     return parameters[index].name;
 }
 
 void CustomExternalForce::setPerParticleParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, parameters);
     parameters[index].name = name;
 }
 
@@ -73,18 +76,22 @@ int CustomExternalForce::addGlobalParameter(const string& name, double defaultVa
 }
 
 const string& CustomExternalForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].name;
 }
 
 void CustomExternalForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].name = name;
 }
 
 double CustomExternalForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].defaultValue;
 }
 
 void CustomExternalForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].defaultValue = defaultValue;
 }
 
@@ -94,11 +101,13 @@ int CustomExternalForce::addParticle(int particle, const vector<double>& paramet
 }
 
 void CustomExternalForce::getParticleParameters(int index, int& particle, std::vector<double>& parameters) const {
+    ASSERT_VALID_INDEX(index, particles);
     particle = particles[index].particle;
     parameters = particles[index].parameters;
 }
 
 void CustomExternalForce::setParticleParameters(int index, int particle, const vector<double>& parameters) {
+    ASSERT_VALID_INDEX(index, particles);
     particles[index].parameters = parameters;
     particles[index].particle = particle;
 }

openmmapi/src/CustomGBForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/CustomGBForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/CustomGBForceImpl.h"
 #include <cmath>
 #include <map>
@@ -71,10 +72,12 @@ int CustomGBForce::addPerParticleParameter(const string& name) {
 }
 
 const string& CustomGBForce::getPerParticleParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, parameters);
     return parameters[index].name;
 }
 
 void CustomGBForce::setPerParticleParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, parameters);
     parameters[index].name = name;
 }
 
@@ -84,18 +87,22 @@ int CustomGBForce::addGlobalParameter(const string& name, double defaultValue) {
 }
 
 const string& CustomGBForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].name;
 }
 
 void CustomGBForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].name = name;
 }
 
 double CustomGBForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].defaultValue;
 }
 
 void CustomGBForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].defaultValue = defaultValue;
 }
 
@@ -105,10 +112,12 @@ int CustomGBForce::addParticle(const vector<double>& parameters) {
 }
 
 void CustomGBForce::getParticleParameters(int index, std::vector<double>& parameters) const {
+    ASSERT_VALID_INDEX(index, particles);
     parameters = particles[index].parameters;
 }
 
 void CustomGBForce::setParticleParameters(int index, const vector<double>& parameters) {
+    ASSERT_VALID_INDEX(index, particles);
     particles[index].parameters = parameters;
 }
 
@@ -118,12 +127,14 @@ int CustomGBForce::addComputedValue(const std::string& name, const std::string&
 }
 
 void CustomGBForce::getComputedValueParameters(int index, std::string& name, std::string& expression, ComputationType& type) const {
+    ASSERT_VALID_INDEX(index, computedValues);
     name = computedValues[index].name;
     expression = computedValues[index].expression;
     type = computedValues[index].type;
 }
 
 void CustomGBForce::setComputedValueParameters(int index, const std::string& name, const std::string& expression, ComputationType type) {
+    ASSERT_VALID_INDEX(index, computedValues);
     computedValues[index].name = name;
     computedValues[index].expression = expression;
     computedValues[index].type = type;
@@ -135,11 +146,13 @@ int CustomGBForce::addEnergyTerm(const std::string& expression, ComputationType
 }
 
 void CustomGBForce::getEnergyTermParameters(int index, std::string& expression, ComputationType& type) const {
+    ASSERT_VALID_INDEX(index, energyTerms);
     expression = energyTerms[index].expression;
     type = energyTerms[index].type;
 }
 
 void CustomGBForce::setEnergyTermParameters(int index, const std::string& expression, ComputationType type) {
+    ASSERT_VALID_INDEX(index, energyTerms);
     energyTerms[index].expression = expression;
     energyTerms[index].type = type;
 }
@@ -149,11 +162,13 @@ int CustomGBForce::addExclusion(int particle1, int particle2) {
     return exclusions.size()-1;
 }
 void CustomGBForce::getExclusionParticles(int index, int& particle1, int& particle2) const {
+    ASSERT_VALID_INDEX(index, exclusions);
     particle1 = exclusions[index].particle1;
     particle2 = exclusions[index].particle2;
 }
 
 void CustomGBForce::setExclusionParticles(int index, int particle1, int particle2) {
+    ASSERT_VALID_INDEX(index, exclusions);
     exclusions[index].particle1 = particle1;
     exclusions[index].particle2 = particle2;
 }
@@ -168,6 +183,7 @@ int CustomGBForce::addFunction(const std::string& name, const std::vector<double
 }
 
 void CustomGBForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
+    ASSERT_VALID_INDEX(index, functions);
     name = functions[index].name;
     values = functions[index].values;
     min = functions[index].min;
@@ -179,6 +195,7 @@ void CustomGBForce::setFunctionParameters(int index, const std::string& name, co
         throw OpenMMException("CustomGBForce: max <= min for a tabulated function.");
     if (values.size() < 2)
         throw OpenMMException("CustomGBForce: a tabulated function must have at least two points");
+    ASSERT_VALID_INDEX(index, functions);
     functions[index].name = name;
     functions[index].values = values;
     functions[index].min = min;

openmmapi/src/CustomHbondForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/CustomHbondForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
 #include <cmath>
 #include <map>
@@ -79,10 +80,12 @@ int CustomHbondForce::addPerDonorParameter(const string& name) {
 }
 
 const string& CustomHbondForce::getPerDonorParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, donorParameters);
     return donorParameters[index].name;
 }
 
 void CustomHbondForce::setPerDonorParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, donorParameters);
     donorParameters[index].name = name;
 }
 
@@ -92,10 +95,12 @@ int CustomHbondForce::addPerAcceptorParameter(const string& name) {
 }
 
 const string& CustomHbondForce::getPerAcceptorParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, acceptorParameters);
     return acceptorParameters[index].name;
 }
 
 void CustomHbondForce::setPerAcceptorParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, acceptorParameters);
     acceptorParameters[index].name = name;
 }
 
@@ -105,18 +110,22 @@ int CustomHbondForce::addGlobalParameter(const string& name, double defaultValue
 }
 
 const string& CustomHbondForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].name;
 }
 
 void CustomHbondForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].name = name;
 }
 
 double CustomHbondForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].defaultValue;
 }
 
 void CustomHbondForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].defaultValue = defaultValue;
 }
 
@@ -126,6 +135,7 @@ int CustomHbondForce::addDonor(int d1, int d2, int d3, const vector<double>& par
 }
 
 void CustomHbondForce::getDonorParameters(int index, int& d1, int& d2, int&  d3, std::vector<double>& parameters) const {
+    ASSERT_VALID_INDEX(index, donors);
     d1 = donors[index].p1;
     d2 = donors[index].p2;
     d3 = donors[index].p3;
@@ -133,6 +143,7 @@ void CustomHbondForce::getDonorParameters(int index, int& d1, int& d2, int&  d3,
 }
 
 void CustomHbondForce::setDonorParameters(int index, int d1, int d2, int d3, const vector<double>& parameters) {
+    ASSERT_VALID_INDEX(index, donors);
     donors[index].p1 = d1;
     donors[index].p2 = d2;
     donors[index].p3 = d3;
@@ -145,6 +156,7 @@ int CustomHbondForce::addAcceptor(int a1, int a2, int a3, const vector<double>&
 }
 
 void CustomHbondForce::getAcceptorParameters(int index, int& a1, int& a2, int& a3, std::vector<double>& parameters) const {
+    ASSERT_VALID_INDEX(index, acceptors);
     a1 = acceptors[index].p1;
     a2 = acceptors[index].p2;
     a3 = acceptors[index].p3;
@@ -152,6 +164,7 @@ void CustomHbondForce::getAcceptorParameters(int index, int& a1, int& a2, int& a
 }
 
 void CustomHbondForce::setAcceptorParameters(int index, int a1, int a2, int a3, const vector<double>& parameters) {
+    ASSERT_VALID_INDEX(index, acceptors);
     acceptors[index].p1 = a1;
     acceptors[index].p2 = a2;
     acceptors[index].p3 = a3;
@@ -163,11 +176,13 @@ int CustomHbondForce::addExclusion(int donor, int acceptor) {
     return exclusions.size()-1;
 }
 void CustomHbondForce::getExclusionParticles(int index, int& donor, int& acceptor) const {
+    ASSERT_VALID_INDEX(index, exclusions);
     donor = exclusions[index].donor;
     acceptor = exclusions[index].acceptor;
 }
 
 void CustomHbondForce::setExclusionParticles(int index, int donor, int acceptor) {
+    ASSERT_VALID_INDEX(index, exclusions);
     exclusions[index].donor = donor;
     exclusions[index].acceptor = acceptor;
 }
@@ -182,6 +197,7 @@ int CustomHbondForce::addFunction(const std::string& name, const std::vector<dou
 }
 
 void CustomHbondForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
+    ASSERT_VALID_INDEX(index, functions);
     name = functions[index].name;
     values = functions[index].values;
     min = functions[index].min;
@@ -193,6 +209,7 @@ void CustomHbondForce::setFunctionParameters(int index, const std::string& name,
         throw OpenMMException("CustomHbondForce: max <= min for a tabulated function.");
     if (values.size() < 2)
         throw OpenMMException("CustomHbondForce: a tabulated function must have at least two points");
+    ASSERT_VALID_INDEX(index, functions);
     functions[index].name = name;
     functions[index].values = values;
     functions[index].min = min;

openmmapi/src/CustomIntegrator.cpp

@@ -32,6 +32,7 @@
 #include "openmm/CustomIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/kernels.h"
 #include <ctime>
@@ -88,8 +89,7 @@ int CustomIntegrator::addGlobalVariable(const string& name, double initialValue)
 }
 
 const string& CustomIntegrator::getGlobalVariableName(int index) const {
-    if (index < 0 || index >= globalNames.size())
-        throw OpenMMException("Index out of range");
+    ASSERT_VALID_INDEX(index, globalNames);
     return globalNames[index];
 }
 
@@ -102,14 +102,12 @@ int CustomIntegrator::addPerDofVariable(const string& name, double initialValue)
 }
 
 const string& CustomIntegrator::getPerDofVariableName(int index) const {
-    if (index < 0 || index >= perDofNames.size())
-        throw OpenMMException("Index out of range");
+    ASSERT_VALID_INDEX(index, perDofNames);
     return perDofNames[index];
 }
 
 double CustomIntegrator::getGlobalVariable(int index) const {
-    if (index < 0 || index >= globalNames.size())
-        throw OpenMMException("Index out of range");
+    ASSERT_VALID_INDEX(index, globalValues);
     if (owner != NULL && !globalsAreCurrent) {
         dynamic_cast<const IntegrateCustomStepKernel&>(kernel.getImpl()).getGlobalVariables(*context, globalValues);
         globalsAreCurrent = true;
@@ -118,8 +116,7 @@ double CustomIntegrator::getGlobalVariable(int index) const {
 }
 
 void CustomIntegrator::setGlobalVariable(int index, double value) {
-    if (index < 0 || index >= globalNames.size())
-        throw OpenMMException("Index out of range");
+    ASSERT_VALID_INDEX(index, globalValues);
     if (owner != NULL && !globalsAreCurrent) {
         dynamic_cast<IntegrateCustomStepKernel&>(kernel.getImpl()).getGlobalVariables(*context, globalValues);
         globalsAreCurrent = true;
@@ -138,8 +135,7 @@ void CustomIntegrator::setGlobalVariableByName(const string& name, double value)
 }
 
 void CustomIntegrator::getPerDofVariable(int index, vector<Vec3>& values) const {
-    if (index < 0 || index >= perDofNames.size())
-        throw OpenMMException("Index out of range");
+    ASSERT_VALID_INDEX(index, perDofValues);
     if (owner == NULL)
         values = perDofValues[index];
     else
@@ -147,8 +143,7 @@ void CustomIntegrator::getPerDofVariable(int index, vector<Vec3>& values) const
 }
 
 void CustomIntegrator::setPerDofVariable(int index, const vector<Vec3>& values) {
-    if (index < 0 || index >= perDofNames.size())
-        throw OpenMMException("Index out of range");
+    ASSERT_VALID_INDEX(index, perDofValues);
     if (values.size() != context->getSystem().getNumParticles())
         throw OpenMMException("setPerDofVariable() called with wrong number of values");
     if (owner == NULL)
@@ -209,8 +204,7 @@ int CustomIntegrator::addUpdateContextState() {
 }
 
 void CustomIntegrator::getComputationStep(int index, ComputationType& type, string& variable, string& expression) const {
-    if (index < 0 || index >= computations.size())
-        throw OpenMMException("Index out of range");
+    ASSERT_VALID_INDEX(index, computations);
     type = computations[index].type;
     variable = computations[index].variable;
     expression = computations[index].expression;

openmmapi/src/CustomNonbondedForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/CustomNonbondedForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/CustomNonbondedForceImpl.h"
 #include <cmath>
 #include <map>
@@ -79,10 +80,12 @@ int CustomNonbondedForce::addPerParticleParameter(const string& name) {
 }
 
 const string& CustomNonbondedForce::getPerParticleParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, parameters);
     return parameters[index].name;
 }
 
 void CustomNonbondedForce::setPerParticleParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, parameters);
     parameters[index].name = name;
 }
 
@@ -92,18 +95,22 @@ int CustomNonbondedForce::addGlobalParameter(const string& name, double defaultV
 }
 
 const string& CustomNonbondedForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].name;
 }
 
 void CustomNonbondedForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].name = name;
 }
 
 double CustomNonbondedForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].defaultValue;
 }
 
 void CustomNonbondedForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].defaultValue = defaultValue;
 }
 
@@ -113,10 +120,12 @@ int CustomNonbondedForce::addParticle(const vector<double>& parameters) {
 }
 
 void CustomNonbondedForce::getParticleParameters(int index, std::vector<double>& parameters) const {
+    ASSERT_VALID_INDEX(index, particles);
     parameters = particles[index].parameters;
 }
 
 void CustomNonbondedForce::setParticleParameters(int index, const vector<double>& parameters) {
+    ASSERT_VALID_INDEX(index, particles);
     particles[index].parameters = parameters;
 }
 
@@ -125,11 +134,13 @@ int CustomNonbondedForce::addExclusion(int particle1, int particle2) {
     return exclusions.size()-1;
 }
 void CustomNonbondedForce::getExclusionParticles(int index, int& particle1, int& particle2) const {
+    ASSERT_VALID_INDEX(index, exclusions);
     particle1 = exclusions[index].particle1;
     particle2 = exclusions[index].particle2;
 }
 
 void CustomNonbondedForce::setExclusionParticles(int index, int particle1, int particle2) {
+    ASSERT_VALID_INDEX(index, exclusions);
     exclusions[index].particle1 = particle1;
     exclusions[index].particle2 = particle2;
 }
@@ -144,6 +155,7 @@ int CustomNonbondedForce::addFunction(const std::string& name, const std::vector
 }
 
 void CustomNonbondedForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
+    ASSERT_VALID_INDEX(index, functions);
     name = functions[index].name;
     values = functions[index].values;
     min = functions[index].min;
@@ -155,6 +167,7 @@ void CustomNonbondedForce::setFunctionParameters(int index, const std::string& n
         throw OpenMMException("CustomNonbondedForce: max <= min for a tabulated function.");
     if (values.size() < 2)
         throw OpenMMException("CustomNonbondedForce: a tabulated function must have at least two points");
+    ASSERT_VALID_INDEX(index, functions);
     functions[index].name = name;
     functions[index].values = values;
     functions[index].min = min;

openmmapi/src/CustomTorsionForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/CustomTorsionForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/CustomTorsionForceImpl.h"
 #include <cmath>
 #include <map>
@@ -60,10 +61,12 @@ int CustomTorsionForce::addPerTorsionParameter(const string& name) {
 }
 
 const string& CustomTorsionForce::getPerTorsionParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, parameters);
     return parameters[index].name;
 }
 
 void CustomTorsionForce::setPerTorsionParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, parameters);
     parameters[index].name = name;
 }
 
@@ -73,18 +76,22 @@ int CustomTorsionForce::addGlobalParameter(const string& name, double defaultVal
 }
 
 const string& CustomTorsionForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].name;
 }
 
 void CustomTorsionForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].name = name;
 }
 
 double CustomTorsionForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
     return globalParameters[index].defaultValue;
 }
 
 void CustomTorsionForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
     globalParameters[index].defaultValue = defaultValue;
 }
 
@@ -94,6 +101,7 @@ int CustomTorsionForce::addTorsion(int particle1, int particle2, int particle3,
 }
 
 void CustomTorsionForce::getTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, std::vector<double>& parameters) const {
+    ASSERT_VALID_INDEX(index, torsions);
     particle1 = torsions[index].particle1;
     particle2 = torsions[index].particle2;
     particle3 = torsions[index].particle3;
@@ -102,6 +110,7 @@ void CustomTorsionForce::getTorsionParameters(int index, int& particle1, int& pa
 }
 
 void CustomTorsionForce::setTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4, const vector<double>& parameters) {
+    ASSERT_VALID_INDEX(index, torsions);
     torsions[index].parameters = parameters;
     torsions[index].particle1 = particle1;
     torsions[index].particle2 = particle2;

openmmapi/src/GBSAOBCForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/GBSAOBCForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/GBSAOBCForceImpl.h"
 
 using namespace OpenMM;
@@ -45,12 +46,14 @@ int GBSAOBCForce::addParticle(double charge, double radius, double scalingFactor
 }
 
 void GBSAOBCForce::getParticleParameters(int index, double& charge, double& radius, double& scalingFactor) const {
+    ASSERT_VALID_INDEX(index, particles);
     charge = particles[index].charge;
     radius = particles[index].radius;
     scalingFactor = particles[index].scalingFactor;
 }
 
 void GBSAOBCForce::setParticleParameters(int index, double charge, double radius, double scalingFactor) {
+    ASSERT_VALID_INDEX(index, particles);
     particles[index].charge = charge;
     particles[index].radius = radius;
     particles[index].scalingFactor = scalingFactor;

openmmapi/src/GBVIForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/GBVIForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/GBVIForceImpl.h"
 #include <sstream>
 
@@ -47,12 +48,14 @@ int GBVIForce::addParticle(double charge, double radius, double gamma) {
 }
 
 void GBVIForce::getParticleParameters(int index, double& charge, double& radius, double& gamma) const {
+    ASSERT_VALID_INDEX(index, particles);
     charge = particles[index].charge;
     radius = particles[index].radius;
     gamma  = particles[index].gamma;
 }
 
 void GBVIForce::setParticleParameters(int index, double charge, double radius, double gamma) {
+    ASSERT_VALID_INDEX(index, particles);
     particles[index].charge = charge;
     particles[index].radius = radius;
     particles[index].gamma  = gamma;
@@ -104,6 +107,7 @@ int GBVIForce::addBond(int particle1, int particle2, double bondLength) {
 }
 
 void GBVIForce::setBondParameters( int index, int particle1, int particle2, double bondLength) {
+    ASSERT_VALID_INDEX(index, bonds);
     bonds[index].particle1  = particle1;
     bonds[index].particle2  = particle2;
     bonds[index].bondLength = bondLength;
@@ -114,6 +118,7 @@ int GBVIForce::getNumBonds( void ) const {
 }
 
 void GBVIForce::getBondParameters(int index, int& bondIndex1, int& bondIndex2, double& bondLength) const {
+    ASSERT_VALID_INDEX(index, bonds);
     bondIndex1 = bonds[index].particle1;
     bondIndex2 = bonds[index].particle2;
     bondLength = bonds[index].bondLength;

openmmapi/src/HarmonicAngleForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/HarmonicAngleForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/HarmonicAngleForceImpl.h"
 
 using namespace OpenMM;
@@ -45,6 +46,7 @@ int HarmonicAngleForce::addAngle(int particle1, int particle2, int particle3, do
 }
 
 void HarmonicAngleForce::getAngleParameters(int index, int& particle1, int& particle2, int& particle3, double& angle, double& k) const {
+    ASSERT_VALID_INDEX(index, angles);
     particle1 = angles[index].particle1;
     particle2 = angles[index].particle2;
     particle3 = angles[index].particle3;
@@ -53,6 +55,7 @@ void HarmonicAngleForce::getAngleParameters(int index, int& particle1, int& part
 }
 
 void HarmonicAngleForce::setAngleParameters(int index, int particle1, int particle2, int particle3, double angle, double k) {
+    ASSERT_VALID_INDEX(index, angles);
     angles[index].particle1 = particle1;
     angles[index].particle2 = particle2;
     angles[index].particle3 = particle3;

openmmapi/src/HarmonicBondForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/HarmonicBondForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/HarmonicBondForceImpl.h"
 
 using namespace OpenMM;
@@ -45,6 +46,7 @@ int HarmonicBondForce::addBond(int particle1, int particle2, double length, doub
 }
 
 void HarmonicBondForce::getBondParameters(int index, int& particle1, int& particle2, double& length, double& k) const {
+    ASSERT_VALID_INDEX(index, bonds);
     particle1 = bonds[index].particle1;
     particle2 = bonds[index].particle2;
     length = bonds[index].length;
@@ -52,6 +54,7 @@ void HarmonicBondForce::getBondParameters(int index, int& particle1, int& partic
 }
 
 void HarmonicBondForce::setBondParameters(int index, int particle1, int particle2, double length, double k) {
+    ASSERT_VALID_INDEX(index, bonds);
     bonds[index].particle1 = particle1;
     bonds[index].particle2 = particle2;
     bonds[index].length = length;

openmmapi/src/NonbondedForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/NonbondedForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include <cmath>
 #include <map>
@@ -88,12 +89,14 @@ int NonbondedForce::addParticle(double charge, double sigma, double epsilon) {
 }
 
 void NonbondedForce::getParticleParameters(int index, double& charge, double& sigma, double& epsilon) const {
+    ASSERT_VALID_INDEX(index, particles);
     charge = particles[index].charge;
     sigma = particles[index].sigma;
     epsilon = particles[index].epsilon;
 }
 
 void NonbondedForce::setParticleParameters(int index, double charge, double sigma, double epsilon) {
+    ASSERT_VALID_INDEX(index, particles);
     particles[index].charge = charge;
     particles[index].sigma = sigma;
     particles[index].epsilon = epsilon;
@@ -125,6 +128,7 @@ int NonbondedForce::addException(int particle1, int particle2, double chargeProd
     return newIndex;
 }
 void NonbondedForce::getExceptionParameters(int index, int& particle1, int& particle2, double& chargeProd, double& sigma, double& epsilon) const {
+    ASSERT_VALID_INDEX(index, exceptions);
     particle1 = exceptions[index].particle1;
     particle2 = exceptions[index].particle2;
     chargeProd = exceptions[index].chargeProd;
@@ -133,6 +137,7 @@ void NonbondedForce::getExceptionParameters(int index, int& particle1, int& part
 }
 
 void NonbondedForce::setExceptionParameters(int index, int particle1, int particle2, double chargeProd, double sigma, double epsilon) {
+    ASSERT_VALID_INDEX(index, exceptions);
     exceptions[index].particle1 = particle1;
     exceptions[index].particle2 = particle2;
     exceptions[index].chargeProd = chargeProd;

openmmapi/src/PeriodicTorsionForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/PeriodicTorsionForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/PeriodicTorsionForceImpl.h"
 
 using namespace OpenMM;
@@ -45,6 +46,7 @@ int PeriodicTorsionForce::addTorsion(int particle1, int particle2, int particle3
 }
 
 void PeriodicTorsionForce::getTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, int& periodicity, double& phase, double& k) const {
+    ASSERT_VALID_INDEX(index, periodicTorsions);
     particle1 = periodicTorsions[index].particle1;
     particle2 = periodicTorsions[index].particle2;
     particle3 = periodicTorsions[index].particle3;
@@ -55,6 +57,7 @@ void PeriodicTorsionForce::getTorsionParameters(int index, int& particle1, int&
 }
 
 void PeriodicTorsionForce::setTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4, int periodicity, double phase, double k) {
+    ASSERT_VALID_INDEX(index, periodicTorsions);
     periodicTorsions[index].particle1 = particle1;
     periodicTorsions[index].particle2 = particle2;
     periodicTorsions[index].particle3 = particle3;

openmmapi/src/RBTorsionForce.cpp

@@ -32,6 +32,7 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/RBTorsionForce.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/RBTorsionForceImpl.h"
 
 using namespace OpenMM;
@@ -45,6 +46,7 @@ int RBTorsionForce::addTorsion(int particle1, int particle2, int particle3, int
 }
 
 void RBTorsionForce::getTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, double& c0, double& c1, double& c2, double& c3, double& c4, double& c5) const {
+    ASSERT_VALID_INDEX(index, rbTorsions);
     particle1 = rbTorsions[index].particle1;
     particle2 = rbTorsions[index].particle2;
     particle3 = rbTorsions[index].particle3;
@@ -58,6 +60,7 @@ void RBTorsionForce::getTorsionParameters(int index, int& particle1, int& partic
 }
 
 void RBTorsionForce::setTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4, double c0, double c1, double c2, double c3, double c4, double c5) {
+    ASSERT_VALID_INDEX(index, rbTorsions);
     rbTorsions[index].particle1 = particle1;
     rbTorsions[index].particle2 = particle2;
     rbTorsions[index].particle3 = particle3;

openmmapi/src/System.cpp

@@ -32,6 +32,8 @@
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/System.h"
+#include "openmm/VirtualSite.h"
+#include "openmm/internal/AssertionUtilities.h"
 
 using namespace OpenMM;
 
@@ -48,6 +50,17 @@ System::~System() {
         delete virtualSites[i];
 }
 
+double System::getParticleMass(int index) const {
+    ASSERT_VALID_INDEX(index, masses);
+    return masses[index];
+}
+
+void System::setParticleMass(int index, double mass) {
+    ASSERT_VALID_INDEX(index, masses);
+    masses[index] = mass;
+}
+
+
 void System::setVirtualSite(int index, VirtualSite* virtualSite) {
     if (index >= virtualSites.size())
         virtualSites.resize(getNumParticles(), NULL);
@@ -66,17 +79,29 @@ int System::addConstraint(int particle1, int particle2, double distance) {
 }
 
 void System::getConstraintParameters(int index, int& particle1, int& particle2, double& distance) const {
+    ASSERT_VALID_INDEX(index, constraints);
     particle1 = constraints[index].particle1;
     particle2 = constraints[index].particle2;
     distance = constraints[index].distance;
 }
 
 void System::setConstraintParameters(int index, int particle1, int particle2, double distance) {
+    ASSERT_VALID_INDEX(index, constraints);
     constraints[index].particle1 = particle1;
     constraints[index].particle2 = particle2;
     constraints[index].distance = distance;
 }
 
+const Force& System::getForce(int index) const {
+    ASSERT_VALID_INDEX(index, forces);
+    return *forces[index];
+}
+
+Force& System::getForce(int index) {
+    ASSERT_VALID_INDEX(index, forces);
+    return *forces[index];
+}
+
 void System::getDefaultPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) const {
     a = periodicBoxVectors[0];
     b = periodicBoxVectors[1];