Merge pull request #1982 from peastman/rmsd

Created RMSDForce

docs-source/usersguide/theory.rst

@@ -762,6 +762,26 @@ where :math:`m_i` and :math:`\mathbf{v}_i` are the mass and velocity of particle
 \ *i*\ .  It then subtracts :math:`\mathbf{v}_\text{CM}` from the velocity of every
 particle.
 
+RMSDForce
+*********
+
+RMSDForce computes the root-mean-squared deviation (RMSD) between the current
+particle positions :math:`\mathbf{x}_i` and a set of reference positions
+:math:`\mathbf{x}_i^\text{ref}`:
+
+.. math::
+   \text{RMSD} = \sqrt{\frac{\sum_{i} \| \mathbf{x}_i - \mathbf{x}_i^\text{ref} \|^2}{N}}
+
+Before computing this, the reference positions are first translated and rotated
+so as to minimize the RMSD.  The computed value is therefore :math:`argmin(\text{RMSD})`,
+where the :math:`argmin` is taken over all possible translations and rotations.
+
+This force is normally used with a CustomCVForce (see Section :ref:`customcvforce`).
+One rarely wants a force whose energy exactly equals the RMSD, but there are many
+situations where it is useful to have a restraining or biasing force that depends
+on the RMSD in some way.
+
+
 Custom Forces
 #############
 
@@ -1158,6 +1178,8 @@ 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
 *************
 

olla/include/openmm/kernels.h

@@ -57,6 +57,7 @@
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/PeriodicTorsionForce.h"
 #include "openmm/RBTorsionForce.h"
+#include "openmm/RMSDForce.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
 #include "openmm/VariableLangevinIntegrator.h"
@@ -981,6 +982,41 @@ public:
     virtual void copyState(ContextImpl& context, ContextImpl& innerContext) = 0;
 };
 
+/**
+ * This kernel is invoked by RMSDForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcRMSDForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcRMSDForce";
+    }
+    CalcRMSDForceKernel(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 RMSDForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const RMSDForce& 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;
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the RMSDForce to copy the parameters from
+     */
+    virtual void copyParametersToContext(ContextImpl& context, const RMSDForce& force) = 0;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

openmmapi/include/OpenMM.h

@@ -65,6 +65,7 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/PeriodicTorsionForce.h"
 #include "openmm/RBTorsionForce.h"
+#include "openmm/RMSDForce.h"
 #include "openmm/State.h"
 #include "openmm/System.h"
 #include "openmm/TabulatedFunction.h"

openmmapi/include/openmm/RMSDForce.h

