Merge pull request #929 from peastman/integrators

Moved the MTS and aMD integrators into the main repository

docs-source/usersguide/application.rst

@@ -1121,6 +1121,25 @@ integrator is:
 The parameter is the integration error tolerance (0.001), whose meaning is the
 same as for the Langevin integrator.
 
+Multiple Time Step Integrator
+-----------------------------
+
+The :class:`MTSIntegrator` class implements the rRESPA multiple time step
+algorithm\ :cite:`Tuckerman1992`.  This allows some forces in the system to be evaluated more
+frequently than others.  For details on how to use it, consult the API
+documentation.
+
+aMD Integrator
+--------------
+
+There are three different integrator types that implement variations of the
+aMD\ :cite:`Hamelberg2007` accelerated sampling algorithm: :class:`AMDIntegrator`,
+:class:`AMDForceGroupIntegrator`, and :class:`DualAMDIntegrator`.  They
+perform integration on a modified potential energy surface to allow much faster
+sampling of conformations.  For details on how to use them, consult the API
+documentation.
+
+
 Temperature Coupling
 ====================
 

docs-source/usersguide/references.bib

@@ -90,6 +90,17 @@
    type = {Journal Article}
 }
 
+@article{Hamelberg2007,
+   author={Hamelberg, Donald and de Oliveira, Cesar Augusto F. and McCammon, J. Andrew},
+   title={Sampling of slow diffusive conformational transitions with accelerated molecular dynamics},
+   journal={Journal of Chemical Physics},
+   volume={127},
+   number={15},
+   pages={155102},
+   year={2007},
+   type = {Journal Article}
+}
+
 @article{Hawkins1995
    author = {Hawkins, Gregory D. and Cramer, Christopher J. and Truhlar, Donald G.},
    title = {Pairwise solute descreening of solute charges from a dielectric medium},
@@ -440,6 +451,17 @@
    type = {Journal Article}
 }
 
