M.F. Ruiz-L ó pez, UMR CNRS-UHP 7565 Université Henri Poincaré, Nancy

1. GDR « Agrégation, Fragmentation et Thermodynamique des Systèmes Complexes Isolés » Modélisation, Dynamique et Thermodynamique Des Systèmes Complexes Isolés Ecole thématique du CNRS Porquerolles, 22 au 27 mai 2005 Modélisation et dynamique de systèmes complexes avec un champ de force mixte QM/MM ou QM/QM' M.F. Ruiz-López, UMR CNRS-UHP 7565 Université Henri Poincaré, Nancy

2. I. Historique, aspects fondamentaux, mise en œuvre • Problèmes posés par l’étude de gros systèmes moléculaires • La méthode QM/MM • Traitement des liaisons frontière, la méthode LSCF • ONIOM, approches QM/QM • Couplage avec la Dynamique Moléculaire • Couplage avec le modèle du Continuum • Exemple : étude d’une molécule d’eau dans l’eau • II. Potentiels d’interaction. • Nécessité d ’une paramétrisation • Stratégies de paramétrisation • Au delà du potentiel de Lennard-Jones • (Termes électrostatiques dans les méthodes semiempiriques) • III. Simulation de réactions chimiques • Difficultés • Calcul de profils d ’énergie libre • Calcul de trajectoires

3. Articles de revue • Gao J.,Methods and Applications of Combined Quantum Mechanical and Molecular Mechanical Potentials, Lipkowitz K.B; Boyd D.B., Eds.; Reviews in Computational Chemistry Vol. 7; VCH, 1996; 119-185 • P. Amara and M. Field Combined Quantum Mechanics and Molecular Mechanics Potentials, Encyclopedia of Computational Chemistry, Vol. 1, Ed. P.v.R. Schleyer, Wiley & sons, 1998, Pg 431-437 • M. F. Ruiz-López and J.L. Rivail, Combined Quantum Mechanics and Molecular Mechanics Approaches to Chemical and Biochemical Reactivity., Encyclopedia of Computational Chemistry, Vol. 1, Ed. P.v.R. Schleyer, Wiley & sons, 1998, Pg 437-448 • G. Monard and K.M. Merz Jr. , Combined Quantum Mechanics and Molecular Mechanics Methodologies Applied toBiomolecular Systems, Acc. Chem. Resd. 1999, 32, 904-911 • THEOCHEM Special Issue « Combined QM/MM calculations in Chemistry and Biochemistry », Vol. 632 (2003)

4. Problèmes posés par l’étude de gros systèmes moléculaires • Objectifs : • Décrire les propriétés électroniques d’un système « complexe »: formation ou rupture des liaisons chimiques, spectroscopies,…. • Difficultés : • Grand nombre d’électrons • Grand nombre de degrés de liberté nucléaires • Approches « traditionnelles » : • Calculs poussés sur un système modèle de petite taille • Calculs approchés (semiempiriques) sur un modèle « réaliste » • Approches actuelles : • Méthodes quantiques à croissance linéaire • Dynamique Moléculaire « Ab initio », Car-Parrinello, …. • Méthodes mixtes QM/MM, ONIOM, QM/QM’…

5. Problèmes posés par l’étude de gros systèmes moléculaires Comparaison du temps de calcul en fonction de la méthode (les calculs ab initio et DFT utilisent une base 6-31G*)

6. La méthode QM/MM • Principe : • Processus chimique (plus ou moins) local • Il est donc possible de diviser le système complet en deux parties : • Sous-système chimique : il est décrit par une approche de mécanique quantique: ab • initio, DFT, semiempirique • Environnement : il est décrit par un champ de force de mécanique moléculaire QM System QM System MM system MM system

