Implementation in ORAC 

From the practical point of view, a metadynamics simulation consists in two steps. In the first one, a set of reaction coordinates is chosen whose dynamics describes the process under study. As we said, such a procedure requires an high degree of chemical and physical intuition for its application to complex molecular system, since these variables are not obviously determined from a molecular structure.

The second step is the metadynamics simulation itself, during which an history-dependent potential is constructed by summing, at regular time intervals, repulsive potential terms centered in the current position of the system in the space of the reaction coordinates. In its standard implementation, the history-dependent potential of metadynamics is given by a sum of small repulsive Gaussian, Eq.7.2. Some variants have been introduced, with the intent of improving the accuracy or the efficiency of the method[136,137]. In the ORAC program we have used Lucy's function[138] as a very efficient alternative to the use of Gaussian distributions.. It is defined as ^7.1

$${\cal L}(s;s_0,h,w) = h \left ( 1 + 2 \frac{\vert s-s_0\vert}{w} ... ...1 - \frac{\vert s-s_0\vert}{w} \right)^2; 0 ~ \textrm{if } \vert s-s_0\vert > w$$ (7.3)

with the origin at $s_0$. The symbols $h$ and $w$ denote the height and the width, respectively. Such a function is normalizable, $\int_{-\infty}^{\infty} {\rm d}s ~ {\cal L}(s;s_0,w) = h w$, has a finite range $w$, has a maximum at the origin and it is differentiable everywhere. A Lucy's function can be compared with a Gaussian function with the same value at the origin and at $\vert s\vert = s_0 + w/2$, such that

$$2\sigma = w / ( 2 \ln 2 )^{1/2}$$ (7.4)

A Lucy's function can be regarded as a Gaussian function with $\sigma$ in Eq.7.4, but without the long tails of the Gaussian, as can be seen in Fig.7.1 where a Lucy's function with $h=w=1$ and a Gaussian function with the same height and $\sigma = w / 2(2 \ln2)^{1/2}$ are shown. The parameters $h$, $w$ and $\tau$ affects the accuracy of the free energy reconstruction in a similar manner to the height and the width of Gaussian functions and a comprehensive review on the analysis of the error during a metadynamics run can be found in [62].

Figure 7.1: Lucy's function $L$ with $h=w=1$, along with a Gaussian function $G$ with the same height and $2\sigma = w/(2 \ln 2)^{1/2}$.

The history dependent potential used during an ORAC simulation can therefore be written as

$$V(s,t) = \sum_{t^{\prime} = \tau, 2\tau, ...} {\cal L}(s;s_{t^{\prime}},h,w)$$ (7.5)

During a simulation, forces from this biasing potential are computed in the shell $n1$ as a sum of derivatives of ${\cal L}$ functions:

$$\frac{\partial {\cal L}(s;s_0,h,w)}{\partial s} = \frac{6h}{w^3} (s-s_0) (\vert s-s_0\vert - w); 0 ~ \textrm{if } \vert s-s_0\vert > w$$ (7.6)

Such a derivative is computationally attractive, since it does not require the evaluation of an exponential function as in the case of the derivative of a Gaussian function. Moreover, since ${\cal L}$ has a finite range by definition, it does not need to be smoothly truncated[137], as there are no contribution to the forces from hills farther than the width $w$.

Using the standard metadynamics approach, during a simulation the algorithm keeps on adding terms to the history-dependent potential (the sum in Eq.7.5) with the same constant rate $\omega = h / \tau$. However, the optimal solution would be to use a faster rate at the beginning of the simulation, so as to produce a rough estimate of the free energy, and then to reduce $\omega$ to refine this estimate[140]. This problem corresponds to finding an optimal protocol for the evolution of the modification factor in the original Wang-Landau algorithm. Various solutions have been proposed[141,133,142,143] in which the energy $h$ in 7.3 is time-dependent. We propose instead to add a term to the biasing potential with a given probability $P_t({\rm add})$, depending parametrically on time. For example, for $P_t({\rm add}) \propto 1/t$, the evolution of the rate would be given by $\omega(t) = P_t({\rm add}) \omega_0 \propto \omega_0/t$. This procedure can be seen on average as an increasing deposition interval $\tau(t)$, such that $\omega(t) = h / \tau(t)$ decreases in time. In the present implementation of ORAC, three different choices are available for the probability $P({\rm add})$: the default one is simply $P({\rm add}) = 1$ and corresponds to the standard metadynamics algorithm. The second one is given by

$$P_t({\rm add}) = e^{-V_{\rm max}(t)/k_BT^{\prime}}$$ (7.7)

where $V_{\rm max}(t)$ is the maximum value of the potential $V(s,t)$ at time $t$. During the simulation, the effective rate $\omega(t)$ decreases as $V_{\rm max}(t)$ increases. As $V_{\rm max} \gg k_B T^{\prime}$, the deposition rate $\omega(t)$ is so slow that the transformation can be considered adiabatic, and the biasing potential converges to the free energy inverted in sign, $A(s) = -V(s,t)$. The slowdown of $\omega$ can be tuned through the parameter $T^{\prime}$. Finally, following the well-tempered metadynamics approach[143], the third choice is given by

$$P_{s,t}({\rm add}) = e^{-V(s,t)/k_BT^{\prime}}$$ (7.8)

where the probability depends parametrically both on time $t$ and on position $s$ of the system along the reaction coordinate through the biasing potential $V$. In this case, the biasing potential does not converge to the free energy inverted in sign as in the previous case, since in general $\omega$ turns out to be coordinate-dependent even when the potential has flatten the free energy profile. However, as shown in [143], the relation

$$A(s) = - \frac{T + T^{\prime}}{T} V(s,t)$$ (7.9)

can be used to recover the original free energy from the biasing potential.

The multiple walkers version of metadynamics algorithm[144] was implemented in the parallel version of the code through the MPI library. This approach is based on running simultaneously multiple replicas of the system, contributing equally to the same history-dependent potential, and therefore to the same free energy surface reconstruction. For $N$ replicas, $V(s,t)$ can be written as a double sum

$$V(s,t) = \sum_{t^{\prime} = \tau, 2\tau,...t} \sum_{i=1,N} {\cal L}(s;s_{i,t^{\prime}},h,\sigma)$$ (7.10)

where $s_{i,t^{\prime}}$ is the position at time $t^{\prime}$ of the $i$-th replica along $s$. In particular, the enhanced efficiency of this algorithm with respect to uncoupled simulations contributes to make the calculation of FESs in high dimensions more accessible.

In the ORAC distribution at http://www.chim.unifi.it/orac we provide some example of metadynamics simulations using Lucy's functions on multi-dimensional surfaces of simple molecules in the gas phase along with some ancillary codes for the analysis of the program output.