---
abstract: 'We present a method to sample reactive pathways via biased molecular dynamics simulations in trajectory space. We show that the use of enhanced sampling techniques enables unconstrained exploration of multiple reaction routes. Time correlation functions are conveniently computed via reweighted averages along a single trajectory and kinetic rates are accessed at no additional cost. These abilities are illustrated analyzing a model potential and the umbrella inversion of NH$_3$. The algorithm allows a parallel implementation and promises to be a powerful tool for the study of rare events.'
author:
- 'Davide Mandelli$^{a}$, Barak Hirshberg$^{b,c}$, Michele Parrinello$^{a,b,c,*}$'
nocite: '[@Salvalaglio2014]'
title: Metadynamics of paths
---

Molecular dynamics (MD) simulations have become an invaluable tool in numerous branches of science. While experiments generally access only spatially and time averaged quantities, atomically detailed MD simulations allow tracking in real time the microscopic mechanisms underlying complex phenomena. Nevertheless, there is a large class of problems where a straightforward application of MD simulations is impractical. Important examples are crystal nucleation, slow diffusion in solids, chemical reactions and conformational changes of large molecules. In all these cases, the presence of large free energy barriers leads to extremely long computational times. Therefore, it is necessary to design efficient algorithms able to accelerate phase space exploration.

In biased MD simulations, external potentials are used to drive transitions between meta-stable states. Among these methods, metadynamics [@Laio2002] (MetaD) has proven to be particularly successful. Sampling efficiency comes at the cost of introducing unphysical forces. Modifying the potential energy surface affects the dynamics of the system and information about the correct (unbiased) time evolution is lost. Methods have been developed aiming at retrieving dynamical information from biased MD simulations [@Tiwary2013; @Wu2014; @Rosta2015; @Wu2016]. In a recent effort, Donati and coworkers [@Donati2017] introduced a scheme to compute time correlation functions from biased Langevin trajectories. The main advantage of this technique is that it applies also to time dependent perturbations [@Donati2018]. However, applications are limited to cases where significant speed-up of phase space exploration can be achieved while keeping the external potential within a range compatible with the reweighting procedure, which otherwise leads to numerical instabilities.

Besides biased MD simulations, there exists an alternative class of techniques that rely on the definition of a functional on the space of all trajectories and tackle the problem of rare events directly in the space of paths [@Elber1987; @Olender1996; @Passerone2001; @Lee2017; @Pratt1986; @Dellago1998a; @Fujisaki2010]. These include the transition path sampling (TPS) algorithm of Chandler and collaborators [@Dellago1998a]. TPS is a Monte Carlo procedure that harvests dynamical RPs connecting two given states. Two practical difficulties arise in its application. First, the computation of rate constants is time consuming [@Dellago1998a; @VanErp2003; @Moroni2004]. Secondly, this approach poses problems in systems with multiple channels, where the presence of energy barriers can prevent exploration of all the relevant reaction routes [@Vlugt2001; @Borrero2016].

In this Letter, we explore the possibility of combining biased MD simulations with action-based algorithms for the sampling of RPs. Specifically, we apply MetaD to MD simulations performed directly in the space of paths. We show that MetaD can be used to harvest RPs without the need of imposing endpoint constraints, allowing exploration of multiple reaction routes in the same run. Furthermore, standard reweighting procedures can be applied straightforwardly and provide a convenient way of computing rate constants and time correlation functions. In the following, we briefly review the theory and formalism behind the algorithm and then present two applications.

The problem of interest here is the time evolution of a molecular system coupled to a thermal bath. The latter acts by introducing random fluctuations in the forces governing the dynamics. In this context, Onsager and Machlup [@Onsager1953] derived an expression for the probability of observing a given path $q(t)$ starting at $q_0$ at $t=0$ and ending in $q_{\tilde t}$ at $t$=$\tilde t$ $$\label{eq.OM}
P[q(t);q_0,q_{\tilde t}]\propto e^{ - \int_{0}^{\tilde t} {\mathcal L} (q(t),\dot{q}(t)){\rm d}t } = e^{-{\mathcal S}[q(t)]}.$$ Here, ${\mathcal L}$ and $\mathcal S$ define the Onsager-Machlup (OM) function and associated OM action. Similar in spirit to Hamilton’s principle of least action, equation provides a route to construct the most probable trajectory by minimizing ${\mathcal S}$ under given endpoints constraints. This observation motivated previous methods that solved the boundary valued problem adopting suitably discretized approximations of $\mathcal S$ [@Elber1987; @Olender1996; @Passerone2001; @Lee2017]. Alternatively, one notices that the integral $$\label{eq.Z}
{\mathcal Z}_{\tilde t} = \int {\mathcal D}{q(t)} e^{-{\mathcal S}[q(t)]},$$ over all paths of given duration, defines a partition function in the space of trajectories so that each of them can be assigned a well-defined probability $$\label{eq.P}
P_{\tilde t}[q(t)] = e^{-{\mathcal S}[q(t)]}/{\mathcal Z}_{\tilde t}.$$ Equation is the starting point of statistical path-sampling techniques [@Pratt1986; @Dellago1998a; @Fujisaki2010] and is central to the present work.

