The Kramers turnover in terms of a macro-state projection on phase space

In order to reproduce the classic Kramers turnover, we selected 25 friction coefficient values $\gamma$ between $0.1\,\mathrm{ps^{-1}}$ and $30.0\,\mathrm{ps^{-1}}$ and we solved the Langevin eq. 1 using the Brünger, Brooks and Karplus (BBK) integrator scheme [19] with an integrator timestep $\Delta t=0.005\,\mathrm{ps}$. For each value of $\gamma$, we ran 500 simulations starting at the bottom of the left well of the potential with an initial momentum randomly drawn from the Boltzmann distribution. After the particle reached the bottom of the right well, the simulations were stopped and we calculated the mean time, i.e. the Mean First Passage Time (MFPT) $\langle\tau_{fp}\rangle$ [20], from which we obtained the transition rate as

The results of this numerical experiment are reported in fig. 2-(a) as black squares. If the friction is very low ($\gamma\approx 0.1\,\mathrm{ps}^{-1}$), the dynamics of the system (eq. 1) is almost deterministic, and the system, unless it has enough initial momentum, is trapped in the well with an extremely low probability of escape. Correspondingly, the MFPT is very large and the value of the rate tends to zero. However, increasing the friction by a small amount ($\gamma\approx 1.5\,\mathrm{ps}^{-1}$) the system gets enough thermal energy through the random force and increases the probability to escape from the well. In fact, for low values of $\gamma$, the stochastic force $\sqrt{2k_{B}T\gamma m}\,\eta_{t}$, which is weighted by the square root of the friction coefficient, dominates the friction force $-\gamma p_{t}$, which is linear with the friction coefficient. Thus, we observe a rapid and linear increase in rates up to a maximum of $k_{A\rightarrow C}\approx 0.02\,\mathrm{ps}^{-1}$. Beyond the threshold of $\gamma\approx 1.5\,\mathrm{ps}^{-1}$, the system enters what is called the moderate friction regime. Here, the friction force dominates the Langevin equation, and the probability of escaping the well, despite the high thermal energy, decays as $k_{A\rightarrow C}\propto\sqrt{1+\nicefrac{{1}}{{\gamma^{2}}}}$. For higher values of the friction coefficient ($\gamma>20\,\mathrm{ps}^{-1}$), the friction term is so strong that the average acceleration of the system tends to zero. The dynamics is overdamped and the transition rate decays as $k_{A\rightarrow C}\propto\nicefrac{{1}}{{\gamma}}$.

The three friction regimes, here qualitatively described, were formalized by Kramers in 1940 [2]. He assumed a two-metastable system governed by the Langevin dynamics with thermal energy $k_{B}T\ll E_{b}^{+}=V(q_{B})-V(q_{A})$, so as to ensure metastability. In addition, he required that the left well and the top of the barrier of the potential $V(q)$ are approximated by harmonic potentials with angular frequencies

Under these conditions, Kramers derived a transition rate formula for the low friction regime ($\gamma<\omega_{B}$)

for the moderate friction regime ($\gamma>\omega_{B}$)

and the high friction regime ($\gamma\gg\omega_{B}$)

Note that Kramers defines the three regimes by comparing the coefficient of friction with the angular frequencies of the harmonic potentials that approximate the potential. In fact, the transition probability also depends on the curvature near the pit and barrier. The prediction of Kramers’ formulas, reported in fig. 2-(a), is excellent, it is only around the threshold separating the low and moderate friction regimes that the model becomes inaccurate. For more details about Kramers theory, we recommend Refs. [21, 22].

