Merge master

plugins/amoeba/openmmapi/include/openmm/AmoebaStretchBendForce.h

@@ -81,15 +81,15 @@ public:
     /**
      * Get the force field parameters for a stretch-bend term.
      * 
-     * @param index         the index of the stretch-bend for which to get parameters
-     * @param particle1     the index of the first particle connected by the stretch-bend
-     * @param particle2     the index of the second particle connected by the stretch-bend
-     * @param particle3     the index of the third particle connected by the stretch-bend
-     * @param lengthAB      the equilibrium length of the stretch-bend in bond ab [particle1, particle2], measured in nm
-     * @param lengthCB      the equilibrium length of the stretch-bend in bond cb [particle3, particle2], measured in nm
-     * @param angle         the equilibrium angle in radians
-     * @param k1            the force constant of the product of bond ab and angle a-b-c
-     * @param k2            the force constant of the product of bond bc and angle a-b-c
+     * @param index              the index of the stretch-bend for which to get parameters
+     * @param[out] particle1     the index of the first particle connected by the stretch-bend
+     * @param[out] particle2     the index of the second particle connected by the stretch-bend
+     * @param[out] particle3     the index of the third particle connected by the stretch-bend
+     * @param[out] lengthAB      the equilibrium length of the stretch-bend in bond ab [particle1, particle2], measured in nm
+     * @param[out] lengthCB      the equilibrium length of the stretch-bend in bond cb [particle3, particle2], measured in nm
+     * @param[out] angle         the equilibrium angle in radians
+     * @param[out] k1            the force constant of the product of bond ab and angle a-b-c
+     * @param[out] k2            the force constant of the product of bond bc and angle a-b-c
      */
     void getStretchBendParameters(int index, int& particle1, int& particle2, int& particle3, double& lengthAB,
                                   double& lengthCB, double& angle, double& k1, double& k2) const;

plugins/amoeba/openmmapi/include/openmm/AmoebaTorsionTorsionForce.h

@@ -89,14 +89,14 @@ public:
     /**
      * Get the force field parameters for a torsion-torsion term.
      *
-     * @param index                     the index of the torsion-torsion for which to get parameters
-     * @param particle1                 the index of the first particle connected by the torsion-torsion
-     * @param particle2                 the index of the second particle connected by the torsion-torsion
-     * @param particle3                 the index of the third particle connected by the torsion-torsion
-     * @param particle4                 the index of the fourth particle connected by the torsion-torsion
-     * @param particle5                 the index of the fifth particle connected by the torsion-torsion
-     * @param chiralCheckAtomIndex      the index of the particle connected to particle3, but not particle2 or particle4 to be used in chirality check
-     * @param gridIndex                 the grid index
+     * @param index                          the index of the torsion-torsion for which to get parameters
+     * @param[out] particle1                 the index of the first particle connected by the torsion-torsion
+     * @param[out] particle2                 the index of the second particle connected by the torsion-torsion
+     * @param[out] particle3                 the index of the third particle connected by the torsion-torsion
+     * @param[out] particle4                 the index of the fourth particle connected by the torsion-torsion
+     * @param[out] particle5                 the index of the fifth particle connected by the torsion-torsion
+     * @param[out] chiralCheckAtomIndex      the index of the particle connected to particle3, but not particle2 or particle4 to be used in chirality check
+     * @param[out] gridIndex                 the grid index
      */
     void getTorsionTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, int& particle5, int& chiralCheckAtomIndex, int& gridIndex) const;
 
@@ -117,7 +117,7 @@ public:
     /**
      * Get the torsion-torsion grid at the specified index
      *
-     * @param gridIndex     the grid index
+     * @param  index        the grid index
      * @return grid         return grid reference
      */
     const std::vector<std::vector<std::vector<double> > >& getTorsionTorsionGrid(int index) const;
@@ -172,7 +172,7 @@ public:
     }
     TorsionTorsionInfo(int particle1, int particle2, int particle3, int particle4, int particle5, int chiralCheckAtomIndex, int gridIndex) :
                        particle1(particle1), particle2(particle2), particle3(particle3),
