Interplay of fast and slow dynamics in rare transition pathways: the disk-to-slab transition in the 2d Ising model

TL;DR

This paper demonstrates that understanding rare transitions in the 2D Ising model requires considering both free energy landscapes and the system's effective dynamics, highlighting the importance of dynamical information for accurate mechanistic insights.

Contribution

It shows how combining free energy landscapes with state-dependent diffusion coefficients via a Smoluchowski equation improves understanding of rare transition mechanisms.

Findings

Static free energy alone can be misleading for rare events.

Dynamical information refines the understanding of transition pathways.

Effective dynamics are essential for accurate rate calculations.

Abstract

Rare transitions between long-lived stable states are often analyzed in terms of free energy landscapes computed as functions of a few collective variables. Here, using transitions between geometric phases as example, we demonstrate that the effective dynamics of a system along these variables are an essential ingredient in the description of rare events and that the static perspective provided by the free energy alone may be misleading. In particular, we investigate the disk-to-slab transition in the two-dimensional Ising model starting with a calculation of a two-dimensional free energy landscape and the distribution of committor probabilities. While at first sight it appears that the committor is incompatible with the free energy, they can be reconciled with each other using a two-dimensional Smoluchowski equation that combines the free energy landscape with state dependent diffusion coefficients. These results illustrate that dynamical information is not only required to calculate rate constants but that neglecting dynamics may also lead to an inaccurate understanding of the mechanism of a given process.

Equations

$$\mathcal{H}=\sum _{<ij,i^{\prime }{j}^{\prime }>}{\sigma }_{ij}{\sigma }_{i^{\prime }{j}^{\prime }},$$