In the second numerical experiment, we estimated the transition rates applying ISOKANN to the same setting of the first experiment. We generated 1000 random initial points $x_{0}=(q_{0},p_{0})$ from a uniform distribution over the phase space defined by the $q$-range $[-2.0,\,2.0]$ nm and the $p$-range $[-10.0,\,10.0]$ amu nm ps^-1, and for each state we simulated $N=100$ trajectories of length $\tau=7$ ps, corresponding to 1400 timesteps using a timestep integrator $\Delta t=0.005$ ps. The ISOKANN algorithm has been applied for 20 iterations using multiquadratic radial basis functions (RBF) [23], which are computationally undemanding and only require a few parameters to optimize during training, resulting in faster convergence. We considered two cases:

• •

  the membership functions $\chi_{A}(q)$ and $\chi_{C}(q)$, and the transition rate $k_{\chi_{A}\rightarrow\chi_{C}}$ between the macro-states of the position space;

• •

  the membership functions $\hat{\chi}_{A}(q,p)$ and $\hat{\chi}_{C}(q,p)$ defined on the two-dimensional phase space, and the transition rate $\hat{k}_{\chi_{A}\rightarrow\chi_{C}}$ between macro-states of the phase space. Note that from now on, each quantity marked with the symbol $\hat{\cdot}$ denotes a quantity measured over the entire phase space.

The two rates, as functions of the friction coefficient $\gamma$, are reported in fig. 2-(b), respectively as blue upside down triangles and red circles. Both curves show a turnover similar to the rate $k_{A\rightarrow C}$ reported in fig. 2-(a): rates have an ascending profile for very low range values, then, having reached the maximum ($k_{\chi_{A}\rightarrow\chi_{C}}\approx 0.4\,\mathrm{ps}^{-1}$ and $\hat{k}_{\chi_{A}\rightarrow\chi_{C}}\approx 0.2\,\mathrm{ps}^{-1}$), descend slowly. However, while the values of the rate $\hat{k}_{\chi_{A}\rightarrow\chi_{C}}$ in phase space are overlapping with those of the Kramers rate $k_{A\rightarrow C}$ (although they are different physical quantities), the rate $k_{\chi_{A}\rightarrow\chi_{C}}$ defined in position space turns out to be higher in the low friction region but converges to the values of $\hat{k}_{\chi_{A}\rightarrow\chi_{C}}$ in the high friction regime. To understand these results, it is useful to take a look at the membership functions obtained from ISOKANN and shown in fig. 3, where figures (d,e,f) on the second row and (g,h,i) on the third row are respectively the membership functions in the phase space and the position space, for low, moderate and high friction.

In fig. 3-(a) (low friction regime), the membership functions of the macro-states only have significant values for those states whose total energy $E=\nicefrac{{p^{2}}}{{2m}}+V(q)$ is less than the height of the barrier. The points with a total energy exceeding the barrier are depicted in white, indicating that they have an equal probability of belonging to either $\chi_{A}$ or $\chi_{C}$, approximately 0.5. This occurs because trajectories originating from this area undergo periodic oscillations in phase space, visiting both wells as depicted in fig. 1-(b) by the blue trajectory. Correspondingly, in fig. 3-(d), we show the membership function values as a projection of $\hat{\chi}_{A}(q,p)$ and $\hat{\chi}_{C}(q,p)$ onto the position space. The apparent noise is due to the fact that the membership functions on phase space are not constant along the axis of momenta. Therefore, when friction is low, the membership values projected to position space are not functions in a strict sense and do not describe position-based macro-states. In fig. 3-(b) (moderate friction regime), the membership functions draw concentric spirals that terminate in the minima of phase space respectively. This may seem counterintuitive, but observing how a trajectory behaves in the moderate friction regime helps to interpret the membership functions correctly. In fig. 1-(b), the red trajectory starts at position $q=-1\,\mathrm{nm}$ and momentum $p=-8\,\mathrm{amu\,nm\,ps^{-1}}$, and reaches the right-hand minimum by following a clockwise trajectory. Similarly, if we start trajectories that are far from the central region of the phase space, we would observe spiral patterns that match with the membership functions. From this figure, we deduce that in the moderate friction regime, the effective barrier, i.e. the transition region, is not at $q=0$, but between the two spirals. In particular, in the central box $[-1,\,1]\times[-5,\,5]$ of the phase space, the barrier corresponds to a diagonal line which is not parallel to the momentum axis. The reason is that if $q=0$ and $p>0$, the system reaches the right well with low probability of recrossing. Conversely, if $q=0$ and $p<0$, the system reaches the left region. Along the white diagonal the system is in an unstable equilibrium, i.e. the system has the same probability of reaching one of the wells, and ISOKANN assigns equal probability of membership to the two macro-states. In fig. 3-(c) (high friction regime), membership functions are independent from momentum space. The two regions $q<0$ and $q>0$ are assigned to the macro-states regardless of the momentum and the transition region is almost a vertical line. The projections $\chi_{A}(q)$ and $\chi_{C}(q)$ onto the position space also appear well defined in fig. 3-(f). This occurs because as the friction is very high, the momentum is quickly damped and it does not provide enough energy to overcome the barrier as shown by the green trajectory in fig. 1-(b). In the high friction regime, only the thermal noise can provide enough energy to jump over the barrier.

