Dynamique moléculaire : modèle de Lennard-Jones à deux dimensions

1. Introduction

Ce document présente la dynamique moléculaire de sphères se déplaçant sur un plan et interagissant avec un potentiel de Lennard-Jones ([1][2]). La méthode utilisée est similaire à celle exposée dans modèle des sphères molles à deux dimensions. Elle consiste à utiliser une grille pour accélérer la recherche des proches voisins. On utilisera aussi une liste de voisins, comme l'a fait L. Verlet ([1]).

2. Mise en équation et méthode numérique

2.a. Forces intermoléculaires

Dans le modèle de Lennard-Jones, les molécules sont sphériques et interagissent par le potentiel suivant :

$$U\left(r\right)=4\epsilon \left[\right(\frac{d}{r}{)}^{12}-\left(\frac{d}{r}{)}^6\right]\text{(1)}$$

La force intermoléculaire correspondante est

$$\overrightarrow{f}\left(r\right)=\frac{48\epsilon }{d}\left[\right(\frac{d}{r}{)}^{13}-\frac{1}{2}\left(\frac{d}{r}{)}^7\right]\overrightarrow{u_r}\text{(2)}$$

Voici le tracé de cette force :

from matplotlib.pyplot import *
import math
import numpy

def force(r):
    return 48.0*(math.pow(r,-13)-0.5*math.pow(r,-7))
r = numpy.arange(0.8,4.0,0.01)
f = r.copy()
for i in range(r.size):
    f[i] = force(r[i])
figure(figsize=(8,6))
plot(r,f)
xlabel("r/d")
ylabel("force")
axis([0,4,-3,10])
grid()
                  

figA.pdf

À courte distance (r<2^1/6d), la force est répulsive, avec une pente très forte. À grande distance, la force attractive en r^-13 modélise la force attractive de van der Waals.

La portée de cette force est en principe infinie, mais comme elle décroît très vite avec la distance, on introduit une distance de coupure r_c au delà de laquelle on considère que la force est nulle. On peut prendre par exemple r_c=3d.

3. Implémentation Python

import math
import numpy
import random
import pygame
            

Comme les sphères se déplacent sur un plan, elles sont représentées par des disques. La classe Disque contient les coordonnées du centre, la vitesse et l'accélérations :

class Disque:
    def __init__(self,x,y,vx,vy):
        self.x = x
        self.y = y
        self.vx  = vx
        self.vy = vy
        self.ax = 0.0
        self.ay = 0.0
            

Il s'agit d'une simple structure de données : on accédera directement aux données, pour les lire ou les modifier.

La classe suivante représente une cellule, qui contient une liste de disques. Les disques seront référencés par un indice sur un tableau. On stocke donc cet indice dans la cellule.

class Cellule:
    def __init__(self):
        self.ensemble_disques = []
        self.decal_x = 0.0
        self.decal_y = 0.0
        self.n = 0
    def ajouter(self,indice):
        self.ensemble_disques.append(indice)
        self.n += 1
    def enlever(self,indice):
        i=0
        cherche = True
        while (i<self.n) and cherche:
            if self.ensemble_disques[i] == indice:
                if self.n == 1 or i==self.n-1:
                    self.ensemble_disques.pop()
                else:
                    self.ensemble_disques[i] = self.ensemble_disques.pop()
                self.n -=1
                cherche = False
            else:
                i += 1
              

La classe Grille contient un tableau bidimensionnel de cellules (une liste de listes). Le constructeur, qui prend en argument la largeur de la grille et le nombre de cellules par dimension, crée le tableau. La fonction obtenir_cellule(i,j) renvoit la cellule dont les indices sont donnés en argument. Lorsque les indices sortent des bornes de la grille, la condition limite périodique est appliquée. Dans ce cas, un décalage doit être appliqué aux coordonnées des disques. Ce décalage est stocké dans la cellule renvoyée.

