Created RMSD force and reference implementation

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/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/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.
  */

platforms/reference/include/ReferenceRMSDForce.h

+/* Portions copyright (c) 2018 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * 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.
+ */
+
+#ifndef __ReferenceRMSDForce_H__
+#define __ReferenceRMSDForce_H__
+
+#include "openmm/RMSDForce.h"
+#include <vector>
+
+namespace OpenMM {
+
+class ReferenceRMSDForce {
+private:
+    std::vector<OpenMM::Vec3> referencePos;
+    std::vector<int> particles;
+
+public:
+    /**
+     * Constructor
+     */
+    ReferenceRMSDForce(std::vector<OpenMM::Vec3>& referencePos, std::vector<int>& particles);
+
+    /**
+     * Destructor
+     */
+    ~ReferenceRMSDForce();
+
+    /**
+     * Calculate the interaction.
+     * 
+     * @param atomCoordinates    atom coordinates
+     * @param forces             the forces are added to this
+     * @return the energy of the interaction
+     */
+   double calculateIxn(std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& forces) const;
+};
+
+} // namespace OpenMM
+
+#endif // __ReferenceRMSDForce_H__

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -80,6 +80,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomCompoundBondForceKernel(name, platform);
     if (name == CalcCustomCVForceKernel::Name())
         return new ReferenceCalcCustomCVForceKernel(name, platform);
+    if (name == CalcRMSDForceKernel::Name())
+        return new ReferenceCalcRMSDForceKernel(name, platform);
     if (name == CalcCustomManyParticleForceKernel::Name())
         return new ReferenceCalcCustomManyParticleForceKernel(name, platform);
     if (name == CalcGayBerneForceKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -57,6 +57,7 @@
 #include "ReferenceMonteCarloBarostat.h"
 #include "ReferenceProperDihedralBond.h"
 #include "ReferenceRbDihedralBond.h"
+#include "ReferenceRMSDForce.h"
 #include "ReferenceStochasticDynamics.h"
 #include "ReferenceTabulatedFunction.h"
 #include "ReferenceVariableStochasticDynamics.h"
@@ -2055,6 +2056,37 @@ void ReferenceCalcCustomCVForceKernel::copyState(ContextImpl& context, ContextIm
         innerContext.setParameter(param.first, context.getParameter(param.first));
 }
 
+void ReferenceCalcRMSDForceKernel::initialize(const System& system, const RMSDForce& force) {
+    particles = force.getParticles();
+    referencePos = force.getReferencePositions();
+    Vec3 center;
+    for (int i : particles)
+        center += referencePos[i];
+    center /= particles.size();
+    for (Vec3& p : referencePos)
+        p -= center;
+}
+
+double ReferenceCalcRMSDForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    vector<Vec3>& posData = extractPositions(context);
+    vector<Vec3>& forceData = extractForces(context);
+    ReferenceRMSDForce rmsd(referencePos, particles);
+    return rmsd.calculateIxn(posData, forceData);
+}
+
+void ReferenceCalcRMSDForceKernel::copyParametersToContext(ContextImpl& context, const RMSDForce& force) {
+    if (referencePos.size() != force.getReferencePositions().size())
+        throw OpenMMException("updateParametersInContext: The number of reference positions has changed");
+    particles = force.getParticles();
+    referencePos = force.getReferencePositions();
+    Vec3 center;
+    for (int i : particles)
+        center += referencePos[i];
+    center /= particles.size();
+    for (Vec3& p : referencePos)
+        p -= center;
+}
+
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferencePlatform.cpp

