CustomCPPForceImpl for writing forces in C++ (#4231)

* Implemented CustomCPPForceImpl

* Documentation for CustomCPPForceImpl

* Attempt at fixing Windows compilation error

* Improved documentation

docs-source/developerguide/09_customcppforceimpl.rst

+.. role:: code
+.. raw:: html
+
+    <style> .code {font-family:monospace;} </style>
+    <style> .caption {text-align:center;} </style>
+
+.. highlight:: c++
+
+
+.. _customcppforceimpl:
+
+Platform-Independent Code with CustomCPPForceImpl
+#################################################
+
+In most cases, Forces are written with platform-specific code that runs as
+efficiently as possible on the available hardware.  There are some situations
+where that is unnecessarily complicated.  For example, you might write a plugin
+that interfaces to an external library through its own public API.  Because all
+expensive calculations are done by the external library, there is little
+opportunity for any kind of platform-specific optimization.  All you can do is
+copy positions over to the external library and copy forces back, hopefully in
+a way that is not too inefficient, and preferably with as little code as
+possible.  CustomCPPForceImpl is a tool for doing that.
+
+To use it, you write your Force subclass in the usual way, providing whatever
+API is appropriate for it.  Then for the corresponding ForceImpl, simply
+subclass CustomCPPForceImpl and override :code:`computeForce()` to perform the
+calculation using platform-independent C++ code.  Nothing more is required.  It
+will automatically work on all platforms.  You do not need to write any GPU
+code, provide a Kernel or KernelFactory, define registration functions, or even
+create plugin libraries to be loaded dynamically.  The single library containing
+the Force and ForceImpl is all you need.
+
+Here is an example of what the code to implement your ForceImpl might look like.
+::
+
+    class ExampleForceImpl : public CustomCPPForceImpl {
+    public:
+        ExampleForceImpl(const ExampleForce& owner) : CustomCPPForceImpl(owner), owner(owner) {
+        }
+
+        double computeForce(ContextImpl& context, const vector<Vec3>& positions, vector<Vec3>& forces) {
+            // Compute the forces and energy here.  Store the forces into the
+            // vector and return the energy.
+        }
+
+        const ExampleForce& getOwner() const {
+            return owner;
+        }
+    private:
+        const ExampleForce& owner;
+    };

docs-source/developerguide/index.rst

@@ -18,6 +18,7 @@ OpenMM Developer's Guide
    06_opencl_platform.rst
    07_cuda_platform.rst
    08_common_compute.rst
+   09_customcppforceimpl.rst
 
 .. only:: html
 

docs-source/developerguide/license.rst

-Portions copyright (c) 2011-2020 Stanford University and the Authors
+Portions copyright (c) 2011-2023 Stanford University and the Authors
 
 Contributors: Peter Eastman
 

olla/include/openmm/kernels.h

@@ -67,6 +67,7 @@
 #include "openmm/NoseHooverIntegrator.h"
 #include "openmm/NoseHooverChain.h"
 #include "openmm/ATMForce.h"
+#include "openmm/internal/CustomCPPForceImpl.h"
 #include <iosfwd>
 #include <set>
 #include <string>
@@ -1684,8 +1685,33 @@ public:
     virtual void copyState(ContextImpl& context, ContextImpl& innerContext0, ContextImpl& innerContext1) = 0;
 };
 
-  
-
+/**
+ * This kernel is invoked by CustomCPPForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCustomCPPForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcCustomCPPForce";
+    }
+    CalcCustomCPPForceKernel(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 CustomCPPForceImpl this kernel will be used for
+     */
+    virtual void initialize(const System& system, CustomCPPForceImpl& force) = 0;
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @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, bool includeForces, bool includeEnergy) = 0;
+};
 
 } // namespace OpenMM
 

openmmapi/include/openmm/internal/CustomCPPForceImpl.h

