3/ Potentiel de Lennard-Jones

Crystal 2D simulé avec le potentiel de Lennard-Jones. La coloration montre le nombre de voisin proche de chaque particule.

Objectifs du TP

  • Implémenter des interactions entre particules basées sur le potentiel de Lennard-Jones
  • Tester l’implémentation en analysant 1) la conservation de l’énergie, 2) les trajectoires des particules (avec Ovito)
  • Détermination d’une structure cristalline stable
  • Calcul de l’interaction par liste de voisin et validation avec le calcul “naïf”
  • Comparaison du temps d’exécution des deux méthodes de calcul d’interaction (faculatif)

1 Implémentation naïve

Dans cette partie du TP, notre but est de réutiliser le code des interactions développé lors du précédent TP, de l’adapter au potentiel de Lennard-Jones, et de tester l’implémentation.

Afin de généraliser le code du TP précédent, on rappelle la formule pour le calcul des forces et l’énergie totale :

\[\begin{align*} \vec f_i & = \sum_{j=1,\ j\neq i}^N U_2'(d_{ij})\vec n_{ij} \\ E_p & = \sum_{i=1}^N\sum_{j=i+1}^N U_2(d_{ij}) = \frac{1}{2} \sum_{i=1}^N \sum_{j=1,j\neq i}^N U_2(d_{ij}) \end{align*}\]

On remarque que le calcul de la distance \(d_{ij}\) et le calcul du vecteur normal \(\vec n_{ij}\) ne dépendent pas du potentiel d’interaction. De même, une fois les forces par paire calculées, la somme se fait indépendamment du potentiel. On peut donc réutiliser une grande partie du code du TP précédent.

Consigne

Vous pouvez donc copier ou déplacer les fonctions compute_dij_nij, elastic_bounce_force et compute_forces_collision dans un fichier potential.py. Dans ce même fichier, sur le modèle du TP précédent, implémentez les deux fonctions :

import numpy as np
import numpy.typing as npt

def lj_force(dij: npt.NDArray[float], epsilon: float, sigma: float):
    """Calcule l'intensité de la force selon le potentiel de Lennard-Jones"""
    return ...

def lj_energy(dij: npt.NDArray[float], epsilon: float, sigma: float):
    """Calcule pour chaque paire l'énergie potentielle de Lennard-Jones"""
    return ...

def compute_forces_energy_lj(positions: npt.NDArray[float], 
                             epsilon: float, sigma: float):
    """Calcule les forces sur chaque particule"""
    # Étape 1: calcul de dij et nij
    # Étape 2: calcul de -U_2'(dij) avec la fonction lj_force
    # Étape 2': calcul de U_2(dij) avec la fonction lj_energy
    # Étape 3: calcul de fij
    # Étape 4: somme pour calculer fi
    # Étape 4': somme pour calculer l'énergie potentielle totale
    return ..., ... # on retourne forces et énergie potentielle

Test : cristal 2D

Pour tester l’implémentation des fonctions ci-dessus, on procède à la simulation d’un cristal hexagonal 2D. Pour générer la structure cristalline, on utilise la fonction ci-dessous :

def make_crystal(lattice_spacing: float,
                 radius: float,
                 angle: float = np.pi / 2):
    """Crée un cristal avec un angle"""
    l = float(lattice_spacing)

    # Nombre de répétitions de la cellule de base du cristal
    nx = int(2 * radius / l) + 1
    ny = int(2 * radius / (l * np.sin(angle))) + 1

    # Génération de la structure cristaline
    i, j = np.mgrid[0:nx, 0:ny]
    ei = np.array([1., 0, 0])
    ej = np.array([np.cos(angle), np.sin(angle), 0])

    lattice = (
        ei[np.newaxis, np.newaxis] * i[..., np.newaxis]
        + ej[np.newaxis, np.newaxis] * j[..., np.newaxis]
    )
    
    lattice = lattice.reshape(lattice.shape[0] * lattice.shape[1], -1).T

    # On passe d'un parallélogramme à un rectangle
    shifts = np.floor(lattice[0] / nx)
    lattice[0] -= shifts * nx

    # On centre le lattice
    center = np.array([radius / l] * 2 + [0])
    lattice -= center[:, np.newaxis]

    # On masque les positions hors du cercle de rayon "radius"
    mask = sum(x**2 for x in lattice) <= (radius / l)**2

    # On retourne le lattice avec la bonne constante
    return lattice[:, mask] * lattice_spacing
positions = make_crystal(lattice_spacing=2, radius=8, angle=np.pi/2)[:2] # On ne garde que x et y

Figure 1: Cristal carré de paramètre \(a\)

La fonction prend trois paramètres : la constante de la structure périodique (paramètre cristallin ou paramètre de maille), notée \(a\) sur la Figure 1, le rayon (approximatif) du cristal, et l’angle de la cellule périodique (un angle de \(\pi / 2\) donne une grille recangulaire comme ci-dessus).

Consigne

Générez une structure cristalline de relativement petite taille (radius = 10 * lattice_spacing). Avec les paramètres \(\sigma = 1\), \(\varepsilon = 1\) et \(m = 1\), calculez l’évolution du système pour environ 2000 pas de temps, avec vitesse initiale nulle. À chaque pas de temps, affichez (avec print) les valeurs numériques de :

  • L’énergie cinétique totale \(E_c = \sum_{i=1}^N \frac{1}{2}m_i v_i^2\)
  • L’énergie potentielle totale (voir équation ci-dessus), calculée avec compute_forces_energy_lj
  • L’énergie totale \(E_t = E_c + E_p\)

Attention : L’énergie totale, bien que non-constante, ne doit pas augmenter de manière significative. Cherchez une valeur convenable du pas de temps. De même, le paramètre de maille \(a\) doit être choisi en fonction du paramètre \(\sigma\) du potentiel de Lennard-Jones pour que la structure soit stable (c’est-à-dire que malgré les vibrations, la structure cristalline se maintient, comme sur l’animation). Expérimentez différentes valeurs pour obtenir un cristal stable.

En sauvegardant les valeurs à chaque pas de temps des trois énergies ci-dessus dans des listes, affichez avec plt.plot l’évolution des énergies en fonction du temps (n’oubliez pas les unités sur le graphe). Expérimentez avec la valeur de dt pour voir comment celle-ci affecte les énergies du système.

En visualisant les trajectoires des atomes dans Ovito, quel est l’angle de la cellule périodique qui donne une structure stable ?

2 Calcul de l’interaction avec liste de voisins

L’algorithme des listes de voisins permet d’accélérer de manière significative le temps de calcul en réduisant la complexité algorithmique. On ne va pas implémenter d’algorithme de calcul de listes de voisins, à la place, on fera appel à une bibliothèque externe : matscipy (Grigorev et al. 2024), qui est une bibliothèque Python spécialisée dans le traitement de données de simulations en dynamique moléculaire.

Consigne

On commence par installer matscipy. Dans Spyder, mettez le curseur dans la console IPython, normalement située en bas à droite de la fenêtre principale (indiquée avec un onglet “Console 1/A”). Dans la console, tapez l’une des deux commandes suivantes :

!pip install --proxy http://proxyweb.upmc.fr:3128 matscipy   # <- pour les ordis de TP
!pip install matscipy                                        # <- pour votre ordi perso

Pour vérifier que l’installation s’est faite correctement, toujours dans la console, exécutez la commande suivante :

from matscipy.neighbours import neighbour_list

S’il ne se passe rien, c’est que l’installation a réussi.

Contrairement à l’algorithme “naïf” que vous avez implémenté ci-dessus, le calcul des forces avec liste de voisins a besoin d’un paramètre supplémentaire : le rayon de troncations des interactions \(r_\mathrm c\) (cutoff en anglais).

Afin de faciliter l’implémentation (la fonction neighbour_list de matscipy prend des paramètres dans un format spécifique), la fonction ci-dessous implémente le calcul des forces (à l’exception de la fonction lj_force que vous avez normalement déjà implémentée).

from matscipy.neighbours import neighbour_list

def compute_forces_energy_neighbours_lj(positions,
                                        epsilon, sigma, cutoff,
                                        domain=None, periodicity=None):
    """Calcule les forces entre particules par liste de voisins"""
    if domain is None:
        domain = np.max(positions, axis=1) - np.min(positions, axis=1)

    domain = np.asanyarray(domain)

    # On complète la matrice du domaine
    if domain.ndim == 1:
        domain = np.diag(domain)

    if domain.shape == (2, 2):
        extended_domain = np.eye(3)
        extended_domain[:2, :2] = domain
        domain = extended_domain

    # Non-périodique par défaut
    if periodicity is None:
        periodicity = np.array([False] * 3)

    periodicity = np.asanyarray(periodicity)

    # On complète la périodicité
    if periodicity.shape[0] == 2:
        periodicity = np.concatenate((periodicity, [False]))

    # On rajoute une coordonnée en 2d
    if positions.shape[0] == 2:
        full_positions = np.vstack((positions, np.zeros(positions.shape[1])))
    else:
        full_positions = positions

    # Calcul des voisins avec matscipy
    i, j, dij, rij = neighbour_list('ijdD', positions=full_positions.T,
                                    cell=domain, pbc=periodicity,
                                    cutoff=float(cutoff))

    # Calcul de la force
    fij = lj_force(dij, epsilon, sigma)
    fij = (fij / dij)[np.newaxis] * rij.T

    # Calcul de l'énergie
    eij = lj_energy(dij, epsilon, sigma) - lj_energy(cutoff, epsilon, sigma)

    # Somme des forces de paires pour chaque particule
    forces = np.zeros_like(positions)
    natoms = forces.shape[1]
    for d in range(positions.shape[0]):
        forces[d] += 0.5 * np.bincount(j, weights=fij[d], minlength=natoms)
        forces[d] -= 0.5 * np.bincount(i, weights=fij[d], minlength=natoms)
    return forces, 0.5 * eij.sum()
Consigne

Utilisez la fonction compute_forces_energy_neighbours_lj à la place de compute_forces_energy_lj, choisissez une valeur pour \(r_\mathrm c\), le rayon de troncation, et vérifiez que le résultat de simulation est (quasiment) identique.

Pour cela, on pourra observer la trajectoire des atomes dans Ovito, mais aussi se fier aux valeurs de l’énergie cinétique et de l’énergie potentielle. On s’attend à voir une (légère) différence dans l’une de ces deux énergies : expliquez pourquoi.

Coût de calcul (facultatif)

La performance est un aspect important de la simulation : certains phénomènes physiques ne peuvent être observés qu’à grande échelle, dans notre cas quand le nombre de particules est suffisamment grand. Il est donc primordial de contrôler le coût de calcul en fonction de la “taille” de la simulation.

On cherche donc à déterminer le temps de calcul en fonction du nombre de particules. Pour calculer le temps d’exécution d’une fonction, on peut se servir du module time de Python, et en particulier de la fonction perf_counter(). Par exemple (à essayer dans un script séparé) :

import time

def heavy_compute():
    # Cette fonction fait des calculs très coûteux que l'on veut mesurer
    time.sleep(1.25) # <- Simule des calculs très coûteux, à remplacer

start = time.perf_counter() # Temps à l'horloge avant le calcul
heavy_compute()
end = time.perf_counter()   # Temps à l'horloge après le calcul

total_time = end - start    # Temps écoulé pendant le calcul
print(f"Il s'est écoulé {total_time} secondes.")
Il s'est écoulé 1.2503020390000188 secondes.
Consigne

Mettez le code de simulation du cristal, avec compute_forces_energy_lj (calcul naïf), dans la fonction suivante :

def time_lennard_jones(radius: float, nsteps: int):
    """Calcul le temps d'exécution d'une simulation"""

    # Générer un cristal et le simuler pour nsteps
    # On utilisera la fonction compute_forces_energy_lj()
    # dans un premier temps

    total_time = ...  # Temps d'exécution à calculer
    natoms = ...      # Nombre total d'atomes
    return natoms, total_time

Pour différents rayons du cristal (de \(R = 5a\) à \(R = 15a\)) et nsteps=2000 calculer le temps d’exécution. Sur un graphe en échelle logarithmique, affichez le temps d’exécution en fonction du nombre de particules. À quelle pente de la courbe doit-on s’attendre ? Est-ce bien la pente observée ?

On souhaite maintenant comparer le coût de calcul de l’implémentation par liste de voisins au coût pour l’implémentation naïve.

Consigne

Faire la même analyse que pour le calcul naïf, mais cette fois avec la fonction compute_forces_energy_neighbours_lj, et afficher les temps de calcul des deux implémentations sur le même graphe en échelle logarithmique. Commentez la complexité algorithmique des deux méthodes à partir du graphe.

Les références

Grigorev, Petr, Lucas Frérot, Fraser Birks, Adrien Gola, Jacek Golebiowski, Jan Grießer, Johannes L. Hörmann, et al. 2024. « Matscipy: Materials Science at the Atomic Scale with Python ». Journal of Open Source Software 9 (93): 5668. https://doi.org/10.21105/joss.05668.