Skip to content

GitLab

Unverified Commit 2e86184b authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2352 from peastman/sampling 

Added simulated tempering and metadynamics scripts
parents ec5f7a18 78804789
Hide whitespace changes
Inline Side-by-side
Showing with 723 additions and 11 deletions
docs-source/usersguide/application.rst
View file @ 2e86184b
......@@ -1101,16 +1101,6 @@ algorithm\ :cite:`Tuckerman1992`. This allows some forces in the system to be e
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.
Compound Integrator
-------------------
......@@ -1371,6 +1361,55 @@ a checkpoint file every 5,000 steps, for example:
Note that the checkpoint reporter will overwrite the last checkpoint file.
Enhanced Sampling Methods
=========================
In many situations, the goal of a simulation is to sample the range of configurations
accessible to a system. It does not matter whether the simulation represents a
single, physically realistic trajectory, only whether it produces a correct distribution
of states. In this case, a variety of methods can be used to sample configuration
space much more quickly and efficiently than a single physical trajectory would.
These are known as enhanced sampling methods. OpenMM offers several that you
can choose from. They are briefly described here. Consult the API documentation
for more detailed descriptions and example code.
Simulated Tempering
-------------------
Simulated tempering\ :cite:`Marinari1992` involves making frequent changes to the
temperature of a simulation. At high temperatures, it can quickly cross energy barriers
and explore a wide range of configurations. At lower temperatures, it more thoroughly
explores each local region of configuration space. This is a powerful method to
speed up sampling when you do not know in advance what motions you want to sample.
Simply specify the range of temperatures to simulate and the algorithm handles
everything for you mostly automatically.
Metadynamics
------------
Metadynamics\ :cite:`Barducci2008` is used when you do know in advance what
motions you want to sample. You specify one or more collective variables, and the
algorithm adds a biasing potential to make the simulation explore a wide range of
values for those variables. It does this by periodically adding "bumps" to the biasing
potential at the current values of the collective variables. This encourages the simulation
to move away from regions it has already explored and sample a wide range of values.
At the end of the simulation, the biasing potential can be used to calculate the
free energy of the system as a function of the collective variables.
Accelerated Molecular Dynamics (aMD)
------------------------------------
aMD\ :cite:`Hamelberg2007` is another method that can be used when you do not know in
advance what motions you want to accelerate. It alters the potential energy surface
by adding a "boost" potential whenever the potential energy is below a threshold.
This makes local minima shallower and allows more frequent transitions between them.
The boost can be applied to the total potential energy, to just a subset of interactions
(typically the dihedral torsions), or both. There are separate integrator classes
for each of these options: :class:`AMDIntegrator`, :class:`AMDForceGroupIntegrator`,
and :class:`DualAMDIntegrator`.
.. _model-building-and-editing:
Model Building and Editing
......
docs-source/usersguide/references.bib
View file @ 2e86184b
......@@ -29,6 +29,18 @@
type = {Journal Article}
}
@article{Barducci2008,
title = {Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method},
author = {Barducci, Alessandro and Bussi, Giovanni and Parrinello, Michele},
journal = {Physical Review Letters},
volume = {100},
issue = {2},
pages = {020603},
year = {2008},
publisher = {American Physical Society},
doi = {10.1103/PhysRevLett.100.020603},
}
@article{Berendsen1987
author = {Berendsen, H. J. C. and Grigera, J. R. and Straatsma, T. P.},
title = {The missing term in effective pair potentials},
......@@ -290,6 +302,18 @@
type = {Journal Article}
}
@article{Marinari1992,
doi = {10.1209/0295-5075/19/6/002},
year = 1992,
publisher = {{IOP} Publishing},
volume = {19},
number = {6},
pages = {451--458},
author = {E Marinari and G Parisi},
title = {Simulated Tempering: A New Monte Carlo Scheme},
journal = {Europhysics Letters ({EPL})},
}
@article{Markland2008
author = {Markland, Thomas E. and Manolopoulos, David E.},
title = {An efficient ring polymer contraction scheme for imaginary time path integral simulations},
......
wrappers/python/simtk/openmm/app/__init__.py
View file @ 2e86184b
......@@ -6,7 +6,7 @@ from __future__ import absolute_import
__docformat__ = "epytext en"
__author__ = "Peter Eastman"
__copyright__ = "Copyright 2016, Stanford University and Peter Eastman"
__copyright__ = "Copyright 2016-2019, Stanford University and Peter Eastman"
__credits__ = []
__license__ = "MIT"
__maintainer__ = "Peter Eastman"
......@@ -32,6 +32,8 @@ from .checkpointreporter import CheckpointReporter
from .charmmcrdfiles import CharmmCrdFile, CharmmRstFile
from .charmmparameterset import CharmmParameterSet
from .charmmpsffile import CharmmPsfFile, CharmmPSFWarning
from .simulatedtempering import SimulatedTempering
from .metadynamics import Metadynamics, BiasVariable
# Enumerated values
......
wrappers/python/simtk/openmm/app/metadynamics.py 0 → 100644
View file @ 2e86184b
"""
metadynamics.py: Well-tempered metadynamics
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) 2018-2019 Stanford University and the Authors.
Authors: Peter Eastman
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.
"""
import simtk.openmm as mm
import simtk.unit as unit
from collections import namedtuple
from functools import reduce
import os
import re
try:
import numpy as np
except:
pass
class Metadynamics(object):
"""Performs metadynamics.
This class implements well-tempered metadynamics, as described in Barducci et al.,
"Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method"
(https://doi.org/10.1103/PhysRevLett.100.020603). You specify from one to three
collective variables whose sampling should be accelerated. A biasing force that
depends on the collective variables is added to the simulation. Initially the bias
is zero. As the simulation runs, Gaussian bumps are periodically added to the bias
at the current location of the simulation. This pushes the simulation away from areas
it has already explored, encouraging it to sample other regions. At the end of the
simulation, the bias function can be used to calculate the system's free energy as a
function of the collective variables.
To use the class you create a Metadynamics object, passing to it the System you want
to simulate and a list of BiasVariable objects defining the collective variables.
It creates a biasing force and adds it to the System. You then run the simulation
as usual, but call step() on the Metadynamics object instead of on the Simulation.
You can optionally specify a directory on disk where the current bias function should
periodically be written. In addition, it loads biases from any other files in the
same directory and includes them in the simulation. It loads files when the
Metqdynamics object is first created, and also checks for any new files every time it
updates its own bias on disk.
This serves two important functions. First, it lets you stop a metadynamics run and
resume it later. When you begin the new simulation, it will load the biases computed
in the earlier simulation and continue adding to them. Second, it provides an easy
way to parallelize metadynamics sampling across many computers. Just point all of
them to a shared directory on disk. Each process will save its biases to that
directory, and also load in and apply the biases added by other processes.
"""
def __init__(self, system, variables, temperature, biasFactor, height, frequency, saveFrequency=None, biasDir=None):
"""Create a Metadynamics object.
Parameters
----------
system: System
the System to simulate. A CustomCVForce implementing the bias is created and
added to the System.
variables: list of BiasVariables
the collective variables to sample
temperature: temperature
the temperature at which the simulation is being run. This is used in computing
the free energy.
biasFactor: float
used in scaling the height of the Gaussians added to the bias. The collective
variables are sampled as if the effective temperature of the simulation were
temperature*biasFactor.
height: energy
the initial height of the Gaussians to add
frequency: int
the interval in time steps at which Gaussians should be added to the bias potential
saveFrequency: int (optional)
the interval in time steps at which to write out the current biases to disk. At
the same time it writes biases, it also checks for updated biases written by other
processes and loads them in. This must be a multiple of frequency.
biasDir: str (optional)
the directory to which biases should be written, and from which biases written by
other processes should be loaded
"""
if not unit.is_quantity(temperature):
temperature = temperature*unit.kelvin
if not unit.is_quantity(height):
height = height*unit.kilojoules_per_mole
if biasFactor < 1.0:
raise ValueError('biasFactor must be >= 1')
if (saveFrequency is None and biasDir is not None) or (saveFrequency is not None and biasDir is None):
raise ValueError('Must specify both saveFrequency and biasDir')
if saveFrequency is not None and (saveFrequency < frequency or saveFrequency%frequency != 0):
raise ValueError('saveFrequency must be a multiple of frequency')
self.variables = variables
self.temperature = temperature
self.biasFactor = biasFactor
self.height = height
self.frequency = frequency
self.biasDir = biasDir
self.saveFrequency = saveFrequency
self._id = np.random.randint(0x7FFFFFFF)
self._saveIndex = 0
self._selfBias = np.zeros(tuple(v.gridWidth for v in variables))
self._totalBias = np.zeros(tuple(v.gridWidth for v in variables))
self._loadedBiases = {}
self._deltaT = temperature*(biasFactor-1)
varNames = ['cv%d' % i for i in range(len(variables))]
self._force = mm.CustomCVForce('table(%s)' % ', '.join(varNames))
for name, var in zip(varNames, variables):
self._force.addCollectiveVariable(name, var.force)
widths = [v.gridWidth for v in variables]
mins = [v.minValue for v in variables]
maxs = [v.maxValue for v in variables]
if len(variables) == 1:
self._table = mm.Continuous1DFunction(self._totalBias.flatten(), mins[0], maxs[0])
elif len(variables) == 2:
self._table = mm.Continuous2DFunction(widths[0], widths[1], self._totalBias.flatten(), mins[0], maxs[0], mins[1], maxs[1])
elif len(variables) == 3:
self._table = mm.Continuous3DFunction(widths[0], widths[1], widths[2], self._totalBias.flatten(), mins[0], maxs[0], mins[1], maxs[1], mins[2], maxs[2])
else:
raise ValueError('Metadynamics requires 1, 2, or 3 collective variables')
self._force.addTabulatedFunction('table', self._table)
self._force.setForceGroup(31)
system.addForce(self._force)
self._syncWithDisk()
def step(self, simulation, steps):
"""Advance the simulation by integrating a specified number of time steps.
Parameters
----------
simulation: Simulation
the Simulation to advance
steps: int
the number of time steps to integrate
"""
stepsToGo = steps
while stepsToGo > 0:
nextSteps = stepsToGo
if simulation.currentStep % self.frequency == 0:
nextSteps = min(nextSteps, self.frequency)
else:
nextSteps = min(nextSteps, simulation.currentStep % self.frequency)
simulation.step(nextSteps)
if simulation.currentStep % self.frequency == 0:
position = self._force.getCollectiveVariableValues(simulation.context)
energy = simulation.context.getState(getEnergy=True, groups={31}).getPotentialEnergy()
height = self.height*np.exp(-energy/(unit.MOLAR_GAS_CONSTANT_R*self._deltaT))
self._addGaussian(position, height, simulation.context)
if self.saveFrequency is not None and simulation.currentStep % self.saveFrequency == 0:
self._syncWithDisk()
stepsToGo -= nextSteps
def getFreeEnergy(self):
"""Get the free energy of the system as a function of the collective variables.
The result is returned as a N-dimensional NumPy array, where N is the number of collective
variables. The values are in kJ/mole. The i'th position along an axis corresponds to
minValue + i*(maxValue-minValue)/gridWidth.
"""
return -((self.temperature+self._deltaT)/self._deltaT)*self._totalBias*unit.kilojoules_per_mole
def getCollectiveVariables(self, simulation):
"""Get the current values of all collective variables in a Simulation."""
return self._force.getCollectiveVariableValues(simulation.context)
def _addGaussian(self, position, height, context):
"""Add a Gaussian to the bias function."""
# Compute a Gaussian along each axis.
axisGaussians = []
for i,v in enumerate(self.variables):
x = (position[i]-v.minValue) / (v.maxValue-v.minValue)
if v.periodic:
x = x % 1.0
dist = np.abs(np.linspace(0, 1.0, num=v.gridWidth) - x)
if v.periodic:
dist = np.min(np.array([dist, np.abs(dist-1)]), axis=0)
axisGaussians.append(np.exp(-dist*dist*v.gridWidth/v.biasWidth))
# Compute their outer product.
if len(self.variables) == 1:
gaussian = axisGaussians[0]
else:
gaussian = reduce(np.multiply.outer, reversed(axisGaussians))
# Add it to the bias.
height = height.value_in_unit(unit.kilojoules_per_mole)
self._selfBias += height*gaussian
self._totalBias += height*gaussian
widths = [v.gridWidth for v in self.variables]
mins = [v.minValue for v in self.variables]
maxs = [v.maxValue for v in self.variables]
if len(self.variables) == 1:
self._table.setFunctionParameters(self._totalBias.flatten(), mins[0], maxs[0])
elif len(self.variables) == 2:
self._table.setFunctionParameters(widths[0], widths[1], self._totalBias.flatten(), mins[0], maxs[0], mins[1], maxs[1])
elif len(self.variables) == 3:
self._table.setFunctionParameters(widths[0], widths[1], widths[2], self._totalBias.flatten(), mins[0], maxs[0], mins[1], maxs[1], mins[2], maxs[2])
self._force.updateParametersInContext(context)
def _syncWithDisk(self):
"""Save biases to disk, and check for updated files created by other processes."""
if self.biasDir is None:
return
# Use a safe save to write out the biases to disk, then delete the older file.
oldName = os.path.join(self.biasDir, 'bias_%d_%d.npy' % (self._id, self._saveIndex))
self._saveIndex += 1
tempName = os.path.join(self.biasDir, 'temp_%d_%d.npy' % (self._id, self._saveIndex))
fileName = os.path.join(self.biasDir, 'bias_%d_%d.npy' % (self._id, self._saveIndex))
np.save(tempName, self._selfBias)
os.rename(tempName, fileName)
if os.path.exists(oldName):
os.remove(oldName)
# Check for any files updated by other processes.
fileLoaded = False
pattern = re.compile('bias_(.*)_(.*)\.npy')
for filename in os.listdir(self.biasDir):
match = pattern.match(filename)
if match is not None:
matchId = int(match.group(1))
matchIndex = int(match.group(2))
if matchId != self._id and (matchId not in self._loadedBiases or matchIndex > self._loadedBiases[matchId].index):
try:
data = np.load(os.path.join(self.biasDir, filename))
self._loadedBiases[matchId] = _LoadedBias(matchId, matchIndex, data)
fileLoaded = True
except IOError:
# There's a tiny chance the file could get deleted by another process between when
# we check the directory and when we try to load it. If so, just ignore the error
# and keep using whatever version of that process' biases we last loaded.
pass
# If we loaded any files, recompute the total bias from all processes.
if fileLoaded:
self._totalBias = self._selfBias
for bias in self._loadedBiases.values():
self._totalBias += bias.bias
class BiasVariable(object):
"""A collective variable that can be used to bias a simulation with metadynamics."""
def __init__(self, force, minValue, maxValue, biasWidth, periodic=False, gridWidth=None):
"""Create a BiasVariable.
Parameters
----------
force: Force
the Force object whose potential energy defines the collective variable
minValue: float
the minimum value the collective variable can take. If it should ever go below this,
the bias force will be set to 0.
maxValue: float
the maximum value the collective variable can take. If it should ever go above this,
the bias force will be set to 0.
biasWidth: float
the width (standard deviation) of the Gaussians added to the bias during metadynamics
periodic: bool
whether this is a periodic variable, such that minValue and maxValue are physical equivalent
gridWidth: int
the number of grid points to use when tabulating the bias function. If this is omitted,
a reasonable value is chosen automatically.
"""
self.force = force
self.minValue = minValue
self.maxValue = maxValue
self.biasWidth = biasWidth
self.periodic = periodic
if gridWidth is None:
self.gridWidth = int(np.ceil(5*(maxValue-minValue)/biasWidth))
else:
self.gridWidth = gridWidth
_LoadedBias = namedtuple('LoadedBias', ['id', 'index', 'bias'])
wrappers/python/simtk/openmm/app/simulatedtempering.py 0 → 100644
View file @ 2e86184b
from __future__ import print_function
"""
simulatedtempering.py: Implements simulated tempering
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.
"""
__author__ = "Peter Eastman"
__version__ = "1.0"
import simtk.unit as unit
import math
import random
from sys import stdout
try:
import bz2
have_bz2 = True
except: have_bz2 = False
try:
import gzip
have_gzip = True
except: have_gzip = False
class SimulatedTempering(object):
"""SimulatedTempering implements the simulated tempering algorithm for accelerated sampling.
It runs a simulation while allowing the temperature to vary. At high temperatures, it can more easily cross
energy barriers to explore a wider area of conformation space. At low temperatures, it can thoroughly
explore each local region. For details, see Marinari, E. and Parisi, G., Europhys. Lett. 19(6). pp. 451-458 (1992).
The set of temperatures to sample can be specified in two ways. First, you can explicitly provide a list
of temperatures by using the "temperatures" argument. Alternatively, you can specify the minimum and
maximum temperatures, and the total number of temperatures to use. The temperatures are chosen spaced
exponentially between the two extremes. For example,
st = SimulatedTempering(simulation, numTemperatures=15, minTemperature=300*kelvin, maxTemperature=450*kelvin)
After creating the SimulatedTempering object, call step() on it to run the simulation.
Transitions between temperatures are performed at regular intervals, as specified by the "tempChangeInterval"
argument. For each transition, a new temperature is selected using the independence sampling method, as
described in Chodera, J. and Shirts, M., J. Chem. Phys. 135, 194110 (2011).
Simulated tempering requires a "weight factor" for each temperature. Ideally, these should be chosen so
the simulation spends equal time at every temperature. You can specify the list of weights to use with the
optional "weights" argument. If this is omitted, weights are selected automatically using the Wang-Landau
algorithm as described in Wang, F. and Landau, D. P., Phys. Rev. Lett. 86(10), pp. 2050-2053 (2001).
To properly analyze the results of the simulation, it is important to know the temperature and weight factors
at every point in time. The SimulatedTempering object functions as a reporter, writing this information
to a file or stdout at regular intervals (which should match the interval at which you save frames from the
simulation). You can specify the output file and reporting interval with the "reportFile" and "reportInterval"
arguments.
"""
def __init__(self, simulation, temperatures=None, numTemperatures=None, minTemperature=None, maxTemperature=None, weights=None, tempChangeInterval=25, reportInterval=1000, reportFile=stdout):
"""Create a new SimulatedTempering.
Parameters
----------
simulation: Simulation
The Simulation defining the System, Context, and Integrator to use
temperatures: list
The list of temperatures to use for tempering, in increasing order
numTemperatures: int
The number of temperatures to use for tempering. If temperatures is not None, this is ignored.
minTemperature: temperature
The minimum temperature to use for tempering. If temperatures is not None, this is ignored.
maxTemperature: temperature
The maximum temperature to use for tempering. If temperatures is not None, this is ignored.
weights: list
The weight factor for each temperature. If none, weights are selected automatically.
tempChangeInterval: int
The interval (in time steps) at which to attempt transitions between temperatures
reportInterval: int
The interval (in time steps) at which to write information to the report file
reportFile: string or file
The file to write reporting information to, specified as a file name or file object
"""
self.simulation = simulation
if temperatures is None:
if unit.is_quantity(minTemperature):
minTemperature = minTemperature.value_in_unit(unit.kelvin)
if unit.is_quantity(maxTemperature):
maxTemperature = maxTemperature.value_in_unit(unit.kelvin)
self.temperatures = [minTemperature*((float(maxTemperature)/minTemperature)**(i/float(numTemperatures-1))) for i in range(numTemperatures)]*unit.kelvin
else:
numTemperatures = len(temperatures)
self.temperatures = [(t.value_in_unit(unit.kelvin) if unit.is_quantity(t) else t)*unit.kelvin for t in temperatures]
if any(self.temperatures[i] >= self.temperatures[i+1] for i in range(numTemperatures-1)):
raise ValueError('The temperatures must be in strictly increasing order')
self.tempChangeInterval = tempChangeInterval
self.reportInterval = reportInterval
self.inverseTemperatures = [1.0/(unit.MOLAR_GAS_CONSTANT_R*t) for t in self.temperatures]
# If necessary, open the file we will write reports to.
self._openedFile = isinstance(reportFile, str)
if self._openedFile:
# Detect the desired compression scheme from the filename extension
# and open all files unbuffered
if reportFile.endswith('.gz'):
if not have_gzip:
raise RuntimeError("Cannot write .gz file because Python could not import gzip library")
self._out = gzip.GzipFile(fileobj=open(reportFile, 'wb', 0))
elif reportFile.endswith('.bz2'):
if not have_bz2:
raise RuntimeError("Cannot write .bz2 file because Python could not import bz2 library")
self._out = bz2.BZ2File(reportFile, 'w', 0)
else:
self._out = open(reportFile, 'w', 1)
else:
self._out = reportFile
# Initialize the weights.
if weights is None:
self._weights = [0.0]*numTemperatures
self._updateWeights = True
self._weightUpdateFactor = 1.0
self._histogram = [0]*numTemperatures
self._hasMadeTransition = False
else:
self._weights = weights
self._updateWeights = False
# Select the initial temperature.
self.currentTemperature = 0
self.simulation.integrator.setTemperature(self.temperatures[self.currentTemperature])
# Add a reporter to the simulation which will handle the updates and reports.
class STReporter(object):
def __init__(self, st):
self.st = st
def describeNextReport(self, simulation):
st = self.st
steps1 = st.tempChangeInterval - simulation.currentStep%st.tempChangeInterval
steps2 = st.reportInterval - simulation.currentStep%st.reportInterval
steps = min(steps1, steps2)
isUpdateAttempt = (steps1 == steps)
return (steps, False, isUpdateAttempt, False, isUpdateAttempt)
def report(self, simulation, state):
st = self.st
if simulation.currentStep%st.tempChangeInterval == 0:
st._attemptTemperatureChange(state)
if simulation.currentStep%st.reportInterval == 0:
st._writeReport()
simulation.reporters.append(STReporter(self))
# Write out the header line.
headers = ['Steps', 'Temperature (K)']
for t in self.temperatures:
headers.append('%gK Weight' % t.value_in_unit(unit.kelvin))
print('#"%s"' % ('"\t"').join(headers), file=self._out)
def __del__(self):
if self._openedFile:
self._out.close()
@property
def weights(self):
return [x-self._weights[0] for x in self._weights]
def step(self, steps):
"""Advance the simulation by integrating a specified number of time steps."""
self.simulation.step(steps)
def _attemptTemperatureChange(self, state):
"""Attempt to move to a different temperature."""
# Compute the probability for each temperature. This is done in log space to avoid overflow.
logProbability = [(self._weights[i]-self.inverseTemperatures[i]*state.getPotentialEnergy()) for i in range(len(self._weights))]
maxLogProb = max(logProbability)
offset = maxLogProb + math.log(sum(math.exp(x-maxLogProb) for x in logProbability))
probability = [math.exp(x-offset) for x in logProbability]
r = random.random()
for j in range(len(probability)):
if r < probability[j]:
if j != self.currentTemperature:
# Rescale the velocities.
scale = math.sqrt(self.temperatures[j]/self.temperatures[self.currentTemperature])
velocities = [v*scale for v in state.getVelocities().value_in_unit(unit.nanometers/unit.picoseconds)]
self.simulation.context.setVelocities(velocities)
# Select this temperature.
self._hasMadeTransition = True
self.currentTemperature = j
self.simulation.integrator.setTemperature(self.temperatures[j])
if self._updateWeights:
# Update the weight factors.
self._weights[j] -= self._weightUpdateFactor
self._histogram[j] += 1
minCounts = min(self._histogram)
if minCounts > 20 and minCounts >= 0.2*sum(self._histogram)/len(self._histogram):
# Reduce the weight update factor and reset the histogram.
self._weightUpdateFactor *= 0.5
self._histogram = [0]*len(self.temperatures)
self._weights = [x-self._weights[0] for x in self._weights]
elif not self._hasMadeTransition and probability[self.currentTemperature] > 0.99:
# Rapidly increase the weight update factor at the start of the simulation to find
# a reasonable starting value.
self._weightUpdateFactor *= 2.0
self._histogram = [0]*len(self.temperatures)
return
r -= probability[j]
def _writeReport(self):
"""Write out a line to the report."""
temperature = self.temperatures[self.currentTemperature].value_in_unit(unit.kelvin)
values = [temperature]+self.weights
print(('%d\t' % self.simulation.currentStep) + '\t'.join('%g' % v for v in values), file=self._out)
wrappers/python/tests/TestMetadynamics.py 0 → 100644
View file @ 2e86184b
import unittest
from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *
class TestMetadynamics(unittest.TestCase):
"""Test the Metadynamics class"""
def testHarmonicOscillator(self):
"""Test running metadynamics on a harmonic oscillator."""
system = System()
system.addParticle(1.0)
system.addParticle(1.0)
force = HarmonicBondForce()
force.addBond(0, 1, 1.0, 100000.0)
system.addForce(force)
cv = CustomBondForce('r')
cv.addBond(0, 1)
bias = BiasVariable(cv, 0.94, 1.06, 0.02)
meta = Metadynamics(system, [bias], 300*kelvin, 2.0, 5.0, 10)
integrator = LangevinIntegrator(300*kelvin, 10/picosecond, 0.001*picosecond)
topology = Topology()
chain = topology.addChain()
residue = topology.addResidue('H2', chain)
topology.addAtom('H1', element.hydrogen, residue)
topology.addAtom('H2', element.hydrogen, residue)
simulation = Simulation(topology, system, integrator, Platform.getPlatformByName('Reference'))
simulation.context.setPositions([Vec3(0, 0, 0), Vec3(1, 0, 0)])
meta.step(simulation, 200000)
fe = meta.getFreeEnergy()
center = bias.gridWidth//2
fe -= fe[center]
# Energies should be reasonably well converged over the central part of the range.
for i in range(center-3, center+4):
r = bias.minValue + i*(bias.maxValue-bias.minValue)/(bias.gridWidth-1)
e = 0.5*100000.0*(r-1.0)**2*kilojoules_per_mole
assert abs(fe[i]-e) < 1.0*kilojoules_per_mole
\ No newline at end of file
wrappers/python/tests/TestSimulatedTempering.py 0 → 100644
View file @ 2e86184b
import unittest
from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *
class TestSimulatedTempering(unittest.TestCase):
"""Test the SimulatedTempering class"""
def testHarmonicOscillator(self):
"""Test running simulated tempering on a harmonic oscillator."""
system = System()
system.addParticle(1.0)
system.addParticle(1.0)
force = HarmonicBondForce()
force.addBond(0, 1, 1.0, 1000.0)
system.addForce(force)
integrator = LangevinIntegrator(300*kelvin, 10/picosecond, 0.001*picosecond)
topology = Topology()
chain = topology.addChain()
residue = topology.addResidue('H2', chain)
topology.addAtom('H1', element.hydrogen, residue)
topology.addAtom('H2', element.hydrogen, residue)
simulation = Simulation(topology, system, integrator, Platform.getPlatformByName('Reference'))
st = SimulatedTempering(simulation, numTemperatures=10, minTemperature=200*kelvin, maxTemperature=400*kelvin, tempChangeInterval=5, reportInterval=10000)
self.assertEqual(10, len(st.temperatures))
self.assertEqual(200*kelvin, st.temperatures[0])
self.assertEqual(400*kelvin, st.temperatures[-1])
simulation.context.setPositions([Vec3(0, 0, 0), Vec3(1, 0, 0)])
# Run for a little while to let the weights stabilize.
st.step(10000)
# Run for a while and record the temperatures and distances.
distances = [[] for i in range(10)]
count = 0
for i in range(7000):
st.step(5)
pos = simulation.context.getState(getPositions=True).getPositions().value_in_unit(nanometers)
r = norm(pos[0]-pos[1])
distances[st.currentTemperature].append(r)
count += 1
# Check that it spent roughly equal time at each temperature, and that the distributions
# are correct.
for d, t in zip(distances, st.temperatures):
n = len(d)
assert count/20 < n < count/5
meanDist = sum(d)/n
assert 0.97 < meanDist < 1.03
meanEnergy = sum([0.5*1000*(r-1)**2 for r in d])/n
expectedEnergy = (0.5*MOLAR_GAS_CONSTANT_R*t).value_in_unit(kilojoules_per_mole)
self.assertAlmostEqual(expectedEnergy, meanEnergy, delta=expectedEnergy*0.3)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment