Ecole ECODOQUI: résumé du cours

Dynamique moléculaire et calcul d’énergie libre

T. Lelièvre

27 novembre 2008

Résumé
Ces notes portent sur l’étude des méthodes numériques pour le calcul de diffé-
rences d’énergie libre en dynamique moléculaire [7, 26]. En dynamique moléculaire,
les modèles consistent à décrire l’état d’un système par la position (et éventuelle-
ment la vitesse) de particules (typiquement les positions des noyaux dans un système
moléculaire), qui interagissent au travers d’un potentiel. L’objectif est de calculer
des moyennes par rapport à la mesure de Boltzmann-Gibbs associée à ce potentiel
(moyennes dans l’ensemble canonique). Mathématiquement, il s’agit d’un problème
d’échantillonnage de mesures métastables (ou multi-modales), en très grande dimen-
sion. La particularité de la dynamique moléculaire est que, bien souvent, on a quelques
informations sur les "directions de métastabilité" au travers de coordonnées de ré-
action. En utilisant cette donnée, de nombreuses méthodes ont été proposées pour
permettre l’échantillonnage de la mesure de Boltzmann-Gibbs. Dans une série de tra-
vaux [29, 9, 24, 15], nous avons analysé les méthodes basées sur des équations différen-
tielles stochastiques avec contraintes (dont les solutions vivent sur des sous-variétés
définies comme des lignes de niveaux de la coordonnée de réaction). Il s’agit en fait
d’analyser des méthodes d’échantillonnage de mesures définies sur des sous-variétés
de grande dimension. Plus récemment, nous avons étudié dans [30, 31, 28, 22] des
méthodes adaptatives qui ont été proposées pour se débarasser des métastabilités.
Mathématiquement, il s’agit de méthodes d’échantillonnage d’importance, avec une
fonction d’importance qui est calculée au cours de la simulation de manière adapta-
tive. Nous avons étudié la vitesse de convergence vers la mesure d’équilibre pour ces
méthodes adaptatives, en utilisant des méthodes d’entropie. Nous avons proposé de
nouvelles méthodes numériques à la communauté appliquée pour utiliser au mieux ces
idées, qui ont été adoptées par des praticiens du domaine [5].
On considère un système moléculaire de N particules, de position (x1 , . . . , xN ) = x ∈
R3N et qui interagissent au travers d’un potentiel V (x1 , . . . , xN ). On note n = 3N la di-
mension totale. Dans l’ensemble canonique NVT, les positions des particules se répartissent
suivant la mesure de probabilité de Boltzmann-Gibbs :

µ(dx) = Z −1 exp(−βV (x)) dx, (1)

où Z = exp(−βV (x)) dx est la fonction de partition et β = (kT )−1 , avec k la constante

R

de Boltzmann et T la température. Un des objectifs essentiels en dynamique moléculaire

est de calculer des moyennes canoniques de la forme
Z
φ(x)µ(dx) (2)

pour certaines fonctionnelles φ, afin d’obtenir des “quantités macroscopiques” comme la

vraisemblance d’une conformation moléculaire, ou l’avancement d’une réaction chimique.

1
Les domaines d’application concernent typiquement la biologie moléculaire (calcul de re-
pliement de protéine, simulation de mécanismes moléculaires d’intérêt biologique) et la
chimie computationnelle (calcul des processus réactionnels : adsorption sur un substrat,
états de transitions, etc...). On renvoie aux monographes [1, 7, 16, 27] pour une présentation
générale.

Remarque 1 On utilise en général des approximations pour le potentiel V , appelées champs

de force. Typiquement, V est la somme de potentiels modélisant l’interaction entre deux
particules, trois particules et quatre particules :
X X X
V = V1 (xi , xj ) + V2 (xi , xj , xk ) + V3 (xi , xj , xk , xl ).
i<j i<j<k i<j<k<l
 
σ 12 σ 6



Par exemple, V1 (xi , xj ) = VLJ (|xi − xj |) où VLJ (r) = 4 r − r est le potentiel de
Lennard-Jones.

Le calcul de moyennes canoniques (2) est un problème en grande dimension (N est

couramment de l’ordre de plusieurs centaines de milliers, jusqu’à plusieurs millions pour
les applications en biologie) qui nécessite le recours à des méthodes de Monte Carlo, typi-
quement basées sur des chaînes de Markov (méthodes de type Markov Chain Monte Carlo,
utilisées également en statistique par exemple).
Par exemple, pour échantillonner µ, on peut utiliser Xt solution de l’équation différen-
tielle stochastique : p
dXt = −∇V (Xt ) dt + 2β −1 dWt . (3)
La dynamique (3) est appelée, dans ce contexte, dynamique de gradient ou dynamique de
Langevin amortie (overdamped Langevin). Elle est à la base de toutes les méthodes que
nous considérons dans cette section. On rappelle que sous de bonnes hypothèses sur le
potentiel V , on a la propriété d’ergodicité : pour µ-presque tout X0 ,
Z T Z
1
lim φ(Xt )dt = φ(x)µ(dx). (4)
T →∞ T 0

Remarque 2 La dynamique stochastique la plus utilisée par les physiciens est plutôt la
dynamique de Langevin :