\[fig1.2Dpot\] ![(a) The double-well potential used in simulations. (b) A polymer configuration corresponding to a RP crossing the lower saddle. Green and blue circles indicate the initial and final position of the particle along the path.](fig_2Dpot.pdf "fig:"){width="0.8\columnwidth"}

We will now specialize equation to the case of a system evolving according to Brownian motion discretized as $$\label{eq.BD}
R^{n+1}_\alpha=R_\alpha^n+\frac{\Delta t}{m_\alpha\nu}F^n_\alpha+\delta r_{\alpha},$$ where ${\bf R}$ is the $D$-dimensional coordinate vector, $F_\alpha^n=-\nabla_{R_\alpha^n}U$ is the force acting on the $\alpha$-[*th*]{} degree of freedom ($\alpha$=1–$D$) at step $n$, $m_\alpha$ is its mass, $\nu$ is the damping coefficient, $\Delta t$ is the time step and $\delta r_\alpha$ is a normally distributed random displacement with variance $\sigma_\alpha^2=2k_BT\Delta t/m_\alpha\nu$ and zero mean. The probability of observing a given trajectory $\{{\bf R}^1 \rightarrow {\bf R}^2 \dots \rightarrow {\bf R}^N \}$ of duration $\tilde t$=$N\Delta t$ can be rewritten as the product $P_{\tilde t}=\prod_{n=1}^{N-1}p({\bf R}^n\rightarrow{\bf R}^{n+1})$ of single-step probabilities. Known analytical results [@Chandrasekhar1943] for $p$ yield the discretized OM action $$\label{eq.S}
{\mathcal S}/\beta=U({\bf R}^1)+\sum_{n=1}^{N-1}\sum_{\alpha=1}^D\frac{m_\alpha\nu}{4\Delta t}\left(R^{n+1}_\alpha-R_\alpha^n-\frac{\Delta t}{m_\alpha\nu}F^n_\alpha\right)^2,$$ where $\beta=1/k_BT$ and the potential energy $U({\bf R}^1)$ has been introduced to account for the Boltzmann weight of the initial configuration.

Using equation , equation becomes formally equivalent to the classical partition function of an open polymer. Hence, one can sample paths distributed according to equation adopting a Hamiltonian approach as done in path integral MD [@Parrinello1984]. Accordingly, we define the effective potential $\mathcal U = \mathcal S/\beta$, introduce auxiliary momenta and masses {${\rm P}_\alpha^n$, ${\rm M}_\alpha$} and solve Hamilton’s equations $$\begin{aligned}
\label{eq4.Hp}
\dot {\rm P}^{n}_\alpha&=&-\nabla_{{R}^{n}_\alpha} \mathcal{U}\\
\dot {R}^{n}_\alpha&=&\frac{{\rm P}^{n}_\alpha}{{\rm M}_\alpha} \label{eq4.Hr}\end{aligned}$$ coupled to a thermostat.

Although this approach was suggested in Ref. [@Dellago1998a], it was not further pursued mainly because of its computational burden. One issue is the presence in equation  of terms containing second derivatives of the potential. This problem is circumvented by using a symmetric finite difference formula [@Putrino2000; @Kapil2016] $$\label{eq5.FD}
\sum_{\alpha=1}^D\frac{\partial F_\alpha^n}{\partial R_k^n}\eta_\alpha^n\approx\frac{F_k^n(R_\alpha^n+\varepsilon \eta_\alpha^n)-F_k^n(R_\alpha^n-\varepsilon \eta_\alpha^n)}{2\varepsilon}$$ where $\eta_\alpha^n=R_\alpha^{n+1}-R_\alpha^n-\frac{\Delta t}{m_\alpha\nu}F_\alpha^n$ and $\varepsilon$ is a number small enough to guarantee energy conservation in microcanonical simulations. Equation  amounts to a modest but necessary increase in computational cost as it avoids direct implementation of the Hessian $\frac{\partial^2  U({\bf R})}{\partial R_\alpha \partial R_\beta}$. Adopting this method, one time step in path space involves 3$\times$$N$$\times$$D$ force evaluations. This has to be compared with the cost of $N$ MD steps in configurational space, which involve $N$$\times$$D$ force evaluations. However, while the standard approach is intrinsically serial, the path approach has the advantage that it can be made highly parallel. Specifically, here we adopt the hyper-parallel scheme of Calhhoun [*et al.*]{} [@Calhoun1996] and implement the algorithm in the LAMMPS [@Plimpton1995] suite of codes.

