Risultati e discussione
L’evento di autoionizzazione è stato studiato usando simulazioni RETIS ab initio come descritto in Materiali e metodi. Per le simulazioni RETIS, abbiamo usato un parametro di ordine di distanza geometrica relativamente semplice, λ, come illustrato in Fig. 1: quando il sistema consiste solo di specie H2O, λ è la più grande distanza di legame covalente O-H, e quando il sistema contiene specie OH- e H3O+, λ è preso come la distanza più breve tra l’ossigeno in OH- e gli atomi di idrogeno in H3O+. Nel seguito, ci riferiamo all’atomo di ossigeno usato per il parametro d’ordine come Oλ. Il tipo di specie (OH-, H2O, o H3O+) è stato identificato assegnando ad ogni idrogeno un singolo legame che lo collega all’ossigeno più vicino. Si noti che la definizione del parametro d’ordine non richiede una soglia per definire un legame chimico né vincola il parametro d’ordine a specifiche molecole d’acqua per tutta la durata della simulazione. Questo significa che calcoliamo il tasso di dissociazione di qualsiasi molecola d’acqua nel sistema invece di un singolo legame O-H mirato o di una molecola d’acqua.
Condizioni d’inizio e variabili collettive locali. (A) Distribuzioni reattive (rλc,λr(ξ)) e non reattive (uλc,λr(ξ)) per ξ={w4,na} e λc=1.16 Å e λr=2.0 Å. Per scopi di visualizzazione, le distribuzioni rappresentate sono normalizzate. Gli inserti in alto e a destra mostrano le proiezioni unidimensionali delle distribuzioni. Una chiara separazione delle due distribuzioni può essere vista lungo la coordinata w4, indicando che le traiettorie reattive sono più compresse rispetto a quelle non reattive. Inoltre, l’atomo di ossigeno usato nel calcolo del parametro di ordine (Oλ) accetta in media un numero maggiore di legami a idrogeno nelle traiettorie reattive, rispetto alle traiettorie non reattive. (B) Istantanea illustrativa di una traiettoria reattiva dove Oλ è mostrato in blu. I quattro atomi di ossigeno circostanti che sono usati per il calcolo del parametro d’ordine tetraedrico q sono mostrati in arancione. Il filo d’acqua è evidenziato con una linea gialla (e sfere grigie trasparenti) e il parametro angolare qcos è indicato. In questa istantanea, il filo d’acqua è compresso, q mostra una deviazione da una struttura tetraedrica, qcos indica che tre atomi di ossigeno sono allineati nel filo, e Oλ accetta tre legami idrogeno e ne dona uno (mostrato con linee verdi).
Se consideriamo la coordinata q, osserviamo che rλc,λr è spostata verso valori q più bassi rispetto a uλc,λr, il che indica che una distorsione da una disposizione tetraedrica intorno alla specie di acqua dissociante può anche iniziare l’evento. Questo risultato è in qualche modo sorprendente, poiché in alcune altre reazioni chimiche in fase acquosa è stato trovato l’effetto opposto (31). Conclusioni simili possono essere tratte per la distribuzione di ξ=(w4,qcos). Qui, c’è un picco lungo la coordinata qcos per la distribuzione reattiva più vicina a una disposizione lineare delle molecole d’acqua. In Fig. 4B mostriamo un’istantanea rappresentativa, ottenuta presto (dopo 3 fs) in una traiettoria reattiva. Nel complesso i risultati mostrati in Fig. 3 riportano che una compressione del filo d’acqua (misurata da w4) e una ipercoordinazione (misurata da na) o distorsione (misurata da q e qcos) sono condizioni di inizio necessarie per l’autoionizzazione. Tuttavia, queste non sono condizioni sufficienti, come mostrato dai valori di TAλc,λr in Fig. 3B: ancora il 60% delle traiettorie che partono all’interno dell’intervallo ideale dei parametri ξ non riescono a stabilire un salto protonico concertato.
Il machine learning (ML) applicato ai dati di path-sampling (33, 34) è un approccio promettente per trovare importanti variabili collettive che possono essere facilmente mancate dall’intuizione umana. Per esplorare questa possibilità, abbiamo costruito modelli ML per prevedere l’esito delle traiettorie dato lo stato del sistema idrico all’inizio delle traiettorie. Ci concentriamo sullo stesso intervallo dell’analisi del potere predittivo e usiamo lo stato del sistema, quando λ>1,15 Å viene raggiunto per la prima volta, per predire il risultato. Abbiamo usato diverse tecniche ML in cui ogni insieme di percorsi dispari è stato incluso nella calibrazione e gli insiemi di percorsi pari sono stati utilizzati per il set di test. Una divisione alternativa in cui i dati all’interno di ogni insieme di percorsi sono stati divisi equamente in due ha dato risultati simili. Inoltre, poiché le distribuzioni fortemente asimmetriche sono difficili da trattare con ML, abbiamo ulteriormente omesso la riponderazione dei set di dati con i pesi statistici degli insiemi di percorsi corrispondenti. Tuttavia, abbiamo applicato le tecniche ML come un approccio qualitativo per trovare nuovi parametri che potrebbero essere testati quantitativamente all’interno del metodo della potenza predittiva (19).
Inoltre, per evitare un potenziale rischio di sovrainterpretazione abbiamo scelto di limitare la complessità del processo decisionale ML e imposto un massimo di quattro parametri di ordine quando si calcola TAλc,λr. Per esempio, eccellenti prestazioni predittive (>90%) sono state ottenute utilizzando le macchine di gradient-boosting basate sull’ensemble (35, 36). Tuttavia, l’interpretazione del modello è problematica poiché viene utilizzato un ensemble di 100-150 alberi decisionali profondi (aggiunti in sequenza). Anche se la performance è migliorata, la possibilità di overfitting con correlazioni accidentali aumenta. Ci siamo quindi limitati ai modelli decisionali a singolo albero basati su alberi decisionali di classificazione e regressione (CART) (20). La restrizione a quattro parametri di ordine per la funzione TAλc,λr è basata su ragioni simili. Aggiungendo più parametri si ottengono matrici più rade che rappresentano le distribuzioni reattive/non reattive e, di conseguenza, l’integrazione numerica per calcolare la sovrapposizione tra queste distribuzioni diventa molto sensibile alla dimensione dei bin e potrebbe sottostimare la sovrapposizione a causa di bin vuoti per statistiche insufficienti.
Abbiamo considerato 138 variabili collettive che consistono nelle distanze ossigeno-ossigeno; le distanze ossigeno-idrogeno per le molecole d’acqua inizialmente legate; tutti gli angoli formati da Oλ e dai suoi quattro vicini di ossigeno più vicini; e i parametri di ordine Steinhardt degli ordini 3, 4 e 6 (32) (vedi Materiali e metodi per maggiori dettagli). Inoltre, sono stati aggiunti i parametri di ordine già considerati. La Fig. 5A mostra l’albero decisionale risultante. Notevolmente, di tutti i parametri di input, il parametro w4 è sia in cima all’albero decisionale che la variabile più importante come misurato dalla riduzione dell’errore di classificazione attribuito ad ogni variabile ad ogni divisione dell’albero decisionale (20) (Appendice SI, Fig. S9). Anche l’ordinamento tetraedrico e il numero di legami a idrogeno accettati appaiono nell’albero decisionale. Per descrivere il primo effetto, l’approccio ML ha dato la priorità al parametro d’ordine q4 di Steinhardt sopra il parametro q simile usato precedentemente da noi. Alcune distanze che appaiono anche nell’albero di decisione come d25, la distanza tra Oλ e il suo 25° ossigeno più vicino, sono molto probabilmente dovute a correlazioni accidentali causate dalla dimensione limitata del dataset. Questo è verificato controllando l’importanza di questa variabile: d25 non appare tra le 20 variabili più importanti (Appendice SI, Fig. S9), e, infatti, altre variabili simili (ad esempio, d24) sono classificate più in alto, anche se con bassa importanza. Un parametro più importante e intuitivamente valido che viene suggerito dall’approccio ML è λ2, la distanza OH tra l’ossigeno più vicino a Oλ e il suo idrogeno con il più grande legame intramolecolare. Il ricalcolo della capacità predittiva utilizzando i parametri dell’albero ML (Fig. 5B) non ha dato prestazioni superiori alla combinazione w4, q, na, e qcos, ma dovrebbe essere concepito come altrettanto buono, considerando le incertezze statistiche.
Risultati dell’analisi di apprendimento automatico. (A) Albero di classificazione e regressione per prevedere l’esito delle traiettorie iniziate. Qui, abbiamo considerato diverse variabili collettive aggiuntive (descrizione in Materiali e Metodi), ma solo un piccolo sottoinsieme è alla fine necessario per costruire l’albero: w4, q4 , λ2 (la lunghezza del legame idrogeno allungato nella molecola d’acqua più vicina alla specie Oλ), di (la distanza da Oλ all’iesimo ossigeno più vicino), e d¯i (la distanza media considerando gli i ossigeni più vicini). La notazione dei nodi è spiegata con il nodo autonomo in alto a sinistra. Questo albero predice che le traiettorie siano reattive, cioè che raggiungano un λ≥2, o non reattive in base alle variabili collettive ottenute al frame nelle traiettorie quando λ è prima ≥1,15. I nodi che predicono traiettorie reattive sono colorati in blu (classe 1) mentre i nodi che predicono traiettorie non reattive sono colorati in verde (classe 0). Si noti che le percentuali nella parte inferiore dei quadrati non riflettono le frazioni fisicamente corrette poiché gli insiemi di percorsi non sono stati ripesati usando i loro pesi statistici. Le regole sono rappresentazioni testuali dell’attraversamento dell’albero; per esempio, la regola 5 (che predice traiettorie reattive) può essere espressa come w4≥7.6 e λ2≥1.1. Queste regole danno diverse condizioni di inizio, e sono elencate nell’Appendice SI, Tabella S1, per la fila inferiore di nodi. (B) Il potere predittivo e la probabilità di attraversamento in funzione di λr per λc=1.16 Å e diverse combinazioni di variabili collettive. Qui confrontiamo il potere predittivo usando le variabili collettive che abbiamo identificato con le variabili segnate come importanti dall’analisi di machine-learning. (C) Distribuzioni reattive (rλc,λr(ξ)) e non reattive (uλc,λr(ξ)) per ξ={λ2,d¯2} e λc=1.16 Å e λr=2.0 Å. Per scopi di visualizzazione, le distribuzioni rappresentate sono normalizzate. Gli inserti in alto e a destra mostrano le proiezioni unidimensionali delle distribuzioni.