Wieloskalowe pobieranie próbek heterogenicznego granicy faz woda/katalizator metalowy przy użyciu teorii funkcjonału gęstości i dynamiki molekularnej pola siłowego

Summary

Celem przedstawionego tutaj protokołu jest generowanie i próbkowanie trajektorii konfiguracji cząsteczek ciekłej wody wokół gatunków katalitycznych na płaskiej powierzchni metalu przejściowego. Próbkowane konfiguracje mogą być używane jako struktury wyjściowe w metodach opartych na mechanice kwantowej.

Abstract

Znaczna liczba heterogenicznie katalizowanych procesów chemicznych zachodzi w warunkach ciekłych, ale symulacja funkcji katalizatora w takich warunkach jest trudna, gdy konieczne jest uwzględnienie cząsteczek rozpuszczalnika. Procesy zrywania i tworzenia wiązań modelowane w tych układach wymagają zastosowania metod chemii kwantowej. Ponieważ cząsteczki w fazie ciekłej znajdują się w ciągłym ruchu termicznym, symulacje muszą również obejmować pobieranie próbek konfiguracyjnych. Oznacza to, że dla każdego interesującego nas gatunku katalitycznego należy symulować wiele konfiguracji cząsteczek cieczy. Celem przedstawionego tutaj protokołu jest wygenerowanie i pobranie próbek trajektorii konfiguracji cząsteczek ciekłej wody wokół form katalitycznych na płaskich powierzchniach metali przejściowych w sposób, który równoważy dokładność chemiczną z kosztami obliczeniowymi. W szczególności symulacje dynamiki molekularnej pola siłowego (FFMD) są wykorzystywane do generowania konfiguracji cząsteczek cieczy, które można następnie wykorzystać w metodach opartych na mechanice kwantowej, takich jak teoria funkcjonału gęstości lub dynamika molekularna ab initio. Aby to zilustrować, w tym manuskrypcie protokół jest używany dla katalitycznych produktów pośrednich, które mogą być zaangażowane w szlak rozkładu glicerolu (_C3H₈O₃). Struktury generowane za pomocą FFMD są modelowane w DFT w celu oszacowania entalpii solwatacji związków katalitycznych i określenia, w jaki sposób cząsteczki_H2Ouczestniczą w rozkładach katalitycznych.

Introduction

Modelowanie zjawisk molekularnych związanych z heterogeniczną katalizą w warunkach ciekłych jest niezbędne do zrozumienia funkcji katalitycznej; jednak pozostaje to wyzwaniem, ponieważ wymaga doskonałej równowagi między dokładnością chemiczną a wydatkami obliczeniowymi. Ogólnie rzecz biorąc, ponieważ kataliza polega na zrywaniu i tworzeniu wiązań chemicznych, mechanika kwantowa musi być stosowana przynajmniej do pewnego stopnia; Jednak długie symulacje są wyzwaniem w mechanice kwantowej, ponieważ wymagają znacznych zasobów komputerowych. Ponieważ cząsteczki w fazie ciekłej znajdują się w ciągłym ruchu termicznym, symulacje muszą również obejmować próbkowanie konfiguracyjne, tj. muszą obejmować wiele układów przestrzennych cząsteczek cieczy, ponieważ każdy inny układ przestrzenny (tj. każda konfiguracja) ma inną energię. Oznacza to, że dla każdego interesującego nas gatunku katalitycznego należy symulować wiele konfiguracji cząsteczek cieczy. Te potrzeby – wykorzystanie mechaniki kwantowej i wykonywanie wielu obliczeń dla każdego rodzaju katalitycznego – mogą sprawić, że modelowanie w katalizie heterogenicznej w fazie ciekłej stanie się obliczeniowo niemożliwe. Celem opisanej w niniejszym dokumencie metody jest umożliwienie wykonalnych obliczeniowo symulacji zjawisk zachodzących w katalizie heterogenicznej w fazie ciekłej.

Jesteśmy szczególnie zainteresowani heterogenicznie katalizowanymi reakcjami, które są przeprowadzane pod ciekłą wodą. Cząsteczki wody mają znaczący wpływ na zjawiska katalityczne, takie jak oddziaływanie z gatunkami katalitycznymi (np. poprzez siły dyspersji i wiązania wodorowe)^{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20, 21,22,23}, uczestniczący w reakcjach katalitycznych1^,7^,8^{,9,15,21,22,24,25,26,27}, oraz wpływające na ścieżki reakcji i/lub szybkość katalityczną1^{,11,12,15,18,23,25,27, 28,29,30,31}. Modelowanie tych zjawisk zostało przeprowadzone przy użyciu QM i/lub ab initio dynamiki molekularnej (AIMD)^{1,2,6,7,14,22,25,27,28,32,33,34}, dynamika molekularna pola siłowego (FFMD)³⁵ oraz mechanika kwantowa/mechanika molekularna (QM/MM)¹⁰. W AIMD i FFMD atomy w układzie są poruszane zgodnie z równaniami ruchu Newtona zgodnie z działającymi na nie siłami. W AIMD energia i siły układu są obliczane za pomocą mechaniki kwantowej, podczas gdy w FFMD energia i siły układu są obliczane za pomocą pól siłowych, które są wyrażeniami algebraicznymi sparametryzowanymi na podstawie danych eksperymentalnych lub QM. W QM/MM część układu, w której następuje zerwanie i tworzenie wiązania, jest obliczana za pomocą QM, a pozostała część systemu jest obliczana za pomocą MM, który wykorzystuje pola siłowe. Ponieważ bezpośrednio wykorzystują QM, AIMD i QM/MM lepiej nadają się do wychwytywania zrywania i tworzenia wiązań, które zachodzą w katalizie heterogenicznej w fazie wodnej; jednak FFMD jest znacznie bardziej wykonalny obliczeniowo, a tym samym lepiej nadaje się do generowania konfiguracji ciekłych cząsteczek H₂O. Metoda przedstawiona w tym protokole równoważy dokładność chemiczną i koszty obliczeniowe poprzez zastosowanie kombinacji QM i FFMD.