7. La méthode QM/MM Un modèle quantique/classique pionnier: le modèle du continuum • le soluté est décrit par une méthode de Chimie Quantique • le solvant est représenté par un milieu continu diélectrique (polarisable) • interaction électrostatique prise en compte dans l’hamiltonien du soluté • interaction non-électrostatique évaluée par des formules paramétrées • D. Rinaldi & J. L. Rivail , Theoret. Chim. Acta, 32 (1973) 57 • O. Tapia & O. Goscinski, Mol. Phys., 29 (1975) 1653 • S. Miertus, E. Scrocco, & J. Tomasi, Chem. Phys., 117 (1981) 117 Diélectrique polarisable Molécule « quantique » H = H° + Vsoluté-solvant

8. La méthode QM/MM Premiers modèles QM/MM Niveau semiempirique Warshel & Levitt J. Mol. Biol., 103 (1976) 227 P. A. Bash, M. J. Field & M. Karplus, JACS, 109 (1987) 8092 M. J. Field, P. A. Bash & M. Karplus, J. Comp. Chem., 11 (1990) 700 Ab initio /MM Alagona, Desmeules, Ghio, Kollman (JACS 1984) Singh, Kollman (JCC,1986) Stanton, Little, Merz (JPC 1995) Freindorf,Gao (JCC 1996) DFT/MM Stanton, Hartsough, Merz (JPC, 1993, 1995), ions en solution Tuñón, Martins-Costa, Millot, Ruiz-López (JMM, CPL, 1995) molécules en solution Parrinello et col. (JCP, 1999), ondes planes Local SCF method Théry, Rinaldi, Rivail, Maigret, Ferenczy (semiempirical, JCC 1994) Assfeld, Rivail (ab initio, CPL 1996)

9. La méthode QM/MM HTOTAL= HQM+ HQM / MM+ HMM HQM Hamiltonien du système QM HQM Hamiltonien du système MM MM QM HQM / MMHamiltonien d’interaction On suppose qu’il n’y a pas de liaison chimique entre les sous-systèmes QM et MM, pour le moment

10. La méthode QM/MM Termes électrostatiques On suppose que le sous-système MM contient des charges qs et des moments dipolaires permanents °s s, t MM centres n QM noyaux i QM électrons Les champs électriques EMM et EQM sont définis par avec

11. La méthode QM/MM Termes de polarisation (cas d’un champ de force MM polarisable) Supposons maintenant que le sous-système MM est polarisable avec inds = moments dipolaires induits : s, t MM centres n QM noyaux i QM électrons (s = polarisabilité isotrope) Énergie nécessaire à la création du dipole induit (théorie de la réponse linéaire) La somme des 2 termes conduit à :

12. La méthode QM/MM Termes de van der Waals Interactions non-électrostatiques. En général = Lennard-Jones. Excellent rapport qualité/prix dans les simulations standard. with répulsion dispersion Principales limitations du potentiel de Lennard-Jones : - pas très flexible (2 paramètres seulement) - indépendant des coordonnées électroniques => correction « post-SCF »

13. La méthode QM/MM Equations Hartree-Fock (Kohn-Sham) l’opérateur de Fock doit contenir tous les termes dépendant explicitement des coordonnées électroniques. Modèle non-polarisable Modèle polarisable L ’opérateur doit inclure en plus la contribution à l ’énergie de polarisation totale (MM + QM/MM) qui contient de termes quadratiques de la matrice densité. La contribution est : La résolution des équations SCF et le calcul des moments induits dans la région MM doivent être faits par un processus itératif auto-cohérent. Coût = environ 5 fois celui d ’un modèle non-polarisable !

14. La méthode QM/MM Energie totale  fonction d ’onde perturbée du sous-système QM Dérivées premières Calcul analytique possible (comme dans la méthode QM) => L’optimisation de géométrie du système QM/MM est possible => forces sur les atomes, dynamique moléculaire

15. QM MM MM Traitement des liaisons frontière, la méthode LSCF Recent reviews Reuter, N.; Dejaegere, A.; Maigret, B.; Karplus, M. J Chem Phys, 2000, 104, 1720 N. Ferré, X. Assfeld, J.L. Rivail, J. Comp. Chem., 2002, 23, 610 G Monard, X. Prat-Resina, A. González-Lafont, J.M. Lluch, Int. J. Quantum Chem., 2003, 93, 229