In order to validate our results, we constructed a reference solution by means of a grid-based technique similar to Ulam’s method [24] which allows to discretize the operator $\mathcal{K}_{\tau}$ defined in eq. 11 into a transition probability matrix $\mathbf{K}(\tau)$, cf. [25]. Given a discretization of the phase space $\Gamma$ into $M$ disjoint subsets $\Gamma_{i}$, with $i=1,\dots,M$, and a set of $N$ simulations of length $\tau$ started in a random position of the subset $\Gamma_{i}$, then the entries of the matrix $\mathbf{K}(\tau)$ are written as

where $\mathbf{1}_{\Gamma_{j}}$ is the indicator function

and $x_{i,n}^{\tau}$ is the final state of the $n$th simulation started in $\Gamma_{i}$. In practice, one counts how many times a simulation starting in $\Gamma_{i}$ ends in $\Gamma_{j}$ and divides by the number of simulations to obtain an estimation of the transition probability. Afterward, an approximation of the infinitesimal generator, sometimes referred to as pseudogenerator, is obtained as

where $\mathbf{I}$ denotes an identity matrix of the same size as $\mathbf{K}{\tau}$. Then the membership functions $\hat{\chi}(q,p)$ are calculated applying PCCA+ to $\widetilde{\mathbf{Q}}$ and the coarse-grained rate matrix between macro-states is calculated as a Galerkin projection of $\widetilde{\mathbf{Q}}$ onto the membership functions:

In eq. 33, $\mathrm{diag}(\pi)$ denotes an $M\times M$ diagonal matrix, whose diagonal entries are the entries of the Boltzmann distribution $\pi(q,p)$ (eq. 7) evaluated at the centers of subsets $\Gamma_{i}$. Assuming a two-metastable system, the rate matrix $\widetilde{\mathbf{Q}}$ has size $2\times 2$:

with $\widetilde{q}_{\chi_{A}\rightarrow\chi_{C}},\,\widetilde{q}_{\chi_{C}\rightarrow\chi_{A}}>0$ representing the transition rates between the macro-states. For the sole case of a bimetastable system, these rates are equivalent to the exit rates defined in eq. 13.

Here, we discretized the $q$-range $[-2.0,\,2.0]$ nm in $80$ intervals of the same length $\Delta q=0.05$ nm, and the $p$-range $[-10.0,\,10.0]$ nm in $70$ intervals of the same length $\Delta p=0.29$ nm. The transition rates estimated by PCCA+ are reported in fig. 2-(b) as black squares, while the membership functions are reported in fig. 3-(a,b,c). For each subset, we ran 500 simulations of length 7 ps, with an integrator timestep of 0.005 ps for a total of 1400 timesteps. There is excellent agreement between ISOKANN and the method based on the discretization of the phase space: both methods recreate the Kramers turnover and show the same patterns for the membership functions.