Implemented center of mass motion removal for OpenMM

olla/include/kernels.h

@@ -274,6 +274,30 @@ public:
     virtual double execute(const Stream& velocities) = 0;
 };
 
+/**
+ * This kernel is invoked to remove center of mass motion from the system.
+ */
+class RemoveCMMotionKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "RemoveCMMotion";
+    }
+    RemoveCMMotionKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel, setting up the atomic masses.
+     * 
+     * @param masses the mass of each atom
+     */
+    virtual void initialize(const std::vector<double>& masses) = 0;
+    /**
+     * Execute the kernel.
+     * 
+     * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
+     */
+    virtual void execute(Stream& velocities) = 0;
+};
+
 } // namespace OpenMM
 
 #endif /*OPENMM_KERNELS_H_*/

openmmapi/include/CMMotionRemover.h

+#ifndef OPENMM_CMMOTIONREMOVER_H_
+#define OPENMM_CMMOTIONREMOVER_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) 2008 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 <string>
+
+namespace OpenMM {
+
+/**
+ * This class prevents the center of mass of a System from drifting.  At each time step, it calculates the
+ * center of mass momentum, then adjusts the individual atom velocities to make it zero.
+ */
+
+class CMMotionRemover : public Force {
+public:
+    /**
+     * Create a CMMotionRemover.
+     */
+    CMMotionRemover();
+protected:
+    ForceImpl* createImpl();
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CMMOTIONREMOVER_H_*/

openmmapi/include/internal/CMMotionRemoverImpl.h

+#ifndef OPENMM_CMMOTIONREMOVERIMPL_H_
+#define OPENMM_CMMOTIONREMOVERIMPL_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) 2008 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 "CMMotionRemover.h"
+#include "Kernel.h"
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of CMMotionRemover.
+ */
+
+class CMMotionRemoverImpl : public ForceImpl {
+public:
+    CMMotionRemoverImpl(CMMotionRemover& owner);
+    void initialize(OpenMMContextImpl& context);
+    CMMotionRemover& getOwner() {
+        return owner;
+    }
+    void updateContextState(OpenMMContextImpl& context);
+    void calcForces(OpenMMContextImpl& context, Stream& forces) {
+        // This force doesn't apply forces to atoms.
+    }
+    double calcEnergy(OpenMMContextImpl& context) {
+        return 0.0; // This force doesn't contribute to the potential energy.
+    }
+    std::map<std::string, double> getDefaultParameters() {
+        return std::map<std::string, double>(); // This force doesn't define any parameters.
+    }
+    std::vector<std::string> getKernelNames();
+private:
+    CMMotionRemover& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CMMOTIONREMOVERIMPL_H_*/

openmmapi/src/CMMotionRemover.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 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 "CMMotionRemover.h"
+#include "internal/CMMotionRemoverImpl.h"
+
+using namespace OpenMM;
+
+CMMotionRemover::CMMotionRemover() {
+}
+
+ForceImpl* CMMotionRemover::createImpl() {
+    return new CMMotionRemoverImpl(*this);
+}

openmmapi/src/CMMotionRemoverImpl.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 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 "internal/CMMotionRemoverImpl.h"
+#include "internal/OpenMMContextImpl.h"
+#include "Integrator.h"
+#include "System.h"
+#include "kernels.h"
+#include <vector>
+
+using namespace OpenMM;
+using std::vector;
+
+CMMotionRemoverImpl::CMMotionRemoverImpl(CMMotionRemover& owner) : owner(owner) {
+}
+
+void CMMotionRemoverImpl::initialize(OpenMMContextImpl& context) {
+    kernel = context.getPlatform().createKernel(RemoveCMMotionKernel::Name());
+    const System& system = context.getSystem();
+    vector<double> masses(system.getNumAtoms());
+    for (int i = 0; i < system.getNumAtoms(); ++i)
+        masses[i] = system.getAtomMass(i);
+    dynamic_cast<RemoveCMMotionKernel&>(kernel.getImpl()).initialize(masses);
+}
+
+void CMMotionRemoverImpl::updateContextState(OpenMMContextImpl& context) {
+    dynamic_cast<RemoveCMMotionKernel&>(kernel.getImpl()).execute(context.getVelocities());
+}
+
+std::vector<std::string> CMMotionRemoverImpl::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(RemoveCMMotionKernel::Name());
+    return names;
+}

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -49,4 +49,6 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceApplyAndersenThermostatKernel(name, platform);
     if (name == CalcKineticEnergyKernel::Name())
         return new ReferenceCalcKineticEnergyKernel(name, platform);
+    if (name == RemoveCMMotionKernel::Name())
+        return new ReferenceRemoveCMMotionKernel(name, platform);
 }

platforms/reference/src/ReferenceKernels.cpp

@@ -432,3 +432,33 @@ double ReferenceCalcKineticEnergyKernel::execute(const Stream& velocities) {
         energy += masses[i]*(velData[i][0]*velData[i][0]+velData[i][1]*velData[i][1]+velData[i][2]*velData[i][2]);
     return 0.5*energy;
 }
+
+void ReferenceRemoveCMMotionKernel::initialize(const vector<double>& masses) {
+    this->masses.resize(masses.size());
+    for (int i = 0; i < masses.size(); ++i)
+        this->masses[i] = masses[i];
+}
+
+void ReferenceRemoveCMMotionKernel::execute(Stream& velocities) {
+    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) velocities.getImpl()).getData();
+    
+    // Calculate the center of mass momentum.
+    
+    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
+    for (int i = 0; i < masses.size(); ++i) {
+        momentum[0] += masses[i]*velData[i][0];
+        momentum[1] += masses[i]*velData[i][1];
+        momentum[2] += masses[i]*velData[i][2];
+    }
+    
+    // Adjust the atom velocities.
+    
+    momentum[0] /= masses.size();
+    momentum[1] /= masses.size();
+    momentum[2] /= masses.size();
+    for (int i = 0; i < masses.size(); ++i) {
+        velData[i][0] -= momentum[0]/masses[i];
+        velData[i][1] -= momentum[1]/masses[i];
+        velData[i][2] -= momentum[2]/masses[i];
+    }
+}

