pyemma.thermo.estimate_umbrella_sampling

pyemma.thermo.estimate_umbrella_sampling(us_trajs, us_dtrajs, us_centers, us_force_constants, md_trajs=None, md_dtrajs=None, kT=None, maxiter=10000, maxerr=1e-15, save_convergence_info=0, estimator='wham', lag=1, dt_traj='1 step', init=None, init_maxiter=10000, init_maxerr=1e-08, **kwargs)

This function acts as a wrapper for tram(), dtram(), and wham() and handles the calculation of bias energies (bias) and thermodynamic state trajectories (ttrajs) when the data comes from umbrella sampling and (optional) unbiased simulations.

Parameters:
  • us_trajs (list of N arrays, each of shape (T_i, d)) – List of arrays, each having T_i rows, one for each time step, and d columns where d is the dimensionality of the subspace in which umbrella sampling was applied. Often d=1, and thus us_trajs will be a list of 1d-arrays.
  • us_dtrajs (list of N int arrays, each of shape (T_i,)) – The integers are indexes in 0,...,n-1 enumerating the n discrete states or the bins the umbrella sampling trajectory is in at any time.
  • us_centers (list of N floats or d-dimensional arrays of floats) – List or array of N center positions. Each position must be a d-dimensional vector. For 1d umbrella sampling, one can simply pass a list of centers, e.g. [-5.0, -4.0, -3.0, ... ].
  • us_force_constants (list of N floats or d- or dxd-dimensional arrays of floats) – The force constants used in the umbrellas, unit-less (e.g. kT per squared length unit). For multidimensional umbrella sampling, the force matrix must be used.
  • md_trajs (list of M arrays, each of shape (T_i, d), optional, default=None) – Unbiased molecular dynamics simulations; format like us_trajs.
  • md_dtrajs (list of M int arrays, each of shape (T_i,)) – The integers are indexes in 0,...,n-1 enumerating the n discrete states or the bins the unbiased trajectory is in at any time.
  • kT (float or None, optional, default=None) – Use this attribute if the supplied force constants are NOT unit-less; kT must have the same energy unit as the force constants.
  • maxiter (int, optional, default=10000) – The maximum number of self-consistent iterations before the estimator exits unsuccessfully.
  • maxerr (float, optional, default=1.0E-15) – Convergence criterion based on the maximal free energy change in a self-consistent iteration step.
  • save_convergence_info (int, optional, default=0) – Every save_convergence_info iteration steps, store the actual increment and the actual loglikelihood; 0 means no storage.
  • estimator (str, optional, default='wham') –

    Specify one of the available estimators

    ‘wham’: use WHAM
    ‘mbar’: use MBAR
    ‘dtram’: use the discrete version of TRAM
    ‘tram’: use TRAM
  • lag (int or list of int, optional, default=1) – Integer lag time at which transitions are counted. Providing a list of lag times will trigger one estimation per lag time.
  • dt_traj (str, optional, default='1 step') –

    Description of the physical time corresponding to the lag. May be used by analysis algorithms such as plotting tools to pretty-print the axes. By default ‘1 step’, i.e. there is no physical time unit. Specify by a number, whitespace and unit. Permitted units are (* is an arbitrary string):

    ‘fs’, ‘femtosecond*’
    ‘ps’, ‘picosecond*’
    ‘ns’, ‘nanosecond*’
    ‘us’, ‘microsecond*’
    ‘ms’, ‘millisecond*’
    ‘s’, ‘second*’
  • init (str, optional, default=None) –

    Use a specific initialization for the self-consistent iteration:

    None: use a hard-coded guess for free energies and Lagrangian multipliers
    ‘wham’: perform a short WHAM estimate to initialize the free energies (only with dtram)
    ‘mbar’: perform a short MBAR estimate to initialize the free energies (only with tram)
  • init_maxiter (int, optional, default=10000) – The maximum number of self-consistent iterations during the initialization.
  • init_maxerr (float, optional, default=1.0E-8) – Convergence criterion for the initialization.
  • **kwargs (dict, optional) – You can use this to pass estimator-specific named parameters to the chosen estimator, which are not already coverd by estimate_umbrella_sampling().
Returns:

estimator_obj – The requested estimator/model object, i.e., WHAM, MBAR, DTRAM or TRAM. If multiple lag times are given, a list of objects is returned (one MEMM per lag time).
Return type:

MEMM or MultiThermModel or list thereof

Example

We look at a 1D umbrella sampling simulation with two umbrellas at 1.1 and 1.3 on the reaction coordinate with spring constant of 1.0; additionally, we have two unbiased simulations.

We start with a joint clustering and use TRAM for the estimation:

>>> from pyemma.coordinates import cluster_regspace as regspace
>>> from pyemma.thermo import estimate_umbrella_sampling as estimate_us
>>> import numpy as np
>>> us_centers = [1.1, 1.3]
>>> us_force_constants = [1.0, 1.0]
>>> us_trajs = [np.array([1.0, 1.1, 1.2, 1.1, 1.0, 1.1]).reshape((-1, 1)), np.array([1.3, 1.2, 1.3, 1.4, 1.4, 1.3]).reshape((-1, 1))]
>>> md_trajs = [np.array([0.9, 1.0, 1.1, 1.2, 1.3, 1.4]).reshape((-1, 1)), np.array([1.5, 1.4, 1.3, 1.4, 1.4, 1.5]).reshape((-1, 1))]
>>> cluster = regspace(data=us_trajs+md_trajs, max_centers=10, dmin=0.15)
>>> us_dtrajs = cluster.dtrajs[:2]
>>> md_dtrajs = cluster.dtrajs[2:]
>>> centers = cluster.clustercenters
>>> tram = estimate_us(us_trajs, us_dtrajs, us_centers, us_force_constants, md_trajs=md_trajs, md_dtrajs=md_dtrajs, estimator='tram', lag=1)
>>> tram.f 
array([ 0.63...,  1.60...,  1.31...])