# Non-equilibrium Molecular Dynamics: Algorithms, Analysis, and Applications, May 2023

- Location
- Alan Walters 111 - Seminar room 2

- Dates
- Thursday 4 May (11:00) - Friday 5 May 2023 (13:00)

- Contact

The workshop will bring together both leading experts and early career researchers on (non)equilibrium molecular dynamics and beyond.

We acknowledge financial support from the Royal Society, Heilbronn Institute for Mathematical Research, QJMAM Fund for Applied Mathematics from the Institute of Mathematics and its Applications, Institute for Interdisciplinary Data Science and AI at the University of Birmingham, and University of Birmingham Digital Chemistry Network+.

**Schedule:**

Day 1 (4 May):

11:00-11:45 **Gabriel Stoltz** (École des Ponts ParisTech): Error estimates and variance reduction for nonequilibrium stochastic dynamics

11:45-12:10 **Xinyi Wu** (University of Birmingham): Computation of transport coefficients in dissipative particle dynamics

12:10-13:45 Lunch Break

13:45-14:30 **Livia Bartók-Pártay** (University of Warwick): Nested sampling: what can we learn about an interatomic potential?

14:30-15:00 Coffee Break

15:00-15:25 **Noé Blassel** (École des Ponts ParisTech): Stochastic Norton dynamics

15:25-15:50 **Renato Spacek** (École des Ponts ParisTech): Extending the regime of linear response with synthetic forcings

15:50-16:15 **Régis Santet** (École des Ponts ParisTech): Optimizing the diffusion for overdamped Langevin dynamics

16:15-16:40 **Shiva Darshan** (École des Ponts ParisTech): Sticky coupling as a control variate for the computation of transport coefficients

Day 2 (5 May):

9:30-10:15 **Urbain Vaes** (Ecole des Ponts): Optimal importance sampling for overdamped Langevin dynamics

10:15-11:00 **Lisa Maria Kreusser** (University of Bath): Interacting particle models: From applied and numerical analysis to applications

11:00-11:30 Coffee Break

11:30-12:15 **Greg Pavliotis** (Imperial College, London): Inference for mean field SDEs

12:15-12:40 **Katerina Karoni** (University of Edinburgh): Friction-adaptive descent: a new family of optimization algorithms

**Titles and abstracts**

**Speaker:** Gabriel Stoltz

**Title:** Error estimates and variance reduction for nonequilibrium stochastic dynamics

**Abstract:** Equilibrium properties in statistical physics are obtained by computing averages with respect to Boltzmann-Gibbs measures, sampled in practice using ergodic dynamics such as the Langevin dynamics. Some quantities however cannot be computed by simply sampling the Boltzmann-Gibbs measure, in particular transport coefficients, which relate the current of some physical quantity of interest with the forcing needed to induce it. For instance, a temperature difference induces an energy current, the proportionality factor between these two quantities being the thermal conductivity. From an abstract point of view, transport coefficients can also be considered as some form of sensitivity analysis with respect to an added forcing to the baseline dynamics.

There are various numerical techniques to estimate transport coefficients, which all suffer from large errors, in particular large statistical errors. I will review the most popular methods, namely the Green-Kubo approach where the transport coefficient is rewritten as some time-integrated correlation function, and the approach based on longtime averages of the stochastic dynamics perturbed by an external driving (so-called nonequilibrium molecular dynamics). I will make precise in each case the various sources of errors, in particular the bias related to the time discretization of the underlying continuous dynamics, and the variance of the associated Monte Carlo estimators. I will also briefly present some recent alternative techniques to estimate transport coefficients.

**Speaker:** Xinyi Wu

**Title:** Computation of transport coefficients in dissipative particle dynamics

**Abstract:** Dissipative particle dynamics (DPD) is a momentum-conserving thermostat that has been widely used in complex fluids and soft matter at mesoscales. In this talk, I will give an overview of DPD and its corresponding numerical methods. In particular, I will compare the performance of various approaches to estimating transport properties (e.g., mobility and self-diffusion), including those based on Green-Kubo formulas, mean squared displacement, and linear response theory.

**Speaker:** Livia Bartók-Pártay

**Title:** Nested sampling: what can we learn about an interatomic potential?

