Added error checking for illegal variable names

openmmapi/include/openmm/internal/CustomCentroidBondForceImpl.h

@@ -41,6 +41,7 @@
 #include "lepton/ParsedExpression.h"
 #include <utility>
 #include <map>
+#include <set>
 #include <string>
 
 namespace OpenMM {
@@ -93,7 +94,7 @@ private:
     class FunctionPlaceholder;
     static Lepton::ExpressionTreeNode replaceFunctions(const Lepton::ExpressionTreeNode& node, std::map<std::string, int> atoms,
             std::map<std::string, std::vector<int> >& distances, std::map<std::string, std::vector<int> >& angles,
-            std::map<std::string, std::vector<int> >& dihedrals);
+            std::map<std::string, std::vector<int> >& dihedrals, std::set<std::string>& variables);
     void addBondsBetweenGroups(int group1, int group2, std::vector<std::pair<int, int> >& bonds) const;
     const CustomCentroidBondForce& owner;
     Kernel kernel;

openmmapi/include/openmm/internal/CustomCompoundBondForceImpl.h

@@ -40,6 +40,7 @@
 #include "lepton/ParsedExpression.h"
 #include <utility>
 #include <map>
+#include <set>
 #include <string>
 
 namespace OpenMM {
@@ -83,7 +84,7 @@ private:
     class FunctionPlaceholder;
     static Lepton::ExpressionTreeNode replaceFunctions(const Lepton::ExpressionTreeNode& node, std::map<std::string, int> atoms,
             std::map<std::string, std::vector<int> >& distances, std::map<std::string, std::vector<int> >& angles,
-            std::map<std::string, std::vector<int> >& dihedrals);
+            std::map<std::string, std::vector<int> >& dihedrals, std::set<std::string>& variables);
     const CustomCompoundBondForce& owner;
     Kernel kernel;
 };

openmmapi/include/openmm/internal/CustomHbondForceImpl.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-2010 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -40,6 +40,7 @@
 #include "lepton/ParsedExpression.h"
 #include <utility>
 #include <map>
+#include <set>
 #include <string>
 
 namespace OpenMM {
@@ -84,7 +85,7 @@ private:
     class FunctionPlaceholder;
     static Lepton::ExpressionTreeNode replaceFunctions(const Lepton::ExpressionTreeNode& node, std::map<std::string, int> atoms,
             std::map<std::string, std::vector<int> >& distances, std::map<std::string, std::vector<int> >& angles,
-            std::map<std::string, std::vector<int> >& dihedrals);
+            std::map<std::string, std::vector<int> >& dihedrals, std::set<std::string>& variables);
     const CustomHbondForce& owner;
     Kernel kernel;
 };

openmmapi/include/openmm/internal/CustomManyParticleForceImpl.h

@@ -40,6 +40,7 @@
 #include "lepton/ParsedExpression.h"
 #include <utility>
 #include <map>
+#include <set>
 #include <string>
 
 namespace OpenMM {
@@ -98,7 +99,7 @@ private:
     class FunctionPlaceholder;
     static Lepton::ExpressionTreeNode replaceFunctions(const Lepton::ExpressionTreeNode& node, std::map<std::string, int> atoms,
             std::map<std::string, std::vector<int> >& distances, std::map<std::string, std::vector<int> >& angles,
-            std::map<std::string, std::vector<int> >& dihedrals);
+            std::map<std::string, std::vector<int> >& dihedrals, std::set<std::string>& variables);
     static void generatePermutations(std::vector<int>& values, int numFixed, std::vector<std::vector<int> >& result);
     const CustomManyParticleForce& owner;
     Kernel kernel;

openmmapi/src/CustomCentroidBondForceImpl.cpp

