Continuing to implement CompoundIntegrator

openmmapi/include/OpenMM.h

@@ -36,6 +36,7 @@
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMAPTorsionForce.h"
 #include "openmm/CMMotionRemover.h"
+#include "openmm/CompoundIntegrator.h"
 #include "openmm/CustomBondForce.h"
 #include "openmm/CustomCentroidBondForce.h"
 #include "openmm/CustomCompoundBondForce.h"

openmmapi/include/openmm/CompoundIntegrator.h

@@ -73,6 +73,7 @@ public:
      * Create a CompoundIntegrator.
      */
     explicit CompoundIntegrator();
+    ~CompoundIntegrator();
     /**
      * Get the number of Integrators that have been added to this CompoundIntegrator.
      */
@@ -103,6 +104,30 @@ public:
      * @param index    the index of the Integrator to use
      */
     void setCurrentIntegrator(int index);
+    /**
+     * Get the size of each time step, in picoseconds.  This method calls getStepSize() on
+     * whichever Integrator has been set as current.
+     * 
+     * @return the step size, measured in ps
+     */
+    double getStepSize() const;
+    /**
+     * Set the size of each time step, in picoseconds.  This method calls setStepSize() on
+     * whichever Integrator has been set as current.
+     * 
+     * @param size    the step size, measured in ps
+     */
+    void setStepSize(double size);
+    /**
+     * Get the distance tolerance within which constraints are maintained, as a fraction of the constrained distance.
+     * This method calls getConstraintTolerance() on whichever Integrator has been set as current.
+     */
+    double getConstraintTolerance() const;
+    /**
+     * Set the distance tolerance within which constraints are maintained, as a fraction of the constrained distance.
+     * This method calls setConstraintTolerance() on whichever Integrator has been set as current.
+     */
+    void setConstraintTolerance(double tol);
    /**
      * Advance a simulation through time by taking a series of time steps.  This method
      * calls step() on whichever Integrator has been set as current.

openmmapi/include/openmm/Integrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -61,30 +61,22 @@ public:
      *
      * @return the step size, measured in ps
      */
-    double getStepSize() const {
-        return stepSize;
-    }
+    virtual double getStepSize() const;
     /**
      * Set the size of each time step, in picoseconds.  If this integrator uses variable time steps,
      * the effect of calling this method is undefined, and it may simply be ignored.
      *
      * @param size    the step size, measured in ps
      */
-    void setStepSize(double size) {
-        stepSize = size;
-    }
+    virtual void setStepSize(double size);
     /**
      * Get the distance tolerance within which constraints are maintained, as a fraction of the constrained distance.
      */
-    double getConstraintTolerance() const {
-        return constraintTol;
-    }
+    virtual double getConstraintTolerance() const;
     /**
      * Set the distance tolerance within which constraints are maintained, as a fraction of the constrained distance.
      */