class Grille:
    def __init__(self,Nc,L):
        self.L = L
        self.Nc = Nc
        self.tableau = []
        for i in range(self.Nc):
            ligne = []
            for j in range(self.Nc):
                ligne.append(Cellule())
            self.tableau.append(ligne)
    def obtenir_cellule(self,i,j):
            decal_x = 0.0
            decal_y = 0.0
            if i<0:
                i += self.Nc
                decal_x = -self.L
            elif i>=self.Nc:
                i -= self.Nc
                decal_x = self.L
            if j<0:
                j += self.Nc
                decal_y = -self.L
            elif j>=self.Nc:
                j -= self.Nc
                decal_y = self.L
            cellule = self.tableau[i][j]
            cellule.decal_x = decal_x
            cellule.decal_y = decal_y
            return cellule
             

La classe Systeme effectue les calculs. Le constructeur prend en argument le nombre de sphères par dimension (la racine du nombre total), la densité, la distance de coupure de la force, et la taille de la couche supplémentaire pour l'utilisation de la liste de voisins (facultatif). Le constructeur calcule le nombre de cellules par dimensions Nc et la largeur des cellules lc, qui doit être supérieure à r_c+Δr.

class Systeme:
    def __init__(self,Nx,densite,rc,deltaR=0):
        self.Nx = Nx
        self.N = Nx*Nx
        self.densite = densite
        self.deltaR = deltaR
        self.rc = rc
        self.rc2 = rc*rc
        self.rv = rc+deltaR
        self.rv2 = self.rv*self.rv
        self.rayon = 0.5
        self.diam = self.rayon*2
        self.L = math.sqrt(math.pi/densite)*0.5*Nx
        self.aire = self.L*self.L
        self.demi_L = self.L*0.5
        self.Nc = int(self.L/self.rv)
        self.lc = self.L/self.Nc
        self.grille = Grille(self.Nc,self.L)
        self.liste_disques = []
        for d in range(self.N):
            self.liste_disques.append(Disque(0,0,0,0))
        self.energie = 0.0 
        self.viriel = 0.0
        self.Ecinetique = 0.0
        self.Epotentielle = 0.0
        self.pression = 0.0
        self.compteur = 0
        self.som_temp = 0
        self.som_pression = 0
        self.som_temp2 = 0
        self.som_pression2 = 0
             

La fonction suivante permet d'initialiser un disque dont l'indice est donné, avec une position et une vitesse données. L'indice est ajouté dans la cellule qui contient le disque.

    def init_disque(self,indice,x,y,vx,vy):
        disque = self.liste_disques[indice]
        disque.x = x
        disque.y = y
        disque.vx = vx
        disque.vy = vy
        i = int(x/self.lc)
        j = int(y/self.lc)
        cellule = self.grille.obtenir_cellule(i,j)
        cellule.ajouter(indice)
             

Pour initialiser le système, on place les sphères régulièrement espacées avec une vitesse de direction aléatoire, dont la norme est donnée. Une correction sur les vitesses est effectuée pour que la quantité de mouvement totale soit nulle. Une fois le système initialisé, on doit calculer les forces et créer la liste de voisins.

    def initialiser(self,vitesse):
        dx = self.L*1.0/(self.Nx)
        print("distance initiale = %f"%dx)
        if dx < self.diam:
            raise Exception("densite trop forte")
        else:
            dy = dx
            x = dx/2
            y = x
            px = 0.0
            py = 0.0
            for k in range(self.N):
                a = random.random()*math.pi*2.0
                vx = vitesse*math.cos(a)
                vy = vitesse*math.sin(a)
                px += vx
                py += vy
                self.init_disque(k,x,y,vx,vy)
                x += dx
                if x > self.L:
                    x = dx/2
                    y += dy
            for k in range(self.N):
                disque = self.liste_disques[k]
                disque.vx -= px/self.N
                disque.vy -= py/self.N 
        self.calculer_forces()
        self.construire_liste_voisins()
             

