Merge pull request #1837 from peastman/cv

Created CustomCVForce

docs-source/usersguide/theory.rst

@@ -1158,6 +1158,21 @@ specified in three ways:
 * Per-donor parameters are defined by specifying a value for each donor group.
 * Per-acceptor parameters are defined by specifying a value for each acceptor group.
 
+CustomCVForce
+*************
+
+CustomCVForce computes an energy as a function of "collective variables".  A
+collective variable may be any scalar valued function of the particle positions
+and other parameters.  Each one is defined by a :code:`Force` object, so any
+function that can be defined via any force class (either standard or custom) can
+be used as a collective variable.  The energy is then computed as
+
+.. math::
+   E=f(...)
+
+where *f*\ (...) is a user supplied mathematical expression of the collective
+variables.  It also may depend on user defined global parameters.
+
 
 .. _writing-custom-expressions:
 

olla/include/openmm/Platform.h

@@ -127,6 +127,14 @@ public:
      * @param properties a set of values for platform-specific properties.  Keys are the property names.
      */
     virtual void contextCreated(ContextImpl& context, const std::map<std::string, std::string>& properties) const;
+    /**
+     * This is called whenever a new Context is created using ContextImpl::createLinkedContext().  It gives the
+     * Platform a chance to initialize the context and store platform-specific data in it.
+     *
+     * @param context          the newly created context
+     * @param originalContext  the original context it is linked to
+     */
+    virtual void linkedContextCreated(ContextImpl& context, ContextImpl& originalContext) const;
     /**
      * This is called whenever a Context is deleted.  It gives the Platform a chance to clean up
      * any platform-specific data that was stored in it.

olla/include/openmm/kernels.h

@@ -39,6 +39,7 @@
 #include "openmm/CustomAngleForce.h"
 #include "openmm/CustomBondForce.h"
 #include "openmm/CustomCentroidBondForce.h"
+#include "openmm/CustomCVForce.h"
 #include "openmm/CustomCompoundBondForce.h"
 #include "openmm/CustomExternalForce.h"
 #include "openmm/CustomGBForce.h"
@@ -943,6 +944,43 @@ public:
     virtual void copyParametersToContext(ContextImpl& context, const GayBerneForce& force) = 0;
 };
 
+/**
+ * This kernel is invoked by CustomCVForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCustomCVForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcCustomCVForce";
+    }
+    CalcCustomCVForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomCVForce this kernel will be used for
+     * @param innerContext   the context created by the CustomCVForce for computing collective variables
+     */
+    virtual void initialize(const System& system, const CustomCVForce& force, ContextImpl& innerContext) = 0;
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param innerContext   the context created by the CustomCVForce for computing collective variables
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    virtual double execute(ContextImpl& context, ContextImpl& innerContext, bool includeForces, bool includeEnergy) = 0;
+    /**
+     * Copy state information to the inner context.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param innerContext   the context created by the CustomCVForce for computing collective variables
+     */
+    virtual void copyState(ContextImpl& context, ContextImpl& innerContext) = 0;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

olla/src/Platform.cpp

@@ -34,6 +34,7 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/Kernel.h"
 #include "openmm/KernelFactory.h"
+#include "openmm/internal/ContextImpl.h"
 #ifdef WIN32
 #include <windows.h>
 #include <sstream>
@@ -113,6 +114,16 @@ void Platform::setPropertyDefaultValue(const string& property, const string& val
 void Platform::contextCreated(ContextImpl& context, const map<string, string>& properties) const {
 }
 
+void Platform::linkedContextCreated(ContextImpl& context, ContextImpl& originalContext) const {
+    // The default implementation just copies over the properties and calls contextCreated().
+    // Subclasses may override this to do something different.
+    
+    map<string, string> properties;
+    for (auto& name : getPropertyNames())
+        properties[name] = getPropertyValue(originalContext.getOwner(), name);
+    contextCreated(context, properties);
+}
+
 void Platform::contextDestroyed(ContextImpl& context) const {
 }
 

openmmapi/include/OpenMM.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) 2009-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -43,6 +43,7 @@
 #include "openmm/CustomAngleForce.h"
 #include "openmm/CustomTorsionForce.h"
 #include "openmm/CustomExternalForce.h"
+#include "openmm/CustomCVForce.h"
 #include "openmm/CustomGBForce.h"
 #include "openmm/CustomHbondForce.h"
 #include "openmm/CustomIntegrator.h"

openmmapi/include/openmm/Context.h

@@ -263,8 +263,11 @@ public:
      */
     const std::vector<std::vector<int> >& getMolecules() const;
 private:
+    friend class ContextImpl;
     friend class Force;
+    friend class ForceImpl;
     friend class Platform;
+    Context(const System& system, Integrator& integrator, ContextImpl& linked);
     ContextImpl& getImpl();
     const ContextImpl& getImpl() const;
     ContextImpl* impl;

openmmapi/include/openmm/CustomCVForce.h