+#ifndef OPENMM_RMSDFORCE_H_
+#define OPENMM_RMSDFORCE_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) 2018 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 "Vec3.h"
+#include <vector>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This is a force whose energy equals the root mean squared deviation (RMSD)
+ * between the current coordinates and a reference structure.  It is intended for
+ * use with CustomCVForce.  You will not normally want a force that exactly equals
+ * the RMSD, but there are many situations where it is useful to have a restraining
+ * or biasing force that depends on the RMSD in some way.
+ * 
+ * The force is computed by first aligning the particle positions to the reference
+ * structure, then computing the RMSD between the aligned positions and the reference.
+ * The computation can optionally be done based on only a subset of the particles
+ * in the system.
+ */
+
+class OPENMM_EXPORT RMSDForce : public Force {
+public:
+    /**
+     * Create an RMSDForce.
+     *
+     * @param referencePositions  the reference positions to compute the deviation
+     *                            from.  The length of this vector must equal the
+     *                            number of particles in the system, even if not
+     *                            all particles are used in computing the RMSD.
+     * @param particles           the indices of the particles to use when computing
+     *                            the RMSD.  If this is empty (the default), all
+     *                            particles in the system will be used.
+     */
+    explicit RMSDForce(const std::vector<Vec3>& referencePositions,
+                       const std::vector<int>& particles=std::vector<int>());
+    /**
+     * Get the reference positions to compute the deviation from.
+     */
+    const std::vector<Vec3>& getReferencePositions() const {
+        return referencePositions;
+    }
+    /**
+     * Set the reference positions to compute the deviation from.
+     */
+    void setReferencePositions(const std::vector<Vec3>& positions);
+    /**
+     * Get the indices of the particles to use when computing the RMSD.  If this
+     * is empty, all particles in the system will be used.
+     */
+    const std::vector<int>& getParticles() const {
+        return particles;
+    }
+    /**
+     * Set the indices of the particles to use when computing the RMSD.  If this
+     * is empty, all particles in the system will be used.
+     */
+    void setParticles(const std::vector<int>& particles);
+    /**
+     * Update the reference positions and particle indices in a Context to match those stored
+     * in this Force object.  This method provides an efficient method to update certain parameters
+     * in an existing Context without needing to reinitialize it.  Simply call setReferencePositions()
+     * and setParticles() to modify this object's parameters, then call updateParametersInContext()
+     * to copy them over to the Context.
+     */
+    void updateParametersInContext(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 {
+        return false;
+    }
+protected:
+    ForceImpl* createImpl() const;
+private:
+    std::vector<Vec3> referencePositions;
+    std::vector<int> particles;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_RMSDFORCE_H_*/

openmmapi/include/openmm/internal/RMSDForceImpl.h

+#ifndef OPENMM_RMSDFORCEIMPL_H_
+#define OPENMM_RMSDFORCEIMPL_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) 2018 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/RMSDForce.h"
+#include "openmm/Kernel.h"
+#include <utility>
+#include <map>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of RMSDForce.
+ */
+
+class RMSDForceImpl : public ForceImpl {
+public:
+    RMSDForceImpl(const RMSDForce& owner);
+    ~RMSDForceImpl();
+    void initialize(ContextImpl& context);
+    const RMSDForce& getOwner() const {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context, bool& forcesInvalid) {
+        // 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() {
+        return std::map<std::string, double>(); // This force doesn't define any parameters.
+    }
+    std::vector<std::string> getKernelNames();
+    void updateParametersInContext(ContextImpl& context);
+private:
+    const RMSDForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_RMSDFORCEIMPL_H_*/

openmmapi/src/RMSDForce.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) 2018 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/RMSDForce.h"
+#include "openmm/internal/RMSDForceImpl.h"
+
+using namespace OpenMM;
+using namespace std;
+
+RMSDForce::RMSDForce(const vector<Vec3>& referencePositions, const vector<int>& particles) :
+        referencePositions(referencePositions), particles(particles) {
+}
+
+void RMSDForce::setReferencePositions(const std::vector<Vec3>& positions) {
+    referencePositions = positions;
+}
+
+void RMSDForce::setParticles(const std::vector<int>& particles) {
+    this->particles = particles;
+}
+
+void RMSDForce::updateParametersInContext(Context& context) {
+    dynamic_cast<RMSDForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
+}
+
+ForceImpl* RMSDForce::createImpl() const {
+    return new RMSDForceImpl(*this);
+}

openmmapi/src/RMSDForceImpl.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) 2018 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/internal/ContextImpl.h"
+#include "openmm/internal/RMSDForceImpl.h"
+#include "openmm/kernels.h"
+#include <set>
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+RMSDForceImpl::RMSDForceImpl(const RMSDForce& owner) : owner(owner) {
+}
+
+RMSDForceImpl::~RMSDForceImpl() {
+}
+
+void RMSDForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcRMSDForceKernel::Name(), context);
+
+    // Check for errors in the specification of particles.
+    const System& system = context.getSystem();
+    int numParticles = system.getNumParticles();
+    if (owner.getReferencePositions().size() != numParticles)
+        throw OpenMMException("RMSDForce: Number of reference positions does not equal number of particles in the System");
+    set<int> particles;
+    for (int i : owner.getParticles()) {
+        if (i < 0 || i >= numParticles) {
+            stringstream msg;
+            msg << "RMSDForce: Illegal particle index for RMSD: ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+        if (particles.find(i) != particles.end()) {
+            stringstream msg;
+            msg << "RMSDForce: Duplicated particle index for RMSD: ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+        particles.insert(i);
+    }
+    kernel.getAs<CalcRMSDForceKernel>().initialize(context.getSystem(), owner);
+}
+
+double RMSDForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&(1<<owner.getForceGroup())) != 0)
+        return kernel.getAs<CalcRMSDForceKernel>().execute(context, includeForces, includeEnergy);
+    return 0.0;
+}
+
+vector<string> RMSDForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcRMSDForceKernel::Name());
+    return names;
+}
+
+void RMSDForceImpl::updateParametersInContext(ContextImpl& context) {
+    kernel.getAs<CalcRMSDForceKernel>().copyParametersToContext(context, owner);
+    context.systemChanged();
+}

platforms/cuda/include/CudaKernels.h

@@ -1278,6 +1278,59 @@ private:
     CUfunction copyStateKernel, copyForcesKernel, addForcesKernel;
 };
 
+/**
+ * This kernel is invoked by RMSDForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcRMSDForceKernel : public CalcRMSDForceKernel {
+public:
+    CudaCalcRMSDForceKernel(std::string name, const Platform& platform, CudaContext& cu) : CalcRMSDForceKernel(name, platform),
+            cu(cu), referencePos(NULL), particles(NULL), buffer(NULL) {
+    }
+    ~CudaCalcRMSDForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the RMSDForce this kernel will be used for
+     */
+    void initialize(const System& system, const RMSDForce& force);
+    /**
+     * Record the reference positions and particle indices.
+     */
+    void recordParameters(const RMSDForce& 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);
+    /**
+     * This is the internal implementation of execute(), templatized on whether we're
+     * using single or double precision.
+     */
+    template <class REAL>
+    double executeImpl(ContextImpl& context);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the RMSDForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const RMSDForce& force);
+private:
+    class ForceInfo;
+    CudaContext& cu;
+    ForceInfo* info;
+    double sumNormRef;
+    CudaArray* referencePos;
+    CudaArray* particles;
+    CudaArray* buffer;
+    CUfunction kernel1, kernel2;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/cuda/src/CudaKernelFactory.cpp

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

