Dynamique Moléculaire Classique Bernard Rousseau Laboratoire de Chimie-Physique, Université Paris-Sud, Orsay

1. Dynamique Moléculaire Classique Bernard Rousseau Laboratoire de Chimie-Physique, Université Paris-Sud, Orsay « Modélisation, dynamique et thermodynamique des agrégats et molécules complexes », Porquerolles, 22-27 mai 2005

2. Quelques ouvrages de référence Computer Simulation of Liquids, M.P. Allen D.J. Tildesley Understanding Molecular Simulation, Daan Frenkel ,B. Smit Molecular Dynamics Simulation, J. M. Haile Molecular Modelling: Principles and Applications, Andrew Leach The Art of Molecular Dynamics Simulation, D. C. Rapaport Theory of Simple Liquids, J-P Hansen, I.R. McDonald Basic Concepts for Simple and Complex Liquids, J-L Barrat, J-P Hansen

3. Simulation de la dynamique moléculaire classique Définition : méthode permettant de calculer les propriétés d’équilibre (P, T, ...) et les propriétés de transport (D, h, ...) d’un système constitué de N particules. • Réalisation moderne (et plus modeste) d’un vieux rêve des sciences, issu du déterminisme classique : « Une intelligence qui connaîtrait à un instant donné l’ensemble des positions et des vitesses de toutes les particules de l’Univers et les forces qui régissent leurs mouvements, pourrait en déterminer le futur et retrouver le passé... » Pierre Simon de Laplace Hypothèse : le mouvement des particules considérées obéit aux lois de la mécanique classique (valide pour la plupart des systèmes, sauf He, H2, D2 si hn > kBT).

4. Dynamique Moléculaire : principes • Modèle pour les particules • atome, qq atomes, molécule • nature des interactions • Modèle pour l’environnement • isolé/périodique, ouvert/fermé, homogène/interface • Energie, température, pression cste • Hors-équilibre • Trajectoire pour chaque particule

5. Dynamique Moléculaire : principes • « Mesurer » des quantités observables • lien entre observable et {rN, vN} : mécanique statistique • Equipartition de l’énergie : • Température « instantannée » • Moyenne temporelle (hyp. ergodicité) • Propriétés dynamiques : corrélations

6. Modélisation moléculaire - Motivations SYSTÈME RÉEL SYSTÈME MODÈLE construction de modèles Expérimentation Simulation Théorie résultats exacts pour le modèle prédictions théoriques résultats expérimentaux comparaison comparaison test des modèles test des théories

7. Equations du mouvement

8. Trajectoire : résoudre les équations du mouvt de Newton

9. Formalisme de Lagrange Lagrangien Equations du mouvement • Système atomique, coordonnées cartésiennes

10. Formalisme de Hamilton • Hamiltonien • Equations du mouvement • Système atomique, coordonnées cartésiennes

11. ri(t+∆t) ri(t+2∆t) ri(t) Intégration des équations du mouvement intégration analytique impossible pour n ≥ 3 (problème à 3 corps ...) Intégration par pas, méthode des « différences finies » r(t+∆t) r(t) v(t) a(t) v(t+∆t) a(t+∆t) v∆t << dimensions particule

12. Algorithme de Verlet Développement en série de Taylor de r(t+∆t) au voisinage de t Positions erreur en t4 Vitesses erreur en t2

13. Algorithme de Verlet - Version « Verlet-vitesse » • Mise en œuvre Etape 1 -> calcul des forces à t+∆t Etape 2

14. Algorithme de Gear (prédicteur-correcteur) Prédiction Calcul des forces -> a(t+∆t) Correction En MD, une seule itération suffit ...

15. Sources d’erreur  2 types d ’erreur - erreur de troncature : approximation « différence finies » - erreur d ’arrondi : implémentation de l’algorithme sur une machine particulière erreur globale arrondi troncature pour un temps de simulation donné pas d’intégration, ∆t

16. Erreurs d’arrondi  Erreurs d ’arrondi Options de compilation Représentation des nombres réels => REAL*8 !!!

17. Erreurs de troncature  Séparation de 2 trajectoires - perturbations initiales : 10-3, 10-6 et 10-9 s Sensibilité aux Conditions Initiales

18. Quel algorithme choisir ? Simplicité et consommation  Verlet : - r(t), r(t-∆t), a(t) => 3 « mots » / particule / deg. de liberté - 2 additions, 1 multiplication  Gear (éq. diff. 1er ordre) : - r(t), v(t), a(t), b(t), F(t) => 5 « mots » / particule /degré de liberté - 10 additions, 10 multiplications Verlet, 106 particules : 3*8 * 3 * 1000000 ≈ 72 Mo calcul des forces ≈ 90% du temps de calcul

19. Quel algorithme choisir ? Conservation de l’énergie • Gear : meilleure conservation à temps court • Verlet : faible dérive d’énergie à temps long • Gear : toutes les coordonnées au même instant, mais pas time-reversible...

20. De l’atome à la molécule Modèles moléculaires

21. De l’atome à la molécule • Particules interagissant via des forces - inter-moléculaires - intra-moléculaires Représenter les degrés de liberté interne : flexibilité moléculaire

22. Potentiels intra-moléculaires élongation torsion pliage interactions « non-liées » • Interactions non-liées : dispersion-répulsion, au-delà de 4 atomes

23. Potentiels intra-moléculaires • Flexibilité  mouvements rapides (élongation) => ∆t petit  mouvements de faible amplitude (par rapport dimensions moléculaires) => découplés de la dynamique Molécule rigide Contraintes géométriques (gel degrés liberté)

