CustomNonbondedForce supports computed values (#3412)

* Reference implementation of computed values for CustomNonbondedForce

* CPU implementation of computed values for CustomNonbondedForce

* Common implementation of computed values for CustomNonbondedForce

* Serialization of computed values

* ForceField supports computed values

docs-source/usersguide/application/05_creating_ffs.rst

@@ -761,6 +761,11 @@ also may depend on an arbitrary list of global or per-particle parameters.  Use
 a :code:`<GlobalParameter>` tag to define a global parameter, and a
 :code:`<PerParticleParameter>` tag to define a per-particle parameter.
 
+The expression may also depend on computed values, each defined with a
+:code:`<ComputedValue>` tag.  The tag should have two attributes, :code:`name`
+with the name of the computed value, and :code:`expression` with the expression
+used to compute it.
+
 Exclusions are created automatically based on the :code:`bondCutoff` attribute.
 After setting the nonbonded parameters for all atoms, the force field calls
 :code:`createExclusionsFromBonds()` on the CustomNonbondedForce, passing in this

openmmapi/include/openmm/CustomNonbondedForce.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-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -88,7 +88,27 @@ namespace OpenMM {
  * above example, the energy only depends on the products sigma1*sigma2 and epsilon1*epsilon2, both of which are unchanged
  * if the labels 1 and 2 are reversed.  In contrast, if it depended on the difference sigma1-sigma2, the results would
  * be undefined, because reversing the labels 1 and 2 would change the energy.
+ * 
+ * The energy also may depend on "computed values".  These are similar to per-particle parameters, but instead of being
+ * specified in advance, their values are computed based on global and per-particle parameters.  For example, the following
+ * code uses a global parameter (lambda) to interpolate between two different sigma values for each particle (sigmaA and sigmaB).
+ * 
+ * \verbatim embed:rst:leading-asterisk
+ * .. code-block:: cpp
  *
+ *    CustomNonbondedForce* force = new CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)");
+ *    force->addComputedValue("sigma", "(1-lambda)*sigmaA + lambda*sigmaB");
+ *    force->addGlobalParameter("lambda", 0);
+ *    force->addPerParticleParameter("sigmaA");
+ *    force->addPerParticleParameter("sigmaB");
+ *    force->addPerParticleParameter("epsilon");
+ *
+ * \endverbatim
+ *
+ * You could, of course, embed the computation of sigma directly into the energy expression, but then it would need to be
+ * repeated for every interaction.  By separating it out as a computed value, it only needs to be computed once for each
+ * particle instead of once for each interaction, thus saving computation time.
+ * 
  * CustomNonbondedForce can operate in two modes.  By default, it computes the interaction of every particle in the System
  * with every other particle.  Alternatively, you can restrict it to only a subset of particle pairs.  To do this, specify
  * one or more "interaction groups".  An interaction group consists of two sets of particles that should interact with
@@ -210,6 +230,12 @@ public:
     int getNumFunctions() const {
         return functions.size();
     }
+    /**
+     * Get the number of per-particle computed values the interaction depends on.
+     */
+    int getNumComputedValues() const {
+        return computedValues.size();
+    }
     /**
      * Get the number of interaction groups that have been defined.
      */
@@ -465,6 +491,36 @@ public:
      * If the specified function is not a Continuous1DFunction, this throws an exception.
      */
     void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max);
