Section: New Results

Computational Statistical Physics

Participants : Grégoire Ferré, Florent Hédin, Frédéric Legoll, Tony Lelièvre, Mouad Ramil, Julien Roussel, Laura Silva Lopes, Gabriel Stoltz, Pierre Terrier.

The objective of computational statistical physics is to compute macroscopic properties of materials starting from a microscopic description, using concepts of statistical physics (thermodynamic ensembles and molecular dynamics). The contributions of the team can be divided into four main topics: (i) the development of methods for sampling the configuration space; (ii) the numerical analysis of such methods; (iii) the efficient computation of dynamical properties which requires to sample metastable trajectories; (iv) coarse-graining techniques to reduce the computational cost of molecular dynamic simulations and gain some insights on the models.

Sampling of the configuration space: new algorithms and applications

New numerical methods in order to sample probability measures on the configuration space have been developed: either measures supported on submanifolds, or stationary states of stochastic dynamics. First, in [51], T. Lelièvre and G. Stoltz, together with M. Rousset (Inria Rennes, France) have studied how to sample probability measures supported on submanifolds, by adding an extra momentum variable to the state of the system, and discretizing the associated Hamiltonian dynamics with some stochastic perturbation in the extra variable. In order to avoid biases in the invariant probability measures sampled by discretizations of these stochastically perturbed Hamiltonian dynamics, a Metropolis rejection procedure can be considered. The so-obtained scheme belongs to the class of generalized Hybrid Monte Carlo (GHMC) algorithms. However, the usual method has to be generalized using a procedure suggested by Goodman, Holmes-Cerfon and Zappa for Metropolis random walks on submanifolds, where a reverse projection check is performed to enforce the reversibility of the algorithm for large timesteps and hence avoid biases in the invariant measure. A full mathematical analysis of such procedures is provided, as well as numerical experiments demonstrating the importance of the reverse projection check on simple toy examples. Second, the work [55] by J. Roussel and G. Stoltz focuses on the use of control variates for non-equilibrium systems. Whereas most variance reduction methods rely on the knowledge of the invariant probability measure, this latter is not explicit out of equilibrium. Control variates offer an attractive alternative in this framework. J. Roussel and G. Stoltz have proposed a general strategy for constructing an efficient control variate, relying on physical simplifications of the dynamics. The authors provide an asymptotic analysis of the variance reduction in a perturbative framework, along with extensive numerical tests on three different systems.

In terms of applications of such sampling techniques, members of the project-team have been working on two different subjects: random matrices models and adaptive techniques to compute large deviation rate functionals. The paper [16] was written by G. Ferré and D. Chafaï (Université Paris Dauphine, France), following the simple idea: the eigenvalues of random matrices are distributed according to Boltzmann–Gibbs measures, but researchers in this field do not use techniques from statistical physics for numerical investigations. The authors therefore used a Hamiltonian Monte Carlo algorithm to investigate numerically conjectures about random matrices and related Coulomb gases. The next step is to add constraints to these systems to understand better the behavior of random matrices with constraints and the large size limit of their spectra (the algorithm mentioned above to sample probability measures supported on submanifolds may be useful in this context). The work [19] focuses on computing free energies and entropy functions, as they arise in large deviations theory, through adaptive techniques. It is actually in the spirit of techniques used in mathematical finance, adapted to the statistical mechanics context, and enriched with new estimators based on variational representations of entropy functions. These tools have been pioneered by H. Touchette (Stellenbosch University, South Africa), with whom the paper was written by G. Ferré.

Sampling of the configuration space: numerical analysis

Concerning the numerical analysis of sampling techniques of probability measures on the configuration space, let us mention three works.

First, in [44], G. Ferré and G. Stoltz study the numerical errors that arise when a stochastic differential equation (SDE) is discretized in order to compute scaled cumulant functions (or free energy) and ergodic properties of Feynman–Kac semigroups. These quantities naturally arise in large deviations theory, for estimating probabilities of rare events. This analysis is made difficult by the nonlinear (mean field) feature of the dynamics at hand. The obtained estimates generalize previous results on the numerical analysis of ergodic properties of discretized SDEs. As a theoretical extension of the previous work, the purpose of the work [43] by G. Ferré and G. Stoltz, in collaboration with M. Rousset (Inria Rennes, France), is to provide further theoretical investigations on the long time behavior of Feynman–Kac semigroups. More precisely, it aims at giving practical criteria for these nonlinear semigroups to have a limit, and makes precise in which sense this limit is to be understood. This was an open problem so far for systems evolving in unbounded configuration spaces, which was addressed through Lyapunov function techniques. Although theoretical, these results are of practical importance since, if these dynamics do not have a well-defined long time behavior, it is hopeless to try to compute rare events.

