Analyse mathematique de methodes numeriques stochastiques en dynamique moleculaire

Introduction à la simulation moléculaire

La simulation moléculaire a été utilisée et développée au cours des 50 dernières années et son nombre d’utilisateurs est en expansion en dépit de ses contraintes intrinsèques à l’échelle spatiale et temporelle. La simulation moléculaire a deux objectifs majeurs:

Premièrement, elle peut être utilisée comme un microscope numérique pour faire des simulations à l’échelle microscopique: de nombreuses théories physiques ont été testées sur des ordinateurs. Historiquement, il n’y avait aucune bonne théorie analytique pour certains systèmes physiques (par exemple la physique de liquides simples), et l’observation des trajectoires simulées informatiquement a été très utile pour guider l’intuition des physiciens sur ce qui se passe dans le système (par exemple les mécanismes conduisant à la diffusion moléculaire). En particulier, les travaux pionniers sur les méthodes de Monte-Carlo [Metropolis et al., 1953], et la première simulation de dynamique moléculaire [Alder and Wainwright, 1956] ont été réalisées avec ces motivations. Aujourd’hui, la compréhension du comportement de la matière à l’échelle microscopique peut encore être difficile du point de vue expérimental parce que nous ne savons tout simplement ce qu’il faut chercher! Des simulations numériques sont alors un bon outil pour tester certaines idées ou obtenir certaines données que l’on peut ensuite traiter et analyser afin d’aider à concevoir les montages expérimentaux.

Le deuxième objectif majeur de la simulation moléculaire, peut-être encore plus important que le précédent, est de calculer les quantités macroscopiques ou les propriétés thermodynamiques, par le biais de moyennes de certains fonctionnelles du système  . Dans ce cas, la simulation moléculaire est un moyen d’obtenir des informations quantitatives sur un système. Parfois, ces propriétés sont accessibles à travers des expériences, mais dans certains cas les calculs numériques sont indispensables parce que les expériences peuvent être infaisables ou trop coûteuses. Plus généralement, la simulation moléculaire est un outil pour explorer les liens entre les propriétés microscopiques et macroscopiques d’un matériau, permettant de répondre à des questions de modélisation telles que « Quels ingrédients microscopiques sont nécessaires pour observer un comportement macroscopique donné ? »

Décrire une transition

Dans cette section, on décrit ce qu’on entend par une transition entre deux états. Une transition doit être comprise ici comme soit

1. Transitions alchimiques: La collection de toutes les configurations possibles microscopiques répartis selon la mesure canonique associé à un hamiltonien en fonction de certains paramètres λ. Le paramètre λ est alors l’indice de l’état. Par exemple, l’immersion d’une molécule dans un solvant, serait caractérisée par les les Hamiltoniens H0 (pour la molécule dans le vide) et H1 (intégrant les interactions soluté-solvant); ou

2. Transitions indexées par une coordonnée de réaction (ou paramètre d’ordre): Dans ce cas, toutes les configurations microscopiques sont distribués selon la mesure canonique πc, et satisfont une contrainte macroscopique ξ(q) = z, où ξ ∶ D→ Rm, avec m petite comparativement à 3N. Dans ce cas, z est l’indice de l’état.

Dans ce qui suit, on définit précisément ces deux approches et on discute comment les différences d’énergie libre sont calculées.

Transformations alchimiques 

Comme mentionné précédemment, la transformation alchimique concerne des transitions indexées par un paramètre λ externe. Ce paramètre est indépendant des configurations microscopiques (q, p). Le paramètre λ peut par exemple modifier les constantes utilisées dans les champs de force empiriques telles que l’énergie ε ou la longueur σ dans le potentiel du Lennard-Jones (1.21), (voir [Chipot and Pohorille, 2007] pour plus d’exemples). Le nom « alchimique » fait référence au fait que la nature chimique des particules peut être modifiée dans la simulation en modifiant les paramètres du potentiel décrivant les interactions moléculaires.

Choisir une bonne coordonnée de réaction 

Le choix d’une «bonne» coordonnée de réaction est un sujet très discuté dans la littérature. Un but de l’analyse mathématique mené ici ou dans des travaux précédents (voir par exemple [lelièvre et al., 2008]) est de quantifier l’efficacité des algorithmes numériques en fonction de la coordonnée de réaction choisie. A minima, une coordonnée de réaction ξ ∶ Rn #→ Rm est une fonction lisse telle que:

1. m ≪ 3N;
2. Les états initiaux et finaux de la transition peuvent être décrits en termes deξ;
3. det(G) ≠ 0, pour que la fonction de partition Zz soit bien définie et plus précisément de sorte que l’ensemble des sous-variétés Σz réalise un feuilletage de l’espace des phases entre les états initiaux et finaux de la transition.

Rappelons qu’en raison de la métastabilité, les moyennes trajectorielles (1.5) convergent très lentement vers leur limite ergodique. De nombreuses méthodes ont été proposées pour surmonter cette difficulté, et on se concentre dans cette thèse sur la méthode adaptative de force biaisante (notée ABF) (voir [Darve and Pohorille, 2001; Hénin and Chipot, 2004]). Pour introduire la méthode ABF (voir la Section 1.4.2 pour plus de détails), un ingrédient majeur est la coordonnée de réaction, ξ = (ξ1, …, ξm) ∶ Rn #→ Rm, ξ(q) = z. La fonction ξ peut être comprise comme une fonction telle que ξ(Xt) est une variable lente par rapport à Xt . On peut dire que ξ décrit les états métastables de la dynamique associée au potentiel V . Quand l’échantillonnage de la mesure canonique e −βV (1.11) est difficile, ξ est un bon choix si la mesure e−β(V −A○ξ) est plus facile à échantillonner que (1.11).

