Merge pull request #2701 from peastman/integrationgroups

Added Integrator.setIntegrationForceGroups()

.travis.yml

@@ -216,7 +216,7 @@ script:
   - python devtools/run-ctest.py --start-time $START_TIME
 
   - if [[ ! -z "${DOCS_DEPLOY}" && "${DOCS_DEPLOY}" = "true" ]]; then
-      pip install sphinx sphinxcontrib-bibtex sphinxcontrib-lunrsearch sphinxcontrib-autodoc_doxygen;
+      pip install sphinx==2.3.1 sphinxcontrib-bibtex sphinxcontrib-lunrsearch sphinxcontrib-autodoc_doxygen;
       make sphinxhtml;
       make sphinxpdf;
       make C++ApiDocs PythonApiDocs;

docs-source/usersguide/theory.rst

@@ -1623,7 +1623,8 @@ Force Groups
 ************
 
 It is possible to split the Force objects in a System into groups.  Those groups
-can then be evaluated independently of each other.  Some Force classes also
+can then be evaluated independently of each other.  This is done by calling
+:code:`setForceGroup()` on the Force.  Some Force classes also
 provide finer grained control over grouping.  For example, NonbondedForce allows
 direct space computations to be in one group and reciprocal space computations
 in a different group.
@@ -1631,8 +1632,15 @@ in a different group.
 The most important use of force groups is for implementing multiple time step
 algorithms with CustomIntegrator.  For example, you might evaluate the slowly
 changing nonbonded interactions less frequently than the quickly changing bonded
-ones.  It also is useful if you want the ability to query a subset of the forces
-acting on the system.
+ones.  This can be done by putting the slow and fast forces into separate
+groups, then using a :class:`MTSIntegrator` or :class:`MTSLangevinIntegrator`
+that evaluates the groups at different frequencies.
+
+Another important use is to define forces that are not used when integrating
+the equations of motion, but can still be queried efficiently.  To do this,
+call :code:`setIntegrationForceGroups()` on the :class:`Integrator`.  Any groups
+omitted will be ignored during simulation, but can be queried at any time by
+calling :code:`getState()`.
 
 Virtual Sites
 *************

openmmapi/include/openmm/Integrator.h

@@ -85,6 +85,18 @@ public:
      * @param steps   the number of time steps to take
      */
     virtual void step(int steps) = 0;
+    /**
+     * Get which force groups to use for integration.  By default, all force groups
+     * are included.  This is interpreted as a set of bit flags: the forces from group i
+     * will be included if (groups&(1<<i)) != 0.
+     */
+    virtual int getIntegrationForceGroups() const;
+    /**
+     * Set which force groups to use for integration.  By default, all force groups
+     * are included.  This is interpreted as a set of bit flags: the forces from group i
+     * will be included if (groups&(1<<i)) != 0.
+     */
+    virtual void setIntegrationForceGroups(int groups);
 protected:
     friend class Context;
     friend class ContextImpl;
@@ -166,6 +178,7 @@ protected:
     }
 private:
     double stepSize, constraintTol;
+    int forceGroups;
 };
 
 } // namespace OpenMM

openmmapi/include/openmm/LocalEnergyMinimizer.h

@@ -43,6 +43,10 @@ namespace OpenMM {
  * force to the potential function.  The strength of the restraining force is steadily increased
  * until the minimum energy configuration satisfies all constraints to within the tolerance
  * specified by the Context's Integrator.
+ * 
+ * Energy minimization is done using the force groups defined by the Integrator.
+ * If you have called setIntegrationForceGroups() on it to restrict the set of forces
+ * used for integration, only the energy of the included forces will be minimized.
  */
 
 class OPENMM_EXPORT LocalEnergyMinimizer {

openmmapi/src/BrownianIntegrator.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) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -80,7 +80,7 @@ void BrownianIntegrator::step(int steps) {
         throw OpenMMException("This Integrator is not bound to a context!");  
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         kernel.getAs<IntegrateBrownianStepKernel>().execute(*context, *this);
     }
 }

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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -39,7 +39,7 @@
 
 using namespace OpenMM;
 
-Integrator::Integrator() : owner(NULL), context(NULL) {
+Integrator::Integrator() : owner(NULL), context(NULL), forceGroups(0xFFFFFFFF) {
 }
 
 Integrator::~Integrator() {
@@ -69,6 +69,14 @@ void Integrator::setConstraintTolerance(double tol) {
     constraintTol = tol;
 }
 
+int Integrator::getIntegrationForceGroups() const {
+    return forceGroups;
+}
+
+void Integrator::setIntegrationForceGroups(int groups) {
+    forceGroups = groups;
+}
+
 std::vector<Vec3> Integrator::getVelocitiesForTemperature(const System &system, double temperature, int randomSeed) const {
     // Generate the list of Gaussian random numbers.
     OpenMM_SFMT::SFMT sfmt;

openmmapi/src/LangevinIntegrator.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) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -76,7 +76,7 @@ void LangevinIntegrator::step(int steps) {
         throw OpenMMException("This Integrator is not bound to a context!");  
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         kernel.getAs<IntegrateLangevinStepKernel>().execute(*context, *this);
     }
 }

openmmapi/src/LangevinMiddleIntegrator.cpp

@@ -80,7 +80,7 @@ void LangevinMiddleIntegrator::step(int steps) {
         throw OpenMMException("This Integrator is not bound to a context!");  
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         kernel.getAs<IntegrateLangevinMiddleStepKernel>().execute(*context, *this);
     }
 }

openmmapi/src/LocalEnergyMinimizer.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-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -79,7 +79,7 @@ struct MinimizerData {
 static double computeForcesAndEnergy(Context& context, const vector<Vec3>& positions, lbfgsfloatval_t *g) {
     context.setPositions(positions);
     context.computeVirtualSites();
-    State state = context.getState(State::Forces | State::Energy);
+    State state = context.getState(State::Forces | State::Energy, false, context.getIntegrator().getIntegrationForceGroups());
     const vector<Vec3>& forces = state.getForces();
     const System& system = context.getSystem();
     for (int i = 0; i < forces.size(); i++) {

openmmapi/src/MonteCarloAnisotropicBarostatImpl.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -72,7 +72,8 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     
     // Compute the current potential energy.
     
-    double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+    int groups = context.getIntegrator().getIntegrationForceGroups();
+    double initialEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
     double pressure;
     
     // Choose which axis to modify at random.
@@ -114,7 +115,7 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     
     // Compute the energy of the modified system.
     
-    double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+    double finalEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
     double kT = BOLTZ*context.getParameter(MonteCarloAnisotropicBarostat::Temperature());
     double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);
     if (w > 0 && SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() > std::exp(-w/kT)) {

openmmapi/src/MonteCarloBarostatImpl.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -68,7 +68,8 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context, bool& forc
 
     // Compute the current potential energy.
 
-    double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+    int groups = context.getIntegrator().getIntegrationForceGroups();
+    double initialEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
 
     // Modify the periodic box size.
 
@@ -83,7 +84,7 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context, bool& forc
 
     // Compute the energy of the modified system.
     
-    double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+    double finalEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
     double pressure = context.getParameter(MonteCarloBarostat::Pressure())*(AVOGADRO*1e-25);
     double kT = BOLTZ*context.getParameter(MonteCarloBarostat::Temperature());
     double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);

openmmapi/src/MonteCarloMembraneBarostatImpl.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-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -70,7 +70,8 @@ void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {
     
     // Compute the current potential energy.
     
-    double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+    int groups = context.getIntegrator().getIntegrationForceGroups();
+    double initialEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
     double pressure = context.getParameter(MonteCarloMembraneBarostat::Pressure())*(AVOGADRO*1e-25);
     double tension = context.getParameter(MonteCarloMembraneBarostat::SurfaceTension())*(AVOGADRO*1e-25);
     
@@ -115,7 +116,7 @@ void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {
     
     // Compute the energy of the modified system.
     
-    double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+    double finalEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
     double kT = BOLTZ*context.getParameter(MonteCarloMembraneBarostat::Temperature());
     double w = finalEnergy-initialEnergy + pressure*deltaVolume - tension*deltaArea - context.getMolecules().size()*kT*std::log(newVolume/volume);
     if (w > 0 && SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() > std::exp(-w/kT)) {

openmmapi/src/NoseHooverIntegrator.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) 2019 Stanford University and the Authors.           *
+ * Portions copyright (c) 2019-2020 Stanford University and the Authors.      *
  * Authors: Andreas Krämer and Andrew C. Simmonett                            *
  * Contributors:                                                              *
  *                                                                            *
@@ -337,7 +337,7 @@ void NoseHooverIntegrator::step(int steps) {
     for (int i = 0; i < steps; ++i) {
         if(context->updateContextState())
             forcesAreValid = false;
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         kernel.getAs<IntegrateNoseHooverStepKernel>().execute(*context, *this, forcesAreValid);
     }
 }

openmmapi/src/VariableLangevinIntegrator.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) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -80,7 +80,7 @@ void VariableLangevinIntegrator::step(int steps) {
         throw OpenMMException("This Integrator is not bound to a context!");  
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         setStepSize(kernel.getAs<IntegrateVariableLangevinStepKernel>().execute(*context, *this, std::numeric_limits<double>::infinity()));
     }
 }
@@ -90,7 +90,7 @@ void VariableLangevinIntegrator::stepTo(double time) {
         throw OpenMMException("This Integrator is not bound to a context!");
     while (time > context->getTime()) {
         context->updateContextState();
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         setStepSize(kernel.getAs<IntegrateVariableLangevinStepKernel>().execute(*context, *this, time));
     }
 }

openmmapi/src/VariableVerletIntegrator.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) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -74,7 +74,7 @@ void VariableVerletIntegrator::step(int steps) {
         throw OpenMMException("This Integrator is not bound to a context!");
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         setStepSize(kernel.getAs<IntegrateVariableVerletStepKernel>().execute(*context, *this, std::numeric_limits<double>::infinity()));
     }
 }