W szczególności, ta metoda wykorzystuje symulacje FFMD do generowania konfiguracji cieczy_H2Oi QM do obliczania energii systemu. FFMD odbywa się przy użyciu LAMMPS. ³⁶ Pola siłowe użyte w FFMD w tej pracy wykorzystują potencjały Lennarda-Jonesa + Coulomba (LJ+C), gdzie parametry LJ zostały zaczerpnięte z modelu TIP3P/CHARMM³⁷ dla H₂O, uniwersalnego pola siłowego³⁸ (UFF) dla Pt, oraz pola siłowego OPLS-AA³⁹ dla form katalitycznych, a parametry Coulomba zostały zaczerpnięte z modelu TIP3P/CHARMM37 dla H₂O i pole siłowe OPLS-AA³⁹ dla gatunków katalitycznych. Parametry kulombowskie dla atomów Pt zostały ustawione na 0. Obliczenia QM są wykonywane przy użyciu kodu VASP ^40,41,42, który jest kodem teorii funkcjonału gęstości (DFT). Insercje cząsteczek wody są wykonywane za pomocą opracowanego przez nas kodu o nazwie Monte Carlo Plug-in for Quantum Methods (MCPliQ). Konwersje plików z VASP do LAMMPS w tym protokole są wykonywane za pomocą oprogramowania Visual Molecular Dynamics (VMD) ⁴³.

Protokół ma na celu generowanie konfiguracji cząsteczek ciekłej wody wokół katalitycznych form na płaskich powierzchniach metalu przejściowego przy niskim pokryciu. Pokrycie jest oznaczone jako θ i zdefiniowane jako liczba adsorbatów na powierzchniowy atom metalu (tj. liczba adsorbatów powierzchniowych znormalizowana przez liczbę atomów metalu w najwyższej warstwie metalowej płyty w modelu katalizatora). W tym manuskrypcie niskie pokrycie zdefiniowano jako θ ≤ 1/9 monowarstwy (ML), gdzie 1 ML oznacza jeden rodzaj katalityczny na atom metalu na powierzchni. Modele katalizatorów należy umieszczać w skrzynkach do symulacji okresowych. Pola symulacji nie muszą być sześcianami. Niniejszy manuskrypt demonstruje zastosowanie protokołu do generowania konfiguracji cieczy_H2O, które można wykorzystać do obliczenia wielkości będących przedmiotem zainteresowania w katalizie heterogenicznej w fazie wodnej.

