6 1/ Balistique de \(N\) corps

Dans ce TP, on mettra en place un code pour prédire la trajectoire dans un champs gravitationel de N particules ayant des vitesses initiales aléatoires.

6.1 Mise en place

On considère un système à deux dimensions de \(N\) particules dont on note les positions \(\{\vec{r}_i\}_{i=1}^N\) et les vitesses \(\{\vec{v}_i\}_{i=1}^N\).

On utilise des arrays NumPy pour représenter ces quantitiés

import numpy as np
import numpy.typing as npt
import matplotlib.pyplot as plt

N = 100
D = 2

positions = np.zeros([D, N])
velocities = np.zeros_like(positions)

Ainsi, la première dimension de l’array positions permet de choisir la composante (\(x\) ou \(y\)) de la position de toutes les particules (positions[0] donne toutes les positions selon \(x\)), alors que la deuxième dimension permet de sélectionner le vecteur position d’une seule particule (positions[:, 9] donne le vecteur position de la 10ᵉ particule), voir Figure 6.1.

Figure 6.1: Disposition des coordonnées de particules dans le tableau `positions`

On devra aussi créer un array pour les masses des particules et les forces que subissent chaque particule.

masses = np.ones([N])
forces = np.zeros_like(positions)

Il est important de noter que le tableau masses est unidimensionnel, puisque la masse de chaque particule n’est pas un vecteur mais un scalaire.

Rappel

Avec Numpy, on évite autant que l’on peut de faire des calculs en utilisant des boucles Python. Par exemple, on évitera :

for i in range(N):
    for d in range(D):
        positions[d, i] += velocities[d, i] * dt

À la place, on écrira simplement :

positions += velocities * dt

Cependant, certains arrays Numpy ont parfois des tailles incompatibles. Par exemple, si on veut appliquer une translation uniforme à toutes les particules :

u = np.array([1., 1.])
positions += u

On obtient l’erreur suivante :

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[15], line 1
----> 1 positions += u

ValueError: operands could not be broadcast together with shapes (2,100) (2,) (2,100)

Pour corriger cela sans écrire de boucle, on doit utiliser une fonctionnalité de Numpy appelée broadcasting (voir cette explication en français). Ici, on a le souci que le vecteur de translation n’a qu’un axe, et que positions en a deux. On peut corriger cela avec np.newaxis pour ajouter un axe au vecteur u:

positions += u[:, np.newaxis]

Numpy comprend alors automatiquement que u doit être rajouté à toutes les colonnes de positions.

6.2 Velocity-Verlet, implémentation

On donne les prototypes des fonctions principales de l’algorithme Velocity-Verlet ci-dessous.

def verlet_predictor(positions: npt.NDArray[float],
                     velocities: npt.NDArray[float],
                     masses: npt.NDArray[float],
                     forces: npt.NDArray[float],
                     dt: float):
    """Étape de prédiction de Velocity-Verlet"""

def verlet_corrector(positions: npt.NDArray[float],
                     velocities: npt.NDArray[float],
                     masses: npt.NDArray[float],
                     forces: npt.NDArray[float],
                     dt: float):
    """Étape de correction de Velocity-Verlet"""

Consigne

Implémentez les étapes de prédiction et corrections de l’algorithme Velocity-Verlet. Pour rappel, l’algorithme de Velocity-Verlet s’écrit :

\(\vec v_{i+1/2} = \vec v_i + \frac{\Delta T}{2} \left( \frac{\vec f_i}{m} \right)\)
\(\vec r_{i+1} = \vec r_i + \Delta T \vec v_{i+1/2}\)
Calcul des forces (ici les forces sont constantes)
\(\vec v_{i+1} = \vec v_{i+1/2} + \frac{\Delta T}{2}\left(\frac{\vec f_{i+1}}{m} \right)\)

Les étapes 1 et 2 sont les étapes de prédiction, l’étape 4 est l’étape de correction. La Figure 6.2 ci-dessous résume les étapes de calcul.

En règle générale, on évite de stocker en mémoire les valeurs de la vitesse et de la position à chaque pas de temps. On peut donc écrire l’algorithme avec des opérations qui modifient les vecteurs:

\(\vec v \gets \vec v + \frac{\Delta T}{2} \left( \frac{\vec f}{m} \right)\)
\(\vec r \gets \vec r + \Delta T \vec v\)
\(\vec f \gets \vec f(\vec r)\)
\(\vec v \gets \vec v + \frac{\Delta T}{2}\left(\frac{\vec f}{m} \right)\)

La flèche \(\gets\) symbolise ici le fait que l’on assigne la valeur d’une variable. On pourra utiliser l’opérateur Python += pour faire cette opération.

graph LR;
pred(Prédiction) --> calc(Calcul des forces)
calc --> corr(Correction)
corr -.-> pred

Figure 6.2: Schématisation de l’algorithme de Verlet. La boucle est répétée jusqu’à ce que le temps simulé \(t = n\Delta t\) soit égal au temps désiré.

6.3 Test

