Modèle des sphères dures à deux dimensions

1. Introduction

Ce document présente une application de l'algorithme de Monte-Carlo Metropolis au modèle des sphères dures à deux dimensions. Metropolis, Rosenbluth et al. [1] ont effectué cette simulation afin d'obtenir l'équation d'état d'un gaz bidimensionnel de sphères dures.

2. Modèle des sphères dures

Le modèle des sphères dures consiste à représenter les atomes d'un fluide (monoatomique) par des sphères de diamètre d. Ce modèle est utilisé en physique statistique pour obtenir une équation d'état ([2][3]).

Lorsque les sphères ne sont pas en contact, elles n'exercent aucune force l'une sur l'autre. Au contact, une force répulsive infinie apparaît, qui empêche l'interpénétration des sphères. L'énergie potentielle d'interaction entre deux sphères a donc la forme suivante :

La distribution statistique des vitesses (distribution de Maxwell) est indépendante de la taille des sphères. Dans une simulation de Monte-Carlo, on peut donc se limiter à considérer les configurations du système dans l'espace et pas les vitesses des atomes. Autrement dit, l'énergie cinétique n'est pas prise en compte et l'énergie E à considérer ne comporte que l'énergie potentielle.

Puisque l'énergie potentielle est nulle (sauf pendant les chocs dont le durée est nulle), l'énergie interne du fluide de sphères dures se réduit à l'énergie cinétique des atomes (comme pour un gaz parfait) :

où N est le nombre de sphères.

L'équation d'état d'un fluide de sphères dures, qui relie la pression à la densité, est en revanche différente de celle du gaz parfait. L'équation d'état s'obtient en utilisant la fonction de distribution radiale g(r), qui représente la probabilité pour que deux sphères (leur centre) soient distantes de r. Plus précisément, si l'on considère qu'une sphère se trouve centrée à l'origine, le nombre moyen de sphères centrées entre les rayons r et r+dr est