Ten protokół wymaga, aby użytkownik miał dostęp do zainstalowanych i działających wersji oprogramowania VASP, MCPliQ, LAMMPS i VMD. Więcej informacji na temat VASP (https://www.vasp.at/), LAMMPS (https://Lammps.sandia.gov/) i VMD (https://www.ks.uiuc.edu/Research/vmd/) można znaleźć na ich stronach internetowych. Oprogramowanie MCPliQ jest udokumentowane na https://github.com/getman-research-group/JoVE_article, wraz ze wszystkimi plikami wejściowymi i skryptami Pythona wymienionymi w tym protokole. Protokół ten zakłada, że wymienione w nim pliki wykonywalne i skrypty będą uruchamiane na komputerze badawczym o wysokiej wydajności i są instalowane w katalogu, który znajduje się w zmiennej $PATH użytkownika. Jeśli plik wykonywalny lub skrypt zostanie umieszczony w lokalizacji, która nie znajduje się w $PATH użytkownika, należy dołączyć ścieżkę do pliku wykonywalnego, aby go wykonać. Pliki wykonywalne i skrypty są wykonywane w krokach 2.1.2, 2.2.1, 2.2.8, 3.1, 4.2, 5.2 i 6.1.2. Na przykład, aby wykonać kod MCPliQ w kroku 2.1.2 z katalogu, który nie znajduje się w $PATH użytkownika, użytkownik powinien wpisać $PATHTOMCPLIQ/mcpliq w interfejsie wiersza poleceń zamiast mcpliq, gdzie $PATHTOMCPLIQ jest lokalizacją, w której został zapisany plik wykonywalny mcpliq (np$PATHTOMCPLIQ. może to być ~/bin). Przed uruchomieniem tego protokołu wszystkie pliki wykonywalne i skrypty powinny mieć uprawnienia do wykonywania plików wykonywalnych (np. w Linuksie można to zrobić, wpisując chmod +x mcpliq w interfejsie wiersza poleceń z katalogu, w którym przechowywany jest plik wykonywalny mcpliq). Ponadto należy załadować wszystkie moduły wymagane przez którekolwiek z programów lub skryptów (zależności te będą specyficzne dla poszczególnych instalacji różnych programów i komputera, na którym zostaną uruchomione symulacje).

Protocol

1. Wygeneruj strukturę adsorbatu

1. Utwórz plik VASP POSCAR składający się z superkomórki z okresowymi warunkami brzegowymi, tak jak w przypadku przeprowadzania symulacji adsorbatów na powierzchniach metalowych w próżni. Superkomórka powinna zawierać wstępne odgadnięcie struktury adsorbatu i powierzchni metalu, a także przestrzeń próżniową nad adsorbatem w celu dodania cząsteczek_H2O. Szczegółowe informacje o tym, jak to zrobić, znajdują się w poprzedniej wersji work^35,44.
   UWAGA: Ważne jest, aby wysokość przestrzeni próżniowej znajdowała się co najmniej 12 A powyżej górnej części adsorbatu.
2. Rozluźnij konstrukcję i zminimalizuj jej energię za pomocą kodu VASP. Szczegółowe informacje o tym, jak to zrobić, znajdują się w poprzedniej wersji work^35,44. Spowoduje to utworzenie pliku o nazwie CONTCAR, który zostanie użyty w następnej sekcji.

2. Dodaj wyraźne cząsteczki_H2O

1. Dodaj N cząsteczek wody do przestrzeni próżniowej w CONTCAR utworzonej w kroku 1.2 przy użyciu kodu MCPliQ, gdzie N = ρV, ρ to gęstość wody, a V to objętość przestrzeni próżniowej nad adsorbatem.
   UWAGA: ρ należy przyjąć jako gęstość wody określoną przez model wody TIP3P/CHARMM w temperaturze symulacji. V zostanie dopracowany w następnym kroku.
   1. Określ następujące informacje w pliku MCPliQ master_input.txt: liczbę cząsteczek H₂O, które należy dodać (N), zmieniając pierwszy argument w wierszu 28, ścieżkę do pliku water.txt, zmieniając drugi argument w wierszu 28 oraz minimalną i maksymalną wysokość superkomórki, która może być zajęta przez cząsteczki wody, zmieniając minimalną współrzędną z w wierszu 11 i maksymalną współrzędną z w wierszu 12.
   2. Wykonaj kod MCPliQ, wpisując mcpliq z interfejsu wiersza poleceń, aby wstawić cząsteczki wody do pliku CONTCAR. Kod wyświetli jeden lub więcej plików z rozszerzeniem pliku . POSCAR.
      UWAGA: Jeśli zostanie wyprodukowany więcej niż jeden plik POSCAR, zostaną one nazwane POSCAR_n.POSCAR. Wybierz plik, w którym n jest największe.
2. Wygeneruj pliki wejściowe LAMMPS dla symulacji NPT i zrównoważ objętość komórki za pomocą FFMD w zespole NPT w LAMMPS.
   1. Wykonaj skrypt lmps_bond_angle.py na . POSCAR wygenerowany w kroku 2.1.2, wpisując lmps_bond_angle.py $filename. POSCAR w interfejsie wiersza poleceń, gdzie $filename jest nazwą pliku . Plik POSCAR wygenerowany w kroku 2.1.2. Ten skrypt tworzy plik o nazwie $filename. POSCAR.bond_angle_info.txt który zawiera listę wiązań i kątów, które zostaną użyte w pliku danych LAMMPS.
   2. Otwórz VMD i wybierz File > New Molecule w oknie głównym, aby otworzyć okno Molecule File Browser. Wybierz VASP_POSCAR z menu rozwijanego Określ typ pliku. Kliknij przycisk Przeglądaj i przejdź do $filename. Plik POSCAR. Kliknij przycisk Załaduj, aby otworzyć $filename. Plik POSCAR.
   3. Otwórz konsolę Tk w VMD, wybierając Rozszerzenia > konsoli Tk w oknie głównym VMD.
   4. Wykonaj następujące polecenie w konsoli Tk: topo writelammpsdata $WDPATH/data.myadsorbate full, gdzie $WDPATH jest katalogiem na komputerze, w którym VMD zapisze plik danych LAMMPS, a data.myadsorbate to nazwa pliku danych LAMMPS.
   5. Usuń sekcję Wiązania i kąty u dołu pliku data.myadsorbate. Następnie dołącz listy wiązań i kątów do pliku $filename. POSCAR.bond_angle_info.text do data.myadsorbate.
      UWAGA: Indeksy dla typu wiązania O-H i typu kąta OH dla cząsteczek wody w $filename. POSCAR.bond_angle_info.txt plik są ustawione na wartość 1. W związku z tym typy wiązań i kątów dla adsorbatów powinny zaczynać się od 2.
   6. Edytuj plik data.myadsorbate, dodając parametry Lennarda-Jonesa do sekcji Współczynniki par i parametry Coulomba do sekcji Atomy. Należy dodać parametry Lennarda-Jonesa i Coulomba dla cząsteczek_H2O, atomów adsorbatu i atomów powierzchni metalu.
      UWAGA: Parametry Lennarda-Jonesa dla atomów Pt, cząsteczek wody i atomów adsorbatu w tym protokole są uzyskiwane odpowiednio z pól siłowych UFF³⁸, TIP3P/CHARMM³⁷ i OPLS-AA³⁹. Parametry kulombowskie dla cząsteczek wody i atomów adsorbatów uzyskuje się odpowiednio z pól siłowych TIP3P/CHARMM³⁷ i OPLS-AA³⁹. Parametry kulombowskie dla atomów Pt są ustawione na 0 w tym protokole. Alternatywnie, obliczone ładunki cząstkowe można wykorzystać do parametrów kulombowskich dla atomów adsorbatu i atomów Pt.
   7. Skopiuj plik wejściowy LAMMPS input.equil do katalogu $WDPATH. Edytuj zmienną grupy w wierszu 34, aby wskazać indeksy typu atomu dla atomów wody, tlenu i wodoru wody, oraz zmienną grupy w wierszu 35, aby wskazać indeksy typu atomu dla atomów Pt i adsorbatu.
   8. Uruchom oprogramowanie LAMMPS, wpisując mpiexec -np XX lmp_mpi < input.equil w interfejsie wiersza poleceń, gdzie XX to liczba rdzeni procesora do użycia, a lmp_mpi to nazwa pliku wykonywalnego LAMMPS. W ten sposób uruchomiona zostanie minimalizacja energii w celu udoskonalenia konfiguracji_H2O, a następnie symulacja FFMD przeprowadzona przy stałej liczbie cząsteczek_H2O(N), objętości (V) i temperaturze (T) w celu doprowadzenia wody do temperatury symulacji, a następnie symulacja FFMD przebieg przy stałym N, ciśnieniu (P) i temperaturze (T), aby określić fizycznie poprawną wysokość pola symulacji. Pliki wyjściowe, które będą używane w sekcji 3, nazywają się data.myadsorbate_npt i log.myadsorbate.
      UWAGA: Czas trwania symulacji NPT powinien być wystarczająco długi, aby obejmował przebieg "równoważenia", w którym objętość superkomórki osiąga stan ustalony, oraz przebieg "produkcyjny", który jest używany do próbkowania średnich zespołu (w tym przypadku wysokości superkomórki). Podczas biegu równoważenia objętość superkomórki wykreślona w funkcji czasu powinna wyrównać się do wartości w stanie ustalonym. Gdy to nastąpi, można powiedzieć, że symulacja NPT jest w fazie produkcji. Sprawdź równowagę symulacji NPT, upewniając się, że wahania wysokości superkomórki (lz) są minimalne lub zbiegły się do stałej wartości. Jeśli wystąpią duże wahania, ponownie wygeneruj konfigurację_H2O, zmniejszając krok czasu w wierszu 92 w pliku input.equil i powtarzając krok 2.2.8 lub zaczynając od kroku 2.1.1 .

3. Wyodrębnij właściwą wysokość superkomórki

1. Wykonaj skrypt get_npt_lz.py w pliku log.myadsorbate, wpisując get_npt_lz.py log.myadsorbate w interfejsie wiersza polecenia. Ten skrypt wyprowadza średnią wysokość superkomórki z części "produkcyjnej" symulacji NPT w pliku avg_lz.txt.
   UWAGA: Skrypt get_npt_lz.py zakłada, że LAMMPS zapisuje długość wymiaru z komórki (lz) w pliku log.myadsorbate co 1000 fs (z możliwością dostosowania w wierszu 20 skryptu get_npt_lz.py), co jest wartością domyślną w podanym pliku wejściowym LAMMPS input.equil. Skrypt get_npt_lz.py wykrywa i odrzuca pierwsze 2 ns (konfigurowalne w wierszu 19 skryptu get_npt_lz.py) wartości lz w pliku log.myadsorbate, ponieważ składają się one na część równoważącą symulację, podczas gdy pozostałe 3 ns składają się na część "produkcyjną" i są zatem używane przez skrypt get_npt_lz.py do obliczenia średniej długości wymiaru z. Oprócz pliku avg_lz.txt skrypt get_npt_lz.py generuje plik o nazwie npt_data.txt, który zawiera wartości lz jako funkcję kroku czasowego, a także plik o nazwie npt_plot.png, który wykreśla te same dane. Wykres może być wykorzystany do weryfikacji równowagi symulacji NPT.
2. Zrekonstruuj superkomórkę, korzystając ze średniej wysokości określonej w NPT.
   1. Skopiuj plik data.myadsorbate_npt do nowego katalogu, nazywanego tutaj $WD 2PATH, i zmień jego nazwę na data.myadsorbate.
   2. Edytuj nowy plik data.myadsorbate tak, aby wysokość lz była równa średniej wartości wyjściowej ze skryptu get_npt_lz.py, zmieniając argumenty zlo i zhi w pliku data.myadsorbate w taki sposób, aby zlo było równe 0,0, a zhi było wartością lz z pliku avg_lz.txt utworzonego w kroku 3.1.

4. Generuj konfiguracje cząsteczek_H2O

1. Skopiuj plik wejściowy LAMMPS input.prod do $WD 2PATH. Edytuj zmienną grupy w wierszu 32, aby wskazać indeksy typu atomu dla atomów wody, tlenu i wodoru wody, oraz zmienną grupy w wierszu 33, aby wskazać indeksy typu atomu dla atomów Pt i adsorbatu.
2. Uruchom oprogramowanie LAMMPS, wpisując mpiexec -np XX lmp_mpi < input.prod w interfejsie wiersza poleceń, gdzie XX to liczba rdzeni procesora do użycia, a lmp_mpi to nazwa pliku wykonywalnego LAMMPS. W ten sposób uruchomiona zostanie ciągła symulacja NVT na cząsteczkach H₂O. Kluczowym plikiem wyjściowym z tej symulacji jest dump.myadsorbate.lammpstrj plik.
   UWAGA: Czas trwania symulacji NVT powinien być wystarczająco długi, aby obejmował przebieg równowagi, w którym energia układu osiąga stan ustalony, oraz przebieg produkcyjny, z którego pobierane są średnie zespoły (w tym przypadku pozycje przestrzenne cząsteczek wody). Podczas przebiegu równowagi energia układu wykreślona w funkcji czasu powinna wyrównać się do wartości w stanie ustalonym. Gdy to nastąpi, można powiedzieć, że symulacja NVT jest w fazie produkcyjnej.

5. Określ czas życia wiązania wodorowego dla prawidłowego pobierania próbek

1. Edytuj skrypt hb_lifetime_dist.py, aby określić: krok czasowy pierwszej klatki pliku dump.myadsorbate.lammpstrj poprzez zmianę zmiennej actualStart w wierszu 22, częstotliwość zapisywania ramek do pliku trajektorii LAMMPS poprzez zmianę zmiennej timestep w wierszu 23, pierwsze i ostatnie kroki czasowe, które skrypt powinien uwzględnić (tj. produkcyjna część trajektorii) poprzez zmianę zmiennych N_first i N_last w wierszach 24 i 25, czy uwzględniane są kolejne klatki, czy też klatki są pomijane przez zmianę zmiennej nevery w wierszu 26, oraz liczbę linii na sekcję klatki pliku trajektorii poprzez zmianę zmiennej frameLine w wierszu 27. Ponadto edytuj wiersze od 31 do 35, aby określić, które typy atomów w pliku data.myadsorbate należą do adsorbatu, a które typy atomów należą do cząsteczek_H2O.
   UWAGA: Skrypt hb_lifetime_dist.py analizuje konfiguracje_H2Ow cyklu produkcyjnym i określa, czy jakiekolwiek cząsteczki_H2Osą związane wodorem z adsorbatem. Następnie zlicza czas symulacji, w którym każde wiązanie wodorowe pozostaje nienaruszone, i podaje tę informację jako rozkład czasów życia wiązań wodorowych w jednostkach ps. Konkretna wersja skryptu, która jest dostarczana z tym protokołem, zakłada, że LAMMPS zapisuje konfigurację cząsteczek H₂O w pliku dump.myadsorbate.lammpstrj co 1000 fs, co jest wartością domyślną w podanym pliku wejściowym LAMMPS input.prod. Wykrywa i odrzuca pierwsze 2 ns konfiguracji w pliku dump.myadsorbate.lammpstrj, ponieważ składają się one na część symulacji dotyczącą równowagi, a pozostałe 3 ns wykorzystuje do obliczenia czasów życia wiązań wodorowych.
2. Wykonaj skrypt hb_lifetime_dist.py w pliku dump.myadsorbate.lammpstrj, wpisując hb_lifetime_dist.py w interfejsie wiersza poleceń. Spowoduje to utworzenie pliku o nazwie distribution_HB_lifetime.dat.
3. Wykreśl dane w pliku distribution_HB_lifetime.dat, aby wyświetlić rozkład czasów życia wiązań wodorowych, które wystąpiły podczas symulacji NVT.
4. Określ przyrost czasu, który ma być używany dla interwału próbkowania czasu, na podstawie obliczonych okresów życia wiązań wodorowych. Najlepszym wyborem jest maksymalna żywotność wiązania wodorowego; Alternatywnie można użyć wartości, która uchwyciłaby 95% przedział ufności.

6. Konfiguracje próbek cieczy_H2O

1. Określ liczbę konfiguracji z serii produkcyjnej trajektorii NVT FFMD w celu dalszych obliczeń. Liczba konfiguracji powinna być tak dobrana, aby minimalny czas między konfiguracjami był równy lub większy niż odstęp czasu od pobierania próbek określony w sekcji 5.
   1. Edytuj wartość domyślną zmiennej num_frames w wierszu 21 skryptu lammps_frames.py, aby określić liczbę konfiguracji do wyodrębnienia.
   2. Wykonaj skrypt lammps_frames.py w pliku dump.myadsorbate.lammpstrj, wpisując lammps_frames.py w interfejsie wiersza poleceń. Spowoduje to wyświetlenie listy czasów symulacji odpowiadających konfiguracjom, które powinny zostać wyodrębnione z pliku dump.myadsorbate.lammpstrj. Te konfiguracje mogą być używane jako struktury wyjściowe w symulacjach AIMD lub QM.
      UWAGA: 1) Skrypt lammps_frames.py automatycznie wykrywa pliki dziennika i zrzutu LAMMPS, a także część produkcyjną trajektorii w pliku zrzutu i dzieli liczbę konfiguracji w pliku zrzutu na 10 grup. Alternatywnie użytkownik może określić plik dziennika, plik zrzutu i liczbę konfiguracji z interfejsu wiersza polecenia, używając odpowiednio opcji -l, -d i -n. Aby to zrobić, użytkownik powinien wpisać lammps_frames.py -n XX -l $logfilename -d $dumpfilename w interfejsie wiersza poleceń, gdzie XX to żądana liczba konfiguracji, $logfilename to nazwa pliku dziennika LAMMPS, a $dumpfilename to nazwa pliku trajektorii (zrzutu) LAMMPS. Czasy symulacji, które są wyprowadzane, odnoszą się do mediany czasów w każdej grupie. 2) Jeśli konfiguracje będą obliczane w VASP z włączoną flagą LDIPOLE, należy dodać niewielką warstwę przestrzeni próżniowej do górnej części superkomórki nad warstwą wody. Ułatwi to konwergencję struktury elektronicznej w obliczeniach VASP. Dodanie dodatkowych 3 A przestrzeni próżniowej powyżej cząsteczek_H2Ozakończyło się sukcesem w symulacjach omówionych poniżej.