Pour une particule dans un champ gravitationnel \(\vec g = -g \vec{e}_y\), ayant un vecteur vitesse initial \(\vec v = v_x\vec e_x + v_y \vec e_y\), calculez en fonction du temps \(t\) la trajectoire analytique.

def analytical_position(initial_velocity: npt.NDArray[float],
                        gravity: npt.NDArray[float],
                        t: float):
    """Calcule la position d'une particule dans un champ gravitationel"""
    final_position = ...
    return final_position

Initialisez les vitesses des particules à des valeurs aléatoires :

# Les lignes suivantes servent uniquement à générer des vitesses initiales aléatoires
# Il n'est pas utile de comprendre le détail de ces fonctions
rng = np.random.default_rng()  # crée un générateur de nombres aleatoires
initial_velocity = rng.uniform(-1, 1, velocities.shape)  # rempli initial_velocity de nombres distribues entre -1 et 1

velocities[:] = initial_velocity

Consigne

À l’aide de Velocity-Verlet et du pas de temps \(\Delta t\) de votre choix, intégrez le mouvement des \(N\) particules jusqu’à \(t = 10\) s. Pour chaque itération, comparez les positions numériques avec la théorie analytique.

# Étape 1 - donner la bonne valeur (en fonction de g) aux forces
forces[:] = ...

# Étape 2 - la boucle d'intégration
for i in range(...):
    ... # on exécute Velocity-Verlet (prédiction et correction)

positions_analytical = ... # appelez la bonne fonction

# Pour tester la position des particules, on peut se servir de np.testing.assert_allclose
tolerance = 1e-14
np.testing.assert_allclose(positions, positions_analytical, atol=tolerance)

Enfin, affichez à l’aide de plt.plot() les trajectoires d’une particule en 2D.

6.4 Input-output

Matplotlib permet facilement de visualiser des trajectoires en 2D pour un faible nombre de particules, mais est assez limité pour un grand nombre de particules (à moins d’animer la trajectoire au lieu de dessiner sa trace), et surtout ses capacités de visualisation en 3D sont très limitées.

Des logiciels spécialisés, tels que Paraview permettent de visualiser des jeux de données plus gros et plus complexe que matplotlib. Dans ce cours, on utilisera un logiciel spécialisé dans le post-traitement de simulations discrètes appelé Ovito, qui a été conçu en particulier pour la dynamique moléculaire.

Astuce

La page de téléchargement d’Ovito est disponible ici. Sur les ordinateurs des salles de TP il n’est pas possible de lancer Ovito depuis le navigateur de fichiers. Il faut donc télécharger le script suivant run_ovito.ps1 et l’exécuter (clic droit > “Exécuter avec Powershell”) après avoir extrait l’archive zip de Ovito dans le dossier “Téléchargements”.

Pour pouvoir traiter les données de notre simulation, il faut les écrire dans des fichiers que Ovito pourra lire. Un format commun pour les simulations en dynamique moléculaire est le format XYZ : c’est un fichier texte contenant le nombre de particules sur la première ligne, une ligne de commentaire (ignorée par Ovito) et un tableau de 4 colonnes : la première est le symbole chimique de l’élément, les trois autres les coordonnées x y z.

La fonction ci-dessous permet d’écrire un fichier XYZ, et suppose que toutes les particules sont des atomes d’hydrogène.

def write_xyz(fname: str, positions: npt.NDArray[float]):
    """Writes an XYZ file from atom positions (2d or 3d)"""
    natoms = positions.shape[1]

    dtype = np.dtype([('symbols', np.unicode_, 2),
                      ('positions', float, 3)])
    
    extended_pos = np.zeros(natoms, dtype=dtype)
    extended_pos['positions'][:, :positions.shape[0]] = positions.T
    extended_pos['symbols'][:] = 'H'
    
    with open(fname, 'w') as fh:
        fh.write(str(natoms) + '\n\n')
        for row in extended_pos:
            s, x = row['symbols'], row['positions']
            fh.write(f"{s} {x[0]:.18e} {x[1]:.18e} {x[2]:.18e}\n")

Consigne

Servez-vous de la fonction ci-dessus pour écrire une séquence de fichiers traj_0000.xyz, traj_0001.xzy, etc. Pour rappel, on peut utiliser les f-strings Python pour formater le numéro de l’itération :

i = 21
print(f"traj_{i:04d}.xyz") # affiche le numéro avec 4 chiffres

traj_0021.xyz

Ensuite, utilisez Ovito pour ouvrir le premier fichier (traj_0000.xyz), et visualisez le mouvement des particules.

Résumé des consignes

Implémenter l’algorithme Velocity-Verlet pour \(N\) particules (en 2D ou en 3D)
Calculer l’expression analytique d’une particule avec vitesse initiale dans un champ gravitationnel uniforme
Comparer l’implémentation à la solution analytique
Afficher les trajectoires des particules
1. En 2D : avec matplotlib
2. En 2D et 3D : écrire des fichiers XYZ et visualiser avec Ovito
(Bonus) générer une video avec Ovito