LocalEnergyMinimizer works with constraints. Also added Context::applyConstraints().

olla/include/openmm/kernels.h

@@ -191,6 +191,31 @@ public:
     virtual void setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) const = 0;
 };
 
+/**
+ * This kernel modifies the positions of particles to enforce distance constraints.
+ */
+class ApplyConstraintsKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "ApplyConstraints";
+    }
+    ApplyConstraintsKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     */
+    virtual void initialize(const System& system) = 0;
+    /**
+     * Update particle positions to enforce constraints.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param tol        the distance tolerance within which constraints must be satisfied.
+     */
+    virtual void apply(ContextImpl& context, double tol) = 0;
+};
+
 /**
  * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
  */

openmmapi/include/openmm/Context.h

@@ -127,7 +127,9 @@ public:
      */
     void setTime(double time);
     /**
-     * Set the positions of all particles in the System (measured in nm).
+     * Set the positions of all particles in the System (measured in nm).  This method simply sets the positions
+     * without checking to see whether they satisfy distance constraints.  If you want constraints to be
+     * enforced, call applyConstraints() after setting the positions.
      * 
      * @param positions   a vector whose length equals the number of particles in the System.  The i'th element
      * contains the position of the i'th particle.
@@ -165,6 +167,12 @@ public:
      * @param c      the vector defining the third edge of the periodic box
      */
     void setPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c);
+    /**
+     * Update the positions of particles so that all distance constraints are satisfied.
+     *
+     * @param tol    the distance tolerance within which constraints must be satisfied.
+     */
+    void applyConstraints(double tol);
     /**
      * When a Context is created, it may cache information about the System being simulated
      * and the Force objects contained in it.  This means that, if the System or Forces are then

openmmapi/include/openmm/LocalEnergyMinimizer.h

@@ -39,6 +39,10 @@ namespace OpenMM {
 /**
  * Given a Context, this class searches for a new set of particle positions that represent
  * a local minimum of the potential energy.  The search is performed with the L-BFGS algorithm.
+ * Distance constraints are enforced during minimization by adding a harmonic restraining
+ * force to the potential function.  The strength of the restraining force is steadily increased
+ * until the minimum energy configuration satisfies all constraints to within the tolerance
+ * specified by the Context's Integrator.
  */
 
 class OPENMM_EXPORT LocalEnergyMinimizer {

openmmapi/include/openmm/internal/ContextImpl.h

@@ -157,6 +157,12 @@ public:
      * @param c      the vector defining the third edge of the periodic box
      */
     void setPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c);