Results

Jednym z zastosowań tego protokołu jest obliczanie energii interakcji między ciekłą wodą a gatunkami katalitycznymi, tj. ΔE_int³⁵:

∆E_int=E_{Gatunki katalityczne+H2O}+E_{Czysta powierzchnia katalizatora-E Gatunki katalityczne-E Oczyść powierzchnię katalizatora+H2O}

gdzie E_{Gatunki katalityczne+H2O} to energia konfiguracji cząsteczek H₂O wokół formy katalitycznej na powierzchni metalu, E_{Czysta powierzchnia katalizatora} to energia czystej powierzchni katalizatora w próżni, E_{Cząstka katalityczna} to energia formy katalitycznej na powierzchni metalu w próżni, a E_{Czysta powierzchnia katalizatora+H2O} jest energią konfiguracji_H2Ona powierzchni katalizatora z usuniętymi formami katalitycznymi. Pozycje cząsteczek_H2Oużywanych do obliczania E_{Gatunki katalityczne + H2O} i E_{Czysta powierzchnia katalizatora + H2O} powinny być identyczne. Wszystkie wartości E są obliczane przy użyciu kodu VASP. Wielkość ΔE_int obejmuje wszystkie fizyczne i chemiczne oddziaływania między wszystkimi cząsteczkami w strukturze ciekłej wody a formami katalitycznymi i daje rozsądne oszacowanie entalpii solwatacji związków katalitycznych, które jest potrzebne do obliczenia jego energii swobodnej solwatacji i całkowitej energii swobodnej. W tabeli 1 przedstawiono wartości ΔE_int obliczone dla gatunków na powierzchni katalizatora Pt(111) o wzorach chemicznych równych C_xH_yO_z w jednostkach eV (1 eV = 96,485 kJ/mol). Wartości zostały obliczone przy zakresach ≤1/9 ML.^35,46 Podane wartości są średnimi przyjętymi dla 10 konfiguracji cieczy_H2O, a niepewności są podawane jako odchylenia standardowe. Wszystkie wartości są ujemne, co wskazuje na korzystne interakcje z wodą.

