Obliczanie stężeń atmosferycznych klastrów molekularnych na podstawie termochemii ab initio

Podsumowanie

Stężenia atmosferyczne słabo związanych klastrów molekularnych można obliczyć na podstawie właściwości termochemicznych struktur o niskiej energii, znalezionych za pomocą wieloetapowej metodologii próbkowania konfiguracyjnego wykorzystującej algorytm genetyczny oraz półempiryczną i ab initio chemię kwantową.

Streszczenie

Obliczeniowe badanie powstawania i wzrostu aerozoli atmosferycznych wymaga dokładnej powierzchni energii swobodnej Gibbsa, którą można uzyskać z obliczeń struktury elektronowej fazy gazowej i częstotliwości drgań. Wielkości te są ważne dla tych klastrów atmosferycznych, których geometria odpowiada minimum na ich powierzchniach energii potencjalnej. Energia swobodna Gibbsa o strukturze minimalnej energii może być wykorzystana do przewidywania stężeń atmosferycznych gromady w różnych warunkach, takich jak temperatura i ciśnienie. Przedstawiamy tanią obliczeniowo procedurę opartą na próbkowaniu konfiguracyjnym opartym na algorytmie genetycznym, po którym następuje seria coraz dokładniejszych obliczeń przesiewowych. Procedura rozpoczyna się od wygenerowania i rozwinięcia geometrii dużego zestawu konfiguracji przy użyciu modeli półempirycznych, a następnie udoskonala powstałe unikalne struktury na szeregu wysokopoziomowych poziomów teorii ab initio. Na koniec obliczane są poprawki termodynamiczne dla wynikowego zestawu struktur o minimalnych energiach energetycznych i wykorzystywane do obliczania energii swobodnych Gibbsa powstawania, stałych równowagi i stężeń atmosferycznych. Przedstawiamy zastosowanie tej procedury do badania uwodnionych klastrów glicyny w warunkach otoczenia.

Wprowadzenie

Najbardziej niepewnym parametrem w badaniach atmosferycznych zmian klimatycznych jest dokładny zakres, w jakim cząstki chmur odbijają nadchodzące promieniowanie słoneczne. Aerozole, które są cząstkami stałymi zawieszonymi w gazie, tworzą cząstki chmur zwane jądrami kondensacji chmur (CCN), które rozpraszają przychodzące promieniowanie, zapobiegając w ten sposób jego absorpcji i późniejszemu podgrzaniu atmosfery¹. Szczegółowe zrozumienie tego efektu chłodzenia netto wymaga zrozumienia wzrostu aerozoli do CCN, co z kolei wymaga zrozumienia wzrostu małych klastrów molekularnych w cząsteczki aerozolu. Ostatnie prace sugerują, że powstawanie aerozoli jest inicjowane przez klastry molekularne o średnicy 3 nm lub mniejszej²; Jednak ten reżim rozmiarów jest trudny do uzyskania przy użyciu technik eksperymentalnych^3,4. W związku z tym pożądane jest podejście do modelowania obliczeniowego w celu przezwyciężenia tego eksperymentalnego ograniczenia.

Korzystając z naszego podejścia do modelowania opisanego poniżej, możemy przeanalizować wzrost dowolnej uwodnionej gromady. Ponieważ interesuje nas rola wody w tworzeniu dużych cząsteczek biologicznych z mniejszych składników w środowiskach prebiotycznych, ilustrujemy nasze podejście glicyną. Napotkane wyzwania i narzędzia potrzebne do odpowiedzi na te pytania badawcze są bardzo podobne do tych, które są związane z badaniem aerozoli atmosferycznych i klastrów prenukleacji^{5,6,7,8,9,10,11,12, 13,14,15}. Tutaj badamy uwodnione klastry glicyny, zaczynając od wyizolowanej cząsteczki glicyny, po której następuje seria stopniowych dodawania do pięciu cząsteczek wody. Ostatecznym celem jest obliczenie stężeń równowagowych klastrów Gly(H2O)_n=0-5 w atmosferze w temperaturze pokojowej na poziomie morza i wilgotności względnej (RH) 100 %.

Niewielka liczba tych sub-nanometrowych klastrów molekularnych rozwija się w metastabilny klaster krytyczny (o średnicy 1-3 nm) albo przez dodanie innych cząsteczek pary, albo przez koagulację na istniejących klastrach. Te krytyczne gromady mają korzystny profil wzrostu prowadzący do powstawania znacznie większych (do 50-100 nm) jąder kondensacji chmur (CCN), które bezpośrednio wpływają na efektywność opadów atmosferycznych chmur, a także na ich zdolność do odbijania padającego światła. W związku z tym dobre zrozumienie termodynamiki klastrów molekularnych i ich rozkładów równowagi powinno prowadzić do dokładniejszych przewidywań wpływu aerozoli na klimat globalny.

Opisowy model powstawania aerozoli wymaga dokładnej termodynamiki powstawania klastrów molekularnych. Obliczenie dokładnej termodynamiki powstawania klastrów molekularnych wymaga identyfikacji najbardziej stabilnych konfiguracji, co wiąże się ze znalezieniem globalnych i lokalnych minimów na powierzchni energii potencjalnej klastra (PES)¹⁶. Proces ten nazywa się próbkowaniem konfiguracyjnym i można go osiągnąć za pomocą różnych technik, w tym opartych na dynamice molekularnej (MD)^17,18,19,20, Monte Carlo (MC)^21,22 oraz algorytmy genetyczne (GA)^23,24,25.

Na przestrzeni lat opracowano różne protokoły, aby uzyskać strukturę i termodynamikę hydratów atmosferycznych na wysokim poziomie teorii. Protokoły te różniły się wyborem (i) metody próbkowania w konfiguracji (ii) charakteru metody niskiego poziomu stosowanej w konfiguracyjnym doborze próbek oraz (iii) hierarchii metod wyższego poziomu stosowanych do doprecyzowania wyników w kolejnych etapach.

Konfiguracyjne metody próbkowania obejmowały intuicję chemiczną²⁶, losowe pobieranie próbek^27,28, dynamika molekularna (MD)^29,30, skakanie po basenie (BH)³¹, oraz algorytm genetyczny (GA)^24,25,32. Najczęstszymi metodami niskopoziomowymi stosowanymi w tych metodach próbkowania są pola sił lub modele półempiryczne, takie jak PM6, PM7 i SCC-DFTB. Po nich często następują obliczenia DFT z coraz większymi zestawami bazowymi i bardziej niezawodnymi funkcjami z wyższych szczebli drabiny Jakuba³³. W niektórych przypadkach są one stosowane przez metody funkcji falowych wyższego poziomu, takie jak MP2, CCSD(T) i ekonomiczne DLPNO-CCSD(T)^34,35.

Kildgaard et al.³⁶ opracowali systematyczną metodę, w której cząsteczki wody są dodawane w punktach sfer Fibonacciego³⁷ wokół mniejszych uwodnionych lub nieuwodnionych klastrów, aby wygenerować kandydatów na większe klastry. Niefizyczni i nadmiarowi kandydaci są usuwani na podstawie progów bliskiego kontaktu i średniej kwadratowej odległości między różnymi konformerami. Kolejne optymalizacje przy użyciu metody półempirycznej PM6 oraz hierarchii metod DFT i funkcji falowych są wykorzystywane do uzyskania zestawu konformerów o niskiej energii na wysokim poziomie teorii.

Algorytm sztucznej kolonii pszczół (ABC)³⁸ to nowe podejście do próbkowania konfiguracyjnego, które zostało niedawno wdrożone przez Zhang et al. do badania klastrów molekularnych w programie o nazwie ABCluster³⁹. Kubecka et al.⁴⁰ użyli ABCluster do próbkowania konfiguracyjnego, a następnie do ponownych optymalizacji niskiego poziomu przy użyciu ściśle wiążącej metody półempirycznej GFN-xTB⁴¹. Następnie udoskonalili struktury i energie za pomocą metod DFT, a następnie energii końcowej przy użyciu DLPNO-CCSD(T).