24. 2 3 1 Méthode des contraintes • Objectif : maintenir les distances entre particules constantes au cours de la dynamique • Méthode : contraintes holonomiques - - bille sur sphère : - si gravité : non-holonomiques !

25. Méthode des contraintes • Modifier les équations du mouvement par addition d’une force « de contrainte : • Contraintes et forces de contrainte : - -

26. Méthode des contraintes - Molécule diatomique • Une contrainte, , un seul multiplicateur l1 • Equations du mouvement • Equation pour la contrainte Note : soit • Solution pour l (Moriss & Evans)

27. Méthode des contraintes • Méthode directe : - cas général => inversion matrice (nc x nc) - problème de « dérive » des longueurs de liaison (la résolution des équations du mouvement est approximative) => fonction de correction • Méthode : SHAKE/RATTLE : Ryckaert et al. (1977) - résolution iterative où les contraintes sont satisfaites explicitement à chaque pas de calcul.

28. Méthode des contraintes - SHAKE • Molécule tri-atomique, algorithme de Verlet (1) (2) • Eq. (1) et (2) au carré, négligeant les termes en l2, et en utilisant les relations : on obtient 2 équations linéaires pour 2 contraintes...

29. Méthode des contraintes - SHAKE • SHAKE : les contraintes sont traitées successivement et itérativement - contrainte 1-2 : Soit : • Convergence lorsque toutes les contraintes sont inférieures à 10-6 / 10-8

30. DM de molécules « rigides »  Molécules rigides - gel des degrés de liberté interne si les tvib << ∆t => gain de temps (benzène, xylènes, CO2, ...)  Séparation de la dynamique : translation CM + rotation Rotation Translation Orientation : angles d ’Euler, quaternions

31. DM de molécules « rigides » - Angles d’Euler singularité pour q = 0, p

32. DM de molécules « rigides » - Quaternions Eq. diff. 1er ordre

33. Modéliser l’environnement

34. Conditions de simulation • Milieu homogène • Agrégat dans le vide • - 10000 atomes • - 50-150 Å • - 1 à 100 ns • corps purs ou • mélanges • Interface explicite (LV, SL) Différents ensembles statistiques {NVT}, {NPT}, {µVT}, ...

35. Ensembles statistiques Dynamique Newtonienne : ensemble {NVE} Autres ensembles : {NVT}, {NPT} - méthodes « stochastiques » : la distribution des vitesses est périodiquement échantillonnée dans une distribution de Maxwell-Boltzmann. Les équations du mouvement sont inchangées. - méthodes « homothétiques » (Berendsen) : les vitesses ou les positions sont modifiées pour introduire un changement de température ou de volume. Les équations du mouvement sont inchangées. - méthodes « étendues » (Nosé, Hoover, Andersen) : le système est en « équilibre » avec un thermostat ou un barostat fictif. Le nombre de degrés de liberté du système change. Les équations du mouvement sont modifiées. - méthode de Gauss (Evans, Morris): on utilise le principe de Gauss pour contraindre une variable thermodynamique

36. Dynamique à température constante • Modifier ou maintenir la température du système • « Recalage » des vitesses

37. Dynamique à température constante • Couplage avec un bain thermique - Berendsen Tbain • décroissance exp. controlée par • t ( ≈ 1 ps) t • variation de température • à chaque pas

38. Dynamique à température constante • Bain thermique de Berendsen • - fluctuations de température ( ≠ recalage) • maintien des différences de température entre les différents • constituants du système (solvent chaud/soluté froid) • - ne génère pas rigoureusement {NVT} • Méthode stochastique (Andersen) • une particule est choisie aléatoirement à intervalles réguliers • sa vitesse est ré-attribuée dans une distribution Gaussienne • chaque changement de vitesse correspond à « une collision » avec un bain fictif... • trajectoire discontinue...

39. Dynamique à température constante - Méthode « étendue » • Système étendu (Nosé-Hoover) • le bain thermique est représenté par un degré de liberté supplémentaire : • énergie potentielle : • énergie cinétique :

40. Dynamique à température constante - Méthode « étendue » • Influence de la « masse fictive » Q (argon supercritique)

41. algorithme d ’intégration espace des phases G Mise en œuvre Configuration initiale - définir {ri(t), vi(t)} + ai(t), ... Positions initiales - aléatoire -> recouvrement - réseau cristallin :  molécules anisotropes - gaz faible densité + compression

42. Mise en œuvre Vitesses initiales - distribution uniforme [-va,va] - distribution Maxwell-Boltzmann - annuler le moment total

43. Stabilisation  Butane : compression isotherme (510 K) de

44. Stabilisation  Butane : compression isotherme (510 K) de

45. Mise en œuvre  Conservation de l’énergie ...

46. Mise en œuvre  Visualisation

47. Optimisation et réduction du temps de calcul  Temps de calcul (estimation) 1 force ≈ 100 opérations élémentaires - ncf = 1000 => FB = 0,5 106par cycle - RISC ≈ 50 MFlops => 1 ns / heure  Evolution des ordinateurs Diminution du temps de calcul  Algorithmique

48. Algorithmique - Listes de voisins (Verlet)  Calcul « brutal »  Avec liste de Verlet

49. Algorithmique - Méthode des cellules  Calcul avec cellules  Calcul avec cellules + liste

50. Algorithmique - Calcul des forces Nombre de forces FB FRC N