Finally, together with C. Andrieu (Univ. Bristol, United-Kingdom), A. Durmus (ENS Saclay, France) and N. Nüsken (Univ. Potsdam, Germany), J. Roussel derived in [32] spectral gap estimates for several Piecewise Deterministic Markov Processes (PDMPs), namely the Randomized Hamiltonian Monte Carlo, the Zig-Zag process and the Bouncy Particle Sampler. The hypocoercivity technique provides estimates with explicit dependence on the parameters of the dynamics. Moreover the general framework considered allows to compare quantitatively the bounds found for the different methods. Such PDMDs are currently more and more used as efficient sampling tools, but their theoretical properties are still not yet well understood.

Sampling of dynamical properties and rare events

The sampling of dynamical properties along molecular dynamics trajectories is crucial to get access to important quantities such as transition rates or reactive paths. This is difficult numerically because of the metastability of trajectories. Members of the project-team are following two numerical approaches to sample metastable trajectories: the accelerated dynamics à la A.F. Voter and the adaptive multilevel splitting (AMS) technique to sample reactive paths between metastable states.

Concerning the mathematical analysis of the accelerated dynamics, in [50], T. Lelièvre reviews the recent mathematical approaches to justify these numerical methods, using the notion of quasi-stationary distribution. Moreover, in [49], T. Lelièvre together with D. Le Peutrec (Université de Paris Saclay, France) and G. Di Gesu and B. Nectoux (TU Wien, Austria) give an overview of the results obtained during the PhD of B. Nectoux. Using the quasi-stationary distribution approach and tools from semi-classical analysis, one can justify the use of kinetic Monte Carlo models parametrized by the Eyring-Kramers formulas to describe exit events from metastable states, for the overdamped Langevin dynamics. Concerning the implementation, in [22], F. Hédin and T. Lelièvre test the Generalized Parallel Replica algorithm to biological systems, and obtain strong linear scalability, providing up to 70% of the maximum possible speedup on several hundreds of CPUs. The “Parallel Replica” (ParRep) dynamics is known for allowing to simulate very long trajectories of metastable Langevin dynamics in the materials science community, but it relies on assumptions that can hardly be transposed to the world of biochemical simulations. The later developed “Generalized ParRep” variant solves those issues, and it had not been applied to significant systems of interest so far. Finally, let us mention the work [27] where T. Lelièvre together with J. Reygner (Ecole des Ponts, France) and L. Pillaud-Vivien (Inria Paris, France) analyze mathematically the Fleming-Viot particle process in the simple case of a finite state space. This Fleming-Viot particle process is a key ingredient of the Generalized ParRep algorithm mentioned above, in order to both approximate the convergence time to the quasi-stationary distribution, and to efficiently sample it.

Concerning the AMS technique, in [36], T. Lelièvre and C.-E. Bréhier (ENS Lyon, France) test new importance functions to compute rare events associated with the law of the solution to a stochastic differential equation at a given fixed time. This can be used for example to estimate the rate functional for large deviation principle applied to time averages.

Coarse-graining

In two related works, members of the project-team have studied the quality of the effective dynamics derived from a high dimension stochastic differential equation on a few degrees of freedom, using a projection approach à la Mori-Zwanzig. More precisely, in [48], F. Legoll, T. Lelièvre and U. Sharma obtain precise error bounds in the case of non reversible dynamics. This analysis also aims at discussing what is a good notion of mean force for non reversible systems. In [53], T. Lelièvre together with W. Zhang (ZIB, Germany) extend previous results on pathwise error estimates for such effective dynamics to the case of nonlinear vectorial reaction coordinates.

Once a good coarse-grained model has been obtained, one can try to use it in order to get a better integrator of the original dynamic in the spirit of a predictor-corrector method. In [52], T. Lelièvre together with G. Samaey and P. Zielinski (KU Leuven, Belgium) analyze such a micro-macro acceleration method for the Monte Carlo simulation of stochastic differential equations with time-scale separation between the (fast) evolution of individual trajectories and the (slow) evolution of the macroscopic function of interest.