Significato
Le reazioni acido-base sono tra i processi chimici più importanti. Eppure ci manca un modo semplice di descrivere questa classe di reazioni in funzione delle coordinate atomiche. Infatti, una volta disciolti in acqua, H+ e il suo anione coniugato OH- hanno una struttura altamente fluttuante e difficile da definire. Qui risolviamo questo problema prendendo il punto di vista di descrivere le reazioni acido-base come un equilibrio tra il soluto e tutto il solvente. Questo permette di identificare descrittori generalmente applicabili. Di conseguenza, è ora possibile eseguire una simulazione quantitativa a campionamento avanzato della reazione acido-base in acqua e in altri ambienti come le cavità della zeolite o le superfici.
Abstract
Le reazioni acido-base sono onnipresenti in natura. La comprensione dei loro meccanismi è cruciale in molti campi, dalla biochimica alla catalisi industriale. Sfortunatamente, gli esperimenti danno solo informazioni limitate senza una grande comprensione del comportamento molecolare. Le simulazioni atomistiche potrebbero integrare gli esperimenti e gettare una luce preziosa sui meccanismi microscopici. Le grandi barriere di energia libera connesse alla dissociazione dei protoni, tuttavia, rendono obbligatorio l’uso di metodi di campionamento avanzati. Qui eseguiamo una simulazione di dinamica molecolare ab initio (MD) e miglioriamo il campionamento con l’aiuto della metadinamica. Ciò è stato reso possibile dall’introduzione di descrittori o variabili collettive (CVs) che si basano su una visione concettualmente diversa degli equilibri acido-base. Testiamo con successo il nostro approccio su tre diverse soluzioni acquose di acido acetico, ammoniaca e bicarbonato. Queste sono rappresentative del comportamento acido, basico e anfotero.
- acido-base
- metadinamica
- variabili collettive
- campionamento potenziato
Le reazioni acido-base hanno un ruolo chiave in molti rami della chimica. Le reazioni di complessazione inorganica, il ripiegamento delle proteine, i processi enzimatici, la polimerizzazione, le reazioni catalitiche e molte altre trasformazioni in diversi settori sono sensibili ai cambiamenti del pH. Comprendere il ruolo del pH in queste reazioni implica avere il controllo sulla loro reattività e cinetica.
L’importanza cruciale del pH ha stimolato la raccolta di una grande quantità di dati sugli equilibri acido-base. Questi sono tipicamente misurati in fasi gassose e condensate, usando tecniche spettroscopiche e potenziometriche. Tuttavia, ci sono limiti pratici alla precisione di questi metodi, specialmente nelle fasi condensate (1). Inoltre è molto difficile estrarre dai dati sperimentali un quadro microscopico dei processi coinvolti. Non è quindi sorprendente che l’equilibrio acido-base sia stato oggetto di un’intensa attività teorica (1⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓-12).
L’acidità di una specie chimica in acqua può essere espressa in termini di pKa, il logaritmo negativo della costante di dissociazione dell’acido. Ci sono due modi di calcolare questi valori, uno statico e l’altro dinamico.
L’approccio più standard è quello statico in cui le energie libere in fase di soluzione, e di conseguenza le pKa, sono ottenute chiudendo un ciclo Born-Haber composto da energie libere di fase gassosa e di solvatazione (1, 3⇓⇓⇓-7). Sebbene sia estremamente efficace in molti casi, l’approccio statico ha alcune limitazioni. Un modello di solvatazione deve essere scelto e i modelli di solvente continuo hanno una precisione limitata. Questo è particolarmente vero in sistemi come le zeoliti o le proteine caratterizzate da cavità irregolari in cui una descrizione implicita del solvente è impegnativa. Ovviamente da un tale approccio non si possono ottenere informazioni dinamiche. Inoltre, ci possono essere reazioni competitive che non possono essere prese in considerazione se non esplicitamente incluse nel modello.
In linea di principio queste limitazioni potrebbero essere eliminate in un approccio più dinamico basato su simulazioni di dinamica molecolare (MD) in cui le molecole del solvente sono trattate esplicitamente. Se si avesse tempo illimitato al computer, tali simulazioni esplorerebbero tutti i percorsi possibili e assegnerebbero il peso statistico relativo ai diversi stati. Sfortunatamente la presenza di colli di bottiglia cinetici vanifica questa possibilità intrappolando il sistema in stati metastabili, poiché i diversi stati di protonazione sono separati da grandi barriere. Inoltre nelle reazioni acido-base i legami chimici vengono rotti e formati. Questo richiede l’uso di MD ab initio in cui le forze interatomiche sono calcolate al volo dalle teorie di struttura elettronica. Questo rende il calcolo più costoso e riduce ulteriormente la scala temporale che può essere esplorata.
Per superare questa difficoltà, l’uso di metodi di campionamento migliorati (13) che accelerano l’esplorazione dello spazio configurazionale diventa obbligatorio. Una classe molto popolare di metodi di campionamento migliorati si basa sull’identificazione dei gradi di libertà che sono coinvolti nella reazione lenta di interesse. Questi gradi di libertà sono di solito indicati come variabili collettive (CV) e sono espressi come funzioni esplicite delle coordinate atomiche R. Il campionamento è quindi migliorato aggiungendo un bias che è una funzione dei CV scelti (14⇓-16). Inoltre, progettare un set adeguato di buoni CV ha anche un significato più profondo. Le CV di successo catturano in modo condensato la fisica del problema, identificano i suoi gradi di libertà lenti e conducono a un’utile descrizione modellistica del processo.
Nelle reazioni chimiche standard, questo è relativamente semplice poiché strutture ben definite possono essere assegnate a reagenti e prodotti (17⇓-19). Questo non è il caso delle reazioni acido-base in cui un protone viene aggiunto o sottratto al soluto. Una volta che questo processo ha avuto luogo, gli ioni d’acqua (H+ o/e OH-) sono solvatati e la loro struttura diventa sfuggente. Infatti, gli ioni d’acqua possono diffondere rapidamente nel mezzo attraverso un meccanismo di Grotthuss (20). Sono diventati altamente fluttuanti e l’identità degli atomi che partecipano alla loro struttura cambia continuamente. La natura di queste specie è quindi difficile da catturare in una funzione analitica esplicita di R. Tuttavia, data la rilevanza delle reazioni acido-base, sono stati fatti molti tentativi per definire queste entità (8⇓⇓⇓-12). Purtroppo queste CV hanno una natura ad hoc e, pur avendo successo in questo o quel caso, non possono essere applicate in generale.
Per costruire CV generali e utili facciamo due passi concettuali. Uno è quello di considerare il processo acido-base come una reazione che coinvolge solo alcune società, cioè l’intero solvente e i residui che reagiscono nella molecola solvata. Per esempio, quando c’è solo un tipo di residuo dissociante, pensiamo all’equilibrio acido-base come a una reazione del tipoA+H2NON⇌Bq0+H2N+q1ONq1,dove N è il numero di molecole d’acqua, A e B sono una generica molecola acido-base in soluzione e la sua specie coniugata, rispettivamente, q0 e q1 sono interi che possono assumere valori +1 e -1 secondo il comportamento acido-base della specie, e q1+q0=0.
Questo implica che non consideriamo il solvente come un insieme di molecole che competono per reagire con le specie acido-base. Piuttosto consideriamo il solvente nella sua interezza come uno dei due addotti. Assumere questo punto di vista è particolarmente rilevante in solventi polari come l’acqua che sono caratterizzati da reti altamente strutturate. In questo caso la presenza di un eccesso o di una carenza di protoni cambia localmente la struttura della rete e questa distorsione si propaga lungo l’intera rete.
Da quando Wicke ed Eigen (21) e Zundel e Metzger (22) hanno iniziato a lottare su quante molecole includere nella definizione della perturbazione (23⇓-25). Data l’assenza di parametri fisici in grado di dare una risposta chiara e inequivocabile a questa domanda, l’idea di considerare il solvente nel suo insieme aggira questo problema. Così, il solvente non è solo un mezzo con un ruolo passivo, ma viene considerato come un insieme di molecole che contribuiscono collettivamente alla formazione della coppia acido-base coniugata. Questo punto di vista è molto più vicino a quello originale proposto da Brønsted e Lowry in cui la reazione può essere vista come un semplice scambio di un catione idrogeno tra una coppia acido-base.
Perché la reazione abbia luogo il centro della perturbazione deve allontanarsi dal soluto. Quindi, il secondo passo importante è monitorare il centro della perturbazione. A causa dei meccanismi tipo Grotthuss, la perturbazione si muove lungo la rete. Questo può portare a diverse definizioni del centro del difetto. Tuttavia, se tasselliamo l’intero spazio usando i poliedri di Voronoi centrati sugli atomi di ossigeno dell’acqua, possiamo assegnare inequivocabilmente ogni atomo di idrogeno a uno e uno solo di questi poliedri. Il sito il cui poliedro di Voronoi contiene un numero anomalo di protoni viene preso come centro della perturbazione (Fig. 1).
Due esempi di partizionamento dello spazio. (A sinistra) Mostriamo un approccio convettivo in cui la distanza dall’atomo di ossigeno è usata per definire i suoi dintorni. Si vedono chiaramente sovrapposizioni artificiali. (Destra) La tassellatura di Voronoi non soffre di questi difetti.
Questo punto di vista dà al metodo una natura molto generale, rendendolo applicabile ad ogni sistema acido-base, senza la necessità di fissare in anticipo le coppie reagenti. Così, è possibile esplorare tutti gli stati di protonazione rilevanti anche in sistemi composti da più di una coppia acido-base.
Questo approccio generale permette di definire le CV senza dover imporre strutture specifiche o selezionare l’identità degli atomi coinvolti. Testiamo il nostro metodo eseguendo simulazioni di metadinamica in un caso di acido debole (acido acetico), in una base debole (ammoniaca), e in una specie anfotera (bicarbonato) scelti come benchmark a causa della loro forza comparabile, ma del loro diverso comportamento acido-base.
Metodi
Come discusso sopra introduciamo due CV, uno relativo allo stato di protonazione e un altro che localizza i difetti di carica e misura la loro distanza relativa. Entrambe queste CV hanno bisogno di una definizione robusta per assegnare gli atomi di idrogeno al rispettivo sito acido-base. Per ottenere questo risultato partizioniamo l’intero spazio in poliedri di Voronoi centrati sui siti acido-base i situati in Ri. I siti includono tutti gli atomi in grado di rompere e formare legami con un protone acido. La partizione standard dello spazio di Voronoi è descritta da un insieme di funzioni indice wi(r) centrate sui diversi Ris tali che wi(r)=1 se l’atomo i è il più vicino a r e wi(r) = 0 altrimenti. Per il loro utilizzo nei metodi di campionamento migliorati le CV devono essere differenziabili. A questo scopo introduciamo una versione liscia delle funzioni indice, wis(r). Queste sono definite usando softmax functionswis(r)=e-λ|Ri-r|∑me-λ|Rm-r|, dove i e m corrono su tutti i siti acido-base e λ controlla la ripidità con cui le curve decadono a 0, cioè la selettività della funzione. Con una scelta appropriata di λ questa definizione raggiunge il risultato desiderato come mostrato in Fig. 2. In questo modo, un atomo di idrogeno in una posizione Rj è assegnato al poliedro centrato sul sito i con il peso wi(Rj). Quindi, il numero totale di atomi di idrogeno assegnati all’iesimo sito acido-base èWi=∑j∈Hwis(Rj),dove la somma su j corre su tutti gli atomi di idrogeno.
Tessitura liscia di uno spazio 2D con celle centrate sui tre atomi di ossigeno della molecola d’acqua. Le regioni piatte blu rappresentano la porzione di spazio in cui la funzione assume un valore di 1 e quelle gialle rappresentano i confini tra le celle. Questa superficie è stata ottenuta con un valore di λ=4.
Si può associare ad ogni sito acido-base un valore di riferimento Wi0 che conta il numero di atomi di idrogeno legati allo stato neutro. La differenza tra il valore istantaneo degli atomi di idrogeno e quello di riferimento èδi=Wi-Wi0.Quando è diverso da zero, δi segnalerà se il sito iesimo ha guadagnato o perso un protone. Nel caso degli atomi di ossigeno dell’acqua, uno ione idronio ha un δi=+1 mentre uno ione idrossido ha δi=-1.
Raggruppiamo poi i siti acido-base in specie. Per esempio, nel caso del più semplice aminoacido glicina in soluzione acquosa il numero di specie Ns sarà uguale a 3. Tutti gli atomi di ossigeno dell’acqua appartengono a una specie; poi si contano in un’altra specie i due atomi di ossigeno carbossilico, e infine si considera come terza specie l’atomo di azoto del gruppo amminico.
Nello spirito di questo lavoro contiamo l’eccesso o il difetto totale di protoni associati a ciascuna specie:qk=∑i∈kδi.Ciò implica che non ci interessa l’identità specifica del sito reagito, ma se la kesima specie nel suo insieme è aumentata (qk=+1), è diminuita (qk=-1), o non ha cambiato il suo numero di protoni. Se consideriamo un soluto con una sola parte reattiva, allora ogni possibile stato del sistema può essere descritto da uno dei tre vettori 2D (0,0), (-1,1), o (1,-1).
Nel caso generale ogni stato di protonazione può essere descritto da un vettore q→=(q0,q1,…qNs-1) con dimensione uguale al numero di siti reattivi ineguali, Ns. Una spiegazione più esaustiva è fornita nell’appendice SI.
Per l’uso nel campionamento avanzato questi vettori devono essere espressi come una funzione scalare f=f(q→) tale che, per ogni q→ fisicamente rilevante, f raggiunga valori in grado di distinguere i diversi stati complessivi di protonazione. Ci sono infiniti modi di costruire uno scalare da un vettore. Forse la scelta più semplice è scrivere f(q→)=X→⋅q→ e, per distinguere i diversi stati di protonazione, scegliere X→=(20,21,22,…2Ns-1).
Questo porta alla seguente definizione per la CV che viene usata per descrivere lo stato di protonazione del sistema,sp=∑k=0Ns-12k⋅qk, dove k sono gli indici usati per etichettare i rispettivi gruppi di siti reattivi. Nell’appendice SI un esempio viene elaborato in dettaglio. Naturalmente il CV è reso continuo dall’uso di Wi nel calcolo del δi necessario per valutare qk nell’Eq. 5.
Il secondo CV è una sommatoria delle distanze tra tutti i siti acido-base moltiplicata per la loro carica parziale δi,sd=∑i,m>i-rim⋅δi⋅δm,dove gli indici i e m corrono su tutti i siti acido-base appartenenti a diversi gruppi k, e rim è la distanza tra i due atomi. In questo modo, solo la coppia acido-base che ha scambiato un protone dà un contributo diverso da zero. L’Eq. 7 è valida solo quando è presente una sola coppia acido-base coniugata. Tuttavia, a causa dell’azione del bias, può accadere occasionalmente che si formino più coppie acido-base. Per evitare il campionamento di questi eventi molto improbabili applichiamo un vincolo sul numero di coppie. Ulteriori dettagli sono forniti nell’appendice SI.
Risultati
Abbiamo applicato il nostro metodo a tre soluzioni acquose di acido acetico, ammoniaca e bicarbonato come rappresentazioni di un acido debole, una base debole e un composto anfotero, rispettivamente. I setup di tutte e tre le simulazioni sono identici tranne che per l’identità delle molecole solvatate. Questo assicura che il risultato riflette la chimica diversa di questi tre sistemi e che non vi è alcuna distorsione a causa della condizione iniziale.
Ogni simulazione dei sistemi è stata eseguita con Born-Oppenheimer simulazioni MD combinato con ben temperato metadinamica (14, 26) utilizzando il pacchetto CP2K (27) patchato con PLUMED 2 (28) e fortemente vincolato e adeguatamente normato funzionale (29) per l’energia xc, Exc. Vedi Appendice SI per i dettagli.
In Fig. 3 tracciamo le superfici di energia libera (FESs) in funzione di sp e sd. Queste FES riproducono vividamente il comportamento atteso. Hanno tutte un minimo a sp=0 che corrisponde allo stato in cui nessuna carica è presente nel solvente. Nella FES dell’acido acetico (Fig. 3A) un secondo minimo vicino a sp=-1 riflette il suo comportamento acido. Al contrario, l’ammoniaca FES (Fig. 3B) mostra un secondo minimo vicino a sp=1. Le forme delle FES dell’ammoniaca e dell’acido acetico sono approssimativamente correlate da una simmetria speculare che riflette il loro comportamento contrastante. Allo stesso modo la FES simmetrica del bicarbonato (Fig. 3C) rispecchia il suo carattere anfotero.
(A-C) Superfici di energia libera lungo sp e sd di acido acetico (A), ammoniaca (B), e bicarbonato (C) in soluzione acquosa. Le barre colorate indicano l’energia libera espressa in unità kJ⋅mol-1 . Il CV sd è espresso in angstroms.
Man mano che la coppia coniugata si forma sd inizia ad assumere valori positivi corrispondenti alla separazione e alla diffusione della coppia coniugata. Rispetto allo stato non associato in cui è permesso solo sd=0, gli stati in cui è presente una coppia coniugata mostrano una forma allungata dei bacini lungo questa variabile. Questo è causato dal comportamento diffusivo degli ioni idronio e idrossido in soluzione che rende accessibile un range continuo di distanze. Inoltre, lungo questo CV possiamo osservare una barriera intorno a 1,5 corrispondente alla rottura del legame covalente tra l’atomo di idrogeno e il sito acido-base.
Conclusioni
L’applicabilità generale di questo metodo a sistemi di diversa natura è un passo importante fatto nella loro comprensione e descrizione. Lo schema può essere esteso per includere effetti nucleari quantistici con l’uso della dinamica molecolare integrale del percorso (30). Questo sarebbe di importanza quantitativa poiché per esempio i valori di pKa sono influenzati dalla deuterizzazione. Inoltre, l’assenza di presupposti o imposizioni sui candidati reattivi o sui percorsi di reazione permette di estendere questo metodo a sistemi di complessità crescente che non possono essere affrontati con i metodi tradizionali. Esempi di domande a cui ora si può rispondere sono gli equilibri tautomerici nei processi biochimici e il comportamento acido nelle zeoliti e sulla superficie degli ossidi esposti all’acqua.
Riconoscimenti
I calcoli sono stati eseguiti sul cluster ETH Euler e sul cluster Mönch presso il Centro nazionale di supercalcolo svizzero. Questa ricerca è stata sostenuta dalla sovvenzione dell’Unione europea ERC-2014-AdG-670227/VARMET.
Pubblicato sotto licenza PNAS.