+    /**
+     * Update the positions of particles so that all distance constraints are satisfied.
+     *
+     * @param tol    the distance tolerance within which constraints must be satisfied.
+     */
+    void applyConstraints(double tol);
     /**
      * Recalculate all of the forces in the system.  After calling this, use getForces() to retrieve
      * the forces that were calculated.
@@ -207,7 +213,7 @@ private:
     mutable std::vector<std::vector<int> > molecules;
     bool hasInitializedForces;
     Platform* platform;
-    Kernel initializeForcesKernel, kineticEnergyKernel, updateStateDataKernel;
+    Kernel initializeForcesKernel, kineticEnergyKernel, updateStateDataKernel, applyConstraintsKernel;
     void* platformData;
 };
 

openmmapi/src/Context.cpp

@@ -127,6 +127,10 @@ void Context::setPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c)
     impl->setPeriodicBoxVectors(a, b, c);
 }
 
+void Context::applyConstraints(double tol) {
+    impl->applyConstraints(tol);
+}
+
 void Context::reinitialize() {
     System& system = impl->getSystem();
     Integrator& integrator = impl->getIntegrator();

openmmapi/src/ContextImpl.cpp

@@ -73,6 +73,8 @@ ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator,
     dynamic_cast<CalcKineticEnergyKernel&>(kineticEnergyKernel.getImpl()).initialize(system);
     updateStateDataKernel = platform->createKernel(UpdateStateDataKernel::Name(), *this);
     dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).initialize(system);
+    applyConstraintsKernel = platform->createKernel(ApplyConstraintsKernel::Name(), *this);
+    dynamic_cast<ApplyConstraintsKernel&>(applyConstraintsKernel.getImpl()).initialize(system);
     Vec3 periodicBoxVectors[3];
     system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
     dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).setPeriodicBoxVectors(*this, periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
@@ -142,6 +144,10 @@ void ContextImpl::setPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3
     dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).setPeriodicBoxVectors(*this, a, b, c);
 }
 
+void ContextImpl::applyConstraints(double tol) {
+    dynamic_cast<ApplyConstraintsKernel&>(applyConstraintsKernel.getImpl()).apply(*this, tol);
+}
+
 void ContextImpl::calcForces() {
     CalcForcesAndEnergyKernel& kernel = dynamic_cast<CalcForcesAndEnergyKernel&>(initializeForcesKernel.getImpl());
     kernel.beginForceComputation(*this);

openmmapi/src/LocalEnergyMinimizer.cpp

@@ -40,32 +40,86 @@
 using namespace OpenMM;
 using namespace std;
 
+typedef struct {
+    Context& context;
+    double k;
+} MinimizerData;
+
 static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsfloatval_t *g, const int n, const lbfgsfloatval_t step) {
-    Context* context = reinterpret_cast<Context*>(instance);
-    int numParticles = context->getSystem().getNumParticles();
+    MinimizerData* data = reinterpret_cast<MinimizerData*>(instance);
+    Context& context = data->context;
+    const System& system = context.getSystem();
+    int numParticles = system.getNumParticles();
+
+    // Compute the force and energy for this configuration.
+
     vector<Vec3> positions(numParticles);
     for (int i = 0; i < numParticles; i++)
         positions[i] = Vec3(x[3*i], x[3*i+1], x[3*i+2]);
-    context->setPositions(positions);
-    State state = context->getState(State::Forces | State::Energy);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
     const vector<Vec3>& forces = state.getForces();
     for (int i = 0; i < numParticles; i++) {
         g[3*i] = -forces[i][0];
         g[3*i+1] = -forces[i][1];
         g[3*i+2] = -forces[i][2];
     }
-    return state.getPotentialEnergy();
+    double energy = state.getPotentialEnergy();
+
+    // Add harmonic forces for any constraints.
+
+    int numConstraints = system.getNumConstraints();
+    double k = data->k;
+    for (int i = 0; i < numConstraints; i++) {
+        int particle1, particle2;
+        double distance;
+        system.getConstraintParameters(i, particle1, particle2, distance);
+        Vec3 delta = positions[particle2]-positions[particle1];
+        double r2 = delta.dot(delta);
+        double r = sqrt(r2);
+        delta *= 1/r;
+        double dr = r-distance;
+        double kdr = k*dr;
+        energy += 0.5*kdr*dr;
+        g[3*particle1] -= kdr*delta[0];
+        g[3*particle1+1] -= kdr*delta[1];
+        g[3*particle1+2] -= kdr*delta[2];
+        g[3*particle2] += kdr*delta[0];
+        g[3*particle2+1] += kdr*delta[1];
+        g[3*particle2+2] += kdr*delta[2];
+    }
+    return energy;
 }
 
 void LocalEnergyMinimizer::minimize(Context& context, double tolerance, int maxIterations) {
-    int numParticles = context.getSystem().getNumParticles();
+    System& system = context.getSystem();
+    int numParticles = system.getNumParticles();
     lbfgsfloatval_t *x = lbfgs_malloc(numParticles*3);
     if (x == NULL)
         throw OpenMMException("LocalEnergyMinimizer: Failed to allocate memory");
+    double constraintTol = context.getIntegrator().getConstraintTolerance();
+    double k = tolerance/constraintTol;
+
+    // Initialize the minimizer.
+
+    lbfgs_parameter_t param;
+    lbfgs_parameter_init(&param);
+    if (!context.getPlatform().supportsDoublePrecision())
+        param.xtol = 1e-7;
+    param.max_iterations = maxIterations;
+    param.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
+
+    // Repeatedly minimize, steadily increasing the strength of the springs until all constraints are satisfied.
+
+    double prevMaxError = 1e10;
+    while (true) {
+        // Make sure the initial configuration satisfies all constraints.
+
+        context.applyConstraints(constraintTol);
 
         // Record the initial positions and determine a normalization constant for scaling the tolerance.
 
-    const vector<Vec3>& positions = context.getState(State::Positions).getPositions();
+        vector<Vec3> positions = context.getState(State::Positions).getPositions();
         double norm = 0.0;
         for (int i = 0; i < numParticles; i++) {
             x[3*i] = positions[i][0];
@@ -75,18 +129,41 @@ void LocalEnergyMinimizer::minimize(Context& context, double tolerance, int maxI
         }
         norm /= numParticles;
         norm = (norm < 1 ? 1 : sqrt(norm));
+        param.epsilon = tolerance/norm;
 
         // Perform the minimization.
 
-    lbfgs_parameter_t param;
-    lbfgs_parameter_init(&param);
-    if (!context.getPlatform().supportsDoublePrecision())
-        param.xtol = 1e-7;
-    param.epsilon = tolerance/norm;
-    param.max_iterations = maxIterations;
-    param.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
         lbfgsfloatval_t fx;
-    lbfgs(numParticles*3, x, &fx, evaluate, NULL, &context, &param);
+        MinimizerData data = (MinimizerData) {context, k};
+        lbfgs(numParticles*3, x, &fx, evaluate, NULL, &data, &param);
+
+        // Check whether all constraints are satisfied.
+
+        positions = context.getState(State::Positions).getPositions();
+        int numConstraints = system.getNumConstraints();
+        double maxError = 0.0;
+        for (int i = 0; i < numConstraints; i++) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(i, particle1, particle2, distance);
+            Vec3 delta = positions[particle2]-positions[particle1];
+            double r = sqrt(delta.dot(delta));
+            double error = fabs(r-distance);
+            if (error > maxError)
+                maxError = error;
+        }
+        if (maxError <= constraintTol)
+            break; // All constraints are satisfied.
+        if (maxError >= prevMaxError) {
+            // Further tightening the springs doesn't seem to be helping, so just to a final
+            // constraint application and return.
+
+            context.applyConstraints(constraintTol);
+            break;
+        }
+        prevMaxError = maxError;
+        k *= 10;
+    }
     lbfgs_free(x);
 }
 

platforms/cuda/src/CudaKernelFactory.cpp

@@ -37,6 +37,8 @@ OPENMMCUDA_EXPORT KernelImpl* CudaKernelFactory::createKernelImpl(std::string na
         return new CudaCalcForcesAndEnergyKernel(name, platform, data);
     if (name == UpdateStateDataKernel::Name())
         return new CudaUpdateStateDataKernel(name, platform, data);
+    if (name == ApplyConstraintsKernel::Name())
+        return new CudaApplyConstraintsKernel(name, platform, data);
     if (name == CalcHarmonicBondForceKernel::Name())
         return new CudaCalcHarmonicBondForceKernel(name, platform, data, context.getSystem());
     if (name == CalcCustomBondForceKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -199,6 +199,13 @@ void CudaUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, cons
     gpuSetConstants(gpu);
 }
 
+void CudaApplyConstraintsKernel::initialize(const System& system) {
+}
+
+void CudaApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
+    kApplyConstraints(data.gpu);
+}
+
 CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
 }
 

platforms/cuda/src/CudaKernels.h

@@ -169,6 +169,30 @@ private:
     CudaPlatform::PlatformData& data;
 };
 
+/**
+ * This kernel modifies the positions of particles to enforce distance constraints.
+ */
+class CudaApplyConstraintsKernel : public ApplyConstraintsKernel {
+public:
+    CudaApplyConstraintsKernel(std::string name, const Platform& platform, CudaPlatform::PlatformData& data) : ApplyConstraintsKernel(name, platform), data(data) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     */
+    void initialize(const System& system);
+    /**
+     * Update particle positions to enforce constraints.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param tol        the distance tolerance within which constraints must be satisfied.
+     */
+    void apply(ContextImpl& context, double tol);
+private:
+    CudaPlatform::PlatformData& data;
+};
+
 /**
  * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/cuda/src/CudaPlatform.cpp

@@ -47,6 +47,7 @@ CudaPlatform::CudaPlatform() {
     CudaKernelFactory* factory = new CudaKernelFactory();
     registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
     registerKernelFactory(UpdateStateDataKernel::Name(), factory);
+    registerKernelFactory(ApplyConstraintsKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomBondForceKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);

platforms/cuda/src/kernels/cudaKernels.h

@@ -69,6 +69,7 @@ extern void kSelectVerletStepSize(gpuContext gpu, float maxTimeStep);
 extern void kBrownianUpdatePart1(gpuContext gpu);
 extern void kBrownianUpdatePart2(gpuContext gpu);
 extern void kScaleAtomCoordinates(gpuContext gpu, float scale, CUDAStream<int>& moleculeAtoms, CUDAStream<int>& moleculeStartIndex);
+extern void kApplyConstraints(gpuContext gpu);
 
 // Extras
 extern void kReduceForces(gpuContext gpu);

platforms/cuda/src/kernels/kApplyConstraints.cu

+/* -------------------------------------------------------------------------- *
+ *                                   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) 2010 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include <stdio.h>
+#include <cuda.h>
+#include <vector_functions.h>
+#include <cstdlib>
+#include <string>
+#include <iostream>
+//#include <fstream>
+using namespace std;
+
+#include "gputypes.h"
+#include "cudaKernels.h"
+
+__global__ void kPrepareConstraints_kernel(int numAtoms, float4* oldPos, float4* posq, float4* posqP) {
+    for (int index = threadIdx.x+blockIdx.x*blockDim.x; index < numAtoms; index += blockDim.x*gridDim.x) {
+        float4 pos = posq[index];
+        oldPos[index] = pos;
+        posqP[index] = make_float4(0.0f, 0.0f, 0.0f, pos.w);
+    }
+}
+
+__global__ void kFinishConstraints_kernel(int numAtoms, float4* posq, float4* posqP) {
+    for (int index = threadIdx.x+blockIdx.x*blockDim.x; index < numAtoms; index += blockDim.x*gridDim.x) {
+        float4 pos = posq[index];
+        float4 delta = posqP[index];
+        posq[index] = make_float4(pos.x+delta.x, pos.y+delta.y, pos.z+delta.z, pos.w);
+    }
+}
+
+void kApplyConstraints(gpuContext gpu)
+{
+    kPrepareConstraints_kernel<<<gpu->sim.blocks, gpu->sim.update_threads_per_block>>>(gpu->natoms, gpu->sim.pOldPosq, gpu->sim.pPosq, gpu->sim.pPosqP);
+    LAUNCHERROR("kPrepareConstraints");
+    kApplyFirstShake(gpu);
+    kApplyFirstSettle(gpu);
+    kApplyFirstCCMA(gpu);
+    kFinishConstraints_kernel<<<gpu->sim.blocks, gpu->sim.update_threads_per_block>>>(gpu->natoms, gpu->sim.pPosq, gpu->sim.pPosqP);
+    LAUNCHERROR("kFinishConstraints");
+}
+

tests/TestLocalEnergyMinimizer.cpp → platforms/cuda/tests/TestCudaLocalEnergyMinimizer.cpp

@@ -30,7 +30,8 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "AssertionUtilities.h"
+#include "../../../tests/AssertionUtilities.h"
+#include "CudaPlatform.h"
 #include "openmm/Context.h"
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/LocalEnergyMinimizer.h"
@@ -62,7 +63,8 @@ void testHarmonicBonds() {
     // Minimize it and check that all bonds are at their equilibrium distances.
 
     VerletIntegrator integrator(0.01);
-    Context context(system, integrator);
+    CudaPlatform platform;
+    Context context(system, integrator, platform);
     context.setPositions(positions);
     LocalEnergyMinimizer::minimize(context, 1e-5);
     State state = context.getState(State::Positions);
@@ -73,7 +75,8 @@ void testHarmonicBonds() {
 }
 
 void testLargeSystem() {
-    const int numParticles = 100;
+    const int numMolecules = 50;
+    const int numParticles = numMolecules*2;
     const double cutoff = 2.0;
     const double boxSize = 5.0;
     const double tolerance = 5;
@@ -84,36 +87,49 @@ void testLargeSystem() {
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     system.addForce(nonbonded);
 
-    // Create a cloud of particles.
+    // Create a cloud of molecules.
 
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++) {
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
         system.addParticle(1.0);
-        nonbonded->addParticle(i%2 == 0 ? 1 : -1, 0.2, 0.2);
-        positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        nonbonded->addParticle(-1.0, 0.2, 0.2);
+        nonbonded->addParticle(1.0, 0.2, 0.2);
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        system.addConstraint(2*i, 2*i+1, 1.0);
     }
 
-    // Minimize it, verify that the energy has decreased, and check that the force magnitude satisfies the requested tolerance.
+    // Minimize it and verify that the energy has decreased.
 
+    CudaPlatform platform;
     VerletIntegrator integrator(0.01);
-    Context context(system, integrator);
+    Context context(system, integrator, platform);
     context.setPositions(positions);
     State initialState = context.getState(State::Forces | State::Energy);
     LocalEnergyMinimizer::minimize(context, tolerance);
-    State finalState = context.getState(State::Forces | State::Energy);
+    State finalState = context.getState(State::Forces | State::Energy | State::Positions);
     ASSERT(finalState.getPotentialEnergy() < initialState.getPotentialEnergy());
-    double initialNorm = 0.0;
-    double finalNorm = 0.0;
-    for (int i = 0; i < numParticles; i++) {
-        initialNorm += initialState.getForces()[i].dot(initialState.getForces()[i]);
-        finalNorm += finalState.getForces()[i].dot(finalState.getForces()[i]);
+
+    // Compute the force magnitude, substracting off any component parallel to a constraint, and
+    // check that it satisfies the requested tolerance.
+
+    double forceNorm = 0.0;
+    for (int i = 0; i < numParticles; i += 2) {
+        Vec3 dir = finalState.getPositions()[i+1]-finalState.getPositions()[i];
+        double distance = sqrt(dir.dot(dir));
+        dir *= 1.0/distance;
+        Vec3 f = finalState.getForces()[i];
+        f -= dir*dir.dot(f);
+        forceNorm += f.dot(f);
+        f = finalState.getForces()[i+1];
+        f -= dir*dir.dot(f);
+        forceNorm += f.dot(f);
     }
-    initialNorm = sqrt(initialNorm/numParticles);
-    finalNorm = sqrt(finalNorm/numParticles);
-    ASSERT(finalNorm < initialNorm);
-    ASSERT(finalNorm < 2*tolerance);
+    forceNorm = sqrt(forceNorm/(4*numMolecules));
+    ASSERT(forceNorm < 3*tolerance);
 }
 
 int main() {
@@ -128,3 +144,4 @@ int main() {
     cout << "Done" << endl;
     return 0;
 }
+

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -38,6 +38,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLCalcForcesAndEnergyKernel(name, platform, cl);
     if (name == UpdateStateDataKernel::Name())
         return new OpenCLUpdateStateDataKernel(name, platform, cl);
+    if (name == ApplyConstraintsKernel::Name())
+        return new OpenCLApplyConstraintsKernel(name, platform, cl);
     if (name == CalcHarmonicBondForceKernel::Name())
         return new OpenCLCalcHarmonicBondForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcCustomBondForceKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -190,6 +190,13 @@ void OpenCLUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, co
     cl.setPeriodicBoxSize(a[0], b[1], c[2]);
 }
 
+void OpenCLApplyConstraintsKernel::initialize(const System& system) {
+}
+
+void OpenCLApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
+    cl.getIntegrationUtilities().applyConstraints(tol);
+}
+
 class OpenCLBondForceInfo : public OpenCLForceInfo {
 public:
     OpenCLBondForceInfo(int requiredBuffers, const HarmonicBondForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {

platforms/opencl/src/OpenCLKernels.h

@@ -164,6 +164,30 @@ private:
     OpenCLContext& cl;
 };
 
+/**
+ * This kernel modifies the positions of particles to enforce distance constraints.
+ */
+class OpenCLApplyConstraintsKernel : public ApplyConstraintsKernel {
+public:
+    OpenCLApplyConstraintsKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyConstraintsKernel(name, platform), cl(cl) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     */
+    void initialize(const System& system);
+    /**
+     * Update particle positions to enforce constraints.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param tol        the distance tolerance within which constraints must be satisfied.
+     */
+    void apply(ContextImpl& context, double tol);
+private:
+    OpenCLContext& cl;
+};
+
 /**
  * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/opencl/src/OpenCLPlatform.cpp

@@ -46,6 +46,7 @@ OpenCLPlatform::OpenCLPlatform() {
     OpenCLKernelFactory* factory = new OpenCLKernelFactory();
     registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
     registerKernelFactory(UpdateStateDataKernel::Name(), factory);
+    registerKernelFactory(ApplyConstraintsKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomBondForceKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);

platforms/opencl/tests/TestOpenCLLocalEnergyMinimizer.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) 2010 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 "../../../tests/AssertionUtilities.h"
+#include "OpenCLPlatform.h"
+#include "openmm/Context.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/LocalEnergyMinimizer.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testHarmonicBonds() {
+    const int numParticles = 10;
+    System system;
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+
+    // Create a chain of particles connected by harmonic bonds.
+
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i] = Vec3(i, 0, 0);
+        if (i > 0)
+            bonds->addBond(i-1, i, 1+0.1*i, 1);
+    }
+
+    // Minimize it and check that all bonds are at their equilibrium distances.
+
+    VerletIntegrator integrator(0.01);
+    OpenCLPlatform platform;
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    LocalEnergyMinimizer::minimize(context, 1e-5);
+    State state = context.getState(State::Positions);
+    for (int i = 1; i < numParticles; i++) {
+        Vec3 delta = state.getPositions()[i]-state.getPositions()[i-1];
+        ASSERT_EQUAL_TOL(1+0.1*i, sqrt(delta.dot(delta)), 1e-4);
+    }
+}
+
+void testLargeSystem() {
+    const int numMolecules = 50;
+    const int numParticles = numMolecules*2;
+    const double cutoff = 2.0;
+    const double boxSize = 5.0;
+    const double tolerance = 5;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        nonbonded->addParticle(-1.0, 0.2, 0.2);
+        nonbonded->addParticle(1.0, 0.2, 0.2);
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        system.addConstraint(2*i, 2*i+1, 1.0);
+    }
+
+    // Minimize it and verify that the energy has decreased.
+
+    OpenCLPlatform platform;
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State initialState = context.getState(State::Forces | State::Energy);
+    LocalEnergyMinimizer::minimize(context, tolerance);
+    State finalState = context.getState(State::Forces | State::Energy | State::Positions);
+    ASSERT(finalState.getPotentialEnergy() < initialState.getPotentialEnergy());
+
+    // Compute the force magnitude, substracting off any component parallel to a constraint, and
+    // check that it satisfies the requested tolerance.
+
+    double forceNorm = 0.0;
+    for (int i = 0; i < numParticles; i += 2) {
+        Vec3 dir = finalState.getPositions()[i+1]-finalState.getPositions()[i];
+        double distance = sqrt(dir.dot(dir));
+        dir *= 1.0/distance;
+        Vec3 f = finalState.getForces()[i];
+        f -= dir*dir.dot(f);
+        forceNorm += f.dot(f);
+        f = finalState.getForces()[i+1];
+        f -= dir*dir.dot(f);
+        forceNorm += f.dot(f);
+    }
+    forceNorm = sqrt(forceNorm/(4*numMolecules));
+    ASSERT(forceNorm < 3*tolerance);
+}
+
+int main() {
+    try {
+        testHarmonicBonds();
+        testLargeSystem();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -42,6 +42,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcForcesAndEnergyKernel(name, platform);
     if (name == UpdateStateDataKernel::Name())
         return new ReferenceUpdateStateDataKernel(name, platform, data);
+    if (name == ApplyConstraintsKernel::Name())
+        return new ReferenceApplyConstraintsKernel(name, platform, data);
     if (name == CalcNonbondedForceKernel::Name())
         return new ReferenceCalcNonbondedForceKernel(name, platform);
     if (name == CalcCustomNonbondedForceKernel::Name())