La fonction suivante déplace un disque à une position donnée, en tenant compte des conditions limites périodiques. Si le disque change de cellule, les deux cellules concernées sont mises à jour.

    def deplacer_disque(self,disque,indice,x1,y1):
        if x1<0:
            x1 += self.L
        elif x1>=self.L:
            x1 -= self.L
        if y1<0:
            y1 += self.L
        elif y1>=self.L:
            y1 -= self.L
        i = int(disque.x/self.lc)
        j = int(disque.y/self.lc)
        i1 = int(x1/self.lc)
        j1 = int(y1/self.lc)
        if i!=i1 or j!=j1:
            cellule = self.grille.obtenir_cellule(i,j)
            cellule.enlever(indice)
            cellule1 = self.grille.obtenir_cellule(i1,j1)
            cellule1.ajouter(indice)
        disque.x = x1
        disque.y = y1
           

La fonction suivante construit la liste de voisins, constituée des paires de sphères dont la distance est inférieure à r_v=r_c+Δr. Pour chaque sphères, les sphères voisines sont recherchées dans la même cellule et dans les 8 cellules voisines.

    def construire_liste_voisins(self):
        self.liste_voisins = []
        self.deplac_max = 0.0
        for i in range(self.Nc):
            for j in range(self.Nc):
                cellule = self.grille.obtenir_cellule(i,j)
                for k in range(-1,2):
                    for l in range(-1,2):
                        cellule1 = self.grille.obtenir_cellule(i+k,j+l)
                        for indice in cellule.ensemble_disques:
                                for indice1 in cellule1.ensemble_disques:
                                    if indice1<indice:
                                        disque = self.liste_disques[indice]
                                        disque1 = self.liste_disques[indice1]
                                        dx = disque1.x+cellule1.decal_x-disque.x
                                        dy = disque1.y+cellule1.decal_y-disque.y
                                        r2 = dx*dx+dy*dy
                                        if r2<= self.rv2:
                                            self.liste_voisins.append([indice1,indice])
           

La fonction suivante calcule les forces sur les disques et stocke le résultat dans les accélérations des disques. La boucle principale se fait sur les cellules de la grille. Pour chaque cellule, les cellules voisines sont explorées afin d'extraire les paires susceptibles d'interagir. Enfin si la distance entre les deux sphères est inférieure à la portée de l'interaction, on calcule la force et on affecte les valeurs des accélérations des deux disques. L'énergie potentielle et le viriel sont aussi calculés.

    def calculer_forces(self):
        for k in range(self.N):
            disque = self.liste_disques[k]
            disque.ax = 0.0
            disque.ay = 0.0
        self.Epotentielle = 0.0
        self.viriel = 0.0
        for i in range(self.Nc):
            for j in range(self.Nc):
                cellule = self.grille.obtenir_cellule(i,j)
                for k in range(-1,2):
                    for l in range(-1,2):
                        cellule1 = self.grille.obtenir_cellule(i+k,j+l)
                        for indice in cellule.ensemble_disques:
                                for indice1 in cellule1.ensemble_disques:
                                    if indice1<indice:
                                        disque = self.liste_disques[indice]
                                        disque1 = self.liste_disques[indice1]
                                        dx = disque1.x+cellule1.decal_x-disque.x
                                        dy = disque1.y+cellule1.decal_y-disque.y
                                        r2 = dx*dx+dy*dy
                                        if r2 < self.rc2:
                                            ir2 = 1.0/r2
                                            ir6 = ir2*ir2*ir2
                                            v = 24.0*ir6*(ir6-0.5)
                                            f = 2.0*v*ir2
                                            fx = f*dx
                                            fy = f*dy
                                            disque1.ax += fx
                                            disque1.ay += fy
                                            disque.ax -= fx
                                            disque.ay -= fy
                                            self.Epotentielle += 4.0*ir6*(ir6-1.0)
                                            self.viriel += v

            

