Integrating Langevin Dynamics equations in computer simulation

Question

In a Langevin Dynamics simulation the following equation is solved numerically : $$ m_i\frac{d^2r_i}{dt^2}=F_{int}-\gamma\frac{dr_i}{dt}+R(t)$$

$$\langle R(t)\rangle=0 \quad \quad \langle R(t)R(t')\rangle=\sqrt{2k_BT\gamma m_i}\delta(t-t') $$ Where, $T$ is the target temperature of the simulation. Can this equation be integrated by using Velocity-Verlet scheme ? I have seen in the literature that there are many integrators which have been specifically developed to integrate this equation of motion.
To cite two examples :
https://doi.org/10.1016/j.cplett.2006.07.086, https://doi.org/10.1063/1.4802990.
Are there any restrictions in using ordinary integrators like Velocity-Verlet, leapfrog, Predictor-corrector which are normally used to integrate ordinary differential equations ?

Answer 1

Because the $R$ term is stochastic, there isn't "a" solution you obtain when numerically integrating this; you actually need to integrate the system a multitude of times for many different paths in order to generate a distribution of possible paths to time $t$. So the focus for solving SDEs is on the stochastic component, finding ways to stabilize integrations/solutions depending on the form of the coefficients (e.g., constants, dependent on $x$, etc).
The other methods you mention, verlet & predictor-corrector, have different focuses that are important to ODEs; for instance, verlet methods are focused on conserving energy during the evolution.

Given a generic stochastic differential equation (SDE), $$\mathrm{d}x=\mu\mathrm{d}t+\sigma\mathrm{d}W,$$ you can actually use the simple Euler scheme for SDEs, called the Euler-Maruyama method, which looks very much like the Euler scheme but with a stochastic term, e.g.

for path = 1 to num_paths
    dt = 1 / num_steps
    time = 0
    x[0] = 0
    for i = 1 to num_steps
        x[i] = x[i-1] + dt * mu + W(dt) * sigma
        time += dt
    // averaging, plotting, etc..

where $W(dt)=\varphi(0, \sqrt{dt})$ is the random normal distribution with 0 mean and variance $dt$.

Other methods generally include additional terms in the stochastic component, such as the Milstein method, which replaces the update of x[i] in the above as,

dW = W(dt)
x[i] = x[i-1] + mu(x,t) * dt + sigma(x,t) * dW 
              + 0.5 * sigma(x,t) * sigma_deriv(x,t) * (dW**2 - dt)

where sigma_deriv is $\partial_x\sigma$. The SDE version for the Runge-Kutta method also is similar to the above. If $\sigma$ is a constant, then these two reduce to the Euler-Maruyama method.

Peter Kloeden has a very good book on the subject of numerical SDEs, it might be worth the time investment in studying from it.