**Abstract:** In recent years we have been working on adapting the Bayesian statistical approach, nested sampling, for studying atomistic systems. Nested sampling automatically generates all the relevant atomic configurations, unhindered by high barriers, and one of its most appealing advantages is that the global partition function can be calculated very easily, thus thermodynamic properties, such as the heat capacity or compressibility becomes accessible. Nested sampling provides an unbiased and exhaustive exploration, starting from the high energy region of the potential energy landscape (gas phase configurations) and progressing towards the ground state structure (crystalline solid) through a series of nested energy levels, estimating the corresponding phase space volume of each. This way the method samples the different basins proportional to their volume, and instead of providing an exhaustive list of the local minima, it identifies the thermodynamically most relevant states without any prior knowledge of the structures or phase transitions.

The use and advantages of the nested sampling method has been demonstrated in sampling the potential energy landscape of several model systems, allowing us to calculate the complete pressure-temperature phase diagram. These often result in the discovery of previously unknown thermodynamically stable solid phases, highlighting that the macroscopic properties of interatomic potential models can be very different from expected or intended.

**Speaker:** Noé Blassel

**Title:** Stochastic Norton dynamics

**Abstract:** Transport coefficients quantify the sensitivity of a flux to the magnitude of a force driving a system out of thermodynamic equilibrium, in the small force regime in which the response is linear. Examples of these include the mobility, thermal conductivity or shear viscosity. Whereas, at the macroscopic level, forces and fluxes play symmetric roles, standard approaches in computational statistical physics apply a constant force and measure the average response flux. In contrast, the Norton method, proposed by Evans et al., proposes a dual approach, which consists in fixing the flux, and measuring the average magnitude of the force needed to induce it. At the dynamical level, this amounts to consider dynamics constrained to remain on a submanifold of the phase space, defined by the constant-flux condition. However, the constraint is typically non-holonomic, thus existing techniques for constrained dynamics cannot be directly applied. The original method is defined in a deterministic setting in which crucial properties such as ergodicity may fail. In this work we propose a more general and stochastic formulation of the Norton method, also providing evidence that it is both valid and sometimes more efficient than the standard approach to transport coefficient computations.

**Speaker:** Renato Spacek

**Title:** Extending the regime of linear response with synthetic forcings

**Abstract:** Transport coefficients, such as the mobility, thermal conductivity and shear viscosity, are quantities of prime interest in statistical physics. At the macroscopic level, transport coefficients relate an external forcing of magnitude η, with η ≪ 1, acting on the system to an average response expressed through some steady-state flux. In practice, steady-state averages involved in the linear response are computed as time averages over a realization of some stochastic differential equation. Variance reduction techniques are of paramount interest in this context, as the linear response is scaled by a factor of 1/η, leading to large statistical error. One way to limit the increase in the variance is to allow for larger values of η by increasing the range of values of the forcing for which the nonlinear part of the response is sufficiently small. In theory, one can add an extra forcing to the physical perturbation of the system, called synthetic forcing, as long as this extra forcing preserves the invariant measure of the reference system. The aim is to find synthetic perturbations allowing to reduce the nonlinear part of the response as much as possible. In this talk, I will present a mathematical framework for quantifying the quality of synthetic forcings, in the context of linear response theory, and discuss various possible choices for them. I will illustrate my analysis with numerical results in low-dimensional systems.

**Speaker:** Régis Santet

**Title:** Optimizing the diffusion for overdamped Langevin dynamics

**Abstract:** Predicting properties of materials and macroscopic physical systems in the framework of statistical physics, or obtaining the distribution of parameter values for Bayesian inference problems, both rely on sampling high dimensional probability measures.

Popular methods rely on stochastic dynamics, in particular Markov Chain Monte Carlo (MCMC) methods. We focus on the overdamped Langevin dynamics, with multiplicative noise (i.e. a position dependent noise), which is ergodic for the target Boltzmann measure. Our objective is to compute, explicitly or numerically, the optimal diffusion function leading to the fastest convergence rate, measured here in terms of the spectral gap of the generator of the dynamics. In particular, we discuss the normalization of the diffusion coefficient, since a too large diffusion coefficient has to be compensated by very small time steps.