+@article{Tuckerman1992,
+   author={Tuckerman, M. and Berne, Bruce J. and Martyna, Glenn J.},
+   title={Reversible multiple time scale molecular dynamics},
+   journal = {Journal of Chemical Physics},
+   volume={97},
+   number={3},
+   pages={1990-2001},
+   year={1992},
+   type = {Journal Article}
+}
+
 @article{Uberuaga2004,
   author =   {Blas P. Uberuaga and Marian Anghel and Arthur
                   F. Voter},

wrappers/python/simtk/openmm/__init__.py

@@ -20,6 +20,8 @@ if sys.platform == 'win32':
 
 from simtk.openmm.openmm import *
 from simtk.openmm.vec3 import Vec3
+from simtk.openmm.mtsintegrator import MTSIntegrator
+from simtk.openmm.amd import AMDIntegrator, AMDForceGroupIntegrator, DualAMDIntegrator
 
 if os.getenv('OPENMM_PLUGIN_DIR') is None and os.path.isdir(version.openmm_library_path):
     pluginLoadedLibNames = Platform.loadPluginsFromDirectory(os.path.join(version.openmm_library_path, 'plugins'))

wrappers/python/simtk/openmm/amd.py

+"""
+amd.py: Implements the aMD integration method.
+
+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) 2012 Stanford University and the Authors.
+Authors: Peter Eastman, Steffen Lindert
+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.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+from simtk.openmm import CustomIntegrator
+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
+        """
+        CustomIntegrator.__init__(self, dt)
+        self.addGlobalVariable("alpha", alpha)
+        self.addGlobalVariable("E", E)
+        self.addPerDofVariable("oldx", 0)
+        self.addUpdateContextState();
+        self.addComputePerDof("v", "v+dt*fprime/m; fprime=f*((1-modify) + modify*(alpha/(alpha+E-energy))^2); modify=step(E-energy)")
+        self.addComputePerDof("oldx", "x")
+        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()
+        E = self.getE()
+        if not is_quantity(energy):
+            energy = energy*kilojoules_per_mole # Assume kJ/mole
+        if (energy > E):
+            return energy
+        return energy+(E-energy)*(E-energy)/(alpha+E-energy)
+
+
+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
+        """
+        CustomIntegrator.__init__(self, dt)
+        self.addGlobalVariable("alphaGroup", alphaGroup)
+        self.addGlobalVariable("EGroup", EGroup)
+        self.addGlobalVariable("groupEnergy", 0)
+        self.addPerDofVariable("oldx", 0)
+        self.addPerDofVariable("fg", 0)
+        self.addUpdateContextState();
+        self.addComputeGlobal("groupEnergy", "energy"+str(group))
+        self.addComputePerDof("fg", "f"+str(group))
+        self.addComputePerDof("v", "v+dt*fprime/m; fprime=fother + fg*((1-modify) + modify*(alphaGroup/(alphaGroup+EGroup-groupEnergy))^2); fother=f-fg; modify=step(EGroup-groupEnergy)")
+        self.addComputePerDof("oldx", "x")
+        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
+        """
+        alphaGroup = self.getAlphaGroup()
+        EGroup = self.getEGroup()
+        if not is_quantity(groupEnergy):
+            groupEnergy = groupEnergy*kilojoules_per_mole # Assume kJ/mole
+        dE = 0.0*kilojoules_per_mole
+        if (groupEnergy < EGroup):
+            dE = dE + (EGroup-groupEnergy)*(EGroup-groupEnergy)/(alphaGroup+EGroup-groupEnergy)
+        return groupEnergy+dE
+
+
+
+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
+        """
+        CustomIntegrator.__init__(self, dt)
+        self.addGlobalVariable("alphaTotal", alphaTotal)
+        self.addGlobalVariable("ETotal", ETotal)
+        self.addGlobalVariable("alphaGroup", alphaGroup)
+        self.addGlobalVariable("EGroup", EGroup)
+        self.addGlobalVariable("groupEnergy", 0)
+        self.addPerDofVariable("oldx", 0)
+        self.addPerDofVariable("fg", 0)
+        self.addUpdateContextState();
+        self.addComputeGlobal("groupEnergy", "energy"+str(group))
+        self.addComputePerDof("fg", "f"+str(group))
+        self.addComputePerDof("v", """v+dt*fprime/m;
+                                      fprime=fprime1 + fprime2;
+                                      fprime2=fg*((1-modifyGroup) + modifyGroup*(alphaGroup/(alphaGroup+EGroup-groupEnergy))^2);
+                                      fprime1=fother*((1-modifyTotal) + modifyTotal*(alphaTotal/(alphaTotal+ETotal-energy))^2);
+                                      fother=f-fg;
+                                      modifyTotal=step(ETotal-energy); modifyGroup=step(EGroup-groupEnergy)""")
+        self.addComputePerDof("oldx", "x")
+        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
+        """
+        alphaTotal = self.getAlphaTotal()
+        ETotal = self.getETotal()
+        alphaGroup = self.getAlphaGroup()
+        EGroup = self.getEGroup()
+        if not is_quantity(totalEnergy):
+            totalEnergy = totalEnergy*kilojoules_per_mole # Assume kJ/mole
+        if not is_quantity(groupEnergy):
+            groupEnergy = groupEnergy*kilojoules_per_mole # Assume kJ/mole
+        dE = 0.0*kilojoules_per_mole
+        if (totalEnergy < ETotal):
+            dE = dE + (ETotal-totalEnergy)*(ETotal-totalEnergy)/(alphaTotal+ETotal-totalEnergy)
+        if (groupEnergy < EGroup):
+            dE = dE + (EGroup-groupEnergy)*(EGroup-groupEnergy)/(alphaGroup+EGroup-groupEnergy)
+        return totalEnergy+dE

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

@@ -4,7 +4,7 @@ OpenMM Application
 __docformat__ = "epytext en"
 
 __author__ = "Peter Eastman"
-__copyright__ = "Copyright 2011, Stanford University and Peter Eastman"
+__copyright__ = "Copyright 2015, Stanford University and Peter Eastman"
 __credits__ = []
 __license__ = "MIT"
 __maintainer__ = "Peter Eastman"

wrappers/python/simtk/openmm/mtsintegrator.py

+"""
+respa.py: Implements the rRESPA multiple time step integration method.
+
+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) 2013-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.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+from simtk.openmm import CustomIntegrator
+
+class MTSIntegrator(CustomIntegrator):
+    """MTSIntegrator implements the rRESPA multiple time step integration algorithm.
+    
+    This integrator allows different forces to be evaluated at different frequencies,
+    for example to evaluate the expensive, slowly changing forces less frequently than
+    the inexpensive, quickly changing forces.
+    
+    To use it, you must first divide your forces into two or more groups (by calling
+    setForceGroup() on them) that should be evaluated at different frequencies.  When
+    you create the integrator, you provide a tuple for each group specifying the index
+    of the force group and the frequency (as a fraction of the outermost time step) at
+    which to evaluate it.  For example:
+    
+    <pre>
+    integrator = MTSIntegrator(4*femtoseconds, [(0,1), (1,2), (2,8)])
+    </pre>
+    
+    This specifies that the outermost time step is 4 fs, so each step of the integrator
+    will advance time by that much.  It also says that force group 0 should be evaluated
+    once per time step, force group 1 should be evaluated twice per time step (every 2 fs),
+    and force group 2 should be evaluated eight times per time step (every 0.5 fs).
+    
+    A common use of this algorithm is to evaluate reciprocal space nonbonded interactions
+    less often than the bonded and direct space nonbonded interactions.  The following
+    example looks up the NonbondedForce, sets the reciprocal space interactions to their
+    own force group, and then creates an integrator that evaluates them once every 4 fs,
+    but all other interactions every 2 fs.
+    
+    <pre>
+    nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
+    nonbonded.setReciprocalSpaceForceGroup(1)
+    integrator = MTSIntegrator(4*femtoseconds, [(1,1), (0,2)])
+    </pre>
+    
+    For details, see Tuckerman et al., J. Chem. Phys. 97(3) pp. 1990-2001 (1992).
+    """
+    
+    def __init__(self, dt, groups):
+        """Create an MTSIntegrator.
+        
+        Parameters:
+         - dt (time) The largest (outermost) integration time step to use
+         - groups (list) A list of tuples defining the force groups.  The first element of each
+           tuple is the force group index, and the second element is the number of times that force
+           group should be evaluated in one time step.
+        """
+        if len(groups) == 0:
+            raise ValueError("No force groups specified")
+        groups = sorted(groups, key=lambda x: x[1])
+        CustomIntegrator.__init__(self, dt)
+        self.addPerDofVariable("x1", 0)
+        self.addUpdateContextState();
+        self._createSubsteps(1, groups)
+        self.addConstrainVelocities();
+            
+    def _createSubsteps(self, parentSubsteps, groups):
+        group, substeps = groups[0]
+        stepsPerParentStep = substeps/parentSubsteps
+        if stepsPerParentStep < 1 or stepsPerParentStep != int(stepsPerParentStep):
+            raise ValueError("The number for substeps for each group must be a multiple of the number for the previous group")
+        if group < 0 or group > 31:
+            raise ValueError("Force group must be between 0 and 31")
+        for i in range(stepsPerParentStep):
+            self.addComputePerDof("v", "v+0.5*(dt/"+str(substeps)+")*f"+str(group)+"/m")
+            if len(groups) == 1:
+                self.addComputePerDof("x1", "x")
+                self.addComputePerDof("x", "x+(dt/"+str(substeps)+")*v")
+                self.addConstrainPositions();
+                self.addComputePerDof("v", "(x-x1)/(dt/"+str(substeps)+")");
+            else:
+                self._createSubsteps(substeps, groups[1:])
+            self.addComputePerDof("v", "v+0.5*(dt/"+str(substeps)+")*f"+str(group)+"/m")