La dynamique moléculaire

Le mouvement atomique, entraînant le mouvement moléculaire est associée à la loi de Newton :

où \(F_i\) représente la somme des forces exercées sur l'atome \(i\) de masse \(m_i\). \(a_i\) est son accélération.

On note bien que la notion de mouvement est ancrée dans cette relation car l'accélération est reliée à la vitesse, elle-même reliée à la position de la particule considérée.

La force sur un atome est obtenue à partir du calcul de l'énergie du système, généralement par une approche de type « champ de force » comme nous l'avons vu dans un chapitre précédent. Notons toutefois que cette énergie peut être une énergie qui contiendrai des termes issus de la mécanique quantique (par exemple la contribution de l'énergie électronique à l'énergie totale peut être obtenue par la résolution de l'équation de Schrödinger électronique).

Connaissant donc l'énergie du système, appelée E, on peut remonter à la force localisée sur un atome i par la relation très classique :

Il s'agit en fait de dériver l'énergie du système par rapport à la position de l'atome i. Pratiquement cela revient à estimer l'évolution de l'énergie lors d'un très petit déplacement de cet atome.

Cette relation peut être associée à la loi de Newton et l'on voit apparaître une relation mêlant l'accélération d'un atome et sa position :

On remarque que dans la mesure où l'accélération est la dérivée seconde de la position par rapport au temps, cette relation fait donc bien intervenir la notion de temps, de position et de force. Il s'agit maintenant d'intégrer cette fonction par rapport au temps pour obtenir une simulation du mouvement des atomes.

En fait, pratiquement, les positions ne peuvent pas être calculées directement avec cette relation dans la mesure ou une solution analytique n'existe pas pour cette relation. Dans un premier temps les accélérations sur les atomes sont calculées avec la loi de Newton. Ensuite, on estime les vitesses en intégrant les accélérations. Les positions sont alors aussi obtenues en intégrant les vitesses :

Les positions au temps \({t}\) sont utilisées pour estimer les positions au temps \({t}+\Delta t\), elles-mêmes ensuite utilisées pour obtenir les positions au temps \(t+2\Delta t\) et ainsi de suite.

Il existe différents algorithmes qui permettent ces estimations de nouvelles positions. L'une des plus simples est l'algorithme dit « leapfrog », traduit par « saut de grenouille ». Ici, les vitesses et les coordonnées ne sont pas déterminés aux même temps \(t\) , mais aux temps \(t+\frac{1}2\Delta t\) et \(t-\frac{1}2\Delta t\) respectivement. Elle est schématisée dans l'animation suivante :

Le calcul d'une trajectoire de dynamique moléculaire utilisant un algorithme de type LeapFrog peut se résumer au protocole suivant :

1. résoudre l'équation permettant d'obtenir les accélérations au temps t :

2. déterminer les vitesses aux temps \(t+\frac{1}2\Delta t\) à partir de leur valeur à \(t-\frac{1}2\Delta t\) :

Le développement de Taylor permettant d'accéder à ces valeurs entraîne nécessairement une approximation dans les résultats obtenus, dans la mesure où il est tronqué au 2ième terme. La validité de ce développement de Taylor est d'autant plus grande qu'il est réalisé sur une valeur de dt petite.

Par ailleurs, dans la mesure où cet algorithme a pour objectif de représenter correctement le mouvement des atomes les uns par rapport aux autres, il doit être en mesure de décrire les mouvements moléculaires les plus rapides. Cette contrainte existe dans tous les cas où un signal (ici un mouvement moléculaire) doit être représenté par un algorithme menant à des valeurs discrètes.

On parle alors de fréquence de Nyquist. La fréquence de Nyquist dans le cadre de systèmes moléculaires correspond à la fréquence de la vibration interatomique la plus rapide. Cette vibration est la vibration entre les atomes lourds et l'atome d'hydrogène. Elle vaut environ 1 fs. C'est pourquoi le pas d'intégration dt dans les algorithmes ne peut pas être supérieur à cette valeur. Une possibilité pour l'augmenter est de figer les liaison X-H dans la molécule (partant du principe qu'elles influences peu le comportement d'une macromolécule, comme une protéine par exemple). La fréquence de Nyquist devient alors celle de la vibration C-C et le pas d'intégration peut alors atteindre 2fs, soit une augmentation de 100 %, doublant ainsi le temps de la simulation pour le même temps de calcul.