-    void setConstraintTolerance(double tol) {
-        constraintTol = tol;
-    }
+    virtual void setConstraintTolerance(double tol);
     /**
      * Advance a simulation through time by taking a series of time steps.
      * 
@@ -94,6 +86,7 @@ public:
 protected:
     friend class Context;
     friend class ContextImpl;
+    friend class CompoundIntegrator;
     ContextImpl* context;
     Context* owner;
     /**

openmmapi/src/CompoundIntegrator.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) 2015 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/CompoundIntegrator.h"
+#include "openmm/Context.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/kernels.h"
+#include <string>
+
+using namespace OpenMM;
+using namespace std;
+
+CompoundIntegrator::CompoundIntegrator() : currentIntegrator(0) {
+}
+
+CompoundIntegrator::~CompoundIntegrator() {
+    for (int i = 0; i < integrators.size(); i++)
+        delete integrators[i];
+}
+
+int CompoundIntegrator::getNumIntegrators() const {
+    return integrators.size();
+}
+
+int CompoundIntegrator::addIntegrator(Integrator* integrator) {
+    if (owner != NULL)
+        throw OpenMMException("addIntegrator() cannot be called after a CustomIntegrator is already bound to a context");
+    integrators.push_back(integrator);
+    return integrators.size()-1;
+}
+
+Integrator& CompoundIntegrator::getIntegrator(int index) {
+    ASSERT_VALID_INDEX(index, integrators);
+    return *integrators[index];
+}
+
+int CompoundIntegrator::getCurrentIntegrator() const {
+    return currentIntegrator;
+}
+
+void CompoundIntegrator::setCurrentIntegrator(int index) {
+    if (index < 0 || index >= integrators.size())
+        throw OpenMMException("Illegal index for setCurrentIntegrator()");
+    currentIntegrator = index;
+}
+
+double CompoundIntegrator::getStepSize() const {
+    return integrators[currentIntegrator]->getStepSize();
+}
+
+void CompoundIntegrator::setStepSize(double size) {
+    integrators[currentIntegrator]->setStepSize(size);
+}
+
+double CompoundIntegrator::getConstraintTolerance() const {
+    return integrators[currentIntegrator]->getConstraintTolerance();
+}
+
+void CompoundIntegrator::setConstraintTolerance(double tol) {
+    integrators[currentIntegrator]->setConstraintTolerance(tol);
+}
+
+void CompoundIntegrator::step(int steps) {
+    integrators[currentIntegrator]->step(steps);
+}
+
+void CompoundIntegrator::initialize(ContextImpl& context) {
+    for (int i = 0; i < integrators.size(); i++)
+        integrators[i]->initialize(context);
+}
+
+void CompoundIntegrator::cleanup() {
+    for (int i = 0; i < integrators.size(); i++)
+        integrators[i]->cleanup();
+}
+
+vector<string> CompoundIntegrator::getKernelNames() {
+    vector<string> kernels;
+    for (int i = 0; i < integrators.size(); i++) {
+        vector<string> integratorKernels = integrators[i]->getKernelNames();
+        kernels.insert(kernels.end(), integratorKernels.begin(), integratorKernels.end());
+    }
+    return kernels;
+}
+
+double CompoundIntegrator::computeKineticEnergy() {
+    return integrators[currentIntegrator]->computeKineticEnergy();
+}

openmmapi/src/Integrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Portions copyright (c) 2013-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -47,3 +47,19 @@ Integrator::~Integrator() {
         context->integratorDeleted();
     }
 }
+
+double Integrator::getStepSize() const {
+    return stepSize;
+}
+
+void Integrator::setStepSize(double size) {
+    stepSize = size;
+}
+
+double Integrator::getConstraintTolerance() const {
+    return constraintTol;
+}
+
+void Integrator::setConstraintTolerance(double tol) {
+    constraintTol = tol;
+}

platforms/reference/tests/TestReferenceCompoundIntegrator.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) 2015 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 "TestCompoundIntegrator.h"
+
+void runPlatformTests() {
+}

tests/TestCompoundIntegrator.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) 2015 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/BrownianIntegrator.h"
+#include "openmm/CompoundIntegrator.h"
+#include "openmm/Context.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testChangingIntegrator() {
+    System system;
+    system.addParticle(2.0);
+    system.addParticle(2.0);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    bonds->addBond(0, 1, 1.5, 1);
+    system.addForce(bonds);
+    CompoundIntegrator integrator;
+    integrator.addIntegrator(new VerletIntegrator(0.01));
+    integrator.addIntegrator(new LangevinIntegrator(300.0, 10.0, 0.01));
+    integrator.addIntegrator(new BrownianIntegrator(300.0, 10.0, 0.01));
+    Context context(system, integrator, platform);
+    ASSERT_EQUAL(0, integrator.getCurrentIntegrator());
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    for (int iteration = 0; iteration < 2; ++iteration) {
+        context.setPositions(positions);
+
+        // First integrate with the Verlet integrator and compare it to the analytical solution.
+
+        const double freq = 1.0;
+        State state = context.getState(State::Energy);
+        const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
+        for (int i = 0; i < 100; ++i) {
+            state = context.getState(State::Positions | State::Velocities | State::Energy);
+            double time = state.getTime();
+            double expectedDist = 1.5+0.5*std::cos(freq*time);
+            ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
+            ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
+            double expectedSpeed = -0.5*freq*std::sin(freq*time);
+            ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
+            ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
+            double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
+            integrator.step(1);
+        }
+        ASSERT_EQUAL_TOL(1.0, context.getState(0).getTime(), 1e-5);
+
+        // Switch to the Langevin integrator and make sure that it heats up.
+
+        integrator.setCurrentIntegrator(1);
+        integrator.step(100);
+        double ke = 0.0;
+        for (int i = 0; i < 1000; ++i) {
+            integrator.step(10);
+            state = context.getState(State::Energy);
+            ke += state.getKineticEnergy();
+        }
+        double expectedKE = 0.5*2*3*BOLTZ*300.0;
+        ASSERT_USUALLY_EQUAL_TOL(expectedKE, ke/1000, 0.1);
+        
+        // Now reinitialize the context and repeat all of these tests to make sure that works correctly.
+        
+        context.reinitialize();
+        integrator.setCurrentIntegrator(0);
+    }
+}
+
+void testChangingParameters() {
+    System system;
+    system.addParticle(1.0);
+    CompoundIntegrator integrator;
+    integrator.addIntegrator(new VerletIntegrator(0.01));
+    integrator.addIntegrator(new LangevinIntegrator(300.0, 10.0, 0.02));
+    integrator.addIntegrator(new BrownianIntegrator(300.0, 10.0, 0.03));
+    
+    // Try getting and setting the step size for different component integrators.
+    
+    for (int i = 0; i < 3; i++) {
+        integrator.setCurrentIntegrator(i);
+        ASSERT_EQUAL_TOL(0.01*(i+1), integrator.getStepSize(), 1e-7);
+    }
+    for (int i = 0; i < 3; i++) {
+        integrator.setCurrentIntegrator(i);
+        integrator.setStepSize(0.02*(i+1));
+        ASSERT_EQUAL_TOL(0.02*(i+1), integrator.getStepSize(), 1e-7);
+    }
+    for (int i = 0; i < 3; i++) {
+        integrator.setCurrentIntegrator(i);
+        ASSERT_EQUAL_TOL(0.02*(i+1), integrator.getStepSize(), 1e-7);
+    }
+    
+    // Try getting and setting the constraint tolerance for different component integrators.
+    
+    for (int i = 0; i < 3; i++) {
+        integrator.setCurrentIntegrator(i);
+        ASSERT_EQUAL_TOL(1e-5, integrator.getConstraintTolerance(), 1e-7);
+    }
+    for (int i = 0; i < 3; i++) {
+        integrator.setCurrentIntegrator(i);
+        integrator.setConstraintTolerance(1e-4*(i+1));
+        ASSERT_EQUAL_TOL(1e-4*(i+1), integrator.getConstraintTolerance(), 1e-7);
+    }
+    for (int i = 0; i < 3; i++) {
+        integrator.setCurrentIntegrator(i);
+        ASSERT_EQUAL_TOL(1e-4*(i+1), integrator.getConstraintTolerance(), 1e-7);
+    }
+}
+
+void runPlatformTests();
+
+int main(int argc, char* argv[]) {
+    try {
+        initializeTests(argc, argv);
+        testChangingIntegrator();
+        testChangingParameters();
+        runPlatformTests();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}