Reverse Non-Equilibrium Molecular Dynamics

• RNEMD is a simulation method that imposes specific fluxes (momentum, energy, or particles) to generate steady-state gradients, enabling the extraction of transport coefficients.
• Various implementations such as Müller–Plathe, VSS-RNEMD, and SPF-RNEMD address heat, momentum, and mass transport across systems like polymers, interfaces, and mixtures.
• The technique’s strength lies in its precise flux control and profile-based analysis, offering a robust alternative to traditional equilibrium and direct nonequilibrium methods.

Reverse non-equilibrium molecular dynamics (RNEMD) is a family of molecular dynamics methods in which a transport flux is imposed directly and the conjugate gradient is measured in the resulting steady state. In this “reverse” formulation, one does not prescribe a velocity gradient, temperature difference, or concentration difference and then observe the flux; instead, one injects momentum, energy, or particles through controlled microscopic operations and extracts viscosity, thermal conductivity, diffusivity, interfacial conductance, or friction from flux–gradient relations. Published implementations include the Müller–Plathe swap scheme for heat and momentum transport, velocity shearing and scaling RNEMD (VSS-RNEMD), and the scaled particle flux RNEMD (SPF-RNEMD) algorithm for mass transport in mixtures (Oishi et al., 2024, Felix et al., 2024, Drisko et al., 2024).

1. Conceptual framework

RNEMD occupies a distinct position relative to equilibrium MD and direct nonequilibrium MD. In equilibrium MD, transport coefficients are inferred from spontaneous fluctuations, as in Green–Kubo or Einstein relations. In direct nonequilibrium MD, one imposes a gradient or boundary driving and measures the resulting flux. RNEMD reverses that logic: a flux is prescribed, the system evolves under otherwise standard dynamics, and the emergent stationary gradient is analyzed with the same constitutive laws used in continuum transport theory (Felix et al., 2024, Drisko et al., 2024).

For shear flow, the relevant constitutive relation is

$j_z(p_x) = -\eta \left(\frac{d v_x}{d z}\right),$

or, equivalently in the polymer formulation,

$\sigma = -\eta\,\dot\gamma.$

For heat transport, Fourier’s law is used,

$$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$

For binary diffusion, the corresponding Fick form is

$\mathbf{J}_1 = -D \nabla c_1.$

In each case, RNEMD replaces fluctuation analysis or externally imposed gradients with an imposed flux and a measured response field (Oishi et al., 2024, Felix et al., 2024, Drisko et al., 2024).

A central feature of classical RNEMD is that no external body forces or background flow fields are required. In the momentum-swap formulation used for polymer shear, the equations of motion are otherwise standard Newtonian dynamics, and total momentum and kinetic energy are conserved by construction for equal masses (Oishi et al., 2024). This suggests why RNEMD is often attractive when one wishes to avoid continuously thermostatted shear algorithms or externally maintained hot and cold reservoirs, although later variants also superimpose compensating heat fluxes or local velocity scaling when interfacial or coupled transport problems demand it (Harless et al., 5 Aug 2025, Shavalier et al., 2023).

2. Algorithmic realizations

The canonical Müller–Plathe heat-flux scheme partitions the simulation cell into slabs, designates one slab as “cold” and another as “hot,” and periodically exchanges kinetic energy between fast particles in the cold slab and slow particles in the hot slab. In the Sun-graphyne study, the system is divided into $n$ slabs along the heat-flow direction, the first slab is chosen as the cold region and the middle slab as the hot region, and swaps are performed every 500 MD steps. The cumulative transferred kinetic energy defines the imposed heat flux, while a symmetric temperature profile develops because periodicity creates two heat-flow directions from the central hot slab toward the boundaries (Felix et al., 2024).

For shear viscosity, the Müller–Plathe idea is applied to momentum exchange rather than energy exchange. In the Kremer–Grest polymer melt study, the box is divided into slabs along the gradient direction $z$; one exchange slab lies at the top of the box and the other at the center. At each exchange event, the particle with the largest positive $v_x$ in the top slab and the particle with the smallest $v_x$ in the center slab are identified, and their $x$-velocity components are swapped. The mean transferred momentum defines a shear stress,

$\sigma = \frac{\langle p_x\rangle}{2L_xL_y},$

and the steady-state velocity profile yields the shear rate $\sigma = -\eta\,\dot\gamma.$0 (Oishi et al., 2024).

