Langevin equation in 4th order Runge-Kutta

Question

I'm trying to figure out how to translate a piece of code from Velocity Verlet to Runge-Kutta, while treating the time step dependence of the thermal noise correctly.

The Langevin equation for my system reads reads

$$ ma = - \gamma v - \frac{dU}{dx} + \xi(t), $$

where $U$ is some interaction potential, $\gamma$ is damping, and $\xi$ is a gaussian noise term with $\mu = 0$ and $\sigma = \sqrt{2\gamma m k_B T}$.

I can use Velocity Verlet with Langevin dynamics as

$$ v_{t+1} = v_t + h(a - \gamma v + \xi(t)). $$

Qualitatively, what the Langevin equation does here is that it models thermal fluctuations by adding random kicks to the acceleration while counteracting them with a constant damping term to stabilize the energy. My question then is, how does this translate to 4th order Runge-Kutta (RK4)?

In RK4 we calculate the velocity as

$$ v_{t+1} = v_t + \frac{h}{6}(a_1 + 2a_2 + 2a_3 + a_4) $$

where $a_i$ are the partial accelerations calculated in the RK4 steps.

It is not obvious to me where to introduce the Langevin dynamics here. My best guess is that it should be applied in every separate RK-step? Meaning e.g. for $a_1$

$$ a_1 = a_t - \gamma v_t + \xi(t). $$

Of course we would have to use the same $\xi(t)$ for all the $a_i$ during one time step for this to make sense. Meaning we generate one $\xi(t)$ at the start of every time step that we then use for every calculation during that time step.

Still, something is missing here... Because now the noise is not dependent on the time step, and it should be! This was not an issue in Velocity Verlet because we just multiplied the noise term with $h$ during every time step, but this is not the case here. It seems to me that the time step has to be included somewhere in the $\sigma$ term of the Langevin equation, but I can't really figure out how...

edit1: changed to a more sensible notation.

edit2: I realized, for RK4 to work, you probably have to add a timestep to the noise term as $\sigma = \sqrt{\frac{2\gamma m k_B T}{h}}$ for the units to come out correctly.

Answer 1

Caveat: This is not a full answer to your question on how to replace VV by RK4. Basically, I recommend another integrator method which also has potentially fourth order accuracy.

You want to integrate the above Langevin equation over time. As an analytical solution is not known, you need to do it numerically. As @Chris is correctly pointing out, this is a stochastic differential equation (SDE). Hence, your time integration consists of two parts: a deterministic part (momentum and position) and a stochastic part (affecting momentum).

One common approach is to use a splitting method: This splits the above Langevin equation (rather the associated Fokker-Planck operator) into distinct parts, each of which can be solved analytically. These analytical solutions are then used as basic building blocks from which numerical integrators of specific accuracy ("order") can be constructed.

I recommend reading the paper by Leimkuhler et al., 2012 on the splitting methods and the resulting basic building blocks -- coined "A", B" and "O" -- from section 2, especially section 2.1 on timestepping methods. (Note that the article focuses on molecular sampling and statistical physics in its terminology but this of no relevance.) I believe that this paper will give you some background on how numerical time integration schemes (for SDEs) are obtained.

So far, you have used Velocity-Verlet with Langevin. Hence, in the notation of the paper you are probably doing "BAO" (B = momentum integration, A = position integration, O = stochastic integration).

If you are only interested in obtaining a more accurate integration scheme, then you could for example use "BAOAB" (see the above paper) which has second order accuracy in the momenta and even up to fourth order in the positions and some other nice features.

The bottom line is: Simply replacing (in the deterministic integration part) VV by RK4 will not help, if you do not treat the stochastic part correctly.