Signification
Les réactions acide-base font partie des processus chimiques les plus importants. Pourtant, il nous manque un moyen simple de décrire cette classe de réactions en fonction des coordonnées atomiques. En effet, une fois dissous dans l’eau, H+ et son anion conjugué OH- ont une structure hautement fluxionnelle difficile à cerner. Nous résolvons ici ce problème en adoptant le point de vue de la description des réactions acide-base comme un équilibre entre le soluté et l’ensemble du solvant. Cela permet d’identifier des descripteurs généralement applicables. En conséquence, il est maintenant possible d’effectuer une simulation quantitative par échantillonnage amélioré de la réaction acide-base dans l’eau et dans d’autres environnements tels que les cavités des zéolithes ou aux surfaces.
Abstract
Les réactions acide-base sont omniprésentes dans la nature. La compréhension de leurs mécanismes est cruciale dans de nombreux domaines, de la biochimie à la catalyse industrielle. Malheureusement, les expériences ne donnent qu’une information limitée sans beaucoup de compréhension du comportement moléculaire. Les simulations atomistiques pourraient compléter les expériences et apporter un éclairage précieux sur les mécanismes microscopiques. Cependant, les importantes barrières d’énergie libre liées à la dissociation des protons rendent obligatoire l’utilisation de méthodes d’échantillonnage améliorées. Nous réalisons ici une simulation ab initio de dynamique moléculaire (MD) et améliorons l’échantillonnage à l’aide de la métadynamique. Ceci a été rendu possible par l’introduction de descripteurs ou de variables collectives (CV) qui sont basés sur une vision conceptuellement différente des équilibres acide-base. Nous testons avec succès notre approche sur trois solutions aqueuses différentes d’acide acétique, d’ammoniac et de bicarbonate. Celles-ci sont représentatives du comportement acide, basique et amphotère.
- acide-base
- métadynamique
- variables collectives
- échantillonnage renforcé
Les réactions acide-base jouent un rôle clé dans de nombreuses branches de la chimie. Les réactions de complexation inorganique, le repliement des protéines, les processus enzymatiques, la polymérisation, les réactions catalytiques et de nombreuses autres transformations dans différents domaines sont sensibles aux changements de pH. Comprendre le rôle du pH dans ces réactions implique d’avoir un contrôle sur leur réactivité et leur cinétique.
L’importance cruciale du pH a stimulé la collecte d’une grande quantité de données sur les équilibres acide-base. Celles-ci sont généralement mesurées dans les phases gazeuses et condensées, à l’aide de techniques spectroscopiques et potentiométriques. Cependant, il existe des limites pratiques à la précision de ces méthodes, en particulier dans les phases condensées (1). En outre, il est très difficile d’extraire des données expérimentales une image microscopique des processus impliqués. Il n’est donc pas surprenant que l’équilibre acide-base ait fait l’objet d’une intense activité théorique (1⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓-12).
L’acidité d’une espèce chimique dans l’eau peut être exprimée en termes de pKa, le logarithme négatif de la constante de dissociation de l’acide. Il existe deux façons de calculer ces valeurs, l’une statique et l’autre dynamique.
L’approche la plus standard est l’approche statique dans laquelle les énergies libres en phase solution, et par conséquent les pKas, sont obtenues en fermant un cycle de Born-Haber composé des énergies libres en phase gazeuse et en solvatation (1, 3⇓⇓⇓-7). Bien qu’extrêmement performante dans de nombreux cas, l’approche statique présente certaines limites. Un modèle de solvatation doit être choisi et les modèles de solvant continuum ont une précision limitée. Ceci est particulièrement vrai dans des systèmes comme les zéolithes ou les protéines caractérisées par des cavités irrégulières dans lesquelles une description implicite du solvant est un défi. Il est évident qu’une telle approche ne permet pas d’obtenir des informations dynamiques. En outre, il peut y avoir des réactions compétitives qui ne peuvent pas être prises en compte si elles ne sont pas explicitement incluses dans le modèle.
En principe, ces limitations pourraient être levées dans une approche plus dynamique basée sur des simulations de dynamique moléculaire (MD) dans lesquelles les molécules de solvant sont traitées explicitement. Si l’on disposait d’un temps informatique illimité, de telles simulations permettraient d’explorer toutes les voies possibles et d’attribuer le poids statistique relatif aux différents états. Malheureusement, la présence de goulots d’étranglement cinétiques empêche cette possibilité en piégeant le système dans des états métastables, puisque les différents états de protonation sont séparés par de grandes barrières. En outre, dans les réactions acide-base, des liaisons chimiques sont formées et brisées. Cela nécessite l’utilisation de la MD ab initio dans laquelle les forces interatomiques sont calculées à la volée à partir des théories de structure électronique. Cela rend le calcul plus coûteux et réduit encore l’échelle de temps qui peut être explorée.
Pour surmonter cette difficulté, l’utilisation de méthodes d’échantillonnage améliorées (13) qui accélèrent l’exploration de l’espace configurationnel devient obligatoire. Une classe très populaire de méthodes d’échantillonnage améliorées est basée sur l’identification des degrés de liberté qui sont impliqués dans la réaction lente d’intérêt. Ces degrés de liberté sont habituellement appelés variables collectives (CV) et sont exprimés comme des fonctions explicites des coordonnées atomiques R. L’échantillonnage est ensuite amélioré en ajoutant un biais qui est une fonction des CV choisies (14⇓-16). En outre, la conception d’un ensemble approprié de bons CV a également une signification plus profonde. Les CVs réussis capturent de manière condensée la physique du problème, identifient ses degrés de liberté lents, et conduisent à une description modéliste utile du processus.
Dans les réactions chimiques standard, ceci est relativement simple puisque des structures bien définies peuvent être assignées aux réactifs et aux produits (17⇓-19). Ce n’est pas le cas pour les réactions acide-base dans lesquelles un proton est ajouté ou soustrait au soluté. Une fois que ce processus a eu lieu, les ions eau (H+ ou/et OH-) sont solvatés et leur structure devient insaisissable. En effet, les ions eau peuvent rapidement diffuser dans le milieu par un mécanisme de Grotthuss (20). Ils deviennent fortement fluxionnels et l’identité des atomes participant à leur structure change continuellement. La nature de ces espèces est donc difficile à capturer dans une fonction analytique explicite de R. Cependant, étant donné la pertinence des réactions acide-base, de nombreuses tentatives ont été faites pour définir ces entités (8⇓⇓⇓⇓-12). Malheureusement, ces CV ont une nature ad hoc et, bien que réussissant dans tel ou tel cas, ne peuvent pas être appliquées de manière générale.
Pour construire des CV généraux et utiles, nous faisons deux étapes conceptuelles. La première consiste à considérer le processus acide-base comme une réaction impliquant seulement quelques fractions, à savoir l’ensemble du solvant et les résidus réagissant dans la molécule solvatée. Par exemple, lorsqu’il n’y a qu’un seul type de résidu dissociant, on pense à l’équilibre acide-base comme une réaction du typeA+H2NON⇌Bq0+H2N+q1ONq1,où N est le nombre de molécules d’eau, A et B sont une molécule acide-base générique en solution et son espèce conjuguée, respectivement, q0 et q1 sont des entiers qui peuvent prendre les valeurs +1 et -1 selon le comportement acide-base de l’espèce, et q1+q0=0.
Cela implique que nous ne regardons pas le solvant comme un ensemble de molécules qui sont en compétition pour réagir avec l’espèce acido-basique. Nous considérons plutôt le solvant dans sa globalité comme l’un des deux adduits. Ce point de vue est particulièrement pertinent dans les solvants polaires comme l’eau qui sont caractérisés par des réseaux hautement structurés. Dans ce cas, la présence d’un excès ou d’un déficit de protons modifie localement la structure du réseau et cette distorsion se propage le long de l’ensemble du réseau.
Depuis les tout débuts de Wicke et Eigen (21) et de Zundel et Metzger (22), les chercheurs se sont battus pour savoir combien de molécules devaient être incluses dans la définition de la perturbation (23⇓-25). Compte tenu de l’absence de paramètres physiques capables de donner une réponse claire et sans équivoque à cette question, l’idée de considérer le solvant dans son ensemble permet de contourner ce problème. Ainsi, le solvant n’est pas seulement un milieu ayant un rôle passif, mais il est regardé comme un ensemble de molécules qui contribuent collectivement à la formation de la paire acide-base conjuguée. Ce point de vue est beaucoup plus proche du point de vue original proposé par Brønsted et Lowry dans lequel la réaction peut être vue comme un simple échange d’un cation hydrogène entre une paire acide-base.
Pour que la réaction ait lieu, le centre de la perturbation doit s’éloigner du soluté. Ainsi, la deuxième étape importante consiste à contrôler le centre de la perturbation. En raison des mécanismes de type Grotthuss, la perturbation se déplace le long du réseau. Cela peut conduire à des définitions différentes du centre du défaut. Cependant, si nous tessellons l’espace entier en utilisant des polyèdres de Voronoï centrés sur les atomes d’oxygène de l’eau, nous pouvons assigner sans équivoque chaque atome d’hydrogène à un et un seul de ces polyèdres. Le site dont le polyèdre de Voronoï contient un nombre anormal de protons est pris comme centre de la perturbation (Fig. 1).
Deux exemples de partition de l’espace. (Gauche) Nous montrons une approche convective dans laquelle la distance de l’atome d’oxygène est utilisée pour définir son environnement. On peut clairement voir des superpositions artificielles. (Droite) La tessellation de Voronoï ne souffre pas de ces défauts.
Ce point de vue confère à la méthode un caractère très général, la rendant applicable à tout système acide-base, sans qu’il soit nécessaire de fixer au préalable les couples réactifs. Ainsi, il est possible d’explorer tous les états de protonation pertinents même dans des systèmes composés de plus d’une paire acide-base.
Cette approche générale permet de définir des CV sans avoir à imposer des structures spécifiques ou à sélectionner l’identité des atomes impliqués. Nous testons notre méthode en effectuant des simulations métadynamiques dans un cas d’acide faible (acide acétique), dans une base faible (ammoniac) et dans une espèce amphotère (bicarbonate) choisis comme points de repère en raison de leur force comparable, mais de leur comportement acide-base différent.
Méthodes
Comme discuté ci-dessus, nous introduisons deux CV, un lié à l’état de protonation et un autre qui localise les défauts de charge et mesure leur distance relative. Ces deux CVs ont besoin d’une définition robuste pour assigner les atomes d’hydrogène au site acide-base respectif. Pour obtenir ce résultat, nous divisons l’espace entier en polyèdres de Voronoï centrés sur les sites acide-base i situés à Ri. Ces sites comprennent tous les atomes capables de rompre et de former des liaisons avec un proton acide. La partition standard de l’espace de Voronoï est décrite par un ensemble de fonctions d’indexation wi(r) centrées sur les différents Ris, telles que wi(r)=1 si le ième atome est le plus proche de r et wi(r) = 0 sinon. Pour leur utilisation dans les méthodes d’échantillonnage amélioré, les CV doivent être différentiables. À cet effet, nous introduisons une version lisse des fonctions d’indexation, wis(r). Celles-ci sont définies en utilisant des fonctions softmaxwis(r)=e-λ|Ri-r|∑me-λ|Rm-r|,où i et m couvrent tous les sites acide-base et λ contrôle la pente avec laquelle les courbes se réduisent à 0, c’est-à-dire la sélectivité de la fonction. Avec un choix approprié de λ, cette définition permet d’obtenir le résultat souhaité, comme le montre la figure 2. De cette manière, un atome d’hydrogène dans une position Rj est assigné au polyèdre centré sur le site i avec le poids wi(Rj). Ensuite, le nombre total d’atomes d’hydrogène affectés au ie site acide-base estWi=∑j∈Hwis(Rj),où la sommation sur j court sur tous les atomes d’hydrogène.
Tessellation lisse d’un espace 2D avec des cellules centrées sur les trois atomes d’oxygène de la molécule d’eau. Les régions bleues plates représentent la portion d’espace dans laquelle la fonction prend une valeur de 1 et les régions jaunes représentent les frontières entre les cellules. Cette surface a été obtenue avec une valeur de λ=4.
On peut associer à chaque site acide-base une valeur de référence Wi0 qui compte le nombre d’atomes d’hydrogène liés à l’état neutre. La différence entre la valeur instantanée des atomes d’hydrogène et celle de référence estδi=Wi-Wi0.Lorsqu’il est différent de zéro, δi signalera si le ie site a gagné ou perdu un proton. Dans le cas des atomes d’oxygène de l’eau, un ion hydronium a un δi=+1 alors qu’un ion hydroxyde a un δi=-1.
On regroupe ensuite les sites acido-basiques en espèces. Par exemple, dans le cas de l’acide aminé le plus simple, la glycine, en solution aqueuse, le nombre d’espèces Ns sera égal à 3. Tous les atomes d’oxygène de l’eau appartiennent à une espèce ; puis on compte dans une autre espèce les deux atomes d’oxygène carboxyliques, et enfin on considère comme troisième espèce l’atome d’azote du groupe amino.
Dans l’esprit de ce travail, nous comptons l’ensemble des excès ou des défauts de protons associés à chaque espèce:qk=∑i∈kδi.Cela implique que nous ne nous intéressons pas à l’identité spécifique du site réagi, mais au fait que la kième espèce dans sa totalité a augmenté (qk=+1), a diminué (qk=-1), ou n’a pas changé son nombre de protons. Si nous considérons un soluté avec une seule fraction réactive, alors chaque état possible du système peut être décrit par l’un des trois vecteurs 2D (0,0), (-1,1), ou (1,-1).
Dans le cas général, chaque état de protonation peut être décrit par un vecteur q→=(q0,q1,…qNs-1) de dimension égale au nombre de sites réactifs non équivalents, Ns. Une explication plus exhaustive est fournie dans l’annexe SI.
Pour être utilisés dans l’échantillonnage amélioré, ces vecteurs doivent être exprimés sous la forme d’une fonction scalaire f=f(q→) telle que, pour chaque q→ physiquement pertinent, f atteigne des valeurs capables de distinguer les différents états de protonation globaux. Il existe une infinité de façons de construire un scalaire à partir d’un vecteur. Le choix le plus simple est peut-être d’écrire f(q→)=X→⋅q→ et, pour distinguer les différents états de protonation, de choisir X→=(20,21,22,…2Ns-1).
Cela conduit à la définition suivante pour le CV qui est utilisé pour décrire l’état de protonation du système,sp=∑k=0Ns-12k⋅qk,où k sont les indices utilisés pour étiqueter les groupes de sites réactifs respectifs. Dans l’annexe SI, un exemple est élaboré en détail. Bien sûr, le CV est rendu continu par l’utilisation de Wi dans le calcul du δi nécessaire pour évaluer qk dans l’équation 5.
Le second CV est une sommation des distances entre tous les sites acido-basiques multipliées par leur charge partielle δi,sd=∑i,m>i-rim⋅δi⋅δm,où les indices i et m courent sur tous les sites acido-basiques appartenant à différents groupes k, et rim est la distance entre les deux atomes. De cette façon, seule la paire acide-base qui a échangé un proton donne une contribution différente de zéro. L’équation 7 n’est valable que lorsqu’une seule paire acide-base conjuguée est présente. Cependant, en raison de l’action du biais, il peut arriver occasionnellement que plusieurs paires acide-base soient formées. Pour éviter d’échantillonner ces événements très improbables, nous appliquons une restriction sur le nombre de paires. De plus amples détails sont fournis dans l’annexe SI.
Résultats
Nous avons appliqué notre méthode à trois solutions aqueuses d’acide acétique, d’ammoniac et de bicarbonate comme représentations d’un acide faible, d’une base faible et d’un composé amphotère, respectivement. Les configurations des trois simulations sont identiques, à l’exception de l’identité des molécules solvatées. Cela permet de s’assurer que le résultat reflète la chimie différente de ces trois systèmes et qu’il n’y a pas de biais dû à la condition initiale.
Chaque simulation des systèmes a été réalisée avec des simulations MD de Born-Oppenheimer combinées à une métadynamique bien tempérée (14, 26) en utilisant le paquet CP2K (27) patché avec PLUMED 2 (28) et une fonctionnelle fortement contrainte et convenablement normée (29) pour l’énergie xc, Exc. Voir l’annexe SI pour plus de détails.
Dans la figure 3, nous traçons les surfaces d’énergie libre (FESs) en fonction de sp et sd. Ces FES reproduisent de manière vivante le comportement attendu. Elles ont toutes un minimum à sp=0 qui correspond à l’état dans lequel aucune charge n’est présente dans le solvant. Dans le FES de l’acide acétique (Fig. 3A), un second minimum proche de sp=-1 reflète son comportement acide. En revanche, le FES de l’ammoniac (Fig. 3B) présente un second minimum proche de sp=1. Les formes des FES de l’ammoniac et de l’acide acétique sont approximativement liées par une symétrie miroir reflétant leur comportement contrasté. De même, le FES symétrique du bicarbonate (Fig. 3C) reflète son caractère amphotère.
(A-C) Surfaces d’énergie libre le long de sp et sd de l’acide acétique (A), de l’ammoniac (B) et du bicarbonate (C) en solution aqueuse. Les barres de couleur indiquent l’énergie libre exprimée en unités kJ⋅mol-1. Le CV sd est exprimé en angströms.
A mesure que la paire conjuguée est formée, sd commence à prendre des valeurs positives correspondant à la séparation et à la diffusion de la paire conjuguée. Par rapport à l’état non dissocié dans lequel seul sd=0 est autorisé, les états où une paire conjuguée est présente montrent une forme allongée des bassins le long de cette variable. Ceci est causé par le comportement diffusif des ions hydronium et hydroxyde en solution qui rend accessible une gamme continue de distances. De plus, le long de cette CV, nous pouvons observer une barrière autour de 1,5 correspondant à la rupture de la liaison covalente entre l’atome d’hydrogène et le site acide-base.
Conclusions
L’applicabilité générale de cette méthode à des systèmes de natures différentes est une étape importante franchie dans leur compréhension et leur description. Le schéma peut être étendu pour inclure les effets nucléaires quantiques avec l’utilisation de la dynamique moléculaire intégrale du chemin (30). Cela aurait une importance quantitative puisque, par exemple, les valeurs de pKa sont affectées par la deutération. De plus, l’absence d’hypothèses ou d’impositions sur les candidats réactifs ou les chemins de réaction permet d’étendre cette méthode à des systèmes de complexité croissante qui ne peuvent être traités par les méthodes traditionnelles. Des exemples de questions auxquelles il est maintenant possible de répondre sont les équilibres tautomères dans les processus biochimiques et le comportement acide dans les zéolites et à la surface des oxydes exposés à l’eau.
Remerciements
Les calculs ont été effectués sur le cluster ETH Euler et sur le cluster Mönch au Centre national suisse de supercalcul. Cette recherche a été soutenue par la subvention de l’Union européenne ERC-2014-AdG-670227/VARMET.
Publiée sous la licence PNAS.
.