@@ -84,7 +84,7 @@ void VariableVerletIntegrator::stepTo(double time) {
         throw OpenMMException("This Integrator is not bound to a context!");  
     while (time > context->getTime()) {
         context->updateContextState();
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         setStepSize(kernel.getAs<IntegrateVariableVerletStepKernel>().execute(*context, *this, time));
     }
 }

openmmapi/src/VerletIntegrator.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) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -73,7 +73,7 @@ void VerletIntegrator::step(int steps) {
         throw OpenMMException("This Integrator is not bound to a context!");
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
-        context->calcForcesAndEnergy(true, false);
+        context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
         kernel.getAs<IntegrateVerletStepKernel>().execute(*context, *this);
     }
 }

platforms/common/src/CommonKernels.cpp

@@ -5709,7 +5709,7 @@ void CommonIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const No
 
     // If the atom reordering has occured, the forces from the previous step are permuted and thus invalid.
     // They need to be either sorted or recomputed; here we choose the latter.
-    if (!forcesAreValid || cc.getAtomsWereReordered()) context.calcForcesAndEnergy(true, false);
+    if (!forcesAreValid || cc.getAtomsWereReordered()) context.calcForcesAndEnergy(true, false, integrator.getIntegrationForceGroups());
 
     const auto& atomList = integrator.getAllThermostatedIndividualParticles();
     const auto& pairList = integrator.getAllThermostatedPairs();
@@ -6770,7 +6770,7 @@ void CommonIntegrateCustomStepKernel::prepareForComputation(ContextImpl& context
 
         // Record the variable names and flags for the force and energy in each step.
 
-        forceGroupFlags.resize(numSteps, -1);
+        forceGroupFlags.resize(numSteps, integrator.getIntegrationForceGroups());
         vector<string> forceGroupName;
         vector<string> energyGroupName;
         for (int i = 0; i < 32; i++) {

platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp

 
-/* Portions copyright (c) 2011-2017 Stanford University and Simbios.
+/* Portions copyright (c) 2011-2020 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -171,7 +171,7 @@ void ReferenceCustomDynamics::initialize(ContextImpl& context, vector<double>& m
 
     // Record the force group flags for each step.
 
-    forceGroupFlags.resize(numSteps, -1);
+    forceGroupFlags.resize(numSteps, integrator.getIntegrationForceGroups());
     for (int i = 0; i < numSteps; i++)
         if (forceGroup[i] > -1)
             forceGroupFlags[i] = 1<<forceGroup[i];

platforms/reference/src/SimTKReference/ReferenceNoseHooverDynamics.cpp

@@ -26,6 +26,7 @@
 #include <iostream>
 #include <sstream>
 
+#include "openmm/Integrator.h"
 #include "openmm/OpenMMException.h"
 #include "SimTKOpenMMUtilities.h"
 #include "openmm/internal/ContextImpl.h"
@@ -55,7 +56,7 @@ void ReferenceNoseHooverDynamics::step1(OpenMM::ContextImpl &context, const Open
                                           double maxPairDistance) {
 
     // first-time-through initialization
-    if (!forcesAreValid) context.calcForcesAndEnergy(true, false);
+    if (!forcesAreValid) context.calcForcesAndEnergy(true, false, context.getIntegrator().getIntegrationForceGroups());
 
     if (getTimeStep() == 0) {
        // invert masses

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.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) 2008-2018 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Mark Friedrichs                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -1782,7 +1782,7 @@ void CudaCalcAmoebaMultipoleForceKernel::ensureMultipolesValid(ContextImpl& cont
         }
     }
     if (!multipolesAreValid)
-        context.calcForcesAndEnergy(false, false, -1);
+        context.calcForcesAndEnergy(false, false, context.getIntegrator().getIntegrationForceGroups());
 }
 
 
@@ -3650,7 +3650,7 @@ void CudaCalcHippoNonbondedForceKernel::ensureMultipolesValid(ContextImpl& conte
         }
     }
     if (!multipolesAreValid)
-        context.calcForcesAndEnergy(false, false, -1);
+        context.calcForcesAndEnergy(false, false, context.getIntegrator().getIntegrationForceGroups());
 }
 
 void CudaCalcHippoNonbondedForceKernel::getLabFramePermanentDipoles(ContextImpl& context, vector<Vec3>& dipoles) {