$$d=\{\begin{array}{ll}{d}_{\text{image}}& \text{for non-spanning clusters}\\ -d_{\text{width}}+3& \text{for spanning clusters}\end{array},$$

$${\tilde{\sigma }}_{jk}=\frac{1}{n^2}\sum _{l,m=1}^n{\sigma }_{lm}\exp \left[\frac{2\pi i}{n}\left(jl+km\right)\right]$$

$$s=\min \left(\frac{\left|{\tilde{\sigma }}_{01}\right|}{\left|{\tilde{\sigma }}_{10}\right|},\frac{\left|{\tilde{\sigma }}_{10}\right|}{\left|{\tilde{\sigma }}_{01}\right|}\right).$$

$${⟨\delta {\tilde{d}}^2⟩}_t={⟨{\left(\tilde{d}(t)-{⟨\tilde{d}⟩}_t\right)}^2⟩}_t.$$

$$C_{ds}(t)=\frac{{⟨(\delta \tilde{d})(\delta \tilde{s})⟩}_t}{\sqrt{{⟨(\delta \tilde{d}{)}^2⟩}_t}\sqrt{{⟨(\delta \tilde{s}{)}^2⟩}_t}},$$

$$\left(D_{ij}\right)\left(\vec{x}\right)=\left(\begin{array}{cc}{D}_d\left(\vec{x}\right)& 0\\ 0& D_s\left(\vec{x}\right)\end{array}\right).$$

$$L=\prod _{\alpha }{\left(e^{t_{\alpha }R}\right)}_{i_{\alpha },j_{\alpha }}.$$

$$R_{ij}=\{\begin{array}{ll}-{\sum }_{l\ne i}{R}_{il}& \text{if }i=j,\\ R_{ji}{P}_i/P_j& \text{if }i<j,\end{array}$$

$$D_{i,i+1}=\Delta X^2{R}_{i+1,i}{\left(\frac{P_i}{P_{i+1}}\right)}^{1/2}.$$

$$\log \tilde{L}$$

$$+\sum _i\frac{{\left(D_{i+1,i}-D_{i,i-1}\right)}^2}{2{\epsilon }^2},$$

$${\tau }_{\text{cross}}\sim \frac{(\Delta d(s){)}^2}{D_d(s)}{e}^{\beta F^{†}(s)}$$

$$1=\beta \Delta F^{†}(s)\approx \beta \frac{\partial F^{†}(s)}{\partial s}\Delta s.$$

$${\tau }_{\text{drift}}\sim \left|\frac{\Delta s}{\dot{s}}\right|\propto \left|\frac{1/\left(\beta \frac{\partial F^{†}(s)}{\partial s}\right)}{D_s\beta \frac{\partial F(s)}{\partial s}}\right|$$

$$\frac{{\tau }_{\text{cross}}}{{\tau }_{\text{drift}}}=\frac{D_s}{D_d}{\beta }^2(\Delta d(s){)}^2\left|\frac{\partial F(s)}{\partial s}\frac{\partial F^{†}(s)}{\partial s}\right|e^{\beta F^{†}(s)},$$

$$\frac{\partial \rho }{\partial t}$$

$$+\sum _{i,j}\frac{{\partial }^2}{\partial x_i\partial x_j}\left[D_{ij}\rho \right],$$

$$\sum _{i,j}\left(\frac{\partial D_{ij}}{\partial x_j}-\beta D_{ij}\frac{\partial F}{\partial x_j}\right)\frac{\partial p_{\mathcal{B}}(\vec{x})}{\partial x_i}$$

$$+\sum _{i,j}{D}_{ij}\frac{{\partial }^2{p}_{\mathcal{B}}(\vec{x})}{\partial x_i\partial x_j}$$

$$\sum _i{n}_i{D}_{ij}\frac{\partial p_{\mathcal{B}}\left(\vec{x}\right)}{\partial x_j}=0,$$

$${p_{\mathcal{B}}(\vec{x})|}_{\partial \mathcal{A}}=0,{p_{\mathcal{B}}(\vec{x})|}_{\partial \mathcal{B}}=1$$

$$\frac{\text{d}{x}_i(t)}{\text{d}t}=\sum _j\left[\frac{\partial D_{ij}}{\partial x_j}-\beta D_{ij}\frac{\partial F}{\partial x_j}\right]+\sum _j{g}_{ij}{\eta }_j(t),$$

Full text

I Introduction

Free energy landscapes are a ubiquitous tool in the study of rare events—such as nucleation events in phase transitions or conformational changes in biomolecules—as they are often used to get a first sense of which regions in configuration space are relevant for a given process. Calculating such a landscape requires the choice of a small number of order parameters that are supposed to capture the important degrees of freedom of the process in question. Examples of such order parameters include the potential energy, the size of the largest cluster in the system (in the case of nucleation), or quantities like the number of native contacts, radii of gyration or root-mean-square deviations from the target structure (in biological systems).

However, care has to be taken in the interpretation of such free energy landscapes, since the choice of coordinates is not unique, even if the relevant degrees of freedom are captured. In particular, any non-linear transformation of the coordinates will change the free energy landscape, including features like the existence and height of free energy barriersFrenkel (2013). Due to this arbitrariness, additional information on the dynamics along a chosen order parameter is required in order to reliably interpret a free energy landscapeChahine et al. (2007); Krivov and Karplus (2008); Hinczewski et al. (2010); Best and Hummer (2010, 2011).

With higher-dimensional free energy landscapes additional complications arise, especially if there is no clear physical relationship between the different order parameters used, as is often the case. For instance, one may be tempted to associate the gradient in such a landscape with a probability flux. However, flow directionsTrinkaus (1983); Berezhkovskii and Szabo (2005); Peters (2009), the channels in which the transition proceedsNorthrup and McCammon (1983); Moro and Cardin (1998); Yang et al. (2007), and isocommittor linesWeinan and Vanden-Eijnden (2006); E et al. (2005) are determined both, by the free energy landscape and the dynamics of the system along the different order parameters and may deviate significantly from naive expectations. In particular, the combination of a fast and a slow degree of freedom which we will encounter in this work, has been discussed using a toy model by Metzner et al. (2006).

Often, especially in the context of rate calculations, one would like to project the system onto a single variable, called a reaction coordinate, that is not only capable of distinguishing between reactant and product state, but also describes the progress of the transformation in between. Due to the effects mentioned above, the free energy landscape may mislead us into a suboptimal choice for this reaction coordinate, which in turn negatively impacts the efficiency and accuracy of common rate calculation methods.

In the following, we use a simple model process to demonstrate some of the issues mentioned above: the disk-to-slab phase transition in the two-dimensional Ising model. Disk and slab phase are so called geometric phasesLeung and Zia (1990), which arise due to periodic boundary conditions used in simulations. They are characterized by different cluster shapes—spherical, cylindrical and slablike in three dimensions—where the stability of the non-spherical shapes is due to the reduction of surface free energy that can be achieved by connecting a cluster to its periodic images. In two dimensions the number of geometries reduces to two: a disk- and a slab phase.

The disk-to-slab transition is an example of a transition between different cluster (or bubble) geometries that play a role in different physical situations Tröster et al. (2005); Singh et al. (2009, 2011); Santos et al. (2010); Kumar et al. (2011); Prestipino et al. (2015), including the dewetting transition observed in volumes that are confined between hydrophobic surfaces Lum and Chandler (1998); Nicolaides and Evans (1989); Evans (1999); Bolhuis and Chandler (2000); Leung and Luzar (2000); Luzar and Leung (2000); Vishnyakov and Neimark (2003). The mechanism of this dewetting has recently been discussed in detail by Remsing et al. (2015), who showed that it may involve the formation of a vapor bubble on one of the surfaces which subsequently changes its shape to a vapor tube that connects the plates. Under some conditions, this mechanism reduces the free energy barrier that has to be overcome, by, in essence, circumventing the states that involve a vapor tube size that is close to critical. The disk-to-slab phase transition can be seen as a simplified version of this process.

In addition, geometric phases pose an obstacle in the sampling of free energy landscapes Neuhaus and Hager (2003); Trebst et al. (2004) and they can be used to calculate surface tensions Tröster and Binder (2011); Tröster et al. (2012).

The remainder of this paper is structured as follows: section II specifies the model that is investigated and introduces two order parameters that are then used to investigate the transition. In section III we discuss the free energy landscape and the distribution of committor probabilities as well as the apparent disconnect between the two. Section IV explains how diffusion coefficients can be calculated reliably in this system, which are then used, in section V, to build a two-dimensional model of the system that predicts the distribution of committor probabilities. In section VI we recapitulate our findings and discuss some of their implications.

II Model and Order Parameters

We investigate the disk-to-slab phase transition in the two-dimensional Ising model with ferromagnetic nearest-neighbor interactions on a square lattice of size $N=100\times 100$ and a vanishing magnetic field. The Hamiltonian is given by

[TABLE]

where the sum extends over all nearest-neighbor pairs and we have set the coupling constant to 1. The spins, denoted by $\sigma_{ij}$ with $i$ being the position in $x$-direction and $j$ the position in $y$-direction, can take on values of $+1$ (up) and $-1$ (down) and the magnetization per spin of the system is given by $m=(\sum_{i,j}\sigma_{ij})/N$. Periodic boundary conditions apply in all directions. Since under these conditions the model is symmetric with respect to a flip of all spins, we restrict the discussion to positive $m$.

In this case, given a fixed temperature $T$ below the critical temperature $T_{\text{c}}$, three distinct phases can be found in order of descending magnetization: a homogeneous phase without a single dominant cluster of down spins, a disk phase where the bulk of the spins pointing down are in a single disk-shaped (non-spanning) cluster and a slab phase where this cluster is in contact with its own periodic copies (spanning). See Fig. 1 for examples of disk- and slab-shaped clusters.

In this study, the dynamics of the system are given by a Metropolis Monte Carlo procedure with moves that keep the magnetization constant by only exchanging spin positions. We distinguish between local (Kawasaki Kawasaki (1966)) dynamics where only exchanges of neighboring spins are allowed and non-local dynamics that allow arbitrary exchanges. Local dynamics, which mimic mass transport in a more realistic system, are used in all simulations where dynamic properties are probed, while non-local moves are used for free energy calculations to enhance sampling. The unit of time $t$ is a Monte Carlo sweep, i.e., $N$ attempted single spin-exchange moves.

In the following discussion the reduced temperature is set to $\theta=T/T_{\text{c}}=0.7$ and the magnetization to $m=0.36$, where the slab state is metastable with respect to the disk state.

We will now introduce the order parameters $d$ and $s$, that we will use to describe transitions between the disk- and the slab state. For reference, Fig. 1 shows examples of configurations together with their respective $d$- and $s$-values.

II.1 The minimum distance/width order parameter $d$

We first introduce an order parameter that is capable of distinguishing between disk- and slab-shaped clusters as well as between intermediate shapes. To do so, we identify the largest (geometricSchmitz et al. (2013); Binder and Virnau (2016)) cluster of spins, i.e. the largest connected cluster of spins that are neighbors of each other and point down. If this cluster is a non-spanning cluster, we characterize its shape by finding the shortest distance, $d_{\text{image}}$, between the surface of the cluster and the surfaces of any of its periodic images (see Fig. 2). If the cluster is a spanning one, we instead characterize constrictions along its length by finding the smallest distance between the surfaces of the delimiting cluster of up spins, which is shown in gray on the right side of Fig 2. This quantity is denoted by $d_{\text{width}}$.

The combined minimum distance/width coordinate $d$ is then given by

[TABLE]

where the number $3$ is added to -$d_{\text{width}}$ in order to close a gap in the possible values of $d$ that is due to the lowest possible $d_{\text{image}}$ and $d_{\text{width}}$ being $\sqrt{2}$ (when the closest spins of the two clusters are offset by one diagonal step) and $2$ (when there is a bridge of exactly one spin left), respectively. Consequently, for spanning clusters $d$ is less than or equal to $1$ and for non-spanning clusters $d$ is larger than $1$. Furthermore, the order parameter $d$ differentiates between different cluster shapes both for spanning and non-spanning clusters and yields information on the local structure of the cluster around the connecting bridge to its periodic images.

The distance order parameter $d$ has units of one lattice constant.

II.2 The shape order parameter $s$

As we will see later, in order to arrive at a complete picture of the dynamics of the transition a second order parameter that characterizes the overall shape of the cluster is required. A straightforward way to do that for non-spanning clusters would be to use the ratio of the main moments of inertia of the cluster. However, this method breaks down as soon as the cluster is a spanning cluster since the main moments of inertia are no longer well defined.

Instead, we will use a proxy in Fourier space, where, for a system with $n\times n$ spins, we first take a Fourier transform of the configuration yielding the components

[TABLE]

and then take the (dimensionless) ratio

[TABLE]

Here, $\tilde{\sigma}_{10}$ and $\tilde{\sigma}_{01}$ are the longest wavelength components in $x$ and $y$ direction, respectively.

The order parameter $s$ takes on values from close to [math]—indicating configurations that are asymmetric with respect to an exchange of the $x$ and $y$ directions—to $1$ indicating that the largest cluster in the system is close to being symmetric. Its values are highly correlated to the ratio of the main moments of inertia for non-spanning clusters (see Fig. 3), but it is also well defined for spanning clusters.

III Free energies and the committor

We start our investigation of the mechanism of the disk-to-slab phase transition by computing the free energy landscape $\beta F(d,s)=-\log\left[P(d,s)\right]$, where $P(d,s)$ is the equilibrium probability density as a function of $d$ and $s$, $\beta=1/k_{\text{B}}T$, and $k_{\text{B}}$ is Boltzmann’s constant. We do so by performing umbrella sampling simulations in the plane spanned by the $d$ and $s$ coordinates. Configurations are biased using a series of harmonic bias potentials that are centered at different points in the $d$-$s$ plane. The resulting probability distributions are then reweighted using the weighted histogram analysis method (WHAM) Ferrenberg and Swendsen (1989); Kumar et al. (1992). Figure 4 shows the result of these calculations. The slab state, $\mathcal{A}$, is metastable with respect to the disk state, $\mathcal{B}$, by $6.4~{}k_{\text{B}}T$ and the two states are separated by a barrier with a height of roughly $8~{}k_{\text{B}}T$ measured from the minimum within state $\mathcal{A}$.

In order to assess the quality of the two chosen order parameters $d$ and $s$, transition path sampling Geissler et al. (1999) (TPS) simulations using aimless shooting Peters and Trout (2006) and variable trajectory lengths have been used to obtain trajectories transitioning from state $\mathcal{A}$ to state $\mathcal{B}$. Configurations are assumed to have reached the basins if the value of $d$ is smaller than $-15$ or greater than $26$, respectively. The trajectories consist of segments with a length of $2000$ MC sweeps each and the shooting point is randomly shifted forward or backward in time by $20$ segments or the same shooting point as in the previous trajectory is chosen. Figure 6 shows example trajectories obtained from these simulations.

As was mentioned in the introduction, a good reaction coordinate is able to characterize the progress of the transition. This means, that solely based on this coordinate one is able to predict the so called committor probability. The committor probability, $p_{\mathcal{B}}$, of a configuration is the probability of reaching state $\mathcal{B}$ before reaching state $\mathcal{A}$ averaged over trajectories started from the configuration. In our case, the averaging is done over trajectories computed with different random number generator seeds.

We computed committor values for configurations from pathways harvested in TPS runs. We estimated $p_{\mathcal{B}}$ by shooting trajectories from a given configuration and counting the number of times these trajectories reach state $\mathcal{B}$ first. The number of shots was chosen such that the statistical uncertainty of the estimated probability is smaller than $10\%$ Dellago et al. (2002).

Figure 5 shows the resulting committor probabilities together with the free energy $\beta F(d,s)$. The shape of the free energy barrier suggests that the committor probability is determined by the value of the $d$-coordinate only. However, an examination of the committor probabilities calculated from simulation paints a different picture, in which there is a significant dependence of the committor on the value of the $s$-coordinate. In particular, some configurations that, based on the free energy landscape, seem to be in the basin of attraction of state $\mathcal{A}$ ($\mathcal{B}$) have a greater (smaller) than $1/2$ chance to evolve into the disk state $\mathcal{B}$. Such configurations, marked with $A$ and $B$, are shown in Fig. 8.

At this point one may suspect that we have not yet captured all the relevant degrees of freedom that determine the progress of the transition. However, Fig. 7—which shows an optimized coordinate obtained using a likelihood maximization method due to Peters and Trout (2006)—suggests that the committor can be fairly accurately described using only $d$ and $s$. In the following we will reconcile this apparent disconnect by including the dynamics along $d$ and $s$ into our considerations.

We start by noting that in the trajectories shown in Fig. 6, the system spends a much larger amount of time close to the free energy barrier than it takes to cross the barrier to the other side. The $s$-coordinate evolves slowly because a large number of spins need to change collectively in order to change the shape parameter $s$ appreciably. Hence, for the dynamics chosen for the system, changing $s$ corresponds to the diffusive rearrangement of spins over long distances. In contrast, close to the barrier region, a change in the $d$-coordinate only requires the movement of relatively few spins and hence changes along the $d$-direction proceed much faster. As we will see later, this difference in relaxation times along the $d$- and $s$-direction strongly influence the distribution of committor probabilities. We will characterize this dynamical behavior by measuring the coordinate dependent diffusion coefficient tensor $D_{ij}(\vec{x})$, where $\vec{x}=(d,s)$.

IV Characterizing System Dynamics

In order to determine $D_{ij}$, we first employ the method of Im and Roux Im and Roux (2002); Pan et al. (2008) to calculate mean square displacements (MSDs) along the two coordinate directions that are corrected (up to first order) for the drift of the order parameters that is due to the underlying free energy landscape. Starting from equilibrium configurations within narrow ranges of $d$- and $s$-values, we generate short, unbiased trajectories. Using the coordinate $d$ as an example, the MSD is then estimated by

[TABLE]

Here $\tilde{d}(t)=d(t)-d(0)$ and the averages $\left<\ldots\right>_{t}$ are taken over the trajectories at time $t$ from their start. This procedure takes the mean drift (given by $\left<\tilde{d}\right>_{t}$) , as well as the slightly different starting points of the trajectories into account.

Figure 9 shows the results of this calculation. In principle, linear fits to the data shown then yield the diffusion coefficients $D_{XX}$. While this is certainly true for the $D_{ss}$ component, the MSDs along the $d$-direction are highly non-linear even for long times, making it impossible to reliably determine the value for $D_{dd}$ this way.

Before we address this problem, note that Fig. 10 shows the correlations

[TABLE]

which are small and, hence, in the following we will treat the fluctuations in $d$- and $s$-direction to be independent of each other by assuming $D_{ij}$ to have the diagonal form

[TABLE]

The values for $D_{s}$ obtained by fitting to the MSDs are shown in the lower panel of Fig. 12.

In order to obtain reliable results for $D_{d}(\vec{x})$, we use a Bayesian method introduced by Best and Hummer Hummer (2005); Best and Hummer (2010), which is based on discretizing order parameter axes into bins followed by the maximization of the likelihood

[TABLE]

Here, $R$ is the matrix of transition rate constants between the bins, $t_{\alpha}$ is the lag time between two measurements of the order parameter and the product runs over all transitions of the system observed in a set of trajectories. The components of $R$ obey the conditions

[TABLE]

where $P_{i}$ is the equilibrium probability of finding the system in bin $i$. Assuming that only the transition rates between neighboring bins are non-zero (i.e., we assume that the system transitioning from bin $i$ to $j$ has to pass through all bins in between), leaves us with $\mathcal{N}-1$ independent components, where $\mathcal{N}$ is the number of bins used. The diffusion coefficient between bin $i$ and $i+1$ along a coordinate $X$ with a bin width of $\Delta X$, is then given by Hummer (2005); Bicout and Szabo (1998)

[TABLE]

The advantage of this method is that we can choose a lag time $\tau=t_{\alpha}$ that is long enough, such that the memory effects seen in the mean square displacements subside while we can still deal with the considerable drift that occurs within this lag time.

Since $D_{d}$ is considered to be independent of the $s$-direction and the evolution along the $s$-direction is slow, we can perform this calculation independently in parallel slices along lines of constant $s$. These slices have a width of $\Delta s=0.025$ and each of these slices is split into bins of size $\Delta d=2$. As suggested in Refs. Hummer, 2005; Best and Hummer, 2010, we count the number of transitions between these bins, $N_{ij}$, for a set of trajectories and subsequently vary the rate parameters $R_{ij}$ in order to maximize the log-likelihood

[TABLE]

where the second sum is a smoothing prior. We take the $P_{i}$ from previous free energy calculations. The matrix exponential $\exp\left(\tau R\right)$ is evaluated using a scaling-and-squaring method Al-Mohy and Higham (2010) as implemented in the python package scipy Jones et al. . The value of $\varepsilon$ was chosen to be $3\cdot 10^{-4}$, which is on the order of the eventual values obtained for $D_{d}$, in order to smooth large statistical fluctuations.

The lag time is chosen by varying the lag time until one observes a plateau in the resulting diffusion coefficients. Figure 11 shows the results of such a calculation as a function of lag time for a single slice. At $\tau=30\,000$ the values of the diffusion coefficients seem to have reached a plateau while statistics are still reasonably good. This lag time is used to calculate the values of $D_{d}$ shown in the upper panel of Fig. 12.

V Understanding the Committor

In this section we will demonstrate that the shape of the committor distribution can be largely explained on the basis of the order parameters we have already introduced, $d$ and $s$, if one includes the dynamics of the system. Based on the free energy landscape and the diffusion coefficients we have obtained, we will first identify and estimate the timescales involved in the transition in order to gain an intuitive understanding, followed by an approach that uses a Smoluchowski equation to model the process.

V.1 A qualitative explanation

We start by noting that the drift towards the stable basins along the $s$-direction is slow compared to the relaxation time along the $d$-direction. Hence, in the following calculation we assume that a local equilibrium along slices of constant $s$ develops and consequently crossing the barrier becomes a one-dimensional problem associated with a timescale $\tau_{\text{cross}}(s)$, similar to a previous discussion of multicomponent nucleation by Trinkaus (1983).

We then compare $\tau_{\text{cross}}$ to the timescale $\tau_{\text{drift}}$ that is associated with the diffusion along the $s$-direction; if $\tau_{\text{cross}}\gg\tau_{\text{drift}}$, the system is likely to stay on the same side of the barrier, whereas if they are comparable, one expects to see a finite probability of crossing to the other side of the barrier.

Given the free energy landscape and the diffusion coefficient $D_{d}(s)$, we can estimate $\tau_{\text{cross}}$ by

[TABLE]

where $F^{\dagger}(s)$ is the height of the free energy barrier found when one looks only at a slice of the free energy landscape along a line of constant $s$ and $\Delta d$ is a length scale given by the width of the basin as seen along the slice (see Fig. 13).

$\tau_{\text{drift}}$ is related to the drift of the system along the $s$-axis. As the system moves towards the stable basin and the barrier height $F^{\dagger}$ becomes larger, $\tau_{\text{cross}}$ grows exponentially as a function of $F^{\dagger}$. This change in barrier height can be used to estimate the time it takes until the system has moved to an area within the free energy landscape, where the rate of transitions has considerably changed relative to where the system has started from. We define a characteristic length $\Delta s$ as the distance at which the barrier height has increased by $\Delta F^{\dagger}=1~{}k_{\text{B}}T$. $\Delta s$ is then defined by setting

[TABLE]

One can now estimate $\tau_{\text{drift}}$ as the time it takes the system on average to undergo a change of $\Delta s$:

[TABLE]

Here we have chosen $\dot{s}$ to be the mean velocity due to the drift caused by the slope in the free energy, i.e. $\dot{s}=-D_{s}\beta\frac{\partial F(s)}{\partial s}$. The ratio of the two timescales, given by

[TABLE]

is shown in Fig. 14.

One can see that in the outer regions with $s\lesssim 0.25$ and $s\gtrsim 0.45$ the crossing time becomes comparable to the drift time. Furthermore, if $\tau_{\text{cross}}/\tau_{\text{drift}}$ is close to $1$ in one direction of the transition, it is orders of magnitude larger in the reverse direction. This suggests that in these regions the probability of crossing from one side of the barrier to the other is non-zero and, on the other hand, after the barrier has been crossed, the probability of going back is small. Hence, we arrive at the explanation of why the committor probabilities in these regions take on the values shown in Fig. 5. It takes a long time to change the overall shape of the cluster and in this time fluctuations can occur that make the system cross the barrier. However, once the barrier has been crossed, the probability of going back the other way is small because the global shape of the cluster disfavors such fluctuations. In other words, there is a high probability that the bridge of spins that can be seen in the right panel of Fig. 8 is severed before the cluster has changed its shape into a slab and, conversely, the likelihood of a bridge forming in the left panel before the cluster becomes disk-shaped is high.

V.2 Smoluchowski analysis

A more accurate analysis of the system’s behavior given the free energy landscape and the diffusion coefficients can be achieved by modeling the dynamics of the system using a two-dimensional Smoluchowski equation. A similar approach has previously been employed by Metzner et al. (2006) in order to predict the behavior of a two-dimensional toy model. We use a generalized ansatz that takes into account state dependent diffusion coefficientsRisken (1984); Gardiner and Others (1985). In this case, the probability density $\rho(\vec{x},t)$ satisfies the equation

[TABLE]

where $\vec{x}=\left(d,s\right)$, $F=F(\vec{x})$ is the free energy shown before, and $D_{ij}(\vec{x})$ are position dependent diffusion coefficients. This form of the Fokker-Planck equation111The Smoluchowski equation (16) is often rewritten in the form $\frac{\partial\rho}{\partial t}=\sum_{i}\frac{\partial}{\partial x_{i}}\left[\sum_{j}D_{ij}e^{-\beta F}\frac{\partial}{\partial x_{j}}\left(e^{\beta F}\rho\right)\right]$, that makes it immediately obvious that the equilibrium distribution is the stationary solution. mimics a random walk on a potential of mean force given by the free energy $F$ and, at the same time, guarantees that the steady state solution is given by the equilibrium distribution $e^{-\beta F}$.

In this model the committor $p_{\mathcal{B}}(\vec{x})$ is then determined by the corresponding backward (or adjoint) equation Gardiner and Others (1985); Metzner et al. (2006)

[TABLE]

The reflecting boundary condition

[TABLE]

where $\vec{n}$ is a normal vector on the domain, and the boundary conditions that surround regions $\mathcal{A}$ and $\mathcal{B}$

[TABLE]

close the set of equations.

This system of equations can be solved using a finite difference scheme similar to the one employed in Ref. Metzner et al., 2006 using the free energy and diffusion coefficients measured from simulation as input. The result of such a calculation is shown in Fig. 15. It recovers the basic features of the committor observed in the brute force simulations. In particular, the $p_{\mathcal{B}}=1/2$ committor line significantly deviates from the ridge of the free energy landscape and the regions close to the barrier that have large or small values of $s$ have a significant probability of evolving to the other side of the free energy barrier.

We can also write down a corresponding Langevin-equation, whose realizations will—if evaluated using Ito’s interpretation (i.e. by evaluating the drift and diffusion terms at the beginning of each step)—evolve according to Equ. (16)Tupper and Yang (2012). It is given by

[TABLE]

where $\vec{x}(t)$ is the position of a random walker, $D_{ij}$, $F$ and their derivatives are evaluated at $\vec{x}(t)$, $\eta$ is delta-correlated white noise with zero mean and unity variance (i.e. $\left<\eta_{i}\right>=0$, $\left<\eta_{i}(t)\eta_{j}(t^{\prime})\right>=\delta(t-t^{\prime})\delta_{ij}$), and $\sum_{k}g_{ik}g_{kj}=2D_{ij}$Risken (1984). The latter condition simplifies to $g_{ij}=\sqrt{2D_{i}}\delta_{ij}$ if we assume the diffusion coefficient is diagonal and given by Equ. (7). Figure 16 shows an example Langevin trajectory compared to a trajectory obtained by TPS and Fig. 17 shows the probability distribution of finding a random walker as a function of time.

These two diagrams demonstrate two things: the trajectories generated using a full simulation and the ones generated from the Langevin equation very much look alike and, secondly, that indeed a relaxation towards equilibrium along the $d$ coordinate occurs before the system has time to evolve to the stable basins.

VI Discussion

We have investigated the disk-to-slab transition in the 2d-Ising model using a variety of methods, including a model based on a two-dimensional Smoluchowski equation. This model allows us to take the different timescales involved in this process into account: the fluctuations of the narrow bridge between the clusters (characterized by the $d$-coordinate) and the change of the overall shape of the cluster (captured by the $s$-coordinate). By comparing committor distributions obtained by brute-force using the full dynamics of the system, to the committor distribution obtained by solving the adjoint equation (17) (Fig. 15) we find that the dynamics of the transition are well described by the Smoluchowski equation (16).

This suggests that, with the coordinates $d$ and $s$, we have indeed captured the degrees of freedom that are most relevant for the disk-to-slab transition. Increased values of $s$, corresponding to disk-like clusters, enhance the likelihood of a trajectory started from a configuration to develop towards the disk state. While stated in this way, this result does not seem surprising, it is nonetheless interesting to note, that this effect only becomes apparent when one includes the dynamics of the system into consideration, while the free energy landscape alone does not capture this feature. In the disk-to-slab transition considered here, this is caused by the large difference in timescales associated with the crossing of the free energy barrier and the relaxation towards the stable basins that causes a local quasi-equilibrium to form where the global shape of the cluster is fairly stable, while its surface fluctuates.

Our analysis highlights the fact that, even in the context of a fairly simple model system, the dynamics of a system are crucial to identifying transition mechanisms, especially when one needs to use a combination of multiple order parameters to obtain a good reaction coordinate. In particular, including dynamics may help to reconcile free energy landscapes and committor probabilities, or their combination may point to the existence of other important degrees of freedom.

Acknowledgements.

The authors thank Georg Menzl and Phillip L. Geissler for many enlightening discussions. C.M. is supported by an uni:docs fellowship of the University of Vienna. A.T. acknowledges support from the Austrian Science Fund (FWF) Project Nr. P27738-N28. The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC).