A Dual Approach to Computing Transport Coefficients

Computing transport coefficients: a persistent challenge

Proportionality constants in constitutive relations.
$Fourier’s law: J \propto \Delta T$ Example: Fourier’s law: $J = \kappa (T_+-T_-)$, $\kappa=$ thermal conductivity.
Diffusion coefficients (Fick’s law), shear viscosity (Newton’s law).
Crucial to parametrize macroscopic models (e.g. evolution PDEs like Navier–Stokes).
Related to sensitivities of fluxes with respect to forcings in systems close to thermal equilibrium.

Nonequilibrium molecular dynamics (NEMD) modelled by perturbed SDEs \[ \mathrm{d}X^\eta_t = b(X^\eta_t)\,\mathrm{d}t+\sigma(X^\eta_t)\,\mathrm{d}W_t +\eta F(X^\eta_t)\,\mathrm{d}t,\qquad X^\eta_t\in\mathcal X. \] Vector fields $b,F:\mathcal X\to T\mathcal X$, (equilibrium drift and nonequilibrium forcing), $\eta>0$ (forcing strength), $\sigma$ (diffusion matrix), $W$ Brownian motion.
Nonequilibrium steady-state (NESS) is an invariant probability measure \[ \psi_\eta\in\mathcal P(\mathcal X),\qquad \mathcal{L}_\eta^\dagger \psi_\eta=0,\qquad \mathcal{L}_\eta = \mathcal{L}_0 + \eta F\cdot\nabla \]
Given a response observable or flux $R:\mathcal X \to\mathbb R$ (assuming w.l.o.g. $\mathbb E_{\psi_0} R = 0$),
linear response is the derivative of the stationary flux with respect to $\eta$ at equilibrium, i.e. \[ \alpha(F,R) = \underset{\eta \to 0}{\lim}\, \frac1\eta\mathbb E_{\psi_\eta}[R]. \]

Classical estimators of the linear response

Finite-difference ergodic estimators [2] \[ \widehat\alpha_{T,\eta}^{\mathrm{NEMD}}(F,R) = \frac1{T\eta}\int_0^T R(X_t^\eta)\,\mathrm{d}t,\quad \eta>0. \]
Equilibrium time-correlation methods based on the Green–Kubo formula [3], [4] \[ \begin{equation} \widehat\alpha_T^{\mathrm{GK}}(F,R) = \int_0^{T} R(X_t^0)S(X_0^0)\,\mathrm{d}t,\qquad X_0^0 \sim \psi_0 \end{equation} \] where $S = (F\cdot\nabla)^* 1$ is the so-called conjugate response.
Transient methods (relaxation to/from equilibrium from/to nonequilibrium [5],[6]), likelihood-ratio based methods [7], tangent dynamics [8], Einstein’s relation [9].

Generally “bad news” error estimates when $\eta\to 0$, $T\to +\infty$ (reviewed in [10],[11] for typical dynamics)
Typical bias \[ \mathbb{E}\left[\left|\widehat\alpha_{T,\eta}^{\mathrm{NEMD}}(F,R)-\alpha(F,R)\right|\right]=\mathcal{O}\left(\eta\right). \] \[ \mathbb{E}\left[\left|\widehat\alpha_T^{\mathrm{GK}}(F,R)-\alpha(F,R)\right|\right]=\mathcal{O}\left(\mathrm{e}^{-c T}\right), \]
Typical variance \[ \mathrm{Var}\left(\widehat\alpha_{T,\eta}^{\mathrm{NEMD}}(F,R)\right)=\mathcal{O}\left(\frac{1}{\eta^2 T}\right), \] \[ \mathrm{Var}\left(\widehat\alpha_T^{\mathrm{GK}}(F,R)\right)=\mathcal{O}\left(T\right). \]
Need for variance reduction strategies.

Dual approach: fixing the flux