platforms/reference/src/ReferenceKernels.h

@@ -302,6 +302,29 @@ private:
     std::vector<double> masses;
 };
 
+/**
+ * This kernel is invoked to remove center of mass motion from the system.
+ */
+class ReferenceRemoveCMMotionKernel : public RemoveCMMotionKernel {
+public:
+    ReferenceRemoveCMMotionKernel(std::string name, const Platform& platform) : RemoveCMMotionKernel(name, platform) {
+    }
+    /**
+     * Initialize the kernel, setting up the atomic masses.
+     * 
+     * @param masses the mass of each atom
+     */
+    void initialize(const std::vector<double>& masses);
+    /**
+     * Execute the kernel.
+     * 
+     * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
+     */
+    void execute(Stream& velocities);
+private:
+    std::vector<double> masses;
+};
+
 } // namespace OpenMM
 
 #endif /*OPENMM_REFERENCEKERNELS_H_*/

platforms/reference/src/ReferencePlatform.cpp

@@ -53,6 +53,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
     registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);
+    registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
 }
 
 bool ReferencePlatform::supportsDoublePrecision() const {

platforms/reference/tests/TestReferenceCMMotionRemover.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the reference implementation of AndersenThermostat.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "CMMotionRemover.h"
+#include "OpenMMContext.h"
+#include "ReferencePlatform.h"
+#include "StandardMMForceField.h"
+#include "System.h"
+#include "VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+Vec3 calcCM(const vector<Vec3>& values, System& system) {
+    Vec3 cm;
+    for (int j = 0; j < system.getNumAtoms(); ++j) {
+        cm[0] += values[j][0]*system.getAtomMass(j);
+        cm[1] += values[j][1]*system.getAtomMass(j);
+        cm[2] += values[j][2]*system.getAtomMass(j);
+    }
+    return cm;
+}
+
+void testMotionRemoval() {
+    const int numAtoms = 8;
+    const double temp = 100.0;
+    const double collisionFreq = 10.0;
+    ReferencePlatform platform;
+    System system(numAtoms, 0);
+    VerletIntegrator integrator(0.01);
+    StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 1, 0, 0, 0);
+    for (int i = 0; i < numAtoms; ++i) {
+        system.setAtomMass(i, i+1);
+        forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    forceField->setBondParameters(0, 2, 3, 2.0, 0.5);
+    system.addForce(forceField);
+    CMMotionRemover* remover = new CMMotionRemover();
+    system.addForce(remover);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(numAtoms);
+    vector<Vec3> velocities(numAtoms);
+    init_gen_rand(0);
+    for (int i = 0; i < numAtoms; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
+    }
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    
+    // Now run it for a while and see if the center of mass remains fixed.
+    
+    Vec3 cmPos = calcCM(context.getState(State::Positions).getPositions(), system);
+    for (int i = 0; i < 1000; ++i) {
+        integrator.step(1);
+        State state = context.getState(State::Positions | State::Velocities);
+        Vec3 pos = calcCM(state.getPositions(), system);
+        ASSERT_EQUAL_VEC(cmPos, pos, 1e-2);
+        Vec3 vel = calcCM(state.getVelocities(), system);
+        ASSERT_EQUAL_VEC(Vec3(0, 0, 0), vel, 1e-2);
+    }
+}
+
+int main() {
+    try {
+        testMotionRemoval();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}