Computing transport coefficients using Molly.jl

\[ \mathrm{d}X_t^\eta = \color{blue}{b(X_t^\eta)\,\mathrm{d}t + \sigma(X_t^\eta)\,\mathrm{d}W_t} + \color{red}{\eta\,F(X_t^\eta)\,\mathrm{d}t} \]

Equilibrium MD (Langevin, dissipative particle dynamics…) + nonequilibrium forcing, of magnitude \(\eta\)

The state \(X^\eta_t\in\mathcal{X}\) (positions and momenta) evolves stochastically

The NESS is an invariant probability measure \(\mu_\eta \in \mathcal{P}(\mathcal{X})\) under \(X^\eta\)

Equilibrium is characterized by an explicit steady-state (Gibbs distribution) \[ \color{blue}{\mu_0(\mathrm{d}x) = \frac1{Z(\beta)}\mathrm{e}^{-\beta E}\,\mathrm{d}x} \]

Given a flux observable \(R:\mathcal{X}\to\mathbb{R}\) with \(\langle R\rangle_{\eta=0} = 0\), the macroscopic flux is \[ \langle R\rangle_\eta = \int_\mathcal{X}R(x)\,\mu_\eta(\mathrm{d}x) \]

The transport coefficient is the sensitivity of the flux to the forcing \[ \alpha_{F,R} = \lim_{\eta\to 0}\,\frac1\eta \langle R\rangle_\eta = \int_{\mathcal{X}}R(x)\,\partial_\eta\left.\mu_\eta(\mathrm{d}x)\right|_{\eta=0} \]

Challenge: derivative of a high-dimensional integral

Fix \(\eta>0\), run the nonequilibrium dynamics and make finite difference + ergodic approximations \[ \widehat\alpha_{T,\eta}^{\,\mathrm{NEMD}} = \frac{1}{\eta}\,\frac1T\int_0^T R(X_t^\eta)\,\mathrm{d}t \]

Variance \(\mathcal O(T^{-1}\eta^{-2})\): \(\eta\gg 1 \rightarrow\) low variance/high bias or \(\eta\ll 1 \rightarrow\) high variance/low bias

Linear response theory gives the Green-Kubo formula \[ \alpha_{F,R} = \int_0^{\infty} \mathbb E_{X_0\sim \mu_0}\left[R(X_t^{0})\,S_F(X_0^{0})\right]\,\mathrm{d}t \]

Integrand is an average over equilibrium paths of the dynamics \(X^0\)

\(S_F:\mathcal{X}\to\mathbb{R}\) is the conjugate flux observable, explicit in terms of \(F\)

Simple estimator \[ \widehat\alpha_T^{\,\mathrm{GK}} = \frac1{K}\sum_{k=1}^K \int_0^T R(X_t^{0,(k)})S_{F,R}(X_0^{0,(k)})\,\mathrm{d}t \]

Unbiased when \(T\to\infty\), but time-truncation \(\rightarrow\) bias/variance tradeoff

Only requires defining a custom interaction (a type, and a method for AtomsCalculators.forces)

struct ShearForcing{F,T} # interaction type
    force_profile::F
    η::T
end

function AtomsCalculators.forces(s::System,
        inter::ShearForcing, args...; kwargs...) # force computation
    f = ustrip_vec.(similar(s.coords)) * s.force_units
    @inbounds for i in eachindex(f)
        fx = inter.η * inter.force_profile(s.coords[i][2]) # force in x according to y-position
        f[i] = SVector(fx, zero(fx), zero(fx))
    end
    return f
end

profile(y) = sin(2π * y / L) * force_unit # sinusoidal force profile
forcing = ShearForcing(profile, 5.0)

sys_nemd = System(sys; general_inters=(forcing,),
    loggers=(R=GeneralObservableLogger(fourier_flux, TR, 1),)) # log the flux

Only requires adding a TimeCorrelationLogger

function shear_conjugate_flux(s::System, args...; kwargs...) # define the conjugate flux S_F
    return β * sum(s.velocities[i][1] * sin(2π * s.coords[i][2] / L)
                   for i in 1:N) * s.force_units
end
const TS = typeof(shear_conjugate_flux(sys))
n_corr = 7500 # correlation window

gk = TimeCorrelationLogger(shear_conjugate_flux, fourier_flux,
                           TS, TR, 1, n_corr)

sys_gk = System(sys; loggers=(gk=gk,)) # log the autocorrelation

Recall the NEMD dynamics \[ \mathrm{d}X_t^\eta =\color{blue}{ b(X_t^\eta)\,\mathrm{d}t + \sigma(X_t^\eta)\,\mathrm{d}W_t} + \color{red}{\eta F(X_t^\eta)\,\mathrm{d}t} \]

Replace the constant forcing magnitude \(\eta\) by a stochastic forcing process \[ \mathrm{d}Y_t^r =\color{blue}{ b(Y_t^r)\,\mathrm{d}t + \sigma(Y_t^r)\,\mathrm{d}W_t} +\,\color{green}{ F(Y_t^r)\mathrm{d}\Lambda_t^r} \]

chosen so that the flux is fixed exactly \[ R(Y_t^r) = r,\qquad \forall\, t\geqslant 0 \]

This determines \(\Lambda^r\) as \[ \Lambda_t^r = \int_0^t \lambda(Y_s^r)\,\mathrm{d}s + \int_0^t \widetilde\lambda(Y_s^r)\,\mathrm{d}W_s \] for explicit observables \(\lambda:\mathcal{X}\to\mathbb{R}\) and \(\widetilde\lambda:\mathcal{X}\to \mathbb{R}^{1\times d}\)

\[ \alpha_{F,R} = \lim_{\eta\to 0}\frac{\langle R\rangle_\eta}{\eta} \qquad \alpha^*_{F,R} = \lim_{r\to 0}\frac{r}{\langle \lambda\rangle_r} \]

Measure a resistance by averaging \(\lambda\) along the constrained trajectory

Dual pair of nonequilibrium ensembles wrt flux/forcing

Fixed exactly vs fixed on average by the steady state (like the microcanonical and canonical duality for energy/temperature)

Averaging \(\lambda\) discards the martingale component \(\widetilde\lambda\) variance reduction with a control variate