-                       particle4(particle4), particle5(particle5), gridIndex(gridIndex), chiralCheckAtomIndex(chiralCheckAtomIndex) {
+                       particle4(particle4), particle5(particle5), chiralCheckAtomIndex(chiralCheckAtomIndex), gridIndex(gridIndex) {
 
     }
 };

plugins/amoeba/openmmapi/include/openmm/AmoebaWcaDispersionForce.h

@@ -75,9 +75,9 @@ public:
     /**
      * Get the force field parameters for a WCA dispersion particle.
      * 
-     * @param particleIndex   the particle index
-     * @param radius          radius
-     * @param epsilon         epsilon 
+     * @param particleIndex        the particle index
+     * @param[out] radius          radius
+     * @param[out] epsilon         epsilon 
      */
     void getParticleParameters(int particleIndex, double& radius, double& epsilon) const;
 

plugins/drude/openmmapi/include/openmm/DrudeForce.h

@@ -44,7 +44,7 @@ namespace OpenMM {
  * it applies: an anisotropic harmonic force connecting each Drude particle to its parent particle; and
  * a screened Coulomb interaction between specific pairs of dipoles.  The latter is typically used between
  * closely bonded particles whose Coulomb interaction would otherwise be fully excluded.
- * 
+ *
  * To use this class, create a DrudeForce object, then call addParticle() once for each Drude particle in the
  * System to define its parameters.  After a particle has been added, you can modify its force field parameters
  * by calling setParticleParameters().  This will have no effect on Contexts that already exist unless you
@@ -91,19 +91,19 @@ public:
     /**
      * Get the parameters for a Drude particle.
      *
-     * @param index           the index of the Drude particle for which to get parameters
-     * @param particle        the index within the System of the Drude particle
-     * @param particle1       the index within the System of the particle to which the Drude particle is attached
-     * @param particle2       the index within the System of the second particle used for defining anisotropic polarizability.
-     *                        This may be set to -1, in which case aniso12 will be ignored.
-     * @param particle3       the index within the System of the third particle used for defining anisotropic polarizability.
-     *                        This may be set to -1, in which case aniso34 will be ignored.
-     * @param particle4       the index within the System of the fourth particle used for defining anisotropic polarizability.
-     *                        This may be set to -1, in which case aniso34 will be ignored.
-     * @param charge          The charge on the Drude particle
-     * @param polarizability  The isotropic polarizability
-     * @param aniso12         The scale factor for the polarizability along the direction defined by particle1 and particle2
-     * @param aniso34         The scale factor for the polarizability along the direction defined by particle3 and particle4
+     * @param index                the index of the Drude particle for which to get parameters
+     * @param[out] particle        the index within the System of the Drude particle
+     * @param[out] particle1       the index within the System of the particle to which the Drude particle is attached
+     * @param[out] particle2       the index within the System of the second particle used for defining anisotropic polarizability.
+     *                             This may be set to -1, in which case aniso12 will be ignored.
+     * @param[out] particle3       the index within the System of the third particle used for defining anisotropic polarizability.
+     *                             This may be set to -1, in which case aniso34 will be ignored.
+     * @param[out] particle4       the index within the System of the fourth particle used for defining anisotropic polarizability.
+     *                             This may be set to -1, in which case aniso34 will be ignored.
+     * @param[out] charge          The charge on the Drude particle
+     * @param[out] polarizability  The isotropic polarizability
+     * @param[out] aniso12         The scale factor for the polarizability along the direction defined by particle1 and particle2
+     * @param[out] aniso34         The scale factor for the polarizability along the direction defined by particle3 and particle4
      */
     void getParticleParameters(int index, int& particle, int& particle1, int& particle2, int& particle3, int& particle4, double& charge, double& polarizability, double& aniso12, double& aniso34) const;
     /**
@@ -135,16 +135,16 @@ public:
     int addScreenedPair(int particle1, int particle2, double thole);
     /**
      * Get the force field parameters for screened pair.
-     * 
-     * @param index      the index of the pair for which to get parameters
-     * @param particle1  the index within this Force of the first particle involved in the interaction
-     * @param particle2  the index within this Force of the second particle involved in the interaction
-     * @param thole      the Thole screening factor
+     *
+     * @param index           the index of the pair for which to get parameters
+     * @param[out] particle1  the index within this Force of the first particle involved in the interaction
+     * @param[out] particle2  the index within this Force of the second particle involved in the interaction
+     * @param[out] thole      the Thole screening factor
      */
     void getScreenedPairParameters(int index, int& particle1, int& particle2, double& thole) const;
     /**
      * Set the force field parameters for screened pair.
-     * 
+     *
      * @param index      the index of the pair for which to get parameters
      * @param particle1  the index within this Force of the first particle involved in the interaction
      * @param particle2  the index within this Force of the second particle involved in the interaction
@@ -156,7 +156,7 @@ public:
      * provides an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
      * Simply call setParticleParameters() and setScreenedPairParameters() to modify this object's parameters, then call
      * updateParametersInContext() to copy them over to the Context.
-     * 
+     *
      * This method has several limitations.  It can be used to modify the numeric parameters associated with a particle or
      * screened pair (polarizability, thole, etc.), but not the identities of the particles they involve.  It also cannot
      * be used to add new particles or screenedPairs, only to change the parameters of existing ones.

plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp

@@ -128,21 +128,15 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
             center *= 1.0/molecules[i].size();
 
             // Find the displacement to move it into the first periodic box.
-
-            int xcell = (int) floor(center[0]/periodicBoxSize[0][0]);
-            int ycell = (int) floor(center[1]/periodicBoxSize[1][1]);
-            int zcell = (int) floor(center[2]/periodicBoxSize[2][2]);
-            double dx = xcell*periodicBoxSize[0][0];
-            double dy = ycell*periodicBoxSize[1][1];
-            double dz = zcell*periodicBoxSize[2][2];
+            Vec3 diff;
+            diff -= periodicBoxSize[0]*static_cast<int>(center[0]/periodicBoxSize[0][0]);
+            diff -= periodicBoxSize[1]*static_cast<int>(center[1]/periodicBoxSize[1][1]);
+            diff -= periodicBoxSize[2]*static_cast<int>(center[2]/periodicBoxSize[2][2]);
 
             // Translate all the particles in the molecule.
-
             for (int j = 0; j < (int) molecules[i].size(); j++) {
                 Vec3& pos = positions[molecules[i][j]];
-                pos[0] -= dx;
-                pos[1] -= dy;
-                pos[2] -= dz;
+                pos -= diff;
             }
         }
 
@@ -170,7 +164,7 @@ double RPMDIntegrator::computeKineticEnergy() {
 
 void RPMDIntegrator::step(int steps) {
     if (context == NULL)
-        throw OpenMMException("This Integrator is not bound to a context!");    
+        throw OpenMMException("This Integrator is not bound to a context!");
     if (!hasSetPosition) {
         // Initialize the positions from the context.
 

serialization/include/openmm/serialization/CompoundIntegratorProxy.h

+#ifndef OPENMM_COMPOUND_INTEGRATOR_PROXY_H_
+#define OPENMM_COMPOUND_INTEGRATOR_PROXY_H_
+
+#include "openmm/serialization/XmlSerializer.h"
+
+namespace OpenMM {
+
+class CompoundIntegratorProxy : public SerializationProxy {
+public:
+    CompoundIntegratorProxy();
+    void serialize(const void* object, SerializationNode& node) const;
+    void* deserialize(const SerializationNode& node) const;
+};
+
+}
+
+#endif /*OPENMM_COMPOUND_INTEGRATOR_PROXY_H_*/
\ No newline at end of file

serialization/include/openmm/serialization/SerializationNode.h

@@ -105,7 +105,7 @@ public:
     /**
      * Determine whether this node has a property with a particular node.
      *
-     * @param the name of the property to check for
+     * @param name  the name of the property to check for
      */
     bool hasProperty(const std::string& name) const;
     /**

serialization/src/CompoundIntegratorProxy.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/serialization/CompoundIntegratorProxy.h"
+#include <OpenMM.h>
+
+using namespace std;
+using namespace OpenMM;
+
+CompoundIntegratorProxy::CompoundIntegratorProxy() : SerializationProxy("CompoundIntegrator") {
+}
+
+void CompoundIntegratorProxy::serialize(const void* object, SerializationNode& node) const {
+    node.setIntProperty("version", 1);
+    const CompoundIntegrator& integrator = *reinterpret_cast<const CompoundIntegrator*>(object);
+    node.setIntProperty("currentIntegrator", integrator.getCurrentIntegrator());
+    for (int i = 0; i < integrator.getNumIntegrators(); i++)
+        node.createChildNode("Integrator", &integrator.getIntegrator(i));
+}
+
+void* CompoundIntegratorProxy::deserialize(const SerializationNode& node) const {
+    if (node.getIntProperty("version") != 1)
+        throw OpenMMException("Unsupported version number");
+    CompoundIntegrator *integrator = new CompoundIntegrator();
+    for (int i = 0; i < node.getChildren().size(); i++)
+        integrator->addIntegrator(node.getChildren()[i].decodeObject<Integrator>());
+    integrator->setCurrentIntegrator(node.getIntProperty("currentIntegrator"));
+    return integrator;
+}
\ No newline at end of file

serialization/src/SerializationProxyRegistration.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-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -33,6 +33,7 @@
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMAPTorsionForce.h"
 #include "openmm/CMMotionRemover.h"
+#include "openmm/CompoundIntegrator.h"
 #include "openmm/CustomAngleForce.h"
 #include "openmm/CustomBondForce.h"
 #include "openmm/CustomCompoundBondForce.h"
@@ -65,6 +66,7 @@
 #include "openmm/serialization/AndersenThermostatProxy.h"
 #include "openmm/serialization/CMAPTorsionForceProxy.h"
 #include "openmm/serialization/CMMotionRemoverProxy.h"
+#include "openmm/serialization/CompoundIntegratorProxy.h"
 #include "openmm/serialization/CustomAngleForceProxy.h"
 #include "openmm/serialization/CustomBondForceProxy.h"
 #include "openmm/serialization/CustomCompoundBondForceProxy.h"
@@ -112,6 +114,7 @@ extern "C" void registerSerializationProxies() {
     SerializationProxy::registerProxy(typeid(BrownianIntegrator), new BrownianIntegratorProxy());
     SerializationProxy::registerProxy(typeid(CMAPTorsionForce), new CMAPTorsionForceProxy());
     SerializationProxy::registerProxy(typeid(CMMotionRemover), new CMMotionRemoverProxy());
+    SerializationProxy::registerProxy(typeid(CompoundIntegrator), new CompoundIntegratorProxy());
     SerializationProxy::registerProxy(typeid(Continuous1DFunction), new Continuous1DFunctionProxy());
     SerializationProxy::registerProxy(typeid(Continuous2DFunction), new Continuous2DFunctionProxy());
     SerializationProxy::registerProxy(typeid(Continuous3DFunction), new Continuous3DFunctionProxy());

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

+/* -------------------------------------------------------------------------- *
+ *                                   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.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);
+    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(100*0.01, 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);
+        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.
+        
+        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 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[]) {
+    try {
+        initializeTests(argc, argv);
+        testChangingIntegrator();
+        testChangingParameters();
+        testDifferentStepSizes();
+        runPlatformTests();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

tests/TestEnforcePeriodicBox.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) 2010-2015 Stanford University and the Authors.      *
+ * Authors: Robert McGibbon                                                   *
+ * 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/Context.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/Platform.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+
+using namespace OpenMM;
+using namespace std;
+
+void testTruncatedOctahedron() {
+    const int numMolecules = 5;
+    const int numParticles = numMolecules*2;
+    const float cutoff = 2.0;
+    Vec3 a(6.7929, 0, 0);
+    Vec3 b(-2.264163559406279, 6.404455775962287, 0);
+    Vec3 c(-2.264163559406279, -3.2019384603140684, 5.54658849047036);
+        
+    System system;
+    system.setDefaultPeriodicBoxVectors(a, b, c);
+    NonbondedForce* force = new NonbondedForce();
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    
+    force->setCutoffDistance(cutoff);
+    force->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        force->addParticle(-1, 0.2, 0.2);
+        force->addParticle(1, 0.2, 0.2);
+        positions[2*i] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        positions[2*i+1] = positions[2*i] + Vec3(1.0, 0.0, 0.0);
+        system.addConstraint(2*i, 2*i+1, 1.0);
+    }
+    system.addForce(force);
+    
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    context.setPositions(positions);
+    State initialState = context.getState(State::Positions | State::Energy, true);
+    double initialEnergy = initialState.getPotentialEnergy();
+
+    context.setState(initialState);
+    State finalState = context.getState(State::Positions | State::Energy, true);
+    double finalEnergy = finalState.getPotentialEnergy();
+
+    ASSERT_EQUAL_TOL(initialEnergy, finalEnergy, 1e-4);
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        testTruncatedOctahedron();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

wrappers/python/simtk/openmm/amd.py

@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
 Authors: Peter Eastman, Steffen Lindert
 Contributors:
 
-Permission is hereby granted, free of charge, to any person obtaining a 
+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,
@@ -37,22 +37,26 @@ from simtk.unit import kilojoules_per_mole, is_quantity
 
 class AMDIntegrator(CustomIntegrator):
     """AMDIntegrator implements the aMD integration algorithm.
-    
+
     The system is integrated based on a modified potential.  Whenever the energy V(r) is less than a
     cutoff value E, the following effective potential is used:
-    
+
     V*(r) = V(r) + (E-V(r))^2 / (alpha+E-V(r))
-    
+
     For details, see Hamelberg et al., J. Chem. Phys. 127, 155102 (2007).
     """
-    
+
     def __init__(self, dt, alpha, E):
         """Create an AMDIntegrator.
-        
-        Parameters:
-         - dt (time) The integration time step to use
-         - alpha (energy) The alpha parameter to use
-         - E (energy) The energy cutoff to use
+
+        Parameters
+        ----------
+        dt : time
+            The integration time step to use
+        alpha : energy
+            The alpha parameter to use
+        E : energy
+            The energy cutoff to use
         """
         CustomIntegrator.__init__(self, dt)
         self.addGlobalVariable("alpha", alpha)
@@ -64,23 +68,23 @@ class AMDIntegrator(CustomIntegrator):
         self.addComputePerDof("x", "x+dt*v")
         self.addConstrainPositions()
         self.addComputePerDof("v", "(x-oldx)/dt")
-    
+
     def getAlpha(self):
         """Get the value of alpha for the integrator."""
         return self.getGlobalVariable(0)*kilojoules_per_mole
-    
+
     def setAlpha(self, alpha):
         """Set the value of alpha for the integrator."""
         self.setGlobalVariable(0, alpha)
-    
+
     def getE(self):
         """Get the energy threshold E for the integrator."""
         return self.getGlobalVariable(1)*kilojoules_per_mole
-    
+
     def setE(self, E):
         """Set the energy threshold E for the integrator."""
         self.setGlobalVariable(1, E)
-    
+
     def getEffectiveEnergy(self, energy):
         """Given the actual potential energy of the system, return the value of the effective potential."""
         alpha = self.getAlpha()
@@ -94,21 +98,26 @@ class AMDIntegrator(CustomIntegrator):
 
 class AMDForceGroupIntegrator(CustomIntegrator):
     """AMDForceGroupIntegrator implements a single boost aMD integration algorithm.
-    
+
     This is similar to AMDIntegrator, but is applied based on the energy of a single force group
     (typically representing torsions).
-    
+
     For details, see Hamelberg et al., J. Chem. Phys. 127, 155102 (2007).
     """
-    
+
     def __init__(self, dt, group, alphaGroup, EGroup):
         """Create a AMDForceGroupIntegrator.
-        
-        Parameters:
-         - dt (time) The integration time step to use
-         - group (int) The force group to apply the boost to
-         - alphaGroup (energy) The alpha parameter to use for the boosted force group
-         - EGroup (energy) The energy cutoff to use for the boosted force group
+
+        Parameters
+        ----------
+        dt : time
+            The integration time step to use
+        group : int
+            The force group to apply the boost to
+        alphaGroup : energy
+            The alpha parameter to use for the boosted force group
+        EGroup : energy
+            The energy cutoff to use for the boosted force group
         """
         CustomIntegrator.__init__(self, dt)
         self.addGlobalVariable("alphaGroup", alphaGroup)
@@ -124,29 +133,35 @@ class AMDForceGroupIntegrator(CustomIntegrator):
         self.addComputePerDof("x", "x+dt*v")
         self.addConstrainPositions()
         self.addComputePerDof("v", "(x-oldx)/dt")
-    
+
     def getAlphaGroup(self):
         """Get the value of alpha for the boosted force group."""
         return self.getGlobalVariable(0)*kilojoules_per_mole
-    
+
     def setAlphaGroup(self, alpha):
         """Set the value of alpha for the boosted force group."""
         self.setGlobalVariable(0, alpha)
-    
+
     def getEGroup(self):
         """Get the energy threshold E for the boosted force group."""
         return self.getGlobalVariable(1)*kilojoules_per_mole
-    
+
     def setEGroup(self, E):
         """Set the energy threshold E for the boosted force group."""
         self.setGlobalVariable(1, E)
-    
+
     def getEffectiveEnergy(self, groupEnergy):
         """Given the actual group energy of the system, return the value of the effective potential.
-        
-        Parameters:
-          - groupEnergy (energy): the actual potential energy of the boosted force group
-        Returns: the value of the effective potential
+
+        Parameters
+        ----------
+        groupEnergy : energy
+            the actual potential energy of the boosted force group
+
+        Returns
+        -------
+        energy
+            the value of the effective potential
         """
         alphaGroup = self.getAlphaGroup()
         EGroup = self.getEGroup()
@@ -161,24 +176,31 @@ class AMDForceGroupIntegrator(CustomIntegrator):
 
 class DualAMDIntegrator(CustomIntegrator):
     """DualAMDIntegrator implements a dual boost aMD integration algorithm.
-    
+
     This is similar to AMDIntegrator, but two different boosts are applied to the potential:
     one based on the total energy, and one based on the energy of a single force group
     (typically representing torsions).
 
     For details, see Hamelberg et al., J. Chem. Phys. 127, 155102 (2007).
     """
-    
+
     def __init__(self, dt, group, alphaTotal, ETotal, alphaGroup, EGroup):
         """Create a DualAMDIntegrator.
 
-        Parameters:
-         - dt (time) The integration time step to use
-         - group (int) The force group to apply the second boost to
-         - alphaTotal (energy) The alpha parameter to use for the total energy
-         - ETotal (energy) The energy cutoff to use for the total energy
-         - alphaGroup (energy) The alpha parameter to use for the boosted force group
-         - EGroup (energy) The energy cutoff to use for the boosted force group
+        Parameters
+        ----------
+        dt : time
+            The integration time step to use
+        group : int
+            The force group to apply the second boost to
+        alphaTotal : energy
+            The alpha parameter to use for the total energy
+        ETotal : energy
+            The energy cutoff to use for the total energy
+        alphaGroup : energy
+            The alpha parameter to use for the boosted force group
+        EGroup : energy
+            The energy cutoff to use for the boosted force group
         """
         CustomIntegrator.__init__(self, dt)
         self.addGlobalVariable("alphaTotal", alphaTotal)
@@ -201,46 +223,53 @@ class DualAMDIntegrator(CustomIntegrator):
         self.addComputePerDof("x", "x+dt*v")
         self.addConstrainPositions()
         self.addComputePerDof("v", "(x-oldx)/dt")
-    
+
     def getAlphaTotal(self):
         """Get the value of alpha for the total energy."""
         return self.getGlobalVariable(0)*kilojoules_per_mole
-    
+
     def setAlphaTotal(self, alpha):
         """Set the value of alpha for the total energy."""
         self.setGlobalVariable(0, alpha)
-    
+
     def getETotal(self):
         """Get the energy threshold E for the total energy."""
         return self.getGlobalVariable(1)*kilojoules_per_mole
-    
+
     def setETotal(self, E):
         """Set the energy threshold E for the total energy."""
         self.setGlobalVariable(1, E)
-    
+
     def getAlphaGroup(self):
         """Get the value of alpha for the boosted force group."""
         return self.getGlobalVariable(2)*kilojoules_per_mole
-    
+
     def setAlphaGroup(self, alpha):
         """Set the value of alpha for the boosted force group."""
         self.setGlobalVariable(2, alpha)
-    
+
     def getEGroup(self):
         """Get the energy threshold E for the boosted force group."""
         return self.getGlobalVariable(3)*kilojoules_per_mole
-    
+
     def setEGroup(self, E):
         """Set the energy threshold E for the boosted force group."""
         self.setGlobalVariable(3, E)
-    
+
     def getEffectiveEnergy(self, totalEnergy, groupEnergy):
         """Given the actual potential energy of the system, return the value of the effective potential.
-        
-        Parameters:
-         - totalEnergy (energy): the actual potential energy of the whole system
-         - groupEnergy (energy): the actual potential energy of the boosted force group
-        Returns: the value of the effective potential
+
+        Parameters
+        ----------
+        totalEnergy : energy
+            the actual potential energy of the whole system
+        groupEnergy : energy
+            the actual potential energy of the boosted force group
+
+        Returns
+        -------
+        energy
+            the value of the effective potential
         """
         alphaTotal = self.getAlphaTotal()
         ETotal = self.getETotal()

wrappers/python/simtk/openmm/app/amberinpcrdfile.py

@@ -60,12 +60,17 @@ class AmberInpcrdFile(object):
     def __init__(self, file, loadVelocities=None, loadBoxVectors=None):
         """Load an inpcrd file.
 
-        An inpcrd file contains atom positions and, optionally, velocities and periodic box dimensions.
-
-        Parameters:
-         - file (string) the name of the file to load
-         - loadVelocities (boolean=None) deprecated. Velocities are loaded automatically if present
-         - loadBoxVectors (boolean=None) deprecated. Box vectors are loaded automatically if present
+        An inpcrd file contains atom positions and, optionally, velocities and
+        periodic box dimensions.
+
+        Parameters
+        ----------
+        file : str
+             The name of the file to load
+        loadVelocities : bool
+             Deprecated. Velocities are loaded automatically if present
+        loadBoxVectors : bool
+            Deprecated. Box vectors are loaded automatically if present
         """
         self.file = file
         if loadVelocities is not None or loadBoxVectors is not None:
@@ -84,8 +89,11 @@ class AmberInpcrdFile(object):
     def getPositions(self, asNumpy=False):
         """Get the atomic positions.
 
-        Parameters:
-         - asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s
+        Parameters
+        ----------
+        asNumpy : bool=False
+            if true, the values are returned as a numpy array instead of a list
+            of Vec3s
         """
         if asNumpy:
             if self._numpyPositions is None:
@@ -97,8 +105,10 @@ class AmberInpcrdFile(object):
     def getVelocities(self, asNumpy=False):
         """Get the atomic velocities.
 
-        Parameters:
-         - asNumpy (boolean=False) if true, the vectors are returned as numpy arrays instead of Vec3s
+        Parameters
+        ----------
+        asNumpy : bool=False
+            if true, the vectors are returned as numpy arrays instead of Vec3s
         """
         if self.velocities is None:
             raise AttributeError('velocities not found in %s' % self.file)
@@ -112,8 +122,11 @@ class AmberInpcrdFile(object):
     def getBoxVectors(self, asNumpy=False):
         """Get the periodic box vectors.
 
-        Parameters:
-         - asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s
+        Parameters
+        ----------
+        asNumpy : bool=False
+            if true, the values are returned as a numpy array instead of a list
+            of Vec3s
         """
         if self.boxVectors is None:
             raise AttributeError('Box information not found in %s' % self.file)