dXt = M −1 Pt dt,
(
p (5)
dPt = −∇V (Xt ) dt − γM −1 Pt dt + 2γβ −1 dWt ,

où M est un tenseur de masse et γ un coefficient de friction. Dans la suite, on se restreint

à la dynamique de Langevin amortie (3), qui est obtenue à partir de la dynamique de
Langevin dans la limite m → 0 ou γ → ∞. La généralisation des résultats ci-dessous à la
dynamique de Langevin n’est pas en général triviale car le générateur infinitésimal associé
à cette dynamique n’est plus elliptique, ce qui rend l’analyse plus compliquée.

Ce qui rend en pratique la méthode de calcul (4) inopérante est la métastabilité du

processus Xt . En effet, la convergence T → ∞ est extrêmement longue à obtenir, du fait
de puits dans le potentiel V dans lesquels le processus Xt reste piégé très longtemps. Pour
donner une image intuitive on peut considérer le potentiel 2-dimensionnel représenté sur la
Figure 1, pour lequel la première coordonnée du processus Xt1 est une variable métastable
(ou variable lente).

2
Fig. 1 – A gauche, les lignes de niveau d’un potentiel métastable, et à droite, une représen-
tation schématique de l’évolution de la première coordonnées d’un processus stochastique
évoluant selon la dynamique (3) dans ce potentiel.

Bien sûr, résoudre ce problème dans toute sa généralité est une tâche impossible (au-
cune méthode numérique ne peut explorer l’ensemble des puits d’un potentiel V en grande
dimension, séparés par des barrières arbitrairement grandes, en un temps de calcul raison-
nable). Nous allons donc ajouter une hypothèse supplémentaire, spécifique au cadre de la
simulation moléculaire. Nous suppons dans la suite que l’on connaît les variables lentes du
système, sous la forme d’une application
ξ : Rn → R ou ξ : Rn → T (6)
appelée coordonnée de réaction, ou variable collective. On note T = R/Z le tore uni-
dimensionel. En pratique, on peut avoir ξ : Rn → Rp (ou ξ : Rn → Tp ) avec p  n, mais
on suppose dans cette présentation p = 1 pour simplifier. L’ensemble des résultats sont
généralisables au cas p > 1. Pour fixer les idées, sur l’exemple de la figure 1, la coordonnée
de réaction serait ξ : (x1 , x2 ) 7→ x1 . En pratique, ξ est typiquement une distance signée
entre deux (groupes de) particules (auquel cas ξ est à valeurs dans R), par exemple pour
indexer l’évolution d’une réaction de synthétisation d’un composé, ou un angle (auquel cas
ξ est à valeurs dans T), par exemple pour caractériser la conformation d’une molécule.
La connaissance de ξ nécessite bien sûr une intuition sur le système considéré. La qualité
des méthodes que nous présentons ci-dessous dépend crucialement du choix d’une bonne
fonction ξ, et déterminer pour un système quelconque une bonne coordonnée de réaction
reste un problème compliqué. Dans la suite, on suppose que ξ est une fonction régulière
telle que |∇ξ| > 0.
A la coordonnée de réaction est associée un potentiel effectif, appelé énergie libre.
Définition 3 L’image de la mesure µ par ξ est Z −1 exp(−βA(z)) dz, où A est appelée
l’énergie libre. On vérifie que
Z 
−1
A(z) = −β ln exp(−βV )|∇ξ| dσΣz = −β −1 ln ZΣz ,
−1
(7)
Σz

où Σz = {x, ξ(x) = z} est une sous-variété (régulière) de Rn , et σΣz est la mesure de

Lebesgue sur Σz .
La formule (7) découle de la formule de la co-aire [2, 14], que l’on peut voir comme une
généralisation de la formule de Fubini (intégration par tranches) : pour toute fonction
régulière φ : Rn → R,
Z Z Z
φ(x)|∇ξ|(x)dx = φ dσΣz dz. (8)
Rn R Σz

3
La mesure |∇ξ|−1 (x)σΣz (dx) est parfois notée δξ(x)−z (dx) dans la littérature.
En utilisant le formule de la co-aire, on voit que l’on peut écrire une moyenne canonique
sous la forme (c’est une formule de conditionnement) :
Z Z Z
φdµ = φdµΣz exp(−βA(z)) dz (9)
R Σz
où la mesure de probabilité
exp(−βV )|∇ξ|−1 dσΣz exp(−βV )|∇ξ|−1 dσΣz
dµΣz = = (10)
ZΣz exp(−βA(z))
est la mesure canonique µ conditionnée à une valeur fixée de la coordonnée de réaction.
Noter que ZΣz = exp(−βA(z)) est la fonction de partition associée à la mesure condition-
née µΣz .
On voit donc que le problème initial (calcul d’une moyenne canonique) est ramené aux
problèmes suivants :
(i) Calculer l’énergie libre A (à une constante additive près qui n’intervient plus après
normalisation, cf. (9)–(10)) ;
(ii) Echantillonner la mesure µΣz . Si ξ a été bien choisi, on peut espérer que cette tâche
soit nettement plus simple que l’échantillonnage de la mesure µ (penser à nouveau
au cas de la figure 1 pour lequel, dans les tranches Σz , le potentiel est convexe, donc
sans métastabilité).
Dans plusieurs applications, la fonctionnelle φ d’intérêt est en fait une fonction de ξ,
auquel cas le seul problème restant est le calcul de l’énergie libre (i). Par ailleurs, le calcul
de l’énergie libre présente un intérêt en soi, car les physiciens utilisent le profil d’énergie
libre pour guider leur intuition sur les mécanismes réactionnels (y compris sur des aspects
dynamiques).
Remarque 4 Il y a ici une question sous-jacente très importante pour la pratique : peut-
on donner un sens dynamique à l’énergie libre A ? De manière plus précise, est-il possible
d’écrire une dynamique effective sur la coordonnée de réaction ξ(Xt ) (où Xt satisfait la
dynamique (3)) en fonction de A ? Cette question est l’objet d’un travail en cours [25].
Dans la suite, on s’intéresse à trois grandes classes de méthodes pour le calcul de
l’énergie libre A (problème (i)) : l’intégration thermodynamique (section 1), les méthodes
hors équilibre (section 2), et les méthodes adaptatives (sections 3, 4 et 5). Pour une revue,
on renvoie à [26] et au livre [7]. On discutera dans la section 1 l’échantillonnage d’une
mesure sur une sous-variété (problème (ii)).

1 Intégration thermodynamique
L’intégration thermodynamique pour le calcul de l’énergie libre [23] est basée sur deux
ingrédients.
Le premier ingrédient est que la dérivée de A par rapport à la coordonnée de réaction
(appelée la force moyenne) peut être obtenue en échantillonnant la probabilité condition-
nelle µΣz . En effet, en utilisant la formule de la co-aire, on montre que
Z   
0 −1 ∇V · ∇ξ −1 ∇ξ
A (z) = ZΣz 2
− β div 2
exp(−βV )|∇ξ|−1 dσΣz (11)
Σz |∇ξ| |∇ξ|
∇ξ 
Z 
−1 −1
= ZΣz 2
· ∇ Ṽ + β H exp(−β Ṽ ) dσΣz
Σz |∇ξ|
Z
= F dµΣz ,
Σz

4
où
Ṽ = V + β −1 ln |∇ξ|, (12)
la force locale F est définie par
 
∇V · ∇ξ ∇ξ
F = − β −1 div (13)
|∇ξ|2 |∇ξ|2
et le vecteur courbure moyenne H est :
 
∇ξ ∇ξ
H = −∇ · . (14)
|∇ξ| |∇ξ|
Le deuxième ingrédient est que la dynamique contrainte :
 p
dXt = −∇Ṽ (Xt ) dt + 2β −1 dWt + ∇ξ(Xt )dΛt ,
(15)
dΛt tel que ξ(Xt ) = z,

admet pour mesure invariante (et, sous de bonnes hypothèses, est ergodique par rapport à)
la probabilité conditionnelle dµΣz = ZΣ−1 z
exp(−β Ṽ )dσΣz . Noter qu’il faut bien utiliser le
potentiel modifié Ṽ (défini par (12)) et non pas le potentiel original V . Plus précisément,
l’équation différentielle stochastique contrainte (15) s’écrit :
 p 
dXt = P (Xt ) −∇Ṽ (Xt ) dt + 2β dWt + β −1 H(Xt ) dt,
−1 (16)

où l’opérateur de projection orthogonale P (x) est :

P (x) = Id − n(x) ⊗ n(x), (17)

∇ξ
où n(x) = |∇ξ| (x) est un vecteur normal à la sous-variété Σξ(x) au point x. Cette dynamique
peut aussi s’écrire en utilisant le produit de Stratonovitch :
p
dXt = −P (Xt )∇Ṽ (Xt ) dt + 2β −1 P (Xt ) ◦ dWt . (18)

La formulation (15) avec un multiplicateur de Lagrange dΛt associé à la contrainte ξ(Xt ) =

z présente l’avantage de suggérer des schémas de discrétisation (cf. le schéma (21) ci-
dessous).
De plus, on vérifie que
dΛt = dΛm f
t + dΛt , (19)
−1 ∇ξ (X ) · dW et
p
avec dΛmt = − 2β |∇ξ|2 t t

∇ξ  −1

dΛft = · ∇ Ṽ + β H (Xt ) dt = F (Xt ) dt (20)
|∇ξ|2
si bien que Z T Z T
0 1 1
A (z) = lim dΛt = lim dΛft .
T →∞ T 0 T →∞ T 0
La valeur de la force moyenne en un point z et donc en pratique obtenue par une moyenne
sur une trajectoire du multiplicateur de Lagrange dΛt associé à la contrainte ξ(x) = z pour
la dynamique contrainte (15) (cf. l’équation (22) ci-dessous).
Le profil d’énergie libre est ensuite obtenu par une intégration en dimension 1 :
Z z K
X
0
A(z) − A(0) = A (z) dz ' ωi A0 (zi ),
0 i=0

5
où zi désigne des points de Gauss d’une formule d’intégration, et ωi les poids associés. Des
deux problèmes (i) et (ii) mentionnés en introduction, on se ramène donc simplement au
problème (ii) : l’échantillonnage d’une mesure sur une sous-variété. Remarquer que le fait
que la partie à variation finie du processus Λt s’exprime en fonction de la force moyenne
locale F (cf. (20)) permet de calculer A0 (z) sans calculer les dérivées compliquées dans
l’expression (13) définissant F .
Dans [9], nous démontrons rigoureusement l’ensemble de ces résultats, qui étaient pour
la plupart connus des physiciens [32, 34, 13, 8]. Nous proposons également une méthode
de réduction de variance pour le schéma numérique associé. Précisons cela. Un schéma de
discrétisation naturel pour (15) est le suivant :
 p
Xn+1 = Xn − ∇Ṽ (Xn )δt + 2β −1 ∆Wn + λn ∇ξ(Xn ),
(21)
où λn ∈ R est tel que ξ(Xn+1 ) = z,
où ∆Wn = W(n+1)δt − Wnδt désigne un incrément brownien. Il s’agit d’un schéma de type
prédicteur correcteur, où l’on fait d’abord un pas du schéma d’Euler sur l’équation diffé-
rentielle stochastique non contrainte (3), puis un pas de projection sur la sous-variété Σz .
On montre que :
T /δt
1 X
lim lim λn = A0 (z). (22)
T →∞ δt→0 T
n=1
PT /δt
On peut réduire grandement la variance de l’estimateur T1 n=1 λn en remarquant que,
conformément à la formule (19), cet estimateur contient une partie fluctuante de moyenne
nulle correspondant à la discrétisation de la martingale Λm t . Une méthode de réduction de
1 PT /δt λn +λ̃n
variance consiste alors à utiliser l’estimateur T n=1 2 , où λ̃n est défini par (remar-
quer le changement de signe devant l’incrément brownien par rapport à (21)) :
 p
X̃n+1 = Xn − ∇Ṽ (Xn )δt − 2β −1 ∆Wn + λ̃n ∇ξ(Xn ),
où λ̃n ∈ R est tel que ξ(X̃n+1 ) = z.

Remarque 5 Il y a bien sûr plusieurs manières de projeter la dynamique (3) sur la sous-
variété Σz , et on a choisi une méthode de projection qui se discrétise de manière très
naturelle. Il est important de remarquer que la statistique (i.e. la mesure invariante) de la
dynamique projetée dépend de la méthode de projection. Par exemple, une projection uti-
lisant le potentiel de pénalisation Vk (x, z) = V (x) + k(z − ξ(x))2 (dans la dynamique (3))
permet d’échantillonner la mesure invariante µΣz , dans la limite k → ∞, sans modification
du potentiel V en Ṽ . Dans le cas d’une équation différentielle ordinaire (sans mouvement
brownien), les deux approches (projection "rigide" (16) et pénalisation) sont équivalentes.
Ceci est relié à une discussion dans la littérature physique sur les différences statistiques
entre l’utilisation d’une barre rigide et l’utilisation d’un ressort infiniment raide pour im-
poser une contrainte (cf. [37, 19]).

Remarque 6 On pourrait penser à utiliser un algorithme de projection différent dans

l’algorithme (21). Par exemple, on pourrait considérer (noter le changement de ∇ξ(Xn ) à
∇ξ(Xn+1 ) par rapport à (21))
 p
Xn+1 = Xn − ∇Ṽ (Xn )δt + 2β −1 ∆Wn + λn ∇ξ(Xn+1 ),
(23)
où λn ∈ R est tel que ξ(Xn+1 ) = z,
qui est l’équation d’Euler Lagrange associée à la discrétisation :
 p
X̃n+1 = Xn − ∇Ṽ (Xn )δt + 2β −1 ∆Wn ,
Xn+1 = arg minY ∈Σz kX̃n+1 − Y k2 .

6
On montre dans [9, 24] que cette méthode de discrétisation, de même que la discrétisa-
tion (21) sont consistantes (dans la limite δt → 0) avec la dynamique projetée (15). Ceci
dit, nous montrons dans [24], qu’une discrétisation du type (23) n’est pas nécessairement
consistante avec la dynamique projetée associée pour une équation différentielle stochas-
tique générale (avec un terme de diffusion devant le mouvement brownien qui n’est pas une
constante fois l’identité). Autrement dit, les opérations de projection de la dynamique et de
discrétisation ne commutent pas en général. On renvoie à [24] pour plus de détails.

Nous terminons cette section en évoquant un résultat de l’article [15]. Une question
naturelle concernant le schéma de discrétisation (21) est de savoir quelle est l’erreur induite
par la discrétisation en temps sur la mesure invariante de la dynamique. Autrement dit,
quel est l’écart entre la mesure µΣz invariante par la dynamique (15) et la mesure µδt Σz
invariante par la dynamique (21), en fonction de δt. Il s’agit donc d’étendre des résultats
bien connus de vitesse de convergence faible pour la discrétisation par le schéma d’Euler des
équations différentielles stochastiques [36, 35], à des équations différentielles stochastiques
projetées. Dans [15], nous montrons, dans un cadre plus général, que l’erreur faible est
d’ordre 1 : il existe une constante positive δt0 telles que, pour tout pas de temps δt < δt0 ,
pour toute fonction régulière g : Σz → R, il existe C > 0 tel que
Z Z
δt
gdµΣz − gdµΣz ≤ Cδt.
Σz Σz

La preuve s’applique à une équation différentielle stochastique générale ergodique pour une
mesure µΣ de support une sous-variété Σ de Rn , et pour tout schéma numérique consistant
tel que le processus discrétisé est ergodique par rapport à une mesure µδt
Σ portée par Σ. La
démonstration repose sur l’utilisation d’un système de coordonnées orthogonales adapté à
la sous-variété Σ. On renvoie à [15] pour des résultats détaillés appliqués à divers schémas
numériques.

2 Méthodes hors équilibre

Les méthodes hors équilibre pour calculer l’énergie libre ont été introduites récem-
ment [21]. Le principe est de forcer la coordonnée de réaction à bouger selon une évolution
déterministe z : [0, T ] 7→ R prescrite a priori pour évaluer la différence d’énergie libre
(A(z(t)) − A(z(0)))0≤t≤T . On parle de dynamique hors équilibre, car on force la coordon-
née de réaction à avancer en temps (en imposant la contrainte ξ(Xt ) = z(t)) à une vitesse
bien supérieure à celle à laquelle elle évoluerait pour la dynamique de base (3). La question
est alors la suivante : comment établir des relations entre des quantités thermodynamiques
à l’équilibre (comme la différence d’énergie libre) et des quantités calculées hors équilibre ?
Dans la littérature de physique, de tels liens ont été découverts par C. Jarzynski dans [21].
D’une point de vue mathématique, il s’agit de formules de Feynman-Kac. Détaillons cela.
Sur le modèle de (15), on introduit la dynamique hors équilibre suivante :
 p
dXt = −∇Ṽ (Xt ) dt + 2β −1 dWt + ∇ξ(Xt )dΛt ,
(24)
dΛt tel que ξ(Xt ) = z(t).

Plus précisément, l’équation différentielle stochastique (24) s’écrit (comparer avec (18)) :
 p
ext
 dXt = −P (Xt )∇Ṽ (Xt ) dt + 2β −1 P (Xt ) ◦ dWt + ∇ξ(Xt )dΛt ,

z 0 (t) (25)
 dΛext
 t = dt,
|∇ξ(Xt )|2

7
où ◦ désigne un produit de Stratonovitch. Comme dans le cas de la dynamique (15), on
vérifie que
dΛt = dΛm f ext
t + dΛt + dΛt ,
−1 ∇ξ (X ) · dW , dΛf = F (X ) dt et dΛext défini dans (25).
p
avec dΛm t = − 2β |∇ξ|2 t t t t t
En utilisant une formule de Feynman-Kac, on montre dans [29] l’égalité suivante :
A(z(t)) − A(z(0)) = −β −1 ln (E (exp(−βW(t)))) ,
où W(t) est un poids (un travail ou une énergie) associé à la trajectoire (Xs )0≤s≤t :
Z t Z t
0
W(t) = F (Xs )z (s) ds = z 0 (s)dΛfs .
0 0

A partir de cette formule, nous proposons un estimateur de A(z(t)) − A(z(0)) de la forme :

M
!
−1 1 X 
m

A(z(t)) − A(z(0)) ' −β ln exp −βWT /δt , (26)
M
m=1

où (Wnm )0≤n≤T /δtest le poids d’une m-ième réplique d’une trajectoire (discrétisée) de (24)
(et 1 ≤ m ≤ M , M désignant le nombre total de répliques). Plus précisément, pour m fixé,
(Wnm )0≤n≤T /δt est défini par W0m = 0 et, ∀n ∈ {0, . . . , T /δt − 1}

m z(tn+1 ) − z(tn ) f,m

Wn+1 = Wnm + λn
tn+1 − tn
avec tn = nδt,
z(tn+1 ) − z(tn ) p −1 ∇ξ
λf,m m
n = λn − + 2β (X m ) · ∆Wnm ,
|∇ξ(Xnm )|2 |∇ξ|2 n
et (sur le modèle du schéma de discrétisation (21))
 m p
Xn+1 = Xnm − ∇V (Xnm )δt + 2β −1 ∆Wnm + λm m
n ∇ξ(Xn ),
m m
où λn ∈ R est tel que ξ(Xn+1 ) = z(tn+1 ).
Dans ces équations, ∆Wnm = W(n+1)δt
m −Wnδtm , où (W m )
t t≥0,1≤m≤M désigne M mouvements
Browniens indépendants. On renvoie à [29] pour des illustrations numériques de cette
approche.
On peut signaler qu’en pratique, l’efficacité de ces méthodes hors équilibre pour le
calcul de l’énergie libre n’a pas été vraiment démontrée. Il apparaît en effet que la va-
riance de l’estimateur (26) explose si la variation en temps de l’application (z(t))0≤t≤T est
trop rapide. On est donc amené à utiliser des variations assez lente, ce qui se rapproche
alors d’une méthode d’intégration thermodynamique (cf. section 1). Mentionnons tout de
même le travail [33], où les auteurs proposent un algorithme comprenant une étape de
reconfiguration stochastique pour mieux contrôler la variance de l’estimateur.

3 Méthodes adaptatives : algorithmes

Nous abordons dans cette section les méthodes qui ont été proposées le plus récemment
pour calculer les différences d’énergie libre. Il s’agit des méthodes adaptatives, du type
Wang-Landau [38], Adaptive Biasing Force [10, 18], ou metadynamics [20], par exemple.
Dans [30], nous proposons un cadre unifié pour toutes ces méthodes, que nous résumons
maintenant.
Le principe des méthodes adaptatives est le suivant. Il s’agit de modifier le potentiel
vu par le processus Xt en fonction des configurations visitées, avec deux objectifs :

8
 (i) explorer de manière efficace l’espace des configurations,
(ii) calculer le profil d’énergie libre.
Le potentiel modifié dépend du temps et est de la forme, pour tout x ∈ Rn ,

Vt (x) = V (x) − At (ξ(x)),

où At est une approximation de l’énergie libre à l’instant t, au vu des configurations visitées

par le système jusqu’à cet instant. On cherche à faire évoluer At de manière à ce que
limt→∞ At = A (cf. objectif (ii) ci dessus).
Pour comprendre comment mettre à jour At , faisons l’hypothèse un instant que le
processus au temps t est immédiatement à l’équilibre par rapport au potentiel Vt , i.e., que

ψ(t, ·) = ψ eq (t, ·) ∝ exp(−βVt )

où ψ désigne la densité de probabilité de la variable aléatoire Xt . On peut alors calculer

une énergie libre observée, en remplaçant V par Vt dans la formule (7) :
Z 
−1 −1
Aobs
t (z) = −β ln exp(−βV t )|∇ξ| dσΣz
Σz

et une force moyenne observée, en remplaçant V par Vt dans la formule (11) :

Z   
∇Vt · ∇ξ −1 ∇ξ
− β div exp(−βVt )|∇ξ|−1 dσΣz
obs 0 Σz |∇ξ|2 |∇ξ|2
(At ) (z) = Z ,
exp(−βVt )|∇ξ|−1 dσΣz
Σz
Z
F exp(−βVt )|∇ξ|−1 dσΣz
= ZΣz − A0t (z),
−1
exp(−βVt )|∇ξ| dσΣz
Σz

où F désigne la force moyenne locale (13). Ici et dans toute la suite, on note par 0 une
dérivée par rapport à la variable z (la valeur de la coordonnée de réaction). On vérifie
facilement que
0 0 0
Aobs
t = A − At et (Aobs
t ) = A − At . (27)
Au vu de (27) et pour que limt→∞ At = A, il semble naturel de proposer les deux types de
mise à jour suivantes (où l’on ne suppose plus cette fois le processus à l’équilibre : ψ(t, ·)
désigne la densité de probabilité de Xt ) :
– Mise à jour du potentiel At (méthode de type Adaptive Biasing Potential (ABP))
Z
∂At 1 −1
(z) = − β ln ψ(t, ·)|∇ξ|−1 dσΣz , (28)
∂t τ
où τ désigne un paramètre numérique homogène à un temps ;
– Ou bien mise à jour de la force moyenne A0t (méthode de type Adaptive Biasing Force
(ABF))
∂A0t F ψ(t, ·)|∇ξ|−1 dσΣz
R 
1 0
(z) = R − At (z) . (29)
∂t τ ψ(t, ·)|∇ξ|−1 dσΣz
En effet, si on suppose dans (28) que ψ = ψ eq , on obtient l’évolution suivante pour At :
∂At 1
∂t = τ (A−At ) et donc At converge bien vers A. Un raisonnement similaire tient pour (29).
Si le processus était à l’équilibre à chaque instant (ψ = ψ eq ), on obtiendrait donc l’énergie

9
libre (ou sa dérivée, la force moyenne) en temps long. On espère que cela sera encore le cas
sans cette hypothèse d’équilibre instantané. C’est l’objet de la section 4.
Il est important de noter que les deux approches (ABP et ABF) sont bien différentes,
au sens où la dérivée de l’équation (28) par rapport à z n’est pas l’équation (29) (ce serait
le cas si ψ = ψ eq ).
Remarquer également que la méthode est construite de telle sorte que si un état sta-
tionnaire est atteint, i.e. si (ψ(t, ·), At ) converge en temps long vers un couple (ψ∞ , A∞ ),
alors, nécessairement, A∞ = A (à une constante additive près). Il suffit pour s’en rendre
compte d’écrire les formules (28) et (29) avec (ψ(t, ·), At ) = (ψ∞ , A∞ ).
Pour le moment, nous n’avons pas précisé la dynamique sous-jacente pour Xt : notre
discussion est valable aussi bien pour une dynamique de Langevin (5), que pour une dyna-
mique de gradient (3). Pour fixer les idées, écrivons le sytème complet pour une dynamique
adaptative de type ABF (29) (en prenant formellement τ = ∞ : mise à jour instantanée
de A0t ), construite à partir d’une dynamique de gradient (3) :
( p
dXt = −∇(V − At ◦ ξ)(Xt ) dt + 2β −1 dWt ,
(30)
A0t (z) = E (F (Xt )|ξ(Xt ) = z) ,

où ◦ désigne ici l’opérateur de composition et E (F (Xt )|ξ(Xt ) = z) désigne une espérance

conditionnelle. En terme de la densité ψ(t, ·) de Xt , on a l’équation de Fokker-Planck
associée :
∇(V − At ◦ ξ)ψ + β −1 ∇ψ ,


 ∂t ψ = div

−1
R
0 Σz F ψ(t, ·)|∇ξ| dσΣz (31)
A
 t
 (z) = R
−1
.
Σz ψ(t, ·)|∇ξ| dσΣz

On vérifie que les équations (30) et (31) sont bien équivalentes, car la mesure conditionnée
ψ|∇ξ|−1 dσ
de la loi de probabilité ψ(x) dx sachant que ξ(x) = z est donnée par R ψ|∇ξ|−1 dσΣz (ceci
Σz
découle de la formule de la co-aire (8)).
Cette dynamique est connue des physiciens [10, 18] sous le nom de Adaptive Biasing
Force method, sous une forme "proche" où on approxime l’espérance conditionnelle par une
moyenne ergodique (cf. la section suivante).
L’intérêt du cadre que nous proposons pour les algorithmes adaptatifs est qu’il permet
de comprendre un grand nombre de méthodes adaptatives dans un formalisme similaire, et
d’imaginer de nouvelles méthodes, notamment en jouant sur la manière de mettre à jour
At et A0t (on a ici donné la méthode la plus simple, pour simplifier l’exposé). On renvoie
à [30] pour plus de détails. De plus, ce cadre permet d’étudier la (vitesse de) convergence
de certains algorithmes adaptatifs (cf. section 4) et de proposer de nouvelles méthodes de
discrétisation (cf. section 5).

Remarque 7 D’un point de vue probabiliste, les méthodes adaptatives que nous avons
décrites peuvent être vues comme des méthodes de réduction de variance par fonction d’im-
portance, avec une fonction d’importance adaptative. En effet, en partant du principe que
pour un bon choix de ξ, la dynamique (30) converge plus vite vers l’équilibre que la dyna-
mique (3) (ce que l’on va démontrer dans la section suivante), on dispose d’une méthode
qui permet d’échantillonner des mesures métastables. Cette remarque donne à ces méthodes
une portée assez générique, du moins pour les applications dans lesquelles on peut intuiter
une bonne "direction de métastabilité". On donne cette "direction", et l’algorithme calcule
automatiquement une bonne fonction d’importance "dans cette direction".

10
4 Méthodes adaptatives : étude théorique
L’objectif de cette section théorique est de prouver, pour la dynamique (30), que A0t
converge vers A0 (objectif (ii) de la section précédente), et de déterminer à quelle vitesse
la convergence a lieu (objectif (i) de la section précédente). Autrement dit, on cherche à
déterminer en quoi est-ce que la dynamique (30) converge plus vite vers un état stationnaire
que la dynamique de gradient (3) dont on a dit qu’elle était inutilisable en pratique car
métastable. Cette section s’appuie sur l’article [31].
Pour montrer que la dynamique (30) va plus vite vers un état stationnaire que (3), il faut
décider d’une mesure de la vitesse de convergence. Il y a plusieurs manières de quantifier un
comportement métastable, comme par exemple le temps de sortie des puits de potentiel, la
variance asymptotique des estimateurs considérés, le temps de “décorrélation” du processus
Xt ou la vitesse de convergence de la loi de Xt vers µ. Dans le suite, c’est ce quatrième
critère que l’on choisit.
Commençons par rappeler une analyse classique de la vitesse de convergence vers l’équi-
libre pour (3). On rappelle qu’associée à la dynamique (3), on peut écrire une équation de
Fokker-Planck (linéaire) sur la densité φ de Xt solution de (3) :

∂t φ = div(∇V φ + β −1 ∇φ). (32)

Une approche classique pour étudier le comportement en temps long de cette équation
consiste à introduire l’entropie de φ par rapport à la solution stationnaire φ∞ = Z −1 exp(−βV )
de (32) : Z  
φ
H(φ|φ∞ ) = ln φ.
Rd φ∞
La fonction φ∞ est la densité de la mesure µ. Remarquer que H ≥ 0 et que H = 0 si et
seulement si φ = φ∞ . On rappelle l’inégalité (de Csiszar-Kullback)
Z p
|φ − φ∞ | ≤ 2H(φ|φ∞ ). (33)
Rd

Contrôler l’entropie permet donc de contrôler la distance (en norme L1 ) entre φ et φ∞ . En

utilisant le fait que φ∞ est une solution stationnaire, on vérifie que
Z  2
dH(φ|φ∞ ) φ
= −β −1 ∇ ln φ. (34)
dt φ∞
Le membre de droite est négatif, et l’entropie H décroît donc au cours du temps.
Pour obtenir la convergence de H vers 0 (et donc de φ vers φ∞ par (33)) on a ensuite
besoin d’une inégalité de Sobolev logarithmique, dont il faut vérifier si elle est satisfaite
pour la fonction φ∞ que l’on considère : il existe une constante R > 0 telle que pour toute
fonction ψ positive et d’intégrale 1,
Z   Z  2
ψ 1 ψ
ln ψ≤ ∇ ln ψ. (35)
φ∞ 2R φ∞
En utilisant (35) dans (34), on en déduit la convergence exponentielle de H vers 0 :

H(φ(t, ·)|φ∞ ) ≤ H(φ(0, ·)|φ∞ ) exp −2β −1 Rt .

Il nous reste donc à expliquer pourquoi φ∞ satisfait génériquement l’inégalité fonction-

nelle (35) dans les applications qui nous intéressent ici.
Avant de rappeler quelques critères bien connus pour établir (35) on donne une défini-
tion précise de l’inégalité de Sobolev logarithmique, pour de futures références.

11
Définition 8 On dit que la mesure ν satisfait une inégalité de Sobolev logarithmique de
constante ρ > 0 (noté par la suite ISL(ρ)) si et seulement si, pour toute mesure µ telle que
µ est absolument continue par rapport à ν (noté µ  ν dans la suite)
1
H(µ|ν) ≤ I(µ|ν)
2ρ
Z  
dµ
où H(µ|ν) = ln dµ est l’entropie de la mesure µ par rapport à la mesure ν et
dν
Z  2
dµ
I(µ|ν) = ∇ ln dµ l’information de Fisher de la mesure µ par rapport à la
dν
mesure ν.

Avec un léger abus de notation, on dit qu’une densité de probabilité satisfait une inégalité de
Sobolev logarithmique si la mesure associée satisfait une inégalité de Sobolev logarithmique
et pour deux densités de probabilité ψ et φ, on note l’entropie des mesures associées H(ψ|φ).
On rappelle le critère dit de Bakry-Emery :

Proposition 9 Soit Π une fonction α-convexe, au sens où pour tout vecteur X, X T ∇2 ΠX ≥

α|X|2 (où ∇2 Π désigne la Hessienne de Π). Alors la densité de probabilité ψ∞ ∝ exp(−Π)
satisfait une inégalité de Sobolev logarithimique pour une constante ρ ≥ α.

On rappelle également le critère perturbatif de Holley-Strook :

Proposition 10 Soit Π une fonction telle que la densité de probabilité ψ∞ ∝ exp(−Π)

satisfait une inégalité de Sobolev logarithmique de constante ρ. Soit une fonction Π̃ bornée
et la densité de probabilité ψ̃∞ ∝ exp(−Π + Π̃). Alors ψ̃∞ satisfait également une inégalité
de Sobolev logarithmique de constante ρ̃ ≥ ρ exp(−osc Π̃) où osc Π̃ = sup Π̃ − inf Π̃.

Donc, dès que le potentiel V est une perturbation bornée d’un potentiel α-convexe, φ∞
satisfait une inégalité de Sobolev logarithmique.
En résumé, si φ∞ satisfait ISL(R), on a :

H(φ(t, ·), φ∞ ) ≤ H(φ(0, ·), φ∞ ) exp(−2β −1 Rt) (36)

R
Par l’inégalité de Csiszar-Kullback (33), on en déduit la convergence vers 0 de |φ(t, ·)−φ∞ |
à vitesse exponentielle, avec un taux β −1 R. On considère donc dans la suite que (3) à
un comportement métastable si et seulement si R est proche de 0. Cette définition est
raisonnable au sens où : (i) on vérifie que R est effectivement très proche de 0 pour des
potentiels comportant des puits séparés par des barrières très hautes (ce qui correspond
à l’image intuitive d’un potentiel métastable) ; (ii) la "vraie" vitesse de convergence de
la dynamique (3) serait donnée par le trou spectral de l’opérateur div(∇V · +β −1 ∇·), et
on peut montrer que la dégénerescence du trou spectral pour cet opérateur implique la
dégénérescence de l’inégalité de Sobolev logarithmique pour la mesure µ (cf. [11, 6]).
Intéressons-nous maintenant à la dynamique (30). Nous allons appliquer une méthode
d’entropie pour analyser la vitesse de convergence de cette dynamique vers son état sta-
tionnaire. La difficulté par rapport à l’étude de (3) et (32) est que l’équation de Fokker-
Planck (31) est non-linéaire. Pour simplifier, nous supposons dans toute la suite que ξ est
à valeurs dans le tore :
ξ : Rn → T.
Les résultats se généralisent au cas où ξ est à valeurs dans un domaine non borné en
ajoutant un potentiel confinant adéquat (cf. [31], un tel potentiel étant d’ailleurs utilisé
en pratique pour les calculs d’énergie libre par de telles méthodes). De plus, nous allons

12
légèrement modifier la dynamique (30) pour simplifier l’étude. Nous considérons dans la
suite de cette section la dynamique :
( p
dXt = −∇ V − At ◦ ξ − β −1 ln(|∇ξ|−2 ) (Xt )|∇ξ|−2 (Xt ) dt + 2β −1 |∇ξ|−1 (Xt )dWt ,

A0t (z) = E (F (Xt )|ξ(Xt ) = z) ,

(37)
associée à l’équation de Fokker-Planck :
|∇ξ|−2 ∇(V − At ◦ ξ)ψ + β −1 ∇ψ ,



 ∂t ψ = div

−1
R
0 Σz F ψ(t, ·)|∇ξ| dσΣz (38)
 At (z) =
 R
−1
.
Σz ψ(t, ·)|∇ξ| dσΣz

Remarquer que dans le cas où |∇ξ| = 1 (i.e. ξ est la distance signée à Σ0 ) la dynamique (37)
est exactement (30). Dans le cas général, nous avons besoin de rajouter les termes fonctions
de |∇ξ| pour simplifier l’analyse. On renvoie à [31] pour des résultats similaires sur la
convergence de la dynamique originale (30).
Nous devons introduire quelques notations. On note ψ ξ la "densité marginale de ψ
selon ξ" : Z
ψ ξ (t, z) = ψ(t, ·)|∇ξ|−1 dσΣz . (39)
Σz

On introduit également µξt (dx|z) la mesure conditionnelle de ψ(t, x) dx sachant que ξ(x) =
z:
ψ(t, ·)|∇ξ|−1 dσΣz
dµξt (dx|z) = . (40)
ψ ξ (t, z)
La formulation "densité marginale" est justifiée par le fait que, par la formule de la co-
aire (8), on vérifie que l’image de la mesure ψ(t, x) dx par la fonction ξ est la mesure
ψ ξ (t, z) dz. Dans le cas ξ(x1 , . . . , xn ) = x1 par exemple, on a simplement (par Fubini)
ψ ξ (t, z) = ψ(t, z, x2 , . . . , xn ) dx2 . . . dxn et dµξt (dx|z) = ψ(t,z,x2 ,...,x n ) dx2 ...dxn
R
ψ ξ (t,z)
.
Une remarque fondamentale est la suivante :
Lemme 11 Soit ψ solution de (38). Alors la marginale ψ ξ de ψ selon ξ définie par (39)
satisfait l’équation de la chaleur :

∂t ψ ξ − β −1 ∂z,z ψ ξ = 0 sur T. (41)

Il est très facile de vérifier ce lemme dans le cas où ξ(x1 , . . . , xn ) = x1 , en intègrant l’équa-
tion (38) par rapport à x2 , . . . , xn . Dans le cas général, la preuve est basée sur la formule
de la co-aire (8). Ce premier résultat montre déjà qualitativement que la dynamique (37)
va plus vite vers son état stationnaire que la dynamique (3), du moins si ξ a été bien choisi.
On voit en effet que l’équation sur ψ ξ est une équation de la chaleur, correspondant à un
potentiel constant : en utilisant la dynamique adaptative, on a "aplati" le potentiel dans
la direction de la coordonnée de réaction. Si cette direction est bien la direction de méta-
stabilité principale pour (3), on a donc très fortement amélioré la vitesse de convergence
en utilisant (37).
Nous allons maintenant quantifier cela. Nous avons besoin de quelques notations sup-
plémentaires. On note ψ∞ la densité stationnaire attendue :

ψ∞ = exp(−β(V − A ◦ ξ)).

A la densité de probabilité ψ∞ , on associe sa marginale selon ξ (de manière consistante

avec (39)) :
ξ
ψ∞ (z) = 1T (z),

13
et la mesure conditionnelle de ψ∞ (x) dx sachant que ξ(x) = z (de manière consistante
avec (40)) :
dµξ∞ (·|z) = ZΣ−1
z
exp(−βV )|∇ξ|−1 dσΣz ,
où on rappelle que (cf. (7))
Z
ZΣz = |∇ξ|−1 exp(−βV )dσΣz = exp(−βA(z)).
Σz

Pour ψ une solution de (38), on introduit l’entropie totale :

E(t) = H(ψ(t, ·)|ψ∞ ),

l’entropie macroscopique :
EM (t) = H(ψ ξ (t, ·)|ψ∞
ξ
),
l’entropie locale à une valeur z de la coordonnée de réaction :
!
ψ(t, ·)|∇ξ|−1 dσΣz
  Z
ψ(t, ·) . ψ∞
em (t, z) = H µξt (·|z) µξ∞ (·|z) = ln ,
Σz
ξ
ψ ξ (t, z) ψ∞ (z) ψ ξ (t, z)

et enfin l’entropie microscopique

Z
Em (t) = em (t, z)ψ ξ (t, z) dz.
M

On vérifie facilement que

E(t) = EM (t) + Em (t).
Noter qu’en utilisant le Lemme (11), on sait déjà que ψ ξ converge vers son état stationnaire
ξ
ψ∞ à vitesse exponentielle. En particulier :

EM (t) ≤ EM (0) exp(−2β −1 rt), (42)

où
r = 4π 2 (43)
est appelé dans la suite taux de convergence macroscopique. Pour étudier la convergence
de ψ vers ψ∞ , il suffit donc d’étudier la convergence de Em vers zéro.
Les dénominations "macroscopique" et "microscopique" font ici référence au fait que
ξ est considéré comme une fonction qui à une configuration associe une grandeur macro-
scopique (variable lente). Les variables "orthogonales à ξ" (dans une feuille Σz ) sont au
contraire des grandeurs microscopiques (variables rapides).
Dans [31], on montre le résultat suivant :

Théorème 12 On suppose :
(H1) ξ : Rn → T est une fonction régulière telle que |∇ξ| > 0.
(H2) V et ξ sont des fonctions régulières telles que k∇ξkL∞ ≤ m < ∞ et k∇Σz F kL∞ ≤
M < ∞ (où ∇Σz = P ∇ (cf. (17)) désigne le gradient surfacique),
(H3) V et ξ sont telles que ∃ρ > 0, ∀z ∈ T, la mesure conditionnelle µξ∞ (·|z) satisfait
une inégalité de Sobolev logarithmique de constante ρ.
Alors, on montre que

14
 1. L’entropie microscopique Em converge à vitesse exponentielle vers zéro : ∃C > 0,
∃λ > 0, ∀t ≥ 0, p
Em (t) ≤ C exp(−λt). (44)
Plus précisément, si ρm−2 6= r, alors

λ = β −1 min(ρm−2 , r). (45)

Dans le cas particulier ρm−2 = r, pour tout λ < β −1 min(ρm−2 , r), il existe C > 0
tel que (44) est vérifié.
√
2. La racine carrée de l’entropie totale E, les quantités kψ(t, ·) − ψ∞ kL1 (Rn ) et kA0t −
A0 k2L2 (T) convergent exponentiellement vite vers zéro au taux λ.

La preuve est basée sur un calcul explicite de dE 0 0

dt et sur une majoration de |At (z) − A (z)|
m

ξ ξ
par l’entropie relative de µt (·|z) par rapport à µ∞ (·|z), basée sur une inégalité de Talagrand.
Ces résultats montrent donc que, essentiellement sous une hypothèse de couplage borné
(H2) (cf. remarque 13 pour une justification de ce terme) et une hypothèse d’ergodicité
des variables microscopiques (H3), A0t converge exponentiellement vite vers A0 au taux
λ = β −1 min(ρm−2 , r), qui doit être comparé à β −1 R dans (36). Le taux de convergence λ
est en fait seulement limité par la constante ρ de l’inégalité de Sobolev logarithmique pour
les mesures conditionnelles µξ∞ (·|z). En effet, on peut rendre en pratique r aussi grand
que l’on veut de deux façons : en restreignant le domaine de calcul de l’énergie libre, ou
ξ
en utilisant une méthode de sélection pour accélérer la convergence de ψ ξ vers ψ∞ (cf.
section 5). La valeur de ρ dépend évidemment du choix de la coordonnée de réaction. Dans
notre langage, une "bonne coordonnée de réaction" est telle que ρ est aussi grand que
possible.
Pour une généralisation de ces résultats (cas où ξ est à valeurs dans Tp ou Rp , cas de la
dynamique ABF non-modifiée (30), discussion d’autres jeux d’hypothèses possibles), nous
renvoyons à [31].

Remarque 13 Il est naturel à ce stade de se poser la question suivante : sous les hypo-
thèses du théorème 12, que peut-on dire de la constante de Sobolev logarithmique R pour
la mesure "totale" µ. Dans [28], nous généralisons un résultat de [17] à une fonction ξ
non-linéaire en montrant le critère suivant (appelé "critère à double-échelle") pour prouver
une inégalité de Sobolev logarithmique pour la mesure µ. Si on suppose : (i) la margi-
nale de µ selon ξ satisfait une ISL(r), (ii) les mesures conditionnelles de µ à ξ(x) = z
fixé satisfont une ISL(ρ) (ce qui est (H3) dans le théorème 12), (iii) V et ξ sont tels que
k∇Σz F kL∞ ≤ M < ∞ et |∇ξ|2 ≥ m > 0 (ce qui est relié à (H2) dans le théorème 12),
alors µ satisfait une ISL de constante R avec
 s 
2
1 M 2m M 2m
R≥ rm + +ρ− rm + + ρ − 4rmρ .
2 ρ ρ

On sait prouver l’égalité pour des lois normales et une fonction ξ linéaire. Noter que le
membre de droite est majoré par min(rm, ρ), qui correspond au cas d’un couplage nul
(M = 0). L’hypothèse (iii) est appelée une hypothèse de couplage borné, car dans le cas
particulier d’une loi normale et d’une fonction ξ linéaire, k∇Σz F kL∞ mesure la covariance
entre les composantes microscopiques et les composantes macroscopiques.
Pour revenir à la question posée au début de cette remarque, on voit que la constante de
Sobolev logarithmique R pour la mesure de Boltzmann-Gibbs µ est typiquement beaucoup
plus petite que λ car la constante de Sobolev logarithmique pour la marginale de µ selon

15
ξ (notée r ci-dessus) est typiquement beaucoup plus petite que la constante de Sobolev
logarithmique pour la marginale de ψ∞ selon ξ (notée r ci-dessus). On peut à nouveau
penser à l’exemple très simple de la Figure 1 pour lequel r et R sont très petits pour une
grande barrière de potentiel, comparativement aux constantes r et ρ.

5 Méthodes adaptatives : discrétisation

Dans cette section, on discute la discrétisation des méthodes adaptatives, en considé-
rant (37) pour fixer les idées. Les schémas de discrétisation proposés sont généralisables à
d’autres dynamiques adaptatives comme celles évoquées dans la section 3, y compris quand
elles sont construites à partir d’une dynamique de Langevin (5).
Dans la communauté "dynamique moléculaire", les physiciens utilisent plus volontiers
des moyennes ergodiques sur de longues trajectoires que des moyennes empiriques sur plu-
sieurs marcheurs pour calculer des moyennes. Ainsi, plutôt qu’une espérance conditionnelle,
une méthode du type (37) utiliserait typiquement une mise à jour :
Rt
0 F (Xs )δ  (ξ(Xs ) − z) ds
At (z) = 0 R t ,

0 δ (ξ(Xs ) − z) ds
où δ  désigne une approximation de l’identité. Dans [30], nous proposons d’utiliser plutôt
des systèmes particulaires, avec plusieurs marcheurs en interaction et des moyennes empi-
riques. Ainsi, une approximation particulaire de (37) s’écrit : pour 1 ≤ n ≤ N , N désignant
le nombre total de marcheurs,
dXtn,N = −∇ V − AN −1
ln(|∇ξ|−2 ) (Xtn,N )|∇ξ|−2 (Xtn,N ) dt



 t ◦ξ−β

+ 2β −1 |∇ξ|−1 (Xtn,N )dWtn ,

 p
(46)
PN m,N  m,N

 N 0 m=1 F (X t )δ (ξ(X t ) − z)
 (At ) (z) = ,


 (ξ(X m,N ) − z)
PN
m=1 δ t
où (Wtn )t≥0,n≥1 désigne une famille de mouvements browniens indépendants. Cette ap-
proche est plus facile à analyser mathématiquement, très facilement parallélisable, et per-
met en outre d’ajouter des mécanismes de sélection pour sélectionner les "bons marcheurs"
et tuer les "mauvais" afin d’accélérer la convergence (cf. ci-dessous). En pratique, on uti-
lise plutôt une discrétisation par histogramme de (AN 0
t ) qu’une approximation de l’identité,
mais la discrétisation (46) reste représentative des méthodes standards, qui comportent ty-
piquement trois paramètres de discrétisation : le pas de temps, le nombre de marcheurs N ,
et un paramètre  donnant la taille caractéristique de discrétisation de la variable coordon-
née de réaction.
Dans [22], nous analysons la vitesse de convergence du système particulaire (46) dans le
cas particulier où ξ(x1 , . . . , xn ) = x1 et l’espace de configuration est Tn ou T × Rn−1 . Nous
montrons notamment que l’erreur comporte un terme C(T,) √
N
, avec une constante C(T, )
qui dépend du temps maximum T et du paramètre . L’expression de C que nous obtenons
est cependant sous-optimale par rapport aux résultats numériques.
Nous terminons cette section en explicitant un exemple de mécanisme de sélection
possible, proposé dans [30]. On modifie la dynamique (38) en considérant : pour une
constante c > 0
∂z,z ψ ξ
  
−2 −1



 ∂t ψ = div |∇ξ| ∇(V − At ◦ ξ)ψ + β ∇ψ + c ◦ ξ ψ,


ψξ

(47)
F ψ(t, .)|∇ξ|−1 dσΣz
R
 0 Σ z
 At (z) = R .


−1
Σz ψ(t, .)|∇ξ| dσΣz

16
On vérifie facilement que cette dynamique modifiée est encore telle que, si un état sta-
tionnaire (ψ∞ , A0∞ ) est atteint, alors A0∞ = A0 . De plus, par un calcul similaire à celui
permettant de prouver le Lemme 11, on montre que la marginale ψ ξ de ψ suivant ξ satis-
fait alors (comparer avec (41)) :

∂t ψ ξ − (β −1 + c)∂z,z ψ ξ = 0 sur T.

L’ajout du terme supplémentaire dans Fokker-Planck permet donc d’améliorer la "vitesse

de convergence macroscopique", qui est celle reliée au paramètre r de la section précédente
(cf. (43)). La diffusion dans l’espace de la coordonnée de réaction est plus rapide.
Le système particulaire associé à (47) consiste à modifier (46) en associant un poids
 Z t
∂z,z ψ ξ

n,N



wn,N (t) = exp c ξ Xs ds
0 ψξ

à la particule (Xtn,N ) et à effectuer une moyenne pondérée par ces poids pour le calcul de
(AN 0
t ) dans (46). En pratique, on effectue de temps en temps une étape de reconfiguration
stochastique, en dupliquant (respectivement tuant) les marcheurs qui ont le poids le plus
élevé (respectivement le plus faible), en gardant un nombre total de marcheurs constant
(cf. [12] pour une présentation de telles méthodes).
∂ ψξ
La quantité z,zψξ
est discrétisée par une méthode d’histogramme pour approximer ψ ξ
et des différences finies :
∂z,z ψ ξ ψ ξ (z + δz) − 2ψ ξ (z) + ψ ξ (z − δz)
(z) ' ,
ψξ ψ ξ (z)δz 2
 ξ
ψ (z + δz) + ψ ξ (z) + ψ ξ (z − δz)

3 ξ
' ξ − ψ (z) .
ψ (z)δz 2 3

Ainsi, la méthode donne plus de poids aux particules qui sont dans des régions moins
explorées que les régions voisines.
On renvoie à [30] pour des expériences numériques montrant l’intérêt de cette mé-
thode sur un exemple concret. En pratique, on observe qu’il existe un optimum pour le
paramètre c : si c devient trop grand, la variance des résultats augmente ce qui dégrade
le résultat (noter qu’augmenter c revient à augmenter la disparité des poids et donc à
diminuer le nombre effectif de marcheurs contribuant à la moyenne donnant (AN 0
t ) (z)).

Références
[1] M.P. Allen and D.J. Tildesley. Computer simulation of liquids. Oxford Science Publi-
cations, 1987.
[2] L. Ambrosio, N. Fusco, and D. Pallara. Functions of bounded variation and free
discontinuity problems. Oxford science publications, 2000.
[3] C. Ané, S. Blachère, D. Chafaï, P. Fougères, I. Gentil, F. Malrieu, C. Roberto, and
G. Scheffer. Sur les inégalités de Sobolev logarithmiques. Société Mathématique de
France, 2000. In French.
[4] A. Arnold, P. Markowich, G. Toscani, and A. Unterreiter. On convex Sobolev inequa-
lities and the rate of convergence to equilibrium for Fokker-Planck type equations.
Comm. Part. Diff. Eq., 26 :43–100, 2001.
[5] V. Babin, C. Roland, and C. Sagui. Adaptively biased molecular dynamics for free
energy calculations. J. Chem. Phys., 128 :134101, 2008.

17
 [6] D. Bakry. On Sobolev and logarithmic Sobolev inequalities for Markov semigroups.
In New trends in stochastic analysis, pages 43–75. World Scientific, 1997.
[7] C. Chipot and A. Pohorille, editors. Free Energy Calculations, volume 86 of Springer
Series in Chemical Physics. Springer, 2007.
[8] G. Ciccotti, R. Kapral, and E. Vanden-Eijnden. Blue moon sampling, vectorial reac-
tion coordinates, and unbiased constrained dynamics. ChemPhysChem, 6(9) :1809–
1814, 2005.
[9] G. Ciccotti, T. Lelièvre, and E. Vanden-Eijnden. Projection of diffusions on sub-
manifolds : Application to mean force computation. Commun. Pur. Appl. Math.,
61(3) :371–408, 2008.
[10] E. Darve and A. Porohille. Calculating free energy using average forces. J. Chem.
Phys., 115 :9169–9183, 2001.
[11] E.B. Davies. Spectral properties of metastable Markov semigroups. J. Funct. Anal.,
52 :315–329, 1983.
[12] P. Del Moral. Feynman-Kac Formulae : Genealogical and Interacting Particle Systems
with Applications. Springer-Verlag, 2004.
[13] W. E and E. Vanden-Eijnden. Metastability, conformation dynamics, and transition
pathways in complex systems. In Multiscale modelling and simulation, volume 39 of
Lect. Notes Comput. Sci. Eng., pages 35–68. Springer, Berlin, 2004.
[14] L.C. Evans and R.F. Gariepy. Measure theory and fine properties of functions. Studies
in Advanced Mathematics. CRC Press, 1992.
[15] E. Faou and T. Lelièvre. Conservative stochastic differential equations : Mathematical
and numerical analysis. Mathematics of computation, 2008. to appear.
[16] D. Frenkel and B. Smit. Understanding molecular simulation : from algorithms to
applications. Academic Press, 2002.
[17] N. Grunewald, F. Otto, C. Villani, and M.G. Westdickenberg. A two-scale approach to
logarithmic Sobolev inequalities and the hydrodynamic limit. Ann. Inst. H. Poincaré
Probab. Statist., 2008. To appear.
[18] J. Hénin and C. Chipot. Overcoming free energy barriers using unconstrained mole-
cular dynamics simulations. J. Chem. Phys., 121 :2904–2914, 2004.
[19] E.J. Hinch. Brownian motion with stiff bonds and rigid constraints. J. Non-Newtonian
Fluid Mech., 271 :219–234, 1994.
[20] M. Iannuzzi, A. Laio, and M. Parrinello. Efficient exploration of reactive poten-
tial energy surfaces using Car-Parrinello molecular dynamics. Phys. Rev. Lett.,
90(23) :238302, 2003.
[21] C. Jarzynski. Equilibrium free energy differences from nonequilibrium measurements :
A master equation approach. Phys. Rev. E, 56(5) :5018–5035, 1997.
[22] B. Jourdain, T. Lelièvre, and R. Roux. Existence, uniqueness and convergence of a
particle approximation for the adaptive biasing force process. in preparation.
[23] J.G. Kirkwood. Statistical mechanics of fluid mixtures. J. Chem. Phys., 3 :300–313,
1935.
[24] C. Le Bris, T. Lelièvre, and E. Vanden Eijnden. Analysis of some discretization
schemes for constrained stochastic differential equations. C.R. Acad. Sci. Paris, Ser.
I., 346(7-8) :471–476, 2008.
[25] F. Legoll and T. Lelièvre. Effective dynamics for the reaction coordinate. In prepara-
tion.

18
[26] F. Legoll, T. Lelièvre, and G. Stoltz. Some remarks on sampling methods in molecular
dynamics. ESAIM Proceedings, 22 :217–233, 2008.
[27] B. Leimkuhler and S. Reich. Simulating Hamiltonian dynamics. Cambridge University
Press, 2004.
[28] T. Lelièvre. A general two-scale criteria for logarithmic Sobolev inequalities. Journal
of Functional Analysis, 2008. to appear.
[29] T. Lelièvre, M. Rousset, and G. Stoltz. Computation of free energy differences through
nonequilibrium stochastic dynamics : The reaction coordinate case. J. Comput. Phys.,
222(2) :624–643, 2007.
[30] T. Lelièvre, M. Rousset, and G. Stoltz. Computation of free energy profiles with
adaptive parallel dynamics. J. Chem. Phys., 126 :134111, 2007.
[31] T. Lelièvre, M. Rousset, and G. Stoltz. Long-time convergence of an adaptive biasing
force method. Nonlinearity, 21 :1155–1181, 2008.
[32] W. K. den Otter and W. J. Briels. The calculation of free-energy differences by
constrained molecular-dynamics simulations. J. Chem. Phys., 109(11) :4139–4146,
1998.
[33] M. Rousset and G. Stoltz. Equilibrium sampling from nonequilibrium dynamics. J.
Stat. Phys., 123(6) :1251–1272, 2006.
[34] M. Sprik and G. Ciccoti. Free energy from constrained molecular dynamics. J. Chem.
Phys., 109(18) :7737–7744, 1998.
[35] D. Talay. Approximation of upper Lyapunov exponents of bilinear stochastic differen-
tial systems. SIAM J. Numer. Anal., 28(4) :1141–1164, 1991.
[36] D. Talay and L. Tubaro. Expansion of the global error for numerical schemes solving
stochastic differential equations. Stochastic Anal. Appl., 8(4) :483–509 (1991), 1990.
[37] N.G. van Kampen. Statistical mechanics for trimers. App. Sci. Res., 37 :67–75, 1981.
[38] F. Wang and D.P. Landau. Determining the density of states for classical statisti-
cal models : A random walk algorithm to produce a flat histogram. Phys. Rev. E,
64 :056101, 2001.

19