Classical Langevin dynamics

Langevin dynamics can be used to sample the canonical ensemble for a classical system. Langevin dynamics are based on classical equations of motion that are modified by an additional drag force and a random force. The Langevin equation of motion can be written as

\[\mathbf{M}\ddot{\mathbf{R}} = - \nabla_R V(\mathbf{R}) + \mathbf{F}(t) - \gamma \dot{\mathbf{R}}\]

where $\mathbf{M}$ are the masses, $\ddot{\mathbf{R}}$ the time-derivative of the positions, i.e., the velocities, $\nabla_R V(\mathbf{R})$ the gradient of the potential and $\mathbf{F}(t)$ the random force that is related to the friction coefficient $\gamma$ by the second fluctuation-dissipation theorem.

Equally the above equation can be written in the form of Ito stochastic differential equations [7]

\[\begin{aligned} d\mathbf{R} &= \mathbf{M}^{-1} \mathbf{P} dt\\ d\mathbf{P} &= [-\nabla V(\mathbf{R}) - \gamma \mathbf{P}] dt + \sigma \mathbf{M}^{1/2} d\mathbf{W} \end{aligned}\]

where $\sigma = \sqrt{2\gamma/\beta}$ and $\mathbf{W}$ is a vector of $N$ independent Wiener processes. As usual, $\mathbf{P}$ is the vector of particle momenta and $\mathbf{M}$ their diagonal mass matrix.

As a stochastic differential equation, these two can be integrated immediately using StochasticDiffEq provided by DifferentialEquations, which offers a variety of stochastic solvers. It is possible to exploit the dynamical structure of the differential equations by splitting the integration steps into parts that can be solved exactly. In this context, it has been shown that the BAOAB method from [7] achieves good accuracy compared to other similar algorithms and this algorithm is used here as the default.

Example

Using Langevin dynamics we can sample the canonical ensemble for a simple harmonic oscillator and investigate the energy expectation values.

Firstly we set up our system parameters. Here, we have two atoms in a harmonic potential at a temperature of 1e-3. We have arbitrarily chosen the dissipation constant $\gamma = 1$, this can be tuned for optimal sampling in more complex systems.

using NQCDynamics
using Unitful

atoms = Atoms([:H, :C])
temperature = 1e-3
sim = Simulation{Langevin}(atoms, Harmonic(m=atoms.masses[1]); γ=1, temperature)

Simulation{Langevin{Float64}}:
  Atoms{Float64}([:H, :C], [1, 6], [1837.4715941070515, 21894.713607956142])
  Harmonic{Float64, Float64, Float64}
  m: Float64 1837.4715941070515
  ω: Float64 1.0
  r₀: Float64 0.0
  dofs: Int64 1

Here we can generate a simple starting configuration with zeros for every degree of freedom.

u = DynamicsVariables(sim, zeros(size(sim)), zeros(size(sim)))

([0.0 0.0], [0.0 0.0])

Running the dynamics proceeds by providing all the parameters along with any extra keywords. This time we have requested both the positions and velocities to be outputted and have selected a timestep dt. Since the default algorithm is a fixed timestep algorithm an error will be thrown if a timestep is not provided.

traj = run_dynamics(sim, (0.0, 2000.0), u; output=(OutputPosition, OutputVelocity), dt=0.5)

3-element Dictionaries.Dictionary{Symbol, Any}
           :Time │ [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5  …  1995.5…
 :OutputPosition │ [[0.0 0.0], [-0.000139940492310738 3.567948941707327e-6], [-…
 :OutputVelocity │ [[0.0 0.0], [-0.0005247768461652675 1.419693745916881e-5], […

Here, we plot the positions of our two atoms throughout the simulation.

using Plots
plot(traj, :OutputPosition, label=["Hydrogen" "Carbon"], legend=true)

We next plot the velocities. Notice how the carbon atom with its heavier mass has a smaller magnitude throughout.

plot(traj, :OutputVelocity, label=["Hydrogen" "Carbon"], legend=true)

Using the configurations from the Langevin simulation we can obtain expectation values along the trajectories. This can be done manually, but we provide the Estimators module to make this as simple as possible.

!!! note Estimators

[Here](@ref `Estimators`) you can find the available quantities that [`Estimators`](@ref) provides.
To add new quantities, you must implement a new function inside `src/Estimators.jl`.

Let's find the expectation for the potential energy during our simulation. This is the potential energy of the final configuration in the simulation:

julia> Estimators.potential_energy(sim, traj[:OutputPosition][end])1.0584545629185378e-6

We could evaluate this for every configuration and average it manually. Fortunately however, we have the @estimate macro that will do this for us:

julia> Estimators.@estimate potential_energy(sim, traj[:OutputPosition])0.0008714360199255442

Similarly, we can evaluate the kinetic energy expectation with:

julia> Estimators.@estimate kinetic_energy(sim, traj[:OutputVelocity])0.0009448229952026531

Again, this takes a similar value since the total energy is evenly split between the kinetic and potential for a classical harmonic system.