More precisely, we formalize the maximization of the convergence rate of the overdamped Langevin dynamics with respect to the diffusion coefficient *D* as a convex optimization program, and we show its well-posedness. We also propose in low dimensional scenarios a numerical procedure combining a finite element parameterization and an optimization algorithm, to compute the optimal diffusion in practice. The numerical procedure is illustrated on simple one-dimensional examples which show the benefits of having a position-dependent diffusion coefficient. We also study the behaviour of the optimal diffusion in the homogenized limit and show it has an analytical expression, proportional to the inverse of the target density, which is in accordance with various previous heuristics. This simple limiting behaviour can be used as an initial guess in an optimization algorithm, or as a proxy for the optimal diffusion which does not require costly convex optimization procedures.

**Speaker:** Shiva Darshan

**Title:** Sticky coupling as a control variate for the computation of transport coefficients

**Abstract:** A standard method to compute transport coefficients is to simulate Langevin dynamics perturbed by a small non-equilibrium forcing up to a time *T* and time-average over the trajectory of a desired observable divided by the magnitude of the forcing, η. Unfortunately, this method suffers from large finite-time sampling bias and variance in the limit of small forcing, on the order of (*T*η)^{-1} and (*T*η^{2})^{-1} respectively.

For overdamped Langevin dynamics, we propose a method to reduce the bias and variance of this computation using a version of the reference (unperturbed) dynamics sticky coupled to the perturbed dynamics as a control variate. We will show that when the potential of the dynamics is strongly convex at infinity, this sticky coupling based estimator reduces the bias and variance by a factor of η^{-1} compared to the standard method. The case of strongly convex at infinity potentials includes commonly used systems such as Lennard-Jones particles confined to box by a quadratic potential.

**Speaker:** Urbain Vaes

**Title:** Optimal importance sampling for overdamped Langevin dynamics

**Abstract:** We present and study an importance sampling approach for calculating averages with respect to multimodal probability distributions. Traditional Markov chain Monte Carlo methods to this end, which are based on time averages along a realization of a Markov process ergodic with respect to the target probability distribution, are usually plagued by a large variance due to the metastability of the process. The estimator we study is based on an ergodic average along a realization of an overdamped Langevin process for a modified potential. We obtain an explicit expression for the optimal perturbation potential in dimension 1 and propose a general numerical approach for approximating the optimal potential in the multi-dimensional setting.

**Speaker:** Lisa Maria Kreusser

**Title:** Interacting particle models: From applied and numerical analysis to applications

**Abstract:** The recent, rapid advances in modern biology and data science have opened up a whole range of challenging mathematical problems which can be tackled by combining applied and numerical analysis. In this talk, I will discuss a class of interacting particle models with anisotropic interaction forces and the associated continuum limit. These models are motivated by the simulation of fingerprint patterns, of which large databases are required in forensic science and biometric applications. The central novelty in the models I consider is an anisotropy induced by an underlying tensor field. This innovation does not only lead to the ability to describe real-world phenomena more accurately but also renders their analysis significantly harder compared to their isotropic counterparts. I will present a positivity-preserving numerical method for simulating realistic fingerprint patterns efficiently with the continuum model. I will also outline a mean-field optimal control algorithm for an inverse problem arising in the context of parameter identification. Going from interaction to transport networks, I will show recent numerical results for a biological network formation problem, and I will demonstrate how similar models can be used in semi-supervised learning and network science.

**Speaker:** Greg Pavliotis

**Title:** Inference for mean field SDEs

**Abstract:** Nonlinear SDEs in the sense of McKean, for which the drift (and possibly the diffusion) coefficient depends on the law of the process, appear in many applications, ranging from plasma physics to social dynamics. These SDEs can be obtained in the mean field limit of systems of interacting particle systems driven by noise. Often one needs to estimate parameters in the drift of the McKean SDE from data. In this talk I will present recently developed techniques for learning parameters in mean field models from observations of particle trajectories. These techniques are based on the stochastic gradient descent algorithm, eigenfunction martingale estimators and the method of moments.

**Speaker:** Katerina Karoni

**Title:** Friction-adaptive descent: a new family of optimization algorithms

**Abstract:** We describe a family of descent algorithms which generalizes common existing schemes used in optimization and neural network training. We create models which can directly address certain situations arising in standard algorithms, for example circumventing barriers delaying approach to minima, exploring rare intermediate states, or stabilizing numerical schemes. Our algorithms are simple to implement and control and convergence can be shown via Lyapunov functions, building on previous works.