+#ifndef OPENMM_CUSTOMCPPFORCEIMPL_H_
+#define OPENMM_CUSTOMCPPFORCEIMPL_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) 2023 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/Kernel.h"
+#include "openmm/Vec3.h"
+#include <string>
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This class provides an easy way of writing Forces that are implemented entirely
+ * in C++, and do not involve any platform-specific code.  To use it, have your
+ * ForceImpl extend this class, implementing computeForce() to compute the forces
+ * and energy in a platform-independent way.  You can optionally override other
+ * methods as well, such as if your Force defines global parameters or bonds
+ * between atoms.
+ * 
+ * Most Forces should not be implemented with this class.  Writing platform-specific
+ * code that can run on a GPU and interact directly with internal data structures
+ * usually leads to much faster performance.  This class is useful for special
+ * situations, such as plugins that interact with external libraries through their
+ * own public APIs.  In that case, there is little opportunity for optimization.
+ * This class allows you to write a single implementation that is automatically
+ * used for all platforms.
+ */
+
+class OPENMM_EXPORT CustomCPPForceImpl : public ForceImpl {
+public:
+    CustomCPPForceImpl(const Force& owner);
+    /**
+     * Subclasses may override this to do their own initialization, but they should
+     * also call CustomCPPForceImpl::initialize(context) to let the parent class
+     * initialize its kernel.
+     * 
+     * @param context    the context in which the system is being simulated
+     */
+    void initialize(ContextImpl& context);
+    /**
+     * Calculate the force on each particle generated by this ForceImpl and/or this ForceImpl's
+     * contribution to the potential energy of the system.
+     * 
+     * @param context    the context in which the system is being simulated
+     * @param positions  the position of every particle
+     * @param forces     store the force on each particle into this array
+     * @return this force's contribution to the potential energy of the system, or 0 if this
+     * force does not contribute to potential energy
+     */
+    virtual double computeForce(ContextImpl& context, const std::vector<Vec3>& positions, std::vector<Vec3>& forces) = 0;
+    /**
+     * Override this if the force updates the context state directly.  In general
+     * it should only make changes by invoking public methods of the ContextImpl.
+     * If you need to access platform-specific data structures, this class is
+     * probably not appropriate.
+     * 
+     * @param context        the context in which the system is being simulated
+     * @param forcesInvalid  if the state was modified in any way that might cause previously
+     *                       calculated forces to no longer be valid (such as modifying
+     *                       positions or parameters), the method should set this to true.
+     */
+    void updateContextState(ContextImpl& context, bool& forcesInvalid) {
+    }
+    /**
+     * Override this if the force defines global parameters.
+     */
+    std::map<std::string, double> getDefaultParameters() {
+        return {};
+    }
+    /**
+     * Override this if the force defines bonds between particles.
+     */
+    std::vector<std::pair<int, int> > getBondedParticles() const {
+        return {};
+    }
+    /**
+     * Override this if the force defines an updateParametersInContext() method.
+     */
+    void updateParametersInContext(ContextImpl& context) {
+    }
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
+    std::vector<std::string> getKernelNames();
+private:
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMCPPFORCEIMPL_H_*/

openmmapi/src/CustomCPPForceImpl.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-2021 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/internal/ContextImpl.h"
+#include "openmm/internal/CustomCPPForceImpl.h"
+#include "openmm/kernels.h"
+
+using namespace OpenMM;
+using namespace std;
+
+CustomCPPForceImpl::CustomCPPForceImpl(const Force& owner) {
+    forceGroup = owner.getForceGroup();
+}
+
+void CustomCPPForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcCustomCPPForceKernel::Name(), context);
+    kernel.getAs<CalcCustomCPPForceKernel>().initialize(context.getSystem(), *this);
+}
+
+double CustomCPPForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&(1<<forceGroup)) != 0)
+        return kernel.getAs<CalcCustomCPPForceKernel>().execute(context, includeForces, includeEnergy);
+    return 0.0;
+}
+
+vector<string> CustomCPPForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCustomCPPForceKernel::Name());
+    return names;
+}

platforms/common/include/openmm/common/CommonKernels.h

@@ -1654,7 +1654,57 @@ private:
     int numParticles;
 };
 
