Dynamique moléculaire

Comme dans toute simulation numérique, le temps évolue de manière discrète : le temps est découpé en une suite d'instants ${\displaystyle t_{init},t_{1},t_{2},\ldots ,t_{fin}}$, séparés par un intervalle très court appelé "pas-de-temps" ou "time-step", noté ${\displaystyle \delta t}$ et tel que ${\displaystyle t_{i+1}=t_{i}+\delta t}$. La simulation consistera alors à calculer la position et la vitesse des particules à chacun des instants, en utilisant les résultats obtenus à l'instant précédent.

Pour cela, on applique la deuxième loi de Newton à chaque instant, qui s'écrit pour une particule :

$${\displaystyle \sum {\vec {f}}_{i}(t)=m\cdot {\vec {a}}(t)}$$

où :

• ${\displaystyle {\vec {f}}_{i}}$ correspond aux différentes forces appliquées à la particule, soit par les particules voisines, soit par une influence extérieure.
• ${\displaystyle {\vec {a}}}$ est l'accélération de la particule dans un référentiel galiléen, assimilée à une masse ponctuelle.
• ${\displaystyle m}$ est la masse de la particule.

Dans le cas où les particules considérées sont des atomes, deux méthodes principales permettent le calcul des forces d'interaction :

• la dynamique moléculaire ab initio, où les forces sont calculées à partir des premiers principes de la mécanique quantique
• la dynamique moléculaire classique, où les forces dérivent d'un potentiel interatomique fixé empiriquement : historiquement, de type Lennard-Jones, aujourd'hui de type EAM (pour Embedded Atom Model).

La première approche a l'avantage d'être particulièrement précise, mais ne peut être envisagée que pour des systèmes de quelques dizaines de particules car elle nécessite un nombre de calculs numériques très important. La seconde est largement utilisée et permet d'étudier facilement des systèmes comportant quelques milliers voire millions de particules, mais pose la difficile question du choix d'un potentiel interatomique.

Le système étudié en dynamique moléculaire est typiquement une "boîte de simulation" parallélépipédique en 3D de dimensions ${\displaystyle L_{x},L_{y},L_{z}}$et contenant ${\displaystyle N}$ particules.

Pour simuler un cristal pseudo-infini dans une, deux ou trois dimensions, on appliquera des conditions aux limites périodiques. Cela signifie que lorsqu'une particule quitte la boîte de simulation par une face, elle est remise dans la boîte avec la même vitesse, au niveau de la face opposée.

Lors du calcul des forces, on devra tenir compte de cette périodicité de l'espace. En pratique, on distinguera dans la force d'interaction des termes à courte portée, qui ne seront pas affectés par la périodicité, c'est-à-dire que seules les particules les plus proches seront prises en compte, et un terme à longue portée, qui devra en tenir compte. Le terme à longue portée est généralement de type coulombien et sera calculé par la somme d'Ewald.

L'application de la seconde loi de Newton donne une équation différentielle du second ordre (dite équation du mouvement) pour chaque particule, qui doivent être résolues. Le but est de déterminer les coordonnées ${\displaystyle r_{i}(t+\delta t)}$ de la particule à l'instant ${\displaystyle t+\delta t}$, connaissant son accélération à l'instant précédent ${\displaystyle {\vec {a}}(t)}$que l'on a calculée connaissant les forces d'interactions ${\textstyle {\vec {f}}_{i}(t)}$.

Plusieurs algorithmes de résolution numérique peuvent être utilisés avec plus ou moins de succès selon le système étudié. Elles fournissent une solution approchée de l'équation, souvent basée sur des développements de Taylor. Un bon algorithme doit permettre une résolution rapide, tout en respectant la loi de conservation d'énergie, puisque l'on considère un système fermé.

Les algorithmes les plus courants sont :