+#ifndef OPENMM_CUSTOMCVFORCE_H_
+#define OPENMM_CUSTOMCVFORCE_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2017 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "Force.h"
+#include "TabulatedFunction.h"
+#include "internal/windowsExport.h"
+#include <string>
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This class supports energy functions that depend on collective variables.  To use it,
+ * you define a set of collective variables (scalar valued functions that depend on the
+ * particle positions), and an algebraic expression for the energy as a function of the
+ * collective variables.  The expression also may involve tabulated functions, and may
+ * depend on arbitrary global parameters.
+ * 
+ * Each collective variable is defined by a Force object.  The Force's potential energy
+ * is computed, and that becomes the value of the variable.  This provides enormous
+ * flexibility in defining collective variables, especially by using custom forces.
+ * Anything that can be computed as a potential function can also be used as a collective
+ * variable.
+ *
+ * To use this class, create a CustomCVForce object, passing an algebraic expression to the
+ * constructor that defines the potential energy.  Then call addCollectiveVariable() to define
+ * collective variables and addGlobalParameter() to define global parameters.  The values
+ * of global parameters may be modified during a simulation by calling Context::setParameter().
+ *
+ * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.
+ * Call addEnergyParameterDerivative() to request that the derivative with respect to a particular parameter be
+ * computed.  You can then query its value in a Context by calling getState() on it.
+ *
+ * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
+ * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.
+ *
+ * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by
+ * creating a TabulatedFunction object.  That function can then appear in the expression.
+ */
+
+class OPENMM_EXPORT CustomCVForce : public Force {
+public:
+    /**
+     * Create a CustomCVForce.
+     *
+     * @param energy   an algebraic expression giving the energy of the system as a function
+     *                 of the collective variables and global parameters
+     */
+    explicit CustomCVForce(const std::string& energy);
+    ~CustomCVForce();
+    /**
+     * Get the number of collective variables that the interaction depends on.
+     */
+    int getNumCollectiveVariables() const {
+        return variables.size();
+    }
+    /**
+     * Get the number of global parameters that the interaction depends on.
+     */
+    int getNumGlobalParameters() const {
+        return globalParameters.size();
+    }
+    /**
+     * Get the number of global parameters with respect to which the derivative of the energy
+     * should be computed.
+     */
+    int getNumEnergyParameterDerivatives() const {
+        return energyParameterDerivatives.size();
+    }
+    /**
+     * Get the number of tabulated functions that have been defined.
+     */
+    int getNumTabulatedFunctions() const {
+        return functions.size();
+    }
+    /**
+     * Get the algebraic expression that gives the energy of the system
+     */
+    const std::string& getEnergyFunction() const;
+    /**
+     * Set the algebraic expression that gives the energy of the system
+     */
+    void setEnergyFunction(const std::string& energy);
+    /**
+     * Add a collective variable that the force may depend on.  The collective variable
+     * is represented by a Force object, which should have been created on the heap with the
+     * "new" operator.  The CustomCVForce takes over ownership of it, and deletes the Force when the
+     * CustomCVForce itself is deleted.
+     * 
+     * @param name      the name of the collective variable, as it will appear in the energy expression
+     * @param variable  the collective variable, represented by a Force object.  The value of the
+     *                  variable is the energy computed by the Force.
+     * @return the index within the Force of the variable that was added
+     */
+    int addCollectiveVariable(const std::string& name, Force* variable);
+    /**
+     * Get the name of a collective variable.
+     *
+     * @param index     the index of the collective variable for which to get the name
+     * @return the variable name
+     */
+    const std::string& getCollectiveVariableName(int index) const;
+    /**
+     * Get a writable reference to the Force object that computes a collective variable.
+     *
+     * @param index     the index of the collective variable to get
+     * @return the Force object
+     */
+    Force& getCollectiveVariable(int index);
+    /**
+     * Get a const reference to the Force object that computes a collective variable.
+     *
+     * @param index     the index of the collective variable to get
+     * @return the Force object
+     */
+    const Force& getCollectiveVariable(int index) const;
+    /**
+     * Add a new global parameter that the interaction may depend on.
+     *
+     * @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);
+    /**
+     * Request that this Force compute the derivative of its energy with respect to a global parameter.
+     * The parameter must have already been added with addGlobalParameter().
+     *
+     * @param name             the name of the parameter
+     */
+    void addEnergyParameterDerivative(const std::string& name);
+    /**
+     * Get the name of a global parameter with respect to which this Force should compute the
+     * derivative of the energy.
+     *
+     * @param index     the index of the parameter derivative, between 0 and getNumEnergyParameterDerivatives()
+     * @return the parameter name
+     */
+    const std::string& getEnergyParameterDerivativeName(int index) const;
+    /**
+     * Add a tabulated function that may appear in the energy expression.
+     *
+     * @param name           the name of the function as it appears in expressions
+     * @param function       a TabulatedFunction object defining the function.  The TabulatedFunction
+     *                       should have been created on the heap with the "new" operator.  The
+     *                       Force takes over ownership of it, and deletes it when the Force itself is deleted.
+     * @return the index of the function that was added
+     */
+    int addTabulatedFunction(const std::string& name, TabulatedFunction* function);
+    /**
+     * Get a const reference to a tabulated function that may appear in the energy expression.
+     *
+     * @param index     the index of the function to get
+     * @return the TabulatedFunction object defining the function
+     */
+    const TabulatedFunction& getTabulatedFunction(int index) const;
+    /**
+     * Get a reference to a tabulated function that may appear in the energy expression.
+     *
+     * @param index     the index of the function to get
+     * @return the TabulatedFunction object defining the function
+     */
+    TabulatedFunction& getTabulatedFunction(int index);
+    /**
+     * Get the name of a tabulated function that may appear in the energy expression.
+     *
+     * @param index     the index of the function to get
+     * @return the name of the function as it appears in expressions
+     */
+    const std::string& getTabulatedFunctionName(int index) const;
+    /**
+     * Get the current values of the collective variables in a Context.
+     *
+     * @param context        the Context for which to get the values
+     * @param[out] values    the values of the collective variables are computed and
+     *                       stored into this
+     */
+    void getCollectiveVariableValues(Context& context, std::vector<double>& values);
+    /**
+     * Get the inner Context used for evaluating collective variables.
+     * 
+     * When you create a Context for a System that contains a CustomCVForce, internally
+     * it creates a new System, adds the Forces that define the CVs to it, creates a new
+     * Context for that System, and uses it to evaluate the variables.  In most cases you
+     * can ignore all of this.  It is just an implementation detail.  However, there are
+     * a few cases where you need to directly access that internal Context.  For example,
+     * if you want to modify one of the Forces that defines a collective variable and
+     * call updateParametersInContext() on it, you need to pass that inner Context to it.
+     * This method returns a reference to it.
+     * 
+     * @param context    the Context containing the CustomCVForce
+     * @return the inner Context used to evaluate the collective variables
+     */
+    Context& getInnerContext(Context& context);
+    /**
+     * Returns whether or not this force makes use of periodic boundary
+     * conditions.
+     *
+     * @returns true if force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions() const;
+protected:
+    ForceImpl* createImpl() const;
+private:
+    class GlobalParameterInfo;
+    class VariableInfo;
+    class FunctionInfo;
+    std::string energyExpression;
+    std::vector<GlobalParameterInfo> globalParameters;
+    std::vector<VariableInfo> variables;
+    std::vector<FunctionInfo> functions;
+    std::vector<int> energyParameterDerivatives;
+};
+
+/**
+ * This is an internal class used to record information about a global parameter.
+ * @private
+ */
+class CustomCVForce::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 tabulated function.
+ * @private
+ */
+class CustomCVForce::VariableInfo {
+public:
+    std::string name;
+    Force* variable;
+    VariableInfo() {
+    }
+    VariableInfo(const std::string& name, Force* variable) : name(name), variable(variable) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a tabulated function.
+ * @private
+ */
+class CustomCVForce::FunctionInfo {
+public:
+    std::string name;
+    TabulatedFunction* function;
+    FunctionInfo() {
+    }
+    FunctionInfo(const std::string& name, TabulatedFunction* function) : name(name), function(function) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMCVFORCE_H_*/

openmmapi/include/openmm/internal/ContextImpl.h

@@ -55,7 +55,8 @@ public:
     /**
      * Create an ContextImpl for a Context;
      */
-    ContextImpl(Context& owner, const System& system, Integrator& integrator, Platform* platform, const std::map<std::string, std::string>& properties);
+    ContextImpl(Context& owner, const System& system, Integrator& integrator, Platform* platform, const std::map<std::string, std::string>& properties,
+            ContextImpl* originalContext=NULL);
     ~ContextImpl();
     /**
      * Get the Context for which this is the implementation.
@@ -264,8 +265,18 @@ public:
      * you should never call it.  It is exposed here because the same logic is useful to other classes too.
      */
     static std::vector<std::vector<int> > findMolecules(int numParticles, std::vector<std::vector<int> >& particleBonds);
+    /**
+     * Create a new Context based on this one.  The new context will use the same Platform, device, and property
+     * values as this one.  With the CUDA and OpenCL platforms, it also shares the same GPU context, allowing data
+     * to be transferred between them without leaving the GPU.
+     * 
+     * This method exists for very specialized purposes.  If you aren't certain whether you should use it, that probably
+     * means you shouldn't.
+     */
+    Context* createLinkedContext(const System& system, Integrator& integrator);
 private:
     friend class Context;
+    void initialize();
     Context& owner;
     const System& system;
     Integrator& integrator;

openmmapi/include/openmm/internal/CustomCVForceImpl.h

+#ifndef OPENMM_CUSTOMCVFORCEIMPL_H_
+#define OPENMM_CUSTOMCVFORCEIMPL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2017 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "ForceImpl.h"
+#include "openmm/Context.h"
+#include "openmm/CustomCVForce.h"
+#include "openmm/Kernel.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include <map>
+#include <string>
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of CustomCVForce.
+ */
+
+class OPENMM_EXPORT CustomCVForceImpl : public ForceImpl {
+public:
+    CustomCVForceImpl(const CustomCVForce& owner);
+    ~CustomCVForceImpl();
+    void initialize(ContextImpl& context);
+    const CustomCVForce& getOwner() const {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context) {
+        // This force field doesn't update the state directly.
+    }
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
+    std::map<std::string, double> getDefaultParameters();
+    std::vector<std::string> getKernelNames();
+    void getCollectiveVariableValues(ContextImpl& context, std::vector<double>& values);
+    Context& getInnerContext();
+private:
+    const CustomCVForce& owner;
+    Kernel kernel;
+    System innerSystem;
+    VerletIntegrator innerIntegrator;
+    Context* innerContext;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMCVFORCEIMPL_H_*/

openmmapi/include/openmm/internal/ForceImpl.h

@@ -32,6 +32,7 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include "openmm/Context.h"
 #include "openmm/internal/windowsExport.h"
 #include <map>
 #include <string>
@@ -104,6 +105,13 @@ public:
     virtual std::vector<std::pair<int, int> > getBondedParticles() const {
         return std::vector<std::pair<int, int> >(0);
     }
+protected:
+    /**
+     * Get the ContextImpl corresponding to a Context.
+     */
+    ContextImpl& getContextImpl(Context& context) {
+        return context.getImpl();
+    }
 };
 
 } // namespace OpenMM

openmmapi/src/Context.cpp

@@ -40,16 +40,25 @@
 using namespace OpenMM;
 using namespace std;
 
+Context::Context(const System& system, Integrator& integrator, ContextImpl& linked) : properties(linked.getOwner().properties) {
+    // This is used by ContextImpl::createLinkedContext().
+    impl = new ContextImpl(*this, system, integrator, &linked.getPlatform(), properties, &linked);
+    impl->initialize();
+}
+
 Context::Context(const System& system, Integrator& integrator) : properties(map<string, string>()) {
     impl = new ContextImpl(*this, system, integrator, 0, properties);
+    impl->initialize();
 }
 
 Context::Context(const System& system, Integrator& integrator, Platform& platform) : properties(map<string, string>()) {
     impl = new ContextImpl(*this, system, integrator, &platform, properties);
+    impl->initialize();
 }
 
 Context::Context(const System& system, Integrator& integrator, Platform& platform, const map<string, string>& properties) : properties(properties) {
     impl = new ContextImpl(*this, system, integrator, &platform, properties);
+    impl->initialize();
 }
 
 Context::~Context() {
@@ -240,6 +249,7 @@ void Context::reinitialize() {
     integrator.cleanup();
     delete impl;
     impl = new ContextImpl(*this, system, integrator, &platform, properties);
+    impl->initialize();
 }
 
 void Context::createCheckpoint(ostream& stream) {

openmmapi/src/ContextImpl.cpp

@@ -53,7 +53,7 @@ using namespace std;
 const static char CHECKPOINT_MAGIC_BYTES[] = "OpenMM Binary Checkpoint\n";
 
 
-ContextImpl::ContextImpl(Context& owner, const System& system, Integrator& integrator, Platform* platform, const map<string, string>& properties) :
+ContextImpl::ContextImpl(Context& owner, const System& system, Integrator& integrator, Platform* platform, const map<string, string>& properties, ContextImpl* originalContext) :
         owner(owner), system(system), integrator(integrator), hasInitializedForces(false), hasSetPositions(false), integratorIsDeleted(false),
         lastForceGroups(-1), platform(platform), platformData(NULL) {
     int numParticles = system.getNumParticles();
@@ -119,8 +119,6 @@ ContextImpl::ContextImpl(Context& owner, const System& system, Integrator& integ
     kernelNames.push_back(VirtualSitesKernel::Name());
     for (int i = 0; i < system.getNumForces(); ++i) {
         forceImpls.push_back(system.getForce(i).createImpl());
-        map<string, double> forceParameters = forceImpls[forceImpls.size()-1]->getDefaultParameters();
-        parameters.insert(forceParameters.begin(), forceParameters.end());
         vector<string> forceKernels = forceImpls[forceImpls.size()-1]->getKernelNames();
         kernelNames.insert(kernelNames.begin(), forceKernels.begin(), forceKernels.end());
     }
@@ -154,7 +152,10 @@ ContextImpl::ContextImpl(Context& owner, const System& system, Integrator& integ
     for (int i = candidatePlatforms.size()-1; i >= 0; i--) {
         try {
             this->platform = platform = candidatePlatforms[i].second;
-            platform->contextCreated(*this, validatedProperties);
+            if (originalContext == NULL)
+                platform->contextCreated(*this, validatedProperties);
+            else
+                platform->linkedContextCreated(*this, *originalContext);
             break;
         }
         catch (...) {
@@ -163,7 +164,9 @@ ContextImpl::ContextImpl(Context& owner, const System& system, Integrator& integ
             throw;
         }
     }
-    
+}
+
+void ContextImpl::initialize() {
     // Create and initialize kernels and other objects.
     
     initializeForcesKernel = platform->createKernel(CalcForcesAndEnergyKernel::Name(), *this);
@@ -177,10 +180,13 @@ ContextImpl::ContextImpl(Context& owner, const System& system, Integrator& integ
     Vec3 periodicBoxVectors[3];
     system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
     updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
-    for (size_t i = 0; i < forceImpls.size(); ++i)
+    for (size_t i = 0; i < forceImpls.size(); ++i) {
         forceImpls[i]->initialize(*this);
+        map<string, double> forceParameters = forceImpls[i]->getDefaultParameters();
+        parameters.insert(forceParameters.begin(), forceParameters.end());
+    }
     integrator.initialize(*this);
-    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, vector<Vec3>(numParticles));
+    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, vector<Vec3>(system.getNumParticles()));
 }
 
 ContextImpl::~ContextImpl() {
@@ -478,3 +484,7 @@ void ContextImpl::loadCheckpoint(istream& stream) {
 void ContextImpl::systemChanged() {
     integrator.stateChanged(State::Energy);
 }
+
+Context* ContextImpl::createLinkedContext(const System& system, Integrator& integrator) {
+    return new Context(system, integrator, *this);
+}

openmmapi/src/CustomCVForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2017 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/OpenMMException.h"
+#include "openmm/CustomCVForce.h"
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/internal/CustomCVForceImpl.h"
+#include <cmath>
+#include <map>
+#include <set>
+#include <utility>
+
+using namespace OpenMM;
+using namespace std;
+
+CustomCVForce::CustomCVForce(const string& energy) : energyExpression(energy) {
+}
+
+CustomCVForce::~CustomCVForce() {
+    for (auto variable : variables)
+        delete variable.variable;
+    for (auto function : functions)
+        delete function.function;
+}
+
+const string& CustomCVForce::getEnergyFunction() const {
+    return energyExpression;
+}
+
+void CustomCVForce::setEnergyFunction(const std::string& energy) {
+    energyExpression = energy;
+}
+
+int CustomCVForce::addCollectiveVariable(const std::string& name, Force* variable) {
+    if (variables.size() >= 32)
+        throw OpenMMException("CustomCVForce cannot have more than 32 collective variables");
+    variables.push_back(VariableInfo(name, variable));
+    return variables.size()-1;
+}
+
+const string& CustomCVForce::getCollectiveVariableName(int index) const {
+    ASSERT_VALID_INDEX(index, variables);
+    return variables[index].name;
+}
+
+Force& CustomCVForce::getCollectiveVariable(int index) {
+    ASSERT_VALID_INDEX(index, variables);
+    return *variables[index].variable;
+}
+
+const Force& CustomCVForce::getCollectiveVariable(int index) const {
+    ASSERT_VALID_INDEX(index, variables);
+    return *variables[index].variable;
+}
+
+int CustomCVForce::addGlobalParameter(const string& name, double defaultValue) {
+    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
+    return globalParameters.size()-1;
+}
+
+const string& CustomCVForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    return globalParameters[index].name;
+}
+
+void CustomCVForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    globalParameters[index].name = name;
+}
+
+double CustomCVForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    return globalParameters[index].defaultValue;
+}
+
+void CustomCVForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    globalParameters[index].defaultValue = defaultValue;
+}
+
+void CustomCVForce::addEnergyParameterDerivative(const string& name) {
+    for (int i = 0; i < globalParameters.size(); i++)
+        if (name == globalParameters[i].name) {
+            energyParameterDerivatives.push_back(i);
+            return;
+        }
+    throw OpenMMException(string("addEnergyParameterDerivative: Unknown global parameter '"+name+"'"));
+}
+
+const string& CustomCVForce::getEnergyParameterDerivativeName(int index) const {
+    ASSERT_VALID_INDEX(index, energyParameterDerivatives);
+    return globalParameters[energyParameterDerivatives[index]].name;
+}
+
+int CustomCVForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
+    functions.push_back(FunctionInfo(name, function));
+    return functions.size()-1;
+}
+
+const TabulatedFunction& CustomCVForce::getTabulatedFunction(int index) const {
+    ASSERT_VALID_INDEX(index, functions);
+    return *functions[index].function;
+}
+
+TabulatedFunction& CustomCVForce::getTabulatedFunction(int index) {
+    ASSERT_VALID_INDEX(index, functions);
+    return *functions[index].function;
+}
+
+const string& CustomCVForce::getTabulatedFunctionName(int index) const {
+    ASSERT_VALID_INDEX(index, functions);
+    return functions[index].name;
+}
+
+void CustomCVForce::getCollectiveVariableValues(Context& context, vector<double>& values) {
+    dynamic_cast<CustomCVForceImpl&>(getImplInContext(context)).getCollectiveVariableValues(getContextImpl(context), values);
+}
+
+ForceImpl* CustomCVForce::createImpl() const {
+    return new CustomCVForceImpl(*this);
+}
+
+Context& CustomCVForce::getInnerContext(Context& context) {
+    return dynamic_cast<CustomCVForceImpl&>(getImplInContext(context)).getInnerContext();
+}
+
+bool CustomCVForce::usesPeriodicBoundaryConditions() const {
+    for (auto& variable : variables)
+        if (variable.variable->usesPeriodicBoundaryConditions())
+            return true;
+    return false;
+}

openmmapi/src/CustomCVForceImpl.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2017 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/NonbondedForce.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/CustomCVForceImpl.h"
+#include "openmm/kernels.h"
+#include "openmm/serialization/XmlSerializer.h"
+#include <map>
+
+using namespace OpenMM;
+using namespace std;
+
+CustomCVForceImpl::CustomCVForceImpl(const CustomCVForce& owner) : owner(owner), innerIntegrator(1.0),
+        innerContext(NULL) {
+}
+
+CustomCVForceImpl::~CustomCVForceImpl() {
+    if (innerContext != NULL)
+        delete innerContext;
+}
+
+void CustomCVForceImpl::initialize(ContextImpl& context) {
+    // Construct the inner system used to evaluate collective variables.
+    
+    const System& system = context.getSystem();
+    for (int i = 0; i < system.getNumParticles(); i++)
+        innerSystem.addParticle(system.getParticleMass(i));
+    for (int i = 0; i < owner.getNumCollectiveVariables(); i++) {
+        Force* variable = XmlSerializer::clone<Force>(owner.getCollectiveVariable(i));
+        variable->setForceGroup(i);
+        NonbondedForce* nonbonded = dynamic_cast<NonbondedForce*>(variable);
+        if (nonbonded != NULL)
+            nonbonded->setReciprocalSpaceForceGroup(-1);
+        innerSystem.addForce(variable);
+    }
+    
+    // Create the inner context.
+    
+    innerContext = context.createLinkedContext(innerSystem, innerIntegrator);
+    vector<Vec3> positions(system.getNumParticles(), Vec3());
+    innerContext->setPositions(positions);
+    
+    // Create the kernel.
+    
+    kernel = context.getPlatform().createKernel(CalcCustomCVForceKernel::Name(), context);
+    kernel.getAs<CalcCustomCVForceKernel>().initialize(context.getSystem(), owner, getContextImpl(*innerContext));
+}
+
+double CustomCVForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&(1<<owner.getForceGroup())) != 0)
+        return kernel.getAs<CalcCustomCVForceKernel>().execute(context, getContextImpl(*innerContext), includeForces, includeEnergy);
+    return 0.0;
+}
+
+vector<string> CustomCVForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCustomCVForceKernel::Name());
+    return names;
+}
+
+map<string, double> CustomCVForceImpl::getDefaultParameters() {
+    map<string, double> parameters;
+    parameters.insert(innerContext->getParameters().begin(), innerContext->getParameters().end());
+    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
+        parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
+    return parameters;
+}
+
+void CustomCVForceImpl::getCollectiveVariableValues(ContextImpl& context, vector<double>& values) {
+    kernel.getAs<CalcCustomCVForceKernel>().copyState(context, getContextImpl(*innerContext));
+    values.clear();
+    for (int i = 0; i < innerSystem.getNumForces(); i++) {
+        double value = innerContext->getState(State::Energy, false, 1<<i).getPotentialEnergy();
+        values.push_back(value);
+    }
+}
+
+Context& CustomCVForceImpl::getInnerContext() {
+    return *innerContext;
+}