\[fig2.trap\] ![Probability distribution of the positions of all beads obtained from (a) dynamical TPS and (b) MetaD in path space. (c) The average value of the external potential, plotted as a function of time, $t=n\Delta t$, along RPs sampled using the two techniques. The MetaD result is obtained without reweighting. Areas indicate standard deviations computed from 5 independent runs.](fig_path-trap.pdf "fig:"){width="0.8\columnwidth"} ![Probability distribution of the positions of all beads obtained from (a) dynamical TPS and (b) MetaD in path space. (c) The average value of the external potential, plotted as a function of time, $t=n\Delta t$, along RPs sampled using the two techniques. The MetaD result is obtained without reweighting. Areas indicate standard deviations computed from 5 independent runs.](fig_2Davgpot.pdf "fig:"){width="0.7\columnwidth"}

As a first test case, we consider the diffusion of a particle in the two dimensional double-well potential of figure \[fig1.2Dpot\](a) [@Dellago1998a]. The latter provides a simple example of a system with multiple RPs connecting meta-stable states. Simulations in path space are carried out using a polymer of $N$=200 beads and auxiliary masses $M$=1. The time step, damping coefficient and mass at the Brownian level are set equal to $\Delta t$=0.15, $\nu$=1, $m$=1. Temperature is controlled via a Nosé-Hoover chains thermostat [@Martina1992] and the equations of motion ,  are propagated adopting a standard velocity-Verlet integrator with time step $\Delta t_{\rm MD}$=0.01. Results are reported as obtained from simulations with the above set of unitless parameters [@SM].

In the first set of simulations, we enforce harmonic constraints [@SM] to pin the initial ($n$=1) and final ($n$=N) replica, respectively close to the left and right minimum (A and B in figure \[fig1.2Dpot\](a)). This setup is equivalent to the dynamical TPS algorithm [@Dellago1998a] that focuses on sampling RPs connecting known initial and final states. We start from an equilibrated polymer configuration realizing a path crossing the lower saddle (see figure \[fig1.2Dpot\](b)) and run 5 unbiased simulations of duration 5$\times$10$^7$ MD steps at temperature $k_BT$=0.05, much smaller than the minimum energy barrier $\Delta E\sim1$ separating the two minima. With this setup, the central maximum of the potential prevents the polymer from visiting path configurations that connect A and B passing via the upper saddle. This is demonstrated in figure \[fig2.trap\](a), showing the histogram of the positions visited by all beads during simulation, averaged over the 5 independent runs. This is a typical example of “path trapping” in a local minimum of the OM action.

MetaD is designed precisely to overcome this type of problems. Moreover, as we will show, it can be applied directly to address the general situation where the final configuration is not known a priori. To this end, we first seek a suitable collective variable (CV) to which the bias is applied. For the present example the choice is straightforward and we select the polymer end-to-end distance $d_{\rm e2e}$=$\lvert {\bf R}^{N}-{\bf R}^1\rvert$ (see figure \[fig1.2Dpot\](b)). Using the enhanced sampling library PLUMED [@Tribello2014], we perform well-tempered MetaD [@Barducci2008; @Dama2014] while constraining the initial replica to reside in the left basin [@SM]. All other beads are free. Also in this case, we perform 5 independent runs of duration 5$\times$10$^7$ MD steps each. Figure \[fig2.trap\](b) reports the histogram of the positions visited by all beads during simulation. While most of the time the polymer remains crumpled within the left basin, the bias enables sampling of elongated configurations corresponding to successful RPs (see also insets of figure \[fig3.FES2D\]). To further illustrate the results, in figure \[fig2.trap\](c) we compare the average value of the potential felt by each replica during the unbiased simulations (red curve) with that computed using the RPs harvested using MetaD (black curve), showing good agreement. This demonstrates that the latter provide a good representation of the transition path ensemble sampled by TPS, even when obtained in presence of the bias.

We now turn to the computation of kinetic rates for transitions between two long-lived states A and B. These can be extracted from simulations computing the time correlation function [@Chandler1987] $$\label{eq6.fAB}
    f_{\rm AB}(t)=\frac{\langle I_{\rm A}(0)I_{\rm B}(t) \rangle}{\langle I_{\rm A}(0)\rangle}.$$ Here, we have introduced the characteristic function $I_{\rm X}(t)$, which takes a value of 1 if ${\bf R}(t)\in$X, and 0 otherwise. Averages, $\langle\cdot\rangle$, are intended over all possible unconstrained paths. Equation expresses the probability of finding the system in state B at time $t$, given that it was in state A at $t=0$. In absence of intermediate states and after a short transient, $f_{\rm AB}$ enters a linear regime and the phenomenological rate constant is given by its slope $$\label{eq.kAB}
 k_{\rm AB}=\frac{{\rm d} f_{\rm AB}(t)}{{\rm d} t}.$$