Innym zastosowaniem tego protokołu jest generowanie struktur startowych dla AIMD. Film 1 to film przedstawiający trajektorię AIMD, która została uruchomiona z konfiguracji wygenerowanej przez ten protokół. Na początku tego filmu adsorbat COH jest pokazany na powierzchni Pt(111) pod strukturą cieczy H₂O. Podkreślono jedną cząsteczkę_H2O, która utworzyła wiązanie wodorowe z COH. W trakcie filmu ta cząsteczka_H2Oabstrahuje proton z adsorbatu COH i osadza drugi atom wodoru na powierzchni Pt(111). Cząsteczka_H2Opomaga zatem katalizować reakcję COH* + * → CO* + H*, gdzie *s wskazują miejsca katalityczne. Symulacja ta uwypukla główną siłę i główny cel opisanej w niniejszym dokumencie wieloskalowej metody pobierania próbek. Liczne konfiguracje cząsteczek H₂O są generowane za pomocą FFMD, ze względu na jego siłę w zakresie wykonalności obliczeniowej. Jednak ograniczeniem FFMD jest to, że nie może uchwycić zrywania i tworzenia wiązań, chyba że zaimplementowane jest pole sił reaktywnych. AIMD wykorzystuje mechanikę kwantową do obliczania energii, a tym samym może uchwycić zrywanie i tworzenie wiązań. Jednak AIMD jest zbyt wymagający obliczeniowo, aby wygenerować wszystkie konfiguracje cząsteczek_H2Oniezbędne do zapewnienia wystarczającej ilości próbkowania. W ten sposób protokół ten łączy te dwie metody.