où V est le volume du système. La fonction de distribution radiale est non seulement un outil de calcul théorique, mais aussi une fonction accessible expérimentalement par diffraction de rayonnement (rayons X, neutrons). Elle apporte une information sur la structure locale du fluide (à l'échelle atomique), dont dépend son équation d'état. Pour un fluide de sphères dures, la théorie statistique ([3]) conduit à l'équation d'état suivante :

Il faut donc déterminer g(d), la valeur de la fonction de distribution radiale pour une distance égale au diamètre des sphères, c'est-à-dire lorsque les sphères sont en contact. Il s'agit d'une valeur moyenne qui peut être obtenue par une simulation de Monte-Carlo, en suivant l'algorithme de Metropolis. En introduisant n le nombre moyen de sphères par unité de volume en contact avec une sphère, l'expression (5) devient :

L'équation d'état peut être développée par rapport à la densité N/V, ce qui donne le développement du viriel :

Le calcul de n en fonction de la densité doit donc permettre de calculer les coefficients du viriel B₂,B₃,....

On se limite au modèle des sphères dures sur un plan (système bidimensionnel). Si A est l'aire du domaine carré, l'équation d'état s'écrit ([1]) :

où n est la densité surfacique moyenne des sphères en contact avec une sphère.

3. Algorithme

4. 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 :

class Disque:
    def __init__(self,x,y):
        self.x = x
        self.y = y
            

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

La classe suivante représente une cellule, qui contient une liste de disques. Le type set (ensemble) est utilisé pour cette liste. Ce type permet d'enlever et d'ajouter facilement un objet. 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 = set()
        self.decal_x = 0.0
        self.decal_y = 0.0
    def ajouter(self,indice):
        self.ensemble_disques.add(indice)
    def enlever(self,indice):
        self.ensemble_disques.remove(indice)
             

La classe Grille contient un tableau bidimensionnel de cellules (une liste de listes). Le constructeur, qui prend en argument le nombre de cellules Nc par dimension et la largeur L de la grille, crée le tableau. La fonction obtenir_cellule(i,j) renvoie 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 disques par dimension et la densité. Une liste est créée pour stocker les disques. Les disques seront référencés par l'indice dans cette liste.

class Systeme:
    def __init__(self,Nx,densite):
        self.N = Nx*Nx
        self.Nx = Nx
        self.rayon = 0.5
        self.diam = self.rayon*2
        self.diam2 = self.diam**2
        self.densite = densite
        self.L = math.sqrt(math.pi/densite)*0.5*Nx
        self.Nc = int(self.L/self.diam)
        self.lc = self.L/self.Nc
        self.deplac = (self.L/Nx-self.diam)*1.0
        self.grille = Grille(self.Nc,self.L)
        self.liste_disques = []
        for d in range(self.N):
            self.liste_disques.append(Disque(0,0))
        self.accept = 0
              

La fonction suivante place un disque, dont l'indice est donné, à une position donnée. Les coordonnées du disque sont modifiées et le disque (son indice) est ajouté dans la cellule correspondante.

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

              

La fonction suivantes initialise le système.

    def initialiser(self):
        dx = self.L*1.0/(self.Nx)
        if dx < self.diam:
            raise Exception("densite trop forte")
        else:
            dy = dx
            x = dx/2
            y = x
            for k in range(self.N):
                self.init_disque(k,x,y)
                x += dx
                if x  > self.L:
                    x = dx/2
                    y += dy

              

La fonction suivante déplace un disque à une position donnée. Le déplacement nécessite éventuellement d'enlever le disque de sa cellule pour le placer dans une nouvelle cellule. Si les coordonnées sortent du domaine, les conditions limites périodiques sont appliquées.

    def deplacer_disque(self,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
        disque = self.liste_disques[indice]
        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 teste si les coordonnées (x,y) d'un disque conduisent à une intersection avec les disques d'une cellule donnée par ces indices (i,j). L'indice du disque est fourni pour éviter de tester le disque avec lui-même.

    def tester_intersection(self,indice,i,j,x,y):
        cellule1 = self.grille.obtenir_cellule(i,j)
        if not cellule1:
            return False
        else:
            for indice1 in cellule1.ensemble_disques:
                if indice!=indice1:
                    disque1 = self.liste_disques[indice1]
                    dx = x-disque1.x-cellule1.decal_x
                    dy = y-disque1.y-cellule1.decal_y
                    if dx*dx+dy*dy < self.diam2:
                        return True
        return False
              

La fonction suivante teste si les coordonnées (x,y) d'un disque conduisent à une intersection avec un des disques du voisinage. Les coordonnées peuvent être hors du domaine, auquel cas les conditions limites périodiques sont appliquées. L'indice du disque est aussi donné, pour être transmis à la fonction tester_intersection.

    def tester_toutes_intersections(self,indice,x,y):
        if x<0:
            x += self.L
        elif x>=self.L:
            x -= self.L
        if y<0:
            y += self.L
        elif y>=self.L:
            y -= self.L
        i = int(x/self.lc)
        j = int(y/self.lc)
        if self.tester_intersection(indice,i,j,x,y):
            return True
        if self.tester_intersection(indice,i-1,j,x,y):
            return True
        if self.tester_intersection(indice,i+1,j,x,y):
            return True
        if self.tester_intersection(indice,i,j-1,x,y):
            return True
        if self.tester_intersection(indice,i,j+1,x,y):
            return True
        if self.tester_intersection(indice,i-1,j-1,x,y):
            return True
        if self.tester_intersection(indice,i-1,j+1,x,y):
            return True
        if self.tester_intersection(indice,i+1,j-1,x,y):
            return True
        if self.tester_intersection(indice,i+1,j+1,x,y):
            return True
        return False
               

La fonction suivante effectue le pas élémentaire de l'algorithme de Metropolis, qui consiste à prendre un disque au hasard, à lui attribuer un déplacement aléatoire et à tester la possibilité de ce déplacement. Si le déplacement est permis, le disque est effectivement déplacé.

    def metropolis(self):
        indice = random.randrange(self.N)
        disque = self.liste_disques[indice]
        x = disque.x
        y = disque.y
        dx = self.deplac*(random.random()-0.5)
        dy = self.deplac*(random.random()-0.5)
        if not self.tester_toutes_intersections(indice,x+dx,y+dy):
            self.deplacer_disque(indice,x+dx,y+dy)
            self.accept += 1

              

La fonction suivante effectue une boucle comportant n pas élémentaires, et renvoie le taux d'acceptation des déplacements. Si ce taux est trop faible, il faudra diminuer le déplacement maximal des disques.

    def boucle_metropolis(self,n):
        self.accept = 0
        for k in range(n):
            self.metropolis()
        return(self.accept*1.0/n)

              

Pour déboguer, il peut être utile de suivre le déplacement d'un disque. La fonction suivante renvoie les coordonnées d'un disque pour n pas de Metropolis :

    def suivi_disque(self,n,indice):
        x = []
        y = []
        disque = self.liste_disques[indice]
        for k in range(n):
            self.metropolis()
            x.append(disque.x)
            y.append(disque.y)
        return(x,y)

             

La représentation graphique des disques sera faite avec la bibliothèque Pygame, qui permettra de visualiser les changements de configuration sous forme d'animation. La fonction suivante trace les disques, avec un facteur d'échelle :

    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)             
             

Voici finalement la fonction de calcul de la distribution radiale, qui devra être appelée une fois que le système sera mis à l'équilibre par un nombre suffisant de déplacements. Les arguments sont : la distance maximale K (un entier), le nombre d'intervalles de l'histogramme p, le nombre de blocs de calculs n₁ et le nombre de pas de Metropolis dans chaque bloc n₂.

    def calculer_distribution_radiale(self,K,p,n1,n2):
        K=int(K)
        delta_r2 = (K**2-1)*self.diam2/p
        imax = int(K*self.diam)
        g = numpy.zeros(p,dtype=numpy.float64)
        r = numpy.zeros(p,dtype=numpy.float64)
        self.accept = 0
        for k in range(n1):
            for l in range(n2):
                self.metropolis()
            for i in range(self.Nc):
                for j in range(self.Nc):
                    cellule = self.grille.obtenir_cellule(i,j)
                    for i1 in range(i-imax,i+imax+1):
                        for j1 in range(j-imax,j+imax+1):
                            cellule1 = self.grille.obtenir_cellule(i1,j1)
                            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
                                        distance2 = dx*dx+dy*dy
                                        k = int((distance2-self.diam2)/delta_r2)
                                        if k<p:
                                            g[k] += 1
        delta_A = math.pi*delta_r2
        b = delta_A*self.N**2/self.L**2*n1
        for k in range(p):
            r[k] = math.sqrt(self.diam2+(k+0.5)*delta_r2)
            g[k] = g[k]/b
        taux_accept = self.accept*1.0/(n1*n2)
        print("taux d'acceptation : %f\n"%taux_accept)
        return (r,g)
                                    

               

5. Exemple

from spheres2dMetropolis import *
import pygame
from matplotlib.pyplot import *
            

On définit un système comportant 49 sphères. On effectue pour commencer 100000 pas de Metropolis pour mettre le système à l'équilibre.

N = 7
densite=0.5
systeme = Systeme(N,densite)
systeme.initialiser()
systeme.boucle_metropolis(100000)
            

Pour contrôler le bon fonctionnement du programme, on effectue une animation qui montre les déplacements des disques. La boucle de calcul prend fin lorsqu'on ferme la fenêtre graphique.

pygame.init()
taille = 500
screen = pygame.display.set_mode([taille,taille])
echelle = taille*1.0/systeme.L
clock = pygame.time.Clock()
fin = False
while not fin:
    clock.tick(30)
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            fin = True
    systeme.dessiner_disques(screen,echelle,(255,0,0))
    systeme.metropolis()
    pygame.display.flip()
pygame.image.save(screen,"../../../../figures/sciphys/physistat/spheres2d/disques.png")
pygame.quit()
             

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

Voyons le calcul de la fonction de distribution radiale. On choisit une distance maximale égale à 4 diamètres de disques. L'histogramme comporte 50 points. On effectuer 10000 blocs de 100 pas élémentaires. Pour chaque bloc, l'histogramme est complété.

figure(figsize=(8,6))
(r,g) = systeme.calculer_distribution_radiale(4,50,1000,100)
plot(r,g,'o')
xlabel('r')
ylabel('g')
title('densite = %f'%densite)
axis([0,4,0,g.max()])
              

La plus grande valeur est obtenue pour r=d=1. Pour cette densité relativement élevée, on remarque que la fonction de distribution tend vers sa valeur limite (égale à 1) avec des oscillations. Le maximum de g pour r=d correspond, pour une sphère donnée, à la couche des sphères en contact avec elle. Le maximum vers r=2.2d correspond à la seconde couche, etc. On voit ainsi qu'il y a un ordre local jusqu'à r=4d environ. Au delà de cette distance, il n'y a plus aucune corrélation des positions et la fonction de distribution radiale est égale à 1.

Faisons le même calcul pour une densité plus faible.

densite = 0.2
systeme = Systeme(N,densite)
systeme.initialiser()
systeme.boucle_metropolis(100000)
figure(figsize=(8,6))
(r,g) = systeme.calculer_distribution_radiale(4,50,10000,100)
plot(r,g,'o')
xlabel('r')
ylabel('g')
title('densite = %f'%densite)
axis([0,4,0,g.max()])
              

Pour cette densité, la décroissance de la fonction de de distribution est monotone. La limite est atteinte à r=2.

Ces calculs durent quelques minutes. Pour comparaison, les calculs effectués par Metropolis et al. en 1953 ([1]) duraient de 4 à 5 heures et se limitaient à une zone proche de r=d (dans le but d'obtenir g(d)).

Pour obtenir g(d) (afin de calculer la pression), on peut se contenter de prendre la permière valeur du tableau G. On remarque toutefois que ce point correspond en fait au rayon r²=d²+Δr²/2. Compte tenu de la variation rapide de g(r) au voisinage de r=d, on peut s'attendre à une forte influence de Δr² sur la première valeur de G. Une solution serait de réduire Δr² pour s'approcher au mieux de r=d, mais cette méthode conduirait à une augmentation de la variance, car le nombre de sphères situées dans la première couche se réduirait fortement. Une meilleure solution, utilisée dans [1], consiste à calculer g(r) au voisinage de d et à effectuer une extrapolation pour obtenir la valeur g(d).

Le calcul de g(d) à différentes densités est fait dans Modèle des sphères dures : équation d'état.