Integral Equations and Monte Carlo molecular simulations
logo_tutelle logo_tutelle 

Equations intégrales et simulations Monte Carlo moléculaires

Toute la nouveauté tient dans le qualificatif "moléculaire" : Il s'agit de calculer les propriétés d'équilibre, thermodynamique et de structure pour des fluides de particules non sphériques. En clair, comment calculer à partir d'un potentiel de paires v(12) entre 2 particules la fonction de distribution de paires (pdf) g(12), l'équation d'état P(r), le facteur de structure S(q), l'intensité diffusée I(q) en diffusion de rayons X ou de neutrons… pour une température T et une densité ρ données ?

La difficulté par rapport au cas "atomique" très classique et bien maîtrisé depuis 30 ans des particules sphériques est que maintenant (12) représente non seulement la distance r entre les 2 particules 1 et 2 mais aussi l'orientation relative Ω de chacune d'elles et du vecteur les joignant, g(12)=g(r,Ω). Plusieurs cas de telles particules anisotropes sont rencontrés dans les recherches effectuées au LIONS. Il peut s'agir de particules linéaires comme les sphères dipolaires ou encore les micelles de tensioactif cylindriques ou ellipsoïdales, les nanotubes d'imogolite… dont l'orientation relative est repérée par 3 angles d'Euler. De manière plus complexe, il peut aussi s'agir de particules non linéaires comme les modèles type SPCE (interaction site-site) de molécules d'eau (H2O vu comme 3 sites à charges partielles) ou comme les "molécules colloïdales" de tétrapodes constituées d'un cœur de silice et de 4 sites de latex en positions tétraédriques. Cette fois-ci, 5 angles d'Euler sont nécessaires !

En mécanique statistique des liquides, deux approches théoriques très complémentaires sont en compétition depuis les années 1960. La 1ère utilise la simulation numérique Monte Carlo (MC) ou la Dynamique Moléculaire, la 2nde résout l'équation Ornstein-Zernike (OZ) avec une équation intégrale, généralement du type HNC (hypernettted chain). La 1ère est exacte, en ce sens que la pdf obtenue au final correspond exactement, en théorie..., au potentiel de paires v donné en entrée. La 2nde apporte une solution approchée, car l'équation HNC n'est qu'une approximation. La 1ère est "facile" à programmer et ne demande aucun artifice (des codes très versatiles sont même proposés sur le web gratuitement), la 2nde peut demander des codes lourds faisant appel à des techniques d'analyse numérique pointues. Mais, les calculs pour obtenir la solution exacte (1ère méthode) sont très longs pour des particules anisotropes, souvent au delà de la limite des capacités des ordinateurs de bureau puissants actuels. Il faut une heure, une nuit, un week-end, une semaine… pour sortir des pdf (pair distribution funciton) ou des intensités diffusées I(q) dont le bruit statistique reste raisonnable.

Au contraire, la 2nde, grâce aux codes très robustes et puissants, est résolue sur ordinateur en quelques secondes ou minutes et les courbes finales sont vierges de bruit statistique. En conclusion, la simulation n'est utilisée que pour quelques cas typiques qui servent de référence et l'équation intégrale est résolue de manière plus systématique pour, par exemple, le dépouillement de spectres expérimentaux. La comparaison des deux techniques permet d'estimer le degré d'approximation de HNC et, si nécessaire, de montrer la direction à suivre pour améliorer la résolution de l'équation intégrale. L'équipe du LIONS a une longue expérience de ces deux techniques pour des particules sphériques (30 ans pour HNC, 15 ans pour MC) et les a développées dans les 5 dernières années pour des particules anisotropes à l'occasion du projet ANR SISCI.

Simulations Monte Carlo (MC)

Cette technique est bien-sûr très classique. Il s'agit de placer N (100-1000) particules dans une boîte cubique de simulation dont la taille est imposée par la densité et de générer une suite (lire plusieurs millions) de configurations en accord avec le potentiel d'interaction et la température. Il "suffit" alors de construire un histogramme des distances (et orientations!) observées pour obtenir la pdf. Seule remarque centrale: la simulation numérique particulaire est un problème essentiellement de bruit statistique qui décroît comme la racine carrée du nombre de configurations observées (indépendantes). Donc, si au bout d'une nuit de calcul on obtient des courbes peu jolies dont on veut diviser le bruit par 10, il faut continuer la simulation pendant … 2 mois!

Le LIONS utilise actuellement MC pour des fluides dipolaires qui sont des modèles simplistes de l'eau (voir page "effets spécifiques ioniques") et sur les molécules colloïdales de tétrapodes. Ce denier cas est illustré ici: les tétrapodes sont fortement imbriqués au sein de la cellule et le spectre de diffusion, donné pour 2 nombres N de particules, est compatible avec une structure diamant de l'organisation.

 

 

Equation intégrale HNC (hypernetted chain equation)

Tout démarre par l'équation exacte OZ qui relie fonctions de corrélation totale h(12)=g(12)-1 et directe c(12) via un produit de convolution (spatial et orientationnel):

γ= h-c = ρ h*c                      (OZ)

Une 2nde relation dite de fermeture, toujours exacte,  s'écrit comme:

g = exp ( -v/kT + h - c + b )

Tout le problème est que la dernière fonction, dite bridge, b(12), ne peut pas s'exprimer simplement en fonction des précédentes. Le graal du physicien de mécanique statistique des liquides est de connaître la meilleure approximation possible pour cette fonction. Dans l'immense majorité des cas, on la néglige purement et simplement, ce qui définit l'équation intégrale HNC:

 g = exp ( -v/kTγ  )           (HNC)

Voilà, tout est là, il faut résoudre l'ensemble des 2 équations OZ, HNC à 2 inconnues g et c (ou g et γ). Pour des particules anisotropes, il est illusoire de résoudre en fonction de r et des 3 ou 5 angles d'Euler. L'idée est plutôt pour chaque r de développer la variation orientationnelle sur une base de fonctions angulaires aux propriétés bien connues et pratiques, les invariants orientationnels ou harmoniques sphériques généralisées définis par 5 indices mnlμν, et de travailler avec les coefficients qui ne dépendent plus que de r:

 g(r,Ω) = Σ gmnlμν(r) Φmnlμν(Ω)

Pour fixer les idées, Φ00000=1 donc g00000 donne la distribution centre de masse-centre de masse; Φ11000=u1u2  donc g11000 dira, par exemple, comment les dipôles sont alignés et donnera la constante diélectrique du solvant; Φ11200=3(u1r)(u2r)-u1u2 qui donne la variation angulaire du potentiel dipôle-dipôle…

En théorie, la base est infinie. En pratique, un nombre fini de termes suffit pour résoudre correctement la variation angulaire de la pdf. Pour des fluides dipolaires, on peut s'en sortir avec 10-20 projections. Pour le modèle SPCE de l'eau, il faut plutôt compter 250-500 projections! La résolution numérique fait appel à des transformées de Hankel (le produit de convolution dans OZ devient un produit algébrique simple dans l'espace de Fourier), à des quadratures de Gauss, à des techniques de convergence Newton-Raphson…

 

Maj : 13/10/2017 (1863)

 

Retour en haut