Struktury ciekłych cząsteczek H₂O generowanych przez tę procedurę zależą od ustawień wejściowych. Niewłaściwe ich ustawienie może mieć niezamierzony wpływ na struktury wodne. Na przykład, gdy odległości międzycząsteczkowe stają się zbyt małe lub gdy inne parametry w plikach wejściowych dynamiki molekularnej są ustawione nieprawidłowo lub przyjmują niefizyczne wartości, struktura wody może stać się nieuzasadniona. W takich okolicznościach struktura wody "eksploduje" w sposób niezamierzony podczas trajektorii FFMD. Rysunek 1 pokazuje przykład tego. Migawka po lewej stronie to struktura początkowa dla przebiegu FFMD, a migawka po prawej stronie to migawka wykonana w ciągu 1 ps od rozpoczęcia symulacji. Jak widać, cząsteczki_H2Oprzesunęły się daleko od powierzchni. Jest to spowodowane nieprawidłowymi ustawieniami wprowadzonymi w plikach wejściowych symulacji i nie jest to struktura, która prawdopodobnie wystąpi w rzeczywistości.

Rysunek 1: Przykład wyniku ujemnego. Symulacja dynamiki molekularnej pola siłowego "wybucha" z powodu niefizycznego ustawienia lub wartości. Obraz po lewej stronie: Początkowa geometria powierzchni Pt(111), adsorbatu i struktury ciekłej wody. Obraz po prawej stronie: Geometria powierzchni Pt(111), adsorbatu i struktury ciekłej wody mniej niż 1 ps później. Na zdjęciu po prawej stronie cząsteczki_H2Ooddzieliły się od powierzchni z powodu niefizycznie dużych sił. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Film 1: Symulacja dynamiki molekularnej ab initio (AIMD) zainicjowana na podstawie konfiguracji wygenerowanej w próbkowaniu wieloskalowym. Cząsteczka H₂O, która jest pierwotnie wodorem związanym z adsorbatem COH na powierzchni Pt(111), abstrahuje proton od COH i osadza drugi wodór na powierzchni Pt(111). To zdarzenie zerwania i tworzenia wiązania może być uchwycone przez AIMD, ale nie za pomocą dynamiki molekularnej pola siłowego (FFMD), chyba że używane jest pole sił reaktywnych. Początkowa konfiguracja cząsteczek_H2Oużytych w tej symulacji AIMD została wygenerowana przy użyciu FFMD, jak opisano w tym manuskrypcie. Kliknij tutaj, aby obejrzeć ten film. (Kliknij prawym przyciskiem myszy, aby pobrać.)

