Nonequilibrium Molecular Dynamics Simulation

• NEMD simulation is an atomistic method that applies external fields to drive systems out of equilibrium and measure transport properties.
• It employs dynamic constraints like Nosé–Hoover and Langevin thermostats to sustain nonequilibrium steady states and manage local thermalization.
• NEMD enables direct extraction of transport coefficients such as thermal conductivity and shear viscosity, informing nanoscale materials design.

Nonequilibrium Molecular Dynamics (NEMD) simulation is an atomistic computational method for investigating the response of molecular systems that are driven out of equilibrium by imposed external fields or constraints. Unlike equilibrium MD, where time averages reflect the equilibrium distribution, NEMD generates nonequilibrium steady states (NESS) by subjecting model systems—typically fluids or solids composed of interacting particles—to gradients of temperature, velocity, chemical potential, electric field, or other external perturbations, while employing dynamic constraints to enforce stationarity of key macroscopic quantities. This approach provides unique access to transport coefficients, nonlinear far-from-equilibrium phenomena, and microscopic signatures of irreversibility across a wide array of materials and processes.

1. Simulation Protocols and Dynamic Constraints

The core of an NEMD simulation is the construction of an atomistic system—often with thousands to millions of interacting particles placed in a simulation cell with periodic or specialized boundary conditions—to which external fields are applied in order to establish nonequilibrium gradients. These perturbations can be steady (such as shear or a constant electric field) or time-dependent (e.g., an oscillatory field) (Hoover et al., 2010, Wang et al., 2013).

Central to achieving and sustaining a stationary NESS is the implementation of thermostats or barostats, which act as dynamic constraints. Typical schemes include:

• Nosé-Hoover (chain) thermostats: Deterministically couple a friction coefficient to velocity, regulating the mean kinetic energy in a region. However, their global nature may lead to nonphysical temperature profiles—especially in nanoscale or asymmetric systems—due to spatial overshoots and poor local thermalization (Li et al., 2019).
• Langevin thermostats: Stochastically reassign momenta with dissipative and random forces, enabling local and efficient temperature control within the thermostat regions. This approach is particularly robust, rapidly establishing local equilibrium and mitigating artifacts like spurious thermal rectification (Li et al., 2019, Wang et al., 2013).

A generic set of equations of motion with a Gaussian isokinetic or Nosé-Hoover thermostat is:

$$\begin{cases} m\,\frac{d^2\mathbf{r}_i}{dt^2} = \mathbf{F}_i - \zeta_i\,\mathbf{p}_i \ \frac{d\zeta_i}{dt} = \dfrac{\left(\mathbf{p}_i^2 / (m k T_i)\right) - 1}{\tau_i^2} \end{cases}$$

where $\zeta_i$ is a time-dependent friction, $\tau_i$ is a relaxation parameter, and $T_i$ is the target temperature (Hoover et al., 2010).

2. Extraction and Interpretation of Transport Properties

NEMD provides direct access to transport coefficients by forcing stationary fluxes and measuring observable gradients. Important examples include:

• Thermal conductivity ($\kappa$): Imposed via source/sink thermostats; computed as $\kappa = Q/|\nabla T|$, where $Q$ is the steady-state heat flux and $\nabla T$ is the temperature gradient (Xu et al., 2018, Li et al., 2019).
• Shear viscosity ($\eta$): Determined from the linear (Newtonian) response to imposed shear flows, either via moving walls or an off-diagonal velocity gradient imposed through boundary conditions (Hoover et al., 2010).
• Interfacial thermal conductance ($G$): Obtained by direct measurement of the energy flux and the temperature drop at an interface, $G = Q / \Delta T$ (Merabia et al., 2012).

In practice, steady-state profiles are often nonlinear near the energy reservoirs due to the Kapitza resistance or other interface-related effects. Recent work demonstrates that for nanoscale systems, the proper evaluation of conductance and conductivity must use the full temperature difference between the thermostatted regions (not just the linear bulk gradient). Neglecting nonlinear interface regions can lead to systematic overestimation of conductance, especially in the ballistic or ballistic–diffusive crossover regimes (Li et al., 2019).

3. Phase-Space Structure and Microscopic Irreversibility

Although NEMD typically implements time-reversible dynamics, the imposition of thermostats and external driving fields enforces the emergence of a stationary nonequilibrium phase-space distribution with distinct fractal properties. The resulting NESS is characterized by multifractal (“strange attractor”) phase-space distributions, reflecting the microscopic foundations of the macroscopic irreversibility embodied in the Second Law of Thermodynamics (Hoover et al., 2010). In these states, the underlying phase-space volume contracts, despite the time-reversible equations of motion—this behavior is a haLLMark of nonequilibrium statistical mechanics.

Moreover, the emergence of key macroscopic phenomena, such as normal stress differences in nonlinear shear flow or the tensorial nature of kinetic temperature in shockwaves, is directly traceable to the microdynamics captured by NEMD (Hoover et al., 2010).

4. Advanced Methods, Boundary Conditions, and Efficiency

To simulate homogeneous flows (shear, elongational, uniaxial, or three-dimensional deformations), special periodic boundary conditions and efficient force calculation algorithms are mandatory.

