Table des matières Python

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(r)=4ε[(dr)12-(dr)6](1)

La force intermoléculaire correspondante est

f(r)=48εd[(dr)13-12(dr)7]ur(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()
                  
figAfigA.pdf

À courte distance (r<21/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 rc au delà de laquelle on considère que la force est nulle. On peut prendre par exemple rc=3d.

2.b. Équations du mouvement

Les équations du mouvement des sphères sont obtenues dans le cadre de la mécanique classique (loi de Newton). Le système comporte N sphères de masse m. On note ri la position de la sphère d'indice i (de son centre). Le système d'équations différentielles s'écrit :

md2ridt2=jif(ri-ri)(3)

où l'indice i varie de 0 à N-1. Pour obtenir une solution de ces équations, il faut fixer les positions et les vitesses initiales des sphères.

L'énergie du système s'écrit :

E=12mivi2+4εi<j(d/rij)12-(d/rij)6(4)

rij désigne la distance entre les sphères i et j. L'énergie du système est constante.

En pratique, les équations sont considérées sous leur forme adimensionnelle, obtenue en posant m=1, ε=1 et d=1.

2.c. Intégration des équations différentielles

L'intégration des équations du mouvement est faite avec la méthode de Stormer-Verlet :

vi(t+h/2)=vi(t)+(h/2)ai(t)(5)ri(t+h)=ri(t)+hvi(t+h/2)(6)ai(t+h)=jif(ri(t+h)-rj(t+h))(7)vi(t+h)=vi(t+h/2)+(h/2)ai(t+h)(8)

Le pas de temps h est fixe.

2.d. Subdivision du domaine

Les sphères de déplacent sur un domaine carré de largeur L. Pour simuler un système de taille infinie, on applique des conditions limites périodiques : une sphère quittant le domaine par un côté réapparaît par le côté opposé.

Pour calculer la force sur une sphère, il faut déterminer rapidement les sphères dont la distance vérifie r<rc, que l'on appelera les sphères voisines. On utilise pour cela une subdivision du domaine carré sous forme de grille. Les cellules sont carrées, avec une largeur lc>rc. Chaque sphère est associée à une cellule de la grille, celle qui contient son centre. Pour une sphère donnée, les interactions ne sont possibles qu'avec les sphères situées soit dans la même cellule, soit dans une des 8 cellules voisines. En effet, les sphères situées en dehors de ces 9 cellules sont à une distance supérieure à rc.

La figure suivante représente la grille, avec des sphères dont on a représenté le rayon d'action sous forme d'un cercle en pointillé de rayon rc (la sphère elle-même est représentée par un cercle de rayon d). Une sphère de centre C est représentée, avec les 9 cellules dans lesquelles il faut rechercher d'éventuelles voisines en interaction.

figureA.svgFigure pleine page

On remarque que pour une sphère située dans une cellule du bord du domaine, des cellules situées près du bord opposé font partie du voisinage à explorer. Une telle cellule est montrée sur la figure (en vert), avec ses 8 proches voisines. Les trois cellules voisines de droite ont un indice i=0, à cause des conditions limites périodiques. Pour ces cellules, les coordonnées des centres des sphères qu'elles contiennent doivent être décalées de la manière suivante :

xd=x+L(9)

2.e. Liste de voisins

La méthode précédente permet d'accélérer considérablement la recherche des voisins (comparée à un parcours systématique de toute les paires). Cependant, la proportion de sphères présentes dans les 9 cellules et dont la distance est effectivement inférieure à rc est en moyenne 0.35. À trois dimensions, cette proportion tombe à 0.16.

Une autre méthode, utilisée par Verlet ([1]), consiste à utiliser une liste de paires de sphères voisines. Cette liste doit être régulièrement générée. Pour éviter d'avoir à le faire trop souvent, on considère une distance de voisinage supérieure à la distance de coupure :

rv=rc+Δr(10)

Δr est choisi de telle sorte que plusieurs pas de calculs (chacun de durée h) puissent se faire sans que les sphères initialement à une distance supérieure à rv n'entrent dans la zone r<rc. On choisit par exemple Δr=0.3 pour rc=2.3. Pour savoir à quel moment la liste de voisins doit être regénérée, on calcule à chaque pas de temps la vitesse maximale vm des sphères et on incrémente une variable deplac, qui représente la borne supérieure du déplacement des sphères, de la manière suivante :

deplacdeplac+2vmh(11)

Lorsque la variable deplac dépasse Δr, la liste de voisins est entièrement regénérée.

Comme la liste de voisins doit être générée assez fréquemment (typiquement tous les 10 pas), il est bon d'utiliser la subdivision en grille pour effectuer la recherche des sphères dont la distance est inférieure à rv ([2]). La largeur des cellules doit pour cela vérifier lc>rv.

2.f. Température et pression

On rappelle ici les résultats démontrés dans modèle des sphères molles à deux dimensions.

L'énergie cinétique moyenne du système est :

<Ec>=NkBT(12)

Les calculs seront fait avec kB=1. La température est alors l'énergie cinétique moyenne par sphère.

La pression multipliée par l'aire A=L2 est obtenue avec le viriel :

PA=NkBT+24ε<i<j((d/rij)12-12(d/rij)6)>(13)

Les fluctuations de pression et de température sont en principe d'autant plus faibles que le nombre de sphères est grand. Lorsqu'on cherche à obtenir l'équation d'état, on souhaite généralement procéder à une température T0 donnée a priori. On peut pour cela choisir la vitesse initiale suivante :

vi=2T0(14)

Dans le cas des sphères molles, cela est suffisant pour obtenir la température voulue car l'énergie potentielle est négligeable. Dans le cas présent, la présence de forces attractives conduit à une variation très importante (et difficile à prévoir) de la température au début de la simulation.

Pour obtenir la température souhaitée, il faut procéder de temps en temps à des ajustements de vitesse ([2]). Lorsque l'équilibre est atteint, on calcule la température moyennes Tm sur un grand nombre de pas de temps (par exemple 1000). On en déduit le facteur multiplicatif à appliquer aux vitesses :

f=T0Tm(15)

Cette correction de vitesse doit être appliquée plusieurs fois avant que la valeur moyenne de la température atteigne la valeur souhaitée (à une tolérance près).

3. Implémentation Python

lennardJones2d.py
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 à rc+Δ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 à rv=rc+Δ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 :

disques

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.

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

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()

           
disques
figure(figsize=(10,4))
plot(energie)
axis([0,energie.size,-2.0,2.0])
title("Energie")
           
figEfigE.pdf
figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           
figFfigF.pdf
figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,0,pression.max()])
title("pression")
           
figGfigG.pdf

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()

           
disques
figure(figsize=(10,4))
plot(energie)
title("Energie")
axis([0,energie.size,-2.0,2.0])
           
figHfigH.pdf
figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           
figIfigI.pdf
figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,0,pression.max()])
title("pression")
           
figJfigJ.pdf

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 N1, 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])
           
figKfigK.pdf
figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           
figLfigL.pdf
figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,-1,1])
title("pression")
           
figMfigM.pdf

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).

Références
[1]  L. Verlet,  Computer experiments on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules,  (Phys Rev 159 p 98, 1967)
[2]  D.C. Rapaport,  The art of molecular dynamics simulation,  (Cambridge University Press, 2004)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.