La fonction suivante calcule les forces en utilisant la liste des voisins.

    def calculer_forces_avec_liste_voisins(self):
        for k in range(self.N):
            disque = self.liste_disques[k]
            disque.ax = 0.0
            disque.ay = 0.0
        self.Epotentielle = 0.0
        self.viriel = 0.0
        for paire in self.liste_voisins:
            disque = self.liste_disques[paire[0]]
            disque1 = self.liste_disques[paire[1]]
            dx = disque1.x-disque.x
            dy = disque1.y-disque.y
            if dx >= self.demi_L:
                dx -= self.L
            elif dx < -self.demi_L:
                dx += self.L
            if dy >= self.demi_L:
                dy -= self.L
            elif dy < -self.demi_L:
                dy += self.L
            r2 = dx*dx+dy*dy
            if r2 < self.rc2:
                ir2 = 1.0/r2
                ir6 = ir2*ir2*ir2
                v = 24.0*ir6*(ir6-0.5)
                f = 2.0*v*ir2
                fx = f*dx
                fy = f*dy
                disque1.ax += fx
                disque1.ay += fy
                disque.ax -= fx
                disque.ay -= fy
                self.Epotentielle += 4.0*ir6*(ir6-1.0)
                self.viriel += v

            

La fonction suivante calcule l'énergie cinétique, la pression (instantanée), et l'énergie totale. Cette dernière est importante pour contrôler la stabilité du schéma numérique d'intégration, car elle doit en principe rester constante. Toutes ces grandeurs sont divisées par le nombre de sphères. On calcule aussi les sommes qui seront utilisées pour le calcul des valeurs moyennes.

    def calculer_cinetique(self):
        self.Ecinetique = 0.0
        for k in range(self.N):
            disque = self.liste_disques[k]
            self.Ecinetique += 0.5*(disque.vx*disque.vx+disque.vy*disque.vy)
        self.pression = (self.Ecinetique+self.viriel)/self.aire
        self.Ecinetique /= self.N
        self.energie = self.Ecinetique+self.Epotentielle/self.N
        self.som_temp += self.Ecinetique
        self.som_temp2 += self.Ecinetique*self.Ecinetique
        self.som_pression += self.pression
        self.som_pression2 += self.pression*self.pression
        self.compteur += 1
             

La fonction suivante initialise le compteur et les sommes pour le calcul des valeurs moyennes.

    def init_moyennes(self):
        self.compteur = 0
        self.som_temp = 0.0
        self.som_pression = 0.0
        self.som_temp2 = 0.0
        self.som_pression2 = 0.0
             

les deux fonctions suivantes calculent la moyenne et l'écart-type de la température et de la pression :

    def temperature_moyenne(self):
        Tm = self.som_temp/self.compteur
        dTm = math.sqrt(self.som_temp2/self.compteur-Tm*Tm)
        return (Tm,dTm)
    def pression_moyenne(self):
        Pm = self.som_pression/self.compteur
        dPm = math.sqrt(self.som_pression2/self.compteur-Pm*Pm)
        return (Pm,dPm)
             

La fonction suivante effectue une ajustement des vitesses pour obtenir une température donnée.

    def ajuster_vitesses(self,T):
        Tm = self.som_temp/self.compteur
        fac = math.sqrt(T/Tm)
        for k in range(self.N):
            disque = self.liste_disques[k]
            disque.vx *= fac
            disque.vy *= fac
             

La fonction suivante effectue le pas élémentaire de la méthode de Verlet, en utilisant les cellules pour calculer les forces.

    def verlet(self,h,hd2):
        for k in range(self.N):
            disque = self.liste_disques[k]
            disque.vx += hd2*disque.ax
            disque.vy += hd2*disque.ay
            self.deplacer_disque(disque,k,disque.x+h*disque.vx,disque.y+h*disque.vy)
        self.calculer_forces()
        for k in range(self.N):
            disque = self.liste_disques[k]
            disque.vx += hd2*disque.ax
            disque.vy += hd2*disque.ay
            

La fonction suivante effectue le pas élémentaire de la méthode de Verlet, en utilisant la liste des voisins. Lors de la seconde évaluation des vitesses, la vitesse maximale (au carré) est aussi calculée, afin d'incrémenter le déplacement maximal des sphères. Si celui-ci dépasse Δr, la liste de voisins est reconstruite.

    def verlet_avec_liste_voisins(self,h,hd2):
        for k in range(self.N):
            disque = self.liste_disques[k]
            disque.vx += hd2*disque.ax
            disque.vy += hd2*disque.ay
            self.deplacer_disque(disque,k,disque.x+h*disque.vx,disque.y+h*disque.vy)
        self.calculer_forces_avec_liste_voisins()
        v2max = 0.0
        for k in range(self.N):
            disque = self.liste_disques[k]
            disque.vx += hd2*disque.ax
            disque.vy += hd2*disque.ay
            v2 = disque.vx*disque.vx+disque.vy*disque.vy
            if v2 > v2max:
                v2max = v2
        self.deplac_max += math.sqrt(v2max)*h
        if self.deplac_max*2.0 > self.deltaR:
            self.construire_liste_voisins()
            self.deplac_max = 0.0
            