• La méthode d'Euler, crée au XVIII^e siècle par Leonhard Euler très peu performant et imprécis, utilisé uniquement à des fins de démonstration et de pédagogie. Les méthodes de Runge-Kutta en sont des améliorations publiées en 1901.
• L'algorithme « Positions de Verlet » (voir Intégration de Verlet) proposé par Loup Verlet en 1967, qui présente de meilleures performances mais donne les positions à l’instant ${\displaystyle t+\delta t}$ et les vitesses à l’instant ${\displaystyle t}$[1].
• L'algorithme « Vitesses de Verlet », dérivé du précédent et donnant positions et vitesses à ${\displaystyle t+\delta t}$.
• L'algorithme « saut de grenouille », qui donne les vitesses à ${\textstyle t+{\frac {\delta t}{2}}}$ et les positions à ${\displaystyle t+\delta t}$.
• Les algorithmes de type « prédicteur correcteur (en) » qui font une approximation de la nouvelle position, vitesse et accélération à ${\displaystyle t+\delta t}$à l'aide d'un développement de Taylor. La nouvelle position permet le calcul des forces d'interactions ${\displaystyle {\vec {f}}_{i}(t+\delta t)}$, et de l'accélération ${\displaystyle {\vec {a}}_{c}(t+\delta t)}$ par application de la 2^e loi de Newton. On peut alors comparer la valeur d'accélération ${\displaystyle {\vec {a}}_{c}}$ avec celle prédite par les développements de Taylor, et appliquer un facteur correctif à cette dernière.

Le calcul des forces est primordial pour une simulation de dynamique moléculaire. Comme évoqué plus haut, deux approches peuvent être mises en œuvre pour le calcul des forces agissant sur les atomes :

• les méthodes dites ab initio prenant en compte la structure électronique des particules et le caractère quantique des interactions mises en jeu. Quel que soit le niveau de théorie appliqué, les forces sont obtenues par l'intermédiaire du Théorème de Hellmann-Feynman :
  ${\displaystyle \mathrm {F} _{\mathrm {X} _{\gamma }}=-{\frac {\partial \mathrm {E} }{\partial \mathrm {X} _{\gamma }}}=-\left\langle \psi \left|{\frac {\partial {\hat {\mathrm {H} }}}{\partial \mathrm {X} _{\gamma }}}\right|\psi \right\rangle .}$
• les méthodes dites classiques faisant intervenir des potentiels empiriques (ou champs de forces).

C'est la partie généralement la plus coûteuse en temps de calcul. Les principales interactions prises en compte sont :

• interactions à longue distance : forces de Coulomb (charges, dipôles, quadrupôles…) ;
• interactions à courte portée : modélisées par des potentiels de paire (ex. potentiel de Lennard-Jones, Buckinngham)
• interaction intra moléculaires : les modèles moléculaires utilisent en général aussi des modèles d'interactions intra moléculaires dits :
  • « 1-2 » élongation (stretching), en général un potentiel harmonique est utilisé,
  • « 1-3 » flexion (bending), en général un potentiel harmonique est utilisé,
  • « 1-4 » torsion, en général un potentiel sinusoïdal est utilisé,
  • …

Si la simulation comprend N particules, il y a ${\textstyle N\cdot {\frac {N-1}{2}}}$ interactions possibles. Le nombre de calculs à effectuer croît donc comme ${\displaystyle N^{2}}$, ce qui explique que les systèmes de grande dimension demandent un temps de calcul conséquent.

Pour réduire le nombre d'interactions à déterminer, on restreint les calculs d'interactions aux particules les plus proches, c'est-à-dire que l'on détermine une distance seuil, ou distance de troncature (cutoff), au-delà de laquelle les interactions sont négligées. Mais cette méthode nécessite de calculer les distances entre particules deux à deux, c'est-à-dire ${\textstyle N\cdot {\frac {N-1}{2}}}$ distances, pour déterminer la liste des proches voisins de chaque particule. Dans la méthode de Verlet[2], on peut augmenter le rayon de troncature, afin que la liste d'interacteurs établie puisse servir à plusieurs étapes de simulation.

D'autres approches permettent de réduire le nombre de distances à calculer pour déterminer les listes d'interacteurs, comme la méthode des listes par cellule (cell-lists, cell linked-lists)[3], une cellule pouvant être la boîte de simulation ou une de ses images :

• on découpe l'espace en petites cellules cubiques dont l'arête est supérieure ou égale à la distance seuil ;
• pour une molécule donnée, on calcule les distances avec les molécules de la même cellule et des cellules voisines.

Si l'interaction se ressent à très longue distance, typiquement si le potentiel est en 1/rⁿ avec n ≤ 2, la limitation par une distance de troncature introduit des erreurs. On a alors recours à des méthodes de type sommation sur un réseau infini, comme la sommation d'Ewald-Kornfeld ou l'expansion multipolaire de Ladd[4].

Le temps de calcul dépend énormément du champ de force utilisé lors de la simulation de la dynamique moléculaire. Les champs de force polarisable (AMOEBA, SIBFA... etc.) sont beaucoup plus coûteux en temps de calcul qu'un champ de force non polarisable (AMBER, CHARMM... etc.)