ABP vs ABF

Il est important de réaliser que les deux dynamiques (1.77) et (1.78) sont différentes. En particulier, il n’y a aucun raison que Ft soit un gradient en (1.78). Comparons maintenant les deux méthodes: la dynamique adaptative ABP et ABF.

1. La principale difficulté pour la dynamique ABP (1.64) est qu’aucun état stationnaire ne peut être atteint «spontanément»: Le potentiel biaisant At évolue en permanence quand de nouveaux échantillons sont ajoutés dans les régions visitées. Ceci est lié au fait que la quantité calculée (At) est définie à une constante additive près. Ces méthodes nécessitent donc un paramètre de réglage de la pénalisation, qui va artificiellement à zéro quand le temps évolue, pour que dAt/dt converge vers zéro dans la limite t → +∞ (cf. le paramètre γ dans (1.65)). Une telle difficulté ne se trouve pas dans la méthode ABF pour laquelle les espérances conditionnelles peuvent en effet atteindre un état stationnaire;

2. Un autre argument en faveur des méthodes ABF est qu’il vaut mieux numériquement estimer un gradient et ensuite l’intégrer, plutôt qu’estimer une fonction, puis la différencier. Pour les méthodes ABP, le biais At est estimé, et son gradient, qui est nettement plus bruté, est ensuite utilisé pour biaiser la dynamique;

3. Un intérêt de la dynamique ABP est qu’elle peut facilement s’appliquer à une coordonnée de réaction discrète, tandis que la dynamique ABF nécessite une coordonnée de réaction continue pour que la force moyenne soit définie.

Table des matières

1 Introduction
1.1 Introduction à la simulation moléculaire
1.1.1 Calcul de moyennes
1.1.1.1 Méthode de Monte Carlo standard
1.1.1.2 Dynamique moléculaire
1.1.2 Ensembles thermodynamiques
1.1.2.1 Ensemble microcanonique
1.1.2.2 Ensemble canonique
1.1.2.3 Ensemble isotherme-isobare
1.1.2.4 Ensemble grand-canonique
1.1.3 Energie potentielle
1.1.3.1 Potentiel d’interaction Ab initio
1.1.3.2 Potentiel empirique
1.2 Méthodes d’échantillonnage
1.2.1 Algorithme de Metropolis-Hastings
1.2.2 Algorithme Hybrid Monte-Carlo (HMC)
1.2.3 Dynamique moléculaire
1.2.3.1 Dynamique Hamiltonienne
1.2.3.2 Dynamique de Langevin
1.2.3.3 Dynamique de Langevin amortie
1.2.3.4 Métastabilité
1.3 Energie libre
1.3.1 Définition de l’énergie libre
1.3.2 Décrire une transition
1.4 Les méthodes adaptatives
1.4.1 Méthode adaptative du potentiel biaisant (ABP)
1.4.2 Méthode adaptative de force biaisante (ABF)
1.4.2.1 Formule de la co-aire et force moyenne
1.4.2.2 La méthode ABF
1.4.3 ABP vs ABF
1.4.4 Barrières d’énergie libre
1.5 Convergence des dynamiques non-adaptatives et adaptatives
1.5.1 Entropie et Inégalité de Sobolev Logarithmique
1.5.2 Convergence de la dynamique de Langevin suramortie
1.5.3 Convergence de la dynamique ABF
1.5.4 Exemple numérique
1.6 Contributions de cette thèse
1.6.1 Méthode adaptative de force biaisante projetée (PABF)
1.6.2 Solution locale et globale pour une équation de Fokker-Planck
nonlocal liée aux méthodes adaptives de force biaisante
1.6.3 Organisation du manuscrit
2 Long-time convergence of an adaptive biasing force method: Variance reduction by Helmholtz projection
2.1 Introduction
2.1.1 The model
2.1.2 Metastability, reaction coordinate and free energy
2.1.3 Adaptive biasing force method (ABF)
2.1.4 Projected adaptive biasing force method (PABF)
2.2 Numerical experiments
2.2.1 Presentation of the model
2.2.1.1 Potential functions
2.2.1.2 Reaction coordinate and physical parameters
2.2.1.3 Numerical methods and numerical parameters
2.2.2 Comparison of the methods
2.2.2.1 Metastability
2.2.2.2 Variance reduction
2.2.2.3 Free energy error
2.2.2.4 Distribution
2.3 Longtime convergence of the projected ABF method
2.3.1 Helmholtz projection
2.3.1.1 Helmholtz decomposition
2.3.1.2 Minimization problem and projection onto a gradient
2.3.2 The projected ABF dynamics
2.3.3 Precise statements of the longtime convergence results
2.3.3.1 Entropy and Fisher information
2.3.3.2 Convergence of the PABF dynamics
2.4 Proofs
2.4.1 Proof of Proposition 2.10
2.4.2 Proof of Corollary 2.13
2.4.3 Proof of Corollary 2.14
2.4.4 Proof of Theorem 2.15
3 Local and global solution for a nonlocal Fokker-Planck
3.1 Introduction
3.2 Notations and preliminaries
3.3 Local existence
3.3.1 Existence of mild solution
3.3.2 Regularity of the solution
3.3.3 Proof of Theorem 3.2
3.4 Diffusion equation and weak maximum principle
3.5 A-priori estimates for solutions and global existence
3.5.1 Polynomial and universal a-priori bounds
3.5.2 Global existence: Proof of Theorem 3.6
4 General Conclusion

Télécharger le rapport complet