VSS-RNEMD replaces discrete pair swaps with continuous, small, deterministic modifications of velocities in two spatial regions. In the ice–water interfacial study, VSS-RNEMD imposes a momentum flux $\sigma = -\eta\,\dot\gamma.$1 and, for interfacial systems, a heat flux $\sigma = -\eta\,\dot\gamma.$2 in the opposite direction to frictional heating. The method modifies particle velocities in two slabs so that prescribed amounts of momentum and energy are exchanged every 2 fs, while conserving total linear momentum in $\sigma = -\eta\,\dot\gamma.$3 and keeping total energy consistent with the imposed fluxes (Harless et al., 5 Aug 2025). In the gold–PEG study, the same broad VSS idea is used to generate a thermal flux in planar slabs or concentric spherical shells, with slab or shell center-of-mass velocities and scaling factors chosen so that total energy and linear momentum are conserved at each step (Shavalier et al., 2023).

SPF-RNEMD extends the RNEMD framework to particle transport in mixtures. Drisko and Gezelter construct a migrating molecule that exists simultaneously in a source region and a sink region, with potential energies and forces linearly combined as

$\sigma = -\eta\,\dot\gamma.$4

$\sigma = -\eta\,\dot\gamma.$5

using $\sigma = -\eta\,\dot\gamma.$6. A continuously increasing scaling variable $\sigma = -\eta\,\dot\gamma.$7 drives the migration, and the imposed particle flux is

$\sigma = -\eta\,\dot\gamma.$8

Local velocity rescaling in the exchange slabs compensates the potential-energy change associated with migration and can simultaneously impose a heat flux (Drisko et al., 2024).

3. Profile construction and transport-coefficient extraction

RNEMD analysis is profile-based. The simulation cell is partitioned into bins or slabs along the transport direction, and time-averaged fields are accumulated once a stationary state has formed. For shear flow, the local mean streaming velocity is measured slabwise, and the shear rate is obtained from a linear fit away from the exchange slabs, where sharp spikes and local nonlinearity occur. In the polymer melt study, the region $\sigma = -\eta\,\dot\gamma.$9 was used for fitting because it was nearly linear and free of exchange-slab artifacts (Oishi et al., 2024).

Thermal RNEMD similarly relies on a spatial temperature profile. In the Sun-graphyne calculations, slab temperatures are computed from the kinetic energy of atoms in each slab, the nonlinear regions near the hot and cold slabs are discarded, and the remaining central region is fit to obtain $$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$0. The thermal conductivity then follows from

$$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$1

Because periodicity creates two equal and opposite heat-flux directions, a factor of $$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$2 enters the heat-flux definition (Felix et al., 2024).

In interfacial systems, RNEMD is also used to extract jump conditions rather than only bulk gradients. For ice–water interfaces, the friction coefficient is defined through the force–velocity relation

$$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$3

where $$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$4 is the transverse velocity drop across the interface. The ice–water study defines the relevant solid-like and liquid-like velocities using the Gibbs dividing surface and the $$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$5 interfacial width obtained from the tetrahedral order parameter profile. This construction is used because classical slip-length-based formulas become ill-defined for negative slip lengths (Harless et al., 5 Aug 2025).

For interfacial heat transport, the gold–PEG work uses the imposed heat flux together with a temperature drop across the interfacial region. In planar geometries the interfacial thermal conductance is

$$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$6

with $$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$7 measured from the last slab containing only gold to the first slab containing only water. In spherical geometries, a series-resistance approximation over concentric shells is used, so that the interfacial Kapitza resistance is built from shellwise temperature differences weighted by $$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$8 (Shavalier et al., 2023).

For diffusion in mixtures, SPF-RNEMD produces a steady concentration gradient from an imposed particle flux. In the binary-mixture formulation used by Drisko and Gezelter,

$$\kappa = -\frac{\langle J\rangle}{\langle \nabla T\rangle}.$$9

so the Fick diffusion coefficient follows directly from the measured concentration profile. When heat and particle fluxes are imposed simultaneously, the same framework yields coupled transport quantities such as the Soret coefficient and a temperature-dependent diffusion coefficient extracted from local gradients along the simulation cell (Drisko et al., 2024).

4. Representative applications

The breadth of RNEMD is illustrated by recent applications spanning polymer rheology, supercooled interfacial liquids, low-dimensional thermal transport, metal–ligand–solvent heat transfer, and coupled mass–heat transport in mixtures.

