Better Than Counting: Density Profiles from Force Sampling

Gauge Invariance of Equilibrium Statistical Mechanics

We identify a recently proposed shifting operation on classical phase space as a gauge transformation for statistical mechanical microstates. The infinitesimal generators of the continuous gauge group form a noncommutative Lie algebra, which induces exact sum rules when thermally averaged. Gauge invariance with respect to finite shifting is demonstrated via Monte Carlo simulation in the transformed phase space which generates identical equilibrium averages. Our results point toward a deeper basis of statistical mechanics than previously known, and they offer avenues for systematic construction of exact identities and of sampling algorithms.

Estimating Free-Energy Surfaces and Their Convergence from Multiple, Independent Static and History-Dependent Biased Molecular-Dynamics Simulations with Mean Force Integration

Addressing the sampling problem is central to obtaining quantitative insight from molecular dynamics simulations. Adaptive biased sampling methods, such as metadynamics, tackle this issue by perturbing the Hamiltonian of a system with a history-dependent bias potential, enhancing the exploration of the ensemble of configurations and estimating the corresponding free energy surface (FES). Nevertheless, efficiently assessing and systematically improving their convergence remains an open problem. Here, building on mean force integration (MFI), we develop and test a metric for estimating the convergence of FESs obtained by combining asynchronous, independent simulations subject to diverse biasing protocols, including static biases, different variants of metadynamics, and various combinations of static and history-dependent biases. The developed metric and the ability to combine independent simulations granted by MFI enable us to devise strategies to systematically improve the quality of FES estimates. We demonstrate our approach by computing FES and assessing the convergence of a range of systems of increasing complexity, including one- and two-dimensional analytical FESs, alanine dipeptide, a Lennard-Jones supersaturated vapor undergoing liquid droplet nucleation, and the model of a colloidal system crystallizing via a two-step mechanism. The methods presented here can be generally applied to biased simulations and are implemented in pyMFI, a publicly accessible, open-source Python library.

Why neural functionals suit statistical mechanics

We describe recent progress in the statistical mechanical description of many-body systems via machine learning combined with concepts from density functional theory and many-body simulations. We argue that the neural functional theory by Sammülleret al(2023Proc. Natl Acad. Sci.120e2312484120) gives a functional representation of direct correlations and of thermodynamics that allows for thorough quality control and consistency checking of the involved methods of artificial intelligence. Addressing a prototypical system we here present a pedagogical application to hard core particle in one spatial dimension, where Percus' exact solution for the free energy functional provides an unambiguous reference. A corresponding standalone numerical tutorial that demonstrates the neural functional concepts together with the underlying fundamentals of Monte Carlo simulations, classical density functional theory, machine learning, and differential programming is available online athttps://github.com/sfalmo/NeuralDFT-Tutorial.

Local measures of fluctuations in inhomogeneous liquids: statistical mechanics and illustrative applications

We show in detail how three one-body fluctuation profiles, namely the local compressibility, the local thermal susceptibility, and the reduced density, can be obtained from a statistical mechanical many-body description of classical particle-based systems. We present several different and equivalent routes to the definition of each fluctuation profile, facilitating their explicit numerical calculation in inhomogeneous equilibrium systems. This underlying framework is used for the derivation of further properties such as hard wall contact theorems and novel types of inhomogeneous one-body Ornstein-Zernike equations. The practical accessibility of all three fluctuation profiles is exemplified by grand canonical Monte Carlo simulations that we present for hard sphere, Gaussian core and Lennard-Jones fluids in confinement.

Noether-Constrained Correlations in Equilibrium Liquids

Liquid structure carries deep imprints of an inherent thermal invariance against a spatial transformation of the underlying classical many-body Hamiltonian. At first order in the transformation field Noether's theorem yields the local force balance. Three distinct two-body correlation functions emerge at second order, namely the standard two-body density, the localized force-force correlation function, and the localized force gradient. An exact Noether sum rule interrelates these correlators. Simulations of Lennard-Jones, Yukawa, soft-sphere dipolar, Stockmayer, Gay-Berne and Weeks-Chandler-Andersen liquids, of monatomic water and of a colloidal gel former demonstrate the fundamental role in the characterization of spatial structure.

Perspective: How to overcome dynamical density functional theory

