Skip to content

GitLab

Commit 9dd49614 authored by Charlles Abreu's avatar Charlles Abreu
Browse files

Tested changes moved from WellTemperedMetadynamics to Metadynamics 

- temporary class WellTemperedMetadynamics deleted
parent 119bcf32
Hide whitespace changes
Inline Side-by-side
Showing with 67 additions and 195 deletions
wrappers/python/simtk/openmm/app/__init__.py
View file @ 9dd49614
......@@ -33,7 +33,7 @@ from .charmmcrdfiles import CharmmCrdFile, CharmmRstFile
from .charmmparameterset import CharmmParameterSet
from .charmmpsffile import CharmmPsfFile, CharmmPSFWarning
from .simulatedtempering import SimulatedTempering
from .metadynamics import Metadynamics, WellTemperedMetadynamics, BiasVariable
from .metadynamics import Metadynamics, BiasVariable
# Enumerated values
......
wrappers/python/simtk/openmm/app/metadynamics.py
View file @ 9dd49614
......@@ -72,7 +72,7 @@ class Metadynamics(object):
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):
def __init__(self, system, variables, temperature, biasFactor, height, frequency, saveFrequency=None, biasDir=None, gridExpansion=20):
"""Create a Metadynamics object.
Parameters
......@@ -100,6 +100,9 @@ class Metadynamics(object):
biasDir: str (optional)
the directory to which biases should be written, and from which biases written by
other processes should be loaded
gridExpansion: int (optional)
the number of extra grid points to be used in periodic directions of multidimensional
tabulated functions. This aims at avoiding boundary discontinuity artifacts.
"""
if not unit.is_quantity(temperature):
temperature = temperature*unit.kelvin
......@@ -120,19 +123,27 @@ class Metadynamics(object):
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))
for v in variables:
v._expanded = v.periodic and len(variables) > 1
v._extraWidth = min(gridExpansion, v.gridWidth - 1) if v._expanded else 0
extraRange = v._extraWidth*(v.maxValue - v.minValue)/(v.gridWidth - 1)
v._actualWidth = v.gridWidth + 2*v._extraWidth
v._actualMin = v.minValue - extraRange
v._actualMax = v.maxValue + extraRange
v._slice = slice(v._extraWidth, v.gridWidth + v._extraWidth)
self._selfBias = np.zeros(tuple(v._actualWidth for v in reversed(variables)))
self._totalBias = np.zeros(tuple(v._actualWidth for v in reversed(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]
widths = [v._actualWidth for v in variables]
mins = [v._actualMin for v in variables]
maxs = [v._actualMax for v in variables]
if len(variables) == 1:
self._table = mm.Continuous1DFunction(self._totalBias.flatten(), mins[0], maxs[0])
self._table = mm.Continuous1DFunction(self._totalBias.flatten(), mins[0], maxs[0], variables[0].periodic)
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:
......@@ -140,7 +151,8 @@ class Metadynamics(object):
else:
raise ValueError('Metadynamics requires 1, 2, or 3 collective variables')
self._force.addTabulatedFunction('table', self._table)
self._force.setForceGroup(31)
freeGroups = set(range(32)) - set(force.getForceGroup() for force in system.getForces())
self._force.setForceGroup(max(freeGroups))
system.addForce(self._force)
self._syncWithDisk()
......@@ -178,7 +190,12 @@ class Metadynamics(object):
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
f = -((self.temperature+self._deltaT)/self._deltaT)*self._totalBias*unit.kilojoules_per_mole
if len(self.variables) == 1:
return f
else:
s = [v._slice for v in self.variables]
return f[s[1], s[0]] if len(self.variables) == 2 else f[s[2], s[1], s[0]]
def getCollectiveVariables(self, simulation):
"""Get the current values of all collective variables in a Simulation."""
......@@ -196,7 +213,11 @@ class Metadynamics(object):
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))
values = np.exp(-0.5*dist*dist/v._scaledVariance)
if v._expanded:
n = v._extraWidth + 1
values = np.hstack((values[-n:-1], values, values[1:n]))
axisGaussians.append(values)
# Compute their outer product.
......@@ -210,10 +231,11 @@ class Metadynamics(object):
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]
widths = [v._actualWidth for v in self.variables]
mins = [v._actualMin for v in self.variables]
maxs = [v._actualMax for v in self.variables]
if len(self.variables) == 1:
self._totalBias[-1] = self._totalBias[0]
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])
......@@ -265,156 +287,6 @@ class Metadynamics(object):
self._totalBias += bias.bias
class WellTemperedMetadynamics(Metadynamics):
"""
Temporary class.
"""
def __init__(self, system, variables, temperature, biasFactor, height, frequency, saveFrequency=None, biasDir=None, gridExpansion=20):
"""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
gridExpansion: int (optional)
the extra number of grid points used in periodic directions for multidimensional
tabulated functions
"""
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
for v in variables:
v._expanded = v.periodic and len(variables) > 1
v._extraWidth = min(gridExpansion, v.gridWidth - 1) if v._expanded else 0
extraRange = v._extraWidth*(v.maxValue - v.minValue)/(v.gridWidth - 1)
v._actualWidth = v.gridWidth + 2*v._extraWidth
v._actualMin = v.minValue - extraRange
v._actualMax = v.maxValue + extraRange
v._slice = slice(v._extraWidth, v.gridWidth + v._extraWidth)
self._selfBias = np.zeros(tuple(v._actualWidth for v in reversed(variables)))
self._totalBias = np.zeros(tuple(v._actualWidth for v in reversed(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._actualWidth for v in variables]
mins = [v._actualMin for v in variables]
maxs = [v._actualMax for v in variables]
if len(variables) == 1:
self._table = mm.Continuous1DFunction(self._totalBias.flatten(), mins[0], maxs[0], variables[0].periodic)
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)
freeGroups = set(range(32)) - set(force.getForceGroup() for force in system.getForces())
self._force.setForceGroup(max(freeGroups))
system.addForce(self._force)
self._syncWithDisk()
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.
"""
f = -((self.temperature+self._deltaT)/self._deltaT)*self._totalBias*unit.kilojoules_per_mole
if len(self.variables) == 1:
return f
else:
s = [v._slice for v in self.variables]
if len(self.variables) == 2:
return f[s[1], s[0]]
else:
return f[s[2], s[1], s[0]]
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)
values = np.exp(-0.5*dist*dist/v._scaledVariance)
if v._expanded:
n = v._extraWidth + 1
values = np.hstack((values[-n:-1], values, values[1:n]))
axisGaussians.append(values)
# 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._actualWidth for v in self.variables]
mins = [v._actualMin for v in self.variables]
maxs = [v._actualMax for v in self.variables]
if len(self.variables) == 1:
self._totalBias[-1] = self._totalBias[0]
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)
class BiasVariable(object):
"""A collective variable that can be used to bias a simulation with metadynamics."""
......@@ -425,17 +297,17 @@ class BiasVariable(object):
----------
force: Force
the Force object whose potential energy defines the collective variable
minValue: float
minValue: float or unit.Quantity
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
maxValue: float or unit.Quantity
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
biasWidth: float or unit.Quantity
the width (standard deviation) of the Gaussians added to the bias during metadynamics
periodic: bool
periodic: bool (optional)
whether this is a periodic variable, such that minValue and maxValue are physical equivalent
gridWidth: int
gridWidth: int (optional)
the number of grid points to use when tabulating the bias function. If this is omitted,
a reasonable value is chosen automatically.
"""
......
wrappers/python/simtk/openmm/app/simulatedtempering.py
View file @ 9dd49614
......@@ -12,7 +12,7 @@ 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
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,
......@@ -55,29 +55,29 @@ except: have_numpy = 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
......@@ -87,7 +87,7 @@ class SimulatedTempering(object):
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
......@@ -108,7 +108,7 @@ class SimulatedTempering(object):
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):
......@@ -143,9 +143,9 @@ class SimulatedTempering(object):
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
......@@ -157,12 +157,12 @@ class SimulatedTempering(object):
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
......@@ -181,11 +181,11 @@ class SimulatedTempering(object):
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))
......@@ -194,7 +194,7 @@ class SimulatedTempering(object):
def __del__(self):
if self._openedFile:
self._out.close()
@property
def weights(self):
return [x-self._weights[0] for x in self._weights]
......@@ -202,10 +202,10 @@ class SimulatedTempering(object):
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))]
......@@ -217,7 +217,7 @@ class SimulatedTempering(object):
if r < probability[j]:
if j != self.currentTemperature:
# Rescale the velocities.
scale = math.sqrt(self.temperatures[j]/self.temperatures[self.currentTemperature])
if have_numpy:
velocities = scale*state.getVelocities(asNumpy=True).value_in_unit(unit.nanometers/unit.picoseconds)
......@@ -226,26 +226,26 @@ class SimulatedTempering(object):
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
......
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