\[fig3.FES2D\] ![Free energy curves as a function of the end-to-end distance. Insets show representative polymer configurations at $d_{\rm e2e}\sim0$ and $d_{\rm e2e}\sim2$.](fig_FES-2Dpot.pdf "fig:"){width="0.8\columnwidth"}

Our method allows computing $f_{\rm AB}(t)$ from a single MD simulation in path space performed without imposing any constraint to the polymer. Alternatively, one notices that equation can be equivalently rewritten as $f_{\rm AB}(t)$=$\langle I_{\rm B}(t) \rangle_A$, where averages are now taken exclusively over paths initiating in A. The latter, in turn, can be sampled adopting a MetaD setup similar to the one described in previous section, adding suitable constraint on the first bead.

We first consider the two dimensional diffusion problem of figure \[fig1.2Dpot\]. Biased simulations are performed using the convergence variant of the recent OPES [@invernizzi2019] method, which turned out to perform better than MetaD. Instead of the bias, OPES focuses on reconstructing the probability distribution, leading to improved convergence speed. Moreover, the method offers a more straightforward reweighting scheme and an efficient importance sampling that avoids visiting uninteresting high free energy regions [@invernizzi2019]. We further exploit the possibility of targeting a specific probability distribution, which, for the problem at hand, we choose to be a uniform distribution of $d_{\rm e2e}$ between 0 and 3 [@SM]. Finally, we do not impose any constraint on the polymer and compute $f_{\rm AB}$ directly via equation . We adopt the same parameters reported above and run simulations in the range $0.125\leq k_BT\leq0.5$.

Figure \[fig3.FES2D\] shows the obtained free energy curves as a function of the end-to-end distance. At all temperatures considered, we observe a global minimum near $d_{\rm e2e}\sim0$, whose population corresponds to paths that never leave the left basin (see left inset). A secondary minimum is found at $d_{\rm e2e}\sim2$, which is the distance between the two minima of the potential. Accordingly, the corresponding polymer configurations represent RPs that successfully reach B after having crossed one of the two equivalent saddles (see right inset). Increasing the temperature lowers the free energy barrier, which reflects the thermally enhanced probability for the particle to cross from A to B. At low temperatures, configurations corresponding to intermediate values of $d_{\rm e2e}\sim1$ are “failed attempts” of the polymer that stretches towards the right basin without reaching it. On the other hand, for $k_BT\geq$0.4 we observe contributions also from configurations recrossing from B to A.

\[fig4.2Drate\] ![(a) The time correlation function, $f_{\rm AB}(t)$, computed at different temperatures. Dashed lines are fit to $k_{\rm AB}t+a$. Labels are scaling factors used for clarity of presentation. (b) Arrhenius plot of the rate constants. The red dashed line is a fit to $-\Delta E_{\rm fit}/k_BT+b$.](fig_rate-2Dpot.pdf "fig:"){width="0.8\columnwidth"}

In figure \[fig4.2Drate\](a) we show the correlation function $f_{\rm AB}(t)$ obtained at various temperatures, displaying the expected short transient, followed by linear growth. Note that for $k_BT\geq$0.4, recrossings of the polymer from B to A are observed at times $t\ge$10, where the conditions for equation to hold cease to apply [@Chandler1987]. Nevertheless, we were able to extract the phenomenological rate constant $k_{\rm AB}$ in the whole range of temperatures by considering the initial linear regime. These are reported in the Arrhenius plot of figure \[fig4.2Drate\](b), showing the expected linear trend of $\log(k_{\rm AB})$. A linear fit of the data yielded an activation energy of $\Delta E_{\rm fit}=1.05\pm0.02$, in agreement with the exact value of the potential barrier $\Delta E\simeq1.08$.

As a second example, we consider the umbrella inversion transformation of ammonia in vacuum. In this configurational change the nitrogen atom passes through the hydrogen plane to reach an equivalent and symmetric position. Thus, the process can be conveniently described in terms of the oriented height $h$ of the NH$_3$ tetrahedron. Here, we describe intra-molecular forces using an empirical model [@Weismiller2010]. Under normal conditions, $h$ takes values of $\sim\pm0.4$Å, the two equivalent configurations being separated by a large barrier of $\Delta E\simeq$120 kJ/mol$\sim$50 $k_BT_{\rm room}$ [@SM].

Simulations in path space are performed with $N$=100, $\Delta t$=0.1 fs and $\nu$=0.14 fs$^{-1}$. Auxiliary masses are set equal to those of the corresponding atom in the replica and the MD integration step is set to $\Delta t_{\rm MD}$=0.25 fs. Instructed by the previous study, we choose as CV the generalized end-to-end distance $\Delta h_{\rm e2e}=(h^{N}-h^{1})$. Biased simulations are performed using the convergence variant of OPES [@invernizzi2019] targeting a uniform distribution of $\Delta h_{\rm e2e}$ between -1Å and 1Å [@SM]. Also in this case, we do not impose any constraint on the polymer and compute $f_{\rm AB}$ directly via equation .