Les deux fonctions suivantes effectuent n pas d'intégration

    def integration(self,h,n):
        hd2 = h/2
        for i in range(n):
            self.verlet(h,hd2)
    def integration_avec_liste_voisins(self,h,n):
        hd2 = h/2
        for i in range(n):
            self.verlet_avec_liste_voisins(h,hd2)
            

La fonction suivante effectue le dessin des disques avec pygame. Le rayon des disques est le paramètre d qui apparaît dans le potentiel de Lennard-Jones, ici égal à 1.

    def dessiner_disques(self,screen,echelle,couleur):
        for k in range(self.N):
            disque = self.liste_disques[k]
            pygame.draw.ellipse(screen,couleur,[(disque.x-self.rayon)*echelle,(disque.y-self.rayon)*echelle,self.rayon*2*echelle,self.rayon*2*echelle],2)
            

4. Exemple

from lennardJones2d import *
import pygame
from matplotlib.pyplot import *
import time
            

On considère 400 sphères. La grille est utilisée pour rechercher les paires en interactions.

Nx = 20
densite = 0.3
rc = 2.5
h = 0.01
sys = Systeme(Nx,densite,rc)
sys.initialiser(1.0)
pygame.init()
taille = 500
screen = pygame.display.set_mode([taille,taille])
echelle = taille*1.0/sys.L
clock = pygame.time.Clock()
done = False
pression = numpy.zeros(0)
ec = numpy.zeros(0)
energie = numpy.zeros(0)
iter = 500
t = time.clock()
while not done and iter > 0:
    iter -= 1
    clock.tick(30)
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            done = True
    sys.dessiner_disques(screen,echelle,(255,0,0))
    sys.integration(h,10)
    sys.calculer_cinetique()
    pression = numpy.append(pression,sys.pression)
    ec = numpy.append(ec,sys.Ecinetique)
    energie = numpy.append(energie,sys.energie)
    pygame.display.flip()
print(time.clock()-t)
pygame.image.save(screen,"../../../../figures/sciphys/dynmol/lennardjones2d/disques1.png")
pygame.quit()

           

print(sys.densite)
--> 0.3

Voici la dernière configuration obtenue, enregistrée sous forme d'image PNG :

Lorsque la fenêtre graphique pygame est fermée, on procède au tracé de l'énergie, de l'énergie cinétique et de la pression instantanées.

figure(figsize=(10,4))
plot(energie)
axis([0,energie.size,-2.0,2.0])
title("Energie")
           

Les fluctuations de l'énergie sont très faibles et sa valeur moyenne reste constante, ce qui montre que le pas de temps est assez petit.

figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Energie cinetique")
           

figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,0,pression.max()])
title("pression")
           

Voici le même calcul, avec utilisation de la liste des voisins pour accélérer le calcul des forces (la grille est aussi utilisée pour construire la liste des voisins).

Nx = 20
densite = 0.3
rc = 2.5
deltaR=0.5
h = 0.01
sys = Systeme(Nx,densite,rc,deltaR)
sys.initialiser(1.0)
pygame.init()
taille = 500
screen = pygame.display.set_mode([taille,taille])
echelle = taille*1.0/sys.L
clock = pygame.time.Clock()
done = False
pression = numpy.zeros(0)
ec = numpy.zeros(0)
energie = numpy.zeros(0)
iter = 500
t = time.clock()
while not done and iter > 0:
    iter -= 1
    clock.tick(30)
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            done = True
    sys.dessiner_disques(screen,echelle,(255,0,0))
    sys.integration_avec_liste_voisins(h,10)
    sys.calculer_cinetique()
    pression = numpy.append(pression,sys.pression)
    ec = numpy.append(ec,sys.Ecinetique)
    energie = numpy.append(energie,sys.energie)
    pygame.display.flip()