Does the forcing cause the flux ? Macroscopically: co-occurence of the forcing and the flux.
Microscopically, dynamics with fixed forcing / fluctuating flux versus dynamics with a fixed flux / fluctuating forcing.
Constant-flux dynamics \[ \mathrm{d}Y^r_t = b(Y^r_t)\,\mathrm{d}t+\sigma(Y^r_t)\,\mathrm{d}W_t + F(Y^r_t)\mathrm{d}\Lambda^r_t,\quad \mathrm{d}\Lambda^r_t = \lambda(Y_t^r)\,\mathrm{d}t + \widetilde{\lambda}(Y_t^r)\,\mathrm{d}W_t. \]
The forcing stochatic process/Lagrange multiplier $\Lambda^r$ is chosen to enforce the constant flux condition: $R(Y_t^{r}) = r \in {\mathbb{{R}}}$ for all $t>0$.
When $\sigma = 0$, this approach is the “Norton” method of Evans & Morriss (for example [12],[13],[14],[15])

Dual approach: fixing the flux

Derivation of the dynamics

Make the ansatz that $\Lambda^r$ satisfies $\mathrm{d}\Lambda^r_t = \lambda(Y_t^r)\mathrm{d}t + \widetilde{\lambda}(Y_t^r)\mathrm{d}W_t$.
Apply Itô’s formula to the constant flux condition $r=R(Y_t^r)$.
\[ 0=\begin{aligned} &\nabla R(Y_t^r)^\top\left(b(Y_t^r) + F(Y_t^r)\lambda(Y_t^r) + \frac12 \left[\nabla^2 R:(\sigma + F\widetilde\lambda)(\sigma + F\widetilde\lambda)^\top\right](Y_t^r)\right)\,\mathrm{d}t\\ &+\nabla R(Y_t^r)^\top\left(\sigma(Y_t^r)+F(Y_t^r)\,\widetilde{\lambda}(Y_t^r)\right)\,\mathrm{d}W_t \end{aligned}. \]
Identify local martingale and finite variation components, this gives \[ \widetilde{\lambda} = -\frac{\nabla R^\top\sigma}{\nabla R^\top F},\qquad \lambda = -\frac{R^\top}{\nabla R^\top F}\left(b+\frac12 \nabla^2R:\left(\mathrm{Id}-\frac{F\nabla R^\top}{F^\top \nabla R}\right)\sigma\sigma^\top\left(\mathrm{Id}-\frac{F\nabla R^\top}{F^\top \nabla R}\right)^\top\right) \]
In terms of projector-valued functions $P_{F,\nabla R}(y)=\frac{F(y)\nabla R(y)^\top}{F(y)^\top\nabla R(y)},\,\overline{P}_{F,\nabla R}(y) = \mathrm{Id}-P_{F,\nabla R}(y)$, \[ \mathrm{d}Y_t^r = \overline{P}_{F,\nabla R}(Y_t^r)\left(b(Y_t^r)\,\mathrm{d}t + \sigma(Y_t^r)\,\mathrm{d}W_t\right) - \frac12\frac{\nabla^2 R:\overline{P}_{F,\nabla R}\sigma\sigma^\top \overline{P}_{F,\nabla R}}{\nabla R^\top F}(Y_t^r)\,\mathrm{d}t \]
This derivation assumes that $\nabla R^\top F$ does not vanish along trajectories of the dynamics: controllability condition.

Oblique projection

Linear response

Denote by $\psi^r \in \mathcal P(R^{-1}\{r\})$ the constant-flux steady-state (or “Norton ensemble”).
Define \[ \alpha^*(F,R) = \underset{r\to 0}{\lim}\, \frac{r}{\mathbb E_{\psi^r}[\lambda]},\qquad \left(\text{recall } \alpha(F,R) = \underset{\eta\to 0}{\lim}\,\frac{\mathbb E_{\psi_\eta}[R]}{\eta}\right). \]
Set $r(\eta):=\mathbb{E}_{\psi_\eta}[R],\,\eta(r):=\mathbb E_{\psi^r}[\lambda]$. If $\eta(r(\eta)) = \eta + o(\eta)$ and $r(\eta(r))=r + o(r)$, then $\alpha^*(F,R)=\alpha(F,R)$.
Reciprocal estimator of the linear response: \[ \begin{equation}\widehat\alpha_{T,r}^{*}(F,R) = \left(\frac1{T r}\int_0^T \lambda(Y_t^r)\,\mathrm{d}t\right)^{-1},\quad r>0. \end{equation} \]
Link with variance reduction: we use $-\frac1{Tr}\int_0^T \widetilde{\lambda}(Y_t^r)\mathrm{d}W_t$ as a control variate.

