Thémadoc : Le mouvement brownien tous azimuts

Comment a-t-on réalisé une simulation pour notre exemple ?

Avertissement
Les outils d'intégration HTML dont nous disposons pour effectuer ces dossiers ne nous permettent pas d'éditer tous les signes mathématiques conventionnels. Nous nous en excusons. Vous trouverez dans la version PDF (791 ko) le document original permettant de visualiser correctement les symboles.

Depuis le développement des ordinateurs, dès les années 1950, les physicochimistes ont tenté de simuler le comportement de la matière à l’état liquide au moyen de calculs où l’on cherche à déduire les propriétés macroscopiques de la dynamique des constituants microscopiques.
Pour réaliser ces simulations, il est nécessaire de bien connaître la nature des interactions entre particules.
Dans un premier temps, ces simulations peuvent donc servir à tester des modèles d’interactions entre particules : on confronte les résultats obtenus par simulations à ceux obtenus expérimentalement. Une fois les interactions entre particules convenablement décrites, les simulations peuvent avoir un caractère prédictif.
Voyons ici les problèmes posés par la mise en œuvre d’une simulation portant sur une solution aqueuse ionique. Nous avons choisi le chlorure de cryptate de potassium (K(C₁₈H₃₆O₆N₂)Cℓ ) dans l’eau lourde à la concentration :
C = 0,15 mol.L^-¹.

Simulation de dynamique moléculaire

Plusieurs méthodes de simulation ont été mises au point. La première à laquelle on peut penser consiste à résoudre les équations de la mécanique classique pour tous les atomes de l’échantillon de matière considéré. On intègre donc numériquement l’ensemble des relations f_i = m_i∙a_i (f_i est le vecteur somme des forces subies par l’atome n° i, m_i la masse de l’atome et a_i son vecteur-accélération). Connaissant les positions et les vitesses à un instant t, les accélérations, calculées à partir des forces subies permettent de calculer les vitesses à un instant t + δt. Ces vitesses conduisent ensuite de la même manière aux nouvelles positions. Cette méthode nécessite l’emploi d’un pas de temps δt extrêmement petit ; on doit en effet pouvoir négliger les éventuelles variations de f_i pendant cet intervalle de temps. Les intervalles utilisés dans cette méthode de simulation – appelée « dynamique moléculaire » – sont de l’ordre de 10^-¹⁵ s. Il est des cas où cette méthode nécessiterait des temps de calculs gigantesques, c’est le cas par exemple lorsqu’on souhaite simuler le comportement d’une solution aqueuse. Dans notre exemple, pour avoir une statistique sur 100 cations et 100 anions, il faut calculer les trajectoires de plus de 100 000 atomes !
Pour notre solution aqueuse de chlorure de cryptate de potassium, on aura donc recours à une méthode différente, appelée « simulation de dynamique brownienne ».

Simulation de dynamique brownienne

Le solvant est traité comme un fond continu caractérisé par une permittivité et une viscosité. Les ions polyatomiques y sont considérés comme des sphères chargées.

En dynamique brownienne, le solvant est considéré comme un fond continu.

Ensuite, le calcul des trajectoires par intégration numérique des équations f_i = m_i∙a_i est remplacé par une approche statistique. Chaque ion suit une marche aléatoire tenant compte des interactions avec les autres ions. À chaque pas de temps, un critère énergétique global permet d’accepter ou de refuser les nouvelles positions envisagées, avec une probabilité p dont le mode de calcul est précisé plus loin.
À partir des positions à l’instant t, les nouvelles positions des ions sont calculées à partir d’une équation stochastique du mouvement :

r(t + δt) = r(t) + R + δr^I

r est le vecteur-position de l’ensemble des ions, si N est le nombre total d’ions considérés, il a 3N composantes.
R est un vecteur-déplacement aléatoire de 3N composantes correspondant au mouvement brownien des ions. Il est de distribution gaussienne, sa moyenne dans le temps vaut <R> = 0. Sa moyenne quadratique < R∙R^T> est proportionnelle à l’intervalle de temps δt et dépend des interactions hydrodynamiques entre les ions et le solvant. δr^I est un vecteur-déplacement calculé grâce aux interactions des ions entre eux et des ions avec le solvant.

Comment a-t-on réalisé une simulation pour notre exemple ?

1. On considère un ensemble de N₊ cations et N_– anions répartis dans un volume cubique d’eau que l’on considère comme un fond continu ; l’arête de ce cube a est choisie de manière à respecter la concentration C. Pour notre exemple N₊ = N_– = 108 et a = 10,6 nm.

2. On choisit un pas de temps δt. Nous avons pris δt = 0,25 ps.

3. Partant de t = 0 et allant de δt en δt :

- On calcule de nouvelles positions pour les ions grâce à l’équation stochastique du mouvement.

- Notons r_n la configuration au temps t et r_n₊₁ la configuration au temps t + δt ; on calcule la probabilité p qu’a le système de passer de l’état r_n à r_n₊₁ :

- si p = 1, la configuration rn+1 est retenue ;

- si p < 1, on n’accepte la nouvelle configuration qu’avec la probabilité p : il suffit pour cela de tirer un nombre au hasard entre 0 et 1 et de n’accepter la configuration rn+1 que si ce nombre est inférieur à p.

4. On arrête le calcul des trajectoires au bout de 16 000 pas de temps environ.

5. On répète ce calcul cinq fois.

Quelles sont les interactions prises en compte ?

Interactions directes

À longue distance, deux ions i et j interagissent par interaction coulombienne :

où e est la charge élémentaire,

_r la permittivité diélectrique relative du solvant,

₀ la permittivité diélectrique du vide et r_ij la distance entre les deux ions.
L’énergie potentielle d’un ion i est alors la somme des énergies potentielles d’interaction :

Pour s’affranchir des effets de bord, la boîte de simulation est périodisée dans les trois directions de l’espace. Cela est analogue à ce que fit Ewald en 1921 sur des cristaux ioniques pour en calculer l’énergie réticulaire.
La figure ci-dessous illustre cette périodisation en deux dimensions seulement.

Périodisation en deux dimensions dans un cristal ionique

Pour tenir compte des interactions coulombiennes, on doit périodiser la boîte de simulation dans les trois directions de l'espace.

À courte distance, les cortèges électroniques des deux ions produisent une forte répulsion, décrite par un potentiel en 1/r⁹:

Elle fait intervenir les rayons des ions a_i et a_j ainsi qu’un paramètre B_ij ajustable. Notons qu’avec deux sortes d’ions il y a au total trois paramètres à ajuster (++, −− et +−).
On a représenté ci-dessous la somme de l’énergie potentielle d’attraction électrostatique d’un cation et d’un anion et de l’énergie potentielle de répulsion à courte portée. On voit par la montée abrupte de cette énergie potentielle que les deux ions ne peuvent s’interpénétrer.

Énergie potentielle d'interactions coulombienne et de répulsion d'un cation et d'un anion en fonction de la distance séparant leurs centres.

La solvatation de deux ions est modifiée lorsqu’ils se touchent

Comme le montre la figure ci-dessous, lorsque deux ions se touchent, ils n’engagent pas le même nombre de molécules d’eau pour leur solvatation. Lorsque deux ions se touchent, certaines molécules d'eau qui participaient à leur solvatation sont exclues.

Certaines molécules d’eau sont exclues de la couche de solvatation.

Cela induit une variation d’énergie potentielle que l’on peut exprimer en fonction de la distance entre les deux ions.

Les ions, en se déplaçant, subissent des forces hydrodynamiques de la part du solvant

La force de friction est donnée pour chaque ion i par la loi de Stokes :

ou a_i est le rayon de l’ion, η la viscosité du solvant et U_i le vecteur-vitesse de l’ion.
Il faut tenir compte également de l’effet d’entraînement des ions présents dans le sillage d’un ion en mouvement dans le solvant. Cela est possible à condition d’exprimer l’ensemble des forces hydrodynamiques subies par les ions (un vecteur F de 3N composantes) en fonction des vitesses de ces ions (un vecteur U de 3N composantes) :

F = − M∙U

où M est une matrice 3N × 3N

Comment est calculée la probabilité p de passer d’une configuration à une autre ?

La probabilité p pour passer d’une configuration r_n à une configuration r_n₊₁ est calculée selon l’expression :

où E_i est l’énergie potentielle du système dans la configuration i, k_B est la constante de Boltzmann et T la température en kelvin. P est un terme dépendant des interactions hydrodynamiques dont l’expression et la justification dépassent le cadre de cet article. On voit cependant que la probabilité de passer d’une configuration à une autre devient quasi nulle si l’énergie potentielle du système augmente considérablement lors de ce passage.

Comment a été calculée la conductivité électrique à partir des simulations de trajectoires ?

Le calcul de la conductivité électrique à partir des trajectoires simulées est un problème de physique statistique dont la solution est donnée par :

S_M est un terme dépendant des interactions hydrodynamiques, V est le volume de la solution, et C(t) est une fonction dite « d’autocorrélation des courants » :

q_i est la charge de l’ion i, U_i(t) est le vecteur-vitesse de l’ion i à l’instant t.
La figure ci-dessous représente l’intégrale de cette fonction en fonction du temps calculée à partir de 5 simulations. On voit que cette intégrale ne fait que fluctuer autour d’une valeur constante à partir de 200 ps environ. On prendra la valeur correspondant à ce plateau comme valeur de l’intégrale de 0 à l’infini.

Intégrales de la fonction d'autocorrélation des courants de cinq simulations et leur moyenne (en rouge)

C’est au prix de toutes ces précisions et de calculs durant en général plusieurs heures que le chimiste, assimilant les ions en solution à des particules browniennes en interaction, peut simuler le comportement de ces ions à l’échelle mésoscopique.
Le résultat de ces simulations, exploité à l’aide de la physique statistique, a fourni des grandeurs comme la pression osmotique, les coefficients de diffusion ou la conductivité électrique en accord avec les valeurs expérimentales.

Thierry OLYNYK