-  
+/**
+ * This kernel is invoked by CustomCPPForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcCustomCPPForceKernel : public CalcCustomCPPForceKernel {
+public:
+    CommonCalcCustomCPPForceKernel(std::string name, const Platform& platform, OpenMM::ContextImpl& contextImpl, ComputeContext& cc) :
+            CalcCustomCPPForceKernel(name, platform), contextImpl(contextImpl), cc(cc), force(NULL) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomCPPForceImpl this kernel will be used for
+     */
+    void initialize(const System& system, CustomCPPForceImpl& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @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, bool includeForces, bool includeEnergy);
+    /**
+     * The is called by the pre-computation to start the calculation running.
+     */
+    void beginComputation(bool includeForces, bool includeEnergy, int groups);
+    /**
+     * This is called by the worker thread to do the computation.
+     */
+    void executeOnWorkerThread(bool includeForces);
+    /**
+     * This is called by the post-computation to add the forces to the main array.
+     */
+    double addForces(bool includeForces, bool includeEnergy, int groups);
+private:
+    class ExecuteTask;
+    class StartCalculationPreComputation;
+    class AddForcesPostComputation;
+    OpenMM::ContextImpl& contextImpl;
+    ComputeContext& cc;
+    CustomCPPForceImpl* force;
+    ComputeArray forcesArray;
+    ComputeKernel addForcesKernel;
+    std::vector<Vec3> positionsVec, forcesVec;
+    std::vector<float> floatForces;
+    int forceGroupFlag;
+    double energy;
+};
+
 } // namespace OpenMM
 
 #endif /*OPENMM_COMMONKERNELS_H_*/

platforms/common/src/CommonKernels.cpp

@@ -7974,3 +7974,111 @@ void CommonCalcATMForceKernel::copyParametersToContext(ContextImpl& context, con
     displ1.upload(displVectorContext1);
     displ0.upload(displVectorContext0);
 }
+
+class CommonCalcCustomCPPForceKernel::StartCalculationPreComputation : public ComputeContext::ForcePreComputation {
+public:
+    StartCalculationPreComputation(CommonCalcCustomCPPForceKernel& owner) : owner(owner) {
+    }
+    void computeForceAndEnergy(bool includeForces, bool includeEnergy, int groups) {
+        owner.beginComputation(includeForces, includeEnergy, groups);
+    }
+    CommonCalcCustomCPPForceKernel& owner;
+};
+
+class CommonCalcCustomCPPForceKernel::ExecuteTask : public ComputeContext::WorkTask {
+public:
+    ExecuteTask(CommonCalcCustomCPPForceKernel& owner, bool includeForces) : owner(owner), includeForces(includeForces) {
+    }
+    void execute() {
+        owner.executeOnWorkerThread(includeForces);
+    }
+    CommonCalcCustomCPPForceKernel& owner;
+    bool includeForces;
+};
+
+class CommonCalcCustomCPPForceKernel::AddForcesPostComputation : public ComputeContext::ForcePostComputation {
+public:
+    AddForcesPostComputation(CommonCalcCustomCPPForceKernel& owner) : owner(owner) {
+    }
+    double computeForceAndEnergy(bool includeForces, bool includeEnergy, int groups) {
+        return owner.addForces(includeForces, includeEnergy, groups);
+    }
+    CommonCalcCustomCPPForceKernel& owner;
+};
+
+void CommonCalcCustomCPPForceKernel::initialize(const System& system, CustomCPPForceImpl& force) {
+    ContextSelector selector(cc);
+    this->force = &force;
+    int numParticles = system.getNumParticles();
+    forcesVec.resize(numParticles);
+    positionsVec.resize(numParticles);
+    floatForces.resize(3*numParticles);
+    int elementSize = (cc.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
+    forcesArray.initialize(cc, 3*numParticles, elementSize, "forces");
+    map<string, string> defines;
+    defines["NUM_ATOMS"] = cc.intToString(numParticles);
+    defines["PADDED_NUM_ATOMS"] = cc.intToString(cc.getPaddedNumAtoms());
+    ComputeProgram program = cc.compileProgram(CommonKernelSources::customCppForce, defines);
+    addForcesKernel = program->createKernel("addForces");
+    addForcesKernel->addArg(forcesArray);
+    addForcesKernel->addArg(cc.getLongForceBuffer());
+    addForcesKernel->addArg(cc.getAtomIndexArray());
+    forceGroupFlag = (1<<force.getOwner().getForceGroup());
+    cc.addPreComputation(new StartCalculationPreComputation(*this));
+    cc.addPostComputation(new AddForcesPostComputation(*this));
+}
+
+double CommonCalcCustomCPPForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    // This method does nothing.  The actual calculation is started by the pre-computation, continued on
+    // the worker thread, and finished by the post-computation.
+    
+    return 0;
+}
+
+void CommonCalcCustomCPPForceKernel::beginComputation(bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&forceGroupFlag) == 0)
+        return;
+    contextImpl.getPositions(positionsVec);
+    
+    // The actual force computation will be done on a different thread.
+    
+    cc.getWorkThread().addTask(new ExecuteTask(*this, includeForces));
+}
+
+void CommonCalcCustomCPPForceKernel::executeOnWorkerThread(bool includeForces) {
+    energy = force->computeForce(contextImpl, positionsVec, forcesVec);
+    if (includeForces) {
+        ContextSelector selector(cc);
+        int numParticles = cc.getNumAtoms();
+        if (cc.getUseDoublePrecision())
+            forcesArray.upload((double*) forcesVec.data());
+        else {
+            for (int i = 0; i < numParticles; i++) {
+                floatForces[3*i] = (float) forcesVec[i][0];
+                floatForces[3*i+1] = (float) forcesVec[i][1];
+                floatForces[3*i+2] = (float) forcesVec[i][2];
+            }
+            forcesArray.upload(floatForces);
+        }
+    }
+}
+
+double CommonCalcCustomCPPForceKernel::addForces(bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&forceGroupFlag) == 0)
+        return 0;
+
+    // Wait until executeOnWorkerThread() is finished.
+    
+    cc.getWorkThread().flush();
+
+    // Add in the forces.
+    
+    if (includeForces) {
+        ContextSelector selector(cc);
+        addForcesKernel->execute(cc.getNumAtoms());
+    }
+    
+    // Return the energy.
+    
+    return energy;
+}
\ No newline at end of file