platforms/cuda/src/CudaKernels.cpp

@@ -51,6 +51,7 @@
 #include "ReferenceTabulatedFunction.h"
 #include "SimTKOpenMMRealType.h"
 #include "SimTKOpenMMUtilities.h"
+#include "jama_eig.h"
 #include <algorithm>
 #include <cmath>
 #include <set>
@@ -6784,6 +6785,200 @@ void CudaCalcCustomCVForceKernel::copyState(ContextImpl& context, ContextImpl& i
         innerContext.setParameter(param.first, context.getParameter(param.first));
 }
 
+class CudaCalcRMSDForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const RMSDForce& force) : force(force) {
+        updateParticles();
+    }
+    void updateParticles() {
+        particles.clear();
+        for (int i : force.getParticles())
+            particles.insert(i);
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        bool include1 = (particles.find(particle1) != particles.end());
+        bool include2 = (particles.find(particle2) != particles.end());
+        return (include1 == include2);
+    }
+private:
+    const RMSDForce& force;
+    set<int> particles;
+};
+
+CudaCalcRMSDForceKernel::~CudaCalcRMSDForceKernel() {
+    if (referencePos != NULL)
+        delete referencePos;
+    if (particles != NULL)
+        delete particles;
+    if (buffer != NULL)
+        delete buffer;
+}
+
+void CudaCalcRMSDForceKernel::initialize(const System& system, const RMSDForce& force) {
+    // Create data structures.
+    
+    bool useDouble = cu.getUseDoublePrecision();
+    int elementSize = (useDouble ? sizeof(double) : sizeof(float));
+    int numParticles = force.getParticles().size();
+    if (numParticles == 0)
+        numParticles = system.getNumParticles();
+    referencePos = new CudaArray(cu, system.getNumParticles(), 4*elementSize, "referencePos");
+    particles = CudaArray::create<int>(cu, numParticles, "particles");
+    buffer = new CudaArray(cu, 13, elementSize, "buffer");
+    recordParameters(force);
+    info = new ForceInfo(force);
+    cu.addForce(info);
+    
+    // Create the kernels.
+
+    CUmodule module = cu.createModule(CudaKernelSources::vectorOps+CudaKernelSources::rmsd);
+    kernel1 = cu.getKernel(module, "computeRMSDPart1");
+    kernel2 = cu.getKernel(module, "computeRMSDForces");
+}
+
+void CudaCalcRMSDForceKernel::recordParameters(const RMSDForce& force) {
+    // Record the parameters and center the reference positions.
+    
+    vector<int> particleVec = force.getParticles();
+    if (particleVec.size() == 0)
+        for (int i = 0; i < cu.getNumAtoms(); i++)
+            particleVec.push_back(i);
+    vector<Vec3> centeredPositions = force.getReferencePositions();
+    Vec3 center;
+    for (int i : particleVec)
+        center += centeredPositions[i];
+    center /= particleVec.size();
+    for (Vec3& p : centeredPositions)
+        p -= center;
+
+    // Upload them to the device.
+
+    particles->upload(particleVec);
+    if (cu.getUseDoublePrecision()) {
+        vector<double4> pos;
+        for (Vec3 p : centeredPositions)
+            pos.push_back(make_double4(p[0], p[1], p[2], 0));
+        referencePos->upload(pos);
+    }
+    else {
+        vector<float4> pos;
+        for (Vec3 p : centeredPositions)
+            pos.push_back(make_float4(p[0], p[1], p[2], 0));
+        referencePos->upload(pos);
+    }
+
+    // Record the sum of the norms of the reference positions.
+
+    sumNormRef = 0.0;
+    for (int i : particleVec) {
+        Vec3 p = centeredPositions[i];
+        sumNormRef += p.dot(p);
+    }
+}
+
+double CudaCalcRMSDForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (cu.getUseDoublePrecision())
+        return executeImpl<double>(context);
+    return executeImpl<float>(context);
+}
+
+template <class REAL>
+double CudaCalcRMSDForceKernel::executeImpl(ContextImpl& context) {
+    // Execute the first kernel.
+
+    int numParticles = particles->getSize();
+    int blockSize = 256;
+    void* args1[] = {&numParticles, &cu.getPosq().getDevicePointer(), &referencePos->getDevicePointer(),
+            &particles->getDevicePointer(), &buffer->getDevicePointer()};
+    cu.executeKernel(kernel1, args1, blockSize, blockSize, blockSize*sizeof(REAL));
+    
+    // Download the results, build the F matrix, and find the maximum eigenvalue
+    // and eigenvector.
+
+    vector<REAL> b;
+    buffer->download(b);
+    Array2D<double> F(4, 4);
+    F[0][0] =  b[0*3+0] + b[1*3+1] + b[2*3+2];
+    F[1][0] =  b[1*3+2] - b[2*3+1];
+    F[2][0] =  b[2*3+0] - b[0*3+2];
+    F[3][0] =  b[0*3+1] - b[1*3+0];
+    F[0][1] =  b[1*3+2] - b[2*3+1];
+    F[1][1] =  b[0*3+0] - b[1*3+1] - b[2*3+2];
+    F[2][1] =  b[0*3+1] + b[1*3+0];
+    F[3][1] =  b[0*3+2] + b[2*3+0];
+    F[0][2] =  b[2*3+0] - b[0*3+2];
+    F[1][2] =  b[0*3+1] + b[1*3+0];
+    F[2][2] = -b[0*3+0] + b[1*3+1] - b[2*3+2];
+    F[3][2] =  b[1*3+2] + b[2*3+1];
+    F[0][3] =  b[0*3+1] - b[1*3+0];
+    F[1][3] =  b[0*3+2] + b[2*3+0];
+    F[2][3] =  b[1*3+2] + b[2*3+1];
+    F[3][3] = -b[0*3+0] - b[1*3+1] + b[2*3+2];
+    JAMA::Eigenvalue<double> eigen(F);
+    Array1D<double> values;
+    eigen.getRealEigenvalues(values);
+    Array2D<double> vectors;
+    eigen.getV(vectors);
+
+    // Compute the RMSD.
+
+    double msd = (sumNormRef+b[9]-2*values[3])/numParticles;
+    if (msd < 1e-20) {
+        // The particles are perfectly aligned, so all the forces should be zero.
+        // Numerical error can lead to NaNs, so just return 0 now.
+        return 0.0;
+    }
+    double rmsd = sqrt(msd);
+    b[9] = rmsd;
+
+    // Compute the rotation matrix.
+
+    double q[] = {vectors[0][3], vectors[1][3], vectors[2][3], vectors[3][3]};
+    double q00 = q[0]*q[0], q01 = q[0]*q[1], q02 = q[0]*q[2], q03 = q[0]*q[3];
+    double q11 = q[1]*q[1], q12 = q[1]*q[2], q13 = q[1]*q[3];
+    double q22 = q[2]*q[2], q23 = q[2]*q[3];
+    double q33 = q[3]*q[3];
+    b[0] = q00+q11-q22-q33;
+    b[1] = 2*(q12-q03);
+    b[2] = 2*(q13+q02);
+    b[3] = 2*(q12+q03);
+    b[4] = q00-q11+q22-q33;
+    b[5] = 2*(q23-q01);
+    b[6] = 2*(q13-q02);
+    b[7] = 2*(q23+q01);
+    b[8] = q00-q11-q22+q33;
+
+    // Upload it to the device and invoke the kernel to apply forces.
+    
+    buffer->upload(b);
+    int paddedNumAtoms = cu.getPaddedNumAtoms();
+    void* args2[] = {&numParticles, &paddedNumAtoms, &cu.getPosq().getDevicePointer(), &referencePos->getDevicePointer(),
+            &particles->getDevicePointer(), &buffer->getDevicePointer(), &cu.getForce().getDevicePointer()};
+    cu.executeKernel(kernel2, args2, numParticles);
+    return rmsd;
+}
+
+void CudaCalcRMSDForceKernel::copyParametersToContext(ContextImpl& context, const RMSDForce& force) {
+    if (referencePos->getSize() != force.getReferencePositions().size())
+        throw OpenMMException("updateParametersInContext: The number of reference positions has changed");
+    int numParticles = force.getParticles().size();
+    if (numParticles == 0)
+        numParticles = context.getSystem().getNumParticles();
+    if (numParticles != particles->getSize()) {
+        // Recreate the particles array.
+        
+        delete particles;
+        particles = NULL;
+        particles = CudaArray::create<int>(cu, numParticles, "particles");
+    }
+    recordParameters(force);
+    
+    // Mark that the current reordering may be invalid.
+    
+    info->updateParticles();
+    cu.invalidateMolecules(info);
+}
+
 CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
 }
 

platforms/cuda/src/CudaPlatform.cpp

@@ -92,6 +92,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
+    registerKernelFactory(CalcRMSDForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);

platforms/cuda/src/kernels/rmsd.cu

+// This file contains kernels to compute the RMSD and its gradient using the algorithm described
+// in Coutsias et al, "Using quaternions to calculate RMSD" (doi: 10.1002/jcc.20110).
+
+/**
+ * Sum a value over all threads.
+ */
+__device__ real reduceValue(real value, volatile real* temp) {
+    const int thread = threadIdx.x;
+    temp[thread] = value;
+    __syncthreads();
+    for (uint step = 1; step < 32; step *= 2) {
+        if (thread+step < blockDim.x && thread%(2*step) == 0)
+            temp[thread] = temp[thread] + temp[thread+step];
+        SYNC_WARPS
+    }
+    for (uint step = 32; step < blockDim.x; step *= 2) {
+        if (thread+step < blockDim.x && thread%(2*step) == 0)
+            temp[thread] = temp[thread] + temp[thread+step];
+        __syncthreads();
+    }
+    return temp[0];
+}
+
+/**
+ * Perform the first step of computing the RMSD.  This is executed as a single work group.
+ */
+extern "C" __global__ void computeRMSDPart1(int numParticles, const real4* __restrict__ posq, const real4* __restrict__ referencePos,
+         const int* __restrict__ particles, real* buffer) {
+    extern __shared__ volatile real temp[];
+
+    // Compute the center of the particle positions.
+    
+    real3 center = make_real3(0);
+    for (int i = threadIdx.x; i < numParticles; i += blockDim.x)
+        center += trimTo3(posq[particles[i]]);
+    center.x = reduceValue(center.x, temp)/numParticles;
+    center.y = reduceValue(center.y, temp)/numParticles;
+    center.z = reduceValue(center.z, temp)/numParticles;
+    
+    // Compute the correlation matrix.
+    
+    real R[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
+    real sum = 0;
+    for (int i = threadIdx.x; i < numParticles; i += blockDim.x) {
+        int index = particles[i];
+        real3 pos = trimTo3(posq[index]) - center;
+        real3 refPos = trimTo3(referencePos[index]);
+        R[0][0] += pos.x*refPos.x;
+        R[0][1] += pos.x*refPos.y;
+        R[0][2] += pos.x*refPos.z;
+        R[1][0] += pos.y*refPos.x;
+        R[1][1] += pos.y*refPos.y;
+        R[1][2] += pos.y*refPos.z;
+        R[2][0] += pos.z*refPos.x;
+        R[2][1] += pos.z*refPos.y;
+        R[2][2] += pos.z*refPos.z;
+        sum += dot(pos, pos);
+    }
+    for (int i = 0; i < 3; i++)
+        for (int j = 0; j < 3; j++)
+            R[i][j] = reduceValue(R[i][j], temp);
+    sum = reduceValue(sum, temp);
+
+    // Copy everything into the output buffer to send back to the host.
+    
+    if (threadIdx.x == 0) {
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++)
+                buffer[3*i+j] = R[i][j];
+        buffer[9] = sum;
+        buffer[10] = center.x;
+        buffer[11] = center.y;
+        buffer[12] = center.z;
+    }
+}
+
+/**
+ * Apply forces based on the RMSD.
+ */
+extern "C" __global__ void computeRMSDForces(int numParticles, int paddedNumAtoms, const real4* __restrict__ posq, const real4* __restrict__ referencePos,
+         const int* __restrict__ particles, const real* buffer, unsigned long long* __restrict__ forceBuffers) {
+    real3 center = make_real3(buffer[10], buffer[11], buffer[12]);
+    real scale = 1 / (real) (buffer[9]*numParticles);
+    for (int i = blockDim.x*blockIdx.x+threadIdx.x; i < numParticles; i += blockDim.x*gridDim.x) {
+        int index = particles[i];
+        real3 pos = trimTo3(posq[index]) - center;
+        real3 refPos = trimTo3(referencePos[index]);
+        real3 rotatedRef = make_real3(buffer[0]*refPos.x + buffer[3]*refPos.y + buffer[6]*refPos.z,
+                                      buffer[1]*refPos.x + buffer[4]*refPos.y + buffer[7]*refPos.z,
+                                      buffer[2]*refPos.x + buffer[5]*refPos.y + buffer[8]*refPos.z);
+        real3 force = (rotatedRef-pos)*scale;
+        atomicAdd(&forceBuffers[index], static_cast<unsigned long long>((long long) (force.x*0x100000000)));
+        atomicAdd(&forceBuffers[index+paddedNumAtoms], static_cast<unsigned long long>((long long) (force.y*0x100000000)));
+        atomicAdd(&forceBuffers[index+2*paddedNumAtoms], static_cast<unsigned long long>((long long) (force.z*0x100000000)));
+    }
+}