Niezależnie od metody, próbkowanie konfiguracyjne rozpoczyna się od losowo lub nielosowo wygenerowanego rozkładu punktów w PES. Każdy punkt odpowiada określonej geometrii danego klastra molekularnego i jest generowany metodą pobierania próbek. Następnie dla każdego punktu znajduje się najbliższe minimum lokalne, podążając w kierunku "w dół" na PES. Znaleziony w ten sposób zbiór minimów odpowiada tym geometriom klastra molekularnego, które są stabilne, przynajmniej przez jakiś czas. W tym przypadku kształt PES i ocena energii w każdym punkcie na powierzchni będą wrażliwe na fizyczny opis systemu, podczas gdy dokładniejszy opis fizyczny skutkuje bardziej kosztownymi obliczeniowo obliczeniami energii. W szczególności użyjemy metody GA zaimplementowanej w programie OGOLEM²⁵, która została z powodzeniem zastosowana do różnych globalnych problemów z optymalizacją i konfiguracyjnym próbkowaniem^42,43,44,45, aby wygenerować początkowy zestaw punktów próbkowania. PES będzie opisany przez model PM7 ⁴⁶ zaimplementowany w programie MOPAC2016⁴⁷. Ta kombinacja jest stosowana, ponieważ generuje większą różnorodność punktów w porównaniu z metodami MD i MC i znajduje lokalne minima szybciej niż bardziej szczegółowe opisy PES.

Zestaw lokalnych minimów zoptymalizowanych pod kątem GA jest traktowany jako geometrie wyjściowe dla serii etapów przesiewania, które prowadzą do zestawu nisko położonych minimalnych energii. Ta część protokołu rozpoczyna się od optymalizacji zestawu unikalnych struktur zoptymalizowanych pod kątem GA przy użyciu teorii funkcjonału gęstości (DFT) z małym zestawem podstawowym. Ten zestaw optymalizacji na ogół daje mniejszy zestaw unikalnych lokalnych struktur minimalnych, które są modelowane bardziej szczegółowo w porównaniu ze strukturami półempirycznymi zoptymalizowanymi pod kątem GA. Następnie kolejna runda optymalizacji DFT jest wykonywana na tym mniejszym zestawie struktur przy użyciu większego zestawu podstawowego. Ponownie, ten krok zazwyczaj daje mniejszy zestaw unikalnych struktur, które są modelowane bardziej szczegółowo w porównaniu z małym krokiem DFT na podstawie. Ostateczny zestaw unikalnych struktur jest następnie optymalizowany pod kątem ściślejszej konwergencji i obliczane są harmoniczne częstotliwości drgań. Po tym kroku mamy wszystko, czego potrzebujemy, aby obliczyć stężenia równowagowe gromad w atmosferze. Ogólne podejście jest podsumowane diagramowo w Rysunek 1. Użyjemy PW91⁴⁸ funkcjonału korelacji wymiany uogólnionego gradientu (GGA) w implementacji DFT Gaussian09⁴⁹ wraz z dwoma wariantami zestawu bazowego Pople50 (6-31+G* dla małego kroku bazowego i 6-311++G** dla dużego kroku bazowego). Ta szczególna kombinacja funkcjonału wymiany i korelacji oraz zestawu bazowego została wybrana ze względu na jej wcześniejsze sukcesy w obliczaniu dokładnych energii swobodnych formowania Gibbsa dla gromad atmosferycznych^51,52.