16. Traitement des liaisons frontière, la méthode LSCF • Approche “Link-atom” • Singh, Kollman (JCC,1986) • M. J. Field, P. A. Bash & M. Karplus (JCC, 1990) • On ajoute des atomes H pour saturer la valence des atomes QM frontière. • Ces atomes n’interagissent pas avec les atomes MM (pas toujours) • Ils sont inclus dans le calcul QM • Avantages : • Méthode intuitive, simple à mettre en œuvre, réaliste si les atomes frontière sont loin. • Inconvénients : • L’expression de l ’énergie total et celle de ses dérivées ne sont pas bien définies • Rôle très important des charges MM proches

17. Traitement des liaisons frontière, la méthode LSCF Méthodes dérivées du « Link-atom » Pseudo-halogen (Hyperchem, 1994) Cummins, Gready (JPC B 2000) Pseudo-bond approach Zhang, Lee, Yang (JCP 1999) Connexion Atoms Antes, Thiel (JPC A, 1999) H  Y (7 electrons) Parameterized atom Effective Core potentials

18. Traitement des liaisons frontière, la méthode LSCF Méthodes utilisant des orbitales gelées LSCF approach (Local SCF) Théry, Rinaldi, Rivail, Maigret, Ferenczy (semiempirical, JCC 1994) Assfeld, Rivail (ab initio, CPL 1996) Ferré, Assfeld, Rivail (JCC 2002) Fragment SCF Náray–Szabó, Surján (CPL 1983) Náray–Szabó (Comp Chem 2000) GHO approach (Generalized Hybrid Orbital) Gao, Amara, Alhambra, Field (JPC, 1998) Frozen orbital QM/MM Philipp, Friesner (JCC 1999)

19. Traitement des liaisons frontière, la méthode LSCF Localized SCF (LSCF) Théry, Rinaldi, Rivail, Maigret, Ferenczy (semiempirical, JCC 1994) Assfeld, Rivail (ab initio, CPL 1996) Ferré, Assfeld, Rivail, J. Comput Chem., 23, 610 (2002) • Principe de la méthode : • Décrire les liaisons frontière par des orbitales hybrides strictement • localisées et gelées pendant le calcul SCF • (SLBOs (= Strictly localized bond orbitals) • Pour ce faire, on doit transformer la base d’AOs • Les SLBO sont obtenues dans un calcul de référence Orbitales hybrides SLBO Base d’orbitales atomiques

20. Traitement des liaisons frontière, la méthode LSCF • Méthode : • 1) Orthogonalisation des SLBO L orbitales gelées orthonormales : • 2) Soustraire la projection des orbitales atomiques sur les SLBOs • 3) Éliminer les L dépendances linéaires (orthogonalisation canonique) • B est rectangulaire Nx(N-L) (équivalente à S-1/2 dans le SCF standard) • On obtient (N-L) fonctions m ’ orthogonales entre elles et aux SLBOs

21. Traitement des liaisons frontière, la méthode LSCF Related method : the GHO approach (Generalized Hybrid Orbital) Gao, Amara, Alhambra, Field (JPC, 1998) - Semiempirical level - Use of 4 hybrid orbitals at the MM frontier atom (“hybrid atom”) - Auxiliary orbitals are not included in the SCF calculation - Auxiliary orbitals + nucleus charge act as an Effective Core Potential (ECP) - The ECP is optimized modifying the original semiempirical parameters of the boundary atom

