Implements concerted rmsd calculation (#3)

* Changes API to deal with groups of particles

* Updated serialization code

* Updated Python API

* Improved particle group handling

* Implemented concerted RMSD calculation proper

* Store positions of group particles only

Author: Charlles Abreu

.github/workflows/Linux.yml

@@ -14,7 +14,7 @@ on:
 jobs:
   unix:
     runs-on: ubuntu-20.04
-    name: Linux (Python ${{ matrix.python-version }}, OpenMM ${{ matrix.openmm-version }}, CUDA ${{ matrix.cuda-version }})
+    name: Linux (Python ${{ matrix.python-version }}, OpenMM ${{ matrix.openmm-version }}, GCC ${{ matrix.gcc-version }})
     strategy:
       fail-fast: false
       matrix:

customcppforces/include/ConcertedRMSDForce.h

@@ -11,18 +11,21 @@
  *  https://github.com/craabreu/customcppforces                               *
  * -------------------------------------------------------------------------- */
 
+#include "internal/windowsExportCustomCPPForces.h"
+
 #include "openmm/Force.h"
 #include "openmm/Vec3.h"
+#include "openmm/internal/AssertionUtilities.h"
 #include <vector>
-#include "internal/windowsExportCustomCPPForces.h"
 
 using namespace OpenMM;
+using namespace std;
 
 namespace CustomCPPForces {
 
 /**
- * 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
+ * This is a force whose energy equals a special type of 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.
@@ -38,43 +41,54 @@ class CUSTOM_CPP_FORCES_EXPORT ConcertedRMSDForce : public Force {
     /**
      * Create an ConcertedRMSDForce.
      *
-     * @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.
+     * @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 Concerted RMSD.
      */
-    explicit ConcertedRMSDForce(const std::vector<Vec3>& referencePositions,
-                       const std::vector<int>& particles=std::vector<int>());
+    explicit ConcertedRMSDForce(const vector<Vec3>& referencePositions);
     /**
      * Get the reference positions to compute the deviation from.
      */
-    const std::vector<Vec3>& getReferencePositions() const {
+    const vector<Vec3>& getReferencePositions() const {
         return referencePositions;
     }
     /**
      * Set the reference positions to compute the deviation from.
+     *
+     * @param positions    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 concerted RMSD.
+     */
+    void setReferencePositions(const vector<Vec3>& positions);
+    /**
+     * Add a group of particles to be included in the concerted RMSD calculation.
+     *
+     * @param particles    the indices of the particles to include
+     *
+     * @return the index of the group that was added
      */
-    void setReferencePositions(const std::vector<Vec3>& positions);
+    int addGroup(const vector<int>& particles);
     /**
-     * Get the indices of the particles to use when computing the RMSD.  If this
-     * is empty, all particles in the system will be used.
+     * Get the number of particle groups included in the concerted RMSD calculation.
      */
-    const std::vector<int>& getParticles() const {
-        return particles;
-    }
+    int getNumGroups() const;
+    /**
+     * Get the particles of a group included in the concerted RMSD calculation.
+     *
+     * @param index    the index of the group whose particles are to be retrieved
+     */
+    const vector<int>& getGroup(int index) const;
     /**
-     * Set the indices of the particles to use when computing the RMSD.  If this
-     * is empty, all particles in the system will be used.
+     * Set the particles of a group included in the concerted RMSD calculation.
      */
-    void setParticles(const std::vector<int>& particles);
+    void setGroup(int index, const 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
+     * Update the reference positions and particle groups in a Context to match those stored
+     * in this Force object.  This method provides an efficient way to update these parameters
      * in an existing Context without needing to reinitialize it.  Simply call setReferencePositions()
-     * and setParticles() to modify this object's parameters, then call updateParametersInContext()
+     * and setGroup() to modify this object's parameters, then call updateParametersInContext()
      * to copy them over to the Context.
      */
     void updateParametersInContext(Context& context);
@@ -90,8 +104,8 @@ class CUSTOM_CPP_FORCES_EXPORT ConcertedRMSDForce : public Force {
 protected:
     ForceImpl* createImpl() const;
 private:
-    std::vector<Vec3> referencePositions;
-    std::vector<int> particles;
+    vector<Vec3> referencePositions;
+    vector<vector<int>> groups;
 };
 
 } // namespace CustomCPPForces

customcppforces/include/internal/ConcertedRMSDForceImpl.h

@@ -36,9 +36,9 @@ class ConcertedRMSDForceImpl : public CustomCPPForceImpl {
     }
     void updateParametersInContext(ContextImpl& context);
 private:
+    void updateParameters(int systemSize);
     const ConcertedRMSDForce& owner;
-    int numParticles;
-    vector<int> particles;
+    vector<vector<int>> groups;
     vector<Vec3> referencePos;
 };
 

customcppforces/src/ConcertedRMSDForce.cpp

@@ -11,20 +11,37 @@
 #include "ConcertedRMSDForce.h"
 #include "internal/ConcertedRMSDForceImpl.h"
 
+#include "openmm/internal/AssertionUtilities.h"
+
 using namespace CustomCPPForces;
 using namespace OpenMM;
 using namespace std;
 
-ConcertedRMSDForce::ConcertedRMSDForce(const vector<Vec3>& referencePositions, const vector<int>& particles) :
-        referencePositions(referencePositions), particles(particles) {
+ConcertedRMSDForce::ConcertedRMSDForce(const vector<Vec3>& referencePositions) :
+        referencePositions(referencePositions) {
 }
 
 void ConcertedRMSDForce::setReferencePositions(const std::vector<Vec3>& positions) {
     referencePositions = positions;
 }
 
-void ConcertedRMSDForce::setParticles(const std::vector<int>& particles) {
-    this->particles = particles;
+int ConcertedRMSDForce::addGroup(const vector<int>& particles) {
+    groups.push_back(particles);
+    return groups.size()-1;
+}
+
+int ConcertedRMSDForce::getNumGroups() const {
+    return groups.size();
+}
+
+const vector<int>& ConcertedRMSDForce::getGroup(int index) const {
+    ASSERT_VALID_INDEX(index, groups);
+    return groups[index];
+}
+
+void ConcertedRMSDForce::setGroup(int index, const std::vector<int>& particles) {
+    ASSERT_VALID_INDEX(index, groups);
+    groups[index] = particles;
 }
 
 void ConcertedRMSDForce::updateParametersInContext(Context& context) {

customcppforces/src/ConcertedRMSDForceImpl.cpp

@@ -21,68 +21,85 @@ using namespace CustomCPPForces;
 using namespace OpenMM;
 using namespace std;
 
-void ConcertedRMSDForceImpl::initialize(ContextImpl& context) {
-    CustomCPPForceImpl::initialize(context);
 
+void ConcertedRMSDForceImpl::updateParameters(int systemSize) {
     // Check for errors in the specification of particles.
-    const System& system = context.getSystem();
-    int systemSize = system.getNumParticles();
     if (owner.getReferencePositions().size() != systemSize)
         throw OpenMMException(
-            "RMSDForce: Number of reference positions does not equal number of particles in the System"
+            "ConcertedRMSDForce: Number of reference positions does not equal number of particles in the System"
         );
 
-    particles = owner.getParticles();
-    numParticles = particles.size();
+    int numGroups = owner.getNumGroups();
+    if (numGroups == 0)
+        throw OpenMMException("ConcertedRMSDForce: No particle groups have been specified");
+
+    groups.resize(numGroups);
+    for (int i = 0; i < numGroups; i++) {
+        groups[i] = owner.getGroup(i);
+        if (groups[i].size() == 0) {
+            groups[i].resize(systemSize);
+            for (int j = 0; j < systemSize; j++)
+                groups[i][j] = j;
+        }
+    }
 
     set<int> distinctParticles;
-    for (int i : particles) {
-        if (i < 0 || i >= systemSize) {
-            stringstream msg;
-            msg << "ConcertedRMSDForce: Illegal particle index for RMSD: ";
-            msg << i;
-            throw OpenMMException(msg.str());
-        }
-        if (distinctParticles.find(i) != distinctParticles.end()) {
-            stringstream msg;
-            msg << "ConcertedRMSDForce: Duplicated particle index for RMSD: ";
-            msg << i;
-            throw OpenMMException(msg.str());
+    for (int k = 0; k < numGroups; k++)
+        for (int i : groups[k]) {
+            if (i < 0 || i >= systemSize) {
+                stringstream msg;
+                msg << "ConcertedRMSDForce: Illegal particle index " << i << " in group " << k;
+                throw OpenMMException(msg.str());
+            }
+            if (distinctParticles.find(i) != distinctParticles.end()) {
+                stringstream msg;
+                msg << "ConcertedRMSDForce: Duplicated particle index " << i << " in group " << k;
+                throw OpenMMException(msg.str());
+            }
+            distinctParticles.insert(i);
         }
-        distinctParticles.insert(i);
+
+    referencePos.resize(0);
+    const vector<Vec3>& positions = owner.getReferencePositions();
+    for (auto& group : groups) {
+        Vec3 center(0.0, 0.0, 0.0);
+        for (int i : group)
+            center += positions[i];
+        center /= group.size();
+        for (int i : group)
+            referencePos.push_back(positions[i] - center);
     }
+}
 
-    referencePos = owner.getReferencePositions();
-    Vec3 center(0.0, 0.0, 0.0);
-    for (int i : particles)
-        center += referencePos[i];
-    center /= numParticles;
-    for (Vec3& p : referencePos)
-        p -= center;
+void ConcertedRMSDForceImpl::initialize(ContextImpl& context) {
+    CustomCPPForceImpl::initialize(context);
+    updateParameters(context.getSystem().getNumParticles());
 }
 
 double ConcertedRMSDForceImpl::computeForce(ContextImpl& context, const vector<Vec3>& positions, vector<Vec3>& forces) {
     // 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.
 
-    Vec3 center(0.0, 0.0, 0.0);
-    for (int i : particles)
-        center += positions[i];
-    center /= numParticles;
+    int numParticles = referencePos.size();
     vector<Vec3> centeredPos(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        centeredPos[i] = positions[particles[i]] - center;
+    int index = 0;
+    for (auto& group : groups) {
+        Vec3 center(0.0, 0.0, 0.0);
+        for (int i : group)
+            center += positions[i];
+        center /= group.size();
+        for (int i : group)
+            centeredPos[index++] = positions[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] += centeredPos[k][i]*referencePos[index][j];
-            }
+            for (int index = 0; index < numParticles; index++)
+                R[i][j] += centeredPos[index][i]*referencePos[index][j];
 
     // Compute the F matrix.
 
@@ -118,10 +135,10 @@ double ConcertedRMSDForceImpl::computeForce(ContextImpl& context, const vector<V
     // Compute the RMSD.
 
     double sum = 0.0;
-    for (int i = 0; i < numParticles; i++) {
-        int index = particles[i];
-        sum += centeredPos[i].dot(centeredPos[i]) + referencePos[index].dot(referencePos[index]);
-    }
+    for (auto& group : groups)
+        for (int index = 0; index < group.size(); index++)
+            sum += centeredPos[index].dot(centeredPos[index]) + referencePos[index].dot(referencePos[index]);
+
     double msd = (sum - 2*values[3])/numParticles;
     if (msd < 1e-20) {
         // The particles are perfectly aligned, so all the forces should be zero.
@@ -143,30 +160,20 @@ double ConcertedRMSDForceImpl::computeForce(ContextImpl& context, const vector<V
 
     // 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]] = -(centeredPos[i] - rotatedRef) / (rmsd*numParticles);
+    index = 0;
+    for (auto& group : groups)
+        for (int i : group) {
+            const Vec3& p = referencePos[index];
+            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[i] = -(centeredPos[index] - rotatedRef) / (rmsd*groups[0].size());
+            index++;
     }
     return rmsd;
 }
 
 void ConcertedRMSDForceImpl::updateParametersInContext(ContextImpl& context) {
-    if (referencePos.size() != owner.getReferencePositions().size())
-        throw OpenMMException("updateParametersInContext: The number of reference positions has changed");
-    particles = owner.getParticles();
-    if (particles.size() == 0)
-        for (int i = 0; i < referencePos.size(); i++)
-            particles.push_back(i);
-    numParticles = particles.size();
-    referencePos = owner.getReferencePositions();
-    Vec3 center(0.0, 0.0, 0.0);
-    for (int i : particles)
-        center += referencePos[i];
-    center /= numParticles;
-    for (Vec3& p : referencePos)
-        p -= center;
+    updateParameters(context.getSystem().getNumParticles());
     context.systemChanged();
 }