A Computational Model: Dimer Isomerization

To show the power of transition path sampling, it is perhaps best to consider a simple example of a two-state system where the states are separated by an energy barrier. This presentation is based on an exercise from the molecular simulation course taught by Berend Smit and Daan Frankel in 2001, and uses a code called tps.

Here, we consider dimer isomerization in a simple 2D liquid sample of $N$ 15 particles confined to a circle of radius $R$. Particles 1 and 2 form the dimer. The total potential energy is given by

$$\mathscr{U} = {\sum_{i<j}}^\prime U_{\rm WCA}\left(r_{ij}\ri... ...sum_i U_{\rm wall}\left(r_i\right) + U_{12}\left(r_{12}\right)$$ (306)

The prime indicates that the 1-2 pair is excluded from this sum. All particles interact with each other via a modified Lennard-Jones pairwise interaction known as the Weeks-Chandler-Andersen (WCA) potential:

$$U_{\rm WCA}\left(r\right) = \left\{\begin{array}{ll} 4\left... ...r < 2^{1/6}$}\\ 0 & \mbox{$r > 2^{1/6}$} \end{array}\right.$$ (307)

The WCA potential is fully repulsive, and cutoff where the force vanishes. The particles interact with the circular wall via another WCA potential:

$$U_{\rm wall}\left({\bf r}\right) = U_{\rm WCA}\left(R+2^{1/6}-r\right)$$ (308)

where ${\bf r}$ is the vector position and $r$ is the radial position (i.e., the circle is centered on the origin). The dimer bond (between particles 1 and 2) is described by a two-well potential:

$$U_{12}\left(r_{12}\right) = h\left[1-\frac{\left(r_{12}-w-2^{1/6}\right)^2}{w^2}\right]^2$$ (309)

This potential has stable minima at $r_{12} = 2^{1/6}$ and $r_{12} = 2^{1/6}+2w$. Here, $h$ defines the height of the potential energy barrier between two stable states as defined by this potential, and $w$ defines the width of this barrier:

The WCA and dimer potentials used in this case study.

For now, we fix $w$ at 0.25, and $R$ at 3.0. This gives an areal particle number density $\rho = 0.53 \sigma^{-2}$.

What is a reasonable order parameter for this system? The bond length, naturally. All configurations for which $r_{12} < 1.30$ are defined to be in region $A$, and all configurations for which $r_{12} > 1.45$ are defined be in region $B$. It is important to note here that $h$ is not a measure of the free energy barrier which determines the kinetic rate in this process. To determine the free energy barrier, we need to construct the free energy profile, $F(r_{12}) = -k_BT\ln Q(r_{12})$. This type of restricted free energy is often termed a Landau free energy.

Let's consider first a small potential barrier height, $h = 2$. For this height, we can compute $P(\lambda,t\rightarrow\infty)$ very accurately from a single long MD simulation (10,000,000 steps; $\Delta t$ = 0.001; total energy 15 $\epsilon$). We see from the linear region of $C(t)$ that the rate constant appears to be $k_{AB} = 4.04\times10^{-2} \tau^{-1}$. Now, let us perform transition path sampling MC on this system.

$C(t)$ computed from a long MD simulation for barrier height $h$ = 2. The total energy is set at 15 $\epsilon$.

First, let's decide how long an interval $T$ we need. Clearly, $C(t)$ appears to be linear up to $t = 2$; following Dellago, we'll choose $T$ = 2.0. ⁴ Now, we need to conduct two types of MC simulations:

1. A single path sampling MC simulation to compute $\left<\dot{h}_B(t)\right>_{AB}$ and $\left<h_B(t^\prime)\right>_{AB}$ for $0 < t,t^\prime < T$. We will conduct a simulation of $10^6$ cycles.
2. A set of umbrella sampling MC simulations on the following intervals for $r_{12}$:

   $i$	min	max
   1	0.00	1.22
   2	1.20	1.26
   3	1.24	1.30
   4	1.28	1.45
   5	1.40	2.45

   When matching the histograms in each of these intervals and renormalizing, we obtain $P\left(r_{12},t=2\right)$. When we then integrate from $r_{12} = 1.45$ to $r_{12} = 2.45$ (region $B$).

We'll use the other default parameter values provided in tps. These include the maximum angle by which velocity vectors are rotated in a shooting move ( $\Delta\phi = 0.6$ rad), and the number of shifting moves per shooting move (95:5), and the number of $r_{12}$-slices per window (200).

We see that $\left<\dot{h}_B(t)\right> = 0.22$ and $\left<h_B(2)\right> = 0.55$.

$\left <h_B(t)\right >$ computed from a transition path sampling MC simulation for barrier height $h$ = 2. The total energy is set at 15 $\epsilon$.

Now, the histogram:

We see that $\left<\dot{h}_B(t)\right> = 0.22$ and $\left<h_B(2)\right> = 0.55$.

$P(r_{12},2)$ computed from a transition path sampling MC simulation and umbrella sampling, for barrier height $h$ = 2. The total energy is set at 15 $\epsilon$.

Integrating over region $B$ yields $C(2)$ = 0.09037. Performing the requisite operations yields $k_{AB} = 0.0360 \tau^{-1}$, which compares well to the MD-calculated value $k_{AB} = 4.04\times10^{-2} \tau^{-1}$.

Repeating this entire procedure for $h$ = 6 yields $k_{AB} = 5.05\times 10^{-4} \tau^{-1}$.

[to be completed; runs currently executing, March 13, 2013]