platforms/common/src/kernels/customCppForce.cc

+KERNEL void addForces(GLOBAL const real* RESTRICT forces, GLOBAL mm_long* RESTRICT forceBuffers, GLOBAL int* RESTRICT atomIndex) {
+    for (int atom = GLOBAL_ID; atom < NUM_ATOMS; atom += GLOBAL_SIZE) {
+        int index = atomIndex[atom];
+        forceBuffers[atom] += (mm_long) (forces[3*index]*0x100000000);
+        forceBuffers[atom+PADDED_NUM_ATOMS] += (mm_long) (forces[3*index+1]*0x100000000);
+        forceBuffers[atom+2*PADDED_NUM_ATOMS] += (mm_long) (forces[3*index+2]*0x100000000);
+    }
+}

platforms/cpu/tests/TestCpuCustomCPPForce.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) 2023 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 "CpuTests.h"
+#include "TestCustomCPPForce.h"
+
+void runPlatformTests() {
+}
+

platforms/cuda/src/CudaKernelFactory.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2023 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -113,6 +113,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcCustomCVForceKernel(name, platform, cu);
     if (name == CalcATMForceKernel::Name())
         return new CudaCalcATMForceKernel(name, platform, cu);
+    if (name == CalcCustomCPPForceKernel::Name())
+        return new CommonCalcCustomCPPForceKernel(name, platform, context, cu);
     if (name == CalcRMSDForceKernel::Name())
         return new CommonCalcRMSDForceKernel(name, platform, cu);
     if (name == CalcCustomManyParticleForceKernel::Name())

platforms/cuda/src/CudaPlatform.cpp

@@ -91,6 +91,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomCPPForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
     registerKernelFactory(CalcATMForceKernel::Name(), factory);
     registerKernelFactory(CalcRMSDForceKernel::Name(), factory);

