Section: New Results

Molecular dynamics and related problems

Participants : Matthew Dobson, Bradley Dickson, Mustapha Hamdi, Claude Le Bris, Frédéric Legoll, Tony Lelièvre, Kimiya Minoukadeh, Stefano Olla, David Pommier, Raphaël Roux, Gabriel Stoltz.

The extremely broad field of Molecular dynamics is a domain where the MICMAC project-team, originally more involved in the quantum chemistry side, has invested a lot of efforts in the recent years. Molecular dynamics may also be termed 'computational statistical physics' since the main aim is to numerically estimate average properties of materials as given by the laws of statistical physics. The project-team studies both deterministic and probabilistic techniques used in the field. On these topics, we benefited from funding from the ARC Hybrid and the ANR MEGAS (“MEthodes Géométriques et échantillonnage : Application à la Simulation moléculaire”). The habilitation thesis [7] of T. Lelièvre has been defended.

Highly oscillatory dynamics.

Constant energy averages are often computed as long time limits of time averages along a typical trajectory of the Hamiltonian dynamics. One difficulty of such a computation is the presence of several time scales in the dynamics: the frequencies of some motions are very high (e.g. for the atomistic bond vibrations), while those of other motions are much smaller. Actually, fast phenomena are only relevant through their mean effect on the slow phenomena, and their precise description is not needed. Consequently, there is a need for time integration algorithms that take into account these fast phenomena only in an averaged way, and for which the time step is not restricted by the highest frequencies.

As a follow-up on a previous study [34] along these lines, M. Dobson, C. Le Bris and F. Legoll continue work that applies homogenization techniques to the Hamilton-Jacobi PDE associated to the Hamiltonian ODE. This yields algorithms whose step sizes are independent of the fastest frequencies. The current goal is to extend the applicability of the previous analysis, by removing certain structural assumptions made on the fast forces. This will improve the generality of the algorithm and enable comparisons with a wide range of related algorithms in the dynamical systems literature.

We collaborate on this subject with F. Castella, P. Chartier and E. Faou from INRIA Rennes, with the funding of ANR MEGAS and ARC Hybrid.

Ergodicity of sampling methods.

In many cases, the dynamics of the system under study is restrained to some submanifold of the whole accessible space. A famous instance is the Hamiltonian dynamics, for which the energy of the system is constant. Hamiltonian dynamics is useful for computing average properties assuming ergodicity, but this dynamics is sometimes not ergodic due to spurious invariants of the motion. It is then appealing to resort to projected stochastic dynamics, which are more likely to destroy these spurious invariants. Such a scheme has been analyzed by E. Faou (INRIA Rennes) and T. Lelièvre in [27] and rates of convergence have also been provided.

Despite the success of stochastic techniques, deterministic sampling methods, such as the Nosé-Hoover dynamics, are still often used in the applied community to compute canonical averages. In collaboration with M. Luskin and R. Moeckel (University of Minnesota), F. Legoll has studied the Nosé-Hoover dynamics when applied to completely integrable systems. Using averaging and KAM techniques, it has been showed that the dynamics is actually not ergodic with respect to the Gibbs measure [37] . This extends a previous work that addressed the simple case of the one-dimensional harmonic oscillator.

Free energy computations.

For large molecular systems, the information of the whole configuration space may be summarized in a few coordinates of interest, called reaction coordinates. An important problem in chemistry or biology is to compute the effective energy felt by those reaction coordinates, called free energy. T. Lelièvre and G. Stoltz, together with M. Rousset, are currently finishing a review book on free energy computations [55] .

Besides this review work, we continued our systematic studies on the new class of adaptive methods:

• The Adaptive Biasing Force (ABF) method is a stochastic algorithm used to compute this free energy. It is based upon a nonlinear dynamics, which uses the mean force as a biasing force to prevent the system from being trapped in metastable regions. The nonlinearity in the dynamics comes from a conditional expectation computed with respect to the solution. The convergence of the algorithm can be accelerated by using multiple walkers, each following similar dynamics but driven by independent Brownian motions. The use of multiple walkers allows for further improvement via a selection mechanism, whereby the walkers are weighted according to a given fitness function and resampled at fixed time intervals. T. Lelièvre and K. Minoukadeh, together with C. Chipot (University of Illinois, on leave from Université Henri Poincaré - Nancy 1), have shown, in [41] , the applicability of the method to realistic biological systems. Their work has in particular highlighted cases in which the standard ABF method, using a single walker, fails to converge within reasonable time scales;