Gatunki katalityczne	∆Eint (eV)
Coh	-0,70 ± 0,07
współ	-0,03 ± 0,03
CH2OH	-0,64 ± 0,12
CHO-CHOH-CH2OH	-0,93 ± 0,22
COH-COH-CH2OH	-0,87 ± 0,23
COH-CHOH-COH	-1,72 ± 0,26
CHOH-COH-CO	-1,57 ± 0,25
CZEKOLADA-CO-CO	-0,31 ± 0,19

Tabela 1: Wyniki energii interakcji między wodą a katalizatorem. Energie oddziaływań w eV obliczone dla ośmiu adsorbatów C_xH_yO_z na Pt(111). Podane wartości są średnimi uzyskanymi dla wielu konfiguracji cieczy_H2O. Niepewności są odchyleniami standardowymi średnich. 1 eV = 96,485 kJ/mol.

Discussion

Przedstawiona metoda została wybrana ze względu na łatwość wdrożenia, ale można dokonać wielu dostosowań. Po pierwsze, pola siłowe używane w symulacjach FFMD mogą być modyfikowane. Zmianę parametrów i/lub potencjałów pola sił można dokonać poprzez edycję plików wejściowych i danych LAMMPS. Podobnie można stosować rozpuszczalniki inne niż_H2O. Aby dokonać tej modyfikacji, pożądana cząsteczka rozpuszczalnika musiałaby zostać wstawiona, począwszy od kroku 2.1.1, a pliki wejściowe LAMMPS musiałyby zostać zmodyfikowane w celu uwzględnienia odpowiednich potencjałów i parametrów. Wprowadzenie nowej cząsteczki rozpuszczalnika wymagałoby również podania wewnętrznych współrzędnych cząsteczki rozpuszczalnika w pliku .txt analogicznym do pliku water.txt.

