9  4/ Transition de phase

Fusion d’une nanoparticule de cuivre

Objectifs du TP

  • Calculer la température d’un système de particules
  • Implémenter et tester une méthode de contrôle de la température (thermostat)
  • Identifier un équilibre thermodynamique
  • Différentier une phase liquide et une phase solide par 1) analyse des trajectoires et 2) déplacement moyen des particules en fonction du temps

9.1 Thermostat

Le fichier cu_particle.txt contient des positions initiales pour une nanoparticule de cuivre.

Consigne

Avec les paramètres \(\varepsilon = 1\), \(\sigma = 2.334\) (Zhang et al. 2022), \(r_\mathrm c = 2.5\sigma\) et \(m = 1\) et vitesse initiale nulle, choisissez \(\Delta t\) pour que la nanoparticule soit stable. Vérifiez que l’énergie totale est conservée.

Implémentez la fonction suivante, qui calcule la température à partir des masses et vitesses instantanées :

def temperature(masses, velocities):
    """Calcule la température instantanée"""
    return ...

On rappelle la formule de la température :

\[ \frac{3}{2} N k_\mathrm{B} T = E_c\]

Puisque l’on travaille dans les unités caractéristiques du potentiel de Lennard-Jones, on prendra la constante de Boltzmann \(k_\mathrm B = 1\).

Dissipation

La configuration de la nanoparticule est une configuration d’équilibre : presque toute l’énergie du système se trouve dans l’énergie potentielle, qui est proche de son minimum théorique. Pour perturber l’équilibre, on va rajouter de l’énergie au système sous forme d’énergie cinétique. Pour cela, on définit des vitesses initiales aléatoires.

Consigne

Pour générer des nombres aléatoires avec NumPy, il faut créer un objet particulier qui va remplir ce rôle :

import numpy as np
import numpy.typing as npt

generator = np.random.default_rng()

On peut alors utiliser l’objet pour générer des séries de données de différentes distributions. Par exemple, pour générer 10 valeurs normalement distribuées avec moyenne de 0.5 et écart-type de 2 :

generator.normal(loc=0.5, scale=2, size=10)
array([ 1.43212404,  2.99489868,  2.19795038,  1.65463977,  0.43499853,
       -2.51798848,  0.3889695 , -0.13012715,  1.05426975,  0.81099987])

L’argument size peut aussi être la shape d’un autre tableau NumPy. En utilisant ces informations, générez des vitesses initiales normalement distribuées, avec une moyenne nulle est un écart-type de 1. Quelle est la température instantanée pour cette distribution de vitesse ?

Pour faciliter la visualisation dans Ovito, il est possible de donner un argument symbol à la fonction write_xyz ci-dessous :

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

    dtype = [('symbols', np.unicode_, 2),
             ('positions', float, 3)]
    dtype = np.dtype(dtype)
    
    extended_pos = np.zeros(natoms, dtype=dtype)
    extended_pos['positions'][:, :positions.shape[0]] = positions.T
    extended_pos['symbols'][:] = 'H' if symbol is None else symbol
    
    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")

Par exemple, dans une boucle de Verlet :

for i in range(1000):
    # ...
    write_xyz(f"copper_{i}.xyz", positions, symbol="Cu")

Visualisez les trajectoires avec vitesse initiale. La nanoparticule est-elle toujours à l’état solide ?

Maintenant que les particules ont une vitesse initiale, on va implémenter une partie du thermostat de Langevin : la force de dissipation \(\vec f^\mathrm d_i = -m\vec v_i / \tau\).

Consigne

Ajoutez la force de dissipation à la force de chaque particule calculée par Lennard-Jones. Enregistrez dans une liste l’énergie cinétique à chaque pas de temps, puis en fin de simulation affichez l’énergie cinétique en fonction du temps (avec plt.plot). Visualisez l’effet de la dissipation sur les trajectoires dans Ovito.

Faites varier \(\tau\). Quelle est son influence sur l’évolution de l’énergie cinétique et de l’énergie potentielle ?

Échange d’énergie

La deuxième partie du thermostat de Langevin consiste à ajouter, pour chaque particule, une force aléatoire dont l’amplitude dépend de la température désirée \(T_0\), de \(\tau\) et de \(\Delta t\) :

\[ \vec f^\mathrm{L}_i = \sqrt{\frac{2 m k_\mathrm B T_0}{\tau \Delta t}} \vec N_i(t),\]

où \(\vec N_i(t)\) est un vecteur dont les composantes suivent une loi normale de moyenne nulle et d’écart-type 1, qui doit être re-généré à chaque boucle temporelle. La constante de Boltzmann \(k_\mathrm B\) est égale à 1 dans les unités de Lennard-Jones. Cette force simule un choc avec un “bain de chaleur” (heat bath en anglais) dont les particules (non-simulées) échangeraient de l’énergie aléatoirement avec les particules du système (simulées).

Consigne

Implémentez la partie aléatoire du thermostat de Langevin. Afin de tester l’implémentation, on posera \(T_0 = 0.4\) avec des vitesses initiales nulles. De même que pour l’énergie cinétique précédemment, calculez la température à chaque pas de temps, puis affichez l’évolution de la température en fin de simulation. Quel est l’effet de \(\tau\) sur l’évolution de la température ? Quand est-ce que le système atteint un équilibre avec son environnement ?

9.2 De vibration à diffusion

On cherche maintenant à déterminer, en fonction de la température, l’éloignement des particules par rapport au centre de gravité. On s’attend à ce que pour les états solide et liquide les atomes restent proches du centre de gravité, mais s’en éloignent pour un gaz.

Pour calculer la position des particules à l’équilibre thermodynamique, on simule l’évolution du système pour un temps \(t_\mathrm{tot} = 75\tau_{LJ}\) avec le thermostat de Langevin. Le temps caractéristique de dissipation est mis à \(\tau = 5 \tau_{LJ}\).

Pour calculer la position d’une particule par rapport au centre de gravité, on peut se servir de la fonction mean de numpy :

r_center = r - r.mean(axis=1)[:, np.newaxis] # valide si toutes les masses sont identiques

Une fois la distance calculée (avec np.linalg.norm, comme vu dans un précédent TP), on peut utiliser np.histogram pour calculer la distribution des distances par rapport au centre. Cette fonction retourne le nombre d’atomes par intervalle de l’histogramme, ainsi que les bornes de chaque intervalle.

Consigne

Implémentez une fonction qui accepte une température en paramètre, fait la simulation d’équilibre thermodynamique décrite ci-dessus et retourne le résultat de la fonction np.histogram.

Pour 4 valeurs de température de \(T = 0.1\) à \(T = 1.0\), affichez l’histogramme dans un graphe avec une légende. Estimez quelle(s) courbe(s) correspondent aux phases liquide, solide et gazeuse.

Pour la phase gazeuse, faites varier \(t_\mathrm{tot}\) (à température constante) et observez l’évolution de l’histogramme.

Zhang, Yuwei, Yunchu Wang, Fei Xia, Zexing Cao, et Xin Xu. 2022. « Accurate and Efficient Estimation of Lennard–Jones Interactions for Coarse-Grained Particles via a Potential Matching Method ». Journal of Chemical Theory and Computation 18 (8): 4879‑90. https://doi.org/10.1021/acs.jctc.2c00513.