openmmtools.multistate.SAMSSampler

class openmmtools.multistate.SAMSSampler(number_of_iterations=1, log_target_probabilities=None, state_update_scheme='global-jump', locality=5, update_stages='two-stage', flatness_criteria='logZ-flatness', flatness_threshold=0.2, weight_update_method='rao-blackwellized', adapt_target_probabilities=False, gamma0=1.0, logZ_guess=None, **kwargs)

Self-adjusted mixture sampling (SAMS), also known as optimally-adjusted mixture sampling.

This class provides a facility for self-adjusted mixture sampling simulations. One or more replicas use the method of expanded ensembles [1] to sample multiple thermodynamic states within each replica, with log weights for each thermodynamic state adapted on the fly [2] to achieve the desired target probabilities for each state.

See also

ReplicaExchangeSampler

References

[1] Lyubartsev AP, Martsinovski AA, Shevkunov SV, and Vorontsov-Velyaminov PN. New approach to Monte Carlo calculation of the free energy: Method of expanded ensembles. JCP 96:1776, 1992 http://dx.doi.org/10.1063/1.462133

[2] Tan, Z. Optimally adjusted mixture sampling and locally weighted histogram analysis, Journal of Computational and Graphical Statistics 26:54, 2017. http://dx.doi.org/10.1080/10618600.2015.1113975

Examples

SAMS simulation of alanine dipeptide in implicit solvent at different temperatures.

Create the system:

>>> import math
>>> from openmm import unit
>>> from openmmtools import testsystems, states, mcmc
>>> testsystem = testsystems.AlanineDipeptideVacuum()
>>> import os
>>> import tempfile

Create thermodynamic states for parallel tempering with exponentially-spaced schedule:

>>> n_replicas = 3  # Number of temperature replicas.
>>> T_min = 298.0 * unit.kelvin  # Minimum temperature.
>>> T_max = 600.0 * unit.kelvin  # Maximum temperature.
>>> temperatures = [T_min + (T_max - T_min) * (math.exp(float(i) / float(n_replicas-1)) - 1.0) / (math.e - 1.0)
...                 for i in range(n_replicas)]
>>> thermodynamic_states = [states.ThermodynamicState(system=testsystem.system, temperature=T)
...                         for T in temperatures]

Initialize simulation object with options. Run with a GHMC integrator:

>>> move = mcmc.GHMCMove(timestep=2.0*unit.femtoseconds, n_steps=50)
>>> simulation = SAMSSampler(mcmc_moves=move, number_of_iterations=2,
...                          state_update_scheme='global-jump', locality=5,
...                          update_stages='two-stage', flatness_criteria='logZ-flatness',
...                          flatness_threshold=0.2, weight_update_method='rao-blackwellized',
...                          adapt_target_probabilities=False)

Create a single-replica SAMS simulation bound to a storage file and run:

>>> storage_path = tempfile.NamedTemporaryFile(delete=False).name + '.nc'
>>> reporter = multistate.MultiStateReporter(storage_path, checkpoint_interval=1)
>>> simulation.create(thermodynamic_states=thermodynamic_states,
...                   sampler_states=[states.SamplerState(testsystem.positions)],
...                   storage=reporter)
Please cite the following:

        Friedrichs MS, Eastman P, Vaidyanathan V, Houston M, LeGrand S, Beberg AL, Ensign DL, Bruns CM, and Pande VS. Accelerating molecular dynamic simulations on graphics processing unit. J. Comput. Chem. 30:864, 2009. DOI: 10.1002/jcc.21209
        Eastman P and Pande VS. OpenMM: A hardware-independent framework for molecular simulations. Comput. Sci. Eng. 12:34, 2010. DOI: 10.1109/MCSE.2010.27
        Eastman P and Pande VS. Efficient nonbonded interactions for molecular dynamics on a graphics processing unit. J. Comput. Chem. 31:1268, 2010. DOI: 10.1002/jcc.21413
        Eastman P and Pande VS. Constant constraint matrix approximation: A robust, parallelizable constraint method for molecular simulations. J. Chem. Theor. Comput. 6:434, 2010. DOI: 10.1021/ct900463w
>>> simulation.run()  # This runs for a maximum of 2 iterations.
>>> simulation.iteration
2
>>> simulation.run(n_iterations=1)
>>> simulation.iteration
2

To resume a simulation from an existing storage file and extend it beyond the original number of iterations.

>>> del simulation
>>> simulation = SAMSSampler.from_storage(reporter)
Please cite the following:

        Friedrichs MS, Eastman P, Vaidyanathan V, Houston M, LeGrand S, Beberg AL, Ensign DL, Bruns CM, and Pande VS. Accelerating molecular dynamic simulations on graphics processing unit. J. Comput. Chem. 30:864, 2009. DOI: 10.1002/jcc.21209
        Eastman P and Pande VS. OpenMM: A hardware-independent framework for molecular simulations. Comput. Sci. Eng. 12:34, 2010. DOI: 10.1109/MCSE.2010.27
        Eastman P and Pande VS. Efficient nonbonded interactions for molecular dynamics on a graphics processing unit. J. Comput. Chem. 31:1268, 2010. DOI: 10.1002/jcc.21413
        Eastman P and Pande VS. Constant constraint matrix approximation: A robust, parallelizable constraint method for molecular simulations. J. Chem. Theor. Comput. 6:434, 2010. DOI: 10.1021/ct900463w
>>> simulation.extend(n_iterations=1)
>>> simulation.iteration
3

You can extract several information from the NetCDF file using the Reporter class while the simulation is running. This reads the SamplerStates of every run iteration.