platforms/cuda/tests/TestCudaRMSDForce.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) 2018 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 "TestRMSDForce.h"
+
+void runPlatformTests() {
+}

platforms/opencl/include/OpenCLKernels.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) 2008-2017 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1256,6 +1256,59 @@ private:
     cl::Kernel copyStateKernel, copyForcesKernel, addForcesKernel;
 };
 
+/**
+ * This kernel is invoked by RMSDForce to calculate the forces acting on the system and the energy of the system.
+ */
+class OpenCLCalcRMSDForceKernel : public CalcRMSDForceKernel {
+public:
+    OpenCLCalcRMSDForceKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcRMSDForceKernel(name, platform),
+            cl(cl), referencePos(NULL), particles(NULL), buffer(NULL) {
+    }
+    ~OpenCLCalcRMSDForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the RMSDForce this kernel will be used for
+     */
+    void initialize(const System& system, const RMSDForce& force);
+    /**
+     * Record the reference positions and particle indices.
+     */
+    void recordParameters(const RMSDForce& 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);
+    /**
+     * This is the internal implementation of execute(), templatized on whether we're
+     * using single or double precision.
+     */
+    template <class REAL>
+    double executeImpl(ContextImpl& context);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the RMSDForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const RMSDForce& force);
+private:
+    class ForceInfo;
+    OpenCLContext& cl;
+    ForceInfo* info;
+    double sumNormRef;
+    OpenCLArray* referencePos;
+    OpenCLArray* particles;
+    OpenCLArray* buffer;
+    cl::Kernel kernel1, kernel2;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

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-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -108,6 +108,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLCalcCustomCompoundBondForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcCustomCVForceKernel::Name())
         return new OpenCLCalcCustomCVForceKernel(name, platform, cl);