System	RNEMD variant	Principal quantities
Unentangled Kremer–Grest polymer melt	Momentum-swap RNEMD	$\mathbf{J}_1 = -D \nabla c_1.$0, $\mathbf{J}_1 = -D \nabla c_1.$1, $\mathbf{J}_1 = -D \nabla c_1.$2, gyration tensor
Ice–water interfaces with solutes	VSS-RNEMD	$\mathbf{J}_1 = -D \nabla c_1.$3, $\mathbf{J}_1 = -D \nabla c_1.$4, interfacial width
Sun-graphyne	Müller–Plathe heat-flux RNEMD	$\mathbf{J}_1 = -D \nabla c_1.$5, $\mathbf{J}_1 = -D \nabla c_1.$6, $\mathbf{J}_1 = -D \nabla c_1.$7
Thiolated PEG on gold interfaces	VSS-RNEMD	Interfacial thermal conductance $\mathbf{J}_1 = -D \nabla c_1.$8
Binary mixtures and nanoporous graphene membranes	SPF-RNEMD	Fick diffusivity, $\mathbf{J}_1 = -D \nabla c_1.$9, permeability

In the polymer-melt study, RNEMD was systematically tested against SLLOD for an unentangled Kremer–Grest chain melt under fast shear. As the shear rate increased, the temperature and density became inhomogeneous, but the average viscosity remained consistent with the SLLOD result under homogeneous temperature and density. The simulations reached $n$0, the zero-shear viscosity was $n$1, and the flow curve was well described by a Carreau-like relation with $n$2 and $n$3, implying $n$4 at high shear. The same study also reported that temperature-density inhomogeneity did not significantly affect polymer conformation, and that the diagonal components of the gyration tensor matched SLLOD across the full $n$5 range (Oishi et al., 2024).

In the ice-active solution study, VSS-RNEMD was used to determine bulk and interfacial viscosities, interfacial widths, and solid–liquid friction coefficients at ice–water interfaces. A momentum flux of $n$6 was imposed, together with a compensating heat flux of $n$7, except for NH$n$8Cl where the heat flux was reduced. The study found a direct correlation between liquid-phase hydrogen-bond jump times and shear viscosity, and reported that only DMSO and sodium formate exhibited increased friction at the ice–water interface (Harless et al., 5 Aug 2025).

In Sun-graphyne, Müller–Plathe RNEMD was used to obtain length-dependent lattice thermal conductivity. The extrapolated intrinsic conductivity was $n$9, with an effective phonon mean free path $z$0. The system displayed ballistic behavior up to about 7 nm, a ballistic–diffusive transition between roughly 7 and 16 nm, and weak length dependence for $z$1. The reduced conductivity relative to graphene was attributed to acetylenic bonds, lower phonon group velocities, and enhanced acoustic–optical scattering (Felix et al., 2024).

In solvated gold interfaces capped with thiolated PEG, VSS-RNEMD was adapted to planar facets and spherical nanoparticles. The interfacial thermal conductance of thiolated PEG capped interfaces was higher than that of pristine gold interfaces. For planar systems, $z$2; for nanospheres, mean $z$3 values were higher, with broad overlap in uncertainties between 10 and 20 Å radii. The largest single temperature drop occurred at the gold–sulfur bond region, indicating that the covalent Au–S linkage was the largest barrier to thermal conduction, while enhanced low-frequency vibrational populations and altered PEG–water coupling helped rationalize the higher conductance of curved systems (Shavalier et al., 2023).

SPF-RNEMD broadened RNEMD to composition-driven transport. In mixtures of identical but distinguishable particles, the method reproduced Fick diffusion coefficients consistent with equilibrium benchmarks except when high fluxes depleted source slabs at very low concentration. In Ar/Kr mixtures, it agreed well with equilibrium and previous RNEMD results; for $z$4, the reported value was $z$5. Under simultaneous heat and particle fluxes, the same framework yielded an Arrhenius activation energy for diffusion of $z$6. In an interfacial membrane application, SPF-RNEMD was used to compute a water diffusive permeability of $z$7 for nanoporous graphene under the stated pore-density assumption (Drisko et al., 2024).

5. Inhomogeneity, artifacts, and methodological limits

RNEMD does not generally produce homogeneous nonequilibrium states. In polymer shear, parabolic temperature and density profiles develop along the gradient direction, together with sharp spikes at the exchange slabs caused by instantaneous velocity swaps. The amplitudes of these inhomogeneities grow with swap rate, and for $z$8 the mapping from swap frequency to shear rate or stress saturates and becomes dependent on the exchange interval $z$9. Nevertheless, the flow curve $v_x$0 collapsed onto a single master curve independent of $v_x$1, $v_x$2, and ensemble, which indicates that the rheology was robust even when the control parameter was not simply related to the shear rate (Oishi et al., 2024).