platforms/cuda/include/CudaContext.h

@@ -77,7 +77,8 @@ public:
     static const int ThreadBlockSize;
     static const int TileSize;
     CudaContext(const System& system, int deviceIndex, bool useBlockingSync, const std::string& precision,
-            const std::string& compiler, const std::string& tempDir, const std::string& hostCompiler, CudaPlatform::PlatformData& platformData);
+            const std::string& compiler, const std::string& tempDir, const std::string& hostCompiler, CudaPlatform::PlatformData& platformData,
+            CudaContext* originalContext);
     ~CudaContext();
     /**
      * This is called to initialize internal data structures after all Forces in the system
@@ -284,6 +285,10 @@ public:
      * Clear all buffers that have been registered with addAutoclearBuffer().
      */
     void clearAutoclearBuffers();
+    /**
+     * Sum the buffer containing energy.
+     */
+    double reduceEnergy();
     /**
      * Get the current simulation time.
      */
@@ -622,6 +627,7 @@ private:
     int numAtomBlocks;
     int numThreadBlocks;
     bool useBlockingSync, useDoublePrecision, useMixedPrecision, contextIsValid, atomsWereReordered, boxIsTriclinic, hasCompilerKernel, isNvccAvailable, forcesValid;
+    bool isLinkedContext;
     std::string compiler, tempDir, cacheDir, gpuArchitecture;
     float4 periodicBoxVecXFloat, periodicBoxVecYFloat, periodicBoxVecZFloat, periodicBoxSizeFloat, invPeriodicBoxSizeFloat;
     double4 periodicBoxVecX, periodicBoxVecY, periodicBoxVecZ, periodicBoxSize, invPeriodicBoxSize;