@@ -150,24 +150,37 @@ ParsedExpression CustomCentroidBondForceImpl::prepareExpression(const CustomCent
     functions["dihedral"] = &dihedral;
     ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
     map<string, int> groups;
+    set<string> variables;
     for (int i = 0; i < force.getNumGroupsPerBond(); i++) {
-        stringstream name;
+        stringstream name, x, y, z;
         name << 'g' << (i+1);
+        x << 'x' << (i+1);
+        y << 'y' << (i+1);
+        z << 'z' << (i+1);
         groups[name.str()] = i;
+        variables.insert(x.str());
+        variables.insert(y.str());
+        variables.insert(z.str());
     }
-    return ParsedExpression(replaceFunctions(expression.getRootNode(), groups, distances, angles, dihedrals)).optimize();
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        variables.insert(force.getGlobalParameterName(i));
+    for (int i = 0; i < force.getNumPerBondParameters(); i++)
+        variables.insert(force.getPerBondParameterName(i));
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), groups, distances, angles, dihedrals, variables)).optimize();
 }
 
 ExpressionTreeNode CustomCentroidBondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> groups,
-        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals, set<string>& variables) {
     const Operation& op = node.getOperation();
+    if (op.getId() == Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("CustomCentroidBondForce: Unknown variable '"+op.getName()+"'");
     if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
     {
         // This is not an angle or dihedral, so process its children.
 
         vector<ExpressionTreeNode> children;
         for (int i = 0; i < (int) node.getChildren().size(); i++)
-            children.push_back(replaceFunctions(node.getChildren()[i], groups, distances, angles, dihedrals));
+            children.push_back(replaceFunctions(node.getChildren()[i], groups, distances, angles, dihedrals, variables));
         return ExpressionTreeNode(op.clone(), children);
     }
     const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);

openmmapi/src/CustomCompoundBondForceImpl.cpp

@@ -136,24 +136,37 @@ ParsedExpression CustomCompoundBondForceImpl::prepareExpression(const CustomComp
     functions["dihedral"] = &dihedral;
     ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
     map<string, int> atoms;
+    set<string> variables;
     for (int i = 0; i < force.getNumParticlesPerBond(); i++) {
-        stringstream name;
+        stringstream name, x, y, z;
         name << 'p' << (i+1);
+        x << 'x' << (i+1);
+        y << 'y' << (i+1);
+        z << 'z' << (i+1);
         atoms[name.str()] = i;
+        variables.insert(x.str());
+        variables.insert(y.str());
+        variables.insert(z.str());
     }
-    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        variables.insert(force.getGlobalParameterName(i));
+    for (int i = 0; i < force.getNumPerBondParameters(); i++)
+        variables.insert(force.getPerBondParameterName(i));
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals, variables)).optimize();
 }
 
 ExpressionTreeNode CustomCompoundBondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
-        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals, set<string>& variables) {
     const Operation& op = node.getOperation();
+    if (op.getId() == Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("CustomCompoundBondForce: Unknown variable '"+op.getName()+"'");
     if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
     {
         // This is not an angle or dihedral, so process its children.
 
         vector<ExpressionTreeNode> children;
         for (int i = 0; i < (int) node.getChildren().size(); i++)
-            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals));
+            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals, variables));
         return ExpressionTreeNode(op.clone(), children);
     }
     const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);

openmmapi/src/CustomHbondForceImpl.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-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -215,19 +215,28 @@ ParsedExpression CustomHbondForceImpl::prepareExpression(const CustomHbondForce&
     atoms["d1"] = 3;
     atoms["d2"] = 4;
     atoms["d3"] = 5;
-    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
+    set<string> variables;
+    for (int i = 0; i < force.getNumPerDonorParameters(); i++)
+        variables.insert(force.getPerDonorParameterName(i));
+    for (int i = 0; i < force.getNumPerAcceptorParameters(); i++)
+        variables.insert(force.getPerAcceptorParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        variables.insert(force.getGlobalParameterName(i));
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals, variables)).optimize();
 }
 
 ExpressionTreeNode CustomHbondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