\[fig5.FES-NH3\] ![Free energy curves as a function of the generalized end-to-end distance $\lvert\Delta h_{\rm e2e}\rvert$. Insets show representative polymer configurations at $\lvert\Delta h_{\rm e2e}\rvert\sim0$ and $\lvert\Delta h_{\rm e2e}\rvert\sim0.8$Å. Replicas have been suitably aligned for clarity of presentation.](fig_FES-NH3.pdf "fig:"){width="0.8\columnwidth"}

Figure \[fig5.FES-NH3\] reports the free energy curves obtained at temperature $T$=300, 500, 700, 900 K. Overall, we observe similar features to those of the previous example, namely, a global minimum at $\lvert \Delta h_{\rm e2e}\rvert\sim0$ Å followed by a plateau and a shallow minimum at $\lvert \Delta h_{\rm e2e}\rvert\sim0.8$ Å. Insets show representative configurations corresponding respectively to a diffusive path where all replicas share the same orientation and a successful RP where ammonia flips between the $h$=$\pm0.4$ Å states. In figure \[fig6.rateNH3\](a) we report the time correlation function $f_{\rm AB}(t)$ extracted from simulations. For the chosen set of parameters, the polymer length adopted is sufficient to observe the onset of the linear regime, from which we extract the phenomenological rates. The latter are displayed in the Arrhenius plot of panel (b), along with results from unbiased MD simulations in configurational space [@SM]. In the latter case, rates are accessible only at high temperatures. At T=900 K, we obtain respectively $k_{\rm AB}^{\rm OPES}\simeq$2.4$\times$10$^{-4}$ and $k_{\rm AB}^{\rm direct}\simeq$4.0$\times$10$^{-4}$ ps$^{-1}$, in fair agreement. On the other hand, the activation barrier of $\Delta E_{\rm fit}$=$118\pm1$ kJ/mol estimated using the results from OPES perfectly matches the exact value of $\Delta E\simeq$120 kJ/mol of the adopted force field.

To conclude, we have presented a method to sample RPs via biased MD simulations in path space. The use of enhanced sampling techniques enables unconstrained exploration of RPs, making this approach more robust against problems like path trapping. Time correlation functions can be computed via straightforward (reweighted) averages along a single MD trajectory. Dynamical information such as kinetic rates are easily accessible at no additional cost. At the same time, adding suitable contraints at the polymer endpoints allows identification of most probable paths via direct minimization of the discretized action, which makes this approach very flexible. The method promises to be a powerful tool with potential applications in many fields including the study of chemical reactions, via implementation within the Car-Parrinello MD approach [@Car1985], and of the kinetics of biological systems.

\[fig6.rateNH3\] ![(a) The time correlation function, $f_{\rm AB}(t)$, computed at different temperatures. Dashed lines are fit to $k_{\rm AB}t+a$. Labels are scaling factors used for clarity of presentation. (b) Arrhenius plot of the rate constants. Circles are results from direct unbiased MD simulations in configurational space. The red dashed line is a fit to $-\Delta E_{\rm fit}/k_BT+b$ of the rate constants obtained from OPES.](fig_rate-NH3.pdf "fig:"){width="0.8\columnwidth"}

In the present work we have adopted MetaD and OPES as biasing schemes, but any other enhanced sampling schemes could be applied as well. As in all biased MD simulations, prior knowledge about the mechanisms underlying the transition of interest is needed in order to build successful CVs. This is crucial to speed up convergence. It is encouraging that in the cases studied here a suitably defined end-to-end distance performed well. This represents the most natural choice. Complex systems will definitely require more fine tuning. However, this should not pose any major problem as one is assisted in the task by the vast literature on the subject [@Valsson2016].

Finally, we note that the proposed path approach effectively realizes parallelization at the level of time [@Rosa-Raices2019]. This, in turn, allows a highly parallel implementation [@Calhoun1996] that takes full advantage of modern massively parallel computer architectures. Given the increasing availability of massive parallel computational resources, we believe that this method will find many successful applications.

Supplemental Materials
======================

Diffusion in a two dimensional double-well potential
----------------------------------------------------

