Résultats et discussion
L’événement d’autoionisation a été étudié à l’aide de simulations ab initio RETIS, comme décrit dans Matériaux et méthodes. Pour les simulations RETIS, nous avons utilisé un paramètre d’ordre de distance géométrique relativement simple, λ, comme illustré sur la figure 1. Lorsque le système est constitué uniquement d’espèces H2O, λ est la plus grande distance de liaison covalente O-H, et lorsque le système contient des espèces OH- et H3O+, λ est pris comme la plus courte distance entre l’oxygène dans OH- et les atomes d’hydrogène dans H3O+. Dans la suite, nous désignons l’atome d’oxygène utilisé pour le paramètre d’ordre par Oλ. Le type d’espèce (OH-, H2O, ou H3O+) a été identifié en attribuant à chaque hydrogène une liaison simple le reliant à l’oxygène le plus proche. Notez que la définition du paramètre d’ordre ne nécessite pas de seuil pour définir une liaison chimique et ne contraint pas le paramètre d’ordre à des molécules d’eau spécifiques pour la durée de la simulation. Cela signifie que nous calculons le taux de dissociation de toute molécule d’eau dans le système au lieu d’une seule liaison O-H ciblée ou d’une seule molécule d’eau.
Conditions d’initiation et variables collectives locales. (A) Distributions réactives (rλc,λr(ξ)) et non réactives (uλc,λr(ξ)) pour ξ={w4,na} et λc=1,16 Å et λr=2,0 Å. Pour des raisons de visualisation, les distributions représentées sont normalisées . Les encarts supérieurs et droits montrent les projections unidimensionnelles des distributions. Une séparation nette des deux distributions peut être observée le long de la coordonnée w4, indiquant que les trajectoires réactives sont plus comprimées par rapport aux trajectoires non réactives. En outre, l’atome d’oxygène utilisé dans le calcul du paramètre d’ordre (Oλ) accepte en moyenne un plus grand nombre de liaisons hydrogène dans les trajectoires réactives, par rapport aux trajectoires non réactives. (B) Instantané illustratif d’une trajectoire réactive où Oλ est représenté en bleu. Les quatre atomes d’oxygène environnants qui sont utilisés pour le calcul du paramètre d’ordre tétraédrique q sont représentés en orange. Le fil d’eau est mis en évidence par une ligne jaune (et des sphères transparentes grises) et le paramètre d’angle qcos est indiqué. Dans cet instantané, le fil d’eau est comprimé, q présente une déviation par rapport à une structure tétraédrique, qcos indique que trois atomes d’oxygène sont alignés dans le fil, et Oλ accepte trois liaisons hydrogène et en donne une (représentée par des lignes vertes).
Si nous considérons la coordonnée q, nous observons que rλc,λr est décalé vers des valeurs q plus faibles par rapport à uλc,λr, ce qui indique qu’une distorsion d’un arrangement tétraédrique autour de l’espèce d’eau dissociée peut également initier l’événement. Cette découverte est quelque peu surprenante car dans certaines autres réactions chimiques en phase aqueuse, l’effet inverse a été trouvé (31). Des conclusions similaires peuvent être tirées pour la distribution de ξ=(w4,qcos). Ici, il y a un pic le long de la coordonnée qcos pour la distribution réactive plus proche d’un arrangement linéaire des molécules d’eau. Dans la Fig. 4B, nous montrons un instantané représentatif, obtenu au début (après 3 fs) d’une trajectoire réactive. Globalement, les résultats présentés dans la Fig. 3 indiquent qu’une compression du fil d’eau (mesurée par w4) et une hypercoordination (mesurée par na) ou une distorsion (mesurée par q et qcos) sont des conditions d’initiation nécessaires à l’autoionisation. Cependant, ce ne sont pas des conditions suffisantes, comme le montrent les valeurs de TAλc,λr dans la figure 3B : encore 60 % des trajectoires démarrant dans la plage de paramètres ξ idéale ne parviennent pas à établir un saut de protons concerté.
L’apprentissage machine (ML) appliqué aux données d’échantillonnage de chemin (33, 34) est une approche prometteuse pour trouver des variables collectives importantes qui peuvent facilement être manquées par l’intuition humaine. Pour explorer cette possibilité, nous avons construit des modèles ML pour prédire le résultat des trajectoires étant donné l’état du système d’eau au début des trajectoires. Nous nous concentrons sur la même plage que dans l’analyse du pouvoir prédictif et nous utilisons l’état du système, lorsque λ>1,15 Å est atteint pour la première fois, pour prédire le résultat. Nous avons utilisé plusieurs techniques ML dans lesquelles chaque ensemble de chemins impairs a été inclus dans la calibration et les ensembles de chemins pairs ont été utilisés pour l’ensemble de test. Une autre répartition dans laquelle les données de chaque ensemble de chemins étaient divisées en deux parties égales a donné des résultats similaires. De plus, comme les distributions fortement asymétriques sont difficiles à traiter avec la méthode ML, nous avons également omis de repondérer les ensembles de données avec les poids statistiques des ensembles de chemins correspondants. Cependant, nous avons appliqué les techniques ML comme une approche qualitative pour trouver de nouveaux paramètres qui pourraient être testés quantitativement dans le cadre de la méthode de puissance prédictive (19).
En outre, pour éviter un risque potentiel de surinterprétation, nous avons choisi de restreindre la complexité du processus de décision ML et avons imposé un maximum de quatre paramètres d’ordre lors du calcul de TAλc,λr. Par exemple, d’excellentes performances prédictives (>90%) ont été obtenues en utilisant les machines à gradient-boosting basées sur un ensemble (35, 36). Cependant, l’interprétation du modèle est problématique puisqu’un ensemble de 100-150 arbres de décision profonds (ajoutés dans une séquence) est utilisé. Bien que la performance soit améliorée, le risque de surajustement avec des corrélations accidentelles augmente. Nous nous sommes donc limités aux modèles de décision à arbre unique basés sur les arbres de décision de classification et de régression (CART) (20). La restriction à quatre paramètres d’ordre pour la fonction TAλc,λr repose sur des raisons similaires. L’ajout de plus de paramètres donne des matrices plus clairsemées représentant les distributions réactives/non réactives et, par conséquent, l’intégration numérique pour calculer le chevauchement entre ces distributions devient très sensible à la taille de la case et pourrait sous-estimer le chevauchement en raison de cases vides par des statistiques insuffisantes.
Nous avons considéré 138 variables collectives constituées des distances oxygène-oxygène ; des distances oxygène-hydrogène pour les molécules d’eau initialement liées ; de tous les angles formés par Oλ et ses quatre voisins les plus proches de l’oxygène ; et des paramètres d’ordre de Steinhardt des ordres 3, 4 et 6 (32) (voir Matériaux et Méthodes pour plus de détails). En outre, les paramètres d’ordre déjà considérés ont été ajoutés. La figure 5A montre l’arbre de décision résultant. Il est remarquable que, parmi tous les paramètres d’entrée, le paramètre w4 soit à la fois au sommet de l’arbre de décision et la variable la plus importante, comme le montre la réduction de l’erreur de classification attribuée à chaque variable à chaque division de l’arbre de décision (20) (annexe SI, figure S9). L’ordre tétraédrique et le nombre de liaisons hydrogène acceptées apparaissent également dans l’arbre de décision. Pour décrire le premier effet, l’approche ML a donné la priorité au paramètre d’ordre q4 de Steinhardt par rapport au paramètre q similaire que nous utilisions précédemment. Certaines distances qui apparaissent également dans l’arbre de décision comme d25, la distance entre Oλ et son 25ème oxygène le plus proche, sont très probablement dues à des corrélations accidentelles causées par la taille limitée de l’ensemble de données. Ceci est vérifié en inspectant l’importance de cette variable : d25 ne figure pas parmi les 20 variables les plus importantes (Annexe SI, Fig. S9), et, en fait, d’autres variables similaires (par exemple, d24) sont classées plus haut, bien qu’avec une faible importance. Un paramètre plus important et intuitif qui est suggéré par l’approche ML est λ2, la distance OH entre l’oxygène le plus proche de Oλ et son hydrogène avec la plus grande liaison intramoléculaire. Le recalcul de la capacité de prédiction en utilisant les paramètres de l’arbre ML (Fig. 5B) n’a pas donné des performances plus élevées que la combinaison w4, q, na, et qcos, mais doit être conçu comme aussi bon, compte tenu des incertitudes statistiques.
Résultats de l’analyse par apprentissage machine. (A) Arbre de classification et de régression pour prédire le résultat des trajectoires initiées. Ici, nous avons considéré plusieurs variables collectives supplémentaires (description dans Matériaux et méthodes), mais seul un petit sous-ensemble est finalement nécessaire pour construire l’arbre : w4, q4, λ2 (la longueur de la liaison hydrogène étirée dans la molécule d’eau la plus proche de l’espèce Oλ), di (la distance de Oλ au ième oxygène le plus proche), et d¯i (la distance moyenne considérant les i oxygènes les plus proches). La notation des nœuds est expliquée avec le nœud autonome dans le coin supérieur gauche. Cet arbre prédit que les trajectoires seront réactives, c’est-à-dire atteignant un λ≥2, ou non réactives sur la base des variables collectives obtenues au cadre des trajectoires lorsque λ est d’abord ≥1,15. Les nœuds prédisant des trajectoires réactives sont colorés en bleu (classe 1) tandis que les nœuds prédisant des trajectoires non réactives sont colorés en vert (classe 0). Notez que les pourcentages au bas des carrés ne reflètent pas les fractions physiquement correctes puisque les ensembles de trajectoires n’ont pas été repondérés en utilisant leurs poids statistiques. Les règles sont des représentations textuelles de la traversée de l’arbre ; par exemple, la règle 5 (qui prédit des trajectoires réactives) peut être exprimée comme w4≥7,6 et λ2≥1,1. Ces règles donnent différentes conditions d’initiation, et elles sont répertoriées dans l’annexe SI, tableau S1, pour la rangée inférieure de nœuds. (B) Le pouvoir prédictif et la probabilité de croisement en fonction de λr pour λc=1,16 Å et différentes combinaisons de variables collectives. Nous comparons ici le pouvoir prédictif en utilisant les variables collectives que nous avons identifiées avec les variables marquées comme importantes par l’analyse par apprentissage automatique. (C) Distributions réactives (rλc,λr(ξ)) et non réactives (uλc,λr(ξ)) pour ξ={λ2,d¯2} et λc=1.16 Å et λr=2.0 Å. Pour des raisons de visualisation, les distributions représentées sont normalisées. Les encarts supérieurs et droits montrent les projections unidimensionnelles des distributions.