>>> reporter = multistate.MultiStateReporter(storage=storage_path, open_mode='r', checkpoint_interval=1)
>>> sampler_states = reporter.read_sampler_states(iteration=3)
>>> len(sampler_states)
1
>>> sampler_states[-1].positions.shape  # Alanine dipeptide has 22 atoms.
(22, 3)

Clean up.

>>> os.remove(storage_path)

Attributes:

log_target_probabilitiesarray-like

log_target_probabilities[state_index] is the log target probability for state state_index

state_update_schemestr

Thermodynamic state sampling scheme. One of [‘global-jump’, ‘local-jump’, ‘restricted-range’]

localityint

Number of neighboring states on either side to consider for local update schemes

update_stagesstr

Number of stages to use for update. One of [‘one-stage’, ‘two-stage’]

weight_update_methodstr

Method to use for updating log weights in SAMS. One of [‘optimal’, ‘rao-blackwellized’]

adapt_target_probabilitiesbool

If True, target probabilities will be adapted to achieve minimal thermodynamic length between terminal thermodynamic states.

gamma0float, optional, default=0.0

Initial weight adaptation rate.

logZ_guessarray-like of shape [n_states] of floats, optional, default=None

Initial guess for logZ for all states, if available.

Methods

Status(iteration, target_error, is_completed)	Attributes:
create(thermodynamic_states, sampler_states, ...)	Create new multistate sampler simulation.
default_options()	dict of all default class options (keyword arguments for __init__ for class and superclasses)
equilibrate(n_iterations[, mcmc_moves])	Equilibrate all replicas.
extend(n_iterations)	Extend the simulation by the given number of iterations.
from_storage(storage)	Constructor from an existing storage file.
minimize([tolerance, max_iterations])	Minimize all replicas.
read_status(storage)	Read the status of the calculation from the storage file.
run([n_iterations])	Run the replica-exchange simulation.

__init__(number_of_iterations=1, log_target_probabilities=None, state_update_scheme='global-jump', locality=5, update_stages='two-stage', flatness_criteria='logZ-flatness', flatness_threshold=0.2, weight_update_method='rao-blackwellized', adapt_target_probabilities=False, gamma0=1.0, logZ_guess=None, **kwargs)

Initialize a SAMS sampler.

Parameters:

log_target_probabilitiesarray-like or None

log_target_probabilities[state_index] is the log target probability for thermodynamic state state_index When converged, each state should be sampled with the specified log probability. If None, uniform probabilities for all states will be assumed.

state_update_schemestr, optional, default=’global-jump’

Specifies the scheme used to sample new thermodynamic states given fixed sampler states. One of [‘global-jump’, ‘local-jump’, ‘restricted-range-jump’] global_jump will allow the sampler to access any thermodynamic state local-jump will propose a move to one of the local neighborhood states, and accept or reject. restricted-range will compute the probabilities for each of the states in the local neighborhood, increasing jump probability

localityint, optional, default=1

Number of neighboring states on either side to consider for local update schemes.

update_stagesstr, optional, default=’two-stage’

One of [‘one-stage’, ‘two-stage’] one-stage will use the asymptotically optimal scheme throughout the entire simulation (not recommended due to slow convergence) two-stage will use a heuristic first stage to achieve flat histograms before switching to the asymptotically optimal scheme

flatness_criteriastring, optiona, default=’logZ-flatness’

Method of assessing when to switch to asymptotically optimal scheme

One of [‘logZ-flatness’,’minimum-visits’,’histogram-flatness’]

flatness_thresholdfloat, optional, default=0.2

Histogram relative flatness threshold to use for first stage of two-stage scheme.

weight_update_methodstr, optional, default=’rao-blackwellized’

Method to use for updating log weights in SAMS. One of [‘optimal’, ‘rao-blackwellized’] rao-blackwellized will update log free energy estimate for all states for which energies were computed optimal will use integral counts to update log free energy estimate of current state only

adapt_target_probabilitiesbool, optional, default=False

If True, target probabilities will be adapted to achieve minimal thermodynamic length between terminal thermodynamic states. (EXPERIMENTAL)

gamma0float, optional, default=0.0

Initial weight adaptation rate.

logZ_guessarray-like of shape [n_states] of floats, optiona, default=None

Initial guess for logZ for all states, if available.

Methods

__init__([number_of_iterations, ...])	Initialize a SAMS sampler.
create(thermodynamic_states, sampler_states, ...)	Create new multistate sampler simulation.
default_options()	dict of all default class options (keyword arguments for __init__ for class and superclasses)
equilibrate(n_iterations[, mcmc_moves])	Equilibrate all replicas.
extend(n_iterations)	Extend the simulation by the given number of iterations.
from_storage(storage)	Constructor from an existing storage file.
minimize([tolerance, max_iterations])	Minimize all replicas.
read_status(storage)	Read the status of the calculation from the storage file.
run([n_iterations])	Run the replica-exchange simulation.

Attributes

is_completed	Check if we have reached any of the stop target criteria (read-only)
is_periodic	Return True if system is periodic, False if not, and None if not initialized
iteration	The integer current iteration of the simulation (read-only).
mcmc_moves	A copy of the MCMCMoves list used to propagate the simulation.
metadata	A copy of the metadata dictionary passed on creation (read-only).
n_replicas	The integer number of replicas (read-only).
n_states	The integer number of thermodynamic states (read-only).
options	dict of all class options (keyword arguments for __init__ for class and superclasses)
sampler_states	A copy of the sampler states list at the current iteration.
online_analysis_interval	interval to carry out online analysis