Finished implementing CompoundIntegrator

serialization/tests/TestSerializeIntegrator.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) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Yutong Zhao                                        *
  * Contributors:                                                              *
  *                                                                            *
@@ -32,6 +32,7 @@
 #include "openmm/internal/AssertionUtilities.h"
 
 #include "openmm/BrownianIntegrator.h"
+#include "openmm/CompoundIntegrator.h"
 #include "openmm/CustomIntegrator.h"
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/VariableLangevinIntegrator.h"
@@ -185,6 +186,29 @@ void testSerializeCustomIntegrator() {
     delete intg2;
 }
 
+void testSerializeCompoundIntegrator() {
+    CompoundIntegrator integ;
+    integ.addIntegrator(new LangevinIntegrator(372.4, 1.234, 0.0018));
+    integ.addIntegrator(new VerletIntegrator(0.002));
+    integ.setCurrentIntegrator(1);
+    stringstream ss;
+    XmlSerializer::serialize<Integrator>(&integ, "CompoundIntegrator", ss);
+    CompoundIntegrator *integ2 = dynamic_cast<CompoundIntegrator*>(XmlSerializer::deserialize<Integrator>(ss));
+    ASSERT_EQUAL(integ.getCurrentIntegrator(), integ2->getCurrentIntegrator());
+    LangevinIntegrator& langevin1 = dynamic_cast<LangevinIntegrator&>(integ.getIntegrator(0));
+    LangevinIntegrator& langevin2 = dynamic_cast<LangevinIntegrator&>(integ2->getIntegrator(0));
+    ASSERT_EQUAL(langevin1.getConstraintTolerance(), langevin2.getConstraintTolerance());
+    ASSERT_EQUAL(langevin1.getStepSize(), langevin2.getStepSize());
+    ASSERT_EQUAL(langevin1.getTemperature(), langevin2.getTemperature());
+    ASSERT_EQUAL(langevin1.getFriction(), langevin2.getFriction());
+    ASSERT_EQUAL(langevin1.getRandomNumberSeed(), langevin2.getRandomNumberSeed());
+    VerletIntegrator& verlet1 = dynamic_cast<VerletIntegrator&>(integ.getIntegrator(1));
+    VerletIntegrator& verlet2 = dynamic_cast<VerletIntegrator&>(integ2->getIntegrator(1));
+    ASSERT_EQUAL(verlet1.getConstraintTolerance(), verlet2.getConstraintTolerance());
+    ASSERT_EQUAL(verlet1.getStepSize(), verlet2.getStepSize());
+    delete integ2;
+}
+
 int main() {
     try {
         testSerializeBrownianIntegrator();
@@ -193,6 +217,7 @@ int main() {
         testSerializeVariableLangevinIntegrator();
         testSerializeVariableVerletIntegrator();
         testSerializeLangevinIntegrator(); 
+        testSerializeCompoundIntegrator();
     }
     catch(const exception& e) {
 		return 1;

tests/TestCompoundIntegrator.h

@@ -55,8 +55,8 @@ void testChangingIntegrator() {
     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));
+    integrator.addIntegrator(new LangevinIntegrator(300.0, 10.0, 0.011));
+    integrator.addIntegrator(new BrownianIntegrator(300.0, 10.0, 0.012));
     Context context(system, integrator, platform);
     ASSERT_EQUAL(0, integrator.getCurrentIntegrator());
     vector<Vec3> positions(2);
@@ -83,7 +83,7 @@ void testChangingIntegrator() {
             ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
             integrator.step(1);
         }
-        ASSERT_EQUAL_TOL(1.0, context.getState(0).getTime(), 1e-5);
+        ASSERT_EQUAL_TOL(100*0.01, context.getState(0).getTime(), 1e-5);
 
         // Switch to the Langevin integrator and make sure that it heats up.
 
@@ -97,6 +97,7 @@ void testChangingIntegrator() {
         }
         double expectedKE = 0.5*2*3*BOLTZ*300.0;
         ASSERT_USUALLY_EQUAL_TOL(expectedKE, ke/1000, 0.1);
+        ASSERT_EQUAL_TOL(100*0.01+10100*0.011, context.getState(0).getTime(), 1e-5);
         
         // Now reinitialize the context and repeat all of these tests to make sure that works correctly.
         
@@ -146,6 +147,67 @@ void testChangingParameters() {
     }
 }
 
+void testDifferentStepSizes() {
+    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.005));
+    integrator.addIntegrator(new VerletIntegrator(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);
+    context.setPositions(positions);
+
+    // Integrate with the first Verlet integrator and compare it to the analytical solution.
+
+    const double freq = 1.0;
+    double expectedTime = 0;
+    for (int i = 0; i < 100; ++i) {
+        State state = context.getState(State::Positions);
+        double time = state.getTime();
+        ASSERT_EQUAL_TOL(expectedTime, time, 1e-5);
+        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);
+        integrator.step(1);
+        expectedTime += 0.005;
+    }
+    
+    // Now switch to the second Verlet integrator which has a different step size.
+
+    integrator.setCurrentIntegrator(1);
+    for (int i = 0; i < 100; ++i) {
+        State state = context.getState(State::Positions);
+        double time = state.getTime();
+        ASSERT_EQUAL_TOL(expectedTime, time, 1e-5);
+        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);
+        integrator.step(1);
+        expectedTime += 0.01;
+    }
+    
+    // Finally, switch back to the first one again.
+
+    integrator.setCurrentIntegrator(0);
+    for (int i = 0; i < 100; ++i) {
+        State state = context.getState(State::Positions);
+        double time = state.getTime();
+        ASSERT_EQUAL_TOL(expectedTime, time, 1e-5);
+        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);
+        integrator.step(1);
+        expectedTime += 0.005;
+    }
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -153,6 +215,7 @@ int main(int argc, char* argv[]) {
         initializeTests(argc, argv);
         testChangingIntegrator();
         testChangingParameters();
+        testDifferentStepSizes();
         runPlatformTests();
     }
     catch(const exception& e) {