Méthodes de simulation

Les simulations moléculaires constituent une approche complémentaire aux approches expérimentale, analytique et théorique pour l’étude de systèmes physiques ou chimiques. En effet, dans de nombreux cas, la réalisation d’expériences et leur interprétation restent limitées par la complexité des phénomènes et le nombre de paramètres entrant en jeu. Les approches analytiques permettent une compréhension complète des phénomènes, mais en raison des approximations qui doivent être faites, restent limitées à des systèmes très simples (gaz parfait, modèle d’Ising, cristal d’Einstein …). Les simulations numériques donnent accès à un certain nombre de grandeurs microscopiques mais aussi macroscopiques, grâce aux lois de la physique statistique. Les méthodes de simulation existantes sont nombreuses et variées. On peut ainsi avoir accès à différents niveaux de description du système, allant de la description fine des liaisons chimiques grâce aux méthodes quantiques jusqu’à une représentation à l’échelle mésoscopique, où même macroscopique. On peut ainsi avoir accès à différentes informations, des informations dynamiques ou thermodynamiques à différentes échelles, depuis l’énergie de liaison d’une liaison chimique, jusqu’à la vitesse d’écoulement d’un fluide, en passant par les coefficients de diffusion et les enthalpies de solvatation des ions. Dans le cas des simulations microscopiques et mésoscopiques, le formalisme de la physique statistique permet de faire le lien entre les grandeurs de sortie de la simulation et les grandeurs macroscopiques. Les méthodes de simulations sont de plus en plus utilisées en complément des approches expérimentales. Lorsque le modèle et la méthode utilisés permettent de reproduire certaines données expérimentales, on a alors accès au détail microscopique souvent inaccessible expérimentalement. De plus, le contrôle de chacun des paramètres permet de découpler l’effet des différents facteurs, de comprendre l’influence de chacun d’entre eux d’extraire les plus importants, et d’apporter une vue nouvelle sur la compréhension des phénomènes étudiés. Elle permet de comparer les caractéristiques d’un modèle à l’expérience et ainsi d’apporter une meilleure compréhension des phénomènes. La simulation moléculaire peut aussi être utilisée pour la prédiction de résultats expérimentaux. Elle permet d’avoir accès à des conditions qui ne sont pas accessibles expérimentalement (domaine de pression ou de température extrêmes) et permet d’obtenir facilement des informations sur le comportement de systèmes réels pour lesquels la conduite d’une série d’expériences serait complexe et fastidieuse, voire impossible. Enfin, elle est de plus en plus souvent utilisée pour l’interprétation de données expérimentales, notamment en RMN où les déplacements chimiques expérimentaux sont parfois attribués en fonction des Méthodes de simulation 34 Chapitre II : Méthodes de simulation résultats issus de simulations [179], dans le cas de l’étude de l’adsorption de molécules dans les matériaux poreux [180] ou de la localisation des cations dans les zéolithes [181]. La simulation est aussi utilisée en complément des approches théoriques ou analytiques. Ce type d’approche est souvent limité à des systèmes très simples pour lesquels les équations de la thermodynamique statistique sont solubles analytiquement (gaz parfait, modèle d’Ising…). L’étude d’un phénomène passe alors par l’introduction d’un grand nombre d’approximations, qui sont souvent trop draconiennes pour espérer obtenir un modèle s’approchant du système réel. La simulation moléculaire permet de résoudre exactement par voie numérique les problèmes de la physique statistique qui ne sont pas solubles analytiquement. On peut alors tester une théorie ou explorer les propriétés de systèmes modèles afin d’obtenir une meilleure compréhension du système via ce que l’on appelle une « expérience » numérique. Je présenterai ici un aperçu des méthodes de simulation moléculaire classiques, puis j’exposerai en détail le principe de la méthode Monte Carlo et de la dynamique moléculaire, qui sont les méthodes que j’ai utilisées au cours de ma thèse. 

Principes des méthodes de simulation 

Calcul de grandeurs thermodynamiques : apports de la physique statistique 