In the first application, we considered the problem of a particle diffusing in the two dimensional potential [@Dellago1998a] $$U(x,y)= 2 + \frac 4 3 x^4 - 2y^2 + y^4 + \frac {10} {3} x^2(y^2-1).$$ Simulations were performed implementing the corresponding Onsager-Machlup (OM) action (equation (5) of the main text) using the following set of parameters at the Brownian level: $m$=1, $\Delta t$=0.15, $\nu$=1. In all simulations, we considered a polymer of size $N$=200 beads inside a square cell of side $L$=10, centered in the origin. Exact analytical expressions of all terms appearing in the definition of the forces governing the dynamics of the polymer (equation (6) of main text) have been hard coded in the LAMMPS [@Plimpton1995] suite of codes. The equations of motion in path space were then solved adopting a standard velocity-Verlet integrator with a time step of $\Delta t_{\rm MD}$=0.01 and auxiliary masses set to $M$=1. Temperature was controlled via a Nosé-Hoover chains thermostat [@Martina1992]. The PLUMED [@Tribello2014] enhanced sampling library was used to introduce harmonic restraints and to perform well-tempered metadynamics [@Barducci2008; @Dama2014] (WT-MetaD) as well as OPES [@invernizzi2019] simulations. Here and in the main text, all quantities are reported as obtained from simulations with the above set of unitless parameters.

In the first set of simulations, we investigated the ability of our metadynamics (MetaD) approach to sample different reactive pathways (RPs) in the same run, and we compared results with those obtained adopting the dynamical transition path sampling (TPS) algorithm of Ref. . TPS simulations were performed at temperature $k_BT$=0.05, applying harmonic constraints (spring constant $K$=100) to the two distances $R_{1,{\rm A}}$=$\lvert {\bf R}^1-{\bf R}_{\rm A} \rvert$ and $R_{N,{\rm B}}$=$\lvert {\bf R}^N-{\bf R}_{\rm B} \rvert$. Here, ${\bf R}^{1,N}$ are the positions of the first and last bead of the polymer, while ${\bf R}_{\rm A,B}\simeq(\mp1,0)$ mark the two minima of the potential. With this choice of parameters, $R_{1,{\rm A}}$ and $R_{N,{\rm B}}$ were bound to values $\leq$0.1. WT-MetaD simulations were performed at the same temperature, applying same harmonic constraint only to $R_{1,{\rm A}}$. The polymer end-to-end distance $d_{\rm e2e}$=$\lvert {\bf R}^1-{\bf R}^N\rvert$ was used as collective variable (CV). The bias potential was built using a bias factor $\gamma$=20, depositing Gaussian kernels (height=1, $\sigma$=0.1) every 2500 MD steps. In both cases, we performed simulations of 5$\times$10$^{7}$ MD steps, sampling configurations every 2500 steps. In the case of TPS, all polymer configurations sampled were used to compute the average potential profile (figure 2(c) of main text). In the case of MetaD, we considered only the second half of the trajectory (well within the asymptotic regime of WT-MetaD) and defined successful RPs those for which $R_{1,{\rm A}}$$\leq$0.1 and $R_{N,{\rm B}}$$\leq$0.1 at the same time. These were used to compute the average potential profile. Each experiment was repeated 5 times. The curves shown in figure 2(c) of the main text report the associated average values and standard deviations.

\[fig.trace2D\] ![(a) Time evolution of $d_{\rm e2e}$ obtained from biased simulation at $k_BT$=0.125 using OPES [@invernizzi2019]. (b) The corresponding biased probability distribution.](fig_trace-and-prob-2Dpot.pdf "fig:"){width="0.9\columnwidth"}

\[fig.converge2D\] ![Time evolution of the correlation function $f_{\rm AB}(t=30)$, obtained from reweighting of 4 independent simulations at $k_BT$=0.125.](fig_fAB-converge-2D.pdf "fig:"){width="0.9\columnwidth"}

In the second set of simulations, we adopted the convergence variant of OPES [@invernizzi2019] to compute the phenomenological kinetic rate, $k_{\rm AB}$, for transitions between A and B. Simulations were performed at temperature $k_BT$=0.125, 0.2, 0.3, 0.4 and 0.5, without applying any constraint to the polymer. We used a bias factor $\gamma$=0 and deposited Gaussian kernels every 2500 MD steps. The standard deviation of the kernels was set equal to $\sigma$=0.1 at $k_BT$=0.125 and to $\sigma$=0.2 at all other temperatures. In OPES, the height of the kernels is automatically adjusted during runtime. The last input parameter is an estimate of the free energy barrier to be overcome, which we set equal to $\Delta F$=1.5 for $k_BT$=0.125, 0.2, 0.3 and to $\Delta F$=1.2 for $k_BT$=0.4, 0.5. These values were estimated from preliminary WT-MetaD simulations. With this setup, we effectively targeted a flat probability distribution of $d_{\rm e2e}$ in the interval 0$\lesssim d_{\rm e2e}\lesssim$3 (see figure \[fig.trace2D\]). At each temperature, we performed a simulation of 5$\times$10$^{8}$ MD steps. The time correlation function $f_{\rm AB}$ (equation (9) of the main text) was computed defining the characteristic functions of the two basins of the potential as $I_{\rm A,B}({\bf R}^n)$=1 if $\lvert {\bf R}^n-{\bf R}_{\rm A,B}\rvert<$0.7 and zero otherwise. Time averages were computed using the reweighting scheme suggested in Ref. [@invernizzi2019], neglecting the initial 5$\times$10$^{4}$ steps (see figure \[fig.converge2D\]). Each experiment was repeated 4 times. Figures 3 and 4(a) of the main text report the associated average values and standard deviations.

