Table des matières Python

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

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

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 B2,B3,....

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]) :

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

3. Algorithme

3.a. Déplacement des sphères

On considère N sphères (ou disques) de diamètre d pouvant se déplacer sur un carré de côté L. Une condition limite périodique est utilisée sur les bords du carré, de manière à simuler un système de taille infinie.

L'algorithme de Metropolis consiste à générer une chaîne de Markov de configurations des sphères. Chaque configuration est générée à partir de la précédente de la manière suivante :

  • Choix aléatoire d'une sphère.
  • Choix d'un déplacement aléatoire pour cette sphère. La nouvelle position est choisie dans un carré centré sur la sphère.
  • Si la nouvelle position de la sphère ne comporte pas d'intersection avec une autre sphère, la sphère est déplacée.
  • Si la nouvelle position de la sphère comporte une intersection avec une autre sphère, le déplacement est refusé. Dans ce cas, la nouvelle configuration est identique à la précédente.

La figure suivante montre un déplacement accepté, un déplacement rejeté (en pointillé), et un déplacement accepté avec franchissement d'un bord du carré. La sphère qui franchit le bord droit se retrouve près du bord gauche.

figureA.svgFigure pleine page

La tentative de déplacement d'une sphère se fait en tirant deux nombres aléatoires et compris entre -1 et 1 :

Cela revient à choisir une nouvelle position aléatoirement dans un carré de côté 2a centré sur l'ancienne position. Le paramètre a sera ajusté en fonction de la densité des sphères, de manière à obtenir un taux de rejet de l'ordre de 50 pour cent. Si la densité est forte, il faudra abaisser ce paramètre, sans quoi trop de mouvements seront rejetés et le système évoluera trop lentement.

3.b. Densité et condition initiale

On définit la densité comme le rapport de l'aire occupée par les sphères sur l'aire du carré :

Cette densité est proportionnelle à N/L2, le nombre de sphères par unité de surface.

Initialement, les sphères sont rangées sur une grille régulière de pas :

3.c. Détection des intersections et subdivision d'espace

À chaque pas de l'algorithme de Metropolis, il faut déterminer si la nouvelle position de la sphère présente une intersection avec une autre sphère. Pour savoir si deux sphères ont une intersection, il suffit de calculer la distance entre leurs centres et de la comparer au diamètre d. En pratique, on évite de calculer une racine carrée : on compare le carré de la distance avec le carré du diamètre.

Pour chaque tentative de déplacement, il y a priori N-1 tests à faire, avec pour chacun le calcul du carré de la distance. Lorsque N est grand, il faut réduire le nombre de tests en les limitant aux sphères voisines. Cela peut être obtenu en utilisant une subdivision du domaine carré sous forme de grille, comme indiqué sur la figure suivante. La largeur des cellules de la grille doit être supérieure au diamètre des sphères. Le diamètre des sphères étant d, la largeur des cellules est lc>d. On note Nc le nombre de cellules par dimension. La largeur du domaine est L=Nclc.

figureC.svgFigure pleine page

Chaque sphère est associée à une cellule de la grille, celle qui contient son centre. Si (x,y) sont les coordonnées du centre de la sphère, les indices de la cellule sont :

E(x) est le plus grand entier inférieur ou égal à x.

À chaque cellule on associe la liste des sphères dont elle contient le centre. Lorsqu'on teste une nouvelle position (x',y'), on commence par calculer les indices (i',j') de sa cellule, puis on effectue les tests d'intersection avec les sphères situées dans cette cellule et dans les 8 cellules voisines. La figure montre un point C à tester avec les 9 cellules concernées. Les distances avec les trois sphères tracées en bleu sont testées.

Lorsque les indices (i',j') de la cellule contenant le point C sont calculés, il est aisé d'obtenir les indices des cellules voisines. Il y a cependant un cas particulier qui nécessite un traitement spécial, celui d'une cellule située près d'un bord. 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 :

3.d. Fonction de distribution radiale

Pour une sphère donnée, le nombre moyen de sphères dont la distance est comprise entre r et r+dr est

g(r) est la fonction de distribution radiale. Pour l'obtenir numériquement, on calcule un histogramme. On choisit tout d'abord une distance maximale

On calcul la taille du voisinage qui permet d'englober ce rayon maximal. Puisque la taille des cellules est lc, il suffit de calculer im=E(rm/lc). Si (i,j) sont les indices de la cellule contenant une sphère, les cellules du voisinage, celles dont les sphères sont à une distance inférieure à rm, sont définies par :

L'intervalle [d2,rm2] est subdivisé en p intervalles égaux de largeur Δr2. Cela donne p couches d'aire πΔr2. Un tableau G à p éléments est utilisé pour stocker l'histogramme, avec des valeurs initiales nulles. À chaque nouvelle configuration du système, on complète l'histogramme de la manière suivante :

  • Pour chacune des N sphères, on considère les sphères situées dans la même cellule et dans les cellules du voisinage.
  • Pour chaque sphère voisine, on calcule la distance au carré r2. Si celle-ci est dans l'intervalle [d2,rm2], on ajoute une unité à l'élément correspondant du tableau G.

La figure suivante montre un exemple pour K=2 avec p=4. Pour chaque sphère, il faut explorer un voisinage constitué de 25 cellules (car im=2).

figureD.svgFigure pleine page

Lorsque n configurations ont été explorées, on divise chaque élément du tableau G par n et par N pour obtenir les valeurs moyennes. On divise aussi par

pour obtenir la fonction g(r).

En principe, ce calcul satistique doit être fait en utilisant toutes les configurations obtenues par l'algorithme de Metropolis. En particulier, deux configurations successives identiques pour cause de rejet d'un déplacement doivent être comptées deux fois dans le calcul.

En pratique, le calcul de l'histogramme consomme beaucoup d'opérations (une boucle sur toutes les sphères avec exploration des cellules voisines pour chacune). On obtient plus rapidement un résultat statisfaisant (une faible variance) en complétant l'histogramme sur quelques configurations seulement, suffisamment espacées. On peut effectuer n2 blocs comportant chacun n1 pas élémentaires (déplacement au maximum d'une sphère). Pour chaque bloc, l'histogramme est complété. À la fin du calcul, on divise par n1 pour obtenir les nombres moyens. Un nombre n2 égal au nombre de sphères donne un bon résultat.

4. Implémentation Python

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

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) 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 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 renvoit 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 renvoit 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 n1 et le nombre de pas de Metropolis dans chaque bloc n2.

    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 :

disques

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()])
              
figA.svgFigure pleine page

La plus grande valeur est obtenue pour r=d=1. Pour cette densité relativement élevée, on remarque que la fonction de distibution 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()])
              
figB.svgFigure pleine page

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 r2=d2+Δr2/2. Compte tenu de la variation rapide de g(r) au voisinage de r=d, on peut s'attendre à une forte influence de Δr2 sur la première valeur de G. Une solution serait de réduire Δr2 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.

Références
[1]  N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller,  Equation of state calculations by fast computing machines,  (J. Chem. Phys. vol. 21 No. 6, 1953)
[2]  J.P. Hansen, I.R. McDonald,  Theory of simple liquids,  (Academic Press, 1990)
[3]  J.L. Bretonnet,  Théorie microscopique des liquides,  (Ellipses, 2010)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.