D’un point de vue microscopique, un système de N particules se caractérise par un ensemble de microétats, l’ensemble de ces microétats constituant l’espace des phases. Cette espace de dimension 6N a pour coordonnées les positions des particules {ri}C et leur quantité de mouvements {qi}C. Le formalisme de la physique statistique permet de relier les grandeurs accessibles à partir de ces microétats et les grandeurs macroscopiques. Ainsi, la valeur d’une observable A macroscopique correspond à la valeur moyenne de cette observable sur l’ensemble des microétats du système en prenant en compte leurs probabilités d’apparition de ce microétat. Cela revient à moyenner sur l’ensemble des configurations microscopiques pondérées par leur probabilité de Boltzmann : ⟨A⟩ = ∑ C AC Pboltz({ri}C, {qi}C) (II.1) où C représente une configuration du système, {ri}C correspond à l’ensemble des positions des atomes, {qi}C leur quantité de mouvements dans cette configuration et AC la valeur de la grandeur microscopique correspondant à A en C. Pboltz correspond à la probabilité de se trouver dans l’état C I . Cette expression ne se résout quasiment jamais analytiquement. L’évaluation de la valeur de AC dans toutes les configurations C possibles du système (ainsi que de leur probabilité de Boltzmann) est impossible pour des systèmes de taille raisonnable (pour 100 atomes, on a I. L’expression de la probabilité d’apparition d’un état dépend des conditions dans lesquelles le système se trouve. Ainsi un système pour lequel on a fixé la température, le nombre de particules et le volume (on parle d’ensemble (N,V,T) ou ensemble canonique), cette probabilité de Boltzmann est égale à exp(−βEC) QN,V,T où QN,V,T est la fonction de partition de l’ensemble canonique. Si on impose d’autres contraintes, l’expression de cette probabilité est différente. L’expression de ces probabilités sera détaillée dans la section II.1.3.Principes des méthodes de simulation 100 × 3 coordonnées dans l’espace et 100 × 6 coordonnées dans l’espace des phases). Il est alors nécessaire de générer un ensemble de configurations représentatives du système en fonction de leurs probabilités de Boltzmann et on approxime la moyenne de A à la valeur de la moyenne sur ces configurations. On s’affranchit ainsi du calcul de ces probabilités de Boltzmann et donc du calcul explicite de la fonction de partition du système, dont le calcul est en général complexe voire impossible. La valeur de A pourra alors être obtenue en faisant la moyenne de AC sur les configurations générées par la simulation. ⟨A⟩ ≃ ∑ C g´en´er´ees suivant Boltzmann AC (II.2) La difficulté de la simulation réside alors dans la génération de cet ensemble de configurations, suffisamment représentatif du système pour que la description de l’espace des phases soit complète sans pour autant que leur génération ne soit trop coûteuse en temps de calcul.

Échantillonnage de l’espace des phases

 Il existe plusieurs méthodes pour explorer l’espace des phases. Je présenterai ici brièvement les deux plus courantes : la dynamique moléculaire et la simulation Monte Carlo. II.1.2.1 Dynamique moléculaire La dynamique moléculaire est la méthode la plus « naturelle » pour générer des configurations. Elle consiste à suivre le mouvement des particules du système en fonction du temps. Pour cela, il est nécessaire d’intégrer les équations de Newton au cours de la simulation, ce qui conduit à introduire un pas de temps dt pour mener à bien cette intégration. On obtient alors pour chacun des pas le couple ({ri}C,{ dri dt }C). Les configurations ainsi générées peuvent alors être exploitées pour accéder à des grandeurs macroscopiques en calculant la moyenne temporelle de la valeur microscopique correspondante : ⟨A⟩ ≃ ∫ t0+∆t t0 AC(t) dt (II.3) où ∆t représente la durée simulée et t0 l’instant initial de la simulation. Il est alors possible d’étudier les propriétés thermodynamiques du système. La simulation de dynamique moléculaire nous apporte également des informations sur la dynamique du système modélisé. Pour obtenir une valeur de A fiable, il faut que la dynamique générée par simulation couvre un domaine de temps suffisant pour permettre au système de parcourir une partie représentative de l’espace des phases. Suivant les phénomènes que l’on désire étudier et leur temps caractéristique, la durée de la simulation à effectuer varie. L’étude des propriétés de diffusion d’un liquide demande un temps de simulation bien moindre que l’étude de sa diffusion dans les canaux d’une zéolithe ou de son processus de cristallisation. Les phénomènes qui ont des temps caractéristiques trop long (> 1 µs) ne pourront jamais être observés avec les moyens de calcul dont nous disposons actuellement. Dans les cas où il existe des grandes barrières de potentiel entre deux états de faible énergie, un échantillonnage trop court par la dynamique moléculaire conduira à l’étude d’un seul des deux puits d’énergie (figure II.1). 36 Chapitre II : Méthodes de simulation Figure II.1 : Représentation d’une surface d’énergie potentielle. La dynamique moléculaire (courbe rouge) explore la surface en suivant les équations du mouvement. Le passage d’un puits de potentiel à un autre se fait en passant par le col d’énergie potentiel les séparant. La simulation Monte Carlo (en bleu) parcourt l’espace des configurations de manière aléatoire en fonction de la probabilité de Boltzmann de chaque état. Le passage entre deux états séparés par des barrières d’énergie importantes peut se faire sans passer par les états intermédiaires de haute énergie. II.1 Principes des méthodes de simulation 37 La méthode de Monte Carlo permet d’échantillonner les états d’énergie faible en traversant les barrières énergétiques sans toutefois obtenir des informations sur leur amplitude. On utilise pour l’étude des phénomènes ayant des temps caractéristiques très grands, d’autres méthodes comme la méthode de transition path sampling (échantillonage de chemin de transition) [182] ou la dynamique biaisée, qui permettent de traverser les cols d’énergie potentielle ou la métadynamique [183] qui permet d’explorer l’ensemble de la surface d’énergie potentielle. 

La simulation Monte Carlo 

La simulation Monte Carlo restreint l’étude de l’espace des phases à l’étude de l’espace des configurations. Elle consiste à générer aléatoirement un ensemble de configurations en respectant leur probabilité de Boltzmann. Chaque configuration est générée à partir de la précédente en utilisant des mouvements type aléatoires (translation, rotation …). Le choix des mouvements est déterminant pour avoir un bon échantillonnage de l’espace des phases. Ces mouvements peuvent être inspirés de mouvements physiques ou être complètement irréels. Par exemple dans le cas des simulations de zéolithes cationiques, il est nécessaire d’introduire un mouvement de « jump », qui consiste en une translation de quelques ångtröms du cation. Ces mouvements non physiques permettent de passer des barrières d’énergie potentielle importantes. Il est nécessaire de bien prendre en compte tous les degrés de liberté du système dans le choix des mouvements Monte Carlo utilisés. La génération aléatoire de configurations par des pas non physiques, permet de passer directement entre deux configurations séparées par une barrière de potentiel importante. On limite alors les problèmes de piégeage dans un état métastable. Toutefois, ce type de problème peut apparaître dans le cas de systèmes très fortement liés, dans le cas où les pas choisis ne sont pas adaptés, ou si leur nombre est insuffisant. La difficulté est de générer un ensemble de configurations représentatives du système. Pour cela on utilise l’algorithme introduit par Metropolis et al. [184] qui permet de générer les configurations en suivant leurs probabilités de Boltzmann I . Cette méthode est très utilisée dans l’étude des phénomènes possédant des grandes barrières d’activation. C’est le cas des phénomènes d’adsorption dans les zéolithes cationiques : le passage des cations à travers les fenêtres hexagonales de la faujasite, la diffusion de molécules à l’intérieur d’une maille de zéolithe afin d’atteindre leurs sites d’adsorption sont des phénomènes très lents qui ne pourraient pas être étudiés par dynamique moléculaire. On perd toutefois avec ce type d’approche toute information sur les propriétés du système dépendantes du temps, seules les propriétés à l’équilibre sont accessibles par cette méthode. II.1.3 Les ensembles statistiques Un ensemble statistique est un ensemble de copies d’un même système (une infinité), chacune étant dans un état possible du système. La notion d’ensemble statistique fait le lien entre les états discrets du système au niveau microscopique décrits par le formalisme de la mécanique quantique (appelé microétats), et les observables macroscopiques. L’ensemble de ces copies peut être placé dans des conditions thermodynamique variées, on définit alors des ensembles statisI. Par exemple, dans l’ensemble canonique la probabilité de Boltzmann d’une configuration C d’énergie EC est proportionnelle à e −βEC où β = 1 kBT .. Ainsi dans l’ensemble microcanonique (N,V ,E) le système évolue de manière isolé, à volume et nombre de particules constant. On peut mettre le système en équilibre avec un thermostat, à volume et nombre de particules constant, on est alors dans l’ensemble canonique. L’expression de la probabilité d’apparition d’un microétat d’un système (ou probabilité de Boltzmann) diffère suivant l’ensemble statistique dans lequel on se place. Suivant les conditions dans lesquelles évolue le système qu’on étudie, il appartient à des ensembles statistiques différents. Je présenterai ici les ensembles statistiques que j’ai utilisés au cours de ma thèse : l’ensemble microcanonique, l’ensemble canonique, l’ensemble grand-canonique, l’ensemble isobare-isotherme et l’ensemble osmotique (cf. figure II.2).