Case of nonequilibrium underdamped Langevin dynamics

Nonequilibrium Langevin dynamics, where $F$ is a non-gradient vector field: \[ \left\{ \begin{aligned} \mathrm{d}q_t^\eta &= M^{-1}p_t^\eta\,\mathrm{d}t \\ \mathrm{d}p_t^\eta &= -\nabla V(q_t^\eta)\,\mathrm{d}t + \eta F(q_t^\eta)\,\mathrm{d}t -\gamma M^{-1}p_t^\eta\,\mathrm{d}t + \sqrt{\frac{2\gamma}{\beta}}\,\mathrm{d}W_t \end{aligned}. \right. \]
Consider fluxes of the form $R(q,p)=G(q)^\top p$: suitable to compute mobility, shear viscosity, thermal conductivity in 1D chains (à la [16]). The constant-flux dynamics is given by \[ \left\{\begin{aligned}\mathrm{d}\mathscr{q}_t^r &= M^{-1}\mathscr{q}_t^r\,\mathrm{d}t \\ \mathrm{d}\mathscr{p}_t^r &= \overline{P}_{F,G}(\mathscr{q}_t^r)\left[-\nabla V(\mathscr{q}_t^r)\,\mathrm{d}t -\gamma M^{-1}\mathscr{p}_t^r\,\mathrm{d}t + \sqrt{\frac{2\gamma}{\beta}}\,\mathrm{d}W_t\right] - \frac{\left(\nabla G(\mathscr{q}_t^r)\mathscr{p}_t^r\right)^\top M^{-1}\mathscr{p}_t^r}{F(\mathscr{q}_t^r)^\top G(\mathscr{q}_t^r)}\,\mathrm{d}t\end{aligned}\right. \]

Time discretization: flux-preserving schemes

The decomposition of the generator of the (nonequilibrium) underdamped Langevin dynamics \[ \mathcal{L}_\eta = \underset{A}{\underbrace{p^\top M^{-1}\nabla_q}} \underset{B_\eta}{\underbrace{-\nabla V(q)^\top \nabla_p +\eta F(q)^\top\nabla_p}}+\gamma\underset{O}{\underbrace{\left(-p^\top M^{-1}\nabla_p + \frac1\beta \Delta_p\right)}} \] gives rise to the (nonequilibrium) “$ABO$” family of splitting schemes, where each of the terms generate analytically integrable Markovian dynamics.
For the constant-flux dynamics, the decomposition of the constant-flux dynamics generator \[ \mathfrak{L} = \underset{\mathfrak A}{\underbrace{p^\top M^{-1}\nabla_q - \frac{(\nabla G(q)p)\top M^{-1}p}{F(q)^\top G(q)}F(q)^\top\nabla_p}}\underset{\mathfrak B}{\underbrace{-(\overline{P}_{F,G}(q)\nabla V(q))^\top\nabla_p}}+\gamma\underset{\mathfrak O}{\underbrace{\left(-(\overline{P}_{F,G}(q)M^{-1}p)^\top \nabla_p + \frac1\beta \overline{P}_{F,G}(q)\overline{P}_{F,G}(q)^\top:\nabla_p^2\right)}} \] leads again to analytically integrable flux-preserving dynamics (which are simply the constant-flux versions of the equilibrium “ABO” dynamics).
Implementation in Molly.jl [17].

Some numerical tests: mobility in a Lennard–Jones fluid

From now on, positions live in a periodic box $(L\mathbb T)^{3N}$. We consider $F(q)=F \in {\mathbb{{R}}}^{3N}$ a constant, and $R_F(q,p)=F^\top M^{-1}p$. In particular, the linear response $\alpha(e_{1,x},R_{e_{1,x}})$ is, by definition, the mobility.
Take the “color-drift” forcing $F = \frac{1}{\sqrt{N}}\sum_{j=1}^N (-1)^j e_{j,x}$, with a potential $V(q)= \frac12\sum_{1\leqslant j< k\leqslant N} v(|q_j-q_k|)$, $v(r) =4(r^{-12}-r^{-6})$.

Linear responses coincide ($N=1000$).

Non-linear responses also coincides. Question of full equivalence of nonequilibrium ensembles (in the regime $N\to +\infty$).

Some numerical tests: mobility in a Lennard–Jones fluid

Take $F = e_{1,x}$, the “single-drift” forcing.
Question of sufficient conditions on $(F,R)$ for equivalence of (non)linear responses.

Cautionary example: take a Langevin particle on $\mathbb T \times {\mathbb{{R}}}$. Set a smooth potential $M=1$, $(F,R)=(1,p)$, and $r>0$. The dynamics writes \[ \left\{\begin{aligned} \mathrm{d}\mathscr{q}_t^r &= p\\ \mathrm{d}\mathscr{p}_t^r & = 0 \end{aligned}\right.,\qquad \mathscr{p}_0^r = r. \]
Therefore, $\psi^r(\mathrm{d}q \mathrm{d}p) = \mathrm{d}q|_{\mathbb T}\otimes \delta_r(\mathrm{d}p)$, and $\lambda(q,p) = V'(q)+\gamma p$.
Integrating, $\mathbb E_{\psi^r}[\lambda]=r\gamma$, therefore $\alpha^*(F,R)=\gamma^{-1}$. 0n the other hand the mobility $\alpha(F,R)$ clearly depends on $V$.

Shear viscosity of a Lennard–Jones fluid

The sinusoidal transverse force method [18], [19].

Fixing $f:L\mathbb T\to {\mathbb{{R}}}$ a periodic profile, define \[ F(q) = \sum_{j=1}^N f(q_{j,y})\mathrm{e}_{j,x},\, R(q,p) = \frac{1}{N}\sum_{j=1}^N \left(M^{-1}p\right)_{j,x}\mathrm{e}^{\frac{2\mathrm{i}\pi q_{j,y}}{L}}. \]
The shear viscosity can be defined, using an analogy with Newton’s law of viscosity, as \[ \mu = \rho\left(\frac{F_1}{U_1}-\gamma\right)\left(\frac{L}{2\pi}\right)^2,\, U_1 = \alpha(F,R),\, F_1 = \frac1L\int_0^L f(y)\mathrm{e}^{\frac{2i\pi y}{L}}\,\mathrm{d}y, \] where $\rho = N/L^3$ is the particle density.

Quantification of the efficiency

Assume that equivalence of the constant-forcing and constant-flux linear responses holds.
Which ensemble to choose to measure the linear response ? The one with lower asymptotic variance for linear response estimators.

Asymptotic variance as a function of system size, at equivalent nonequilibrium state points (shear viscosity of a Lennard-Jones fluid).

Anomalous scaling of the stationary variance: $\mathbb E_{\pi^r_{N}}[(\lambda - \langle \lambda_{\pi_N^r} \rangle)^2]\propto N^{-5/3}$

Some perspectives

Theoretical questions. Identify sufficient conditions on $(F,R)$ for the following.
- Ergodicity of the constant-flux dynamics.
- Equivalence of linear responses in the limit $N\to +\infty$.
- Equivalence of nonequilibrium ensembles.
- Asymptotics as $N\to\infty$ for typical fluctuations and correlation time of the forcing process.
Some ongoing work by Shiva Darshan & Gabriel Stoltz in a mean-field setting (next talk).
Similar results have been obtained recently for constant-flux dissipative particle dynamics [20].

Some bibliography

1.

Blassel, N., Stoltz, G.: Fixing the flux: A dual approach to computing transport coefficients. Journal of Statistical Physics. 191, (2024)

2.

Ciccotti, G., Jacucci, G.: Direct computation of dynamical response by molecular dynamics: The mobility of a charged Lennard–Jones particle. Physical Review Letters. 35, 789–792 (1975)

3.

Green, M.S.: Markoff random processes and the statistical mechanics of time‐dependent phenomena. II. Irreversible processes in fluids. The Journal of Chemical Physics. 22, 1281–1295 (1954)

4.

Kubo, R.: Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems. Journal of the Physical Society of Japan. 12, 570–586 (1957)

5.

Arya, G., Maginn, E.J.J., Chang, H.-C.: Efficient viscosity estimation from molecular dynamics simulation via momentum impulse relaxation. The Journal of Chemical Physics. 113, 2079–2087 (2000)

6.

Morriss, G.P., Evans, D.J.: Application of transient correlation functions to shear flow far from equilibrium. Physical Review A. 35, 792–797 (1987)

7.

Plecháč, P., Stoltz, G., Wang, T.: Convergence of the likelihood ratio method for linear response of non-equilibrium stationary states. ESAIM: Mathematical Modelling and Numerical Analysis. 55, S593–S623 (2021)

8.

Assaraf, R., Jourdain, B., Lelievre, T., Roux, R.: Computation of sensitivities for the invariant measure of a parameter dependent diffusion. Stochastics and Partial Differential Equations: Analysis and Computations. 6, 125–183 (2018)

9.

Rodenhausen, H.: Einstein’s relation between diffusion constant and mobility for a diffusion model. Journal of Statistical Physics. 55, 1065–1088 (1989)

10.

Stoltz, G.: Error estimates and variance reduction for nonequilibrium stochastic dynamics. In: Hinrichs, A., Kritzer, P., and Pillichshammer, F. (eds.) Monte carlo and quasi-monte carlo methods (MCQMC 2022). pp. 163–187 (2024)

11.

Blassel, N., Carillo, L., Darshan, S., Gastaldello, R., Iacobucci, A., Marini, E., Santet, R., Shang, X., Stoltz, G., Vaes, U.: Mathematical analysis and numerical methods for the computation of transport coefficients in molecular dynamics. In preparation. (2025)

12.

Evans, D.J., Hoover, W.G., Failor, B.H., Moran, B., Ladd, A.J.C.: Nonequilibrium molecular dynamics via Gauss’s principle of least constraint. Physical Review A. 28, (1983)

13.

Evans, D.J., Ely, J.F.: Viscous flow in the stress ensemble. Molecular Physics. 59, 1043–1048 (1986)

14.

Evans, D.J.: The equivalence of Norton and Thévenin ensembles. Molecular Physics. 80, (1993)

15.

Evans, D.J., Morriss, G.P.: Statistical mechanics of nonequilbrium liquids. ANU Press (2007)

16.

Lefevere, R.: On the local space–time structure of non-equilibrium steady states. Journal of Statistical Mechanics: Theory and Experiment. 2007, P01004 (2007)

17.

Blassel, N.: Github repository, (2023)

18.

Gosling, E.M., McDonald, I.R., Singer, K.: On the calculation by molecular dynamics of the shear viscosity of a simple fluid. Molecular Physics. 26, 1475–1484 (1973)

19.

Joubaud, R., Stoltz, G.: Nonequilibrium shear viscosity computations with Langevin dynamics. Multiscale Modeling & Simulation. 10, 191–216 (2012)

20.

Wu, X., Shang, X.: Stochastic Norton dynamics: An alternative approach for the computation of transport coefficients in Dissipative Particle Dynamics. Journal of Computational Physics. 541, 114316:1–19 (2025)