Kolejną modyfikacją, którą można wprowadzić, jest modyfikacja powierzchni płyty. W wyniku omówionego w tym manuskrypcie wykorzystano płyty powierzchniowe 3 Pt x 3 Pt lub 4 Pt x 4 Pt, których powierzchnia jest mniejsza niż 120 A². Wraz ze wzrostem powierzchni płyty wzrasta również koszt obliczeniowy. Koszt obliczeniowy ma największy wpływ na sekcję 5 tego protokołu. Jeżeli etapy przetwarzania danych w sekcji 5 staną się niewystarczające obliczeniowo, można zastosować strategie przetwarzania końcowego dużych zbiorów danych, takie jak te omówione w Li i in. 2018⁴⁵ .

Możliwe źródła niepewności dla tej procedury obejmują zastosowane pole siłowe, metodę próbkowania i częstotliwość próbkowania. Struktura wody jest określana przez użyte pole siłowe, co oznacza, że wybór pola siłowego może wpływać na określone konfiguracje cząsteczek_H2O. Nasza grupa oceniła, w jaki sposób wybór pola sił dla cząsteczek_H2Oi atomów Pt wpływa na energie oddziaływań obliczone w FFMD i odkryła, że wybór pola siłowego przyczynia się do mniej niż 0,1 eV do tej energii oddziaływania. Innym źródłem niepewności jest metoda doboru próby, która wpływa na określone konfiguracje używane do obliczania interesującej nas wielkości. Nasza grupa porównała wydajność metody "próbkowania w czasie" przedstawionej w tym protokole z metodą "próbkowania energii", która jest stronnicza w stosunku do niższych konfiguracji energetycznych cząsteczek_H2O, na energiach oddziaływań obliczonych w DFT i stwierdziła, że obie te metody próbkowania dają statystycznie równe wartości^35,46. Częstotliwość próbkowania może również mieć wpływ na wyniki. Oceniliśmy, w jaki sposób zwiększenie liczby konfiguracji od 10 do 30 000 wpływa na średnie energie oddziaływań obliczone w FFMD dla 40 różnych adsorbatów_C3H_xO3 i stwierdziliśmy, że częstotliwość próbkowania przyczynia się mniej niż 0,1 eV do średniej energii oddziaływania⁴⁴.

Głównym ograniczeniem tej metody jest to, że adsorbaty są aproksymowane przez struktury w próżni podczas symulacji FFMD. W rzeczywistości adsorbaty wykazywałyby zmiany konformacyjne (rozciąganie wiązań, zgięcia kątowe, ruchy skręcające itp.) spowodowane normalnymi ruchami termicznymi, w tym interakcjami z cząsteczkami rozpuszczalnika. Próby włączenia zmian konformacyjnych adsorbatów do symulacji FFMD wymagałyby szczegółowego opracowania pól siłowych dla katalitycznych adsorbatów powierzchniowych, tj. które zawierają między innymi terminy opisujące rozciąganie wiązań, zagięcia kątowe i terminy skręcania. Jako przyszły kierunek tego protokołu, opracowujemy takie pola siłowe dla adsorbatów na powierzchniach stałych, których użyjemy do określenia zakresu, w jakim użycie sztywnych adsorbatów wpływa na wyniki.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Oprogramowanie VASP	Computational Materials Physics, Wydział Fizyki, Uniwersytet Wiedeński	vasp.5.4.4	Standardowy równoległy plik wykonywalny VASP w najnowszej wersji.
Oprogramowanie LAMMPS	Sandia National Laboratory	31Mar17-dp	Plik wykonywalny LAMMPS o podwójnej precyzji od 31 marca 2017 r.
Oprogramowanie VMD	Theoretical and Computational Biophysics Group, University of Illinois at Urbana-Champaign	1.9.3	Standardowy plik wykonywalny VMD w najnowszej wersji.
Oprogramowanie MCPliQ	Getman Research Group, Dept. of Chemical and Biomolecular Engineering, Clemson University		Pliki wykonywalne i wejściowe dla oprogramowania MCPliQ dostępne na stronie GitHub Getman Research Group.
Skrypty artykułów JoVE	Getman Research Group, Dept. of Chemical and Biomolecular Engineering, Clemson University		Skrypty Pythona dla tego manuskryptu JoVE dostępne na stronie GitHub Getman Research Group.
Plik PDB H2O	Getman Research Group, Dept. of Chemical and Biomolecular Engineering, Clemson University lub RCSB Protein Data Bank		plik PDB dla cząsteczki wody, dostępny na stronie GitHub Getman Research Group lub pod adresem http://www.rcsb.org/ligand/HOH.