@@ -636,6 +642,7 @@ private:
     CUfunction clearFourBuffersKernel;
     CUfunction clearFiveBuffersKernel;
     CUfunction clearSixBuffersKernel;
+    CUfunction reduceEnergyKernel;
     CUfunction setChargesKernel;
     std::vector<CudaForceInfo*> forces;
     std::vector<Molecule> molecules;
@@ -647,6 +654,7 @@ private:
     CudaArray* velm;
     CudaArray* force;
     CudaArray* energyBuffer;
+    CudaArray* energySum;
     CudaArray* energyParamDerivBuffer;
     CudaArray* atomIndexDevice;
     CudaArray* chargeBuffer;

platforms/cuda/include/CudaKernels.h

@@ -38,6 +38,7 @@
 #include "openmm/internal/CompiledExpressionSet.h"
 #include "openmm/internal/CustomIntegratorUtilities.h"
 #include "lepton/CompiledExpression.h"
+#include "lepton/ExpressionProgram.h"
 #include <cufft.h>
 
 namespace OpenMM {
@@ -1229,6 +1230,54 @@ private:
     CUevent event;
 };
 
+/**
+ * This kernel is invoked by CustomCVForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcCustomCVForceKernel : public CalcCustomCVForceKernel {
+public:
+    CudaCalcCustomCVForceKernel(std::string name, const Platform& platform, CudaContext& cu) : CalcCustomCVForceKernel(name, platform),
+            cu(cu), hasInitializedListeners(false), invAtomOrder(NULL), innerInvAtomOrder(NULL) {
+    }
+    ~CudaCalcCustomCVForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomCVForce this kernel will be used for
+     * @param innerContext   the context created by the CustomCVForce for computing collective variables
+     */
+    void initialize(const System& system, const CustomCVForce& force, ContextImpl& innerContext);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param innerContext   the context created by the CustomCVForce for computing collective variables
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, ContextImpl& innerContext, bool includeForces, bool includeEnergy);
+    /**
+     * Copy state information to the inner context.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param innerContext   the context created by the CustomCVForce for computing collective variables
+     */
+    void copyState(ContextImpl& context, ContextImpl& innerContext);
+private:
+    class ReorderListener;
+    CudaContext& cu;
+    bool hasInitializedListeners;
+    Lepton::ExpressionProgram energyExpression;
+    std::vector<std::string> variableNames, paramDerivNames, globalParameterNames;
+    std::vector<Lepton::ExpressionProgram> variableDerivExpressions;
+    std::vector<Lepton::ExpressionProgram> paramDerivExpressions;
+    std::vector<CudaArray*> cvForces;
+    CudaArray* invAtomOrder;
+    CudaArray* innerInvAtomOrder;
+    CUfunction copyStateKernel, copyForcesKernel, addForcesKernel;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/cuda/include/CudaPlatform.h

