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#

Parameter

Value

Potential

\( a = 2.0 \)

\( \gamma \)

1.0

\( T \)

0.25–1.0

\( dt \)

0.01

\( t_{\text{total}} \)

\( 10^4 \)–\( 10^6 \) steps