La dynamique moléculaire
Introduction
Depuis les premières de simulations de liquides simples effectuées par Metropolis et ses collaborateurs (Metropolis et al., 1953) et Wood et Parker (Wood et Parker, 1957), les techniques de « physique et de chimie computationelle » ont connu un essor considérable et se révèlent des outils de choix, de nos jours, pour étudier au niveau atomique les propriétés des phases condensées (comme, nous le verrons plus loin). Depuis une quinzaine d’années, avec l’augmentation des moyens informatiques et l’existence d’algorithmes de simulation de plus en plus performants, de nombreuses techniques de physique computationnelle, méthodes stochastiques (Monte Carlo (Frenkel et Smit, 2002), Brownienne (Ermak et McCammon, 1978)) ou déterministes (telle que la dynamique moléculaire (DM) (van Gunsteren et al., 2006) se révèlent de plus en plus présentes et accessibles dans les laboratoires pour interpréter les expériences et peuvent être confrontées avec les résultats expérimentaux. Dans ce qui suit, nous parlerons uniquement de la dynamique moléculaire (DM). Nous verrons comment, à partir d’une structure atomique initiale et d’une fonction d’énergie potentielle, il est possible de modéliser l’évolution temporelle d’un système moléculaire en solution et de reproduire des propriétés structurales de systèmes de plus en plus complexes contenant plusieurs dizaines à plusieurs centaines de milliers d’atomes (avec une limite, à l’heure actuelle, de l’ordre de ∼2.7.106 atomes13). Nous discuterons ensuite quelques techniques de dynamiques moléculaires disponibles dans le code de dynamique moléculaire ORAC, que nous avons utilisé pour ce travail de thèse. Enfin, nous discuterons les différentes approches, permettant de contrôler les paramètres physiques comme la température ou la pression.Le comportement d’un système de particules composé de noyaux, d’atomes, d’électrons est en principe décrit par les lois de la mécanique quantique et devrait donc être étudié par des méthodes ab-initio. Les moyens informatiques actuels permettent d’étudier par des méthodes abinitio, que des systèmes contenant, au plus, quelques centaines d’atomes, le temps de calcul étant proportionnel au cube du nombre d’électrons. Par conséquent, pour des systèmes de grande taille, comme ceux étudiés dans cette thèse, l’approche la plus viable est d’utiliser la dynamique (ou mécanique) moléculaire classique. Cette approche est acceptable, si on fait quelques approximations. La première est de considérer que le mouvement des électrons est beaucoup plus rapide que celui des noyaux (approximation de Born-Oppenheimer) et donc de considérer que l’évolution d’un système moléculaire peut être uniquement décrit par l’énergie du système, en fonction des coordonnées des noyaux et de quelques paramètres (comme l’étirement des liaisons, l’ouverture et la fermeture des angles de valence, la rotation autour des liaisons etc.). Le système moléculaire est donc décrit à l’aide de la mécanique classique et des lois de Newton. A cause de ces approximations, la DM ne peut pas être utilisée directement pour étudier les réactions chimiques, car le système garde tout le temps son identité chimique. Dans ce cas, il est possible d’utiliser une approche « mixte » mêlant les principes de la mécanique quantique et de la mécanique moléculaire ou QM/MM (Gao, 1996). En utilisant cette approche, la partie « réaction chimique » du système est traitée avec la mécanique quantique alors que l’evolution du reste du système est traitée avec la 13 Comme les simulations explicites récentes de ribosome (2.64 x 106 atomes) (Sanbonmatsu, 2005) ou du virus de la mosaïque du Tabac ~1 x 106 atomes) (Freddolino, 2006). II. La dynamique moléculaire mécanique moléculaire classique.
Hamiltonien et fonction d’interaction du système
La dynamique moléculaire est utilisable pour l’étude d’un système atomique ou moléculaire, si une modélisation préalable du système a été faite. Cette condition implique qu’on ait défini le potentiel qui agit entre les molécules et les degrés de libertés (DDL) du système. Tous les degrés de liberté représentés explicitement doivent être traités de façon classique. Ensuite, il faut choisir la fonction qui décrit les interactions en fonction des degrés de liberté explicites du modèle. Enfin, il faut sélectionner une méthode de simulation pour échantillonner statistiquement les configurations du système. Puisque que la DM est basée sur la mécanique classique, toutes les méthodes qui la concernent sont dérivées de l’Hamiltonien H du système, qui a la forme : H K (r,p;m,s) (p;m) + U(r;s) = (II.1) où r, p, m et s sont les coordonnées, les impulsions, la masse et les paramètres du champ de force, respectivement. L’énergie cinétique, 2 n n p i 1 2 (p;m) = = m vi i i=1 i=1 2m 2 i K ∑ ∑ (II.2) dépend uniquement des masses atomiques mi et des vitesses vi . Elle est indépendante des coordonnées r de la particule. L’énergie potentielle qui décrit l’énergie d’interaction en fonction des coordonnées r des particules : U(r;s) = U(r ,r ,…,r ,s ,s ,…,s ) 1 2 1 2 n n (II.3) La fonction d’interaction est une somme de termes qui modélise les différents types d’interactions entre les atomes d’une molécule ou d’un ensemble d’atomes. Lors des simulations de dynamique moléculaire, les trajectoires sont générées à partir des forces agissant sur les atomes, qui sont définis par : U f = – (r , r ,…, r ,s ,s ,…,s ) i r 1 2 1 2 M M i ∂ ∂ (II.4)
La fonction d’énergie potentielle ou champs de force
Les informations sur les interactions entre tous les atomes d’un système moléculaire sont donc contenues dans une fonction d’énergie potentielle U(r;s) 14. Les paramètres de cette fonction sont spécifiques du modèle moléculaire étudié et sont une des clefs de voûte des calculs effectués en DM. Pour un système moléculaire, l’expression de l’énergie potentielle, en fonction des positions r des N particules (usuellement d’atomes), U(r N ), est donnée généralement par la somme de deux termes liés aux interactions intramoléculaires, Uintra et intermoléculaires, Uinter : N U(r ) = U + U intra inter (II.5) 14 On simplifiera par la suite U(r N ). II. La dynamique moléculaire – 43 – Ces termes seront détaillés dans les sous-sections suivantes. II.4.1 Les interactions intramoléculaires ou « liantes » Le terme « intra » désigne les interactions mettant en jeu les atomes reliés par une ou plusieurs liaisons covalentes (au maximum 3). Dans les champs de force les plus courants, comme AMBER (Case et al., 2005) ou CHARMM (MacKerell et al., 1998b), les interactions intramoléculaires sont décrites par la somme des contributions suivantes : U = U + U + U + U propre liaison angle impropre (II.6) où Uliaison, Uangle, Upropre et Uimpropre représentent les contributions à l’énergie potentielle d’interaction des liaisons covalentes, des angles de valences et des torsions propres et impropres. Ces potentiels seront détaillés dans les paragraphes suivants et une représentation est donnée sur la Fig. II.3 II
Les potentiels d’interaction de liaison et d’angle
La modélisation des potentiels d’interaction de liaison entre deux atomes successifs ou « liaisons » et de torsion des angles de valence entre 2 liaisons covalentes successives ou « angle » est décrite par des fonctions harmoniques de forme : 2 U = K (r – r ) liaison b 0 liaison ∑ (II.7) 2 U = K (θ – θ ) angle θ 0 angle ∑ (II.8) avec Kb et Kθ, les constantes de force associées aux liaisons et aux angles de valence, respectivement, et r0 et θ0 les valeurs des liaisons et des angles de valences de référence (i.e. à l’équilibre), respectivement. Dans plusieurs cas, les liaisons covalentes, dont les fréquences de vibrations sont très élevées, comme celles impliquant les atomes d’hydrogène, sont modélisées à l’aide de liaisons rigides à leurs valeurs d’équilibre (i.e. r=ro), par la méthode de contraintes (Ciccotti et Ryckaert., 1986).
Les potentiels d’interaction de dièdre
La variation de l’énergie liée à la déformation des angles dièdres pour 4 atomes successifs fait intervenir un potentiel de torsion Upropre. La description de Upropre est faite à partir d’un développement en série de Fourier : = + φ − γ ∑ [ ] φ U K 1 cos(n ) propre propre Avec n=1, 2, 3… (II.9) avec Kφ, n et γ sont respectivement la barrière de torsions, la périodicité et la phase respectivement. Le choix de n dépend du type de liaison dièdre qu’on veut décrire (Figure. II-1). Le premier ordre de la série de Fourier, n=1, agit sur la stabilisation des formes de conformation trans ou cis, le deuxième ordre, n=2, sur la stabilisation planaire des doubles liaisons entre atomes de type sp2 et le troisième ordre, n=3, sur la stabilisation des formes de conformation décalée ou éclipsée. Figure II-1 : Ordres de n pour la modélisation des dièdres. Ainsi, si on prend le cas de deux carbones tétraédriques sp3 reliés par une liaison covalente (R-CH2-CH3), on a une fonction Upropre avec n = 3 et γ = 0o . Pour le propane (R = H), on aura 3 seuils de rotation avec 3 minima aux angles 60o , 180o et 300o correspondant aux conformations, dites décalées, et 3 maxima aux angles 0o , 120o et 240o , correspondant aux conformations, dites éclipsées des différents rotamers (Voir la Fig. II-2). Figure II-2 : Variation de l’énergie de torsion du dièdre R-CH2-CH3 du propane en fonction d eφ, n et γ. Parfois, il est nécessaire que certains groupes atomiques d’une molécule soient fixés dans une certaine conformation particulière vis à vis de ces voisins, par exemple, le caractère « plan » de l’atome de carbone hybridé sp2 (i.e. CH2=) (voir Fig. II-3). On ajoute alors à la somme II.6 un terme d’interaction dit impropre, Uimpropre. Ce type d’interaction est modélisé, dans le champ de force CHARMM, par une fonction harmonique avec une constante de force, Kω et une valeur de référence, ω0 : 2 U = K ( ω ω-ω ) impropre o impropre ∑ (II.10) Dans le cas du champs de force AMBER (Case et al., 2005), le traitement des angles dièdres (propre et impropre) se fait à l’aide d’une seule série de Fourier, ayant la même forme que II.9. Généralement, pour modéliser l’énergie liée à la déformation des dièdres, le développement en série de II.9 ne comporte qu’un seul terme, mais il n’est pas rare, dans certains cas qu’il en faille plus (c’est le cas, par exemple du dièdre O-C-C-O dans le champ de force AMBER où deux termes sont nécessaires, pour que ce dièdre adopte sa conformation gauche (Leach, 2001)). L’ensemble de ces potentiels d’interaction est représentés sur la Fig. II-3. Figure II-3 : Illustration des différents termes de la fonction d’énergie potentielle U(rN ). Les schémas (a-d) montrent les termes de valence du champ de force. Le terme (d) est le terme dit de torsion « impropre », garantissant que l’atome central A reste dans le même plan que les atomes B, C et D. La Figure (e) représente les interactions non liantes dites «1-4 » qui seront décrites dans la partie II.4.2. Les valeurs de référence pour les longueurs de liaison et les angles et γ, peuvent être déduites à partir des données expérimentales de la géométrie d’une molécule, obtenue par diffraction des rayons X (RX), de neutrons. Quant aux valeurs des constantes de force Kb, Kθ, Kω, elles sont calculées soit à partir des modes vibrationnels, obtenus par spectroscopie (infrarouge, Raman) ou par des calculs de modes normaux (Leach, 2001).
Introduction générale |