22. ONIOM et méthodes QM/QM • Principe de la méthode ONIOM • Choix d’un sous-système, le « modèle » • 3 calculs: • système modèle au niveau de calcul HAUT => E1 • système modèle au niveau de calcul BAS => E2 • système complet au niveau de calcul BAS => E3 • L’énergie du système complet au niveau HAUT, E4, • est estimée à l’aide de l’équation : • E4 = E3 + (E1 - E2) • Liaison frontière : link-atom Low level « modèle » High level IMOMM F. Maseras and K. Morokuma, J. Comp. Chem. 16 (1995) 1170-1179. IMOMO S. Humbel, S. Sieber and K. Morokuma, J. Chem. Phys. 105 (1996) 1959-1967. ONIOM M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber, and K. Morokuma, J. Phys. Chem., 100, 19357-19363 (1996).

23. ONIOM et méthodes QM/QM Réaction PM3 B3LYP PhOH  PhO¯ + H+ 331.1 350.9 DMP + H+ DMPH+ -194.7 -238.5 PhOH + DMP  PhO¯ + DMPH+136.4 112.4 E(B3LYP) = E (PM3) - 24 kcal/mol Équation valable pour d ’autres systèmes proches ?

24. ONIOM et méthodes QM/QM E(B3LYP) = 25.2 + 0.989 E(PM3) R= 0.999 Correlation between PM3 and B3LYP results for the energies of (p-R-PhO¯ + HDMP+) and (p-R-PhO¯··HDMP+) (relative energies with respect to the separated neutral molecules).

25. ONIOM et méthodes QM/QM

26. ONIOM et méthodes QM/QM Méthode QM/SE-I Cui, Guo & Karplus, J. Chem. Phys., 117, 5617 (2003)) Système SE (semiempirique) • Principe: • Approximation NDDO pour toutes les interactions SE et QM/SE • Matrice densité QM/SE => Pmn = 0 • Seules les interactions de Coulomb QM/SE sont donc prises en compte au niveau de l’opérateur de Fock • Les interaction non-électrostatiques sont estimées par un potentiel de van der Waals (comme dans QM/MM) • Les interactions coulombiennes sont calculées par l’une de ces 2 méthodes : • base auxiliaire de gaussiennes pour décrire  dans SE • distribution multipolaire pour représenter  dans SE (charges de Mulliken dans la pratique) Système QM (ab initio ou DFT)

27. Core b Core a Buffer a Buffer b ONIOM et méthodes QM/QM Méthode basée sur l ’approche Divide & Conquer V. Gogonea, L.W. Westerhoff and K.M. Merz Jr., J. Chem. Phys., 113, 5604 (2000) Principe: T valeur fictive (1000 K, valeur typique) eF determiné itérativement, Nb total é constant

28. ONIOM et méthodes QM/QM Méthodes basées sur la DFT Gordon-Kim, JCP, 1972 Principe: Supposons que l’on connaît la densité électronique des systèmes A et B. L ’énergie du système AB, EAB est calculée à partir de la densité AB = A +B en utilisant la fonctionnelle de la densité : EAB = F [AB] = T [AB] + Vne [AB] + Vee [AB] T kinetic term Vne Coulomb Z - e- Vee Coulomb e- - e- + exchange-correlation Système A Système B

29. ONIOM et méthodes QM/QM Implémentations : Frozen-density: Wesolowski, Warshel JPC, 1993 Méthode itérative: Jeanvoine, Ruiz-Lopez, 1994 Extensions possibles pour un calcul QM/QM général Carter et col., JCP 1999. B A Calcul SCF de A en présence du potentiel crée par B (charges ponctuelles, développement multipolaire,…). Calcul itératif nécessaire si B n’est pas gelée. Principale difficulté : T[] Thomas-Fermi CF ∫5/3 (r) dr (LDA) von Weizsacker (correction 2ème ordre)

30. ONIOM et méthodes QM/QM Energie d’interaction (HCl)2 Y. Jeanvoine, D.E.A., Nancy 1994 Calcul DFT utilisant F[a+b] Calcul DFT standard