+    /**
+     * Add a computed value to calculate for each particle.
+     *
+     * @param name        the name of the value
+     * @param expression  an algebraic expression to evaluate when calculating the computed value.  It may
+     *                    depend on the values of per-particle and global parameters, but not one other
+     *                    computed values.
+     * @return the index of the computed value that was added
+     */
+    int addComputedValue(const std::string& name, const std::string& expression);
+    /**
+     * Get the properties of a computed value.
+     *
+     * @param index            the index of the computed value for which to get parameters
+     * @param[out] name        the name of the value
+     * @param[out] expression  an algebraic expression to evaluate when calculating the computed value.  It may
+     *                         depend on the values of per-particle and global parameters, but not one other
+     *                         computed values.
+     */
+    void getComputedValueParameters(int index, std::string& name, std::string& expression) const;
+    /**
+     * Set the properties of a computed value.
+     *
+     * @param index       the index of the computed value for which to set parameters
+     * @param name        the name of the value
+     * @param expression  an algebraic expression to evaluate when calculating the computed value.  It may
+     *                    depend on the values of per-particle and global parameters, but not one other
+     *                    computed values.
+     */
+    void setComputedValueParameters(int index, const std::string& name, const std::string& expression);
     /**
      * Add an interaction group.  An interaction will be computed between every particle in set1 and every particle in set2.
      *
@@ -520,6 +576,7 @@ private:
     class GlobalParameterInfo;
     class ExclusionInfo;
     class FunctionInfo;
+    class ComputedValueInfo;
     class InteractionGroupInfo;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance, switchingDistance;
@@ -530,6 +587,7 @@ private:
     std::vector<ParticleInfo> particles;
     std::vector<ExclusionInfo> exclusions;
     std::vector<FunctionInfo> functions;
+    std::vector<ComputedValueInfo> computedValues;
     std::vector<InteractionGroupInfo> interactionGroups;
     std::vector<int> energyParameterDerivatives;
 };
@@ -603,6 +661,19 @@ public:
     }
 };
 
+/**
+ * This is an internal class used to record information about a computed value.
+ * @private
+ */
+class CustomNonbondedForce::ComputedValueInfo {
+public:
+    std::string name, expression;
+    ComputedValueInfo() {
+    }
+    ComputedValueInfo(const std::string& name, const std::string& expression) : name(name), expression(expression) {
+    }
+};
+
 /**
  * This is an internal class used to record information about an interaction group.
  * @private

openmmapi/include/openmm/internal/CustomNonbondedForceImpl.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-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -78,7 +78,8 @@ public:
     static void calcLongRangeCorrection(const CustomNonbondedForce& force, LongRangeCorrectionData& data, const Context& context, double& coefficient, std::vector<double>& derivatives, ThreadPool& threads);
 private:
     static double integrateInteraction(Lepton::CompiledExpression& expression, const std::vector<double>& params1, const std::vector<double>& params2,
-            const CustomNonbondedForce& force, const Context& context, const std::vector<std::string>& paramNames);
+            const std::vector<double>& computedValues1, const std::vector<double>& computedValues2, const CustomNonbondedForce& force, const Context& context,
+            const std::vector<std::string>& paramNames, const std::vector<std::string>& computedValueNames);
     const CustomNonbondedForce& owner;
     Kernel kernel;
 };
@@ -87,10 +88,10 @@ class CustomNonbondedForceImpl::LongRangeCorrectionData {
 public:
     CustomNonbondedForce::NonbondedMethod method;
     std::vector<std::vector<double> > classes;
-    std::vector<std::string> paramNames;
+    std::vector<std::string> paramNames, computedValueNames;
     std::map<std::pair<int, int>, long long int> interactionCount;
     Lepton::CompiledExpression energyExpression;
-    std::vector<Lepton::CompiledExpression> derivExpressions;
+    std::vector<Lepton::CompiledExpression> derivExpressions, computedValueExpressions;
 };
 
 } // namespace OpenMM

openmmapi/src/CustomNonbondedForce.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-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -63,6 +63,7 @@ CustomNonbondedForce::CustomNonbondedForce(const CustomNonbondedForce& rhs) {
     useLongRangeCorrection = rhs.useLongRangeCorrection;
     parameters = rhs.parameters;
     globalParameters = rhs.globalParameters;
+    computedValues = rhs.computedValues;
     energyParameterDerivatives = rhs.energyParameterDerivatives;
     particles = rhs.particles;
     exclusions = rhs.exclusions;
@@ -282,6 +283,23 @@ void CustomNonbondedForce::setFunctionParameters(int index, const std::string& n
     function->setFunctionParameters(values, min, max);
 }
 
+int CustomNonbondedForce::addComputedValue(const string& name, const string& expression) {
+    computedValues.push_back(ComputedValueInfo(name, expression));
+    return computedValues.size()-1;
+}
+
+void CustomNonbondedForce::getComputedValueParameters(int index, string& name, string& expression) const {
+    ASSERT_VALID_INDEX(index, computedValues);
+    name = computedValues[index].name;
+    expression = computedValues[index].expression;
+}
+
+void CustomNonbondedForce::setComputedValueParameters(int index, const string& name, const string& expression) {
+    ASSERT_VALID_INDEX(index, computedValues);
+    computedValues[index].name = name;
+    computedValues[index].expression = expression;
+}
+
 int CustomNonbondedForce::addInteractionGroup(const std::set<int>& set1, const std::set<int>& set2) {
     for (set<int>::iterator it = set1.begin(); it != set1.end(); ++it)
         ASSERT(*it >= 0);

openmmapi/src/CustomNonbondedForceImpl.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-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -131,6 +131,8 @@ void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
                 throw OpenMMException(msg.str());
             }
     }
+    if (owner.getNumEnergyParameterDerivatives() > 0 && owner.getNumComputedValues() > 0)
+        throw OpenMMException("CustomNonbondedForce: Cannot compute parameter derivatives for a force that uses computed values.");
 
     kernel.getAs<CalcCustomNonbondedForceKernel>().initialize(context.getSystem(), owner);
 }
@@ -231,11 +233,16 @@ CustomNonbondedForceImpl::LongRangeCorrectionData CustomNonbondedForceImpl::prep
     for (int k = 0; k < force.getNumEnergyParameterDerivatives(); k++)
         data.derivExpressions.push_back(Lepton::Parser::parse(force.getEnergyFunction(), functions).differentiate(force.getEnergyParameterDerivativeName(k)).createCompiledExpression());
     for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
-        stringstream name1, name2;
-        name1 << force.getPerParticleParameterName(i) << 1;
-        name2 << force.getPerParticleParameterName(i) << 2;
-        data.paramNames.push_back(name1.str());
-        data.paramNames.push_back(name2.str());
+        string name = force.getPerParticleParameterName(i);
+        data.paramNames.push_back(name+"1");
+        data.paramNames.push_back(name+"2");
+    }
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        string name, exp;
+        force.getComputedValueParameters(i, name, exp);
+        data.computedValueNames.push_back(name+"1");
+        data.computedValueNames.push_back(name+"2");
+        data.computedValueExpressions.push_back(Lepton::Parser::parse(exp, functions).createCompiledExpression());
     }
     return data;
 }
@@ -245,10 +252,31 @@ void CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedForc
         coefficient = 0.0;
         return;
     }
+    
+    // Calculate the computed values for all atom classes.
+    
+    int numClasses = data.classes.size();
+    vector<vector<double> > computedValues(numClasses, vector<double>(force.getNumComputedValues()));
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        Lepton::CompiledExpression& expression = data.computedValueExpressions[i];
+        const set<string>& variables = expression.getVariables();
+        for (int j = 0; j < force.getNumGlobalParameters(); j++) {
+            const string& name = force.getGlobalParameterName(j);
+            if (variables.find(name) != variables.end())
+                expression.getVariableReference(name) = context.getParameter(name);
+        }
+        for (int j = 0; j < numClasses; j++) {
+            for (int k = 0; k < force.getNumPerParticleParameters(); k++) {
+                const string& name = force.getPerParticleParameterName(k);
+                if (variables.find(name) != variables.end())
+                    expression.getVariableReference(name) = data.classes[j][k];
+            }
+            computedValues[j][i] = expression.evaluate();
+        }
+    }
 
     // Compute the coefficient.  Use multiple threads to compute the integrals in parallel.
 
-    int numClasses = data.classes.size();
     double nPart = (double) context.getSystem().getNumParticles();
     double numInteractions = (nPart*(nPart+1))/2;
     vector<double> threadSum(threads.getNumThreads(), 0.0);
@@ -260,7 +288,8 @@ void CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedForc
             if (i >= numClasses)
                 break;
             for (int j = i; j < numClasses; j++)
-                threadSum[threadIndex] += data.interactionCount.at(make_pair(i, j))*integrateInteraction(expression, data.classes[i], data.classes[j], force, context, data.paramNames);
+                threadSum[threadIndex] += data.interactionCount.at(make_pair(i, j))*integrateInteraction(expression, data.classes[i], data.classes[j],
+                        computedValues[i], computedValues[j], force, context, data.paramNames, data.computedValueNames);
         }
     });
     threads.waitForThreads();
@@ -284,7 +313,8 @@ void CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedForc
                 if (i >= numClasses)
                     break;
                 for (int j = i; j < numClasses; j++)
-                    threadSum[threadIndex] += data.interactionCount.at(make_pair(i, j))*integrateInteraction(expression, data.classes[i], data.classes[j], force, context, data.paramNames);
+                    threadSum[threadIndex] += data.interactionCount.at(make_pair(i, j))*integrateInteraction(expression, data.classes[i], data.classes[j],
+                            computedValues[i], computedValues[j], force, context, data.paramNames, data.computedValueNames);
             }
         });
         threads.waitForThreads();
@@ -297,7 +327,8 @@ void CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedForc
 }
 
 double CustomNonbondedForceImpl::integrateInteraction(Lepton::CompiledExpression& expression, const vector<double>& params1, const vector<double>& params2,
-        const CustomNonbondedForce& force, const Context& context, const vector<string>& paramNames) {
+        const vector<double>& computedValues1, const vector<double>& computedValues2, const CustomNonbondedForce& force, const Context& context,
+        const vector<string>& paramNames, const vector<string>& computedValueNames) {
     const set<string>& variables = expression.getVariables();
     for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
         if (variables.find(paramNames[2*i]) != variables.end())
@@ -305,6 +336,12 @@ double CustomNonbondedForceImpl::integrateInteraction(Lepton::CompiledExpression
         if (variables.find(paramNames[2*i+1]) != variables.end())
             expression.getVariableReference(paramNames[2*i+1]) = params2[i];
     }
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        if (variables.find(computedValueNames[2*i]) != variables.end())
+            expression.getVariableReference(computedValueNames[2*i]) = computedValues1[i];
+        if (variables.find(computedValueNames[2*i+1]) != variables.end())
+            expression.getVariableReference(computedValueNames[2*i+1]) = computedValues2[i];
+    }
     for (int i = 0; i < force.getNumGlobalParameters(); i++) {
         const string& name = force.getGlobalParameterName(i);
         if (variables.find(name) != variables.end())

platforms/common/include/openmm/common/CommonKernels.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-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -591,7 +591,7 @@ private:
 class CommonCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
 public:
     CommonCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, ComputeContext& cc, const System& system) : CalcCustomNonbondedForceKernel(name, platform),
-            cc(cc), params(NULL), forceCopy(NULL), system(system), hasInitializedKernel(false) {
+            cc(cc), params(NULL), computedValues(NULL), forceCopy(NULL), system(system), hasInitializedKernel(false) {
     }
     ~CommonCalcCustomNonbondedForceKernel();
     /**
@@ -625,13 +625,16 @@ private:
     ComputeContext& cc;
     ForceInfo* info;
     ComputeParameterSet* params;
+    ComputeParameterSet* computedValues;
     ComputeArray globals, interactionGroupData, filteredGroupData, numGroupTiles;
-    ComputeKernel interactionGroupKernel, prepareNeighborListKernel, buildNeighborListKernel;
+    ComputeKernel interactionGroupKernel, prepareNeighborListKernel, buildNeighborListKernel, computedValuesKernel;
     std::vector<void*> interactionGroupArgs;
     std::vector<std::string> globalParamNames;
     std::vector<float> globalParamValues;
     std::vector<ComputeArray> tabulatedFunctionArrays;
     std::map<std::string, const TabulatedFunction*> tabulatedFunctions;
+    std::vector<std::string> paramNames, computedValueNames;
+    std::vector<ComputeParameterInfo> paramBuffers, computedValueBuffers;
     double longRangeCoefficient;
     std::vector<double> longRangeCoefficientDerivs;
     bool hasInitializedLongRangeCorrection, hasInitializedKernel, hasParamDerivs, useNeighborList;

platforms/common/src/CommonKernels.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-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1852,6 +1852,8 @@ CommonCalcCustomNonbondedForceKernel::~CommonCalcCustomNonbondedForceKernel() {
     ContextSelector selector(cc);
     if (params != NULL)
         delete params;
+    if (computedValues != NULL)
+        delete computedValues;
     if (forceCopy != NULL)
         delete forceCopy;
 }
@@ -1868,7 +1870,7 @@ void CommonCalcCustomNonbondedForceKernel::initialize(const System& system, cons
     int numParticles = force.getNumParticles();
     int paddedNumParticles = cc.getPaddedNumAtoms();
     int numParams = force.getNumPerParticleParameters();
-    params = new ComputeParameterSet(cc, numParams, paddedNumParticles, "customNonbondedParameters");
+    params = new ComputeParameterSet(cc, numParams, paddedNumParticles, "customNonbondedParameters", true);
     if (force.getNumGlobalParameters() > 0)
         globals.initialize<float>(cc, force.getNumGlobalParameters(), "customNonbondedGlobals");
     vector<vector<float> > paramVector(paddedNumParticles, vector<float>(numParams, 0));
@@ -1895,6 +1897,7 @@ void CommonCalcCustomNonbondedForceKernel::initialize(const System& system, cons
     vector<pair<string, string> > functionDefinitions;
     vector<const TabulatedFunction*> functionList;
     vector<string> tableTypes;
+    stringstream tableArgs;
     tabulatedFunctionArrays.resize(force.getNumTabulatedFunctions());
     for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
         functionList.push_back(&force.getTabulatedFunction(i));
@@ -1912,6 +1915,10 @@ void CommonCalcCustomNonbondedForceKernel::initialize(const System& system, cons
             tableTypes.push_back("float");
         else
             tableTypes.push_back("float"+cc.intToString(width));
+        tableArgs << ", GLOBAL const float";
+        if (width > 1)
+            tableArgs << width;
+        tableArgs << "* RESTRICT " << arrayName;
     }
 
     // Record information for the expressions.
@@ -1932,6 +1939,26 @@ void CommonCalcCustomNonbondedForceKernel::initialize(const System& system, cons
     forceExpressions["real customEnergy = "] = energyExpression;
     forceExpressions["tempForce -= "] = forceExpression;
 
+    // Record which per-particle parameters and computed values appear in the energy expression.
+
+    if (force.getNumComputedValues() > 0)
+        computedValues = new ComputeParameterSet(cc, force.getNumComputedValues(), paddedNumParticles, "customNonbondedComputedValues", true);
+    for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
+        string name = force.getPerParticleParameterName(i);
+        if (usesVariable(energyExpression, name+"1") || usesVariable(energyExpression, name+"2")) {
+            paramNames.push_back(name);
+            paramBuffers.push_back(params->getParameterInfos()[i]);
+        }
+    }
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        string name, expression;
+        force.getComputedValueParameters(i, name, expression);
+        if (usesVariable(energyExpression, name+"1") || usesVariable(energyExpression, name+"2")) {
+            computedValueNames.push_back(name);
+            computedValueBuffers.push_back(computedValues->getParameterInfos()[i]);
+        }
+    }
+
     // Create the kernels.
 
     vector<pair<ExpressionTreeNode, string> > variables;
@@ -1939,10 +1966,13 @@ void CommonCalcCustomNonbondedForceKernel::initialize(const System& system, cons
     variables.push_back(make_pair(rnode, "r"));
     variables.push_back(make_pair(ExpressionTreeNode(new Operation::Square(), rnode), "r2"));
     variables.push_back(make_pair(ExpressionTreeNode(new Operation::Reciprocal(), rnode), "invR"));
-    for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
-        const string& name = force.getPerParticleParameterName(i);
-        variables.push_back(makeVariable(name+"1", prefix+"params"+params->getParameterSuffix(i, "1")));
-        variables.push_back(makeVariable(name+"2", prefix+"params"+params->getParameterSuffix(i, "2")));
+    for (int i = 0; i < paramNames.size(); i++) {
+        variables.push_back(makeVariable(paramNames[i]+"1", prefix+"params"+cc.intToString(i+1)+"1"));
+        variables.push_back(makeVariable(paramNames[i]+"2", prefix+"params"+cc.intToString(i+1)+"2"));
+    }
+    for (int i = 0; i < computedValueNames.size(); i++) {
+        variables.push_back(makeVariable(computedValueNames[i]+"1", prefix+"values"+cc.intToString(i+1)+"1"));
+        variables.push_back(makeVariable(computedValueNames[i]+"2", prefix+"values"+cc.intToString(i+1)+"2"));
     }
     for (int i = 0; i < force.getNumGlobalParameters(); i++) {
         const string& name = force.getGlobalParameterName(i);
@@ -1973,16 +2003,69 @@ void CommonCalcCustomNonbondedForceKernel::initialize(const System& system, cons
         initInteractionGroups(force, source, tableTypes);
     else {
         cc.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
-        for (int i = 0; i < (int) params->getParameterInfos().size(); i++) {
-            ComputeParameterInfo& parameter = params->getParameterInfos()[i];
-            cc.getNonbondedUtilities().addParameter(ComputeParameterInfo(parameter.getArray(), prefix+"params"+cc.intToString(i+1),
-                    parameter.getComponentType(), parameter.getNumComponents()));
-        }
+        for (int i = 0; i < paramBuffers.size(); i++)
+            cc.getNonbondedUtilities().addParameter(ComputeParameterInfo(paramBuffers[i].getArray(), prefix+"params"+cc.intToString(i+1),
+                    paramBuffers[i].getComponentType(), paramBuffers[i].getNumComponents()));
+        for (int i = 0; i < computedValueBuffers.size(); i++)
+            cc.getNonbondedUtilities().addParameter(ComputeParameterInfo(computedValueBuffers[i].getArray(), prefix+"values"+cc.intToString(i+1),
+                    computedValueBuffers[i].getComponentType(), computedValueBuffers[i].getNumComponents()));
         if (globals.isInitialized()) {
             globals.upload(globalParamValues);
             cc.getNonbondedUtilities().addArgument(ComputeParameterInfo(globals, prefix+"globals", "float", 1));
         }
     }
+    if (force.getNumComputedValues() > 0) {
+        // Create the kernel to calculate computed values.
+
+        stringstream valuesSource, args;
+        for (int i = 0; i < computedValues->getParameterInfos().size(); i++) {
+            ComputeParameterInfo& buffer = computedValues->getParameterInfos()[i];
+            string valueName = "values"+cc.intToString(i+1);
+            if (i > 0)
+                args << ", ";
+            args << "GLOBAL " << buffer.getType() << "* RESTRICT global_" << valueName;
+            valuesSource << buffer.getType() << " local_" << valueName << ";\n";
+        }
+        if (force.getNumGlobalParameters() > 0)
+            args << ", GLOBAL const float* globals";
+        for (int i = 0; i < params->getParameterInfos().size(); i++) {
+            ComputeParameterInfo& buffer = params->getParameterInfos()[i];
+            string paramName = "params"+cc.intToString(i+1);
+            args << ", GLOBAL const " << buffer.getType() << "* RESTRICT " << paramName;
+        }
+        map<string, string> variables;
+        for (int i = 0; i < force.getNumPerParticleParameters(); i++)
+            variables[force.getPerParticleParameterName(i)] = "params"+params->getParameterSuffix(i, "[index]");
+        for (int i = 0; i < force.getNumGlobalParameters(); i++)
+            variables[force.getGlobalParameterName(i)] = "globals["+cc.intToString(i)+"]";
+        for (int i = 0; i < force.getNumComputedValues(); i++) {
+            string name, expression;
+            force.getComputedValueParameters(i, name, expression);
+            variables[name] = "local_values"+computedValues->getParameterSuffix(i);
+            map<string, Lepton::ParsedExpression> valueExpressions;
+            valueExpressions["local_values"+computedValues->getParameterSuffix(i)+" = "] = Lepton::Parser::parse(expression, functions).optimize();
+            valuesSource << cc.getExpressionUtilities().createExpressions(valueExpressions, variables, functionList, functionDefinitions, "value"+cc.intToString(i)+"_temp");
+        }
+        for (int i = 0; i < (int) computedValues->getParameterInfos().size(); i++) {
+            string valueName = "values"+cc.intToString(i+1);
+            valuesSource << "global_" << valueName << "[index] = local_" << valueName << ";\n";
+        }
+        map<string, string> replacements;
+        replacements["PARAMETER_ARGUMENTS"] = args.str()+tableArgs.str();
+        replacements["COMPUTE_VALUES"] = valuesSource.str();
+        map<string, string> defines;
+        defines["NUM_ATOMS"] = cc.intToString(cc.getNumAtoms());
+        ComputeProgram program = cc.compileProgram(cc.replaceStrings(CommonKernelSources::customNonbondedComputedValues, replacements), defines);
+        computedValuesKernel = program->createKernel("computePerParticleValues");
+        for (auto& value : computedValues->getParameterInfos())
+            computedValuesKernel->addArg(value.getArray());
+        if (globals.isInitialized())
+            computedValuesKernel->addArg(globals);
+        for (auto& parameter : params->getParameterInfos())
+            computedValuesKernel->addArg(parameter.getArray());
+        for (auto& function : tabulatedFunctionArrays)
+            computedValuesKernel->addArg(function);
+    }
     info = new ForceInfo(force);
     cc.addForce(info);
     
@@ -2180,21 +2263,21 @@ void CommonCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
     const string suffixes[] = {"x", "y", "z", "w"};
     stringstream localData;
     int localDataSize = 0;
-    vector<ComputeParameterInfo>& buffers = params->getParameterInfos(); 
-    for (int i = 0; i < (int) buffers.size(); i++) {
-        if (buffers[i].getNumComponents() == 1)
-            localData<<buffers[i].getComponentType()<<" params"<<(i+1)<<";\n";
-        else {
-            for (int j = 0; j < buffers[i].getNumComponents(); ++j)
-                localData<<buffers[i].getComponentType()<<" params"<<(i+1)<<"_"<<suffixes[j]<<";\n";
-        }
-        localDataSize += buffers[i].getSize();
+    for (int i = 0; i < paramBuffers.size(); i++) {
+        localData<<paramBuffers[i].getComponentType()<<" params"<<(i+1)<<";\n";
+        localDataSize += paramBuffers[i].getSize();
+    }
+    for (int i = 0; i < computedValueBuffers.size(); i++) {
+        localData<<computedValueBuffers[i].getComponentType()<<" values"<<(i+1)<<";\n";
+        localDataSize += computedValueBuffers[i].getSize();
     }
     replacements["ATOM_PARAMETER_DATA"] = localData.str();
     stringstream args;
-    for (int i = 0; i < (int) buffers.size(); i++)
-        args<<", GLOBAL const "<<buffers[i].getType()<<"* RESTRICT global_params"<<(i+1);
-    for (int i = 0; i < (int) tabulatedFunctionArrays.size(); i++)
+    for (int i = 0; i < paramBuffers.size(); i++)
+        args<<", GLOBAL const "<<paramBuffers[i].getType()<<"* RESTRICT global_params"<<(i+1);
+    for (int i = 0; i < computedValueBuffers.size(); i++)
+        args<<", GLOBAL const "<<computedValueBuffers[i].getType()<<"* RESTRICT global_values"<<(i+1);
+    for (int i = 0; i < tabulatedFunctionArrays.size(); i++)
         args << ", GLOBAL const " << tableTypes[i]<< "* RESTRICT table" << i;
     if (globals.isInitialized())
         args<<", GLOBAL const float* RESTRICT globals";
@@ -2202,34 +2285,22 @@ void CommonCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNon
         args << ", GLOBAL mixed* RESTRICT energyParamDerivs";
     replacements["PARAMETER_ARGUMENTS"] = args.str();
     stringstream load1;
-    for (int i = 0; i < (int) buffers.size(); i++)
-        load1<<buffers[i].getType()<<" params"<<(i+1)<<"1 = global_params"<<(i+1)<<"[atom1];\n";
+    for (int i = 0; i < paramBuffers.size(); i++)
+        load1<<paramBuffers[i].getType()<<" params"<<(i+1)<<"1 = global_params"<<(i+1)<<"[atom1];\n";
+    for (int i = 0; i < computedValueBuffers.size(); i++)
+        load1<<computedValueBuffers[i].getType()<<" values"<<(i+1)<<"1 = global_values"<<(i+1)<<"[atom1];\n";
     replacements["LOAD_ATOM1_PARAMETERS"] = load1.str();
     stringstream loadLocal2;
-    for (int i = 0; i < (int) buffers.size(); i++) {
-        if (buffers[i].getNumComponents() == 1)
-            loadLocal2<<"localData[LOCAL_ID].params"<<(i+1)<<" = global_params"<<(i+1)<<"[atom2];\n";
-        else {
-            loadLocal2<<buffers[i].getType()<<" temp_params"<<(i+1)<<" = global_params"<<(i+1)<<"[atom2];\n";
-            for (int j = 0; j < buffers[i].getNumComponents(); ++j)
-                loadLocal2<<"localData[LOCAL_ID].params"<<(i+1)<<"_"<<suffixes[j]<<" = temp_params"<<(i+1)<<"."<<suffixes[j]<<";\n";
-        }
-    }
+    for (int i = 0; i < paramBuffers.size(); i++)
+        loadLocal2<<"localData[LOCAL_ID].params"<<(i+1)<<" = global_params"<<(i+1)<<"[atom2];\n";
+    for (int i = 0; i < computedValueBuffers.size(); i++)
+        loadLocal2<<"localData[LOCAL_ID].values"<<(i+1)<<" = global_values"<<(i+1)<<"[atom2];\n";
     replacements["LOAD_LOCAL_PARAMETERS"] = loadLocal2.str();
     stringstream load2;
-    for (int i = 0; i < (int) buffers.size(); i++) {
-        if (buffers[i].getNumComponents() == 1)
-            load2<<buffers[i].getType()<<" params"<<(i+1)<<"2 = localData[localIndex].params"<<(i+1)<<";\n";
-        else {
-            load2<<buffers[i].getType()<<" params"<<(i+1)<<"2 = make_"<<buffers[i].getType()<<"(";
-            for (int j = 0; j < buffers[i].getNumComponents(); ++j) {
-                if (j > 0)
-                    load2<<", ";
-                load2<<"localData[localIndex].params"<<(i+1)<<"_"<<suffixes[j];
-            }
-            load2<<");\n";
-        }
-    }
+    for (int i = 0; i < paramBuffers.size(); i++)
+        load2<<paramBuffers[i].getType()<<" params"<<(i+1)<<"2 = localData[localIndex].params"<<(i+1)<<";\n";
+    for (int i = 0; i < computedValueBuffers.size(); i++)
+        load2<<computedValueBuffers[i].getType()<<" values"<<(i+1)<<"2 = localData[localIndex].values"<<(i+1)<<";\n";
     replacements["LOAD_ATOM2_PARAMETERS"] = load2.str();
     stringstream initDerivs, saveDerivs;
     const vector<string>& allParamDerivNames = cc.getEnergyParamDerivNames();
@@ -2303,6 +2374,8 @@ double CommonCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool
         else
             hasInitializedLongRangeCorrection = false;
     }
+    if (computedValues != NULL)
+        computedValuesKernel->execute(cc.getNumAtoms());
     if (interactionGroupData.isInitialized()) {
         if (!hasInitializedKernel) {
             hasInitializedKernel = true;
@@ -2315,8 +2388,10 @@ double CommonCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool
             interactionGroupKernel->addArg((int) useNeighborList);
             for (int i = 0; i < 5; i++)
                 interactionGroupKernel->addArg(); // Periodic box information will be set just before it is executed.
-            for (auto& parameter : params->getParameterInfos())
-                interactionGroupKernel->addArg(parameter.getArray());
+            for (auto& buffer : paramBuffers)
+                interactionGroupKernel->addArg(buffer.getArray());
+            for (auto& buffer : computedValueBuffers)
+                interactionGroupKernel->addArg(buffer.getArray());
             for (auto& function : tabulatedFunctionArrays)
                 interactionGroupKernel->addArg(function);
             if (globals.isInitialized())

platforms/common/src/kernels/customNonbondedComputedValues.cc

+/**
+ * Calculate per-particle computed values for a CustomNonbondedForce.
+ */
+
+KERNEL void computePerParticleValues(PARAMETER_ARGUMENTS) {
+    for (int index = GLOBAL_ID; index < NUM_ATOMS; index += GLOBAL_SIZE) {
+        COMPUTE_VALUES
+    }
+}

