Merge pull request #2119 from peastman/offsetparams

Global parameters for NonbondedForce are explicitly defined

openmmapi/include/openmm/NonbondedForce.h

@@ -80,8 +80,9 @@ namespace OpenMM {
  * In some applications, it is useful to be able to inexpensively change the parameters of small groups of particles.
  * Usually this is done to interpolate between two sets of parameters.  For example, a titratable group might have
  * two states it can exist in, each described by a different set of parameters for the atoms that make up the
- * group.  You might then want to smoothly interpolate between the two states.  This is done by calling
- * addParticleParameterOffset() to create a "parameter offset".  Each offset defines the following:
+ * group.  You might then want to smoothly interpolate between the two states.  This is done by first calling
+ * addGlobalParameter() to define a Context parameter, then addParticleParameterOffset() to create a "parameter offset"
+ * that depends on the Context parameter.  Each offset defines the following:
  * 
  * <ul>
  * <li>A Context parameter used to interpolate between the states.</li>
@@ -158,6 +159,12 @@ public:
     int getNumExceptions() const {
         return exceptions.size();
     }
+    /**
+     * Get the number of global parameters that have been added.
+     */
+    int getNumGlobalParameters() const {
+        return globalParameters.size();
+    }
     /**
      * Get the number of particles parameter offsets that have been added.
      */
@@ -393,11 +400,49 @@ public:
      *                        multiplied by this factor
      */
     void createExceptionsFromBonds(const std::vector<std::pair<int, int> >& bonds, double coulomb14Scale, double lj14Scale);
+    /**
+     * Add a new global parameter that parameter offsets may depend on.  The default value provided to
+     * this method is the initial value of the parameter in newly created Contexts.  You can change
+     * the value at any time by calling setParameter() on the Context.
+     * 
+     * @param name             the name of the parameter
+     * @param defaultValue     the default value of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addGlobalParameter(const std::string& name, double defaultValue);
+    /**
+     * Get the name of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getGlobalParameterName(int index) const;
+    /**
+     * Set the name of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setGlobalParameterName(int index, const std::string& name);
+    /**
+     * Get the default value of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the default value
+     * @return the parameter default value
+     */
+    double getGlobalParameterDefaultValue(int index) const;
+    /**
+     * Set the default value of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the default value
+     * @param defaultValue   the default value of the parameter
+     */
+    void setGlobalParameterDefaultValue(int index, double defaultValue);
     /**
      * Add an offset to the per-particle parameters of a particular particle, based on a global parameter.
      * 
-     * @param parameter       the name of the global parameter.  Its value is initially set to 0, and can be modified
-     *                        at any time by calling Context::setParameter().
+     * @param parameter       the name of the global parameter.  It must have already been added with addGlobalParameter().
+     *                        Its value can be modified at any time by calling Context::setParameter().
      * @param particleIndex   the index of the particle whose parameters are affected
      * @param chargeScale     this value multiplied by the parameter value is added to the particle's charge
      * @param sigmaScale      this value multiplied by the parameter value is added to the particle's sigma
@@ -409,8 +454,7 @@ public:
      * Get the offset added to the per-particle parameters of a particular particle, based on a global parameter.
      * 
      * @param index           the index of the offset to query, as returned by addParticleParameterOffset()
-     * @param parameter       the name of the global parameter.  Its value is initially set to 0, and can be modified
-     *                        at any time by calling Context::setParameter().
+     * @param parameter       the name of the global parameter
      * @param particleIndex   the index of the particle whose parameters are affected
      * @param chargeScale     this value multiplied by the parameter value is added to the particle's charge
      * @param sigmaScale      this value multiplied by the parameter value is added to the particle's sigma
@@ -421,8 +465,8 @@ public:
      * Set the offset added to the per-particle parameters of a particular particle, based on a global parameter.
      * 
      * @param index           the index of the offset to modify, as returned by addParticleParameterOffset()
-     * @param parameter       the name of the global parameter.  Its value is initially set to 0, and can be modified
-     *                        at any time by calling Context::setParameter().
+     * @param parameter       the name of the global parameter.  It must have already been added with addGlobalParameter().
+     *                        Its value can be modified at any time by calling Context::setParameter().
      * @param particleIndex   the index of the particle whose parameters are affected
      * @param chargeScale     this value multiplied by the parameter value is added to the particle's charge
      * @param sigmaScale      this value multiplied by the parameter value is added to the particle's sigma
@@ -432,8 +476,8 @@ public:
     /**
      * Add an offset to the parameters of a particular exception, based on a global parameter.
      * 
-     * @param parameter       the name of the global parameter.  Its value is initially set to 0, and can be modified
-     *                        at any time by calling Context::setParameter().
+     * @param parameter       the name of the global parameter.  It must have already been added with addGlobalParameter().
+     *                        Its value can be modified at any time by calling Context::setParameter().
      * @param exceptionIndex  the index of the exception whose parameters are affected
      * @param chargeProdScale this value multiplied by the parameter value is added to the exception's charge product
      * @param sigmaScale      this value multiplied by the parameter value is added to the exception's sigma
@@ -445,8 +489,7 @@ public:
      * Get the offset added to the parameters of a particular exception, based on a global parameter.
      * 
      * @param index           the index of the offset to query, as returned by addExceptionParameterOffset()
-     * @param parameter       the name of the global parameter.  Its value is initially set to 0, and can be modified
-     *                        at any time by calling Context::setParameter().
+     * @param parameter       the name of the global parameter
      * @param exceptionIndex  the index of the exception whose parameters are affected
      * @param chargeProdScale this value multiplied by the parameter value is added to the exception's charge product
      * @param sigmaScale      this value multiplied by the parameter value is added to the exception's sigma
@@ -457,8 +500,8 @@ public:
      * Set the offset added to the parameters of a particular exception, based on a global parameter.
      * 
      * @param index           the index of the offset to modify, as returned by addExceptionParameterOffset()
-     * @param parameter       the name of the global parameter.  Its value is initially set to 0, and can be modified
-     *                        at any time by calling Context::setParameter().
+     * @param parameter       the name of the global parameter.  It must have already been added with addGlobalParameter().
+     *                        Its value can be modified at any time by calling Context::setParameter().
      * @param exceptionIndex  the index of the exception whose parameters are affected
      * @param chargeProdScale this value multiplied by the parameter value is added to the exception's charge product
      * @param sigmaScale      this value multiplied by the parameter value is added to the exception's sigma
@@ -530,6 +573,7 @@ protected:
 private:
     class ParticleInfo;
     class ExceptionInfo;
+    class GlobalParameterInfo;
     class ParticleOffsetInfo;
     class ExceptionOffsetInfo;
     NonbondedMethod nonbondedMethod;
@@ -537,8 +581,10 @@ private:
     bool useSwitchingFunction, useDispersionCorrection;
     int recipForceGroup, nx, ny, nz, dnx, dny, dnz;
     void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
+    int getGlobalParameterIndex(const std::string& parameter) const;
     std::vector<ParticleInfo> particles;
     std::vector<ExceptionInfo> exceptions;
+    std::vector<GlobalParameterInfo> globalParameters;
     std::vector<ParticleOffsetInfo> particleOffsets;
     std::vector<ExceptionOffsetInfo> exceptionOffsets;
     std::map<std::pair<int, int>, int> exceptionMap;
@@ -576,20 +622,33 @@ public:
     }
 };
 
+/**
+ * This is an internal class used to record information about a global parameter.
+ * @private
+ */
+class NonbondedForce::GlobalParameterInfo {
+public:
+    std::string name;
+    double defaultValue;
+    GlobalParameterInfo() {
+    }
+    GlobalParameterInfo(const std::string& name, double defaultValue) : name(name), defaultValue(defaultValue) {
+    }
+};
+
 /**
  * This is an internal class used to record information about a particle parameter offset.
  * @private
  */
 class NonbondedForce::ParticleOffsetInfo {
 public:
-    std::string parameter;
-    int particle;
+    int particle, parameter;
     double chargeScale, sigmaScale, epsilonScale;
     ParticleOffsetInfo() {
-        particle = -1;
+        particle = parameter = -1;
         chargeScale = sigmaScale = epsilonScale = 0.0;
     }
-    ParticleOffsetInfo(const std::string& parameter, int particle, double chargeScale, double sigmaScale, double epsilonScale) :
+    ParticleOffsetInfo(int parameter, int particle, double chargeScale, double sigmaScale, double epsilonScale) :
         parameter(parameter), particle(particle), chargeScale(chargeScale), sigmaScale(sigmaScale), epsilonScale(epsilonScale) {
     }
 };
@@ -600,14 +659,13 @@ public:
  */
 class NonbondedForce::ExceptionOffsetInfo {
 public:
-    std::string parameter;
-    int exception;
+    int exception, parameter;
     double chargeProdScale, sigmaScale, epsilonScale;
     ExceptionOffsetInfo() {
-        exception = -1;
+        exception = parameter = -1;
         chargeProdScale = sigmaScale = epsilonScale = 0.0;
     }
-    ExceptionOffsetInfo(const std::string& parameter, int exception, double chargeProdScale, double sigmaScale, double epsilonScale) :
+    ExceptionOffsetInfo(int parameter, int exception, double chargeProdScale, double sigmaScale, double epsilonScale) :
         parameter(parameter), exception(exception), chargeProdScale(chargeProdScale), sigmaScale(sigmaScale), epsilonScale(epsilonScale) {
     }
 };

openmmapi/src/NonbondedForce.cpp

@@ -256,14 +256,46 @@ void NonbondedForce::addExclusionsToSet(const vector<set<int> >& bonded12, set<i
     }
 }
 
