This is an implementation of a Langevin dynamics simulator developed as part of my masters thesis. All figures, code are a work of my own.

Introduction

Consider a physical system \(\mathcal{S}\) with \(N\) molecules. Each of this molecule can occupy some portion of the Euclidian space \(q \in \Omega \subset \mathbb{R}^{3N}\) and have a velocity \(p \in \mathcal{T}\Omega \subset \mathbb{R}^{3N}\). These quantities are sometimes called the canonical coordinates of \(\mathcal{S}\). It is often true, that one is interested more in the average properties of this molecular system, such as the temperature, pressure, enthalpy; and less in the actual trajectories of these molecules. It suffices to say, if there exists a probability distribution over \(\Omega \times \mathcal{T}\Omega\), say \(P(q,p)\) then any macroscopic property of interest \(\Phi\), can be obtained from the following integral

\[\begin{align} \Phi := \mathbb{E}\left[\phi(q)\right] = \int P(q,p) \: \phi(q) \: dp \: dq \end{align}\]

where \(\phi(q)\) is a constitutive relation, that we know from prior experience. (More or less)

But the question of how one arrives at \(P(q,p)\) still remains to be answered. Typically, this is done by sampling sufficiently many realizations from this distribution using some Monte Carlo scheme. Solving the Langevin SDE leads to one such Monte Carlo method, whose trajectories sufficently samples \(P(q,p)\). In its most generic form, the Langevin equation is as follows.

\[\begin{align} d \begin{bmatrix} q \\ p \end{bmatrix} = \begin{bmatrix} M^{-1}p \: dt \\ -\nabla_q U(q) \: dt- \gamma p \: dt+ \sigma M^{-1/2} \: dW(t) \end{bmatrix} \end{align}\]

\(U : q \mapsto \mathbb{R}\) is the potential energy function. \(M, \gamma, \sigma\) are parameters used in tuning the simulation. \(W(t)\) is a Wiener process. Sometimes referred to as Brownian motion; for \(t,s,k,l \in [0,T]\) a Wiener process is defined as follows:

\[\begin{align} W(0) = 0 \\ \forall s>t \quad W(s)-W(t) \sim \sqrt{s-t}\:\mathcal{N}(0,1) \\ W(s)-W(t) \textrm{ and } W(l)-W(k) \textrm{ are indepenedent of each other.} \end{align}\]

For most practical purposes the Langevin SDE needs to be numerically simulated. This entails making time discretizations. It is usually true that the equations are stiff and small time discretization steps need to be taken for accurately simulating \(P(q,p)\). On the other hand, it could also be that smaller time steps, would require significantly longer computational runs before the entire distribution can be sampled. This “chicken and egg” dilemma, presents the need for numerical integrators that can efficiently simulate the Langevin equation without compromising accuracy.

This article discusses some of these numerical schemes, provides implementations of these methods and illustrates them on the classical harmonic and double well potentials.

Operator splitting

Operator splitting is a numerical approximation technique, where the vector field of an ordinary differential equation is complicated, but can be decomposed into an additive sum of terms that are trivial when treated individually. In the context of the Langevin SDE, this is

\[\begin{align} d \begin{bmatrix} q \\ p \end{bmatrix} = \begin{bmatrix} M^{-1}p \\ 0 \end{bmatrix} \: dt + \begin{bmatrix} 0 \\ -\nabla_q U(q) \end{bmatrix} \: dt + \begin{bmatrix} 0 \\ -\gamma p \: dt + \sigma M^{-1/2} \: dW(t) \end{bmatrix} \end{align}\]

Here, each of the terms on the right hand side, is tractable – some analytically and the one with \(U\) numerically.

Formally let the solution to

\[\begin{align} d \begin{bmatrix} q \\ p \end{bmatrix} = \begin{bmatrix} M^{-1}p \\ 0 \end{bmatrix} \: dt \end{align}\]

be

\[\begin{align} \phi_A(t) := \begin{bmatrix} q(t) \\ p(t) \end{bmatrix}_A = \begin{bmatrix} e^{M^{-1}t}q(0) \\ p(0) \end{bmatrix} \end{align}\]

And the solution to

\[\begin{align} d \begin{bmatrix} q \\ p \end{bmatrix} = \begin{bmatrix} 0 \\ -\gamma p \: dt + \sigma M^{-1/2} \: dW(t)\end{bmatrix} \end{align}\]

be

\[\begin{align} \phi_O(t) := \begin{bmatrix} q(t) \\ p(t) \end{bmatrix}_O = \begin{bmatrix} q(0) \\ e^{-\gamma t} p(0) + \frac{\sigma}{\sqrt{2 \gamma}} \sqrt{1-e^{-2\gamma t}} M^{-1/2} R(t)\end{bmatrix} \end{align}\]

where \(R(t) \sim \mathcal{N}(0,1)\).

Finally, there exists no analytical solution to

\[\begin{align} d \begin{bmatrix} q \\ p \end{bmatrix} = \begin{bmatrix} 0 \\ -\nabla_q U(q)\end{bmatrix} \: dt \end{align}\]

in the general case and needs to be simulated numerically. We call this solution

\[\begin{align} \phi_B(t) := \begin{bmatrix} q(t) \\ p(t) \end{bmatrix}_B \end{align}\]

The idea behind operator splitting is to compose these individual solutions in a specific sequence that results in minimum error. For this Langevin problem, Leimkuhler considered several sequences. We will discuss two of them.

BAOAB
\[\begin{align} \phi(t+\tau) := \phi_B \left(t+\frac{\tau}{2}\right) \circ \phi_A \left(t+\frac{\tau}{2}\right) \circ \phi_O(t + \tau) \circ \phi_A \left(t+\frac{\tau}{2}\right) \circ \phi_B \left(t+\frac{\tau}{2}\right) \end{align}\]
ABOBA
\[\begin{align} \phi(t+\tau) := \phi_A \left(t+\frac{\tau}{2}\right) \circ \phi_B \left(t+\frac{\tau}{2}\right) \circ \phi_O(t + \tau) \circ \phi_B \left(t+\frac{\tau}{2}\right) \circ \phi_A \left(t+\frac{\tau}{2}\right) \end{align}\]

Examples

In the following, we consider one dimensional molecular systems. That is \(q,p \in \mathbb{R}^N\). With these assumptions, the simulated trajectories and their distributions, along with the energies are shown below.

One body system

Here \(N=1\)

Simple Harmonic oscillator
\[\begin{align} U(q) = \frac{1}{2} k(q-c)^2 \end{align}\]
Double Well potential oscillator
\[\begin{align} U(q) = \frac{1}{4} k(q-a)^2(q+a)^2 \end{align}\]

Two body systems

Here \(N=2\), that is \(q := \begin{bmatrix} q_1 \\ q_2 \end{bmatrix}\)

Double well with pairwise harmonic potential
\[\begin{align} U(q_1) = \frac{1}{4} k(q_1-a)^2(q_1-b)^2 + \frac{1}{2} k(q_1-q_2)^2 \end{align}\]
Double well with soft-bound potential
\[\begin{align} U(q_1) = \frac{1}{4} k(q_1-a)^2(q_1-b)^2 + V(q_1, q_2) \\ V(q_1, q_2) = \begin{cases} cos\left(1 + \frac{ (q_2 - q_1) \pi }{r_c}\right) \quad r<r_c \\ 0 \quad r \geq r_c \end{cases} \end{align}\]

Extending these simulations to higher dimensions is trivial, but remains out of the scope of this article.

Code