Thermal RNEMD likewise exhibits localized nonlinearity near the exchange regions. In Sun-graphyne, the temperature profile was explicitly reported to be nonlinear near the hot and cold slabs due to finite-size effects. This is why only the interior linear region was used for the thermal-gradient fit, and why length scaling was essential for extracting an intrinsic conductivity from finite systems (Felix et al., 2024).

Interfacial applications add further complications. In ice–water systems, a small heat flux was required to compensate frictional heating, which introduced a weak temperature gradient across the box and necessitated a modified hyperbolic-tangent fit for the tetrahedral order parameter. Because negative slip lengths occurred, the usual relation $v_x$3 was rejected and replaced by the flux–velocity drop relation $v_x$4 (Harless et al., 5 Aug 2025). In gold–PEG systems, the definition of the interfacial region mattered because the interface was a finite-width mixed zone rather than a sharp plane, and spherical systems required a shell-based series-resistance construction instead of a simple planar area formula (Shavalier et al., 2023).

SPF-RNEMD introduces a different class of limitations. Because the method perturbs positions rather than only velocities, potential-energy changes must be compensated by local velocity scaling. The algorithm therefore restricts scaling to within about $v_x$5 of the original velocities and rejects attempted $v_x$6 increments when the required scaling would be excessive or $v_x$7 or $v_x$8 would become negative. High particle fluxes can deplete source slabs of dilute components, break linear response, and produce failure modes that do not arise in the same form in heat- or momentum-flux RNEMD (Drisko et al., 2024).

A common misconception is that energy and momentum conservation by themselves guarantee absence of artifacts. The published studies do not support that conclusion. They show instead that RNEMD can generate physically useful steady states while still producing exchange-slab spikes, parabolic $v_x$9 and $v_x$0 profiles, finite-size effects, and geometry-dependent analysis issues that must be evaluated case by case (Oishi et al., 2024, Felix et al., 2024).

6. Scope, best practices, and evolving directions

RNEMD has expanded from a technique for thermal conductivity and shear viscosity into a broader framework for transport in heterogeneous and coupled systems. Current applications include homogeneous liquids, polymer melts, two-dimensional crystals, ice–water interfaces, ligand-functionalized nanoparticles, binary mixtures, and nanoporous membranes (Oishi et al., 2024, Shavalier et al., 2023, Drisko et al., 2024). This suggests a methodological trajectory in which RNEMD is less a single algorithm than a design principle: impose a conserved microscopic flux, let the system build the conjugate field, and infer transport coefficients from the stationary profile.

Several best practices recur across the published implementations. One is to measure gradients away from source, sink, hot, cold, or exchange regions, because those regions exhibit spikes or nonlinearity. Another is to verify that the steady-state profiles are linear or otherwise well fit in the intended analysis window. For size-dependent heat transport, multiple lengths should be simulated and fit to a scaling relation such as

$v_x$1

rather than assuming a single finite system is already in the diffusive regime (Felix et al., 2024). For coupled transport or interfacial friction, structural order parameters and carefully defined dividing surfaces are often required to identify the appropriate fitting region or jump condition (Harless et al., 5 Aug 2025).

The recent polymer and mixture studies also clarify the importance of flux selection. At moderate driving, stress and shear rate or concentration gradient may scale simply with the imposed RNEMD rate, making parameter tuning straightforward. At stronger driving, that mapping can become nonlinear even when the extracted constitutive relation remains stable. In such cases, the transport coefficient can still be meaningful, but only if steady-state profiles are well resolved and the measurement region remains representative (Oishi et al., 2024, Drisko et al., 2024).

The present literature also marks clear boundaries. For the polymer system studied under shear, RNEMD remained accurate because the box dimensions were much larger than the molecular size and the induced deviations in the measurement region were small enough that viscosity and conformation were not significantly altered. The same study cautioned that if molecular dimensions approach the scale of temperature or density inhomogeneity, or if interactions are long-ranged, swapping could distort dynamics at domain boundaries; entangled polymers or very long chains may therefore require additional scrutiny (Oishi et al., 2024). In a parallel way, interfacial RNEMD requires explicit attention to local structural gradients, negative slip, or finite-width mixed layers before continuum transport formulas can be applied without ambiguity (Harless et al., 5 Aug 2025, Shavalier et al., 2023).

Taken together, these developments establish RNEMD as a technically versatile route to transport coefficients in regimes where equilibrium fluctuation methods can be noisy and direct gradient-imposing methods can be cumbersome. Its distinctive strength is precise control of flux with direct access to the resulting spatial response; its enduring challenge is that the steady states it creates are inherently structured, and their inhomogeneities are part of the method rather than incidental numerical noise (Felix et al., 2024, Drisko et al., 2024).