@@ -53,6 +53,7 @@ public:
     const std::string& getPropertyValue(const Context& context, const std::string& property) const;
     void setPropertyValue(Context& context, const std::string& property, const std::string& value) const;
     void contextCreated(ContextImpl& context, const std::map<std::string, std::string>& properties) const;
+    void linkedContextCreated(ContextImpl& context, ContextImpl& originalContext) const;
     void contextDestroyed(ContextImpl& context) const;
     /**
      * This is the name of the parameter for selecting which CUDA device or devices to use.
@@ -130,7 +131,7 @@ class OPENMM_EXPORT_CUDA CudaPlatform::PlatformData {
 public:
     PlatformData(ContextImpl* context, const System& system, const std::string& deviceIndexProperty, const std::string& blockingProperty, const std::string& precisionProperty,
             const std::string& cpuPmeProperty, const std::string& compilerProperty, const std::string& tempProperty, const std::string& hostCompilerProperty,
-            const std::string& pmeStreamProperty, const std::string& deterministicForcesProperty, int numThreads);
+            const std::string& pmeStreamProperty, const std::string& deterministicForcesProperty, int numThreads, ContextImpl* originalContext);
     ~PlatformData();
     void initializeContexts(const System& system);
     void syncContexts();

platforms/cuda/src/CudaContext.cpp

@@ -106,9 +106,9 @@ static int executeInWindows(const string &command) {
 #endif
 
 CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlockingSync, const string& precision, const string& compiler,
-        const string& tempDir, const std::string& hostCompiler, CudaPlatform::PlatformData& platformData) : system(system), currentStream(0),
+        const string& tempDir, const std::string& hostCompiler, CudaPlatform::PlatformData& platformData, CudaContext* originalContext) : system(system), currentStream(0),
         time(0.0), platformData(platformData), stepCount(0), computeForceCount(0), stepsSinceReorder(99999), contextIsValid(false), atomsWereReordered(false), hasCompilerKernel(false), isNvccAvailable(false),
-        pinnedBuffer(NULL), posq(NULL), posqCorrection(NULL), velm(NULL), force(NULL), energyBuffer(NULL), energyParamDerivBuffer(NULL), atomIndexDevice(NULL), chargeBuffer(NULL),
+        pinnedBuffer(NULL), posq(NULL), posqCorrection(NULL), velm(NULL), force(NULL), energyBuffer(NULL), energySum(NULL), energyParamDerivBuffer(NULL), atomIndexDevice(NULL), chargeBuffer(NULL),
         integration(NULL), expression(NULL), bonded(NULL), nonbonded(NULL), thread(NULL) {
     // Determine what compiler to use.
     
@@ -173,40 +173,49 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
     cacheDir = cacheDir+"/";
 #endif
     contextIndex = platformData.contexts.size();
-    int numDevices;
     string errorMessage = "Error initializing Context";
-    CHECK_RESULT(cuDeviceGetCount(&numDevices));
-    if (deviceIndex < -1 || deviceIndex >= numDevices)
-        throw OpenMMException("Illegal value for DeviceIndex: "+intToString(deviceIndex));
-
-    vector<int> devicePrecedence;
-    if (deviceIndex == -1) {
-        devicePrecedence = getDevicePrecedence();
-    } else {
-        devicePrecedence.push_back(deviceIndex);
-    }
-
-    this->deviceIndex = -1;
-    for (int i = 0; i < static_cast<int>(devicePrecedence.size()); i++) {
-        int trialDeviceIndex = devicePrecedence[i];
-        CHECK_RESULT(cuDeviceGet(&device, trialDeviceIndex));
-        defaultOptimizationOptions = "--use_fast_math";
-        unsigned int flags = CU_CTX_MAP_HOST;
-        if (useBlockingSync)
-            flags += CU_CTX_SCHED_BLOCKING_SYNC;
-        else
-            flags += CU_CTX_SCHED_SPIN;
+    if (originalContext == NULL) {
+        isLinkedContext = false;
+        int numDevices;
+        CHECK_RESULT(cuDeviceGetCount(&numDevices));
+        if (deviceIndex < -1 || deviceIndex >= numDevices)
+            throw OpenMMException("Illegal value for DeviceIndex: "+intToString(deviceIndex));
+
+        vector<int> devicePrecedence;
+        if (deviceIndex == -1) {
+            devicePrecedence = getDevicePrecedence();
+        } else {
+            devicePrecedence.push_back(deviceIndex);
+        }
 
-        if (cuCtxCreate(&context, flags, device) == CUDA_SUCCESS) {
-            this->deviceIndex = trialDeviceIndex;
-            break;
+        this->deviceIndex = -1;
+        for (int i = 0; i < static_cast<int>(devicePrecedence.size()); i++) {
+            int trialDeviceIndex = devicePrecedence[i];
+            CHECK_RESULT(cuDeviceGet(&device, trialDeviceIndex));
+            defaultOptimizationOptions = "--use_fast_math";
+            unsigned int flags = CU_CTX_MAP_HOST;
+            if (useBlockingSync)
+                flags += CU_CTX_SCHED_BLOCKING_SYNC;
+            else
+                flags += CU_CTX_SCHED_SPIN;
+
+            if (cuCtxCreate(&context, flags, device) == CUDA_SUCCESS) {
+                this->deviceIndex = trialDeviceIndex;
+                break;
+            }
         }
+        if (this->deviceIndex == -1)
+            if (deviceIndex != -1)
+                throw OpenMMException("The requested CUDA device could not be loaded");
+            else
+                throw OpenMMException("No compatible CUDA device is available");
+    }
+    else {
+        isLinkedContext = true;
+        context = originalContext->context;
+        this->deviceIndex = originalContext->deviceIndex;
+        this->device = originalContext->device;
     }
-    if (this->deviceIndex == -1)
-        if (deviceIndex != -1)
-            throw OpenMMException("The requested CUDA device could not be loaded");
-        else
-            throw OpenMMException("No compatible CUDA device is available");
 
     int major, minor;
     CHECK_RESULT(cuDeviceComputeCapability(&major, &minor, device));
@@ -298,6 +307,7 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
     clearFourBuffersKernel = getKernel(utilities, "clearFourBuffers");
     clearFiveBuffersKernel = getKernel(utilities, "clearFiveBuffers");
     clearSixBuffersKernel = getKernel(utilities, "clearSixBuffers");
+    reduceEnergyKernel = getKernel(utilities, "reduceEnergy");
     setChargesKernel = getKernel(utilities, "setCharges");
 
     // Set defines based on the requested precision.
@@ -411,6 +421,8 @@ CudaContext::~CudaContext() {
         delete force;
     if (energyBuffer != NULL)
         delete energyBuffer;
+    if (energySum != NULL)
+        delete energySum;
     if (energyParamDerivBuffer != NULL)
         delete energyParamDerivBuffer;
     if (atomIndexDevice != NULL)
@@ -428,7 +440,7 @@ CudaContext::~CudaContext() {
     if (thread != NULL)
         delete thread;
     string errorMessage = "Error deleting Context";
-    if (contextIsValid) {
+    if (contextIsValid && !isLinkedContext) {
         cuProfilerStop();
         CHECK_RESULT(cuCtxDestroy(context));
     }
@@ -441,16 +453,19 @@ void CudaContext::initialize() {
     int numEnergyBuffers = max(numThreadBlocks*ThreadBlockSize, nonbonded->getNumEnergyBuffers());
     if (useDoublePrecision) {
         energyBuffer = CudaArray::create<double>(*this, numEnergyBuffers, "energyBuffer");
+        energySum = CudaArray::create<double>(*this, 1, "energySum");
         int pinnedBufferSize = max(paddedNumAtoms*4, numEnergyBuffers);
         CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(double), 0));
     }
     else if (useMixedPrecision) {
         energyBuffer = CudaArray::create<double>(*this, numEnergyBuffers, "energyBuffer");
+        energySum = CudaArray::create<double>(*this, 1, "energySum");
         int pinnedBufferSize = max(paddedNumAtoms*4, numEnergyBuffers);
         CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(double), 0));
     }
     else {
         energyBuffer = CudaArray::create<float>(*this, numEnergyBuffers, "energyBuffer");
+        energySum = CudaArray::create<float>(*this, 1, "energySum");
         int pinnedBufferSize = max(paddedNumAtoms*6, numEnergyBuffers);
         CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(float), 0));
     }
@@ -870,6 +885,18 @@ void CudaContext::clearAutoclearBuffers() {
     }
 }
 
+double CudaContext::reduceEnergy() {
+    int bufferSize = energyBuffer->getSize();
+    int workGroupSize  = 512;
+    void* args[] = {&energyBuffer->getDevicePointer(), &energySum->getDevicePointer(), &bufferSize, &workGroupSize};
+    executeKernel(reduceEnergyKernel, args, workGroupSize, workGroupSize, workGroupSize*energyBuffer->getElementSize());
+    energySum->download(pinnedBuffer);
+    if (getUseDoublePrecision() || getUseMixedPrecision())
+        return *((double*) pinnedBuffer);
+    else
+        return *((float*) pinnedBuffer);
+}
+
 void CudaContext::setCharges(const vector<double>& charges) {
     if (chargeBuffer == NULL)
         chargeBuffer = new CudaArray(*this, numAtoms, useDoublePrecision ? sizeof(double) : sizeof(float), "chargeBuffer");

platforms/cuda/src/CudaKernelFactory.cpp

@@ -108,6 +108,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcCustomCentroidBondForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcCustomCompoundBondForceKernel::Name())
         return new CudaCalcCustomCompoundBondForceKernel(name, platform, cu, context.getSystem());
+    if (name == CalcCustomCVForceKernel::Name())
+        return new CudaCalcCustomCVForceKernel(name, platform, cu);
     if (name == CalcCustomManyParticleForceKernel::Name())
         return new CudaCalcCustomManyParticleForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcGayBerneForceKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -115,21 +115,8 @@ double CudaCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bo
     for (auto computation : cu.getPostComputations())
         sum += computation->computeForceAndEnergy(includeForces, includeEnergy, groups);
     cu.getIntegrationUtilities().distributeForcesFromVirtualSites();
-    if (includeEnergy) {
-        CudaArray& energyArray = cu.getEnergyBuffer();
-        if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
-            double* energy = (double*) cu.getPinnedBuffer();
-            energyArray.download(energy);
-            for (int i = 0; i < energyArray.getSize(); i++)
-                sum += energy[i];
-        }
-        else {
-            float* energy = (float*) cu.getPinnedBuffer();
-            energyArray.download(energy);
-            for (int i = 0; i < energyArray.getSize(); i++)
-                sum += energy[i];
-        }
-    }
+    if (includeEnergy)
+        sum += cu.reduceEnergy();
     if (!cu.getForcesValid())
         valid = false;
     return sum;
@@ -6609,6 +6596,176 @@ void CudaCalcGayBerneForceKernel::sortAtoms() {
     exclusionStartIndex->upload(startIndexVec);
 }
 
+class CudaCalcCustomCVForceKernel::ReorderListener : public CudaContext::ReorderListener {
+public:
+    ReorderListener(CudaContext& cu, CudaArray& invAtomOrder) : cu(cu), invAtomOrder(invAtomOrder) {
+    }
+    void execute() {
+        vector<int> invOrder(cu.getPaddedNumAtoms());
+        const vector<int>& order = cu.getAtomIndex();
+        for (int i = 0; i < order.size(); i++)
+            invOrder[order[i]] = i;
+        invAtomOrder.upload(invOrder);
+    }
+private:
+    CudaContext& cu;
+    CudaArray& invAtomOrder;
+};
+
+CudaCalcCustomCVForceKernel::~CudaCalcCustomCVForceKernel() {
+    for (auto force : cvForces)
+        delete force;
+    if (invAtomOrder != NULL)
+        delete invAtomOrder;
+    if (innerInvAtomOrder != NULL)
+        delete innerInvAtomOrder;
+}
+
+void CudaCalcCustomCVForceKernel::initialize(const System& system, const CustomCVForce& force, ContextImpl& innerContext) {
+    int numCVs = force.getNumCollectiveVariables();
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+    
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < (int) force.getNumTabulatedFunctions(); i++)
+        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
+
+    // Create the expressions.
+
+    Lepton::ParsedExpression energyExpr = Lepton::Parser::parse(force.getEnergyFunction(), functions);
+    energyExpression = energyExpr.createProgram();
+    for (int i = 0; i < numCVs; i++) {
+        string name = force.getCollectiveVariableName(i);
+        variableNames.push_back(name);
+        variableDerivExpressions.push_back(energyExpr.differentiate(name).optimize().createProgram());
+    }
+    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
+        string name = force.getEnergyParameterDerivativeName(i);
+        paramDerivNames.push_back(name);
+        paramDerivExpressions.push_back(energyExpr.differentiate(name).optimize().createProgram());
+        cu.addEnergyParameterDerivative(name);
+    }
+
+    // Delete the custom functions.
+
+    for (auto& function : functions)
+        delete function.second;
+        
+    // Copy parameter derivatives from the inner context.
+
+    CudaContext& cu2 = *reinterpret_cast<CudaPlatform::PlatformData*>(innerContext.getPlatformData())->contexts[0];
+    for (auto& param : cu2.getEnergyParamDerivNames())
+        cu.addEnergyParameterDerivative(param);
+    
+    // Create arrays for storing information.
+    
+    int elementSize = (cu.getUseDoublePrecision() || cu.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
+    for (int i = 0; i < numCVs; i++)
+        cvForces.push_back(CudaArray::create<long long>(cu, 3*cu.getPaddedNumAtoms(), "cvForce"));
+    invAtomOrder = CudaArray::create<int>(cu, cu.getPaddedNumAtoms(), "invAtomOrder");
+    innerInvAtomOrder = CudaArray::create<int>(cu, cu.getPaddedNumAtoms(), "innerInvAtomOrder");
+    
+    // Create the kernels.
+    
+    stringstream args, add;
+    for (int i = 0; i < numCVs; i++) {
+        args << ", long long* __restrict__ force" << i << ", real dEdV" << i;
+        add << "forces[i] += (long long) (force" << i << "[i]*dEdV" << i << ");\n";
+    }
+    map<string, string> replacements;
+    replacements["PARAMETER_ARGUMENTS"] = args.str();
+    replacements["ADD_FORCES"] = add.str();
+    CUmodule module = cu.createModule(cu.replaceStrings(CudaKernelSources::vectorOps+CudaKernelSources::customCVForce, replacements));
+    copyStateKernel = cu.getKernel(module, "copyState");
+    copyForcesKernel = cu.getKernel(module, "copyForces");
+    addForcesKernel = cu.getKernel(module, "addForces");
+}
+
+double CudaCalcCustomCVForceKernel::execute(ContextImpl& context, ContextImpl& innerContext, bool includeForces, bool includeEnergy) {
+    copyState(context, innerContext);
+    int numCVs = variableNames.size();
+    int numAtoms = cu.getNumAtoms();
+    int paddedNumAtoms = cu.getPaddedNumAtoms();
+    CudaContext& cu2 = *reinterpret_cast<CudaPlatform::PlatformData*>(innerContext.getPlatformData())->contexts[0];
+    vector<double> cvValues;
+    vector<map<string, double> > cvDerivs(numCVs);
+    void* copyForcesArgs[] = {NULL, &invAtomOrder->getDevicePointer(), &cu2.getForce().getDevicePointer(), &cu2.getAtomIndexArray().getDevicePointer(), &numAtoms, &paddedNumAtoms};
+    for (int i = 0; i < numCVs; i++) {
+        cvValues.push_back(innerContext.calcForcesAndEnergy(true, true, 1<<i));
+        copyForcesArgs[0] = &cvForces[i]->getDevicePointer();
+        cu.executeKernel(copyForcesKernel, copyForcesArgs, numAtoms);
+        innerContext.getEnergyParameterDerivatives(cvDerivs[i]);
+    }
+    
+    // Compute the energy and forces.
+    
+    map<string, double> variables;
+    for (auto& name : globalParameterNames)
+        variables[name] = context.getParameter(name);
+    for (int i = 0; i < numCVs; i++)
+        variables[variableNames[i]] = cvValues[i];
+    double energy = energyExpression.evaluate(variables);
+    int bufferSize = cu.getForce().getSize();
+    vector<void*> addForcesArgs;
+    addForcesArgs.push_back(&cu.getForce().getDevicePointer());
+    addForcesArgs.push_back(&bufferSize);
+    vector<double> dEdV(numCVs);
+    vector<float> dEdVFloat(numCVs);
+    for (int i = 0; i < numCVs; i++) {
+        dEdV[i] = variableDerivExpressions[i].evaluate(variables);
+        dEdVFloat[i] = (float) dEdV[i];
+        addForcesArgs.push_back(&cvForces[i]->getDevicePointer());
+        if (cu.getUseDoublePrecision())
+            addForcesArgs.push_back(&dEdV[i]);
+        else
+            addForcesArgs.push_back(&dEdVFloat[i]);
+    }
+    cu.executeKernel(addForcesKernel, &addForcesArgs[0], numAtoms);
+    
+    // Compute the energy parameter derivatives.
+    
+    map<string, double>& energyParamDerivs = cu.getEnergyParamDerivWorkspace();
+    for (int i = 0; i < paramDerivExpressions.size(); i++)
+        energyParamDerivs[paramDerivNames[i]] += paramDerivExpressions[i].evaluate(variables);
+    for (int i = 0; i < numCVs; i++) {
+        double dEdV = variableDerivExpressions[i].evaluate(variables);
+        for (auto& deriv : cvDerivs[i])
+            energyParamDerivs[deriv.first] += dEdV*deriv.second;
+    }
+    return energy;
+}
+
+void CudaCalcCustomCVForceKernel::copyState(ContextImpl& context, ContextImpl& innerContext) {
+    int numAtoms = cu.getNumAtoms();
+    CudaContext& cu2 = *reinterpret_cast<CudaPlatform::PlatformData*>(innerContext.getPlatformData())->contexts[0];
+    if (!hasInitializedListeners) {
+        hasInitializedListeners = true;
+        
+        // Initialize the listeners.
+        
+        ReorderListener* listener1 = new ReorderListener(cu, *invAtomOrder);
+        ReorderListener* listener2 = new ReorderListener(cu2, *innerInvAtomOrder);
+        cu.addReorderListener(listener1);
+        cu2.addReorderListener(listener2);
+        listener1->execute();
+        listener2->execute();
+    }
+    CUdeviceptr posCorrection = (cu.getUseMixedPrecision() ? cu.getPosqCorrection().getDevicePointer() : 0);
+    CUdeviceptr posCorrection2 = (cu2.getUseMixedPrecision() ? cu2.getPosqCorrection().getDevicePointer() : 0);
+    void* copyStateArgs[] = {&cu.getPosq().getDevicePointer(), &posCorrection, &cu.getVelm().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(),
+        &cu2.getPosq().getDevicePointer(), &posCorrection2,& cu2.getVelm().getDevicePointer(), &innerInvAtomOrder->getDevicePointer(), &numAtoms};
+    cu.executeKernel(copyStateKernel, copyStateArgs, numAtoms);
+    Vec3 a, b, c;
+    context.getPeriodicBoxVectors(a, b, c);
+    innerContext.setPeriodicBoxVectors(a, b, c);
+    innerContext.setTime(context.getTime());
+    map<string, double> innerParameters = innerContext.getParameters();
+    for (auto& param : innerParameters)
+        innerContext.setParameter(param.first, context.getParameter(param.first));
+}
+
 CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
 }