• In [32] , B. Jourdain, T. Lelièvre and R. Roux study the convergence of a particle approximation of the Adaptive Biasing Force method;

• in addition, a current work in progress by B. Dickson, T. Lelièvre, G. Stoltz and F. Legoll, in collaboration with P. Fleurat-Lessard (Département de chimie, ENS Lyon) is to improve adaptive methods where the biasing potential is updated (in contrast to the update of the derivative of this potential as for ABF), see [52] . This study is supported by the ANR Sire;

• In [51] , N. Chopin (CREST, ENSAE), T. Lelièvre and G. Stoltz are studying the relevance of adaptive methods to sample posterior distributions in Bayesian statistics. The preliminary results are very encouraging.

• On a more applicative side, M. Hamdi currently studies the biological relevance of functionalized carbon nanotubes for targetted drug delivery, using the ABF methodology to compute the free energy profile and the transition mechanism associated with the penetration of the nanotube in the cell lipidic membrane.

The thermodynamic integration method, and nonequilibrium methods in the Jarzynski fashion for Langevin dynamics, have been extended in [54] .

As a by-product of these free-energy studies, T. Lelièvre has obtained in [38] a new result for proving a logarithmic Sobolev inequality on a measure defined on , assuming that a logarithmic Sobolev inequality holds for the marginals * and the conditional measures (·|) associated to some function (with m<n ). This theoretical result has practical interest in molecular dynamics, where is the reaction coordinate and where the above assumptions are often met in practice.

The free energy completely describes the statistics of the reaction coordinates. F. Legoll and T. Lelièvre have been working on the definition of a dynamics closed in these reaction coordinates. The problem hence amounts to reducing the dimension of a set of SDEs, from the full set of degrees of freedom to only a small subset of them. Encouraging numerical results have been obtained [36] , along with estimates on the accuracy on the proposed effective dynamics (again using entropy techniques). Extensions of this work to more realistic molecular systems, as well as to solid materials, are currently under study.

Search for reaction paths.

The microscopic dynamics used to sample the configurations of the system are often trapped in metastable states. A major numerical issue is therefore the search for transition paths connecting metastable states. E. Cancès, F. Legoll, K. Minoukadeh and two of their collaborators at CEA Saclay proposed in [21] an improvement to an existing eigenvector following method, the Activation-Relaxation Technique nouveau (ARTn), for searches of saddle points and transition pathways on a given potential energy surface. Local convergence and robustness of the algorithm have been established, and the new method has been successfully tested on point defects in body centered cubic iron. The possible application of this algorithm to other large scale problems is currently studied, in collaboration with D. Wales (Cambridge).

In addition, D. Pommier and T. Lelièvre are currently studying various methods numerically to compute reactive paths. Mathematically, this amounts to sampling trajectories which start from a metastable region, and end in another metastable region.

Nonequilibrium systems.

The study of nonequilibrium systems in the case of mechanical shock waves has been pursued. A first work dealt with the computation of the so-called Hugoniot curve, which is the locus of all the thermodynamic states which can be attained from a given initial state with shocks of increasing strengths [40] . The numerical method employed is an original adaptive scheme which improves on previous iterative techniques, and which can be used more generally to sample systems with a constraint fixed in average. In [39] , we additionally investigated whether release waves behind shock fronts are isentropic. This was assessed by comparing the numerical results from equilibrium methods based on techniques for free-energy computations.

This year, a work has also started on the thermal conductivity of one-dimensional chains, benefiting from the presence of S. Olla (on leave from Université Paris Dauphine), see [53] . In this classical setting, we can use nonlinear interaction potentials and perturb the dynamics of the system by some energy and momentum preserving noise. The aim is to determine whether this is enough to achieve a finite conductivity. This study is somehow the classical counterpart of the quantum computations [42] , [43] .

Finally, we mention a PhD thesis starting under the cosupervision of T. Lelièvre and A. Ern, on the modelling of clays (funding from ANDRA Agence Nationale pour la gestion des Déchets RadioActifs). The PhD student is R. Joubaud. One part of the project contains a modelling at the molecular level of the diffusion of ions in clays, and we plan to use nonequilibrium methods to understand this process. This work involves F. Legoll, T. Lelièvre and G. Stoltz, within a collaboration with B. Rotenberg and P. Turq (Paris 6, Laboratoire PECSA).