\[fig.pes\] ![Potential energy of NH$_3$ as a function of its oriented heigth, $h$.](fig_PES-NH3.pdf "fig:"){width="0.9\columnwidth"}

Umbrella inversion of NH$_3$ in vacuum
--------------------------------------

In the second application, we considered the umbrella inversion of ammonia in vacuum. Intra-molecular interactions were described using the ReaxFF force field of Ref. , neglecting electrostatics. We used a cubic box of side $L$=30 Å, larger than twice the cutoff of 14 Å adopted for long range van der Waals interactions, and we enforced periodic boundary conditions. Figure \[fig.pes\] shows the potential energy profile obtained via a sequence of geometry optimizations at fixed values of the oriented height, $h$, of the NH$_3$ tetrahedron, which we used to describe the transformation. The two symmetric minimum energy confingurations at $h$$\sim\pm$0.4 Å are separated by a barrier of $\sim$120 kJ/mol.

\[fig2.trace-e2e\] ![(a) Time evolution of $\Delta h_{\rm e2e}$ obtained from biased simulation at $T$=500 K using OPES [@invernizzi2019]. (b) The corresponding biased probability distribution.](fig_trace-and-prob-NH3.pdf "fig:"){width="0.9\columnwidth"}

\[fig.convergeNH3\] ![Time evolution of the correlation function $f_{\rm AB}(t=8\,{\rm ps})$, obtained from reweighting of 4 independent simulations at $T$=500 K.](fig_fAB-converge-NH3.pdf "fig:"){width="0.9\columnwidth"}

Simulations in path space were performed setting $\Delta t$=0.1 fs, $\nu$=0.14 fs$^{-1}$ in the OM action and using a polymer of length $N$=100 beads. The equations of motion have been implemented in LAMMPS [@Plimpton1995] and solved adopting a velocity-Verlet algorithm with an MD integration step of $\Delta t_{\rm MD}$=0.25 fs. The auxiliary masses $M_{\alpha}$ were set equal to those of the corresponding atom in the ammonia replica. Terms in the forces containing second derivatives of the potential energy were estimated adopting the finite difference expression of equation (8) of the main text. Biased MD simulations were performed adopting the convergence variant of OPES [@invernizzi2019], as implemented in PLUMED [@Tribello2014]. We considered temperatures of $T$=300, 500, 700, 900 K and we did not apply any constraint to the polymer. We chose as CV the generalized polymer end-to-end distance $\Delta h_{\rm e2e}=(h^N-h^1)$, equal to the difference between the oriented height of the last and first ammonia replica. We used a bias factor $\gamma$=0 and deposited Gaussian kernels every 500 MD steps. The standard deviation of the kernels was set equal to $\sigma$=0.0032, 0.0043, 0.0051, 0.0059 Å, respectively for $T$=300, 500, 700, 900 K. In OPES, the height of the kernels is automatically adjusted during runtime. The last input parameter is an estimate of the free energy barrier to be overcome. This was estimated from preliminary WT-MetaD simulations and set equal to $\Delta F$=140 kJ/mol. With this setup, we effectively targeted a flat probability distribution of $\Delta h_{\rm e2e}$ in the interval -1 Å$\lesssim \Delta h_{\rm e2e}\lesssim$1 Å (see figure \[fig2.trace-e2e\]). At each temperature, we ran a simulation of $\sim$3$\times$10$^8$ MD steps. The time correlation function $f_{\rm AB}(t)$ was computed defining the characteristic function of the two basins $I_{\rm A,B}(h)$=1 if $\lvert h\rvert>0.02$ Å and zero otherwise. Time averages were computed using the reweighting scheme suggested in Ref. [@invernizzi2019], skipping the initial $\sim$10$^5$ MD steps (see figure \[fig.convergeNH3\]). Each experiment was repeated 4 times. Figures 5 and 6(a) of the main text report the associated average values and standard deviations.

We further performed simulations of a single NH$_3$ molecule at high temperatures $T=$900-1700 K, solving directly equation (4) of the main text ($\Delta t$=0.1 fs, $\nu$=0.14 fs$^{-1}$). At each temperature, we ran a simulation of 2$\times$$10^9$ MD steps. Following Ref. [@Salvalaglio2014], the phenomenological rate constants were obtained from a fit of the cumulative distribution function of first passage times (see figure \[fig.CDF\]). Each experiment was repeated 5 times. Black circles in figure 6(b) of the main text report the associated average values and standard deviations.