• Specialized boundary conditions: Lees–Edwards for shear, and Kraynik–Reinelt or “rotating box” approaches for elongational flows, prevent the simulation box from excessive distortion over long timescales (Dobson et al., 2021). The most recent advances introduce time-periodic remappings (modulo rotation) based on integer automorphisms with complex eigenvalues, improving minimum image spacing and computational stability.
• Cell list algorithms: In deforming simulation cells, dynamic cell list methods (either congruent with box deformation or with fixed geometry and dynamic offsets) enable $O(N)$ scaling for short-range force computations even as the box deforms with the flow (Dobson et al., 2014).

These schemes, often deployed in tandem, support the simulation of very large systems and long time scales necessary for the reliable extraction of nonlinear response and rare-event statistics.

5. Methodological Comparisons and Limitations

NEMD is commonly evaluated alongside equilibrium MD (EMD) and homogeneous NEMD (HNEMD) for the computation of transport coefficients:

Method	Key Strengths	Limitations
NEMD	Direct access to nonlinear, far-from-equilibrium response; intuitive setup for imposed fluxes; enables local analysis and spatial profiles	May suffer from finite-size effects and increased computational cost for slow-relaxing systems; typically needs multiple simulations per parameter set
EMD/HNEMD	All response functions from a single run (Green–Kubo or Einstein–Helfand formalism); more efficient for full-spectrum viscoelasticity or for systems with multiple coupled fluxes	Requires accurate sampling of autocorrelation functions, suffers higher statistical noise at long times, less efficient for single property/perturbation studies (Mangaud et al., 2020, Adeyemi et al., 2022)

For certain nanoscale problems (including ionic, electro-osmotic, or mixed-gradient transport in membranes), NEMD offers a unified simulation framework that can directly probe Onsager off-diagonal coefficients and capture cross-phenomena (e.g., streaming currents, diffusio-osmosis), achieving full symmetry of the dynamic Onsager matrix and facilitating quantitative modeling of multicomponent flows (Monet et al., 2023).

6. Applications, Frontier Developments, and Theoretical Insights

NEMD underpins quantitative atomistic modeling across a wide span of emerging material systems and physical regimes:

• Nanofluidics and filtration: Direct simulation of water and ion transport through nanoporous membranes—where the methodology reveals the dominance of entrance effects and nearly frictionless internal flow in thin carbon-based membranes, crucial to their exceptional selectivity and permeance (Monet et al., 2023).
• Phonon transport and size effects: Spectral-resolved approaches (e.g., TDDDM and FDDDM) decompose the modal contributions to thermal conductivity, quantifying the truncation of long-MFP phonons and the impact of finite-size scaling, and revealing precise relationships between simulated and theoretical conductance (Zhou et al., 2015, Zhou et al., 2015, Merabia et al., 2012, 2207.13405).
• Polymer viscoelasticity: Oscillatory NEMD with optimized pre-averaging and DFT signal analysis facilitates the extraction of storage/loss moduli over broad frequency ranges with significantly reduced simulation time (Adeyemi et al., 2022).
• Large-scale and high-accuracy materials design: The incorporation of machine-learning force fields (e.g., atomic cluster expansion with 4-body terms) now enables DFT-level accuracy for NEMD thermal transport simulations in systems exceeding $10^5$ atoms, including the correct reproduction of phonon dispersion and vibrational spectra for predictive device-scale modeling (Araki et al., 2023).
• Coarse-grained and integrable system connections: Recent NEMD studies on SWCNTs demonstrate that the coarse-grained potential in the NESS regime adopts a Toda lattice form, supporting solitonic dynamics and providing a link between anomalous thermal transport in low-dimensional materials and integrable system theory (Koh et al., 24 Jul 2025).

7. Limitations, Uncertainties, and Emerging Challenges

While NEMD is a powerful and flexible method, significant considerations include:

• Sensitivity to simulation parameters: The accuracy of predictions depends strongly on thermostatting protocols, boundary condition handling, and choice of interatomic potentials. Artifacts may arise—e.g., spurious rectification or overestimated conductance—if global thermostats are used or interface regions are mishandled (Li et al., 2019).
• Finite-size/finite-time effects: Especially acute in systems where the relevant mean free path is comparable to the system dimensions, necessitating careful extrapolation and, in some cases, surrogate modeling with uncertainty quantification via polynomial chaos and sensitivity analysis (Vohra et al., 2018).
• Computational cost: For systems requiring many-body accuracy (e.g., large-scale MLFFs) or suffering from slow-relaxing dynamics (entangled polymers, rare event transitions), the computational burden remains substantial despite advances in linear-scaling algorithms.

Continued methodological improvements in force fields, thermostatting, efficient boundary remappings, and statistical error reduction are critical for the ongoing extension of NEMD into new domains—ranging from biomolecular machines under chemical fuel (chemostatted NEMD), to strongly correlated quantum systems, to macroscale multiscale modeling where atomistic data informs continuum and field-theoretic descriptions.

---

NEMD stands as an indispensable tool for probing the origins and the full spectrum of nonequilibrium phenomena at the molecular scale, providing microscopic insight into macroscopic irreversibility, yielding quantitative transport coefficients, and facilitating the rational engineering of nanoscale materials and devices in both fundamental and applied contexts (Hoover et al., 2010, Li et al., 2019, Monet et al., 2023, Araki et al., 2023, Koh et al., 24 Jul 2025).