-        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals, set<string>& variables) {
     const Operation& op = node.getOperation();
+    if (op.getId() == Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("CustomHBondForce: Unknown variable '"+op.getName()+"'");
     if (op.getId() != Operation::CUSTOM || op.getNumArguments() < 2)
     {
         // This is not an angle or dihedral, so process its children.
 
         vector<ExpressionTreeNode> children;
         for (int i = 0; i < (int) node.getChildren().size(); i++)
-            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals));
+            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals, variables));
         return ExpressionTreeNode(op.clone(), children);
     }
     const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);

openmmapi/src/CustomManyParticleForceImpl.cpp

@@ -165,24 +165,40 @@ ParsedExpression CustomManyParticleForceImpl::prepareExpression(const CustomMany
     functions["dihedral"] = &dihedral;
     ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
     map<string, int> atoms;
+    set<string> variables;
     for (int i = 0; i < force.getNumParticlesPerSet(); i++) {
-        stringstream name;
+        stringstream name, x, y, z;
         name << 'p' << (i+1);
+        x << 'x' << (i+1);
+        y << 'y' << (i+1);
+        z << 'z' << (i+1);
         atoms[name.str()] = i;
+        variables.insert(x.str());
+        variables.insert(y.str());
+        variables.insert(z.str());
+        for (int j = 0; j < force.getNumPerParticleParameters(); j++) {
+            stringstream param;
+            param << force.getPerParticleParameterName(j) << (i+1);
+            variables.insert(param.str());
+        }
     }
-    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        variables.insert(force.getGlobalParameterName(i));
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals, variables)).optimize();
 }
 
 ExpressionTreeNode CustomManyParticleForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
-        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals, set<string>& variables) {
     const Operation& op = node.getOperation();
+    if (op.getId() == Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("CustomManyParticleForce: Unknown variable '"+op.getName()+"'");
     if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
     {
         // This is not an angle or dihedral, so process its children.
 
         vector<ExpressionTreeNode> children;
         for (int i = 0; i < (int) node.getChildren().size(); i++)
-            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals));
+            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals, variables));
         return ExpressionTreeNode(op.clone(), children);
     }
     const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);

platforms/cpu/src/CpuKernels.cpp

@@ -47,6 +47,7 @@
 #include "RealVec.h"
 #include "lepton/CompiledExpression.h"
 #include "lepton/CustomFunction.h"
+#include "lepton/Operation.h"
 #include "lepton/Parser.h"
 #include "lepton/ParsedExpression.h"
 
@@ -83,6 +84,17 @@ static ReferenceConstraints& extractConstraints(ContextImpl& context) {
     return *(ReferenceConstraints*) data->constraints;
 }
 