We argue in favour of developing a comprehensive dynamical theory for rationalizing, predicting, designing, and machine learning nonequilibrium phenomena that occur in soft matter. To give guidance for navigating the theoretical and practical challenges that lie ahead, we discuss and exemplify the limitations of dynamical density functional theory (DDFT). Instead of the implied adiabatic sequence of equilibrium states that this approach provides as a makeshift for the true time evolution, we posit that the pending theoretical tasks lie in developing a systematic understanding of the dynamical functional relationships that govern the genuine nonequilibrium physics. While static density functional theory gives a comprehensive account of the equilibrium properties of many-body systems, we argue that power functional theory is the only present contender to shed similar insights into nonequilibrium dynamics, including the recognition and implementation of exact sum rules that result from the Noether theorem. As a demonstration of the power functional point of view, we consider an idealized steady sedimentation flow of the three-dimensional Lennard-Jones fluid and machine-learn the kinematic map from the mean motion to the internal force field. The trained model is capable of both predicting and designing the steady state dynamics universally for various target density modulations. This demonstrates the significant potential of using such techniques in nonequilibrium many-body physics and overcomes both the conceptual constraints of DDFT as well as the limited availability of its analytical functional approximations.

Reduced-variance orientational distribution functions from torque sampling

We introduce a method to sample the orientational distribution function in computer simulations. The method is based on the exact torque balance equation for classical many-body systems of interacting anisotropic particles in equilibrium. Instead of the traditional counting of events, we reconstruct the orientational distribution function via an orientational integral of the torque acting on the particles. We test the torque sampling method in two- and three-dimensions, using both Langevin dynamics and overdamped Brownian dynamics, and with two interparticle interaction potentials. In all cases the torque sampling method produces profiles of the orientational distribution function with better accuracy than those obtained with the traditional counting method. The accuracy of the torque sampling method is independent of the bin size, and hence it is possible to resolve the orientational distribution function with arbitrarily small angular resolutions.

Comparative study of force-based classical density functional theory

We reexamine results obtained with the recently proposed density functional theory framework based on forces (force-DFT) [S. M. Tschopp et al., Phys. Rev. E 106, 014115 (2022)2470-004510.1103/PhysRevE.106.014115]. We compare inhomogeneous density profiles for hard sphere fluids to results from both standard density functional theory and from computer simulations. Test situations include the equilibrium hard sphere fluid adsorbed against a planar hard wall and the dynamical relaxation of hard spheres in a switched harmonic potential. The comparison to grand canonical Monte Carlo simulation profiles shows that equilibrium force-DFT alone does not improve upon results obtained with the standard Rosenfeld functional. Similar behavior holds for the relaxation dynamics, where we use our event-driven Brownian dynamics data as benchmark. Based on an appropriate linear combination of standard and force-DFT results, we investigate a simple hybrid scheme which rectifies these deficiencies in both the equilibrium and the dynamical case. We explicitly demonstrate that although the hybrid method is based on the original Rosenfeld fundamental measure functional, its performance is comparable to that of the more advanced White Bear theory.

Inhomogeneous steady shear dynamics of a three-body colloidal gel former

We investigate the stationary flow of a colloidal gel under an inhomogeneous external shear force using adaptive Brownian dynamics simulations. The interparticle forces are derived from the Stillinger-Weber potential, where the three-body term is tuned to enable network formation and gelation in equilibrium. When subjected to the shear force field, the system develops remarkable modulations in the one-body density profile. Depending on the shear magnitude, particles accumulate either in quiescent regions or in the vicinity of maximum net flow, and we deduce this strong non-equilibrium response to be characteristic of the gel state. Studying the components of the internal force parallel and perpendicular to the flow direction reveals that the emerging flow and structure of the stationary state are driven by significant viscous and structural superadiabatic forces. Thereby, the magnitude and nature of the observed non-equilibrium phenomena differ from the corresponding behavior of simple fluids. We demonstrate that a simple power functional theory reproduces accurately the viscous force profile, giving a rationale of the complex dynamical behavior of the system.

Kirkwood-Buff integrals: From fluctuations in finite volumes to the thermodynamic limit

The Kirkwood-Buff theory is a cornerstone of the statistical mechanics of liquids and solutions. It relates volume integrals over the radial distribution function, so-called Kirkwood-Buff integrals (KBIs), to particle number fluctuations and thereby to various macroscopic thermodynamic quantities such as the isothermal compressibility and partial molar volumes. Recently, the field has seen a strong revival with breakthroughs in the numerical computation of KBIs and applications to complex systems such as bio-molecules. One of the main emergent results is the possibility to use the finite volume KBIs as a tool to access finite volume thermodynamic quantities. The purpose of this Perspective is to shed new light on the latest developments and discuss future avenues.

Force density functional theory in- and out-of-equilibrium