Ten protokół zakłada, że użytkownik ma dostęp do klastra obliczeń o wysokiej wydajności (HPC) z przenośnym systemem wsadowym⁵³ (PBS), MOPAC2016 (http://openmopac.net/MOPAC2016.html)⁴⁷, OGOLEM (https://www.ogolem.org)²⁵, Gaussian 09 (https://gaussian.com)⁴⁹, oraz OpenBabel⁵⁴ (http://openbabel.org/wiki/Main_Page) oprogramowanie zainstalowane zgodnie z instrukcjami instalacji. Każdy krok w tym protokole wykorzystuje również zestaw wewnętrznych skryptów powłoki i Pythona 2.7, które muszą być zapisane w katalogu zawartym w zmiennej środowiskowej $PATH użytkownika. Wszystkie niezbędne moduły środowiskowe i uprawnienia do wykonywania do uruchomienia wszystkich powyższych programów muszą być również załadowane do sesji użytkownika. Użycie dysku i pamięci przez kod GA (OGOLEM) i kody półempiryczne (MOPAC) są bardzo małe jak na współczesne standardy zasobów komputerowych. Ogólne użycie pamięci i dysku przez OGOLEM/MOPAC zależy od tego, ile wątków chce się użyć, a nawet wtedy zużycie zasobów będzie niewielkie w porównaniu z możliwościami większości systemów HPC. Potrzeby w zakresie zasobów metod zarządzania jakością zależą od wielkości klastrów i zastosowanego poziomu teorii. Zaletą korzystania z tego protokołu jest to, że można zmieniać poziom teorii, aby móc obliczyć końcowy zestaw struktur niskoenergetycznych, pamiętając, że zwykle szybsze obliczenia prowadzą do większej niepewności co do dokładności wyników.

Dla jasności, lokalny komputer użytkownika będzie określany jako "komputer lokalny", podczas gdy klaster HPC, do którego mają dostęp, będzie określany jako "zdalny klaster".

Protokół

1. Znalezienie minimalnej struktury energetycznej izolowanej glicyny i wody

UWAGA: Cel jest tutaj dwojaki: (i) uzyskanie minimalnych struktur energetycznych izolowanych cząsteczek wody i glicyny do wykorzystania w próbkowaniu konfiguracji algorytmu genetycznego, (ii) oraz obliczenie poprawek termodynamicznych do energii fazy gazowej tych cząsteczek do wykorzystania w obliczeniach stężeń atmosferycznych.

1. Na komputerze lokalnym otwórz nową sesję Avogadro.
   1. Kliknij opcję Build > Insert > Peptide i wybierz opcję Gly w oknie Insert Peptide, aby wygenerować monomer glicyny w oknie wizualizacji.
   2. Kliknij pozycję Rozszerzenia > Gaussa i edytuj pierwszy wiersz w polu tekstowym, aby brzmiał '# pw91pw91/6-311++G** int(Acc2E=12,UltraFine) scf(conver=12) opt(tight,maxcyc=300) freq'. Kliknij przycisk Generuj i zapisz plik wejściowy w glycine.com.
   3. Należy pamiętać, że jeśli cząsteczka ma znaczną elastyczność konformacyjną, tak jak glicyna⁵⁵, bardzo ważne jest przeprowadzenie analizy konformacyjnej w celu zidentyfikowania globalnej struktury minimalnej i innych nisko położonych konformerów. OpenBabel⁵⁴ dostarcza solidnych narzędzi do wyszukiwania konformacyjnego, wykorzystujących różne algorytmy i szybkie pola siłowe. Podczas gdy konformiści mogą się rozluźniać i konwertować podczas GA i kolejnych obliczeń, czasami konieczne jest uruchomienie wielu obliczeń GA, z których każde zaczyna się od innego konformera.
2. Na komputerze lokalnym otwórz nową sesję w Avogadro.
   1. Kliknij pozycję Build > Insert > Fragment i wyszukaj "water" w oknie Insert Fragment, aby uzyskać współrzędne wody.
   2. Kliknij pozycję Rozszerzenia > Gaussa i edytuj pierwszy wiersz w polu tekstowym, aby brzmiał '# pw91pw91/6-311++G** int(Acc2E=12,UltraFine) scf(conver=12) opt(tight,maxcyc=300) freq'. Kliknij przycisk Generuj i zapisz plik wejściowy jako water.com.
3. Przenieś dwa pliki .com do klastra zdalnego. Po zalogowaniu się do klastra zdalnego wywołaj Gaussian 09 w skrypcie przesyłania wsadowego, aby rozpocząć obliczenia. Po zakończeniu obliczeń wyodrębnij współrzędne kartezjańskie (pliki .xyz) struktur o minimalnej energii, wywołując OpenBabel. W przypadku glicyny polecenie do wykonania to:
   obabel -ig09 glycine.log -oxyz > glycine.xyz
   Te dwa pliki >.xyz będą używane przez próbkowanie konfiguracyjne ogólnie dostępne w następnym kroku.

2. Oparte na algorytmie genetycznym próbkowanie konfiguracyjne klastrów Gly(H₂O)_n=1-5

UWAGA: Celem tutaj jest uzyskanie zestawu niskoenergetycznych struktur dla Gly(H₂O)_n=1-5 na niedrogim półempirycznym poziomie teorii, przy użyciu modelu PM7⁴⁶ zaimplementowanego w MOPAC⁴⁷. Konieczne jest, aby katalog roboczy miał dokładnie taką organizację i strukturę, jak pokazano na Rysunek 2. Ma to na celu zapewnienie, że niestandardowa powłoka i skrypty języka Python działają bez błędów.

1. Skopiuj wszystkie niezbędne skrypty do klastra zdalnego i dodaj ich lokalizację do $PATH
   1. Umieść wszystkie skrypty i pliki szablonów w folderze (np. skrypty) i skopiuj je do zdalnego klastra
   2. Upewnij się, że wszystkie skrypty są wykonywalne
   3. Dodaj lokalizację katalogu scripts do zmiennej środowiskowej $PATH, wprowadzając następujące polecenia w terminalu. Domyślna lokalizacja skryptów to $HOME/JoVE-demo/scripts, jednak można zdefiniować zmienną środowiskową o nazwie $SCRIPTS_HOME wskazującą na katalog zawierający skrypty i dodać $SCRIPTS_HOME do swojej ścieżki
      1. Powłoka bash:
         export SCRIPTS_HOME=/ścieżka/do/skryptów
         export PATH=${SCRIPTS_HOME}:${PATH}
      2. Powłoka Tcsh / Csh:
         setenv SCRIPTS_HOME /ścieżka/do/skryptów
         setenv PATH ${SCRIPTS_HOME}:${PATH}
2. W klastrze zdalnym skonfiguruj i uruchom obliczenie ogólnie dostępne:
   1. Utwórz katalog o nazwie gly-h2o-n, gdzie n to liczba cząsteczek wody.
   2. Utwórz podkatalog o nazwie GA w katalogu gly-h2o-n, aby uruchomić obliczenia algorytmu genetycznego.
   3. Skopiuj pliki wejściowe OGOLEM (np. pm7.ogo), monomery, współrzędne kartezjańskie (np. glycine.xyz, water.xyz) i skrypt przesyłania wsadowego PBS (np. run.pbs) do katalogu GA.
   4. Wprowadź niezbędne zmiany w pliku wejściowym OGOLEM i pliku przesyłania partii.
   5. Prześlij obliczenie. Po rozpoczęciu obliczeń OGOLEM utworzy nowy katalog o nazwie prefiks pliku wejściowego OGOLEM (Np. pm7) w katalogu GA i zapisze w nim nowo wygenerowane współrzędne.
3. Po zakończeniu obliczeń skompiluj energie i stałe rotacji i wykorzystaj te informacje do określenia, które są unikalnymi strukturami niskoenergetycznymi:
   1. Zmień katalog na gly-h2o-n/GA/pm7 i
   2. Wyodrębnij energie i oblicz stałe rotacji klastrów zoptymalizowanych pod kątem GA za pomocą polecenia:
      getRotConsts-GA.csh N 0 99
      gdzie N to liczba atomów w klastrze molekularnym, a "0 99" wskazuje, że wielkość puli GA wynosi 100, ze wskaźnikami od 0 do 99. Spowoduje to wygenerowanie pliku o nazwie rotConstsData_C, który zawiera posortowaną listę wszystkich konfiguracji klastrów zoptymalizowanych pod kątem GA, ich energii i stałych obrotowych.
   3. Wykonaj polecenie:
      similarityAnalysis.py PM7 rotConstsData_C
      gdzie pm7 będzie używany jako etykieta nazwy plików, aby znaleźć i zapisać unikalne klastry zoptymalizowane pod kątem GA. Spowoduje to wygenerowanie pliku o nazwie uniqueStructures-pm7.data, który zawiera posortowaną listę unikalnych konfiguracji zoptymalizowanych pod kątem GA. Jest to lista unikalnych lokalnych struktur minimalnych dla klastra Gly(H2O)_n zoptymalizowanych na poziomie teorii PM7, a struktury te są teraz gotowe do udoskonalenia przy użyciu DFT.
4. Przejdź do katalogu gly-h2o-n/GA i połącz wyniki z wielu porównywalnych przebiegów ogólnie dostępnych przy użyciu skryptu combine-GA.csh. Składnia jest następująca:
   combine-GA.csh
   W tym konkretnym przypadku polecenie:
   kombinować-GA.csh pm7 pm7
   wygeneruje nową listę unikalnych struktur o nazwie 'uniqueStructures-pm7.data' w katalogu gly-h2o-n/GA.

3. Udoskonalanie za pomocą metody QM z małym zestawem podstaw

UWAGA: Celem tutaj jest udoskonalenie próbkowania konfiguracyjnego klastrów Gly(H₂O)_n=1-5 przy użyciu lepszego opisu kwantowo-mechanicznego, aby uzyskać mniejszy, ale dokładniejszy zestaw struktur klastrowych Gly(H₂O)_n=1-5. Strukturami wyjściowymi dla tego kroku są dane wyjściowe kroku 2.

1. Przygotuj i uruchom obliczenie DFT dla małego zestawu podstawowego:
   1. Utwórz podkatalog o nazwie QM w katalogu gly-h2o-n. W katalogu QM utwórz kolejny podkatalog o nazwie pw91-sb.
   2. Skopiuj listę unikalnych struktur (uniqueStructures-pm7.data) z katalogu gly-h2o-n/GA do katalogu QM/pw91-sb.
   3. Zmień katalog na ten gly-h2o-n/QM/pw91-sb.
   4. Uruchom skrypt próbkowania konfiguracyjnego DFT z małym zestawem podstawowym za pomocą polecenia:
      run-pw91-sb.csh uniqueStructures-pm7.data sb KOLEJKA 10
      gdzie sb jest etykietą dla tego zestawu obliczeń, QUEUE jest preferowaną kolejką w klastrze obliczeniowym, a 10 oznacza, że 10 obliczeń ma zostać zgrupowanych w jedno zadanie wsadowe. Ten skrypt automatycznie wygeneruje dane wejściowe dla Gaussian 09 i prześle wszystkie obliczenia. Wprowadź "test" dla "QUEUE", aby wykonać próbny przebieg.
2. Po zakończeniu przesłanych obliczeń wyodrębnij i przeanalizuj wyniki.
   1. Wyodrębnij energie i oblicz stałe rotacji klastrów zoptymalizowanych pod kątem małych baz za pomocą polecenia:
      getRotConsts-dft-sb.csh pw91 N
      gdzie pw91 wskazuje, że użyto funkcjonału gęstości PW91, a N to liczba atomów w klastrze. Spowoduje to utworzenie pliku o nazwie rotConstsData_C.
   2. Teraz zidentyfikuj unikalne struktury za pomocą polecenia:
      similarityAnalysis.py SB rotConstsData_C
      Gdzie sb jest używany jako etykieta nazwy pliku. W pliku uniqueStructures-sb.data zostanie teraz zapisana lista unikalnych konfiguracji zoptymalizowanych na poziomie teorii PW91/6-31+G*.
3. Przejdź do katalogu gly-h2o-n/QM i połącz wyniki z wielu porównywalnych przebiegów QM za pomocą skryptu combine-QM.csh. Składnia jest następująca:
   combine-QM.csh
   W tym konkretnym przypadku polecenie:
   combine-QM.csh sb pw91-sb
   wygeneruje nową unikalną listę struktur o nazwie 'uniqueStructures-sb.data' w katalogu gly-h2o-n/QM.

4. Dalsze udoskonalanie za pomocą metody QM z dużym zestawem podstaw

UWAGA: Celem tutaj jest dalsze udoskonalenie próbkowania konfiguracyjnego klastrów Gly(H₂O)_n=1-5 przy użyciu lepszego opisu kwantowo-mechanicznego. Strukturami wyjściowymi dla tego kroku są dane wyjściowe kroku 3.

1. Przesyłaj bardziej wiarygodne obliczenia przy użyciu większego zestawu podstawowego.
   1. Utwórz podkatalog o nazwie pw91-lb w katalogu QM.
   2. Skopiuj listę unikalnych struktur (uniqueStructures-sb.data) z katalogu gly-h2o-n/QM do katalogu gly-h2o-n/QM/pw91-lb i przejdź do tego katalogu.
   3. Uruchom skrypt próbkowania konfiguracyjnego DFT o dużej podstawie za pomocą polecenia:
      run-pw91-lb.csh uniqueStructures-sb.data lb KOLEJKA 10
      gdzie lb jest etykietą dla tego zestawu obliczeń, QUEUE jest preferowaną kolejką w klastrze obliczeniowym, a 10 oznacza, że 10 obliczeń ma zostać zgrupowanych w jedno zadanie wsadowe. Ten skrypt automatycznie wygeneruje dane wejściowe dla Gaussian 09 i prześle wszystkie obliczenia. Wprowadź wartość "test" dla pola "QUEUE", aby przeprowadzić test na sucho.
2. Po zakończeniu przesłanych obliczeń wyodrębnij i przeanalizuj dane
   1. Obliczanie stałych obrotowych klastrów zoptymalizowanych pod kątem dużych podstaw za pomocą polecenia:
      getRotConsts-dft-lb.csh pw91 N
      gdzie pw91 wskazuje, że użyto funkcjonału gęstości PW91, a N to liczba atomów w klastrze.
   2. Teraz zidentyfikuj unikalne struktury za pomocą polecenia:
      similarityAnalysis.py funty rotConstsData_C
      gdzie lb jest używany jako etykieta nazwy pliku. Masz teraz listę unikalnych konfiguracji zoptymalizowanych na poziomie teorii PW91/6-311++G** zapisaną w pliku uniqueStructures-lb.data.

5. Obliczenia końcowej energii i korekcji termodynamicznej

UWAGA: Celem tutaj jest uzyskanie struktury wibracyjnej i energii klastrów Gly(H₂O)_n=1-5 przy użyciu dużego zestawu bazowego i ultradrobnej siatki integracji w celu obliczenia pożądanych poprawek termochemicznych.

1. Zaczynając od wyników z poprzedniego kroku, prześlij bardziej wiarygodne obliczenia.
   1. Utwórz podkatalog o nazwie ultrafine w katalogu QM/pw91-lb. Następnie skopiuj listę unikalnych struktur (uniqueStructures-lb.data) z katalogu QM/pw91-lb do katalogu QM/pw91-lb/ultrafine i przejdź do tego katalogu.
   2. Prześlij ultradrobny, duży skrypt DFT za pomocą polecenia:
      run-pw91-lb-ultrafine.csh uniqueStructures-lb.data uf KOLEJKA 10
      gdzie uf jest etykietą dla tego zestawu obliczeń, QUEUE jest preferowaną kolejką w klastrze obliczeniowym, a 10 oznacza, że 10 obliczeń ma zostać zgrupowanych w jedno zadanie wsadowe. Ten skrypt automatycznie wygeneruje dane wejściowe dla Gaussian 09 i prześle wszystkie obliczenia. Wprowadź wartość "test" dla pola "QUEUE", aby przeprowadzić test na sucho.
2. Po zakończeniu przesłanych obliczeń wyodrębnij i przeanalizuj dane
   1. Wyodrębnij energie i oblicz stałe rotacji klastrów zoptymalizowanych pod kątem dużej bazy za pomocą polecenia:
      getRotConsts-dft-lb-ultrafine.csh pw91 N
      gdzie pw91 wskazuje, że użyto funkcjonału gęstości PW91, a N to liczba atomów w klastrze.
   2. Teraz zidentyfikuj unikalne struktury za pomocą polecenia:
      similarityAnalysis.py uf rotConstsData_C
      gdzie UF jest używane jako etykieta nazwy pliku. Masz teraz listę unikalnych konfiguracji zoptymalizowanych na poziomie teorii PW91/6-311++G** zapisaną w pliku uniqueStructures-uf.data.
3. Wykonaj końcową ekstrakcję informacji potrzebnych do obliczenia poprawek termodynamicznych. Skorzystaj z tych informacji, aby obliczyć poprawki termodynamiczne.
   1. Wyodrębnij końcowe energie elektronowe, stałe obrotowe i częstotliwości drgań i użyj ich do obliczenia poprawek termodynamicznych za pomocą polecenia:
      run-thermo-pw91.csh uniqueStructures-uf.data
   2. Skopiuj/wklej dane wyjściowe wiersza poleceń do arkusza "Raw_Energies" arkusza kalkulacyjnego Excel o nazwie "gly-h2o-n.xlsx". Musiałbyś to zrobić dla monomerów (glicyny i wody), a także dla najniższego składnika energetycznego każdego hydratu (gly-h2o-n, gdzie n = 1,2, ...).
   3. Gdy surowe energie są dodawane do pierwszego arkusza arkusza kalkulacyjnego "gly-h2o-n.xlsx", kolejne arkusze "Binding_Energies" i "Hydrate_Distribution" są automatycznie aktualizowane. W szczególności arkusz "Hydrate_Distribution" podaje stężenie równowagowe hydratów w różnych temperaturach (np. 298,15 K), wilgotności względnej (20%, 50%, 100%) i początkowych stężeniach wody ([H₂O]) i glicyny ([Glicyna]). Teoria stojąca za tymi obliczeniami jest opisana w następnym kroku.

6. Obliczanie stężeń atmosferycznych klastrów Gly(H₂O)_n=0-5 w temperaturze pokojowej na poziomie morza

UWAGA: Osiąga się to najpierw kopiując dane termodynamiczne wygenerowane w poprzednim kroku do arkusza kalkulacyjnego i obliczając swobodne energie Gibbsa sekwencyjnego nawadniania. Następnie swobodne energie Gibbsa są wykorzystywane do obliczania stałych równowagi dla każdego sekwencyjnego uwodnienia. Na koniec rozwiązywany jest zestaw równań liniowych, aby uzyskać stężenie równowagowe hydratów dla danego stężenia monomerów, temperatury i ciśnienia.

1. Zacznij od ustanowienia systemu równowagi chemicznej dla sekwencyjnego uwodnienia glicyny, jak pokazano poniżej:
   
2. Oblicz stałe równowagi K_n za pomocą K_n = ^{e-ΔG_n/(k_BT),} gdzie n jest poziomem hydratacji, ΔG_n jest zmianą energii swobodnej Gibbsa w^n-tej reakcji hydratacji, k_B jest stałą Boltzmanna, a T jest temperaturą.
   
3. Ustaw równanie zachowania masy, przyjmując założenie, że suma stężeń równowagowych uwodnionych i nieuwodnionych klastrów glicyny jest równa początkowemu stężeniu wyizolowanej glicyny [Gly]₀. Przepisz ten układ sześciu równoczesnych równań, używając pewnego algebraicznego przegrupowania wyrażeń stałych równowagi, jako
   
4. Rozwiąż układ równań pokazany powyżej, aby uzyskać stężenia równowagowe Gly(H₂O)_{n = 0-5} przy użyciu wartości eksperymentalnej^56,57,58 dla stężenia glicyny w atmosferze, [Gly]₀ = 2,9 x 106 ^cm-3, oraz stężenia wody w atmosferze przy 100% wilgotności względnej i temperaturze 298,15 K⁵⁹, [H₂O] = 7,7 x 10^{17 cm-3}.

Wyniki

Pierwszym zestawem wyników z tego protokołu powinien być zestaw niskoenergetycznych struktur Gly(H₂O)_n=1-5 znalezionych w procedurze konfiguracyjnego próbkowania. Struktury te zostały zoptymalizowane na poziomie teorii PW91/6-311++G** i zakłada się, że są dokładne na potrzeby tego artykułu. Nie ma dowodów sugerujących, że PW91/6-311++G** konsekwentnie zaniża lub przecenia energię wiązania tych gromad. Jego zdolność do przewidywania energii wiązania w stosunku do szacunków MP2/CBS³² i [DLPNO-]CCSD(T)/CBS^60,61 i experiment⁵² wykazuje wiele wahań. To samo dotyczy większości innych funkcjonałów gęstości. Ogólnie rzecz biorąc, każda wartość n = 1 – 5 powinna dać garść struktur niskoenergetycznych w granicach około 5 kcal ^mol-1 struktury o najniższej energii. W tym miejscu skupimy się na pierwszej strukturze utworzonej przez skrypt run-thermo-pw91.csh dla zwięzłości. Rysunek 3 pokazuje najniższe izomery energii elektronowej klastrów Gly(H₂O)_n=0-5. Można zauważyć, że sieć wiązań wodorowych staje się coraz bardziej złożona wraz ze wzrostem liczby cząsteczek wody, a nawet przechodzi z sieci przeważnie płaskiej do trójwymiarowej struktury przypominającej klatkę przy n = 5. W dalszej części tego tekstu posługujemy się energiami i wielkościami termodynamicznymi odpowiadającymi tym pięciu specyficznym gromadom.

Tabela 1 zawiera wielkości termodynamiczne niezbędne do przeprowadzenia protokołu. Tabela 2 przedstawia przykład danych wyjściowych skryptu run-thermo-pw91.csh, w którym drukowane są energie elektronowe, drgania w punkcie zerowym i poprawki termodynamiczne w trzech różnych temperaturach. Dla każdego klastra (rzędu) E[PW91/6-311++G**] odpowiada energiom elektronowym fazy gazowej na poziomie teorii PW91/6-311++G** obliczonym na ultradrobnych siatkach całkowania w jednostkach Hartree, a także energii drgań punktu zerowego (ZPVE) w jednostkach kcal ^mol-1. W każdej temperaturze, 216,65 K, 273,15 K i 298,15 K, wymienione są poprawki termodynamiczne, ∆H entalpia tworzenia w jednostkach kcal ^mol-1, S entropia tworzenia w jednostkach cal ^mol-1, a ∆G energia swobodna Gibbsa w jednostkach kcal ^mol-1. Tabela 3 przedstawia przykładowe obliczenia całkowitej zmiany energii swobodnej Gibbsa podczas hydratacji, a także dla sekwencyjnego nawodnienia. Przykładowe obliczenie całkowitej zmiany energii swobodnej Gibbsa hydratacji dla reakcji

rozpoczyna się od obliczenia energii elektronowej E_PW91 jako

gdzie E_PW91[Gly∙(H₂O)] jest zaczerpnięte z Tabeli 2 kolumna C, a E_PW91[Gly] i E_PW91[H₂O] są pobrane z Tabeli 1 kolumna B. Następnie obliczamy całkowitą zmianę energii fazy gazowej ΔE(0), uwzględniając zmianę energii drgań w punkcie zerowym reakcji jako

aby uzyskać kolumnę D. Tutaj ΔE_{PW91/6−311++G**} jest pobierane z tabeli 3 kolumna C, E_ZPVE[Gly ∙ (H₂O)] z tabeli 2 kolumna D, a E_ZPVE[Gly] i E_ZPVE[H₂O] z tabeli 1 kolumna C. Dla zwięzłości przejdziemy do klastrów temperatury pokojowej, więc pominiemy dane 216.65 K i 273.15 K. W temperaturze pokojowej obliczamy następnie zmianę entalpii reakcji ΔH, korygując zmianę energii fazy gazowej jako

gdzie ΔE(0) jest zaczerpnięte z Tabeli 3 kolumna D, ΔH[Gly∙(H₂O)] jest pobrane z Tabeli 2 kolumna K, a ΔH[Gly] i ΔH[H₂O] są pobrane z Tabeli 1 kolumna J. Na koniec obliczamy zmianę energii swobodnej Gibbsa reakcji ΔG jako

gdzie ΔH jest zaczerpnięte z Tabeli 3 kolumny I, S[Gly∙(H₂O)] jest pobrane z Tabeli 2 kolumny L, a S[Gly] i S[H₂O] są pobrane z Tabeli 1 kolumny K. Zauważ tutaj, że wartości entropii muszą zostać przeliczone na jednostki kcal ^{mol-1 K-1} podczas tego kcal.

Mamy teraz niezbędne ilości do obliczenia atmosferycznych stężeń uwodnionej glicyny, jak pokazano w kroku 6. Wyniki powinny być zbliżone do danych przedstawionych w tabeli 4, ale należy spodziewać się niewielkich różnic liczbowych. Tabela 4 przedstawia stężenia hydratów równowagi stwierdzone w wyniku sformułowania układu sześciu równań w kroku 6.2 w jedno równanie macierzowe i jego późniejsze rozwiązanie. Zaczynamy od uznania faktu, że układ równań można zapisać jako

gdzie K_n jest stałą równowagi dla^n-tego sekwencyjnego uwodnienia glicyny, w jest stężeniem wody w atmosferze, g jest początkowym stężeniem izolowanej glicyny w atmosferze, a g_n jest stężeniem równowagowym Gly(H₂O)_n. Jeśli przepiszemy powyższe równanie jako Ax = b, otrzymamy x = A⁻¹b, gdzie A⁻¹ jest odwrotnością macierzy A. Tę odwrotność można łatwo obliczyć za pomocą wbudowanych funkcji arkusza kalkulacyjnego, jak pokazano w tabeli 4, aby uzyskać ostateczne wyniki.

Rysunek 4 pokazuje stężenie równowagowe uwodnionej glicyny obliczone w Tabeli 4 jako funkcję temperatury przy 100% wilgotności względnej i 1 ciśnienia atmosferycznego. Wynika z niego, że wraz ze spadkiem temperatury z 298,15 K do 216,65 K stężenie glicyny nieuwodnionej (n=0) maleje, a wzrasta glicyny uwodnionej. W szczególności dwuwodzian glicyny (n=2) gwałtownie wzrasta wraz ze spadkiem temperatury, podczas gdy zmiana stężenia innych hydratów jest mniej zauważalna. Ta odwrotna korelacja między temperaturą a stężeniem hydratów jest zgodna z oczekiwaniami, że niższe energie swobodne Gibbsa hydratacji w niższych temperaturach sprzyjają tworzeniu się hydratów.

Rysunek 5 ilustruje zależność wilgotności względnej od stężenia równowagowego hydratów glicyny przy ciśnieniu 298,15 K i 1 atmosferze. Pokazuje to jednoznacznie, że wraz ze wzrostem wilgotności względnej z 20% do 100%, stężenie hydratów (n>0) wzrasta kosztem glicyny nieuwodnionej (n=0). Po raz kolejny bezpośrednia korelacja między wilgotnością względną a stężeniem hydratów jest zgodna z ideą, że obecność większej liczby cząsteczek wody przy wyższej wilgotności względnej sprzyja tworzeniu się hydratów.

Jak przedstawiono, ten protokół daje jakościowe zrozumienie uwodnionych populacji glicyny w atmosferze. Zakładając początkowe stężenie wyizolowanej glicyny wynoszące 2,9 miliona cząsteczek na centymetr sześcienny, widzimy, że glicyna nieuwodniona (n = 0) jest najliczniejszym gatunkiem w większości warunków z wyjątkiem T = 216,65 K i RH = 100%. Dihydrat (n=2), który ma najniższą sekwencyjną energię swobodną Gibbsa hydratacji we wszystkich trzech temperaturach, jest najobfitszym hydratem w rozważanych tutaj warunkach. Przewiduje się, że monohydrat (n=1) i większe hydraty (n≥3) zostaną znalezione w znikomych ilościach. Po przyjrzeniu się klasie Rysunek 3, obfitość klastrów n = 1–4 można powiązać ze stabilnością i odkształceniem w sieci wiązań wodorowych klastrów. Klastry te mają cząsteczki wody związane z wodorem związanym z ugrupowaniem kwasu karboksylowego glicyny w geometrii bardzo przypominającej różne struktury pierścieniowe związane z wodorem, co czyni je szczególnie stabilnymi.

Rysunek 1: Schematyczny opis bieżącej procedury. Duża pula struktur przypuszczenia generowanych przez algorytm genetyczny (GA) jest udoskonalana przez serię optymalizacji geometrii PW91, aż do uzyskania zestawu zbieżnych struktur. Częstotliwości drgań tych struktur są obliczane i wykorzystywane do obliczania swobodnej energii powstawania Gibbsa, która z kolei jest wykorzystywana do obliczania stężeń równowagowych gromad w warunkach otoczenia. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Rysunek 2: Reprezentatywna struktura katalogów dla każdego klastra. Wewnętrzne skrypty zawarte w tym protokole wymagają struktury katalogów pokazanej powyżej, gdzie n jest liczbą cząsteczek wody. Dla każdego n w gly-h2o-n istnieją następujące podkatalogi: GA dla algorytmu genetycznego z katalogiem GA/pm7, QM dla mechaniki kwantowej z QM/pw91-sb dla PW91/6-31+G*, QM/pw91-lb dla PW91/6-311++G** oraz QM/pw91-lb/ultrafine do optymalizacji i końcowych obliczeń wibracyjnych na ultradrobnych siatkach całkowania. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Rysunek 3: Reprezentatywne struktury niskoenergetyczne Gly(_H2O)_n=0-5. Klastry te były globalnymi minimami energii elektronowej zoptymalizowanymi na poziomie teorii PW91/6-311++G**. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Rysunek 4: Zależność Gly(_H2O)_n=0-5 od temperatury przy wilgotności względnej 100% i ciśnieniu 1 atm. Stężenie hydratów podawane jest w jednostkach cząsteczek ^cm-3. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Rysunek 5: Zależność wilgotności względnej Gly(_H2O)_n=0-5 przy ciśnieniu 298,15 K i 1 atm. Stężenie hydratów podawane jest w jednostkach cząsteczek ^cm-3. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

. . . Papierów TGL pkt. zł TGL zł pkt. pkt. pkt. pkt. pkt. pkt.

	E[PW91/6-311++G**]		216,65 tys			273.15 tys			298,15 tys
	LB-UF	ZPVE (Organizacja Wyzwolenia	∆H	S	∆G	∆H	S	∆G	∆H	S	∆G
Woda	-76.430500	Godzina 13.04	1,72	Z dnia 42,59	Klasa 5,54	Informacja o tym, że 2,17	Z dnia 44,44	Pytanie 3,08	Godzina 2,37	Kategoria 45,14	1,96
glicyna	-284.434838	Z dnia 48,55	Godzina 2,65	69,53	Rejon 36,14	Godzina 3,70	73,81	Rozdział 32.09	Rozdział 4,22	75,61	Dnia 30,22

Tabela 1: Energie monomerów. Energie elektronowe są w jednostkach Hartree, podczas gdy wszystkie inne wielkości są w jednostkach kcal ^mol-1. Wodę i glicynę zoptymalizowano na poziomie teorii PW91/6-311++G** i obliczono częstotliwości drgań. Poprawki termodynamiczne dla ciśnienia 1 atm i temperatury 298,15 K obliczono za pomocą skryptu thermo.pl.

. . . . Papierów powiedział: pkt. pkt. pkt. pkt. pkt. pkt. powiedział: pkt. pkt. pkt. pkt. powiedział: TGL pkt. pkt. pkt. pkt. powiedział: zł pkt. pkt. pkt. pkt. pkt. pkt. powiedział: pkt. pkt. pkt. pkt. pkt. pkt.

E[PW91/6-311++G**]

0 tys

216,65 tys

273.15 tys

298,15 tys

n

nazwa

LB-UF

ZPVE (Organizacja Wyzwolenia

∆H

S

∆G

∆H

S

∆G

∆H

S

∆G

1

GLY-H2O-1

-360.88481

63,96

Pytanie 3,61

Nr 80,12

Klasa 50,22

Klasa 5.12

86,27

45,52

Klasa 5,85

88,83

Z dnia 43,33

cyfra arabska

GLY-H2O-2

-437.33763

79,33

Godzina 4,53

90,86

64,17

Godzina 6,46

98,78

Klasa 58,81

godz. 7.40

102.06

Godzina 56,30

3

GLY-H2O-3

-513.78620

Nr 94,52

Norma 5,67

Godzina 105.08

Nr kat. 77,42

Godzina 8.08

114,94

Nr 71,19

Godzina 9,23

godz. 119.00

Ocena z dnia 68,27 ust.

4

GLY-H2O-4

-590.23667

Godzina 109,80

6.03

104,98

O godz. 91,30

8,78

116,21

84,40

Dnia 10.11

120,87

81,14

5

GLY-H2O-5

-666.68845

125,80 jeniec

7,26

121,70

106,69

Godzina 10.47

134,83

99,44

Godzina 12.01

140,24

o godz. 96.00

Tabela 2: Energie klastrów. Energie struktur o najniższej energii Gly(H₂O)_n=1-5 znalezione przy użyciu naszej procedury opisanej w Rysunek 1. Energie elektronowe są w jednostkach Hartree, podczas gdy wszystkie inne wielkości są w jednostkach kcal ^mol-1.

powiedział: pkt. jedn. pkt. jedni pkt. jezdnego pkt. pkt. pkt. jedn. pkt. jedni pkt. jezdnego pkt. pkt. powiedział: zł pkt. pkt. pkt. pkt. jezdnego pkt. pkt. pkt. jezdny powiedział: pkt. pkt. pkt. pkt. pkt. pkt. pkt. jedni pkt. pkt. jedni pkt. pkt. jedni jedni powiedział: pkt. pkt. pkt. pkt. pkt. pkt. pkt. pkt. jedni pkt. pkt. pkt. jezdnego pkt. pkt. powiedział: pkt. pkt. pkt. pkt. pkt. ZŁ pkt. jedni pkt. pkt. zł pkt. zł

		Całkowite nawodnienie: Gly + nH2O <-> Gly(H2O)n								Sekwencyjne nawodnienie: Gly(H2O)n-1 + H2O <-> Gly(H2O)n
		E[PW91/6-311++G**]		216.65		273.15		298.15				216.65		273.15		298.15
n	Nazwa systemu	LB-UF	∆E(0)	∆H(T)	∆G(T)	∆H(T)	∆G(T)	∆H(T)	∆G(T)	LB-UF	∆E(0)	∆H(T)	∆G(T)	H(T)	∆G(T)	∆H(T)	∆G(T)
1	GLY-H2O-1	-12,22	-9,85	-10,61	-3,68	-10,61	-1,87	-10,59	-1,07	-12,22	-9,85	-10,61	-3,68	-10,61	-1,87	-10,59	-1,07
cyfra arabska	GLY-H2O-2	-26,22	-21,53	-23,10	-9,27	-23,11	-5,66	-23.09	-4.06	-godz. 14.00	-11,68	-12,49	-5,59	-12.50	-3,79	-12.50	-2,99
3	GLY-H2O-3	-37,56	-30,72	-32,88	-12.90	-32,87	-7,69	-32,82	-5,38	-11,34	-9,19	-9,78	-3,63	-9,76	-2,03	-9,73	-1,32
4	GLY-H2O-4	-50,10	-40,34	-43,48	-15,87	-43,54	-8,71	-43,51	-5,55	-12,54	-9,62	-10.60	-2,97	-10,67	-1,02	-10,69	-0,17
5	GLY-H2O-5	-63,45	-51,41	-55,42	-20,58	-55,51	-11,48	-55,48	-7,45	-13,35	-11.07	-11,94	-4,71	-11,97	-2,77	-11,97	-1,90

Tabela 3: Energie nawodnienia. Całkowita energia hydratacji i energia sekwencyjnego uwodnienia dla Gly(H₂O)_n=1-5 w jednostkach kcal ^mol-1. Tutaj E[PW91/6-311++G**] to zmiana energii elektronowej, ∆E(0) to skorygowana zmiana energii energii w punkcie zerowym (ZPVE), ∆H(T) to zmiana entalpii w temperaturze T, a ∆G(T) to zmiana energii swobodnej Gibbsa w hydratacji każdego klastra Gly(H₂O)_n=1-5.

TGL TGL TGL TGL TGL TGL TGL

	Rozkład hydratu równowagi jako funkcja temperatury i wilgotności względnej
	T=298.15K			T=273.15K			T=216.65K
Gly(H2O)n	RH=100%	RH=50%	RH=20%	RH=100%	RH=50%	RH=20%	RH=100%	RH=50%	RH=20%
0	1.3E+06	2.2E+06	2.7E+06	1.1E+06	2.0E+06	2.7E+06	6.1E+05	1,5E+06	2,5E+06
1	2.3E+05	1.9E+05	9,5E+04	2.0E+05	1.9E+05	9.9E+04	1.2E+05	1,5E+05	9,5E+04
cyfra arabska	1.0E+06	4.3E+05	8.4E+04	1.3E+06	6.1E+05	1.3E+05	1.8E+06	1.1E+06	3.0E+05
3	2.8E+05	5.8E+04	4,5E+03	3.2E+05	7.4E+04	6.3E+03	3.1E+05	9,6E+04	1.0E+04
4	1.1E+04	1.1E+03	3.4E+01	1.3E+04	1.5E+03	5.0E+01	1.1E+04	1.8E+03	7,5E+01
5	7,5E+03	3.9E+02	4.9E+00	1.2E+04	7.2E+02	9.7E+00	2.4E+04	1.9E+03	3.1E+01

Tabela 4: Stężenie hydratu równowagi Gly(H₂O)_n=0-5 jako funkcja temperatura (T=298.15K, 273.15K, 216.65K) i wilgotność względna (RH=100%, 50%, 20%). Stężenie hydratów jest podane w jednostkach cząsteczek ^cm-3 przy założeniu wartości eksperymentalnych^56,57,58, of [Gly]₀ = 2,9 x 10^{6 cm-3} i [H₂O] = 7,7 x 10^{17 cm-3}, 1,6 x 10^{17 cm-3} i 9,9 x 10^{14 cm-3} przy 100% wilgotności względnej i T = 298,15 K, 273,15 K i 216,65 K, odpowiednio⁵⁹.

Pliki uzupełniające. Kliknij tutaj, aby pobrać te pliki.

Dyskusja

Dokładność danych generowanych przez ten protokół zależy głównie od trzech rzeczy: (i) różnorodności konfiguracji próbkowanych w Kroku 2, (ii) dokładności struktury elektronicznej systemu, (iii) oraz dokładności poprawek termodynamicznych. Każdy z tych czynników można rozwiązać, modyfikując metodę poprzez edycję dołączonych skryptów. Pierwszy czynnik można łatwo przezwyciężyć za pomocą większej początkowej puli losowo generowanych struktur, liczniejszych iteracji GA i luźniejszej definicji kryteriów związanych z GA. Ponadto można zastosować inną metodę półempiryczną, taką jak model samospójnego ścisłego wiązania gęstości ładunku (SCC-DFTB)⁶² i model efektywnego potencjału fragmentów (EFP)⁶³ w celu zbadania skutków różnych opisów fizycznych. Głównym ograniczeniem jest tutaj niezdolność metody do tworzenia lub zrywania wiązań kowalencyjnych, co oznacza, że monomery są zamrożone. Procedura GA znajduje tylko najbardziej stabilne względne pozycje tych zamrożonych monomerów zgodnie z opisem półempirycznym.

Dokładność struktury elektronicznej systemu można poprawić na wiele sposobów, z których każdy wiąże się z kosztem obliczeniowym. Można wybrać lepszy funkcjonał gęstości, taki jak M06-2X⁶⁴ i wB97X-V⁶⁵, lub metodę mechaniki kwantowej (QM), taką jak teorie perturbacji Møllera-Plesseta 66,67,68 (MPn) i metody sprzężonego klastra ⁶⁹ (CC) w celu poprawy fizycznego opisu układu. W hierarchii funkcjonałów wydajność generalnie poprawia się po przejściu od funkcjonałów z uogólnionym gradientem aproksymacji (GGA), takich jak PW91, do oddzielonych zakresem funkcjonałów hybrydowych, takich jak wB97X-D i hybrydowych funkcji meta-GGA, takich jak M06-2X.

Wadą metod DFT jest to, że systematyczna zbieżność w kierunku dokładnej wartości nie jest możliwa; jednak metody DFT są niedrogie obliczeniowo i istnieje szeroka gama funkcji dla szerokiej gamy zastosowań.

Energie obliczane przy użyciu metod funkcji falowych, takich jak MP2 i CCSD(T) w połączeniu z korelacją spójnych zbiorów bazowych rosnącej liczby kardynalnej ([aug-]cc-pV[D,T,Q,...]Z) systematycznie zbliżają się do pełnego limitu zestawu podstawowego, ale koszt obliczeniowy każdego obliczenia staje się zaporowy wraz ze wzrostem rozmiaru systemu. Dalsze udoskonalenie struktury elektronicznej można osiągnąć przy użyciu wyraźnie skorelowanych zestawów podstawowych⁷⁰ i ekstrapolacji do limitu pełnego zbioru podstawowego (CBS)⁷¹ . Nasze ostatnie prace sugerują, że wyraźnie skorelowane podejście perturbacyjne drugiego rzędu Møllera-Plesseta (DF-MP2-F12) z dopasowaną gęstością daje energie zbliżone do obliczeń MP2/CBS³². Modyfikacja obecnego protokołu w celu wykorzystania różnych metod struktury elektronicznej obejmuje dwa kroki: (i) przygotowanie pliku wejściowego szablonu zgodnie ze składnią podaną przez oprogramowanie, (ii) i zredagowanie skryptów run-pw91-sb.csh, run-pw91-lb.csh i run-pw91-lb-ultrafine.csh w celu wygenerowania poprawnej składni pliku wejściowego, a także poprawnego skryptu przesyłania dla oprogramowania.

Wreszcie, dokładność poprawek termodynamicznych zależy od metody struktury elektronowej, a także od opisu PES wokół globalnego minimum. Dokładny opis PES wymaga obliczenia pochodnych PES trzeciego i wyższego rzędu w odniesieniu do przemieszczeń w jądrowych stopniach swobody, takich jak kwartowe pole sił^72,73 (QFF), co jest wyjątkowo kosztownym zadaniem. Obecny protokół wykorzystuje przybliżenie oscylatora harmonicznego do częstotliwości drgań, co powoduje konieczność obliczenia tylko do sekundowych pochodnych PES. Takie podejście staje się problematyczne w układach o wysokiej anharmoniczności, takich jak bardzo wiotkie cząsteczki i symetryczne potencjały podwójnej studni ze względu na dużą różnicę między prawdziwym PES a harmonicznym PES. Co więcej, koszt posiadania wysokiej jakości PES z wymagającej obliczeniowo metody struktury elektronicznej tylko potęguje problem kosztów obliczeń częstotliwości drgań. Jednym ze sposobów przezwyciężenia tego problemu jest wykorzystanie energii elektronicznej z wysokiej jakości obliczeń struktury elektronicznej wraz z częstotliwościami drgań obliczonymi na PES o niższej jakości, co skutkuje równowagą między kosztami a dokładnością. Obecny protokół może zostać zmodyfikowany w celu użycia różnych opisów publicznych służb zatrudnienia, jak opisano w poprzednim akapicie; Można jednak również edytować słowa kluczowe częstotliwości drgań w skryptach i szablonach, aby obliczyć anharmoniczne częstotliwości drgań.

Dwie kluczowe kwestie dla każdego konfiguracyjnego protokołu pobierania próbek to wstępna metoda próbkowania powierzchni energii potencjalnej oraz kryteria stosowane do identyfikacji każdego klastra. W naszej poprzedniej pracy szeroko korzystaliśmy z różnych metod. W pierwszym zagadnieniu, wstępnej metodzie próbkowania powierzchni energii potencjalnej, zdecydowaliśmy się na użycie GA z metodami półempirycznymi opartymi na tych czynnikach. Próbkowanie konfiguracyjne przy użyciu intuicji chemicznej²⁶, losowego próbkowania i dynamiki molekularnej (MD)^29,30 nie pozwala na regularne znajdowanie przypuszczalnych globalnych minimów dla klastrów większych niż 10 monomerów, jak zaobserwowaliśmy w naszych badaniach klastrów wodnych¹⁸. Z powodzeniem zastosowaliśmy skakanie po basenie (BH) do badania złożonego PES (_H2O)₁₁⁷⁴, ale wymagało to ręcznego włączenia niektórych potencjalnych izomerów o niskiej energii, których algorytm BH nie znalazł. Porównanie wydajności BH i GA w znajdowaniu globalnego minimum klastrów wody, (H₂O)_n=10-20 wykazało, że GA konsekwentnie znajdowało globalne minimum szybciej niż BH⁷⁵. GA zaimplementowany w OGOLEM i CLUSTER jest bardzo wszechstronny, ponieważ może być stosowany do dowolnego klastra molekularnego i może łączyć się z ogromną liczbą pakietów z klasycznym polem siłowym, półempirycznym, funkcjonałem gęstości i możliwościami ab initio. Wybór PM7 zależy od jego szybkości i rozsądnej dokładności. Praktycznie każda inna metoda półempiryczna wiązałaby się ze znacznie wyższymi kosztami obliczeniowymi.

Jeśli chodzi o drugą kwestię, badaliśmy ją przy użyciu różnych kryteriów, aby zidentyfikować unikalne struktury, począwszy od energii elektronowej, momentów dipolowych, nakładających się RMSD i stałych rotacyjnych. Wykorzystanie momentów dipolowych okazało się trudne, ponieważ obie składowe momentu dipolowego były zależne od orientacji cząsteczki, a całkowity moment dipolowy był bardzo wrażliwy na różnice geometryczne w taki sposób, że trudno było ustalić progi określające, czy struktury są takie same lub unikalne. Kombinacja energii elektronowych i stałych rotacyjnych okazała się najbardziej użyteczna.

Obecne kryteria uznania dwóch struktur za unikalne opierają się na progu różnicy energii wynoszącym 0,10 kcal ^mol-1 i różnicy stałej obrotowej wynoszącej 1%. Dlatego dwie struktury uważa się za różne, jeśli ich energie różnią się o więcej niż 0,10 kcal ^mol-1 (~0,00015 a.u.) I każda z ich trzech stałych obrotowych (A, B, C) różni się o więcej niż 1%. W obszernych wewnętrznych testach porównawczych na przestrzeni lat stwierdzono, że progi te są rozsądnym wyborem. Nasze konfiguracyjne podejście do pobierania próbek i metodologia badań przesiewowych zostały zastosowane do bardzo słabo związanych klastrów, takich jak węglowodory poliaromatyczne skompleksowane z wodą^76,77, a także silnie związane trójskładnikowe hydraty siarczanowe zawierające amoniak i aminy³². W przypadku klastrów, w których należy wziąć pod uwagę różne stany protonacji, najlepszym podejściem jest przeprowadzenie różnych obliczeń GA, z których każde zaczyna się od monomerów w różnych stanach protonacji. Gwarantuje to, że struktury o różnych stanach protonacji są dokładnie przemyślane. Jednak niskopoziomowe obliczenia DFT często pozwalają na zmianę stanów protonacji w trakcie optymalizacji geometrii, uzyskując w ten sposób najbardziej stabilny stan protonacji niezależnie od geometrii początkowej.

Nasze metody próbkowania konfiguracyjnego GA powinny działać dobrze nawet w przypadku cząsteczek dyskietkowych, o ile kody GA są połączone z ogólnymi, niesparametryzowanymi metodami, które pozwalają monomerom na przyjęcie różnych konfiguracji w trakcie przebiegu GA. Na przykład, połączenie GA z PM7 pozwoliłoby na zmianę struktur monomerów, ale jeśli ich wiązania zostaną zerwane, co miałoby miejsce, gdy zmieniają się stany protonacji, struktury mogą zostać odrzucone jako nieakceptowalni kandydaci.

Rozważaliśmy różne sposoby korygowania niedociągnięć w przybliżeniu harmonicznych, zwłaszcza tych wynikających z niskich częstotliwości drgań. Włączenie aproksymacji quasi-harmonicznej do obecnej metodologii nie jest trudne. Nadal jednak pojawiają się pytania dotyczące metody quasi-harmonicznej, zwłaszcza jeśli chodzi o częstotliwość graniczną, poniżej której będzie ona stosowana. Nie ma również rygorystycznych prac porównawczych oceniających wiarygodność aproksymacji quasi-RRHO, mimo że konwencjonalna mądrość sugeruje, że powinna to być poprawa w stosunku do aproksymacji RRHO.

Przedstawiony w ten sposób protokół można uogólnić na dowolny układ niekowalencyjnie związanych klastrów molekularnych fazy gazowej. Można go również uogólnić na użycie dowolnej metody półempirycznej, metody i oprogramowania do struktury elektronicznej oraz metody i oprogramowania do analizy wibracyjnej poprzez edycję skryptów i szablonów. Zakłada się, że użytkownik czuje się komfortowo z interfejsem wiersza poleceń systemu Linux, skryptami języka Python i obliczeniami o wysokiej wydajności. Nieznana składnia i wygląd systemu operacyjnego Linux oraz brak doświadczenia w pisaniu skryptów to największa pułapka w tym protokole i to właśnie w niej nowi studenci mają największe problemy. Protokół ten jest od lat z powodzeniem stosowany w różnych implementacjach w naszej grupie, koncentrując się głównie na wpływie kwasu siarkowego i amoniaku na tworzenie aerozoli. Dalsze ulepszenia tego protokołu będą obejmowały bardziej niezawodny interfejs do większej liczby oprogramowania do budowy struktur elektronowych, alternatywne implementacje algorytmu genetycznego i być może wykorzystanie nowszych metod do szybszych obliczeń energii elektronowej i wibracyjnej. Nasze obecne zastosowania tego protokołu badają znaczenie aminokwasów we wczesnych stadiach tworzenia aerozoli w obecnej atmosferze oraz w tworzeniu większych cząsteczek biologicznych w środowiskach prebiotycznych.

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
Avogadro	https://avogadro.cc		Program do wizualizacji molekularnej o otwartym kodzie źródłowym
Gaussian [09/16] Oprogramowanie	http://www.gaussian.com/		Komercyjny program do struktury elektronicznej ab initio
MOPAC 2016	http://openmopac.net/MOPAC2016.html		Półempiryczny program o otwartym kodzie źródłowym
OGOLEM Software	https://www.ogolem.org		Globalny program optymalizacyjny oparty na algorytmach genetycznych
OpenBabel	http://openbabel.org/wiki/Main_Page		Biblioteka chemiczno-komputerowa typu open source
calcRotConsts.py	Shields Group, Wydział Chemii, Uniwersytet		Skrypt Python do obliczania stałych obrotowych
calcSymmetry.csh	Grupa Shields, Wydział Chemii, Uniwersytet		FurmanaSkrypt powłoki do obliczania liczby symetrii cząsteczki o podanych współrzędnych kartezjańskich
combine-GA.csh	Grupa Shields, Wydział Chemii, Uniwersytet		Skrypt powłoki do łączenia stałych energetycznych i rotacyjnych z różnych katalogów GA
combine-QM.csh	Grupa Shields, Wydział Chemii, Uniwersytet		Skrypt powłoki do łączenia stałych energetycznych i rotacyjnych z różnych katalogów
gaussianE.csh	Shields Group, Wydział Chemii, Uniwersytet		Skrypt powłoki do ekstrakcji energii Gaussian 09
gaussianFreqs.csh	Shields Group, Wydział Chemii, Uniwersytet		Skrypt powłoki do wyodrębniania częstotliwości drgań Gaussa 09
getrotconsts	Grupa Shields, Wydział Chemii, Uniwersytet		FurmanaPlik wykonywalny do obliczania stałych rotacyjnych przy danych współrzędnych kartezjańskich cząsteczki
getRotConsts-dft-lb.csh	Shields Group, Wydział Chemii, Uniwersytet		Skrypt powłoki do obliczania stałych rotacyjnych dla partii struktur zoptymalizowanych pod kątem DFT o dużej podstawie
getRotConsts-dft-lb-ultrafine.csh	Grupa Shields, Wydział Chemii, Uniwersytet		FurmanaSkrypt powłoki do obliczania stałych obrotowych dla partii ultradrobnych struktur zoptymalizowanych pod kątem DFT
getRotConsts-dft-sb.csh	Grupa Shields, Wydział Chemii, Uniwersytet		Skrypt powłoki do obliczania stałych obrotowych dla partii struktur zoptymalizowanych pod kątem DFT o małej podstawie
getRotConsts-GA.csh	Grupa Shields, Wydział Chemii, Uniwersytet		FurmanaSkrypt powłoki do obliczania stałych rotacyjnych dla partii struktur zoptymalizowanych pod kątem algorytmu genetycznego
global-minimum-coords.xyz	Grupa Shields, Wydział Chemii, Uniwersytet		Współrzędne kartezjańskie globalnych struktur minimalnych gly-(h2o)n, gdzie n=0-5
make-thermo-gaussian.csh	Grupa Shields, Wydział Chemii, Uniwersytet		FurmanaSkrypt powłoki do wyodrębniania danych z gaussowskich plików wyjściowych i tworzenia plików wejściowych dla skryptu thermo.pl
ogolem-input-file.ogo	Grupa Shields, Wydział Chemii, Uniwersytet		Ogolem przykładowy plik wejściowy
ogolem-submit-script.pbs	Shields Group, Wydział Chemii, Furman University		PBS plik wsadowy do obliczeń Ogolem
README.docx	Shields Group, Wydział Chemii, Uniwersytet		Wyjaśnienia pomagające czytelnikom efektywnie korzystać ze skryptów
runogolem.csh	Shields Group, Wydział Chemii, Uniwersytet		Skrypt powłoki do uruchamiania OGOLEM
run-pw91-lb.csh	Shields Group, Wydział Chemii, Uniwersytet		Skrypt powłoki do uruchamiania partii obliczeń optymalizacyjnych DFT o
dużej podstawierun-pw91-lb-ultrafine.csh	Grupa Shields, Wydział Chemii, Uniwersytet		Skrypt powłoki do uruchamiania partii ultradrobnych obliczeń optymalizacyjnych DFT
run-pw91-sb.csh	Grupa Shields, Wydział Chemii, Uniwersytet		Skrypt powłoki do uruchamiania partii obliczeń optymalizacyjnych DFT o małej podstawie
run-thermo-pw91.csh	Grupa Shields, Wydział Chemii, Uniwersytet		FurmanaSkrypt powłoki do obliczania poprawek termodynamicznych dla partii struktur zoptymalizowanych pod kątem DFT
similarityAnalysis.py	Shields Group, Wydział Chemii, Uniwersytet		Skrypt Pythona do określania unikalnych struktur na podstawie stałych rotacyjnych i
symetrii	Grupa Shields, Wydział Chemii, Uniwersytet		Plik wykonywalny do obliczania symetrii molekularnej o współrzędnych
symmetry.c	(C) 1996, 2003 S. Patchkovskii, Serguei.Patchkovskii@sympatico.ca		Kod C do określania symmstrii molekularnej cząsteczki o danych współrzędnych kartezjańskich
template-marcy.pbs	Grupa Shields, Wydział Chemii, Uniwersytet		Szablon skryptu przesyłania PBS, który wykorzystuje OGOLEM
template-pw91.com	Shields Group, Wydział Chemii, Uniwersytet		FurmanaSzablon wejściowy Gaussian 09
template-pw91-HL.com	Shields Group, Wydział Chemii, Uniwersytet		Szablon wejściowy Gaussian 09 do ultradrobnej optymalizacji DFT
thermo.pl https://www.nist.gov/mml/csd/chemical-informatics-research-group/products-and-services/program-computing-ideal-gas			Skrypt open source Perl do obliczania poprawek
termodynamicznych gazu doskonałegogly-h2o-n.xlsx	Shields Group, Wydział Chemii, Uniwersytet		Arkusz kalkulacyjny Excel dla pełnego protokołu
table-1.xlsx	Grupa Shields, Wydział Chemii, Uniwersytet		Arkusz kalkulacyjny Excel
table-2.xlsx	Grupa Shields, Wydział Chemii, Uniwersytet Furmana		Arkusz kalkulacyjny Excel
table-3.xlsx	Grupa Shields, Wydział Chemia, Arkusz	kalkulacyjny Excel
table-4.xlsx	Grupa Shields, Wydział Chemii,	Arkusz kalkulacyjny Excel
water.xyz	Grupa Shields, Wydział Chemii, Uniwersytet		Współrzędne kartezjańskie wody
glycine.xyz	Grupa Shields, Wydział Chemii, Uniwersytet		Współrzędne kartezjańskie glicyna