\[fig.CDF\] ![The cumulative distribution function of first passage times, extracted from direct simulations at $T$=1100 K. The red dashed line is a fit to $1-e^{-k_{\rm AB}t}$, used to extract the phenomenological rate constant.](fig_CDF.pdf "fig:"){width="0.9\columnwidth"}

[36]{} A. Laio and M. Parrinello, PNAS [**99**]{}, 12562 (2002). P. Tiwary and M. Parrinello, Phys. Rev. Lett. [**111**]{}, 230602 (2013). H. Wu, A. S. J. S. Mey, E. Rosta, and F. Noè, J. Chem. Phys. [**141**]{}, 214106 (2014). E. Rosta and G. Hummer, J. Chem. Theory Comput. [**11**]{}, 276 (2015). H. Wu, F. Paul, C. Wehmeyer, and F. Noè, PNAS [**113**]{}, E3221 (2016). L. Donati, C. Hartmann, and B. G. Keller, J. Chem. Phys. [**146**]{}, 244112 (2017). L. Donati and B. G. Keller, J. Chem. Phys. [**149**]{}, 072335 (2018). R. Elber and M. Karplus, Chem. Phys. Lett. [**139**]{}, 375 (1987). R. Olender and R. Elber, J. Chem. Phys. [**105**]{}, 9299 (1996). D. Passerone and M. Parrinello, Phys. Rev. Lett. [**87**]{}, 108302 (2001). J. Lee, I.-H. Lee, I. Joung, J. Lee, and B. R. Brooks, Nat. Commun. [**8**]{}, 15443 (2017). L. R. Pratt, J. Chem. Phys. [**85**]{}, 5045 (1986). C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, J. Chem. Phys. [**108**]{}, 1964 (1998). H. Fujisaki, M. Shiga, and A. Kidera, J. Chem. Phys. [**132**]{}, 134101 (2010). T. S. van Erp, D. Moroni, and P. G. Bolhuis, J. Chem. Phys. [**118**]{}, 7762 (2003). D. Moroni, P. G. Bolhuis, and T. S. van Erp, J. Chem. Phys. [**120**]{}, 4055 (2004). T. J. H. Vlugt and B. Smit, Phys. Chem. Comm. [**4**]{}, 11(2001). E. E. Borrero and C. Dellago, Eur. Phys. J. Spec. Top. [**225**]{}, 1609 (2016). L. Onsager and S. Machlup, Phys. Rev. [**91**]{}, 1505 (1953). S. Chandrasekhar, Rev. Mod. Phys. [**15**]{}, 1 (1943). M. Parrinello and A. Rahman, J. Chem. Phys. [**80**]{}, 860 (1984). A. Putrino, D. Sebastiani, and M. Parrinello, J. Chem. Phys. [**113**]{}, 7102 (2000). V. Kapil, J. Behler, and M. Ceriotti, J. Chem. Phys. [**145**]{}, 234103 (2016). A. Calhoun, M. Pavese, and G. A. Voth, Chem. Phys. Lett. [**262**]{}, 415 (1996). S. Plimpton, J. Comput. Phys. [**117**]{}, 1 (1995). G. J. Martyna, M. L. Klein, and M. Tuckerman, J. Chem. Phys. [**97**]{}, 2635 (1992). See Supplemental Material at URL for a detailed discussion of all simulation setup, which includes Ref. [@Salvalaglio2014]. M. Salvalaglio, P. Tiwary, and M. Parrinello, J. Chem. Theory Comput. [**10**]{}, 1420 (2014). G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, and G. Bussi, Comput. Phys. Commun. [**185**]{}, 604 (2014). A. Barducci, G. Bussi, and M. Parrinello, Phys. Rev. Lett. [**100**]{}, 020603 (2008). J. F. Dama, M. Parrinello, and G. A. Voth, Phys. Rev. Lett. [**112**]{}, 240602 (2014). D. Chandler, [*Introduction to Modern Statistical Mechanics*]{}, edited by O. U. Press (New York, 1987). M. Invernizzi and M. Parrinello, arXiv:1909.07250 (2019). M. R. Weismiller, A. C. T. van Duin, J. Lee, and R. A. Yetter, J. Phys. Chem. A [**114**]{}, 5485 (2010). R. Car and M. Parrinello, Phys. Rev. Lett. [**55**]{}, 2471 (1985). O. Valsson, P. Tiwary, and M. Parrinello, Annu. Rev. Phys. Chem. [**67**]{}, 159 (2016). J. L. Rosa-Raı́ces, B. Zhang, and T. F. Miller, J. Chem. Phys. [**151**]{}, 164120 (2019).
