3 Dynamique moléculaire
Dans le chapitre précédent, nous avons évoqué que le potentiel d’interaction \(U\) contient toute la physique d’un problème donné. Dans ce chapitre, nous allons essayer de mettre en place un potentiel \(U\) pour des interactions interatomiques simples, mais dont on va voir qu’il contient une richesse suffisante pour simuler gaz, liquide et solide.
4 Physique classique, physique quantique
À petite échelle, les interactions entre atomes et molécules sont décrites par la mécanique quantique. Moyennant quelques approximations, il est possible de calculer l’énergie potentielle d’un système de particules à partir des principes de la physique quantique (par exemple l’équation de Schrödinger). Les méthodes basées sur ce genre de calcul (par exemple la théorie fonctionnelle de la densité) sont appelées ab-initio : elles ne font pas d’hypothèse autre que la présence d’un certain nombre de protons, neutrons et électrons.
On pourrait penser que ce type d’approche est idéal pour la simulation de systèmes d’atomes et de molécules, mais l’emploi de ces méthodes se heurte à un impératif : le coût de calcul. La résolution, même approchée, de l’équation de Schrödinger est trop coûteuse pour simuler plus qu’une centaine de particules, même sur des supercalculateurs. Il est donc souvent nécessaire de formaliser des formes approximées des interactions entre particules qui évitent la résolution de l’équation de Schrödinger, au prix d’hypothèse supplémentaires.
5 Dynamique moléculaire
Puisque la résolution des équations de la mécanique quantique n’est plus nécessaire, on retourne dans le domaine de la physique dite classique : les énergies, quantité de mouvement, etc., sont des quantités continues. Ce type de méthode est appelé Dynamique Moléculaire. L’objet central de la méthode est (encore et toujours) le potentiel d’interaction, qui établit un compromis entre trois caractéristiques :
- La précision : c’est une mesure d’à quel point le potentiel est représentatif d’une interaction que l’on peut mesurer expérimentalement (ou calculer ab-initio) entre les atomes / molécules que l’on souhaite simuler. C’est en quelque sorte une erreur par rapport à la réalité pour une caractéristique particulière ;
- La transférabilité : c’est une mesure de la capacité d’un potentiel à faire des prédictions pour un large éventail de caractéristiques. Par exemple, un potentiel capable de prédire la distance entre les atomes dans une phase cristalline, le module d’élasticité de Young et l’énergie de surface est plus transférable qu’un potentiel qui ne prédit que l’une de ces quantités. Il n’existe pas de potentiel parfaitement transférable : une nouvelle situation doit toujours être testée ;
- Le coût de calcul : en général, plus un potentiel est précis, plus il est coûteux à calculer, et plus il est difficile de faire des simulations avec un grand nombre de particules.
C’est au scientifique qui fait la simulation de trouver le potentiel faisant le compromis le plus acceptable entre ces trois facteurs et le choix d’un potentiel fait très souvent l’objet de débat dans la littérature de la science des matériaux.
5.1 Interactions de paires
On a vu au précédent chapitre quelques exemples d’interactions de paires de particules. À l’échelle atomique, pour des particules / molécules électriquement neutres, on doit régulièrement modéliser deux effets importants qui sont décrits par des interactions de paires :
- La répulsion à courte distance : le prince d’exclusion de Pauli interdit à des électrons d’occuper le même état. En pratique, cela empêche deux atomes de se chevaucher.
- L’attraction à courte distance : même des atomes non-chargés peuvent s’attirer, grâce aux forces de dispersion de London, qui font partie des interactions de van der Waals. Le couple proton-électron, dû aux fluctuations quantiques, peut spontanément former un dipôle, qui induit un dipôle opposé dans un autre atome, et donc une force d’attraction. London a montré, avec une solution approchée de l’équation de Schrödinger, que cet effet était proportionnel à \(1 / r_{ij}^6\).
Ce dernier effet est à distinguer des liaisons covalentes (les liaisons chimiques) : les énergies d’interaction typiques des forces de van der Waals sont beaucoup plus faibles que les énergies des liaisons covalentes, de sortes qu’elles peuvent être brisées par la température seule. Il faut aussi distinguer les interactions ci-dessus de la force électrostatique de Coulomb (qui peut être attractive ou répulsive). Cette dernière s’exerce sur des distances bien plus grandes que l’attraction de London.
5.2 Potentiel de Lennard-Jones
Le potentiel de Lennard-Jones (1931), que l’on va noter \(U_\mathrm{LJ}\), modélise des interactions répulsives et attractives, et est donné par l’expression suivante :
\[ U_\mathrm{LJ}(r) = 4\varepsilon \left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right] \]
La Figure 5.1 montre l’énergie potentielle en fonction de \(r\).
On note plusieurs choses sur la Figure 5.1 :
- l’énergie potentielle est décroissante jusqu’à un minimum de \(-\varepsilon\) pour \(r = 2^\frac{1}{6}\sigma\), après quoi l’énergie augmente et tend vers zéro,
- le paramètre \(\varepsilon\) est homogène à une énergie (typiquement exprimé en eV, électron-Volt, soit \(\approx 1.6\cdot 10^{-19}\) J) et le paramètre \(\sigma\) est homogène à une longueur (typiquement exprimé en Å, Ångstrom, soit 0.1 nm)
- L’énergie tend vers \(+\infty\) en \(r = 0\), elle est nulle en \(r = \sigma\)
On a tiré de la première observation la zone (bleue) où la force est répulsive et la zone (rouge) où la force est attractive. La distance \(2^\frac{1}{6}\sigma\) où l’énergie est minimale, et donc la dérivée est nulle, est la distance d’équilibre. On affiche maintenant la force dérivée du potentiel, avec la convention qu’une force positive est une force d’attraction.
On remarque que la force montrée dans la Figure 5.2 est exprimée avec le ratio \(\frac{\varepsilon}{\sigma}\), homogène à une énergie sur une longueur, c’est-à-dire une force. En fait, les paramètres \(\varepsilon\) et \(\sigma\), avec une masse \(m\), permettent de définir un jeu d’unités naturelles :
| Quantité | Unité |
|---|---|
| Énergie | \(\varepsilon\) |
| Distance | \(\sigma\) |
| Force | \(\varepsilon \sigma^{-1}\) |
| Pression | \(\varepsilon \sigma^{-3}\) |
| Température1 | \(\varepsilon k_\mathrm B^{-1}\) |
| Masse | \(m\) |
| Densité | \(m\sigma^{-3}\) |
| Temps | \(\sqrt{{m\sigma^2}{\varepsilon^{-1}}}\) |
Ces unités peuvent donc servir à adimensionaliser les quantités d’intérêt (position, vitesse, etc.).
Paramétrisation
La paramétrisation d’un potentiel consiste à choisir les valeurs de paramètres (ici \(\sigma\) et \(\varepsilon\)) pour représenter un système physique réel. Le potentiel de Lennard-Jones, puisqu’il décrit des particules ne s’attirant que par des forces de van der Waals, est bien adapté à la description des gaz nobles.
Par exemple, pour modéliser l’argon, on pourra choisir \(\varepsilon\) = 0.01 eV et \(\sigma\) = 3 Å.
Pas de temps critique
On notera que si on utilise Velocity-Verlet, il faut choisir un pas de temps convenable. L’unité de temps \(\tau_\mathrm{LJ} = \sqrt{{m\sigma^2}{\varepsilon^{-1}}}\) est en fait établie à partir de la pulsation propre d’un système masse-ressort \(\omega = \sqrt{k / m}\). Puisque \(k\) est une force par longueur, \(\omega\) est homogène à \(\sqrt{\varepsilon m^{-1}\sigma^{-2}} = \tau_{LJ}^{-1}\). En pratique, on choisira un pas de temps dans l’intervalle \([10^{-3}\tau_\mathrm{LJ}, 10^{-2}\tau_\mathrm{LJ}]\).
Pour chaque système, il faut bien vérifier que le pas de temps choisi permet la conservation (en moyenne) de l’énergie totale du système.
Distance limite
L’énergie potentielle des forces de dispersion de London (la partie attractive du potentiel LJ) décroit vers zéro en \(1 / r^6\). Comparée à la décroissance de l’énergie potentielle de l’interaction électrostatique, en \(1 / r\), c’est une décroissance très rapide. En fait, la décroissance est si rapide que l’on peut approximer l’énergie d’une paire de particule par zéro si la distance séparant les particules est plus grande qu’une distance caractéristique \(r_\mathrm c\), qui est un rayon de troncation (cutoff radius en anglais), au-delà duquel on considère que les particules n’interagissent plus. On peut donc définir :
\[ U_\mathrm{LJ,cut}(r) = \begin{cases} U_\mathrm{LJ}(r) - U_\mathrm{LJ}(r_\mathrm c) & r \leq r_\mathrm c\\ 0 & r > r_\mathrm c \end{cases}\]
La figure Figure 5.3 montre une comparaison entre la version tronquée et la version complète du potentiel de Lennard-Jones.
On remarque sur la Figure 5.3 l’effet du terme \(-U_\mathrm{LJ}(r_\mathrm c\)) dans l’expression de \(U_\mathrm{LJ,cut}\) : l’énergie du potentiel tronqué est plus élevée, de façon à ce que l’énergie soit continue en \(r = r_\mathrm c\). Cependant, comme ce changement est dû à une constante, la position d’équilibre, et par conséquent les zones d’attraction et de répulsion, restent inchangées.
Coût de calcul
La formulation tronquée du potentiel LJ a un avantage majeur : le coût computationnel. Rappelons la formule de l’énergie potentielle totale d’un système à \(N\) particules :
\[ U(\{\vec r_k\}) = \sum_{i=1}^N \sum_{j = i + 1}^N U_2(r_{ij}) \]
Puisque l’on a deux sommes imbriquées, on doit effectuer \(O(N^2)\) opérations pour calculer l’énergie totale. Si, en revanche, on sait que les paires de particules plus éloignées que la distance \(r_\mathrm c\) ne contribuent pas à l’énergie, on peut réécrire la somme :
\[ U(\{\vec r_k\}) = \sum_{i=1}^N \sum_{j = i + 1\ |\ r_{ij}\leq r_\mathrm{c}}^N U_2(r_{ij}) \]
Supposons qu’une particule ait en moyenne \(\bar n\) voisins proches (qui satisfont le critère \(r_{ij} \leq r_\mathrm c\)), alors le nombre d’opérations requis pour calculer l’énergie est \(O(N\bar n)\) : si la densité (dont dépend \(\bar n\)) est constante, le coût de calcul croît donc linéairement par rapport aux nombres de particules, contrairement au calcul naïf où le coût croît quadratiquement. C’est une différence très significative que l’on va explorer au TP 3.
Recherche des voisins
Pour maintenir la complexité linéaire du calcul des interactions, il faut pouvoir établir la liste des voisins (neighbor list en anglais), soit l’ensemble des couples \(\{ (i, j)\ |\ r_{ij} \leq r_\mathrm c\}\). On ne peut pas simplement faire la liste de toutes les paires de particules puisqu’il y a de l’ordre de \(N^2\) paires.
On va donc utiliser l’algorithme des listes de cellules :
On divise l’espace de la simulation en une grille dont les cellules ont une taille \(b > r_\mathrm c\). À partir de la position de chaque particule, on peut déterminer sa cellule, et ainsi construire, pour chaque cellule, la liste des particules qu’elle contient. Pour trouver les voisins d’une particule, on doit uniquement chercher dans sa cellule et dans les cellules directement voisines. Puisque toutes les cellules sont plus grandes que \(r_\mathrm c\), deux cellules qui ne sont pas directement voisines (la zone hachurée de la Figure 5.4) ne contribuent pas à l’énergie potentielle.
Cette façon de calculer la liste des voisins permet de conserver un coût de calcul linéaire. En revanche, l’implémentation efficace de cet algorithme n’est pas aisée (surtout en Python), c’est pourquoi nous utiliserons une libraire externe lors du TP 3.
\(k_\mathrm{B}\) est la constante de Boltzmann qui intervient dans l’équation de la température \(\frac{3}{2} k_\mathrm B T = {\frac{1}{2}m\langle v^2 \rangle}\)↩︎