When a fluid is subject to an external field, as is the case near an interface or under spatial confinement, then the density becomes spatially inhomogeneous. Although the one-body density provides much useful information, a higher level of resolution is provided by the two-body correlations. These give a statistical description of the internal microstructure of the fluid and enable calculation of the average interparticle force, which plays an essential role in determining both the equilibrium and dynamic properties of interacting fluids. We present a theoretical framework for the description of inhomogeneous (classical) many-body systems, based explicitly on the two-body correlation functions. By consideration of local Noether-invariance against spatial distortion of the system we demonstrate the fundamental status of the Yvon-Born-Green (YBG) equation as a local force-balance within the fluid. Using the inhomogeneous Ornstein-Zernike equation we show that the two-body correlations are density functionals and, thus, that the average interparticle force entering the YBG equation is also a functional of the one-body density. The force-based theory we develop provides an alternative to standard density functional theory for the study of inhomogeneous systems both in- and out-of-equilibrium. We compare force-based density profiles to the results of the standard potential-based (dynamical) density functional theory. In-equilibrium, we confirm both analytically and numerically that the standard approach yields profiles that are consistent with the compressibility pressure, whereas the force-density functional gives profiles consistent with the virial pressure. For both approaches we explicitly prove the hard-wall contact theorem that connects the value of the density profile at the hard-wall with the bulk pressure. The structure of the theory offers deep insights into the nature of correlation in dense and inhomogeneous systems.

Adaptive Brownian Dynamics

A framework for performant Brownian Dynamics (BD) many-body simulations with adaptive timestepping is presented. Contrary to the Euler-Maruyama scheme in common non-adaptive BD, we employ an embedded Heun-Euler integrator for the propagation of the overdamped coupled Langevin equations of motion. This enables the derivation of a local error estimate and the formulation of criteria for the acceptance or rejection of trial steps and for the control of optimal stepsize. Introducing erroneous bias in the random forces is avoided by rejection sampling with memory due to Rackauckas and Nie, which makes use of the Brownian bridge theorem and guarantees the correct generation of a specified random process even when rejecting trial steps. For test cases of Lennard-Jones fluids in bulk and in confinement, it is shown that adaptive BD solves performance and stability issues of conventional BD, already outperforming the latter even in standard situations. We expect this novel computational approach to BD to be especially helpful in long-time simulations of complex systems, e.g., in non-equilibrium, where concurrent slow and fast processes occur.

Reduced variance analysis of molecular dynamics simulations by linear combination of estimators

Building upon recent developments of force-based estimators with a reduced variance for the computation of densities, radial distribution functions, or local transport properties from molecular simulations, we show that the variance can be further reduced by considering optimal linear combinations of such estimators. This control variates approach, well known in statistics and already used in other branches of computational physics, has been comparatively much less exploited in molecular simulations. We illustrate this idea on the radial distribution function and the one-dimensional density of a bulk and confined Lennard-Jones fluid, where the optimal combination of estimators is determined for each distance or position, respectively. In addition to reducing the variance everywhere at virtually no additional cost, this approach cures an artifact of the initial force-based estimators, namely, small but non-zero values of the quantities in regions where they should vanish. Beyond the examples considered here, the present work highlights, more generally, the underexplored potential of control variates to estimate observables from molecular simulations.

Convergence behaviour of solvation shells in simulated liquids

A convenient way to analyse solvent structure around a solute is to use solvation shells, whereby solvent position around the solute is discretised by the size of a solvent molecule, leading to multiple shells around the solute. The two main ways to define multiple shells around a solute are either directly with respect to the solute, called solute-centric, or locally for both solute and solvent molecules alike. It might be assumed that both methods lead to solvation shells with similar properties. However, our analysis suggests otherwise. Solvation shells are analysed in a series of simulations of five pure liquids of differing polarity. Shells are defined locally working outwards from each molecule treated as a reference molecule using two methods: the cutoff at the first minimum in the radial distribution function and the parameter-free Relative Angular Distance method (RAD). The molecular properties studied are potential energy, coordination number and coordination radius. Rather than converging to bulk values, as might be expected for pure solvents, properties are found to deviate as a function of shell index. This behaviour occurs because molecules with larger coordination numbers and radius have more neighbours, which make them more likely to be connected to the reference molecule via fewer shells. The effect is amplified for RAD because of its more variable coordination radii and for water with its more open structure and stronger interactions. These findings indicate that locally defined shells should not be thought of as directly comparable to solute-centric shells or to distance. As well as showing how box size and cutoff affect the non-convergence, to restore convergence we propose a hybrid method by defining a new set of shells with boundaries at the uppermost distance of each locally derived shell.

