Metadynamics Simulation: history-dependent algorithms in Non-Boltzmann sampling

If we are studying a prototypical elementary reaction, in which two stable states are separated by a high free energy barrier $\Delta A^* \gg k_B T$ along the reaction coordinate $s$, configurations corresponding to the free energy maximum (the transition state $s^*$) can be sampled by adding a restraining potential to the original Hamiltonian of the system, so as to obtain a frequency histogram for the value of the reaction coordinate $s$ centered around the transition state itself. If we were good enough in locating the transition state and matching the curvature of the potential, this distribution will overlap with the two distributions obtained starting two different simulations from the two metastable states. The free energy difference between the metastable states, as well as the height of the free energy barrier at the transition state, can then be computed using the sampling from this ``bridging'' distribution. This solution is known as Umbrella Sampling[57]. More generally, if the transition state can be identified and located at some value of the reaction coordinate, the procedure of modifying the energetics of a system in order to balance the activation barrier and flatten the free energy profile is known with the name of Non-Boltzmann sampling. The original free energy can be computed from the free energy of the modified ensemble through the formula

$$A(s) \sim A^{\prime}(s) - V(s)$$ (7.1)

where $A^{\prime}(s)$ denotes the free energy computed by simulating the modified ergodic system. As in the Umbrella Sampling algorithm, the hardest part of the Non-Boltzmann sampling approach is the construction of a good biasing potential, since this task can be performed only iteratively. Given a rough (because of some free energy barrier) estimate of $A(s)$ from an old simulation, the simplest way to know how good this estimate is consists in performing a new simulation using this estimate, inverted in sign, as a bias potential. If the free energy profile of the modified system is flat, $A^{\prime} = {\rm constant}$, then $A(s) \sim - V(s)$ is the free energy inverted in sign. Otherwise, from this simulation we can compute an improved estimate for $A(s)$ through Eq. 7.1. The effectiveness of this tedious approach is due to the fact that each correction to the biasing potential makes the system more ergodic, and therefore each successive simulation is statistically more accurate than the former.

This iterative approach to the problem[130,131] led to the development of adaptive biasing potential methods that improve the potential ``on the fly'' [132,60,58,133], i.e., while the simulation is performed. All these methods share all the common basic idea, namely, ``to introduce the concept of memory''[132] during a simulation by changing the potential of mean force perceived by the system, in order to penalize conformations that have been already sampled before. The potential becomes history-dependent since it is now a functional of the past trajectory along the reaction coordinate. Among these algorithms, the Wang-Landau [60] and the metadynamics[58] algorithms have received most attention in the fields of the Monte Carlo (MC) and Molecular Dynamics (MD) simulations, respectively. This success is mainly due to the clearness and the ease of implementation of the algorithm, that is basically the same for the two methods. The Wang-Landau algorithm was initially proposed as a method to compute the density of states $g(E)$, and therefore the entropy $S(E) = \ln g(E)$, of a simulated discrete system. During a Wang-Landau MC simulation, $S(E)$ is estimated as an histogram, increasing by a fixed quantity the frequency of the visited energy levels, while moves are generated randomly and accepted with a Metropolis probability ${\rm acc}(E \rightarrow E^{\prime}) = {\rm min} \left \{ 1, \exp (- \Delta S ) \right \}$, where $\Delta S = S(E^{\prime}) - S(E)$ is the current estimate of the entropy change after the move. While for a random walk in energy the system would have been trapped in entropy maxima, the algorithm, that can be easily extended to the computation of any entropy-related thermodynamic potential along a generic collective variable, helps the system in escaping from these maxima and reconstructs the entropy $S(E)$. The metadynamics algorithm extends this approach to off-lattice systems and to Molecular Dynamics. Metadynamics has been successfully applied in the computation of free energy profiles in disparate fields, ranging from chemical physics to biophysics and material sciences. For a system in the canonical ensemble, metadynamics reconstructs the free energy along some reaction coordinate $s$ as a sum of Gaussian functions deposed along the trajectory of the system. This sum inverted in sign is used during the simulation as a biasing potential $V(s,t)$ that depends explicitly on time $s$:

$$V(s,t) = \sum_{t^{\prime} = \tau, 2\tau,...t} {\cal G}(s;s_{t^{\prime}},h,\sigma)$$ (7.2)

where ${\cal G}(s;s_t,h,\sigma) = h \exp \left (-(s-s_t)^2/2\sigma^2 \right )$ is a Gaussian function centered in $s_t$ with height $h$ and variance $\sigma^2$. During a metadynamics simulation, the potential $V(s,t)$ will grow faster for states with an higher probability, pushing out the system from minima in the free energy landscape. If the rate of deposition $\omega = h / \tau$ is sufficiently slow, the system can be considered in equilibrium with the biased Hamiltonian $H^{\prime}(x,t) = H(x) + V(s,t)$, and therefore the probability of visiting state $s$ at time $t$ is the equilibrium canonical distribution $p(s,t) \propto \exp [-\beta(A(s)+V(s,t)]$. Once all the free energy minima have been ``filled'' by the biasing potential, and therefore $V(s,t) = - A(s)$, such a probability is uniform along $s$ and the potential will grow uniformly.

The thermodynamic work spent in changing the potential from the original Hamiltonian $H(x)$ to $H^{\prime}(x,t) = H(x) + V(s,t)$ can be computed through the relation $W = \int_0^t {\rm d} \tau \left ( \frac{\partial H}{\partial t} \right )_\tau$. In the limit of an adiabatic transformation, this quantity is equal to the Helmholtz free energy difference $\Delta A = A^{\prime} - A_0$ between two systems with energy functions $H^{\prime}$ and $H$, where $A^{\prime} = \int dx \exp(-\beta H^{\prime})$ and $A_0 = \int dx \exp(-\beta H)$[134]. However, if the process is too fast with respect to the ergodic time scale, a part of the work spent during the switching will be dissipated in the system, resulting in an non-equilibrium, non-canonical distribution, and in a systematic error in the free energy estimate. In particular, it is assumed that during a metadynamics simulation all the microscopic variables different from the macroscopic reaction coordinate $s$ are always in the equilibrium state corresponding to the value of $s$[135]. This property is known with the name of Markov property, and it summarizes the main assumption of the algorithm: all the slow modes of the system coupled to the reaction under study have to be known a priori and they have to be included in the number of the reaction coordinates. Therefore, at variance with the methods presented in the previous chapters, metadynamics should be considered a quasi-equilibrium method, in which the knowledge about the variables that capture the mechanism of a reaction is exploited to gain insight on the transition states and more generally to compute the free energy landscape along the relevant reaction coordinates.