Barrier Crossing and Kramers’ Rate from Langevin Dynamics

The model: a particle in a double well¶

Consider a 1D reaction coordinate with potential

V(x) = a (x^2 - 1)^2, \qquad V'(x) = 4 a x (x^2 - 1).

(1)

The two minima at $x = \pm 1$ are reactant and product states; the barrier at $x = 0$ has height $\Delta V = a$ . Motion is governed by the Langevin equation

m \ddot x = -\gamma\, \dot x - V'(x) + \sqrt{2 \gamma\, k_B T}\, \eta(t),

(2)

with $\langle \eta(t) \eta(t') \rangle = \delta(t-t')$ . We integrate with the BAOAB splitting (Leimkuhler & Matthews), which is accurate for configurational averages even at large $\Delta t$ .

Kramers’ high-friction rate (moderate-to-large $\gamma$ ):

k_{\text{Kramers}} = \frac{\omega_0\, \omega_b}{2 \pi \gamma}\, e^{-\beta \Delta V},

(3)

where $\omega_0^2 = V''(x=\pm 1)/m = 8a/m$ and $\omega_b^2 = |V''(0)|/m = 4a/m$ .

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)

BAOAB Langevin integrator¶

def langevin_baoab(x0=-1.0, v0=0.0, T=0.5, gamma=1.0, m=1.0,
                   dt=0.01, n_steps=200_000, a=2.0, seed=0):
    r = np.random.default_rng(seed)
    x, v = float(x0), float(v0)
    traj = np.empty(n_steps)
    c1 = np.exp(-gamma * dt)
    c2 = np.sqrt((1.0 - c1 * c1) * T / m)
    force = lambda q: -4.0 * a * q * (q * q - 1.0)
    F = force(x)
    for i in range(n_steps):
        v = v + 0.5 * dt * F / m              # B
        x = x + 0.5 * dt * v                   # A
        v = c1 * v + c2 * r.standard_normal()  # O
        x = x + 0.5 * dt * v                   # A
        F = force(x)
        v = v + 0.5 * dt * F / m              # B
        traj[i] = x
    return traj

A trajectory with many barrier crossings¶

At $T = 0.5$ , $a = 2$ we have $\beta \Delta V = 4$ , so a few dozen transitions in 10⁵ steps — long enough to get a decent rate estimate but fast enough to run in seconds.

T, gamma, a, dt, n_steps = 0.5, 1.0, 2.0, 0.01, 200_000
traj = langevin_baoab(x0=-1.0, T=T, gamma=gamma, a=a, dt=dt, n_steps=n_steps, seed=1)
t = np.arange(n_steps) * dt

fig, ax = plt.subplots(figsize=(9, 3.5))
ax.plot(t, traj, lw=0.5, color='steelblue')
ax.axhline(0, color='k', lw=0.5)
ax.set(xlabel='time', ylabel='$x(t)$',
       title=f'Langevin trajectory ($T={T}$, $\\gamma={gamma}$, $a={a}$)')
ax.grid(True, ls=':', alpha=0.5)
fig.tight_layout()

Counting transitions, estimating $k$ ¶

We label each time step with its most recent committed well (reached $|x| > 0.8$ ) to avoid double-counting fast recrossings at the barrier top. For a symmetric well, the population is $1/2$ in each state, so the rate of escape equals the total transition rate:

k = \frac{\text{number of committed transitions}}{\text{total time}}.

(4)

Kramers’ high-friction formula $k_{\text{hi-}\gamma} = \omega_0 \omega_b / (2 \pi \gamma)\, e^{-\beta \Delta V}$ is only accurate when $\gamma \gg \omega_b$ . At moderate $\gamma$ the correct expression is

k_{\text{full}} = \frac{\lambda_+}{\omega_b} \cdot \frac{\omega_0}{2\pi} e^{-\beta \Delta V}, \qquad \lambda_+ = \sqrt{\gamma^2/4 + \omega_b^2} - \gamma/2,

(5)

which reduces to the high- $\gamma$ form when $\gamma \gg \omega_b$ . At $\gamma = 1$ we are in the intermediate regime and the two differ noticeably — Problem 2 explores the full turnover.

def count_transitions(x, commit=0.8):
    state = np.zeros(len(x), dtype=np.int8)
    state[x > commit] = 1
    state[x < -commit] = -1
    last = 0
    labels = np.zeros_like(state)
    for i in range(len(state)):
        if state[i] != 0:
            last = state[i]
        labels[i] = last
    diffs = np.diff(labels)
    return int(np.sum(diffs != 0))

n_trans = count_transitions(traj)
total_time = n_steps * dt
k_emp = n_trans / total_time          # counts both L->R and R->L; symmetric well

omega0 = np.sqrt(8 * a)               # at x=+-1
omegab = np.sqrt(4 * a)               # at x=0
dV = a
k_hi = omega0 * omegab / (2 * np.pi * gamma) * np.exp(-dV / T)          # high friction
lam = np.sqrt(0.25 * gamma**2 + omegab**2) - 0.5 * gamma
k_full = (lam / omegab) * (omega0 / (2 * np.pi)) * np.exp(-dV / T)      # full moderate-friction

print(f"transitions          : {n_trans}")
print(f"k (empirical)        : {k_emp:.4f}")
print(f"k (Kramers high-gam) : {k_hi:.4f}")
print(f"k (Kramers full)     : {k_full:.4f}")

Problems¶

Arrhenius plot. Run simulations at $T \in \{0.3, 0.4, 0.5, 0.6, 0.8\}$ (same $\gamma = 1$ , $a = 2$ ). Estimate $k(T)$ from transition counts and plot $\ln k$ vs. $1/T$ . Fit the slope and confirm it recovers the barrier $\Delta V = a = 2$ .
Kramers’ turnover. Fix $T = 0.5$ , $a = 2$ , and vary $\gamma \in \{0.05, 0.2, 1, 5, 20\}$ . Plot $k(\gamma)$ on log axes. You should see a non-monotonic curve: at small $\gamma$ the bath is too weak to reset momentum between crossings (energy-diffusion limited); at large $\gamma$ friction damps the motion (spatial-diffusion limited).
Asymmetric well. Add a tilt $-\delta x$ to the potential (so the right well is lower by $2\delta$ ). Estimate forward and backward rates separately. Verify detailed balance: $k_{+\to-}/k_{-\to+} = e^{-2\beta\delta}$ .

Reference¶

H. A. Kramers, Brownian motion in a field of force and the diffusion model of chemical reactions, Physica 7, 284 (1940). Kramers (1940)