31. Couplage avec la dynamique moléculaire Approximations de base: Traitement « classique » des noyaux Born-Oppenheimer Algorithmes possibles: 1) Calculs QM/MM sur un ensemble de configurations issues d’une simulation de dynamique moléculaire classique 1 simulation MD classique + quelques centaines de calculs QM/MM Limitations : celles du champ de force (réactions chimiques) 2) Calcul auto-cohérent QM/MM en présence d’un champ électrostatique « moyen » obtenu par un simulation de dynamique classique. (Méthode ASEP = average solvent electrostatic potential) plusieurs simulations MD classiques + quelques calculs QM Limitations : pas de dynamique explicite pour le soluté 3) Simulations avec forces QM/MM calculées « on the fly » 1 simulation QM/MM impliquant plusieurs milliers (millions) de calculs QM/MM

32. Couplage avec la dynamique moléculaire • Simulations QM/MM « on the fly » : • Pour une configuration initiale donnée, il faut résoudre les équations SCF du système QM en interaction avec le système MM • La fonction d’onde est stockée (« guess » de l’itération suivante) • Calcul des forces sur les noyaux QM et MM • Résoudre les équations du mouvement, déduire les nouvelles positions des noyaux au temps t + Dt. • Itérer le processus pour la nouvelle configuration. • Dt : suffisamment petit pour un échantillonnage correct du système, fonction des modes de vibration les plus rapides => typiquement les liaisons XH…. (Dt = 0.2-1.0 fs )

33. Couplage avec la dynamique moléculaire • Conditions initiales • Configuration initiale : peut être déduite d’une simulation purement classique. Cependant, le système doit être re-équilibré car le champ de force change !!! • Fonction  d ’essai : à chaque pas de la dynamique, la fonction d ’onde est calculée : le choix d ’une bonne fonction d ’essai est donc très important afin de réduite au minimum le nombre d ’itérations SCF : • Option standard : utiliser  calculée dans l’itération précédente • Plus efficace : utiliser une technique d’extrapolation (matrice densité, densité électronique) Percentage of SCF calculations requiring N iterations. Comparison between two simulations of 300 fs using the standard density guess (full line) and an extrapolated guess (dashed line)

34. Continuum diélectrique polarisable, constante diélectrique o MM system QM system Couplage avec le modèle du continuum S. Chalmet, D. Rinaldi and M.F. Ruiz-López Int. J. Quantum Chem., 84, 559 (2001) Continuum (eo) Continuum (eo) Etotal = EQM + EQM/MM + EMM + Esolv

35. Exemple : étude d’une molécule d’eau dans l’eau Première étude DFT/MM d’une molécule en solution Tuñón et al, J. Mol. Model., 1, 196 (1995) Modèle 1 molécule QM + 128 molécules MM (molécule QM rigide, géométrie =TIP3P) Niveau de calcul QM : DFT(BP) O(7111/411/1) H(41/1) MM : champ de force TIP3P Interactions QM/MM vdW = TIP3P Simulations MD Conditions périodiques Intégration 1 fs Durée : 70 ps après thermalisation ensemble NVE

36. Exemple : étude d’une molécule d’eau dans l’eau (c) O(QM) - H (MM) (d) Q(MM) - H (QM)

37. Exemple : étude d’une molécule d’eau dans l’eau Internal energy in kcal/mol Ei = 1/2(EQM+QM/MM-E°QM) Experimental -9.92 Theoretical Standard TIP3P MD -9.86 0.03 DFT/MM MD -9.72 2.1 Induced dipole moment « Experimental »= 0.75 D (lower limit) Theoretical AIMD Car-Parrinello = 1.08 D DFT/MM MD = 0.82 D Chalmet et al, JCP 2001

38. Exemple : étude d’une molécule d’eau dans l’eau Moment dipolaire instantané dans la simulation QM/MM MD

39. Exemple : étude d’une molécule d’eau dans l’eau Comparaison de différents modèles Modèle 1: Continuum 1 molécule QM + continuum e=80 Modèle 2: QM/MM/Continuum 1 molécule QM + 4 molécules MM TIP3P + continuum Modèle 3: QM/MM MD 1 molécule QM + 215 molécules MM TIP3P + conditions périodiques + simulation MD 25 ps.

