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 potentielle8 3/ Potentiel de Lennard-Jones
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)
8.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 :
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_spacingpositions = make_crystal(lattice_spacing=2, radius=8, angle=np.pi/2)[:2] # On ne garde que x et y
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 8.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 ?
8.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 persoPour vérifier que l’installation s’est faite correctement, toujours dans la console, exécutez la commande suivante :
from matscipy.neighbours import neighbour_listS’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.250301496999981 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_timePour 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.
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.