+    if (name == CalcRMSDForceKernel::Name())
+        return new OpenCLCalcRMSDForceKernel(name, platform, cl);
     if (name == CalcCustomManyParticleForceKernel::Name())
         return new OpenCLCalcCustomManyParticleForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcGayBerneForceKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -51,6 +51,7 @@
 #include "ReferenceTabulatedFunction.h"
 #include "SimTKOpenMMRealType.h"
 #include "SimTKOpenMMUtilities.h"
+#include "jama_eig.h"
 #include <algorithm>
 #include <cmath>
 #include <set>
@@ -7064,6 +7065,207 @@ void OpenCLCalcCustomCVForceKernel::copyState(ContextImpl& context, ContextImpl&
         innerContext.setParameter(param.first, context.getParameter(param.first));
 }
 
+class OpenCLCalcRMSDForceKernel::ForceInfo : public OpenCLForceInfo {
+public:
+    ForceInfo(const RMSDForce& force) : OpenCLForceInfo(0), force(force) {
+        updateParticles();
+    }
+    void updateParticles() {
+        particles.clear();
+        for (int i : force.getParticles())
+            particles.insert(i);
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        bool include1 = (particles.find(particle1) != particles.end());
+        bool include2 = (particles.find(particle2) != particles.end());
+        return (include1 == include2);
+    }
+private:
+    const RMSDForce& force;
+    set<int> particles;
+};
+
+OpenCLCalcRMSDForceKernel::~OpenCLCalcRMSDForceKernel() {
+    if (referencePos != NULL)
+        delete referencePos;
+    if (particles != NULL)
+        delete particles;
+    if (buffer != NULL)
+        delete buffer;
+}
+
+void OpenCLCalcRMSDForceKernel::initialize(const System& system, const RMSDForce& force) {
+    // Create data structures.
+    
+    bool useDouble = cl.getUseDoublePrecision();
+    int elementSize = (useDouble ? sizeof(cl_double) : sizeof(cl_float));
+    int numParticles = force.getParticles().size();
+    if (numParticles == 0)
+        numParticles = system.getNumParticles();
+    referencePos = new OpenCLArray(cl, system.getNumParticles(), 4*elementSize, "referencePos");
+    particles = OpenCLArray::create<cl_int>(cl, numParticles, "particles");
+    buffer = new OpenCLArray(cl, 13, elementSize, "buffer");
+    recordParameters(force);
+    info = new ForceInfo(force);
+    cl.addForce(info);
+    
+    // Create the kernels.
+
+    cl::Program program = cl.createProgram(OpenCLKernelSources::rmsd);
+    kernel1 = cl::Kernel(program, "computeRMSDPart1");
+    kernel2 = cl::Kernel(program, "computeRMSDForces");
+}
+
+void OpenCLCalcRMSDForceKernel::recordParameters(const RMSDForce& force) {
+    // Record the parameters and center the reference positions.
+    
+    vector<int> particleVec = force.getParticles();
+    if (particleVec.size() == 0)
+        for (int i = 0; i < cl.getNumAtoms(); i++)
+            particleVec.push_back(i);
+    vector<Vec3> centeredPositions = force.getReferencePositions();
+    Vec3 center;
+    for (int i : particleVec)
+        center += centeredPositions[i];
+    center /= particleVec.size();
+    for (Vec3& p : centeredPositions)
+        p -= center;
+
+    // Upload them to the device.
+
+    particles->upload(particleVec);
+    if (cl.getUseDoublePrecision()) {
+        vector<mm_double4> pos;
+        for (Vec3 p : centeredPositions)
+            pos.push_back(mm_double4(p[0], p[1], p[2], 0));
+        referencePos->upload(pos);
+    }
+    else {
+        vector<mm_float4> pos;
+        for (Vec3 p : centeredPositions)
+            pos.push_back(mm_float4(p[0], p[1], p[2], 0));
+        referencePos->upload(pos);
+    }
+
+    // Record the sum of the norms of the reference positions.
+
+    sumNormRef = 0.0;
+    for (int i : particleVec) {
+        Vec3 p = centeredPositions[i];
+        sumNormRef += p.dot(p);
+    }
+}
+
+double OpenCLCalcRMSDForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (cl.getUseDoublePrecision())
+        return executeImpl<double>(context);
+    return executeImpl<float>(context);
+}
+
+template <class REAL>
+double OpenCLCalcRMSDForceKernel::executeImpl(ContextImpl& context) {
+    // Execute the first kernel.
+
+    int numParticles = particles->getSize();
+    int blockSize = min(256, (int) kernel1.getWorkGroupInfo<CL_KERNEL_WORK_GROUP_SIZE>(cl.getDevice()));
+    kernel1.setArg<cl_int>(0, numParticles);
+    kernel1.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(2, referencePos->getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(3, particles->getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(4, buffer->getDeviceBuffer());
+    kernel1.setArg(5, blockSize*sizeof(REAL), NULL);
+    cl.executeKernel(kernel1, blockSize, blockSize);
+    
+    // Download the results, build the F matrix, and find the maximum eigenvalue
+    // and eigenvector.
+
+    vector<REAL> b;
+    buffer->download(b);
+    Array2D<double> F(4, 4);
+    F[0][0] =  b[0*3+0] + b[1*3+1] + b[2*3+2];
+    F[1][0] =  b[1*3+2] - b[2*3+1];
+    F[2][0] =  b[2*3+0] - b[0*3+2];
+    F[3][0] =  b[0*3+1] - b[1*3+0];
+    F[0][1] =  b[1*3+2] - b[2*3+1];
+    F[1][1] =  b[0*3+0] - b[1*3+1] - b[2*3+2];
+    F[2][1] =  b[0*3+1] + b[1*3+0];
+    F[3][1] =  b[0*3+2] + b[2*3+0];
+    F[0][2] =  b[2*3+0] - b[0*3+2];
+    F[1][2] =  b[0*3+1] + b[1*3+0];
+    F[2][2] = -b[0*3+0] + b[1*3+1] - b[2*3+2];
+    F[3][2] =  b[1*3+2] + b[2*3+1];
+    F[0][3] =  b[0*3+1] - b[1*3+0];
+    F[1][3] =  b[0*3+2] + b[2*3+0];
+    F[2][3] =  b[1*3+2] + b[2*3+1];
+    F[3][3] = -b[0*3+0] - b[1*3+1] + b[2*3+2];
+    JAMA::Eigenvalue<double> eigen(F);
+    Array1D<double> values;
+    eigen.getRealEigenvalues(values);
+    Array2D<double> vectors;
+    eigen.getV(vectors);
+
+    // Compute the RMSD.
+
+    double msd = (sumNormRef+b[9]-2*values[3])/numParticles;
+    if (msd < 1e-20) {
+        // The particles are perfectly aligned, so all the forces should be zero.
+        // Numerical error can lead to NaNs, so just return 0 now.
+        return 0.0;
+    }
+    double rmsd = sqrt(msd);
+    b[9] = rmsd;
+
+    // Compute the rotation matrix.
+
+    double q[] = {vectors[0][3], vectors[1][3], vectors[2][3], vectors[3][3]};
+    double q00 = q[0]*q[0], q01 = q[0]*q[1], q02 = q[0]*q[2], q03 = q[0]*q[3];
+    double q11 = q[1]*q[1], q12 = q[1]*q[2], q13 = q[1]*q[3];
+    double q22 = q[2]*q[2], q23 = q[2]*q[3];
+    double q33 = q[3]*q[3];
+    b[0] = q00+q11-q22-q33;
+    b[1] = 2*(q12-q03);
+    b[2] = 2*(q13+q02);
+    b[3] = 2*(q12+q03);
+    b[4] = q00-q11+q22-q33;
+    b[5] = 2*(q23-q01);
+    b[6] = 2*(q13-q02);
+    b[7] = 2*(q23+q01);
+    b[8] = q00-q11-q22+q33;
+
+    // Upload it to the device and invoke the kernel to apply forces.
+    
+    buffer->upload(b);
+    kernel2.setArg<cl_int>(0, numParticles);
+    kernel2.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
+    kernel2.setArg<cl::Buffer>(2, referencePos->getDeviceBuffer());
+    kernel2.setArg<cl::Buffer>(3, particles->getDeviceBuffer());
+    kernel2.setArg<cl::Buffer>(4, buffer->getDeviceBuffer());
+    kernel2.setArg<cl::Buffer>(5, cl.getForceBuffers().getDeviceBuffer());
+    cl.executeKernel(kernel2, numParticles);
+    return rmsd;
+}
+
+void OpenCLCalcRMSDForceKernel::copyParametersToContext(ContextImpl& context, const RMSDForce& force) {
+    if (referencePos->getSize() != force.getReferencePositions().size())
+        throw OpenMMException("updateParametersInContext: The number of reference positions has changed");
+    int numParticles = force.getParticles().size();
+    if (numParticles == 0)
+        numParticles = context.getSystem().getNumParticles();
+    if (numParticles != particles->getSize()) {
+        // Recreate the particles array.
+        
+        delete particles;
+        particles = NULL;
+        particles = OpenCLArray::create<cl_int>(cl, numParticles, "particles");
+    }
+    recordParameters(force);
+    
+    // Mark that the current reordering may be invalid.
+    
+    info->updateParticles();
+    cl.invalidateMolecules(info);
+}
+
 OpenCLIntegrateVerletStepKernel::~OpenCLIntegrateVerletStepKernel() {
 }
 

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-2017 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -83,6 +83,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
+    registerKernelFactory(CalcRMSDForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);

platforms/opencl/src/kernels/rmsd.cl

+// This file contains kernels to compute the RMSD and its gradient using the algorithm described
+// in Coutsias et al, "Using quaternions to calculate RMSD" (doi: 10.1002/jcc.20110).
+
+/**
+ * Sum a value over all threads.
+ */
+real reduceValue(real value, __local volatile real* temp) {
+    const int thread = get_local_id(0);
+    temp[thread] = value;
+    barrier(CLK_LOCAL_MEM_FENCE);
+    for (uint step = 1; step < 32; step *= 2) {
+        if (thread+step < get_local_size(0) && thread%(2*step) == 0)
+            temp[thread] = temp[thread] + temp[thread+step];
+        SYNC_WARPS
+    }
+    for (uint step = 32; step < get_local_size(0); step *= 2) {
+        if (thread+step < get_local_size(0) && thread%(2*step) == 0)
+            temp[thread] = temp[thread] + temp[thread+step];
+        barrier(CLK_LOCAL_MEM_FENCE);
+    }
+    return temp[0];
+}
+
+/**
+ * Perform the first step of computing the RMSD.  This is executed as a single work group.
+ */
+__kernel void computeRMSDPart1(int numParticles, __global const real4* restrict posq, __global const real4* restrict referencePos,
+        __global const int* restrict particles, __global real* buffer, __local volatile real* restrict temp) {
+    // Compute the center of the particle positions.
+    
+    real3 center = (real3) 0;
+    for (int i = get_local_id(0); i < numParticles; i += get_local_size(0))
+        center += posq[particles[i]].xyz;
+    center.x = reduceValue(center.x, temp)/numParticles;
+    center.y = reduceValue(center.y, temp)/numParticles;
+    center.z = reduceValue(center.z, temp)/numParticles;
+    
+    // Compute the correlation matrix.
+    
+    real R[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
+    real sum = 0;
+    for (int i = get_local_id(0); i < numParticles; i += get_local_size(0)) {
+        int index = particles[i];
+        real3 pos = posq[index].xyz - center;
+        real3 refPos = referencePos[index].xyz;
+        R[0][0] += pos.x*refPos.x;
+        R[0][1] += pos.x*refPos.y;
+        R[0][2] += pos.x*refPos.z;
+        R[1][0] += pos.y*refPos.x;
+        R[1][1] += pos.y*refPos.y;
+        R[1][2] += pos.y*refPos.z;
+        R[2][0] += pos.z*refPos.x;
+        R[2][1] += pos.z*refPos.y;
+        R[2][2] += pos.z*refPos.z;
+        sum += dot(pos, pos);
+    }
+    for (int i = 0; i < 3; i++)
+        for (int j = 0; j < 3; j++)
+            R[i][j] = reduceValue(R[i][j], temp);
+    sum = reduceValue(sum, temp);
+
+    // Copy everything into the output buffer to send back to the host.
+    
+    if (get_local_id(0) == 0) {
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++)
+                buffer[3*i+j] = R[i][j];
+        buffer[9] = sum;
+        buffer[10] = center.x;
+        buffer[11] = center.y;
+        buffer[12] = center.z;
+    }
+}
+
+/**
+ * Apply forces based on the RMSD.
+ */
+__kernel void computeRMSDForces(int numParticles, __global const real4* restrict posq, __global const real4* restrict referencePos,
+        __global const int* restrict particles, __global const real* buffer, __global real4* restrict forceBuffers) {
+    real3 center = (real3) (buffer[10], buffer[11], buffer[12]);
+    real scale = 1 / (real) (buffer[9]*numParticles);
+    for (int i = get_global_id(0); i < numParticles; i += get_global_size(0)) {
+        int index = particles[i];
+        real3 pos = posq[index].xyz - center;
+        real3 refPos = referencePos[index].xyz;
+        real3 rotatedRef = (real3) (buffer[0]*refPos.x + buffer[3]*refPos.y + buffer[6]*refPos.z,
+                                    buffer[1]*refPos.x + buffer[4]*refPos.y + buffer[7]*refPos.z,
+                                    buffer[2]*refPos.x + buffer[5]*refPos.y + buffer[8]*refPos.z);
+        forceBuffers[index].xyz -= (pos-rotatedRef)*scale;
+    }
+}

platforms/opencl/tests/TestOpenCLRMSDForce.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) 2018 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 "TestRMSDForce.h"
+
+void runPlatformTests() {
+}

platforms/reference/include/ReferenceKernels.h

@@ -1045,6 +1045,41 @@ private:
     std::vector<std::string> globalParameterNames, energyParamDerivNames;
 };
 
+/**
+ * This kernel is invoked by RMSDForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcRMSDForceKernel : public CalcRMSDForceKernel {
+public:
+    ReferenceCalcRMSDForceKernel(std::string name, const Platform& platform) : CalcRMSDForceKernel(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the RMSDForce this kernel will be used for
+     */
+    void initialize(const System& system, const RMSDForce& 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);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the RMSDForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const RMSDForce& force);
+private:
+    std::vector<Vec3> referencePos;
+    std::vector<int> particles;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */