Computing transport coefficients with Molly.jl

Transport coefficients: a computational challenge

• Relate to sensitivities of fluxes in response to nonequilibrium forcings (linear response).
• Dynamical properties of molecular systems (thermal transport, diffusion, shear viscosity…)
• Parametrize macroscopic evolution equations (e.g. Navier–Stokes).
• Standard methods require long simulations to get decent estimators of these quantities.
• Need for innovative variance reduction techniques (e.g. coupling methods, synthetic forcings).
• Molly is a promising tool for testing new methods, because the total cost of learning + implementing/debugging + running simulations is low.

Setting: nonequilibrium kinetic Langevin dynamics

• State space: \(\mathcal X = \mathcal D \times {\mathbb{{R}}}^{dN}\) (positions and momenta).
• Underdamped Langevin with external forcing: \[ \begin{equation} \left\{\begin{aligned} \mathrm{d}q_t^\eta & = \color{red}{M^{-1}p_t^\eta \mathrm{d}t}, \\ \mathrm{d}p_t^\eta &= \color{blue}{-\nabla V(q_t^\eta)\,\mathrm{d}t} \color{green}{-\gamma M^{-1}p_t^\eta\,\mathrm{d}t +\sqrt{\frac{2\gamma}{\beta}}\mathrm{d}W_t} +\color{purple}{\eta F(q_t^\eta)\,\mathrm{d}t}, \end{aligned}\right. \end{equation} \]
• Generator: \(\mathcal{L}_\eta = \mathcal{L}+ \color{purple}{\eta F\cdot \nabla_p}\) with \(\mathcal{L}= \color{red}{M^{-1}p\cdot \nabla_q} \color{blue}{-\nabla V\cdot \nabla_p} \color{green}{-\gamma M^{-1}p\cdot \nabla_p +\frac{\gamma}{\beta}\Delta_p}\)
• Parameter: \(\eta\) = strength of the perturbation.
• Steady state density: solution to the Fokker–Planck equation \(\mathcal{L}_\eta^\dagger \psi_\eta = 0\).
• Forcing: \(F\) not the gradient of a potential function on \(\mathcal D\).
• No explicit invariant measure: \(\psi_\eta(q,p)\,\mathrm{d}q\,\mathrm{d}p \neq Z_{\beta,\eta}^{-1}\mathrm{e}^{-\beta H_\eta(q,p)}\,\mathrm{d}q\,\mathrm{d}p\) except for \(\eta=0\).

Linear response

• Response observable or flux: \[R:\mathcal X \to\mathbb R.\]

• Assume: \[ \int_{\mathcal X}R(q,p)\,\mathrm{e}^{-\beta H(q,p)}\,\mathrm{d}q\,\mathrm{d}p =0 \] (No flux at equilibrium)

• Linear response: \[ \begin{equation}\rho(F,R) = \underset{\eta \to 0}{\lim}\, \frac1\eta\int_{\mathcal X} R(q,p)\psi_\eta(q,p)\,\mathrm{d}q\,\mathrm{d}p \end{equation} \]

• Derivative of \(\eta \mapsto \psi_\eta(q,p) \,\mathrm{d}q\, \mathrm{d}p\) at \(0\), tested against \(R\).

• For materials applications, \(F\) and \(R\) have very specific forms informed by physics.

• Running example: mobility. \(F\) constant, response \(R(q,p)=F\cdot M^{-1}p\), particle flux in direction \(F\).

Discretization in time

• Nonequilibrium generator: \[ \mathcal{L}_{\eta} = \mathcal{L}^{\mathrm{A}} + \mathcal{L}^{\mathrm{B}_\eta} +\gamma\mathcal{L}^{\mathrm{O}},\\ \begin{equation} \left\{ \begin{aligned} \mathcal{L}^{\mathrm{A}}&=M^{-1}p\cdot \nabla_q,\\ \mathcal{L}^{\mathrm{B}_\eta}&=\left(-\nabla V+\eta F\right)\cdot \nabla_p,\\ \mathcal{L}^{\mathrm{O}}&=-M^{-1}p\cdot \nabla_p +\frac{1}{\beta}\Delta_p,\\ \end{aligned}\right. \end{equation} \]
• As in the equilibrium case, use splitting schemes.
• Individually, these are the generators of analytically integrable dynamics.
• Example: nonequilibrium BAOAB scheme \[ \mathrm{e}^{\frac{\Delta t}2 \mathcal{L}^{\mathrm{B}_\eta}}\mathrm{e}^{\frac{\Delta t}2 \mathcal{L}^{\mathrm{A}}}\mathrm{e}^{\Delta t\mathcal{L}^{\mathrm{O}}}\mathrm{e}^{\frac{\Delta t}2\mathcal{L}^{\mathrm{A}}}\mathrm{e}^{\frac{\Delta t}2 \mathcal{L}^{\mathrm{B}_\eta}}. \]
• Simply “ignore” \(F\) is non-gradient and add it to the force: possible to abuse equilibrium integrators based on this splitting of the generator.

NEMD method

• Finite-difference estimator of the linear response, based on an ergodic average: \[ \widehat{\rho}(R,F)_{\eta,T} = \frac1{\eta T}\int_0^T R(q_t^\eta,p_t^\eta)\,\mathrm{d}t. \]
• In practice, replace time integral by average over \(N_{\rm sim} = \lceil T/\Delta t\rceil\) steps of a numerical trajectory.
• Variance scales asymptotically as \(C(R)T^{-1}\eta^{-2} \implies T > C(R)\varepsilon^{-2}\eta^{-2}\) for CIs of order \(\varepsilon\) .

Green–Kubo method

• Equilibrium formula: \[\rho(R,F) = \int_0^\infty \mathbb{E}_{\psi_0}\left[R\left(q_t,p_t\right)S\left(q_0,p_0\right)\right]\,\mathrm{d}t,\quad S = (F\cdot \nabla_p)^* {\bf 1}_{\mathcal X}.\]
• Conjugate response: \(S(q,p) = \beta F(q)\cdot M^{-1}p\) (integration by parts).
• Proof by perturbative expansion \((\mathcal{L}+\eta F\cdot \nabla_p)^\dagger\left[\psi_0\left(1+\eta r_1 + \eta^2 r_2 + \dots\right)\right] = 0\).
• Estimate of correlation \(\mathbb{E}_{\psi_0}\left[R\left(q_{n\Delta t},p_{n\Delta t}\right)S\left(q_0,p_0\right)\right]\) from a single trajectory: \[\frac{1}{N_{\rm sim}-n+1}\sum_{k=0}^{N_{\rm sim}}R(q^{n+k},p^{n+k})S(q^k,p^k), \] for \(0\leqslant n\leqslant N_{\rm corr} \ll N_{\rm sim}\) and equilibrium \((q^0,p^0)\).
• Truncate integral in time and compute \(\rho(R,F)\) with numerical quadrature (trapezoid rule).
  
• Extrapolate the correlations to \(t\to\infty\) from short-time correlations by using cepstral analysis, c.f. Bertossa, Grasselli, Ercole & Baroni, 2019, and integrate analytically.

Green–Kubo method: mobility

• Mobility: Forcing \(F(q) = F\), flux \(R(q,p)= F\cdot M^{-1}p\), conjugate response \(S(q,p) = \beta F\cdot M^{-1}p\).
• \[ \rho(F,R) = \beta\int_0^\infty \mathbb{E}_{\psi_0}\left[(F\cdot M^{-1}p_t)(F\cdot M^{-1}p_0)\right]=\beta F^\intercal M^{-1} \left(\int_0^{\infty}\mathbb{E}_{\psi_0}\left[p_t\otimes p_0\right]\,\mathrm{d}t\right)M^{-1}F. \]
• Estimated \(\rho(R,F)\approx 0.1218\).

Molly: setting up nonequilibrium system

• Setup basic Lennard–Jones fluid system

using Molly

n_atoms = 100
boundary = CubicBoundary(2.0u"nm")
temp = 298.0u"K"
atom_mass = 10.0u"g/mol"

atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]
coords = place_atoms(n_atoms, boundary; min_dist=0.3u"nm")
velocities = [random_velocity(atom_mass, temp) for i in 1:n_atoms]
pairwise_inters = (LennardJones(cutoff=ShiftedForceCutoff(0.8u"nm")),);

Molly: setting up nonequilibrium system

To add a nonequilibrium forcing, define a custom interaction type (say MyInterType), and provide a method for Molly.forces(inter::MyInterType,s::System,neighbors) . . . -

struct MyConstantForcing{F}
  forcing::F
end

MyConstantForcing(η::Number,F) = MyConstantForcing(η * F)

Molly.forces(inter::MyConstantForcing,s::System,neighbors=nothing;kwargs...) = inter.forcing
Molly.potential_energy(inter::MyConstantForcing,s::System,neighbors=nothing;kwargs...) = zero(Unitful.numtype(eltype(sys.velocities[1])))*sys.energy_units

η = 10.0u"ps^-1" # very large forcing

const F = zero(velocities) * u"g * mol^-1 "
F[1] = SVector(1.0,0.0,0.0) * u"nm * g * mol^-1 * ps^-1"

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    velocities=velocities,
    pairwise_inters=pairwise_inters,
    general_inters =(MyConstantForcing(η,F),),
    loggers=(temp=TemperatureLogger(10),),
)

System with 100 atoms, boundary CubicBoundary{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}(Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}[2.0 nm, 2.0 nm, 2.0 nm])

Molly: Langevin integrators

sim1 = LangevinSplitting(dt = 2.0u"fs", temperature = temp, friction = 10.0u"g * mol^-1 * ps^-1", splitting = "BAOAB", remove_CM_motion = false)
sim2 = Langevin(dt = 2.0u"fs",temperature = temp, friction = 10.0u"ps^-1",remove_CM_motion = false)

simulate!(sys,sim1,1500)
plot(values(sys.loggers.temp),label="LangevinSplitting (BAOAB)")
empty!(sys.loggers.temp.history) # clear logger
simulate!(sys,sim2,1500);
plot!(values(sys.loggers.temp),label="Langevin")

• Langevin: BAOA splitting with \(\gamma \propto M\).
• LangevinSplitting: General A,B,O splittings, scalar \(\gamma\).

Molly: tracking the flux over trajectories

• Create a function my_response_observable(sys::System,neighbors;n_threads) that returns the microscopic flux.

flux_observable(sys,args...;kwargs...) = F[1][1]*sys.velocities[1][1] # velocity in direction of forcing
sys = System(sys; loggers = (flux = GeneralObservableLogger(flux_observable,typeof(flux_observable(sys)),10),))
simulate!(sys,sim1,10_000)
println(Unitful.dimension(first(values(sys.loggers.flux))))
plot(values(sys.loggers.flux),ylabel="",label=L"R")
hline!([mean(values(sys.loggers.flux)),zero(first(values(sys.loggers.flux)))],color=[:red,:blue],label="")