40. Exemple : étude d’une molécule d’eau dans l’eau Potentiel électrostatique créé par le solvant QM/MM MD QM/MM/Continuum Continuum Moment dipolaire induit Dm « Expérimental » 0.75 D Calculs AIMD Car-Parrinello : 1.08 D Continuum 0.60 D QM/MM/Continuum 0.80 D QM/MM MD 0.82 D

41. II. Potentiels d’interaction. • Nécessité d ’une paramétrisation • Stratégies de paramétrisation • Au delà du potentiel de Lennard-Jones • Termes électrostatiques dans les méthodes semiempiriques

42. Nécessité d ’une paramétrisation Pour un niveau QM et un champ de force MM donnés, seul l’hamiltonien non électrostatique HvdWQM/MM n’est pas complètement défini. s MM sites n QM nuclei n , n ?

43. Nécessité d ’une paramétrisation Molécule d ’eau dans l ’eau : simulation DFT/TIP3P MD utilisant des paramètres de Lennard-Jones standard TIP3P pour les interactions QM/MM gOO(R)

44. Nécessité d ’une paramétrisation Molécule d ’eau dans l ’eau : simulation AM1/TIP3P MD utilisant des paramètres de Lennard-Jones standard TIP3P pour les interactions QM/MM gOO(R) expérimentale AM1/TIP3P

45. Nécessité d ’une paramétrisation Donc : Une paramétrisation est en général nécessaire Les paramètres n , n devraientdépendre de la méthode QM utilisée Pour les approches QM ab initio ou DFT, les paramètres peuvent dépendre de la base d ’orbitales Dans un traitement QM semiempirique, les paramètres peuvent dépendre de la méthode (AM1, PM3, ....)

46. Stratégies de paramétrisation • Paramétrisation : • Que veut-on reproduire ? • On peut envisager plusieurs cas : • Type de données : reproduire des calculs ab initio poussés, des calculs au même niveau QM ou encore des données expérimentales • Type de propriété : reproduire la structure (ou d’autres propriétés, ex surface d’énergie potentielle) d’un complexe, d’un liquide, d’une solution, etc

47. Stratégies de paramétrisation • Paramétrisation à partir de calculs « ab initio » sur des complexes de référence • Recette : • calculer la surface d’énergie potentielle ab initio du complexe AB • calculer l’énergie d’interaction QM/MM électrostatique • la différence représente l’interaction non-électrostatique • faire un « fit » de cette énergie à partir de l’expression de Lennard-Jones, c’est-à-dire, minimiser la quantité :

48. Stratégies de paramétrisation Paramétrisation des interactions non électrostatiques NH3 (DFT) - H2O (TIP3P) QMMM QMMM

49. QM/MM, vdW TIP3P classical experimental QM/MM vdWdimer optimized Stratégies de paramétrisation Paramétrisation pour des simulations en phase condensée. la stratégie précédente n’est pas valable : effets non-additifs ! Exemple Tu &Laaksonen (JCP 1999) Eau liquide Calculs QM/MM au niveau HF/6-311G(d,p)/TIP3P => Résultats très mauvais lorsque la simulation est faite en utilisant les paramètres LJ optimisés pour le dimère de l’eau. Moment dipolaire induit : 0.35 D (exp 0.75 D) Energie interne -6.74 kcal/mol (exp -9.92 kcal/mol)

50. Stratégies de paramétrisation • Paramétrisation pour des simulations en phase condensée • Martín et al, Chem. Phys. 284 (2001) 607 • La paramétrisation doit se faire en ajustant une propriété du liquide (ou de la solution) • => processus itératif (plusieurs simulations QM/MM !) • On peut remplacer les simulations QM/MM par des simulations ASEP, moins coûteuses, (ou par une simulation classique avec des charges atomiques appropriées) • Les différences avec les paramètres du champ de force classique peuvent être significatives O(solute)-O(solvent) radial distribution function of H2O Experimental (full line) ASEP/MD (dashed line) DFT/MM/MD (dotted line).