+int NonbondedForce::addGlobalParameter(const string& name, double defaultValue) {
+    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
+    return globalParameters.size()-1;
+}
+
+const string& NonbondedForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    return globalParameters[index].name;
+}
+
+void NonbondedForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    globalParameters[index].name = name;
+}
+
+double NonbondedForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    return globalParameters[index].defaultValue;
+}
+
+void NonbondedForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    globalParameters[index].defaultValue = defaultValue;
+}
+
+int NonbondedForce::getGlobalParameterIndex(const std::string& parameter) const {
+    for (int i = 0; i < globalParameters.size(); i++)
+        if (globalParameters[i].name == parameter)
+            return i;
+    throw OpenMMException("There is no global parameter called '"+parameter+"'");
+}
+
 int NonbondedForce::addParticleParameterOffset(const std::string& parameter, int particleIndex, double chargeScale, double sigmaScale, double epsilonScale) {
-    particleOffsets.push_back(ParticleOffsetInfo(parameter, particleIndex, chargeScale, sigmaScale, epsilonScale));
+    particleOffsets.push_back(ParticleOffsetInfo(getGlobalParameterIndex(parameter), particleIndex, chargeScale, sigmaScale, epsilonScale));
     return particleOffsets.size()-1;
 }
 
 void NonbondedForce::getParticleParameterOffset(int index, std::string& parameter, int& particleIndex, double& chargeScale, double& sigmaScale, double& epsilonScale) const {
     ASSERT_VALID_INDEX(index, particleOffsets);
-    parameter = particleOffsets[index].parameter;
+    parameter = globalParameters[particleOffsets[index].parameter].name;
     particleIndex = particleOffsets[index].particle;
     chargeScale = particleOffsets[index].chargeScale;
     sigmaScale = particleOffsets[index].sigmaScale;
@@ -272,7 +304,7 @@ void NonbondedForce::getParticleParameterOffset(int index, std::string& paramete
 
 void NonbondedForce::setParticleParameterOffset(int index, const std::string& parameter, int particleIndex, double chargeScale, double sigmaScale, double epsilonScale) {
     ASSERT_VALID_INDEX(index, particleOffsets);
-    particleOffsets[index].parameter = parameter;
+    particleOffsets[index].parameter = getGlobalParameterIndex(parameter);
     particleOffsets[index].particle = particleIndex;
     particleOffsets[index].chargeScale = chargeScale;
     particleOffsets[index].sigmaScale = sigmaScale;
@@ -280,13 +312,13 @@ void NonbondedForce::setParticleParameterOffset(int index, const std::string& pa
 }
 
 int NonbondedForce::addExceptionParameterOffset(const std::string& parameter, int exceptionIndex, double chargeProdScale, double sigmaScale, double epsilonScale) {
-    exceptionOffsets.push_back(ExceptionOffsetInfo(parameter, exceptionIndex, chargeProdScale, sigmaScale, epsilonScale));
+    exceptionOffsets.push_back(ExceptionOffsetInfo(getGlobalParameterIndex(parameter), exceptionIndex, chargeProdScale, sigmaScale, epsilonScale));
     return exceptionOffsets.size()-1;
 }
 
 void NonbondedForce::getExceptionParameterOffset(int index, std::string& parameter, int& exceptionIndex, double& chargeProdScale, double& sigmaScale, double& epsilonScale) const {
     ASSERT_VALID_INDEX(index, exceptionOffsets);
-    parameter = exceptionOffsets[index].parameter;
+    parameter = globalParameters[exceptionOffsets[index].parameter].name;
     exceptionIndex = exceptionOffsets[index].exception;
     chargeProdScale = exceptionOffsets[index].chargeProdScale;
     sigmaScale = exceptionOffsets[index].sigmaScale;
@@ -295,7 +327,7 @@ void NonbondedForce::getExceptionParameterOffset(int index, std::string& paramet
 
 void NonbondedForce::setExceptionParameterOffset(int index, const std::string& parameter, int exceptionIndex, double chargeProdScale, double sigmaScale, double epsilonScale) {
     ASSERT_VALID_INDEX(index, exceptionOffsets);
-    exceptionOffsets[index].parameter = parameter;
+    exceptionOffsets[index].parameter = getGlobalParameterIndex(parameter);
     exceptionOffsets[index].exception = exceptionIndex;
     exceptionOffsets[index].chargeProdScale = chargeProdScale;
     exceptionOffsets[index].sigmaScale = sigmaScale;

openmmapi/src/NonbondedForceImpl.cpp

@@ -111,24 +111,10 @@ double NonbondedForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includ
 }
 
 map<string, double> NonbondedForceImpl::getDefaultParameters() {
-    map<string, double> params;
-    for (int i = 0; i < owner.getNumParticleParameterOffsets(); i++) {
-        string parameter;
-        int particle;
-        double charge, sigma, epsilon;
-        owner.getParticleParameterOffset(i, parameter, particle, charge, sigma, epsilon);
-        if (params.find(parameter) == params.end())
-            params[parameter] = 0.0;
-    }
-    for (int i = 0; i < owner.getNumExceptionParameterOffsets(); i++) {
-        string parameter;
-        int exception;
-        double chargeProd, sigma, epsilon;
-        owner.getExceptionParameterOffset(i, parameter, exception, chargeProd, sigma, epsilon);
-        if (params.find(parameter) == params.end())
-            params[parameter] = 0.0;
-    }
-    return params;
+    map<string, double> parameters;
+    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
+        parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
+    return parameters;
 }
 
 std::vector<std::string> NonbondedForceImpl::getKernelNames() {

serialization/src/NonbondedForceProxy.cpp

@@ -65,6 +65,9 @@ void NonbondedForceProxy::serialize(const void* object, SerializationNode& node)
     node.setIntProperty("ljny", ny);
     node.setIntProperty("ljnz", nz);
     node.setIntProperty("recipForceGroup", force.getReciprocalSpaceForceGroup());
+    SerializationNode& globalParams = node.createChildNode("GlobalParameters");
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParams.createChildNode("Parameter").setStringProperty("name", force.getGlobalParameterName(i)).setDoubleProperty("default", force.getGlobalParameterDefaultValue(i));
     SerializationNode& particleOffsets = node.createChildNode("ParticleOffsets");
     for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
         int particle;
@@ -124,6 +127,9 @@ void* NonbondedForceProxy::deserialize(const SerializationNode& node) const {
         }
         force->setReciprocalSpaceForceGroup(node.getIntProperty("recipForceGroup", -1));
         if (version >= 3) {
+            const SerializationNode& globalParams = node.getChildNode("GlobalParameters");
+            for (auto& parameter : globalParams.getChildren())
+                force->addGlobalParameter(parameter.getStringProperty("name"), parameter.getDoubleProperty("default"));
             const SerializationNode& particleOffsets = node.getChildNode("ParticleOffsets");
             for (auto& offset : particleOffsets.getChildren())
                 force->addParticleParameterOffset(offset.getStringProperty("parameter"), offset.getIntProperty("particle"), offset.getDoubleProperty("q"), offset.getDoubleProperty("sig"), offset.getDoubleProperty("eps"));

serialization/tests/TestSerializeNonbondedForce.cpp

@@ -61,6 +61,8 @@ void testSerialization() {
     force.addParticle(-0.5, 0.3, 0.03);
     force.addException(0, 1, 2, 0.5, 0.1);
     force.addException(1, 2, 0.2, 0.4, 0.2);
+    force.addGlobalParameter("scale1", 1.0);
+    force.addGlobalParameter("scale2", 2.0);
     force.addParticleParameterOffset("scale1", 2, 1.5, 2.0, 2.5);
     force.addExceptionParameterOffset("scale2", 1, -0.1, -0.2, -0.3);
 
@@ -83,6 +85,7 @@ void testSerialization() {
     ASSERT_EQUAL(force.getUseDispersionCorrection(), force2.getUseDispersionCorrection());
     ASSERT_EQUAL(force.getNumParticles(), force2.getNumParticles());
     ASSERT_EQUAL(force.getNumExceptions(), force2.getNumExceptions());
+    ASSERT_EQUAL(force.getNumGlobalParameters(), force2.getNumGlobalParameters());
     ASSERT_EQUAL(force.getNumParticleParameterOffsets(), force2.getNumParticleParameterOffsets());
     ASSERT_EQUAL(force.getNumExceptionParameterOffsets(), force2.getNumExceptionParameterOffsets());
     double alpha2;
@@ -98,7 +101,11 @@ void testSerialization() {
     ASSERT_EQUAL(dalpha, dalpha2);
     ASSERT_EQUAL(dnx, dnx2);
     ASSERT_EQUAL(dny, dny2);
-    ASSERT_EQUAL(dnz, dnz2);    
+    ASSERT_EQUAL(dnz, dnz2);
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        ASSERT_EQUAL(force.getGlobalParameterName(i), force2.getGlobalParameterName(i));
+        ASSERT_EQUAL(force.getGlobalParameterDefaultValue(i), force2.getGlobalParameterDefaultValue(i));
+    }
     for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
         int index1, index2;
         string param1, param2;

tests/TestNonbondedForce.h

@@ -745,6 +745,8 @@ void testParameterOffsets() {
     force->addParticle(0.5, 2.0, 0.8);
     force->addException(0, 1, 1.0, 1.5, 1.0);
     force->addException(2, 3, 0.5, 1.0, 1.5);
+    force->addGlobalParameter("p1", 0.0);
+    force->addGlobalParameter("p2", 1.0);
     force->addParticleParameterOffset("p1", 0, 3.0, 0.5, 0.5);
     force->addParticleParameterOffset("p2", 1, 1.0, 1.0, 2.0);
     force->addExceptionParameterOffset("p1", 1, 0.5, 0.5, 1.5);
@@ -757,7 +759,7 @@ void testParameterOffsets() {
     context.setPositions(positions);
     ASSERT_EQUAL(2, context.getParameters().size());
     ASSERT_EQUAL(0.0, context.getParameter("p1"));
-    ASSERT_EQUAL(0.0, context.getParameter("p2"));
+    ASSERT_EQUAL(1.0, context.getParameter("p2"));
     context.setParameter("p1", 0.5);
     context.setParameter("p2", 1.5);