+/**
+ * Make sure an expression doesn't use any undefined variables.
+ */
+static void validateVariables(const Lepton::ExpressionTreeNode& node, const set<string>& variables) {
+    const Lepton::Operation& op = node.getOperation();
+    if (op.getId() == Lepton::Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("Unknown variable in expression: "+op.getName());
+    for (int i = 0; i < (int) node.getChildren().size(); i++)
+        validateVariables(node.getChildren()[i], variables);
+}
+
 /**
  * Compute the kinetic energy of the system, possibly shifting the velocities in time to account
  * for a leapfrog integrator.
@@ -737,6 +749,14 @@ void CpuCalcCustomNonbondedForceKernel::initialize(const System& system, const C
         globalParameterNames.push_back(force.getGlobalParameterName(i));
         globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
     }
+    set<string> variables;
+    variables.insert("r");
+    for (int i = 0; i < numParameters; i++) {
+        variables.insert(parameterNames[i]+"1");
+        variables.insert(parameterNames[i]+"2");
+    }
+    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    validateVariables(expression.getRootNode(), variables);
 
     // Delete the custom functions.
 
@@ -950,6 +970,18 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
     vector<vector<Lepton::CompiledExpression> > valueGradientExpressions(force.getNumComputedValues());
     vector<Lepton::CompiledExpression> valueExpressions;
     vector<Lepton::CompiledExpression> energyExpressions;
+    set<string> particleVariables, pairVariables;
+    pairVariables.insert("r");
+    particleVariables.insert("x");
+    particleVariables.insert("y");
+    particleVariables.insert("z");
+    for (int i = 0; i < numPerParticleParameters; i++) {
+        particleVariables.insert(particleParameterNames[i]);
+        pairVariables.insert(particleParameterNames[i]+"1");
+        pairVariables.insert(particleParameterNames[i]+"2");
+    }
+    particleVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    pairVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
     for (int i = 0; i < force.getNumComputedValues(); i++) {
         string name, expression;
         CustomGBForce::ComputationType type;
@@ -958,15 +990,21 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
         valueExpressions.push_back(ex.createCompiledExpression());
         valueTypes.push_back(type);
         valueNames.push_back(name);
-        if (i == 0)
+        if (i == 0) {
             valueDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
+            validateVariables(ex.getRootNode(), pairVariables);
+        }
         else {
             valueGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
             valueGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
             valueGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
             for (int j = 0; j < i; j++)
                 valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).createCompiledExpression());
+            validateVariables(ex.getRootNode(), particleVariables);
         }
+        particleVariables.insert(name);
+        pairVariables.insert(name+"1");
+        pairVariables.insert(name+"2");
     }
 
     // Parse the expressions for energy terms.
@@ -988,10 +1026,12 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
                 energyGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
                 energyGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
                 energyGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
+                validateVariables(ex.getRootNode(), particleVariables);
             }
             else {
                 energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").createCompiledExpression());
                 energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").createCompiledExpression());
+                validateVariables(ex.getRootNode(), pairVariables);
             }
         }
     }

platforms/reference/src/ReferenceKernels.cpp

@@ -77,6 +77,7 @@
 #include "openmm/OpenMMException.h"
 #include "SimTKOpenMMUtilities.h"
 #include "lepton/CustomFunction.h"
+#include "lepton/Operation.h"
 #include "lepton/Parser.h"
 #include "lepton/ParsedExpression.h"
 #include <cmath>
@@ -146,6 +147,17 @@ static ReferenceConstraints& extractConstraints(ContextImpl& context) {
     return *(ReferenceConstraints*) data->constraints;
 }
 
+/**
+ * Make sure an expression doesn't use any undefined variables.
+ */
+static void validateVariables(const Lepton::ExpressionTreeNode& node, const set<string>& variables) {
+    const Lepton::Operation& op = node.getOperation();
+    if (op.getId() == Lepton::Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("Unknown variable in expression: "+op.getName());
+    for (int i = 0; i < (int) node.getChildren().size(); i++)
+        validateVariables(node.getChildren()[i], variables);
+}
+
 /**
  * Compute the kinetic energy of the system, possibly shifting the velocities in time to account
  * for a leapfrog integrator.
@@ -422,6 +434,11 @@ void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const
         parameterNames.push_back(force.getPerBondParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
+    set<string> variables;
+    variables.insert("r");
+    variables.insert(parameterNames.begin(), parameterNames.end());
+    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    validateVariables(expression.getRootNode(), variables);
 }
 
 double ReferenceCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -536,6 +553,11 @@ void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const
         parameterNames.push_back(force.getPerAngleParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
+    set<string> variables;
+    variables.insert("theta");
+    variables.insert(parameterNames.begin(), parameterNames.end());
+    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    validateVariables(expression.getRootNode(), variables);
 }
 
 double ReferenceCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -780,6 +802,11 @@ void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, con
         parameterNames.push_back(force.getPerTorsionParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
+    set<string> variables;
+    variables.insert("theta");
+    variables.insert(parameterNames.begin(), parameterNames.end());
+    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    validateVariables(expression.getRootNode(), variables);
 }
 
 double ReferenceCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1038,6 +1065,14 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
         globalParameterNames.push_back(force.getGlobalParameterName(i));
         globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
     }
+    set<string> variables;
+    variables.insert("r");
+    for (int i = 0; i < numParameters; i++) {
+        variables.insert(parameterNames[i]+"1");
+        variables.insert(parameterNames[i]+"2");
+    }
+    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    validateVariables(expression.getRootNode(), variables);
 
     // Delete the custom functions.
 
@@ -1314,6 +1349,18 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
 
     valueDerivExpressions.resize(force.getNumComputedValues());
     valueGradientExpressions.resize(force.getNumComputedValues());
+    set<string> particleVariables, pairVariables;
+    pairVariables.insert("r");
+    particleVariables.insert("x");
+    particleVariables.insert("y");
+    particleVariables.insert("z");
+    for (int i = 0; i < numPerParticleParameters; i++) {
+        particleVariables.insert(particleParameterNames[i]);
+        pairVariables.insert(particleParameterNames[i]+"1");
+        pairVariables.insert(particleParameterNames[i]+"2");
+    }
+    particleVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    pairVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
     for (int i = 0; i < force.getNumComputedValues(); i++) {
         string name, expression;
         CustomGBForce::ComputationType type;
@@ -1322,15 +1369,21 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
         valueExpressions.push_back(ex.createProgram());
         valueTypes.push_back(type);
         valueNames.push_back(name);
-        if (i == 0)
+        if (i == 0) {
             valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
+            validateVariables(ex.getRootNode(), pairVariables);
+        }
         else {
             valueGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
             valueGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
             valueGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
             for (int j = 0; j < i; j++)
                 valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
+            validateVariables(ex.getRootNode(), particleVariables);
         }
+        particleVariables.insert(name);
+        pairVariables.insert(name+"1");
+        pairVariables.insert(name+"2");
     }
 
     // Parse the expressions for energy terms.
@@ -1352,10 +1405,12 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
                 energyGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
                 energyGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
                 energyGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
+                validateVariables(ex.getRootNode(), particleVariables);
             }
             else {
                 energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
                 energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").optimize().createProgram());
+                validateVariables(ex.getRootNode(), pairVariables);
             }
         }
     }
@@ -1475,6 +1530,13 @@ void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, co
         parameterNames.push_back(force.getPerParticleParameterName(i));
     for (int i = 0; i < force.getNumGlobalParameters(); i++)
         globalParameterNames.push_back(force.getGlobalParameterName(i));
+    set<string> variables;
+    variables.insert("x");
+    variables.insert("y");
+    variables.insert("z");
+    variables.insert(parameterNames.begin(), parameterNames.end());
+    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    validateVariables(expression.getRootNode(), variables);
 }
 
 double ReferenceCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {

tests/TestCustomAngleForce.h

@@ -126,12 +126,33 @@ void testAngles() {
     }
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    CustomAngleForce* force = new CustomAngleForce("theta+none");
+    vector<double> params;
+    force->addAngle(0, 1, 2, params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testAngles();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomBondForce.h

@@ -130,6 +130,25 @@ void testManyParameters() {
     ASSERT_EQUAL_TOL(f*2.5, state.getPotentialEnergy(), TOL);
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    CustomBondForce* force = new CustomBondForce("r+none");
+    vector<double> params;
+    force->addBond(0, 1, params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -137,6 +156,7 @@ int main(int argc, char* argv[]) {
         initializeTests(argc, argv);
         testBonds();
         testManyParameters();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomCentroidBondForce.h

@@ -252,6 +252,32 @@ void testCustomWeights() {
     ASSERT_EQUAL_VEC(Vec3(0, -2*9*(1.0/3.0), 0), state.getForces()[3], TOL);
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "1+none");
+    vector<int> particles;
+    particles.push_back(0);
+    force->addGroup(particles);
+    force->addGroup(particles);
+    vector<int> groups;
+    groups.push_back(0);
+    groups.push_back(1);
+    vector<double> params;
+    force->addBond(groups, params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -260,6 +286,7 @@ int main(int argc, char* argv[]) {
         testHarmonicBond();
         testComplexFunction();
         testCustomWeights();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomCompoundBondForce.h

@@ -309,6 +309,28 @@ void testMultipleBonds() {
     ASSERT_EQUAL_VEC(Vec3(0.776, 0.0, -0.084), forces[3], 1e-3);
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    CustomCompoundBondForce* force = new CustomCompoundBondForce(2, "1+none");
+    vector<int> particles;
+    particles.push_back(0);
+    particles.push_back(1);
+    vector<double> params;
+    force->addBond(particles, params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -319,6 +341,7 @@ int main(int argc, char* argv[]) {
         testContinuous2DFunction();
         testContinuous3DFunction();
         testMultipleBonds();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomExternalForce.h

@@ -165,7 +165,24 @@ void testPeriodic() {
         ASSERT_EQUAL_TOL(delta.dot(delta), state.getPotentialEnergy(), TOL);
         integrator.step(1);
     }
+}
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    CustomExternalForce* force = new CustomExternalForce("x+none");
+    vector<double> params;
+    force->addParticle(0, params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
 }
 
 void runPlatformTests();
@@ -176,6 +193,7 @@ int main(int argc, char* argv[]) {
         testForce();
         testManyParameters();
         testPeriodic();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomGBForce.h

@@ -470,6 +470,25 @@ void testExclusions() {
     }
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    CustomGBForce* force = new CustomGBForce();
+    force->addComputedValue("a", "r+none", CustomGBForce::ParticlePair);
+    vector<double> params;
+    force->addParticle(params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -483,6 +502,7 @@ int main(int argc, char* argv[]) {
         testMultipleChainRules();
         testPositionDependence();
         testExclusions();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomHbondForce.h

@@ -225,6 +225,26 @@ void testCustomFunctions() {
     ASSERT_EQUAL_TOL(0.1*2+0.1*2, state.getPotentialEnergy(), TOL);
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    CustomHbondForce* force = new CustomHbondForce("1+none");
+    vector<double> params;
+    force->addDonor(0, -1, -1, params);
+    force->addAcceptor(1, -1, -1, params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -234,6 +254,7 @@ int main(int argc, char* argv[]) {
         testExclusions();
         testCutoff();
         testCustomFunctions();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomManyParticleForce.h

@@ -707,6 +707,24 @@ void testCentralParticleModeLargeSystem() {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    CustomManyParticleForce* force = new CustomManyParticleForce(2, "x1+y2+none");
+    vector<double> params;
+    force->addParticle(params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -725,6 +743,7 @@ int main(int argc, char* argv[]) {
         testCentralParticleModeNoCutoff();
         testCentralParticleModeCutoff();
         testCentralParticleModeLargeSystem();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomNonbondedForce.h

@@ -993,6 +993,24 @@ void testMultipleCutoffs() {
     }
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    CustomNonbondedForce* force = new CustomNonbondedForce("r+none");
+    vector<double> params;
+    force->addParticle(params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -1017,6 +1035,7 @@ int main(int argc, char* argv[]) {
         testLargeInteractionGroup();
         testInteractionGroupLongRangeCorrection();
         testMultipleCutoffs();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomTorsionForce.h

@@ -166,6 +166,27 @@ void testRange() {
     ASSERT(maxAngle <= M_PI);
 }
 
+void testIllegalVariable() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    CustomTorsionForce* force = new CustomTorsionForce("theta+none");
+    vector<double> params;
+    force->addTorsion(0, 1, 2, 3, params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    bool threwException = false;
+    try {
+        Context(system, integrator, platform);
+    }
+    catch (const exception& e) {
+        threwException = true;
+    }
+    ASSERT(threwException);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -173,6 +194,7 @@ int main(int argc, char* argv[]) {
         initializeTests(argc, argv);
         testTorsions();
         testRange();
+        testIllegalVariable();
         runPlatformTests();
     }
     catch(const exception& e) {