print(time.clock()-t)
pygame.image.save(screen,"../../../../figures/sciphys/dynmol/lennardjones2d/disques2.png")
pygame.quit()

           

figure(figsize=(10,4))
plot(energie)
axis([0,energie.size,-2.0,2.0])
title("Energie")
           

figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           

figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,0,pression.max()])
title("pression")
           

L'utilisation de la liste de voisins avec Δr=0.5 permet de gagner environ 30 pourcent en temps de calcul. Si la valeur de Δr est trop faible, la reconstruction de la liste de voisins est trop fréquente. Si la valeur de Δr est trop grande (supérieure à 1), le nombre de sphères dans la liste est trop grand. La plage optimale est entre 0.3 et 1.

Voyons le même calcul avec deux fois plus de sphères (même densité).

Nx = 40
densite = 0.3
rc = 2.5
deltaR=0.5
h = 0.01
sys = Systeme(Nx,densite,rc,deltaR)
sys.initialiser(1.0)
pygame.init()
taille = 500
screen = pygame.display.set_mode([taille,taille])
echelle = taille*1.0/sys.L
clock = pygame.time.Clock()
done = False
pression = numpy.zeros(0)
ec = numpy.zeros(0)
energie = numpy.zeros(0)
iter = 500
t = time.clock()
while not done and iter > 0:
    iter -= 1
    clock.tick(30)
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            done = True
    sys.dessiner_disques(screen,echelle,(255,0,0))
    sys.integration_avec_liste_voisins(h,10)
    sys.calculer_cinetique()
    pression = numpy.append(pression,sys.pression)
    ec = numpy.append(ec,sys.Ecinetique)
    energie = numpy.append(energie,sys.energie)
    pygame.display.flip()
print(time.clock()-t)
pygame.image.save(screen,"../../../../figures/sciphys/dynmol/lennardjones2d/disques3.png")
pygame.quit()

           

figure(figsize=(10,4))
plot(energie)
title("Energie")
axis([0,energie.size,-2.0,2.0])
           

figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           

figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,0,pression.max()])
title("pression")
           

Les fluctuations de pression et de température sont plus faibles (en théorie d'un facteur 2). Pour comparaison, les calculs de L. verlet [1] ont été faits avec 864 sphères. On en a ici 1600.

Le temps de calcul est un peu plus de 4 fois plus grand, ce qui montre que la complexité est proche de N¹, c'est-à-dire que le temps de calcul est proportionnel au nombre de sphères.

Voyons à présent un exemple de calcul avec ajustements des vitesses pour une obtenir une température donnée.

Nx = 40
densite = 0.3
rc = 2.5
deltaR=0.5
h = 0.01
T = 0.5
vinit = math.sqrt(2*T)
pression = numpy.zeros(0)
ec = numpy.zeros(0)
energie = numpy.zeros(0)
sys = Systeme(Nx,densite,rc,deltaR)
sys.initialiser(vinit)
sys.integration_avec_liste_voisins(h,1000)
sys.init_moyennes()
iter = 1000
while iter>0:
    iter -= 1
    sys.integration_avec_liste_voisins(h,10)
    sys.calculer_cinetique()
    pression = numpy.append(pression,sys.pression)
    ec = numpy.append(ec,sys.Ecinetique)
    energie = numpy.append(energie,sys.energie)
    if iter%100==0:
        sys.ajuster_vitesses(T)
            

figure(figsize=(10,4))
plot(energie)
title("Energie")
axis([0,energie.size,-2.0,2.0])
           

figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           

figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,-1,1])
title("pression")
           

Il y a un ajustement de vitesse tous les 1000 pas (une itération fait 10 pas). Il faut environ 5 ajustement pour que la température atteigne la valeur voulue. Lorsque la température est atteinte, on peut arrêter les ajustements et procéder à des calculs statistiques sur les vitesses (distribution, corrélations, etc).