platforms/cuda/tests/TestCudaCustomCPPForce.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) 2023 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 "CudaTests.h"
+#include "TestCustomCPPForce.h"
+
+void runPlatformTests() {
+}
+

platforms/opencl/src/OpenCLKernelFactory.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2023 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -110,7 +110,9 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
     if (name == CalcCustomCVForceKernel::Name())
         return new OpenCLCalcCustomCVForceKernel(name, platform, cl);
     if (name == CalcATMForceKernel::Name())
-        return new OpenCLCalcATMForceKernel(name, platform, cl);    
+        return new OpenCLCalcATMForceKernel(name, platform, cl);
+    if (name == CalcCustomCPPForceKernel::Name())
+        return new CommonCalcCustomCPPForceKernel(name, platform, context, cl);
     if (name == CalcRMSDForceKernel::Name())
         return new CommonCalcRMSDForceKernel(name, platform, cl);
     if (name == CalcCustomManyParticleForceKernel::Name())

platforms/opencl/src/OpenCLPlatform.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-2018 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2023 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -82,6 +82,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomCPPForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
     registerKernelFactory(CalcATMForceKernel::Name(), factory);
     registerKernelFactory(CalcRMSDForceKernel::Name(), factory);

platforms/opencl/tests/TestOpenCLCustomCPPForce.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) 2023 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 "OpenCLTests.h"
+#include "TestCustomCPPForce.h"
+
+void runPlatformTests() {
+}
+

platforms/reference/include/ReferenceKernels.h

@@ -34,6 +34,7 @@
 
 #include "ReferencePlatform.h"
 #include "openmm/kernels.h"
+#include "openmm/internal/CustomCPPForceImpl.h"
 #include "openmm/internal/CustomNonbondedForceImpl.h"
 #include "SimTKOpenMMRealType.h"
 #include "ReferenceNeighborList.h"
@@ -1680,6 +1681,34 @@ private:
     std::vector<Vec3> displ0;
 };
 
+/**
+ * This kernel is invoked by CustomCPPForceImpl to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcCustomCPPForceKernel : public CalcCustomCPPForceKernel {
+public:
+    ReferenceCalcCustomCPPForceKernel(std::string name, const Platform& platform) : CalcCustomCPPForceKernel(name, platform), force(NULL) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomCPPForceImpl this kernel will be used for
+     */
+    void initialize(const System& system, CustomCPPForceImpl& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @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, bool includeForces, bool includeEnergy);
+private:
+    CustomCPPForceImpl* force;
+    std::vector<Vec3> forces;
+};
+
 } // namespace OpenMM
 
 #endif /*OPENMM_REFERENCEKERNELS_H_*/

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -82,6 +82,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomCVForceKernel(name, platform);
     if (name == CalcATMForceKernel::Name())
         return new ReferenceCalcATMForceKernel(name, platform);
+    if (name == CalcCustomCPPForceKernel::Name())
+        return new ReferenceCalcCustomCPPForceKernel(name, platform);
     if (name == CalcRMSDForceKernel::Name())
         return new ReferenceCalcRMSDForceKernel(name, platform);
     if (name == CalcCustomManyParticleForceKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -3032,3 +3032,18 @@ void ReferenceCalcATMForceKernel::copyParametersToContext(ContextImpl& context,
         displ0[i] = displacement0;
     }
 }
+
+void ReferenceCalcCustomCPPForceKernel::initialize(const System& system, CustomCPPForceImpl& force) {
+    this->force = &force;
+    forces.resize(system.getNumParticles());
+}
+
+double ReferenceCalcCustomCPPForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    vector<Vec3>& posData = extractPositions(context);
+    vector<Vec3>& forceData = extractForces(context);
+    double energy = force->computeForce(context, posData, forces);
+    if (includeForces)
+        for (int i = 0; i < forces.size(); i++)
+            forceData[i] += forces[i];
+    return energy;
+}

platforms/reference/src/ReferencePlatform.cpp

@@ -61,6 +61,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomCPPForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
     registerKernelFactory(CalcATMForceKernel::Name(), factory);
     registerKernelFactory(CalcRMSDForceKernel::Name(), factory);