Comment utiliser le potentiel de Lennard Jones en dynamique moléculaire?

Je vous avais présenté un outil de dynamique moléculaire, développé en C++. Je vais maintenant vous en expliquer le fonctionnement, le paramétrage et présenter une nouvelle vidéo de simulation numérique.

Principes

Soit Φij(r) le potentiel d'interaction de paire entre la particule i avec la particule j. La somme des gradients de potentiel correspondants donne la force d'interaction appliquée sur chaque particule i :

Si l'on fige le temps et que l'on fait les calculs, on connait alors la force d'interaction appliquée à chaque particule, ce qui nous donne l'accélération a=F/m. En appliquant cette force un temps dt, on connait la variation de vitesse dv=a*dt. On a donc la nouvelle vitesse au point r+dr.

Il y a plusieurs algorithmes possibles, moi j'ai utilisé l'algorithme de Beeman. Pour ceux qui sont familiers avec le développement de Taylor d'une fonction, cet algorithme est une approximation à l'ordre 2. Le choix de l'algorithme impacte sur la précision et le temps de calcul. On peut par exemple aller à l'ordre 4 avec l'algorithme de Runge-Kutta.

On itère par pas de temps dt:

- Au temps t, on calcule la nouvelle position en t+dt, à l'aide des vitesses calculées au coup d'avant : r(t+dt)=r(t)+v(t)*dt+a(t)*dt².
- Les nouvelles positions calculées, on en déduit les forces à l'aide du potentiel comme expliqué précédemment.
- On calcule la variation de vitesse dv=a(t)*dt et on en déduit la nouvelle vitesse v(t+dt)=v(t)+a(t)*dt.

Potentiel de Lennard-Jones

Le potentiel utilisé est celui de Lennard-Jones qui se rapproche de la modélisation des Gaz de Van der Waals.

Il est composé d'une partie attractive de longue portée en puissance 6 et d'une partie répulsive à courte portée en puissance 12. Par la méthode du col ou en utilisant le lagrangien, on voit que, entre deux particules, il y aura un minimum à ce potentiel et donc un point d'équilibre ou une zone d'oscillation. Avec un grand nombre de particules, le comportement devient chaotique. Cependant cela permet de comprendre l'équilibre statistique de la distribution des vitesses. Sur le schéma ci-contre j'ai tracé la droite d'énergie nulle et colorié la zone d'énergie potentielle dans laquelle une particule seule et d'énergie initiale nulle oscillerait.

Conditions aux bords

J'ai modélisé deux types de conditions aux bords : un espace vide infini et un espace fini. Pour délimiter l'espace fini, on pourrait entourer le tout de molécules fixes interagissant avec celles à l'intérieur : boite non déformable. Mais de manière plus simple, on a défini un espace carré de taille fixe dont le bas communique avec le haut et la droite communique avec la gauche.

Variantes

Il est alors intéressant de réaliser des variantes. Par exemple, en enlevant la partie attractive du potentiel. Cela permet de simuler le comportement de boules de billards qui ne font que s'entrechoquer.

D'ailleurs, si vous pensez que cette approximation est moins bonne que l'approche mécanique consistant à évaluer leur distance puis à calculer les quantités de mouvement après un choc potentiel, pensez-y à nouveau et demandez-vous quelle est l'interaction physique qui est à l'origine des chocs ? Car la nature me semble plus proche de cette approche utilisant un potentiel que celle définissant une loi mécanique à respecter.

Nouvelle vidéo