Final Project: Simulating Chemical Reactions with Langevin Dynamics

Scientific Motivation**¶

Many chemical and biological reactions involve crossing energy barriers due to thermal fluctuations. These processes can be understood as stochastic dynamics over a potential energy landscape, where the system hops from one stable state (reactant) to another (product).

In this project, you’ll simulate a 1D model reaction using Langevin dynamics:

• The particle represents a reactant molecule.

• The potential energy landscape represents the reaction coordinate, e.g., reactant → transition state → product.

• The particle’s crossings of the barrier represent reaction events.

This is a numerical experiment in thermal activation and a testbed for rate theory.

Project Objectiv¶

Simulate a 1D chemical reaction using Langevin dynamics, estimate the reaction rate constant using multiple approaches, and compare your results to Kramers’ theory.

Tasks¶

1. Define the Potential Energy Surface¶

Use a double-well potential to represent a reaction:

$V(x) = a(x^2 - 1)^2$

• Minima at $x = \pm 1$ (reactant and product)

• Barrier at $x = 0$, height $\Delta V = a$

• Example: $a = 2$ gives $\Delta V = 2$

You can also try asymmetric wells or multi-barrier landscapes later.

2. Implement Langevin Dynamics¶

Use the Langevin equation for a particle in a potential:

$m \ddot{x} = -\gamma \dot{x} - \frac{dV}{dx} + \sqrt{2 \gamma k_B T} \, \eta(t)$

• $\eta(t)$ is Gaussian white noise

• Use BAOAB scheme (or Euler-Maruyama for simplicity)

Parameters:

• Mass $m = 1$

• Friction $\gamma = 1.0$

• Temperature $T = 0.25 \text{ to } 1.0$

• Time step $dt = 0.01$

• Simulate long enough to observe many transitions

3. Label States and Count Transitions¶

• Define:

  • Left well: $x < 0$

  • Right well: $x > 0$

• Count the number of times the particle switches wells

• Compute:

$k = \frac{\text{Number of transitions}}{\text{Total time}}$

This is your empirical rate constant.

4. Compare to Kramers’ Theory¶

Kramers’ estimate of the rate constant is:

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

• $\omega_0$: curvature at the well bottom (from second derivative of $V(x)$)

• $\omega_b$: curvature at the barrier top

• $\Delta V$: barrier height

Compare your simulated $k$ to this analytical expression at various temperatures.

5. Correlation Function Method¶

Define the population indicator:

$h(x) = \begin{cases} 1, & x > 0 \\ 0, & x < 0 \end{cases}$

Compute the autocorrelation function:

$C(t) = \langle h(0) h(t) \rangle \sim e^{-kt}$

Fit $C(t)$ to extract the rate constant.

Challenges (Optional)¶

• Try asymmetric double-well: $V(x) = a(x^2 - 1)^2 + \delta x$

• Compute the mean first passage time from one well to the other

• Vary friction $\gamma$ to observe turnover behavior in Kramers’ theory

• Simulate in 2D or on rugged landscapes

Learning Outcomes¶

• Understand the connection between stochastic dynamics and chemical reaction rates

• Explore rare event sampling and its challenges

• Learn to extract kinetic information from trajectories

• See how temperature and friction affect barrier crossing

💡 Starter Parameters¶