Skip to content

Commit

Permalink
Implements concerted rmsd calculation (#3)
Browse files Browse the repository at this point in the history 
* 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
  • Loading branch information
@craabreu
craabreu authored Jan 24, 2024
1 parent 70379c6 commit ed4ec69
Show file tree
Hide file tree
Showing 9 changed files with 219 additions and 137 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/Linux.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
68 changes: 41 additions & 27 deletions customcppforces/include/ConcertedRMSDForce.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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);
Expand All @@ -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
Expand Down
4 changes: 2 additions & 2 deletions customcppforces/include/internal/ConcertedRMSDForceImpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};

Expand Down
25 changes: 21 additions & 4 deletions customcppforces/src/ConcertedRMSDForce.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
129 changes: 68 additions & 61 deletions customcppforces/src/ConcertedRMSDForceImpl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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.
Expand All @@ -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();
}
Loading

0 comments on commit ed4ec69

Please # to comment.