Znaczenie
Reakcje kwasowo-zasadowe należą do najważniejszych procesów chemicznych. Brakuje nam jednak prostego sposobu na opisanie tej klasy reakcji w funkcji współrzędnych atomowych. W rzeczywistości, po rozpuszczeniu w wodzie, H+ i jego sprzężony anion OH- mają bardzo fluksyjną strukturę, trudną do określenia. Tutaj rozwiązujemy ten problem przyjmując punkt widzenia opisujący reakcje kwasowo-zasadowe jako równowagę pomiędzy solutem a całym rozpuszczalnikiem. Pozwala to na zidentyfikowanie ogólnie stosowanych deskryptorów. W konsekwencji możliwe jest obecnie przeprowadzenie ilościowej symulacji reakcji kwasowo-zasadowych w wodzie oraz w innych środowiskach, takich jak wnęki zeolitów lub na powierzchniach.
Abstract
Reakcje kwasowo-zasadowe są wszechobecne w przyrodzie. Zrozumienie ich mechanizmów jest kluczowe w wielu dziedzinach, od biochemii po katalizę przemysłową. Niestety, eksperymenty dostarczają jedynie ograniczonych informacji, bez większego wglądu w zachowanie cząsteczek. Symulacje atomistyczne mogą uzupełnić eksperymenty i rzucić cenne światło na mechanizmy mikroskopowe. Jednak duże bariery swobodnej energii związane z dysocjacją protonów sprawiają, że konieczne jest zastosowanie metod rozszerzonego próbkowania. Tutaj przeprowadzamy symulację ab initio dynamiki molekularnej (MD) i zwiększamy próbkowanie za pomocą metadynamiki. Stało się to możliwe dzięki wprowadzeniu deskryptorów lub zmiennych kolektywnych (CV), które bazują na koncepcyjnie odmiennym spojrzeniu na równowagę kwasowo-zasadową. Nasze podejście testujemy z powodzeniem na trzech różnych wodnych roztworach kwasu octowego, amoniaku i wodorowęglanu. Są one reprezentatywne dla kwasowych, zasadowych i amfoterycznych zachowań.
- kwasowo-zasadowe
- metadynamika
- zmienne kolektywne
- wzmocnione próbkowanie
Reakcje kwasowo-zasadowe odgrywają kluczową rolę w wielu gałęziach chemii. Reakcje kompleksowania nieorganicznego, fałdowania białek, procesy enzymatyczne, polimeryzacja, reakcje katalityczne i wiele innych przemian w różnych dziedzinach są wrażliwe na zmiany pH. Zrozumienie roli pH w tych reakcjach implikuje kontrolę nad ich reaktywnością i kinetyką.
Kluczowe znaczenie pH pobudziło zbieranie dużej ilości danych na temat równowagi kwasowo-zasadowej. Są one zazwyczaj mierzone w gazie i fazach skondensowanych, przy użyciu technik spektroskopowych i potencjometrycznych. Istnieją jednak praktyczne ograniczenia dokładności tych metod, zwłaszcza w fazach skondensowanych (1). Ponadto, bardzo trudno jest uzyskać z danych eksperymentalnych mikroskopowy obraz zachodzących procesów. Nie jest więc zaskakujące, że równowaga kwasowo-zasadowa była przedmiotem intensywnej działalności teoretycznej (1⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓-12).
Kwasowość gatunku chemicznego w wodzie może być wyrażona w kategoriach pKa, ujemny logarytm stałej dysocjacji kwasu. Istnieją dwa sposoby obliczania tych wartości, jeden statyczny, a drugi dynamiczny.
Najbardziej standardowym podejściem jest podejście statyczne, w którym energie swobodne fazy roztworu, a w konsekwencji pKa, uzyskuje się przez zamknięcie cyklu Borna-Habera złożonego z energii swobodnych fazy gazowej i solwatacji (1, 3⇓⇓⇓⇓-7). Podejście statyczne, choć bardzo skuteczne w wielu przypadkach, ma pewne ograniczenia. Konieczne jest wybranie modelu solwatacji, a modele rozpuszczalników continuum mają ograniczoną dokładność. Jest to szczególnie prawdziwe w układach takich jak zeolity lub białka charakteryzujące się nieregularnymi wgłębieniami, w których opis rozpuszczalnika jest trudny. Oczywiście w takim podejściu nie można uzyskać informacji dynamicznej. Ponadto, mogą zachodzić reakcje konkurencyjne, które nie mogą być brane pod uwagę, chyba że są wyraźnie uwzględnione w modelu.
W zasadzie te ograniczenia mogłyby być zniesione w bardziej dynamicznym podejściu opartym na symulacjach dynamiki molekularnej (MD), w których cząsteczki rozpuszczalnika są traktowane w sposób jawny. Gdybyśmy mieli nieograniczony czas komputerowy, takie symulacje zbadałyby wszystkie możliwe ścieżki i przypisałyby względną wagę statystyczną różnym stanom. Niestety, obecność kinetycznych „wąskich gardeł” niweczy tę możliwość poprzez uwięzienie układu w stanach metastabilnych, ponieważ różne stany protonacji oddzielone są od siebie dużymi barierami. Ponadto w reakcjach kwasowo-zasadowych wiązania chemiczne są zrywane i tworzone. Wymaga to zastosowania ab initio MD, w którym siły międzyatomowe są obliczane w locie na podstawie teorii struktury elektronicznej. To sprawia, że obliczenia są droższe i jeszcze bardziej ogranicza skalę czasową, którą można zbadać.
Aby przezwyciężyć tę trudność, obowiązkowe staje się zastosowanie metod wzmocnionego próbkowania (13), które przyspieszają eksplorację przestrzeni konfiguracyjnej. Bardzo popularna klasa metod zwiększonego próbkowania oparta jest na identyfikacji stopni swobody, które są zaangażowane w interesującą nas powolną reakcję. Te stopnie swobody są zwykle nazywane zmiennymi kolektywnymi (CV) i są wyrażone jako jawne funkcje współrzędnych atomowych R. Próbkowanie jest następnie wzmocnione przez dodanie stronniczości, która jest funkcją wybranych CV (14⇓-16). Ponadto, zaprojektowanie odpowiedniego zestawu dobrych CV ma również głębsze znaczenie. Udane CV ujmują w skondensowany sposób fizykę problemu, identyfikują jego wolne stopnie swobody i prowadzą do użytecznego modelowego opisu procesu.
W standardowych reakcjach chemicznych jest to stosunkowo proste, ponieważ dobrze zdefiniowane struktury mogą być przypisane do reaktantów i produktów (17⇓-19). Nie jest tak w przypadku reakcji kwasowo-zasadowych, w których proton jest dodawany do lub odejmowany od solutu. Po zajściu tego procesu, jony wody (H+ lub/i OH-) są rozpuszczane, a ich struktura staje się nieuchwytna. W rzeczywistości, jony wody mogą szybko dyfundować w medium poprzez mechanizm Grotthussa (20). Stają się one wysoce fluksyjne, a tożsamość atomów biorących udział w ich strukturze ulega ciągłym zmianom. Natura tych gatunków jest więc trudna do uchwycenia w jednoznacznej funkcji analitycznej R. Jednakże, biorąc pod uwagę znaczenie reakcji kwasowo-zasadowych, podjęto wiele prób zdefiniowania tych jednostek (8⇓⇓⇓⇓-12). Niestety te CV mają charakter ad hoc i, choć udane w tym czy innym przypadku, nie mogą być powszechnie stosowane.
Aby zbudować ogólne i użyteczne CV robimy dwa kroki koncepcyjne. Jednym z nich jest spojrzenie na proces kwasowo-zasadowy jako na reakcję, w którą zaangażowanych jest tylko kilka cząsteczek, mianowicie cały rozpuszczalnik i reagujące reszty w rozpuszczanej cząsteczce. Na przykład, gdy istnieje tylko jeden typ dysocjującej reszty, myślimy o równowadze kwasowo-zasadowej jako o reakcji typuA+H2NON⇌Bq0+H2N+q1ONq1,gdzie N jest liczbą cząsteczek wody, A i B są odpowiednio ogólną cząsteczką kwasowo-zasadową w roztworze i jej sprzężonym gatunkiem, q0 i q1 są liczbami całkowitymi, które mogą przyjmować wartości +1 i -1 zgodnie z zachowaniem kwasowo-zasadowym gatunku, a q1+q0=0.
Z tego wynika, że nie patrzymy na rozpuszczalnik jako na zbiór cząsteczek, które konkurują ze sobą w reakcji z gatunkami kwasowo-zasadowymi. Raczej rozważamy rozpuszczalnik w całości jako jeden z dwóch adduktów. Ten punkt widzenia jest szczególnie istotny w rozpuszczalnikach polarnych, takich jak woda, które charakteryzują się wysoce uporządkowanymi sieciami. W tym przypadku obecność nadmiaru lub niedoboru protonów zmienia lokalnie strukturę sieci i to zniekształcenie propaguje się wzdłuż całej sieci.
Od bardzo wczesnych dni Wicke i Eigen (21) oraz Zundel i Metzger (22), badacze zmagali się z tym, ile cząsteczek powinno być uwzględnionych w definicji perturbacji (23⇓-25). Ze względu na brak parametrów fizycznych, które mogłyby dać jasną i jednoznaczną odpowiedź na to pytanie, koncepcja rozpatrywania rozpuszczalnika jako całości omija ten problem. Rozpuszczalnik nie jest więc tylko medium pełniącym bierną rolę, ale jest postrzegany jako zespół cząsteczek, które wspólnie przyczyniają się do powstania sprzężonej pary kwas-zasada. Ten punkt widzenia jest znacznie bliższy oryginalnemu zaproponowanemu przez Brønsteda i Lowry’ego, w którym reakcja może być postrzegana jako prosta wymiana kationu wodorowego pomiędzy parą kwas-zasada.
Aby reakcja mogła zajść, centrum perturbacji musi oddalić się od solutu. Dlatego drugim ważnym krokiem jest monitorowanie centrum perturbacji. Ze względu na mechanizmy podobne do Grotthussa, perturbacja przemieszcza się wzdłuż sieci. Może to prowadzić do różnych definicji centrum defektu. Jeśli jednak teselujemy całą przestrzeń za pomocą wielościanów Voronoi skupionych na atomach tlenu wody, możemy jednoznacznie przypisać każdy atom wodoru do jednego i tylko jednego z tych wielościanów. Miejsce, którego wielościan Voronoi zawiera anomalną liczbę protonów, przyjmujemy za centrum perturbacji (Rys. 1).
Dwa przykłady podziału przestrzeni. (Po lewej) Pokazujemy podejście konwekcyjne, w którym odległość od atomu tlenu jest używana do definiowania jego otoczenia. Wyraźnie widać sztuczne superpozycje. (Po prawej) Teselacja Voronoi nie ma tych wad.
Ten punkt widzenia nadaje metodzie bardzo ogólny charakter, czyniąc ją możliwą do zastosowania w każdym układzie kwasowo-zasadowym, bez potrzeby wcześniejszego ustalania reagujących par. Dzięki temu możliwe jest badanie wszystkich istotnych stanów protonacyjnych nawet w układach złożonych z więcej niż jednej pary kwas-zasada.
Takie ogólne podejście pozwala na definiowanie CV bez konieczności narzucania specyficznych struktur lub wybierania tożsamości atomów biorących udział w reakcji. Testujemy naszą metodę wykonując symulacje metadynamiki w przypadku słabego kwasu (kwas octowy), słabej zasady (amoniak) i amfoterycznego gatunku (wodorowęglan) wybranych jako wzorce ze względu na ich porównywalną siłę, ale różne zachowanie kwasowo-zasadowe.
Metody
Jak omówiono powyżej wprowadzamy dwa CV, jedno związane ze stanem protonacji i drugie, które lokalizuje defekty ładunkowe i mierzy ich względną odległość. Oba te CV wymagają solidnej definicji, aby przypisać atomy wodoru do odpowiedniego miejsca kwasowo-zasadowego. Aby osiągnąć ten rezultat, dzielimy całą przestrzeń na wielościany Voronoi skupione wokół miejsc kwasowo-zasadowych i zlokalizowanych w Ri. Miejsca te obejmują wszystkie atomy zdolne do zrywania i tworzenia wiązań z protonem kwasu. Standardowy podział przestrzeni Voronoi opisany jest przez zbiór funkcji indeksowych wi(r) skupionych na różnych Ris, takich, że wi(r)=1 jeśli i-ty atom jest najbliżej r, a wi(r) = 0 w przeciwnym przypadku. W celu wykorzystania ich w metodach rozszerzonego próbkowania CV muszą być różniczkowalne. W tym celu wprowadzamy gładką wersję funkcji indeksowych, wis(r). Są one zdefiniowane przy użyciu funkcji softmaxswis(r)=e-λ|Ri-r|∑me-λ|Rm-r|, gdzie i oraz m przebiegają przez wszystkie miejsca kwasowo-zasadowe, a λ kontroluje stromość, z jaką krzywe opadają do 0, czyli selektywność funkcji. Przy odpowiednim doborze λ definicja ta pozwala uzyskać pożądany rezultat przedstawiony na Rys. 2. W ten sposób atom wodoru w pozycji Rj jest przypisywany do wielościanu skupionego na miejscu i z wagą wi(Rj). Wówczas całkowita liczba atomów wodoru przypisanych do i-tego miejsca kwasowo-zasadowego wynosiWi=∑j∈Hwis(Rj),gdzie sumowanie na j przebiega po wszystkich atomach wodoru.
Gładka teselacja przestrzeni 2D z komórkami skupionymi na trzech atomach tlenu cząsteczki wody. Płaskie niebieskie obszary reprezentują część przestrzeni, w której funkcja przyjmuje wartość 1, a żółte reprezentują granice między komórkami. Powierzchnia ta została uzyskana przy wartości λ=4.
Do każdego miejsca kwasowo-zasadowego można przypisać wartość odniesienia Wi0, która liczy liczbę związanych atomów wodoru w stanie obojętnym. Różnica pomiędzy chwilową wartością atomów wodoru a wartością referencyjną wynosiδi=Wi-Wi0.Gdy jest różna od zera, δi sygnalizuje, czy i-te miejsce zyskało czy straciło proton. W przypadku atomów tlenu wody, jon hydroniowy ma δi=+1, natomiast jon wodorotlenkowy ma δi=-1.
Następnie grupujemy miejsca kwasowo-zasadowe w gatunki. Na przykład, w przypadku najprostszego aminokwasu glicyny w roztworze wodnym liczba gatunków Ns będzie równa 3. Do jednego gatunku należą wszystkie atomy tlenu wody, do drugiego gatunku zaliczamy dwa karboksylowe atomy tlenu, wreszcie za trzeci gatunek uważamy atom azotu grupy aminowej.
W duchu tej pracy liczymy całkowity nadmiar lub ubytek protonów związanych z każdym gatunkiem:qk=∑i∈kδi.Oznacza to, że nie interesuje nas konkretna tożsamość miejsca reakcji, lecz to, czy k-ty gatunek w całości zwiększył (qk=+1), zmniejszył (qk=-1), czy też nie zmienił swojej liczby protonów. Jeśli rozważamy solut z tylko jedną reaktywną cząsteczką, to każdy możliwy stan układu może być opisany przez jeden z trzech wektorów 2D (0,0), (-1,1) lub (1,-1).
W ogólnym przypadku każdy stan protonacji może być opisany przez wektor q→=(q0,q1,…qNs-1) o wymiarze równym liczbie nierównoważnych miejsc reaktywnych, Ns. Bardziej wyczerpujące wyjaśnienie znajduje się w Dodatku SI.
Do zastosowania w próbkowaniu rozszerzonym wektory te muszą być wyrażone jako funkcja skalarna f=f(q→) taka, że dla każdego fizycznie istotnego q→, f osiąga wartości zdolne do rozróżnienia różnych ogólnych stanów protonacji. Istnieje nieskończenie wiele sposobów na skonstruowanie skalara z wektora. Prawdopodobnie najprostszym wyborem jest zapisanie f(q→)=X→⋅q→ oraz, w celu rozróżnienia różnych stanów protonacji, wybranie X→=(20,21,22,…2Ns-1).
Prowadzi to do następującej definicji CV, która jest używana do opisu stanu protonacji układu,sp=∑k=0Ns-12k⋅qk,gdzie k są indeksami używanymi do oznaczenia odpowiednich grup miejsc reaktywnych. W Dodatku SI przykład jest szczegółowo opracowany. Oczywiście CV jest ciągłe poprzez użycie Wi do obliczenia δi potrzebnego do oszacowania qk w równaniu 5.
Drugie CV jest sumą odległości pomiędzy wszystkimi miejscami kwasowo-zasadowymi pomnożoną przez ich ładunek cząstkowy δi,sd=∑i,m>i-rim⋅δi⋅δm,gdzie indeksy i oraz m biegną przez wszystkie miejsca kwasowo-zasadowe należące do różnych grup k, a rim jest odległością pomiędzy dwoma atomami. W ten sposób tylko ta para kwas-zasada, która wymieniła proton, daje wkład różny od zera. Równanie 7 jest słuszne tylko wtedy, gdy obecna jest jedna pojedyncza sprzężona para kwasowo-zasadowa. Jednakże, ze względu na działanie nierówności, może się zdarzyć, że od czasu do czasu powstanie kilka par kwas-zasada. Aby uniknąć próbkowania tych bardzo mało prawdopodobnych zdarzeń, stosujemy ograniczenie na liczbę par. Dalsze szczegóły podano w Dodatku SI.
Wyniki
Naszą metodę zastosowaliśmy do trzech wodnych roztworów kwasu octowego, amoniaku i wodorowęglanu jako reprezentacji odpowiednio słabego kwasu, słabej zasady i związku amfoterycznego. Układy wszystkich trzech symulacji są identyczne, z wyjątkiem tożsamości rozpuszczonych cząsteczek. Zapewnia to, że wyniki odzwierciedlają różną chemię tych trzech układów i że nie ma uprzedzeń wynikających z warunków początkowych.
Każda symulacja układów została przeprowadzona za pomocą symulacji MD Borna-Oppenheimera w połączeniu z dobrze wytłumioną metadynamiką (14, 26) przy użyciu pakietu CP2K (27) z poprawką PLUMED 2 (28) i silnie ograniczoną i odpowiednio znormalizowaną funkcją (29) dla energii xc, Exc. Zobacz Dodatek SI dla szczegółów.
Na Rys. 3 wykreślamy powierzchnie swobodnej energii (FES) jako funkcję sp i sd. Te FES-y żywo odtwarzają oczekiwane zachowanie. Wszystkie mają minimum przy sp=0, które odpowiada stanowi, w którym nie ma ładunków w rozpuszczalniku. W przypadku kwasu octowego FES (Rys. 3A) drugie minimum w pobliżu sp=-1 odzwierciedla jego kwasowe zachowanie. Natomiast w przypadku amoniaku FES (Rys. 3B) drugie minimum znajduje się w pobliżu sp=1. Kształty FES amoniaku i kwasu octowego są w przybliżeniu powiązane przez lustrzaną symetrię odzwierciedlającą ich kontrastowe zachowanie. Podobnie symetryczny FES wodorowęglanu (Rys. 3C) odzwierciedla jego amfoteryczny charakter.
(A-C) Powierzchnie swobodnej energii wzdłuż sp i sd kwasu octowego (A), amoniaku (B) i wodorowęglanu (C) w roztworze wodnym. Kolorowe paski oznaczają energię swobodną wyrażoną w jednostkach kJ⋅mol-1. CV sd wyrażone jest w angstremach.
W miarę tworzenia się pary sprzężonej sd zaczyna przyjmować wartości dodatnie odpowiadające separacji i dyfuzji pary sprzężonej. W porównaniu ze stanem niezdysocjowanym, w którym dopuszczalne jest tylko sd=0, stany, w których obecna jest para sprzężona wykazują wydłużony kształt basenów wzdłuż tej zmiennej. Spowodowane jest to dyfuzyjnym zachowaniem jonów hydroniowych i wodorotlenkowych w roztworze, które udostępnia kontinuum odległości. Ponadto, wzdłuż tego CV możemy zaobserwować barierę około 1.5 odpowiadającą rozerwaniu wiązania kowalencyjnego pomiędzy atomem wodoru a miejscem kwasowo-zasadowym.
Wnioski
Ogólne zastosowanie tej metody do układów o różnej naturze jest ważnym krokiem poczynionym w ich zrozumieniu i opisie. Schemat ten może być rozszerzony o kwantowe efekty jądrowe z wykorzystaniem dynamiki molekularnej z całką ścieżkową (30). Miałoby to znaczenie ilościowe, ponieważ na przykład wartości pKa są zależne od deuteracji. Co więcej, brak założeń lub narzuceń dotyczących kandydatów na reaktantów lub ścieżek reakcji pozwala na rozszerzenie tej metody na układy o rosnącej złożoności, które nie mogą być badane tradycyjnymi metodami. Przykładami pytań, na które można teraz odpowiedzieć są równania tautomeryczne w procesach biochemicznych oraz zachowanie kwasów w zeolitach i na powierzchni tlenków wystawionych na działanie wody.
Podziękowania
Obliczenia przeprowadzono na klastrze ETH Euler oraz na klastrze Mönch w Swiss National Supercomputing Center. Badania te były wspierane przez Grant Unii Europejskiej ERC-2014-AdG-670227/VARMET.
Publikacja na licencji PNAS.
.