𝐋^2 𝐌 𝐍^-1 𝐓^-2

Molly: computing time correlations

const β = temp*sys.k
conj_flux_observable(sys,args...;kwargs...) = β*F[1][1]*sys.velocities[1][1]

sys = System(sys; general_inters=(),loggers=(corr = TimeCorrelationLogger(flux_observable,conj_flux_observable,typeof(flux_observable(sys)),typeof(conj_flux_observable(sys)),1,1000),))
simulate!(sys,sim1,10_000;run_loggers=false) # thermalize
simulate!(sys,sim1,10_000)
println(Unitful.dimension(first(values(sys.loggers.corr;normalize=false))^2*sim1.dt))
plot(values(sys.loggers.corr;normalize=false),xlabel="n",label=L"\widehat{C}(n\Delta t)")

𝐋^12 𝐌^6 𝐍^-6 𝐓^-11

• Define a function my_conjugate_response_observable(sys::System,neighbors;n_threads)
• Use the TimeCorrelationLogger logger to estimate the correlation function over a single trajectory.
• Dot-product time correlations for better estimates (also with TimeCorrelationLogger)

Experiments with Molly, past and future:

• Transient substraction technique (momentum kick relaxation rate) + couplings.
• Stochastic Norton method (fix nonequilibrium flux and measure fluctuating forcing).
• DPD transport coefficients (standard methods + Norton).
• Sticky coupling.
• Importance sampling for Green–Kubo.
• Tangent dynamics with ForwardDiff?