platforms/cpu/include/CpuCustomNonbondedForce.h

 
-/* Portions copyright (c) 2009-2018 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2022 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -49,7 +49,9 @@ class CpuCustomNonbondedForce {
 
        CpuCustomNonbondedForce(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
                                const std::vector<std::string>& parameterNames, const std::vector<std::set<int> >& exclusions,
-                               const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions, ThreadPool& threads);
+                               const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions,
+                               const std::vector<std::string>& computedValueNames, const std::vector<Lepton::CompiledExpression> computedValueExpressions,
+                               ThreadPool& threads);
 
       /**---------------------------------------------------------------------------------------
 
@@ -137,9 +139,10 @@ private:
     ThreadPool& threads;
     const std::vector<std::set<int> > exclusions;
     std::vector<ThreadData*> threadData;
-    std::vector<std::string> paramNames;
+    std::vector<std::string> paramNames, computedValueNames;
     std::vector<std::pair<int, int> > groupInteractions;
     std::vector<double> threadEnergy;
+    std::vector<std::vector<double> > atomComputedValues;
     // The following variables are used to make information accessible to the individual threads.
     int numberOfAtoms;
     float* posq;
@@ -178,14 +181,16 @@ private:
 class CpuCustomNonbondedForce::ThreadData {
 public:
     ThreadData(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression, const std::vector<std::string>& parameterNames,
-            const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions);
+            const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions, const std::vector<std::string>& computedValueNames,
+            const std::vector<Lepton::CompiledExpression> computedValueExpressions, std::vector<std::vector<double> >& atomComputedValues);
     Lepton::CompiledExpression energyExpression;
     Lepton::CompiledExpression forceExpression;
-    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+    std::vector<Lepton::CompiledExpression> computedValueExpressions, energyParamDerivExpressions;
     CompiledExpressionSet expressionSet;
-    std::vector<double> particleParam;
+    std::vector<double> particleParam, computedValues;
     double r;
     std::vector<double> energyParamDerivs; 
+    std::vector<std::vector<double> >& atomComputedValues;
 };
 
 } // namespace OpenMM

platforms/cpu/include/CpuKernels.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) 2013-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -331,7 +331,7 @@ private:
     CustomNonbondedForceImpl::LongRangeCorrectionData longRangeCorrectionData;
     std::map<std::string, double> globalParamValues;
     std::vector<std::set<int> > exclusions;
-    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
+    std::vector<std::string> parameterNames, globalParameterNames, computedValueNames, energyParamDerivNames;
     std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
     std::vector<double> longRangeCoefficientDerivs;
     std::map<std::string, const TabulatedFunction*> tabulatedFunctions;

platforms/cpu/src/CpuCustomNonbondedForce.cpp

 
-/* Portions copyright (c) 2009-2018 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2022 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -33,11 +33,17 @@ using namespace OpenMM;
 using namespace std;
 
 CpuCustomNonbondedForce::ThreadData::ThreadData(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
-            const vector<string>& parameterNames, const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions) :
-            energyExpression(energyExpression), forceExpression(forceExpression), energyParamDerivExpressions(energyParamDerivExpressions) {
+            const vector<string>& parameterNames, const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions,
+            const vector<string>& computedValueNames, const vector<Lepton::CompiledExpression> computedValueExpressions,
+            vector<vector<double> >& atomComputedValues) :
+            energyExpression(energyExpression), forceExpression(forceExpression), energyParamDerivExpressions(energyParamDerivExpressions),
+            computedValueExpressions(computedValueExpressions), atomComputedValues(atomComputedValues) {
+    // Prepare for passing variables to expressions.
+
     map<string, double*> variableLocations;
     variableLocations["r"] = &r;
     particleParam.resize(2*parameterNames.size());
+    computedValues.resize(2*computedValueNames.size());
     for (int i = 0; i < (int) parameterNames.size(); i++) {
         for (int j = 0; j < 2; j++) {
             stringstream name;
@@ -45,6 +51,13 @@ CpuCustomNonbondedForce::ThreadData::ThreadData(const Lepton::CompiledExpression
             variableLocations[name.str()] = &particleParam[i*2+j];
         }
     }
+    for (int i = 0; i < (int) computedValueNames.size(); i++) {
+        for (int j = 0; j < 2; j++) {
+            stringstream name;
+            name << computedValueNames[i] << (j+1);
+            variableLocations[name.str()] = &computedValues[i*2+j];
+        }
+    }
     energyParamDerivs.resize(energyParamDerivExpressions.size());
     this->energyExpression.setVariableLocations(variableLocations);
     this->forceExpression.setVariableLocations(variableLocations);
@@ -54,14 +67,26 @@ CpuCustomNonbondedForce::ThreadData::ThreadData(const Lepton::CompiledExpression
         expression.setVariableLocations(variableLocations);
         expressionSet.registerExpression(expression);
     }
+
+    // Prepare for passing variables to the computed value expressions.
+
+    map<string, double*> valueVariableLocations;
+    for (int i = 0; i < parameterNames.size(); i++)
+        valueVariableLocations[parameterNames[i]] = &particleParam[i];
+    for (auto& expression : this->computedValueExpressions) {
+        expression.setVariableLocations(valueVariableLocations);
+        expressionSet.registerExpression(expression);
+    }
 }
 
 CpuCustomNonbondedForce::CpuCustomNonbondedForce(const Lepton::CompiledExpression& energyExpression,
             const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, const vector<set<int> >& exclusions,
-            const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions, ThreadPool& threads) :
-            cutoff(false), useSwitch(false), periodic(false), useInteractionGroups(false), paramNames(parameterNames), exclusions(exclusions), threads(threads) {
+            const vector<Lepton::CompiledExpression> energyParamDerivExpressions, const vector<string>& computedValueNames,
+            const vector<Lepton::CompiledExpression> computedValueExpressions, ThreadPool& threads) :
+            cutoff(false), useSwitch(false), periodic(false), useInteractionGroups(false), paramNames(parameterNames), exclusions(exclusions),
+            computedValueNames(computedValueNames), threads(threads) {
     for (int i = 0; i < threads.getNumThreads(); i++)
-        threadData.push_back(new ThreadData(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions));
+        threadData.push_back(new ThreadData(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions, computedValueNames, computedValueExpressions, atomComputedValues));
 }
 
 CpuCustomNonbondedForce::~CpuCustomNonbondedForce() {
@@ -133,12 +158,15 @@ void CpuCustomNonbondedForce::calculatePairIxn(int numberOfAtoms, float* posq, v
     this->includeForce = includeForce;
     this->includeEnergy = includeEnergy;
     threadEnergy.resize(threads.getNumThreads());
+    atomComputedValues.resize(computedValueNames.size(), vector<double>(numberOfAtoms));
     atomicCounter = 0;
     
     // Signal the threads to start running and wait for them to finish.
     
     threads.execute([&] (ThreadPool& threads, int threadIndex) { threadComputeForce(threads, threadIndex); });
-    threads.waitForThreads();
+    threads.waitForThreads(); // Computed values
+    threads.resumeThreads();
+    threads.waitForThreads(); // Interactions
     
     // Combine the energies from all the threads.
     
@@ -157,15 +185,31 @@ void CpuCustomNonbondedForce::calculatePairIxn(int numberOfAtoms, float* posq, v
 }
 
 void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int threadIndex) {
+    int numThreads = threads.getNumThreads();
+    ThreadData& data = *threadData[threadIndex];
+    for (auto& param : *globalParameters)
+        data.expressionSet.setVariable(data.expressionSet.getVariableIndex(param.first), param.second);
+
+    // Process computed values for this thread's subset of interactions.
+
+    int numComputedValues = atomComputedValues.size();
+    if (numComputedValues > 0) {
+        int start = threadIndex*numberOfAtoms/numThreads;
+        int end = (threadIndex+1)*numberOfAtoms/numThreads;
+        for (int i = start; i < end; i++) {
+            for (int j = 0; j < paramNames.size(); j++)
+                data.particleParam[j] = atomParameters[i][j];
+            for (int j = 0; j < numComputedValues; j++)
+                atomComputedValues[j][i] = data.computedValueExpressions[j].evaluate();
+        }
+    }
+    threads.syncThreads();
+
     // Compute this thread's subset of interactions.
 
-    int numThreads = threads.getNumThreads();
     threadEnergy[threadIndex] = 0;
     double& energy = threadEnergy[threadIndex];
     float* forces = &(*threadForce)[threadIndex][0];
-    ThreadData& data = *threadData[threadIndex];
-    for (auto& param : *globalParameters)
-        data.expressionSet.setVariable(data.expressionSet.getVariableIndex(param.first), param.second);
     for (auto& deriv : data.energyParamDerivs)
         deriv = 0.0;
     fvec4 boxSize(periodicBoxVectors[0][0], periodicBoxVectors[1][1], periodicBoxVectors[2][2], 0);
@@ -178,10 +222,14 @@ void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int thread
         for (int i = start; i < end; i++) {
             int atom1 = groupInteractions[i].first;
             int atom2 = groupInteractions[i].second;
-            for (int j = 0; j < (int) paramNames.size(); j++) {
+            for (int j = 0; j < paramNames.size(); j++) {
                 data.particleParam[j*2] = atomParameters[atom1][j];
                 data.particleParam[j*2+1] = atomParameters[atom2][j];
             }
+            for (int j = 0; j < computedValueNames.size(); j++) {
+                data.computedValues[j*2] = atomComputedValues[j][atom1];
+                data.computedValues[j*2+1] = atomComputedValues[j][atom2];
+            }
             calculateOneIxn(atom1, atom2, data, forces, energy, boxSize, invBoxSize);
         }
     }
@@ -198,13 +246,17 @@ void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int thread
             const auto& exclusions = neighborList->getBlockExclusions(blockIndex);
             for (int i = 0; i < (int) neighbors.size(); i++) {
                 int first = neighbors[i];
-                for (int j = 0; j < (int) paramNames.size(); j++)
+                for (int j = 0; j < paramNames.size(); j++)
                     data.particleParam[j*2] = atomParameters[first][j];
+                for (int j = 0; j < computedValueNames.size(); j++)
+                    data.computedValues[j*2] = atomComputedValues[j][first];
                 for (int k = 0; k < blockSize; k++) {
                     if ((exclusions[i] & (1<<k)) == 0) {
                         int second = blockAtom[k];
-                        for (int j = 0; j < (int) paramNames.size(); j++)
+                        for (int j = 0; j < paramNames.size(); j++)
                             data.particleParam[j*2+1] = atomParameters[second][j];
+                        for (int j = 0; j < computedValueNames.size(); j++)
+                            data.computedValues[j*2+1] = atomComputedValues[j][second];
                         calculateOneIxn(first, second, data, forces, energy, boxSize, invBoxSize);
                     }
                 }
@@ -220,10 +272,14 @@ void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int thread
                 break;
             for (int jj = ii+1; jj < numberOfAtoms; jj++) {
                 if (exclusions[jj].find(ii) == exclusions[jj].end()) {
-                    for (int j = 0; j < (int) paramNames.size(); j++) {
+                    for (int j = 0; j < paramNames.size(); j++) {
                         data.particleParam[j*2] = atomParameters[ii][j];
                         data.particleParam[j*2+1] = atomParameters[jj][j];
                     }
+                    for (int j = 0; j < computedValueNames.size(); j++) {
+                        data.computedValues[j*2] = atomComputedValues[j][ii];
+                        data.computedValues[j*2+1] = atomComputedValues[j][jj];
+                    }
                     calculateOneIxn(ii, jj, data, forces, energy, boxSize, invBoxSize);
                 }
             }

platforms/cpu/src/CpuKernels.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-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -930,20 +930,35 @@ void CpuCalcCustomNonbondedForceKernel::createInteraction(const CustomNonbondedF
         globalParameterNames.push_back(force.getGlobalParameterName(i));
         globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
     }
-    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+    set<string> particleVariables, pairVariables;
+    particleVariables.insert("r");
+    pairVariables.insert("r");
+    for (auto& name : parameterNames) {
+        particleVariables.insert(name);
+        pairVariables.insert(name+"1");
+        pairVariables.insert(name+"2");
+    }
+    particleVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    pairVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    vector<Lepton::CompiledExpression> computedValueExpressions, energyParamDerivExpressions;
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        string name, exp;
+        force.getComputedValueParameters(i, name, exp);
+        Lepton::ParsedExpression parsed = Lepton::Parser::parse(exp, functions);
+        validateVariables(parsed.getRootNode(), particleVariables);
+        computedValueNames.push_back(name);
+        computedValueExpressions.push_back(parsed.createCompiledExpression());
+    }
     for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
         string param = force.getEnergyParameterDerivativeName(i);
         energyParamDerivNames.push_back(param);
         energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
     }
-    set<string> variables;
-    variables.insert("r");
-    for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
-        variables.insert(parameterNames[i]+"1");
-        variables.insert(parameterNames[i]+"2");
+    for (auto& name : computedValueNames) {
+        pairVariables.insert(name+"1");
+        pairVariables.insert(name+"2");
     }
-    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
-    validateVariables(expression.getRootNode(), variables);
+    validateVariables(expression.getRootNode(), pairVariables);
 
     // Delete the custom functions.
 
@@ -952,7 +967,8 @@ void CpuCalcCustomNonbondedForceKernel::createInteraction(const CustomNonbondedF
 
     // Create the object that computes the interaction.
 
-    nonbonded = new CpuCustomNonbondedForce(energyExpression, forceExpression, parameterNames, exclusions, energyParamDerivExpressions, data.threads);
+    nonbonded = new CpuCustomNonbondedForce(energyExpression, forceExpression, parameterNames, exclusions, energyParamDerivExpressions,
+            computedValueNames, computedValueExpressions, data.threads);
     if (interactionGroups.size() > 0)
         nonbonded->setInteractionGroups(interactionGroups);
 }

platforms/reference/include/ReferenceCustomNonbondedIxn.h

 
-/* Portions copyright (c) 2009-2018 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2022 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -47,10 +47,10 @@ class ReferenceCustomNonbondedIxn {
       double cutoffDistance, switchingDistance;
       Lepton::CompiledExpression energyExpression;
       Lepton::CompiledExpression forceExpression;
-      std::vector<std::string> paramNames;
-      std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
+      std::vector<std::string> paramNames, computedValueNames;
+      std::vector<Lepton::CompiledExpression> computedValueExpressions, energyParamDerivExpressions;
       CompiledExpressionSet expressionSet;
-      std::vector<int> particleParamIndex;
+      std::vector<int> particleParamIndex, computedValueIndex;
       int rIndex;
       std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
 
@@ -79,7 +79,8 @@ class ReferenceCustomNonbondedIxn {
          --------------------------------------------------------------------------------------- */
 
        ReferenceCustomNonbondedIxn(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
-                                   const std::vector<std::string>& parameterNames, const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions);
+                                   const std::vector<std::string>& parameterNames, const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions,
+                                   const std::vector<std::string>& computedValueNames, const std::vector<Lepton::CompiledExpression> computedValueExpressions);
 
       /**---------------------------------------------------------------------------------------
 

platforms/reference/include/ReferenceKernels.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-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -692,8 +692,8 @@ private:
     std::map<std::string, double> globalParamValues;
     std::vector<std::set<int> > exclusions;
     Lepton::CompiledExpression energyExpression, forceExpression;
-    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
-    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
+    std::vector<Lepton::CompiledExpression> computedValueExpressions, energyParamDerivExpressions;
+    std::vector<std::string> parameterNames, globalParameterNames, computedValueNames, energyParamDerivNames;
     std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
     std::vector<double> longRangeCoefficientDerivs;
     std::map<std::string, const TabulatedFunction*> tabulatedFunctions;

platforms/reference/src/ReferenceKernels.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-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1211,6 +1211,8 @@ void ReferenceCalcCustomNonbondedForceKernel::createExpressions(const CustomNonb
     parameterNames.clear();
     globalParameterNames.clear();
     globalParamValues.clear();
+    computedValueNames.clear();
+    computedValueExpressions.clear();
     energyParamDerivNames.clear();
     energyParamDerivExpressions.clear();
     for (int i = 0; i < force.getNumPerParticleParameters(); i++)
@@ -1219,19 +1221,34 @@ void ReferenceCalcCustomNonbondedForceKernel::createExpressions(const CustomNonb
         globalParameterNames.push_back(force.getGlobalParameterName(i));
         globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
     }
+    set<string> particleVariables, pairVariables;
+    particleVariables.insert("r");
+    pairVariables.insert("r");
+    for (auto& name : parameterNames) {
+        particleVariables.insert(name);
+        pairVariables.insert(name+"1");
+        pairVariables.insert(name+"2");
+    }
+    particleVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    pairVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        string name, exp;
+        force.getComputedValueParameters(i, name, exp);
+        Lepton::ParsedExpression parsed = Lepton::Parser::parse(exp, functions);
+        validateVariables(parsed.getRootNode(), particleVariables);
+        computedValueNames.push_back(name);
+        computedValueExpressions.push_back(parsed.createCompiledExpression());
+    }
     for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
         string param = force.getEnergyParameterDerivativeName(i);
         energyParamDerivNames.push_back(param);
         energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
     }
-    set<string> variables;
-    variables.insert("r");
-    for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
-        variables.insert(parameterNames[i]+"1");
-        variables.insert(parameterNames[i]+"2");
+    for (auto& name : computedValueNames) {
+        pairVariables.insert(name+"1");
+        pairVariables.insert(name+"2");
     }
-    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
-    validateVariables(expression.getRootNode(), variables);
+    validateVariables(expression.getRootNode(), pairVariables);
 
     // Delete the custom functions.
 
@@ -1244,7 +1261,7 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
     vector<Vec3>& forceData = extractForces(context);
     Vec3* boxVectors = extractBoxVectors(context);
     double energy = 0;
-    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions);
+    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions, computedValueNames, computedValueExpressions);
     bool periodic = (nonbondedMethod == CutoffPeriodic);
     if (nonbondedMethod != NoCutoff) {
         computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp

 
-/* Portions copyright (c) 2009-2018 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2022 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -45,13 +45,17 @@ using namespace OpenMM;
 
 ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::CompiledExpression& energyExpression,
         const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames,
-        const vector<Lepton::CompiledExpression> energyParamDerivExpressions) :
+        const vector<Lepton::CompiledExpression> energyParamDerivExpressions,
+        const std::vector<std::string>& computedValueNames, const std::vector<Lepton::CompiledExpression> computedValueExpressions) :
             cutoff(false), useSwitch(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression),
-            paramNames(parameterNames), energyParamDerivExpressions(energyParamDerivExpressions) {
+            paramNames(parameterNames), energyParamDerivExpressions(energyParamDerivExpressions), computedValueNames(computedValueNames),
+            computedValueExpressions(computedValueExpressions) {
     expressionSet.registerExpression(this->energyExpression);
     expressionSet.registerExpression(this->forceExpression);
     for (int i = 0; i < this->energyParamDerivExpressions.size(); i++)
         expressionSet.registerExpression(this->energyParamDerivExpressions[i]);
+    for (int i = 0; i < this->computedValueExpressions.size(); i++)
+        expressionSet.registerExpression(this->computedValueExpressions[i]);
     rIndex = expressionSet.getVariableIndex("r");
     for (auto& param : paramNames) {
         for (int j = 1; j < 3; j++) {
@@ -60,6 +64,13 @@ ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::CompiledE
             particleParamIndex.push_back(expressionSet.getVariableIndex(name.str()));
         }
     }
+    for (auto& value : computedValueNames) {
+        for (int j = 1; j < 3; j++) {
+            stringstream name;
+            name << value << j;
+            computedValueIndex.push_back(expressionSet.getVariableIndex(name.str()));
+        }
+    }
 }
 
 /**---------------------------------------------------------------------------------------
@@ -159,7 +170,25 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn(int numberOfAtoms, vector<Vec
 
     for (auto& param : globalParameters)
         expressionSet.setVariable(expressionSet.getVariableIndex(param.first), param.second);
-    if (interactionGroups.size() > 0) { 
+    
+    // Calculate computed values.
+    
+    vector<int> paramIndex, valueIndex;
+    for (auto& name : paramNames)
+        paramIndex.push_back(expressionSet.getVariableIndex(name));
+    for (auto& name : computedValueNames)
+        valueIndex.push_back(expressionSet.getVariableIndex(name));
+    vector<vector<double> > computedValues(computedValueNames.size(), vector<double>(numberOfAtoms));
+    for (int i = 0; i < numberOfAtoms; i++) {
+        for (int j = 0; j < paramNames.size(); j++)
+            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
+        for (int j = 0; j < computedValueNames.size(); j++)
+            computedValues[j][i] = computedValueExpressions[j].evaluate();
+    }
+    
+    // Compute the interactions.
+
+    if (interactionGroups.size() > 0) {
         // The user has specified interaction groups, so compute only the requested interactions.
         
         for (auto& group : interactionGroups) {
@@ -175,6 +204,10 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn(int numberOfAtoms, vector<Vec
                         expressionSet.setVariable(particleParamIndex[j*2], atomParameters[*atom1][j]);
                         expressionSet.setVariable(particleParamIndex[j*2+1], atomParameters[*atom2][j]);
                     }
+                    for (int j = 0; j < (int) computedValueNames.size(); j++) {
+                        expressionSet.setVariable(computedValueIndex[j*2], computedValues[j][*atom1]);
+                        expressionSet.setVariable(computedValueIndex[j*2+1], computedValues[j][*atom2]);
+                    }
                     calculateOneIxn(*atom1, *atom2, atomCoordinates, forces, totalEnergy, energyParamDerivs);
                 }
             }
@@ -188,6 +221,10 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn(int numberOfAtoms, vector<Vec
                 expressionSet.setVariable(particleParamIndex[j*2], atomParameters[pair.first][j]);
                 expressionSet.setVariable(particleParamIndex[j*2+1], atomParameters[pair.second][j]);
             }
+            for (int j = 0; j < (int) computedValueNames.size(); j++) {
+                expressionSet.setVariable(computedValueIndex[j*2], computedValues[j][pair.first]);
+                expressionSet.setVariable(computedValueIndex[j*2+1], computedValues[j][pair.second]);
+            }
             calculateOneIxn(pair.first, pair.second, atomCoordinates, forces, totalEnergy, energyParamDerivs);
         }
     }
@@ -201,6 +238,10 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn(int numberOfAtoms, vector<Vec
                         expressionSet.setVariable(particleParamIndex[j*2], atomParameters[ii][j]);
                         expressionSet.setVariable(particleParamIndex[j*2+1], atomParameters[jj][j]);
                     }
+                    for (int j = 0; j < (int) computedValueNames.size(); j++) {
+                        expressionSet.setVariable(computedValueIndex[j*2], computedValues[j][ii]);
+                        expressionSet.setVariable(computedValueIndex[j*2+1], computedValues[j][jj]);
+                    }
                     calculateOneIxn(ii, jj, atomCoordinates, forces, totalEnergy, energyParamDerivs);
                 }
             }

serialization/src/CustomNonbondedForceProxy.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-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -42,7 +42,7 @@ CustomNonbondedForceProxy::CustomNonbondedForceProxy() : SerializationProxy("Cus
 }
 
 void CustomNonbondedForceProxy::serialize(const void* object, SerializationNode& node) const {
-    node.setIntProperty("version", 2);
+    node.setIntProperty("version", 3);
     const CustomNonbondedForce& force = *reinterpret_cast<const CustomNonbondedForce*>(object);
     node.setIntProperty("forceGroup", force.getForceGroup());
     node.setStringProperty("name", force.getName());
@@ -60,6 +60,12 @@ void CustomNonbondedForceProxy::serialize(const void* object, SerializationNode&
     for (int i = 0; i < force.getNumGlobalParameters(); i++) {
         globalParams.createChildNode("Parameter").setStringProperty("name", force.getGlobalParameterName(i)).setDoubleProperty("default", force.getGlobalParameterDefaultValue(i));
     }
+    SerializationNode& computedValues = node.createChildNode("ComputedValues");
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        string name, expression;
+        force.getComputedValueParameters(i, name, expression);
+        computedValues.createChildNode("Value").setStringProperty("name", name).setStringProperty("expression", expression);
+    }
     SerializationNode& energyDerivs = node.createChildNode("EnergyParameterDerivatives");
     for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
         energyDerivs.createChildNode("Parameter").setStringProperty("name", force.getEnergyParameterDerivativeName(i));
@@ -103,7 +109,7 @@ void CustomNonbondedForceProxy::serialize(const void* object, SerializationNode&
 
 void* CustomNonbondedForceProxy::deserialize(const SerializationNode& node) const {
     int version = node.getIntProperty("version");
-    if (version < 1 || version > 2)
+    if (version < 1 || version > 3)
         throw OpenMMException("Unsupported version number");
     CustomNonbondedForce* force = NULL;
     try {
@@ -121,6 +127,11 @@ void* CustomNonbondedForceProxy::deserialize(const SerializationNode& node) cons
         const SerializationNode& globalParams = node.getChildNode("GlobalParameters");
         for (auto& parameter : globalParams.getChildren())
             force->addGlobalParameter(parameter.getStringProperty("name"), parameter.getDoubleProperty("default"));
+        if (version > 2) {
+            const SerializationNode& computedValues = node.getChildNode("ComputedValues");
+            for (auto& value : computedValues.getChildren())
+                force->addComputedValue(value.getStringProperty("name"), value.getStringProperty("expression"));
+        }
         if (version > 1) {
             const SerializationNode& energyDerivs = node.getChildNode("EnergyParameterDerivatives");
             for (auto& parameter : energyDerivs.getChildren())

serialization/tests/TestSerializeCustomNonbondedForce.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-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -53,6 +53,8 @@ void testSerialization() {
     force.addGlobalParameter("y", 2.221);
     force.addPerParticleParameter("z");
     force.addEnergyParameterDerivative("y");
+    force.addComputedValue("c1", "x+y");
+    force.addComputedValue("c2", "x*z+2");
     vector<double> params(1);
     params[0] = 1.0;
     force.addParticle(params);
@@ -97,6 +99,14 @@ void testSerialization() {
         ASSERT_EQUAL(force.getGlobalParameterName(i), force2.getGlobalParameterName(i));
         ASSERT_EQUAL(force.getGlobalParameterDefaultValue(i), force2.getGlobalParameterDefaultValue(i));
     }
+    ASSERT_EQUAL(force.getNumComputedValues(), force2.getNumComputedValues());
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        string name1, name2, exp1, exp2;
+        force.getComputedValueParameters(i, name1, exp1);
+        force2.getComputedValueParameters(i, name2, exp2);
+        ASSERT_EQUAL(name1, name2);
+        ASSERT_EQUAL(exp1, exp2)
+    }
     ASSERT_EQUAL(force.getNumEnergyParameterDerivatives(), force2.getNumEnergyParameterDerivatives());
     for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++)
         ASSERT_EQUAL(force.getEnergyParameterDerivativeName(i), force2.getEnergyParameterDerivativeName(i));

tests/TestCustomNonbondedForce.h

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2018 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2022 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1445,6 +1445,69 @@ void testEnergyParameterDerivativesWithGroups() {
     }
 }
 
+void testComputedValues(int mode) {
+    // Create a box of particles.
+
+    int gridSize = 5;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = gridSize*0.7;
+    double cutoff = boxSize/3;
+    System system;
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* nb1 = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2); sigma1=(1-lambda)*sigmaA1+lambda*sigmaB1; sigma2=(1-lambda)*sigmaA2+lambda*sigmaB2");
+    CustomNonbondedForce* nb2 = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
+    nb2->addComputedValue("sigma", "(1-lambda)*sigmaA+lambda*sigmaB");
+    for (int model = 0; model < 2; model++) {
+        CustomNonbondedForce* force = (model == 0 ? nb1 : nb2);
+        force->addGlobalParameter("lambda", 0);
+        force->addPerParticleParameter("sigmaA");
+        force->addPerParticleParameter("sigmaB");
+        force->addPerParticleParameter("eps");
+        if (mode == 1) {
+            // Test with a cutoff.
+            
+            force->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
+            force->setCutoffDistance(cutoff);
+            force->setUseLongRangeCorrection(true);
+            force->setUseSwitchingFunction(true);
+            force->setSwitchingDistance(0.8*cutoff);
+        }
+        if (mode == 2) {
+            // Test with interaction groups.
+            
+            force->addInteractionGroup({0, 1, 2, 3, 4, 5}, {0, 3, 10, 15, 20, 25, 30});
+        }
+        force->setForceGroup(model);
+        for (int i = 0; i < numParticles; i++) {
+            if (i%2 == 0)
+                force->addParticle({1.1, 1.6, 0.5});
+            else
+                force->addParticle({1.0, 1.0, 1.0});
+        }
+        system.addForce(force);
+    }
+    vector<Vec3> positions(numParticles);
+    int index = 0;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                system.addParticle(1.0);
+                positions[index++] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
+            }
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+
+    // Compute the force and energy for a few values of lambda and make sure both versions agree.
+
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    for (double lambda : {0.0, 0.3, 0.7, 1.0}) {
+        context.setParameter("lambda", lambda);
+        double e1 = context.getState(State::Energy, false, 1<<0).getPotentialEnergy();
+        double e2 = context.getState(State::Energy, false, 1<<1).getPotentialEnergy();
+        ASSERT_EQUAL_TOL(e1, e2, 1e-5);
+    }
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -1479,6 +1542,9 @@ int main(int argc, char* argv[]) {
         testEnergyParameterDerivatives();
         testEnergyParameterDerivatives2();
         testEnergyParameterDerivativesWithGroups();
+        testComputedValues(0); // No cutoff
+        testComputedValues(1); // Cutoff, periodic
+        testComputedValues(2); // Interaction groups
         runPlatformTests();
     }
     catch(const exception& e) {

wrappers/python/openmm/app/forcefield.py

@@ -2925,6 +2925,7 @@ class CustomNonbondedGenerator(object):
         self.bondCutoff = bondCutoff
         self.globalParams = {}
         self.perParticleParams = []
+        self.computedValues = {}
         self.functions = []
 
     @staticmethod
@@ -2935,6 +2936,8 @@ class CustomNonbondedGenerator(object):
             generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
         for param in element.findall('PerParticleParameter'):
             generator.perParticleParams.append(param.attrib['name'])
+        for value in element.findall('ComputedValue'):
+            generator.computedValues[value.attrib['name']] = value.attrib['expression']
         generator.params = ForceField._AtomTypeParameters(ff, 'CustomNonbondedForce', 'Atom', generator.perParticleParams)
         generator.params.parseDefinitions(element)
         generator.functions += _parseFunctions(element)
@@ -2953,6 +2956,8 @@ class CustomNonbondedGenerator(object):
             force.addGlobalParameter(param, self.globalParams[param])
         for param in self.perParticleParams:
             force.addPerParticleParameter(param)
+        for name in self.computedValues:
+            force.addComputedValue(name, self.computedValues[name])
         _createFunctions(force, self.functions)
         for atom in data.atoms:
             values = self.params.getAtomParameters(atom, data)