@@ -65,6 +65,7 @@ ReferencePlatform::ReferencePlatform() {
     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/reference/src/SimTKReference/ReferenceRMSDForce.cpp

+/* Portions copyright (c) 2018 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * 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 "ReferenceRMSDForce.h"
+#include "jama_eig.h"
+
+using namespace OpenMM;
+using namespace std;
+
+ReferenceRMSDForce::ReferenceRMSDForce(vector<Vec3>& referencePos, vector<int>& particles) :
+        referencePos(referencePos), particles(particles) {
+}
+
+ReferenceRMSDForce::~ReferenceRMSDForce() {
+}
+
+double ReferenceRMSDForce::calculateIxn(vector<Vec3>& atomCoordinates, vector<Vec3>& forces) const {
+    // Compute the RMSD and its gradient using the algorithm described in Coutsias et al,
+    // "Using quaternions to calculate RMSD" (doi: 10.1002/jcc.20110).  First subtract
+    // the centroid from the atom positions.  The reference positions have already been centered.
+    
+    int numParticles = particles.size();
+    Vec3 center;
+    for (int i : particles)
+        center += atomCoordinates[i];
+    center /= numParticles;
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = atomCoordinates[particles[i]]-center;
+    
+    // Compute the correlation matrix.
+    
+    double R[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
+    for (int i = 0; i < 3; i++)
+        for (int j = 0; j < 3; j++)
+            for (int k = 0; k < numParticles; k++) {
+                int index = particles[k];
+                R[i][j] += positions[k][i]*referencePos[index][j];
+            }
+
+    // Compute the F matrix.
+
+    Array2D<double> F(4, 4);
+    F[0][0] =  R[0][0] + R[1][1] + R[2][2];
+    F[1][0] =  R[1][2] - R[2][1];
+    F[2][0] =  R[2][0] - R[0][2];
+    F[3][0] =  R[0][1] - R[1][0];
+
+    F[0][1] =  R[1][2] - R[2][1];
+    F[1][1] =  R[0][0] - R[1][1] - R[2][2];
+    F[2][1] =  R[0][1] + R[1][0];
+    F[3][1] =  R[0][2] + R[2][0];
+
+    F[0][2] =  R[2][0] - R[0][2];
+    F[1][2] =  R[0][1] + R[1][0];
+    F[2][2] = -R[0][0] + R[1][1] - R[2][2];
+    F[3][2] =  R[1][2] + R[2][1];
+
+    F[0][3] =  R[0][1] - R[1][0];
+    F[1][3] =  R[0][2] + R[2][0];
+    F[2][3] =  R[1][2] + R[2][1];
+    F[3][3] = -R[0][0] - R[1][1] + R[2][2];
+
+    // Find the maximum eigenvalue and eigenvector.
+
+    JAMA::Eigenvalue<double> eigen(F);
+    Array1D<double> values;
+    eigen.getRealEigenvalues(values);
+    Array2D<double> vectors;
+    eigen.getV(vectors);
+
+    // Compute the RMSD.
+    
+    double sum = 0.0;
+    for (int i = 0; i < numParticles; i++) {
+        int index = particles[i];
+        sum += positions[i].dot(positions[i]) + referencePos[index].dot(referencePos[index]);
+    }
+    double rmsd = sqrt((sum-2*values[3])/numParticles);
+
+    // 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];
+    double U[3][3] = {{q00+q11-q22-q33, 2*(q12-q03), 2*(q13+q02)},
+                      {2*(q12+q03), q00-q11+q22-q33, 2*(q23-q01)},
+                      {2*(q13-q02), 2*(q23+q01), q00-q11-q22+q33}};
+
+    // Rotate the reference positions and compute the forces.
+    
+    for (int i = 0; i < numParticles; i++) {
+        const Vec3& p = referencePos[particles[i]];
+        Vec3 rotatedRef(U[0][0]*p[0] + U[1][0]*p[1] + U[2][0]*p[2],
+                        U[0][1]*p[0] + U[1][1]*p[1] + U[2][1]*p[2],
+                        U[0][2]*p[0] + U[1][2]*p[1] + U[2][2]*p[2]);
+        forces[particles[i]] -= (positions[i]-rotatedRef) / (rmsd*numParticles);
+    }
+    return rmsd;
+}

platforms/reference/tests/TestReferenceRMSDForce.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 "ReferenceTests.h"
+#include "TestRMSDForce.h"
+
+void runPlatformTests() {
+}

tests/TestRMSDForce.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 "openmm/internal/AssertionUtilities.h"
+#include "openmm/RMSDForce.h"
+#include "openmm/Context.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testRMSD() {
+    const int numParticles = 20;
+    System system;
+    vector<Vec3> referencePos(numParticles);
+    vector<Vec3> positions(numParticles);
+    vector<int> particles;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(1.0);
+        referencePos[i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*10;
+        positions[i] = referencePos[i] + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.2;
+        if (i%5 != 0)
+            particles.push_back(i);
+    }
+    RMSDForce* force = new RMSDForce(referencePos, particles);
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // Estimate the RMSD.  For simplicity we omit the orientation alignment, but they should
+    // already be almost perfectly aligned.
+    
+    Vec3 center1, center2;
+    for (int i : particles) {
+        center1 += referencePos[i];
+        center2 += positions[i];
+    }
+    center1 /= particles.size();
+    center2 /= particles.size();
+    double estimate = 0.0;
+    for (int i : particles) {
+        Vec3 delta = (referencePos[i]-center1) - (positions[i]-center2);
+        estimate += delta.dot(delta);
+    }
+    estimate = sqrt(estimate/particles.size());
+    
+    // Have the force compute the RMSD.  It should be very slightly less than
+    // what we calculated above (since that omitted the rotation).
+    
+    State state1 = context.getState(State::Energy);
+    double rmsd = state1.getPotentialEnergy();
+    ASSERT(rmsd <= estimate);
+    ASSERT(rmsd > 0.9*estimate);
+
+    // Translate and rotate all the particles.  This should have no effect on the RMSD.
+
+    double cs = cos(1.1), sn = sin(1.1);
+    for (int i = 0; i < numParticles; i++) {
+        Vec3 p = positions[i];
+        positions[i] = Vec3( cs*p[0] + sn*p[1] + 0.1,
+                            -sn*p[0] + cs*p[1] - 11.3,
+                            p[2] + 1.5);
+    }
+    context.setPositions(positions);
+    state1 = context.getState(State::Energy | State::Forces);
+    ASSERT_EQUAL_TOL(rmsd, state1.getPotentialEnergy(), 1e-5);
+
+    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+    const vector<Vec3>& forces = state1.getForces();
+    double norm = 0.0;
+    for (int i = 0; i < (int) forces.size(); ++i)
+        norm += forces[i].dot(forces[i]);
+    norm = std::sqrt(norm);
+    const double stepSize = 1e-3;
+    double step = 0.5*stepSize/norm;
+    vector<Vec3> positions2(numParticles), positions3(numParticles);
+    for (int i = 0; i < (int) positions.size(); ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = forces[i];
+        positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
+    }
+    context.setPositions(positions2);
+    State state2 = context.getState(State::Energy);
+    context.setPositions(positions3);
+    State state3 = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
+}
+
+void runPlatformTests();
+
+int main(int argc, char* argv[]) {
+    try {
+        initializeTests(argc, argv);
+        testRMSD();
+        runPlatformTests();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}