Parallelized SETTLE on CPU platform

platforms/cpu/include/CpuSETTLE.h

+#ifndef OPENMM_CPUSETTLE_H_
+#define OPENMM_CPUSETTLE_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) 2013 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 "ReferenceSETTLEAlgorithm.h"
+#include "openmm/System.h"
+#include "openmm/internal/ThreadPool.h"
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This class uses multiple ReferenceSETTLEAlgorithm objects to execute the algorithm in parallel.
+ */
+class OPENMM_EXPORT CpuSETTLE : public ReferenceConstraintAlgorithm {
+public:
+    class ApplyToPositionsTask;
+    class ApplyToVelocitiesTask;
+    CpuSETTLE(const System& system, const ReferenceSETTLEAlgorithm& settle, ThreadPool& threads);
+    ~CpuSETTLE();
+
+    /**
+     * Apply the constraint algorithm.
+     * 
+     * @param atomCoordinates  the original atom coordinates
+     * @param atomCoordinatesP the new atom coordinates
+     * @param inverseMasses    1/mass
+     * @param tolerance        the constraint tolerance
+     */
+    void apply(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
+
+    /**
+     * Apply the constraint algorithm to velocities.
+     * 
+     * @param atomCoordinates  the atom coordinates
+     * @param atomCoordinatesP the velocities to modify
+     * @param inverseMasses    1/mass
+     * @param tolerance        the constraint tolerance
+     */
+    void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
+private:
+    std::vector<ReferenceSETTLEAlgorithm*> threadSettle;
+    ThreadPool& threads;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CPUSETTLE_H_*/

platforms/cpu/src/CpuPlatform.cpp

@@ -32,6 +32,8 @@
 #include "CpuPlatform.h"
 #include "CpuKernelFactory.h"
 #include "CpuKernels.h"
+#include "CpuSETTLE.h"
+#include "ReferenceConstraints.h"
 #include "openmm/internal/hardware.h"
 
 using namespace OpenMM;
@@ -77,6 +79,12 @@ void CpuPlatform::contextCreated(ContextImpl& context, const map<string, string>
     ReferencePlatform::contextCreated(context, properties);
     PlatformData* data = new PlatformData(context.getSystem().getNumParticles());
     contextData[&context] = data;
+    ReferenceConstraints& constraints = *(ReferenceConstraints*) reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData())->constraints;
+    if (constraints.settle != NULL) {
+        CpuSETTLE* parallelSettle = new CpuSETTLE(context.getSystem(), *(ReferenceSETTLEAlgorithm*) constraints.settle, data->threads);
+        delete constraints.settle;
+        constraints.settle = parallelSettle;
+    }
 }
 
 void CpuPlatform::contextDestroyed(ContextImpl& context) const {

platforms/cpu/src/CpuSETTLE.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) 2013 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 "CpuSETTLE.h"
+
+using namespace OpenMM;
+using namespace std;
+
+class CpuSETTLE::ApplyToPositionsTask : public ThreadPool::Task {
+public:
+    ApplyToPositionsTask(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses,
+            RealOpenMM tolerance, vector<ReferenceSETTLEAlgorithm*>& threadSettle) : atomCoordinates(atomCoordinates), atomCoordinatesP(atomCoordinatesP),
+            inverseMasses(inverseMasses), tolerance(tolerance), threadSettle(threadSettle) {
+    }
+    void execute(ThreadPool& threads, int threadIndex) {
+        threadSettle[threadIndex]->apply(atomCoordinates, atomCoordinatesP, inverseMasses, tolerance);
+    }
+    vector<OpenMM::RealVec>& atomCoordinates;
+    vector<OpenMM::RealVec>& atomCoordinatesP;
+    vector<RealOpenMM>& inverseMasses;
+    RealOpenMM tolerance;
+    vector<ReferenceSETTLEAlgorithm*>& threadSettle;
+};
+
+class CpuSETTLE::ApplyToVelocitiesTask : public ThreadPool::Task {
+public:
+    ApplyToVelocitiesTask(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses,
+            RealOpenMM tolerance, vector<ReferenceSETTLEAlgorithm*>& threadSettle) : atomCoordinates(atomCoordinates), velocities(velocities),
+            inverseMasses(inverseMasses), tolerance(tolerance), threadSettle(threadSettle) {
+    }
+    void execute(ThreadPool& threads, int threadIndex) {
+        threadSettle[threadIndex]->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
+    }
+    vector<OpenMM::RealVec>& atomCoordinates;
+    vector<OpenMM::RealVec>& velocities;
+    vector<RealOpenMM>& inverseMasses;
+    RealOpenMM tolerance;
+    vector<ReferenceSETTLEAlgorithm*>& threadSettle;
+};
+
+CpuSETTLE::CpuSETTLE(const System& system, const ReferenceSETTLEAlgorithm& settle, ThreadPool& threads) : threads(threads) {
+    int numThreads = threads.getNumThreads();
+    int numClusters = settle.getNumClusters();
+    vector<RealOpenMM> mass(system.getNumParticles());
+    for (int i = 0; i < system.getNumParticles(); i++)
+        mass[i] = system.getParticleMass(i);
+    for (int i = 0; i < numThreads; i++) {
+        int start = i*numClusters/numThreads;
+        int end = (i+1)*numClusters/numThreads;
+        if (start != end) {
+            int numThreadClusters = end-start;
+            vector<int> atom1(numThreadClusters), atom2(numThreadClusters), atom3(numThreadClusters);
+            vector<RealOpenMM> distance1(numThreadClusters), distance2(numThreadClusters);
+            for (int j = 0; j < numThreadClusters; j++)
+                settle.getClusterParameters(start+j, atom1[j], atom2[j], atom3[j], distance1[j], distance2[j]);
+            threadSettle.push_back(new ReferenceSETTLEAlgorithm(atom1, atom2, atom3, distance1, distance2, mass));
+        }
+    }
+}
+
+CpuSETTLE::~CpuSETTLE() {
+    for (int i = 0; i < (int) threadSettle.size(); i++)
+        delete threadSettle[i];
+}
+
+void CpuSETTLE::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
+    ApplyToPositionsTask task(atomCoordinates, atomCoordinatesP, inverseMasses, tolerance, threadSettle);
+    threads.execute(task);
+    threads.waitForThreads();
+}
+
+void CpuSETTLE::applyToVelocities(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
+    ApplyToVelocitiesTask task(atomCoordinates, velocities, inverseMasses, tolerance, threadSettle);
+    threads.execute(task);
+    threads.waitForThreads();
+}

platforms/cpu/tests/TestCpuSettle.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-2013 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CPU implementation of the SETTLE algorithm.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CpuPlatform.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testConstraints() {
+    const int numMolecules = 10;
+    const int numParticles = numMolecules*3;
+    const int numConstraints = numMolecules*3;
+    const double temp = 100.0;
+    CpuPlatform platform;
+    System system;
+    LangevinIntegrator integrator(temp, 2.0, 0.001);
+    integrator.setConstraintTolerance(1e-5);
+    NonbondedForce* forceField = new NonbondedForce();
+    for (int i = 0; i < numMolecules; ++i) {
+        system.addParticle(16.0);
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        forceField->addParticle(-0.82, 0.317, 0.65);
+        forceField->addParticle(0.41, 1.0, 0.0);
+        forceField->addParticle(0.41, 1.0, 0.0);
+        system.addConstraint(i*3, i*3+1, 0.1);
+        system.addConstraint(i*3, i*3+2, 0.1);
+        system.addConstraint(i*3+1, i*3+2, 0.163);
+    }
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numMolecules; ++i) {
+        positions[i*3] = Vec3((i%4)*0.4, (i/4)*0.4, 0);
+        positions[i*3+1] = positions[i*3]+Vec3(0.1, 0, 0);
+        positions[i*3+2] = positions[i*3]+Vec3(-0.03333, 0.09428, 0);
+        velocities[i*3] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+        velocities[i*3+1] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+        velocities[i*3+2] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+    }
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+
+    // Simulate it and see whether the constraints remain satisfied.
+
+    for (int i = 0; i < 1000; ++i) {
+        integrator.step(1);
+        State state = context.getState(State::Positions | State::Forces);
+        for (int j = 0; j < numConstraints; ++j) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(j, particle1, particle2, distance);
+            Vec3 p1 = state.getPositions()[particle1];
+            Vec3 p2 = state.getPositions()[particle2];
+            double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
+            ASSERT_EQUAL_TOL(distance, dist, 1e-5);
+        }
+    }
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        testConstraints();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/reference/include/ReferenceConstraints.h

@@ -33,8 +33,6 @@
  * -------------------------------------------------------------------------- */
 
 #include "ReferenceConstraintAlgorithm.h"
-#include "ReferenceCCMAAlgorithm.h"
-#include "ReferenceSETTLEAlgorithm.h"
 #include "openmm/System.h"
 
 namespace OpenMM {
@@ -47,7 +45,7 @@ namespace OpenMM {
 class OPENMM_EXPORT ReferenceConstraints : public ReferenceConstraintAlgorithm {
 public:
     ReferenceConstraints(const System& system);
-    ~ReferenceConstraints();
+    virtual ~ReferenceConstraints();
 
     /**
      * Apply the constraint algorithm.
@@ -68,9 +66,8 @@ public:
      * @param tolerance        the constraint tolerance
      */
     void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
-private:
-    ReferenceCCMAAlgorithm* ccma;
-    ReferenceSETTLEAlgorithm* settle;
+    ReferenceConstraintAlgorithm* ccma;
+    ReferenceConstraintAlgorithm* settle;
 };
 
 } // namespace OpenMM

platforms/reference/include/ReferenceSETTLEAlgorithm.h

@@ -46,6 +46,21 @@ public:
     ReferenceSETTLEAlgorithm(const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<int>& atom3,
             const std::vector<RealOpenMM>& distance1, const std::vector<RealOpenMM>& distance2, std::vector<RealOpenMM>& masses);
 
+    /**
+     * Get the number of clusters.
+     */
+    int getNumClusters() const;
+    /**
+     * Get the parameters describing one cluster.
+     * 
+     * @param index       the index of the cluster to get
+     * @param atom1       the index of the first atom in the cluster
+     * @param atom2       the index of the second atom in the cluster
+     * @param atom3       the index of the third atom in the cluster
+     * @param distance1   the distance between atoms 1 and 2
+     * @param distance2   the distance between atoms 2 and 3
+     */
+    void getClusterParameters(int index, int& atom1, int& atom2, int& atom3, RealOpenMM& distance1, RealOpenMM& distance2) const;
     /**
      * Apply the constraint algorithm.
      * 

platforms/reference/src/SimTKReference/ReferenceConstraints.cpp

@@ -30,6 +30,8 @@
  * -------------------------------------------------------------------------- */
 
 #include "ReferenceConstraints.h"
+#include "ReferenceCCMAAlgorithm.h"
+#include "ReferenceSETTLEAlgorithm.h"
 #include "openmm/HarmonicAngleForce.h"
 #include "openmm/OpenMMException.h"
 #include <map>

platforms/reference/src/SimTKReference/ReferenceSETTLEAlgorithm.cpp

@@ -39,6 +39,18 @@ ReferenceSETTLEAlgorithm::ReferenceSETTLEAlgorithm(const vector<int>& atom1, con
         atom1(atom1), atom2(atom2), atom3(atom3), distance1(distance1), distance2(distance2), masses(masses) {
 }
 
+int ReferenceSETTLEAlgorithm::getNumClusters() const {
+    return atom1.size();
+}
+
+void ReferenceSETTLEAlgorithm::getClusterParameters(int index, int& atom1, int& atom2, int& atom3, RealOpenMM& distance1, RealOpenMM& distance2) const {
+    atom1 = this->atom1[index];
+    atom2 = this->atom2[index];
+    atom3 = this->atom3[index];
+    distance1 = this->distance1[index];
+    distance2 = this->distance2[index];
+}
+
 void ReferenceSETTLEAlgorithm::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
     for (int index = 0; index < (int) atom1.size(); ++index) {
         RealVec apos0 = atomCoordinates[atom1[index]];