Shear-induced deconfinement of hard disks

Using Brownian dynamics simulations, we investigate the response to shear of a two-dimensional system of quasi-hard disks that are confined in the velocity gradient direction by a smooth external potential. Shearing the confined system leads to a homogenization of the one-body density profile. In order to rationalize this deconfinement effect, we split the internal one-body force field into adiabatic and superadiabatic contributions. We demonstrate that the superadiabatic force field consists of viscous and of structural contributions. We give an empirical scaling law that yields results for the superadiabatic force profiles both in the flow and in the gradient direction, in excellent agreement with the simulation data.

Kirkwood-Buff Integrals Using Molecular Simulation: Estimation of Surface Effects

Kirkwood-Buff (KB) integrals provide a connection between microscopic properties and thermodynamic properties of multicomponent fluids. The estimation of KB integrals using molecular simulations of finite systems requires accounting for finite size effects. In the small system method, properties of finite subvolumes with different sizes embedded in a larger volume can be used to extrapolate to macroscopic thermodynamic properties. KB integrals computed from small subvolumes scale with the inverse size of the system. This scaling was used to find KB integrals in the thermodynamic limit. To reduce numerical inaccuracies that arise from this extrapolation, alternative approaches were considered in this work. Three methods for computing KB integrals in the thermodynamic limit from information of radial distribution functions (RDFs) of finite systems were compared. These methods allowed for the computation of surface effects. KB integrals and surface terms in the thermodynamic limit were computed for Lennard–Jones (LJ) and Weeks–Chandler–Andersen (WCA) fluids. It was found that all three methods converge to the same value. The main differentiating factor was the speed of convergence with system size L. The method that required the smallest size was the one which exploited the scaling of the finite volume KB integral multiplied by L. The relationship between KB integrals and surface effects was studied for a range of densities.

Custom flow in overdamped Brownian dynamics

When an external field drives a colloidal system out of equilibrium, the ensuing colloidal response can be very complex, and obtaining a detailed physical understanding often requires case-by-case considerations. To facilitate systematic analysis, here we present a general iterative scheme for the determination of the unique external force field that yields prescribed inhomogeneous stationary or time-dependent flow in an overdamped Brownian many-body system. The computer simulation method is based on the exact one-body force balance equation and allows to specifically tailor both gradient and rotational velocity contributions, as well as to freely control the one-body density distribution. Hence, compressibility of the flow field can be fully adjusted. The practical convergence to a unique external force field demonstrates the existence of a functional map from both velocity and density to external force field, as predicted by the power functional variational framework. In equilibrium, the method allows to find the conservative force field that generates a prescribed target density profile, and hence implements the Mermin-Evans classical density functional map from density distribution to external potential. The conceptual tools developed here enable one to gain detailed physical insight into complex flow behaviour, as we demonstrate in prototypical situations.

Homogeneous and heterogeneous populations of active rods in two-dimensional channels

Active swarms, consisting of individual agents which consume energy to move or produce work, are known to generate a diverse range of collective behaviors. Many examples of active swarms are biological in nature (e.g., fish shoals and bird flocks) and have been modeled extensively by numerical simulations. Such simulations of swarms usually assume that the swarm is homogeneous; that is, every agent has exactly the same dynamical properties. However, many biological swarms are highly heterogeneous, such as multispecies communities of micro-organisms in soil, and individual species may have a wide range of different physical properties. Here we explore heterogeneity by developing a simple model for the dynamics of a swarm of motile heterogeneous rodlike bacteria in the absence of hydrodynamic effects. Using molecular dynamics simulations of active rods confined within a two-dimensional rectangular channel, we first explore the case of homogeneous swarms and show that the key parameter governing both dynamics is ratio of the motility force to the steric force. Next we explore heterogeneous or mixed swarms in which the constituent self-propelled rods have a range of motilities and steric interactions. Our results show that the confining boundaries play a strong role in driving the segregation of mixed populations.

Structural Nonequilibrium Forces in Driven Colloidal Systems

We identify a structural one-body force field that sustains spatial inhomogeneities in nonequilibrium overdamped Brownian many-body systems. The structural force is perpendicular to the local flow direction, it is free of viscous dissipation, it is microscopically resolved in both space and time, and it can stabilize density gradients. From the time evolution in the exact (Smoluchowski) low-density limit, Brownian dynamics simulations, and a novel power functional approximation, we obtain a quantitative understanding of viscous and structural forces, including memory and shear migration.