Artykuł metodologiczny

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

9.4K wyświetleń

⸱

DOI:

10.3791/60964

⸱

8 kwietnia 2020

 ,  ,  ,  , 
1Department of Chemistry, Furman University, 2College of Charleston

W tym artykule

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 atmosfery1. 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 mniejszej2; Jednak ten reżim rozmiarów jest trudny do uzyskania przy użyciu technik eksperymentalnych3,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 prenukleacji5,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)16. 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ą26, losowe pobieranie próbek27,28, dynamika molekularna (MD)29,30, skakanie po basenie (BH)31, 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 Jakuba33. 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.36 opracowali systematyczną metodę, w której cząsteczki wody są dodawane w punktach sfer Fibonacciego37 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)38 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 ABCluster39. Kubecka et al.40 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-xTB41. 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 OGOLEM25, która została z powodzeniem zastosowana do różnych globalnych problemów z optymalizacją i konfiguracyjnym próbkowaniem42,43,44,45, aby wygenerować początkowy zestaw punktów próbkowania. PES będzie opisany przez model PM7 46 zaimplementowany w programie MOPAC201647. 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 PW9148 funkcjonału korelacji wymiany uogólnionego gradientu (GGA) w implementacji DFT Gaussian0949 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 atmosferycznych51,52.

Ten protokół zakłada, że użytkownik ma dostęp do klastra obliczeń o wysokiej wydajności (HPC) z przenośnym systemem wsadowym53 (PBS), MOPAC2016 (http://openmopac.net/MOPAC2016.html)47, OGOLEM (https://www.ogolem.org)25, Gaussian 09 (https://gaussian.com)49, oraz OpenBabel54 (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. Wyznaczanie struktury o minimalnej energii dla izolowanej glicyny i wody

UWAGA: Celem jest tutaj realizacja dwóch zadań: (i) uzyskanie struktur o minimalnej energii dla odizolowanych cząsteczek wody i glicyny do wykorzystania w próbkowaniu konfiguracyjnym algorytmu genetycznego oraz (ii) obliczenie poprawek termodynamicznych do energii tych cząsteczek w fazie gazowej do wykorzystania w obliczeniach stężeń atmosferycznych.

  1. Na komputerze lokalnym otwórz nową sesję programu Avogadro.
    1. Kliknij Build > Insert > Peptide i w oknie Insert Peptide wybierz Gly, aby wygenerować monomer glicyny w oknie wizualizacji.
    2. Kliknij Extensions > Gaussian i edytuj pierwszą linię w polu tekstowym tak, aby brzmiała: '# pw91pw91/6-311++G** int(Acc2E=12,UltraFine) scf(conver=12) opt(tight,maxcyc=300) freq'. Kliknij Generate i zapisz plik wejściowy jako glycine.com.
    3. Należy pamiętać, że jeśli cząsteczka wykazuje znaczną elastyczność konformacyjną, co ma miejsce w przypadku glicyny55, kluczowe jest przeprowadzenie analizy konformacyjnej w celu zidentyfikowania struktury minimum globalnego oraz innych niskoenergetycznych konformerów. OpenBabel54 oferuje zaawansowane narzędzia do przeszukiwania konformacyjnego, wykorzystujące różne algorytmy i szybkie pola siłowe. Choć konformery mogą ulegać relaksacji i interkonwersji podczas GA oraz późniejszych obliczeń, czasami konieczne jest przeprowadzenie wielu obliczeń GA, z których każde rozpoczyna się od innego konformera.
  2. Na komputerze lokalnym otwórz nową sesję w programie Avogadro.
    1. Kliknij Build > Insert > Fragment i w oknie Insert Fragment wyszukaj frazę "water", aby uzyskać współrzędne wody.
    2. Kliknij Extensions > Gaussian i edytuj pierwszą linię w polu tekstowym tak, aby brzmiała: '# pw91pw91/6-311++G** int(Acc2E=12,UltraFine) scf(conver=12) opt(tight,maxcyc=300) freq'. Kliknij Generate i zapisz plik wejściowy jako water.com.
  3. Prześlij oba pliki .com na zdalny klaster. Po zalogowaniu się do klastra zdalnego wywołaj program Gaussian 09 w skrypcie wysyłającym zadanie w trybie wsadowym, aby rozpocząć obliczenia. Po zakończeniu obliczeń wyodrębnij współrzędne kartezjańskie (pliki .xyz) struktur o minimalnej energii, korzystając z programu OpenBabel. W przypadku glicyny należy wykonać następującą komendę:
    obabel -ig09 glycine.log -oxyz > glycine.xyz
    Te dwa pliki .xyz zostaną wykorzystane do próbkowania konfiguracyjnego GA w następnym kroku.

2. Próbkowanie konfiguracyjne klastrów Gly(H2O)n=1-5 w oparciu o algorytm genetyczny

UWAGA: Celem jest tutaj uzyskanie zestawu struktur o niskiej energii dla Gly(H2O)n=1-5 na niedrogim półempirycznym poziomie teorii, przy użyciu modelu PM746 zaimplementowanego w programie MOPAC47. Jest niezbędne, aby katalog roboczy miał dokładnie taką organizację i strukturę, jak pokazano na Rysunku 2. Ma to zapewnić poprawne działanie dostosowanej powłoki oraz skryptów Python.

  1. Skopiuj wszystkie niezbędne skrypty na zdalny klaster i dodaj ich lokalizację do $PATH
    1. Umieść wszystkie skrypty i pliki szablonów w folderze (np. scripts) i skopiuj go na zdalny klaster
    2. Upewnij się, że wszystkie skrypty mają uprawnienia do wykonywania
    3. Dodaj lokalizację katalogu ze skryptami do zmiennej środowiskowej $PATH, wpisując następujące polecenia w terminalu. Domyślna lokalizacja skryptów jest ustawiona na $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 ścieżki dostępu
      1. Powszechna powłoka Bash:
        export SCRIPTS_HOME=/path/to/scripts
        export PATH=${SCRIPTS_HOME}:${PATH}
      2. Powłoka Tcsh/Csh:
        setenv SCRIPTS_HOME /path/to/scripts
        setenv PATH ${SCRIPTS_HOME}:${PATH}
  2. Na zdalnym klastrze skonfiguruj i uruchom obliczenia GA:
    1. Utwórz katalog o nazwie gly-h2o-n, gdzie n oznacza liczbę cząsteczek wody.
    2. Wewnątrz katalogu gly-h2o-n utwórz podkatalog o nazwie GA w celu uruchomienia obliczeń algorytmu genetycznego.
    3. Skopiuj pliki wejściowe OGOLEM (np. pm7.ogo), współrzędne kartezjańskie monomerów (np. glycine.xyz, water.xyz) oraz skrypt przesyłania zadań PBS (np. run.pbs) do katalogu GA.
    4. Wprowadź niezbędne zmiany w pliku wejściowym OGOLEM oraz pliku przesyłania zadań.
    5. Prześlij obliczenia do realizacji. Po uruchomieniu obliczeń program OGOLEM utworzy w katalogu GA nowy katalog o nazwie odpowiadającej przedrostkowi pliku wejściowego OGOLEM (np. pm7) i zapisze w nim nowo wygenerowane współrzędne.
  3. Po zakończeniu obliczeń zestaw energie i stałe rotacyjne, a następnie wykorzystaj te informacje do określenia unikalnych struktur o niskiej energii:
    1. Zmień katalog na gly-h2o-n/GA/pm7 a następnie
    2. Wyodrębnij energie i oblicz stałe rotacyjne klastrów zoptymalizowanych metodą GA za pomocą polecenia:
      getRotConsts-GA.csh N 0 99
      gdzie N to liczba atomów w klastrze molekularnym, a „0 99” oznacza, że rozmiar puli GA wynosi 100, z indeksami od 0 do 99. Wygeneruje to plik o nazwie rotConstsData_C, który zawiera posortowaną listę wszystkich konfiguracji klastrów zoptymalizowanych metodą GA, ich energie oraz stałe rotacyjne.
    3. Wykonaj polecenie:
      similarityAnalysis.py pm7 rotConstsData_C
      gdzie pm7 zostanie użyte jako etykieta nazewnictwa plików, aby znaleźć i zapisać unikalne klastry zoptymalizowane metodą GA. Wygeneruje to plik o nazwie uniqueStructures-pm7.data, który zawiera posortowaną listę unikalnych konfiguracji zoptymalizowanych metodą GA. Jest to lista unikalnych struktur lokalnych minimów dla klastra Gly(H2O)n zoptymalizowanego na poziomie teorii PM7; struktury te są teraz gotowe do doprecyzowania przy użyciu DFT.
  4. Przejdź do katalogu gly-h2o-n/GA i połącz wyniki z wielu porównywalnych uruchomień GA za pomocą skryptu combine-GA.csh. Składnia to:
    combine-GA.csh <etykieta> <lista katalogów z uruchomieniami GA>
    W tym konkretnym przypadku polecenie:
    combine-GA.csh pm7 pm7
    wygeneruje nową listę unikalnych struktur o nazwie 'uniqueStructures-pm7.data' w katalogu gly-h2o-n/GA.

3. Optymalizacja z wykorzystaniem metody QM z małym zestawem funkcji bazowych

UWAGA: Celem tego etapu jest udoskonalenie próbkowania konfiguracyjnego klastrów Gly(H2O)n=1-5 przy użyciu lepszego opisu kwantowo-mechanicznego, aby uzyskać mniejszy, ale dokładniejszy zestaw struktur klastrów Gly(H2O)n=1-5. Strukturami wyjściowymi dla tego kroku są wyniki z kroku 2.

  1. Przygotuj i przeprowadź obliczenia DFT z małym zestawem funkcji bazowych:
    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. Przejdź do katalogu gly-h2o-n/QM/pw91-sb.
    4. Uruchom skrypt próbkowania konfiguracyjnego DFT z małym zestawem funkcji bazowych za pomocą polecenia:
      run-pw91-sb.csh uniqueStructures-pm7.data sb QUEUE 10
      gdzie sb jest etykietą dla tego zestawu obliczeń, QUEUE to preferowana kolejka w klastrze obliczeniowym, a 10 oznacza, że 10 obliczeń zostanie zgrupowanych w jednym zadaniu wsadowym. Skrypt ten automatycznie wygeneruje pliki wejściowe dla Gaussian 09 i prześle wszystkie obliczenia. Wprowadź 'test' w miejscu 'QUEUE', aby wykonać uruchomienie próbne.
  2. Po zakończeniu przesłanych obliczeń wyodrębnij i przeanalizuj wyniki.
    1. Wyodrębnij energie i oblicz stałe rotacyjne dla klastrów zoptymalizowanych z użyciem małego zestawu funkcji bazowych za pomocą polecenia:
      getRotConsts-dft-sb.csh pw91 N
      gdzie pw91 wskazuje na użycie funkcjonału gęstości PW91, a N to liczba atomów w klastrze. Operacja ta utworzy plik o nazwie rotConstsData_C.
    2. Następnie zidentyfikuj unikalne struktury za pomocą polecenia:
      similarityAnalysis.py sb rotConstsData_C
      gdzie sb służy jako etykieta nazewnictwa plików. Lista unikalnych konfiguracji zoptymalizowanych na poziomie teorii PW91/6-31+G* zostanie zapisana w pliku uniqueStructures-sb.data.
  3. Przejdź do katalogu gly-h2o-n/QM i połącz wyniki z wielu porównywalnych uruchomień QM za pomocą skryptu combine-QM.csh. Składnia brzmi:
    combine-QM.csh <etykieta> <lista katalogów z obliczeniami QM>
    W tym konkretnym przypadku polecenie:
    combine-QM.csh sb pw91-sb
    wygeneruje nową listę unikalnych struktur o nazwie 'uniqueStructures-sb.data' w katalogu gly-h2o-n/QM.

4. Dalsze doprecyzowanie z wykorzystaniem metody QM z dużym zestawem funkcji bazowych

UWAGA: Celem tego etapu jest dalsze udoskonalenie próbkowania konfiguracyjnego klastrów Gly(H2O)n=1-5 przy użyciu lepszego opisu mechaniki kwantowej. Strukturami wyjściowymi dla tego kroku są wyniki z Kroku 3.

  1. Prześlij bardziej wiarygodne obliczenia przy użyciu większego bazy funkcji.
    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 z dużą bazą funkcji za pomocą polecenia:
      run-pw91-lb.csh uniqueStructures-sb.data lb QUEUE 10
      gdzie lb jest etykietą dla tego zestawu obliczeń, QUEUE to preferowana kolejka w klastrze obliczeniowym, a 10 oznacza, że 10 obliczeń ma zostać zgrupowanych w jedno zadanie wsadowe. Skrypt ten automatycznie wygeneruje pliki wejściowe dla Gaussian 09 i prześle wszystkie obliczenia. Aby przeprowadzić próbne uruchomienie testowe, wpisz 'test' w miejscu 'QUEUE'.
  2. Po zakończeniu przesłanych obliczeń wyodrębnij i przeanalizuj dane
    1. Oblicz stałe rotacyjne klastrów zoptymalizowanych z dużą bazą funkcji 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. Następnie zidentyfikuj unikalne struktury za pomocą polecenia:
      similarityAnalysis.py lb rotConstsData_C
      gdzie lb jest używane jako etykieta nazewnictwa plików. Otrzymasz teraz listę unikalnych konfiguracji zoptymalizowanych na poziomie teorii PW91/6-311++G**, zapisaną w pliku uniqueStructures-lb.data.

5. Obliczenia końcowej energii i poprawek termodynamicznych

UWAGA: Celem jest tutaj uzyskanie struktury wibracyjnej i energii klastrów Gly(H2O)n=1-5 przy użyciu dużego zbioru funkcji bazowych i ultrawąskiej siatki całkowania, aby obliczyć wymagane poprawki termochemiczne.

  1. Wykorzystując wyniki z poprzedniego kroku, przeprowadź bardziej precyzyjne 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. Uruchom skrypt DFT dla ultrafine z dużą bazą funkcji za pomocą polecenia:
      run-pw91-lb-ultrafine.csh uniqueStructures-lb.data uf QUEUE 10
      gdzie uf jest etykietą dla tego zestawu obliczeń, QUEUE to preferowana kolejka w klastrze obliczeniowym, a 10 oznacza, że 10 obliczeń zostanie zgrupowanych w jedno zadanie wsadowe. Skrypt ten automatycznie wygeneruje pliki wejściowe dla programu Gaussian 09 i uruchomi wszystkie obliczenia. Aby przeprowadzić test próbny, wpisz 'test' w miejscu 'QUEUE'.
  2. Po zakończeniu obliczeń wyodrębnij i przeanalizuj dane
    1. Wyodrębnij energie i oblicz stałe rotacyjne dla klastrów zoptymalizowanych z użyciem dużej bazy funkcji za pomocą polecenia:
      getRotConsts-dft-lb-ultrafine.csh pw91 N
      gdzie pw91 wskazuje na użycie funkcjonału gęstości PW91, a N to liczba atomów w klastrze.
    2. Zidentyfikuj teraz unikalne struktury za pomocą polecenia:
      similarityAnalysis.py uf rotConstsData_C
      gdzie uf służy jako etykieta nazewnictwa plików. Otrzymasz listę unikalnych konfiguracji zoptymalizowanych na poziomie teorii PW91/6-311++G**, zapisaną w pliku uniqueStructures-uf.data.
  3. Przeprowadź końcowe wyodrębnienie informacji niezbędnych do obliczenia poprawek termodynamicznych. Wykorzystaj te informacje do wyznaczenia poprawek termodynamicznych.
    1. Wyodrębnij końcowe energie elektronowe, stałe rotacyjne oraz częstotliwości wibracyjne i wykorzystaj je do obliczenia poprawek termodynamicznych za pomocą polecenia:
      run-thermo-pw91.csh uniqueStructures-uf.data
    2. Skopiuj i wklej wynik z wiersza poleceń do arkusza 'Raw_Energies' w pliku Excel o nazwie 'gly-h2o-n.xlsx'. Należy to zrobić dla monomerów (glicyny i wody), a także dla struktury o najniższej energii z każdego hydratu (gly-h2o-n, gdzie n=1,2, …).
    3. W miarę dodawania surowych energii do pierwszego arkusza pliku 'gly-h2o-n.xlsx', kolejne arkusze 'Binding_Energies' oraz 'Hydrate_Distribution' będą aktualizowane automatycznie. W szczególności arkusz 'Hydrate_Distribution' dostarcza równowagowe stężenia hydratów w różnych temperaturach (np. 298.15K), przy różnej wilgotności względnej (20%, 50%, 100%) oraz początkowych stężeniach wody ([H2O]) i glicyny ([Glycine]). Teoria stojąca za tymi obliczeniami jest opisana w następnym kroku.

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

UWAGA: Osiąga się to poprzez skopiowanie danych termodynamicznych wygenerowanych w poprzednim kroku do arkusza kalkulacyjnego i obliczenie energii swobodnej Gibbsa dla sekwencyjnej hydratacji. Następnie energie swobodnej Gibbsa są wykorzystywane do obliczenia stałych równowagi dla każdego etapu sekwencyjnej hydratacji. Na koniec rozwiązuje się układ równań liniowych, aby uzyskać równowagowe stężenie hydratów dla określonego stężenia monomerów, temperatury i ciśnienia.

  1. Zacznij od ustalenia układu równowag chemicznych dla sekwencyjnej hydratacji glicyny, jak pokazano poniżej:
    figure-protocol-1
  2. Oblicz stałe równowagi Kn korzystając ze wzoru Kn = e-ΔGn/(kBT), gdzie n oznacza stopień hydratacji, ΔGn to zmiana energii swobodnej Gibbsa dla ntej reakcji hydratacji, kB to stała Boltzmanna, a T to temperatura.
    figure-protocol-2
  3. Sformułuj równanie zachowania masy, przyjmując, że suma równowagowych stężeń klastrów glicyny hydratowanej i niehydratowanej jest równa stężeniu początkowemu izolowanej glicyny [Gly]0. Zapisz ten układ sześciu równań jednoczesnych, wykorzystując przekształcenia algebraiczne wyrażeń na stałą równowagi, w postaci:
    figure-protocol-3
  4. Rozwiąż powyższy układ równań, aby wyznaczyć stężenia równowagowe Gly(H2O)n = 0-5, wykorzystując wartość eksperymentalną56,57,58 dla stężenia glicyny w atmosferze [Gly]0 = 2.9 x 106 cm-3 oraz stężenie wody w atmosferze przy 100% wilgotności względnej i temperaturze 298.15 K59, [H2O] = 7.7 x 1017 cm-3.

Wyniki

Pierwszy zestaw wyników z niniejszego protokołu powinien stanowić zbiór struktur o niskiej energii kompleksów Gly(H2O)n=1-5 wyznaczonych w procedurze próbkowania konfiguracyjnego. Struktury te zostały zoptymalizowane na poziomie teorii PW91/6-311++G** i przyjmuje się, że są one wystarczająco dokładne na potrzeby niniejszej pracy. Nie ma dowodów sugerujących, aby metoda PW91/6-311++G** systematycznie zaniżała lub zawyżała energię wiązania tych klastrów. Jej zdolność do przewidywania energii wiązania w odniesieniu do szacunków MP2/CBS32 oraz [DLPNO-]CCSD(T)/CBS60,61 i wyników eksperymentalnych52 wykazuje duże fluktuacje. To samo dotyczy większości innych funkcjonałów gęstości. Generalnie dla każdej wartości n = 1 – 5 powinno zostać otrzymanych kilka struktur o niskiej energii, mieszczących się w granicach około 5 kcal mol-1 od struktury o najniższej energii. Dla zwięzłości skupiamy się tutaj na pierwszej strukturze wygenerowanej przez skrypt run-thermo-pw91.csh. Rysunek 3 przedstawia izomery o najniższej energii elektronowej dla klastrów Gly(H2O)n=0-5. Można zauważyć, że sieć wiązań wodorowych staje się bardziej złożona wraz ze wzrostem liczby cząsteczek wody i przechodzi od sieci głównie planarnej do trójwymiarowej struktury klatkowej przy n = 5. W dalszej części tekstu wykorzystano energie i wielkości termodynamiczne odpowiadające tym pięciu konkretnym klastrom.

Tabela 1 zawiera wielkości termodynamiczne niezbędne do przeprowadzenia protokołu. Tabela 2 przedstawia przykład wyniku działania skryptu run-thermo-pw91.csh, w którym wydrukowane są energie elektronowe, poprawki energii punktu zerowego drgań oraz poprawki termodynamiczne w trzech różnych temperaturach. Dla każdego klastra (wiersza), E[PW91/6-311++G**] odpowiada energiom elektronowym w fazie gazowej na poziomie teorii PW91/6-311++G** obliczonym na ultrafine siatkach całkowania w jednostkach Hartree, a także energii wibracyjnej punktu zerowego (ZPVE) w jednostkach kcal mol-1. Dla każdej temperatury: 216.65 K, 273.15 K oraz 298.15 K, wymieniono poprawki termodynamiczne: ∆H – entalpię tworzenia w jednostkach kcal mol-1, S – entropię tworzenia w jednostkach cal mol-1 oraz ∆G – energię swobodną Gibbsa tworzenia w jednostkach kcal mol-1. Tabela 3 przedstawia przykładowe obliczenie całkowitej zmiany energii swobodnej Gibbsa hydratacji, a także dla hydratacji sekwencyjnej. Przykładowe obliczenie całkowitej zmiany energii swobodnej Gibbsa hydratacji dla reakcji

figure-results-1

zaczyna się od obliczenia energii elektronicznej EPW91 jako

figure-results-2

gdzie EPW91[Gly∙(H2O)] pochodzi z kolumny C Tabeli 2, a EPW91[Gly] oraz EPW91[H2O] pochodzą z kolumny B Tabeli 1. Następnie obliczamy całkowitą zmianę energii w fazie gazowej ΔE(0), uwzględniając zmianę energii drgań zerowych reakcji jako

figure-results-3

aby otrzymać kolumnę D. Tutaj ΔEPW91/6−311++G** pochodzi z Tabeli 3, kolumna C, EZPVE[Gly ∙ (H2O)] z Tabeli 2, kolumna D, a EZPVE[Gly] oraz EZPVE[H2O] z Tabeli 1, kolumna C. Dla zwięzłości przejdziemy do klastrów w temperaturze pokojowej, pomijając dane dla 216,65 K i 273,15 K. W temperaturze pokojowej obliczamy następnie zmianę entalpii reakcji ΔH, korygując zmianę energii w fazie gazowej w następujący sposób:

figure-results-4

gdzie ΔE(0) pochodzi z kolumny D Tabeli 3, ΔH[Gly∙(H2O)] pochodzi z kolumny K Tabeli 2, a ΔH[Gly] oraz ΔH[H2O] pochodzą z kolumny J Tabeli 1. Na koniec obliczamy zmianę energii swobodnej Gibbsa reakcji ΔG w następujący sposób:

figure-results-5

gdzie ΔH pochodzi z kolumny I Tabeli 3, S[Gly∙(H2O)] pochodzi z kolumny L Tabeli 2, a S[Gly] i S[H2O] pochodzą z kolumny K Tabeli 1. Należy tutaj zauważyć, że na tym etapie wartości entropii muszą zostać przeliczone na jednostki kcal mol-1 K-1.

Mamy teraz niezbędne ilości do obliczenia stężeń atmosferycznych uwodnionej glicyny, zgodnie z opisem w Kroku 6. Wyniki powinny być zbliżone do danych przedstawionych w Tabeli 4, jednak należy spodziewać się niewielkich różnic liczbowych. Tabela 4 przedstawia równowagowe stężenia uwodnionych form otrzymane z przekształcenia układu sześciu równań z Kroku 6.2 w jedno równanie macierzowe i jego późniejszego rozwiązania. Zaczynamy od przyjęcia założenia, że układ równań można zapisać jako

figure-results-6

gdzie Kn jest stałą równowagi dla n-tej sekwencyjnej hydratacji glicyny, w to stężenie wody w atmosferze, g to początkowe stężenie wyizolowanej glicyny w atmosferze, a gn to stężenie równowagowe Gly(H2O)n. Jeśli zapiszemy powyższe równanie jako Ax = b, otrzymamy x = A−1b, gdzie A−1 jest macierzą odwrotną do macierzy A. Macierz odwrotną można łatwo obliczyć za pomocą wbudowanych funkcji arkusza kalkulacyjnego, jak pokazano w Tabeli 4, aby uzyskać końcowe wyniki.

Rysunek 4 przedstawia równowagowe stężenie uwodnionej glicyny obliczone w Tabeli 4 w funkcji temperatury przy 100% wilgotności względnej i ciśnieniu 1 atmosfery. Widać na nim, że wraz ze spadkiem temperatury od 298,15K do 216,65K stężenie nieuwodnionej glicyny (n=0) maleje, natomiast stężenia uwodnionej glicyny rosną. W szczególności stężenie dwuwodnika glicyny (n=2) gwałtownie wzrasta wraz ze spadkiem temperatury, podczas gdy zmiana stężeń pozostałych hydratów jest mniej zauważalna. Ta odwrotna korelacja między temperaturą a stężeniem hydratów jest zgodna z przewidywaniem, że niższe wartości energii swobodnej Gibbsa uwodnienia w niższych temperaturach sprzyjają tworzeniu się hydratów.

Rysunek 5 przedstawia zależność stężenia równowagowego hydratów glicyny od wilgotności względnej w temperaturze 298.15K i pod ciśnieniem 1 atmosfery. Wyraźnie pokazuje on, że wraz ze wzrostem RH od 20% do 100% stężenie hydratów (n>0) wzrasta kosztem niehydratowanej glicyny (n=0). Ponownie bezpośrednia korelacja między wilgotnością względną a stężeniem hydratów jest zgodna z tezą, że obecność większej liczby cząsteczek wody przy wyższej RH sprzyja tworzeniu się hydratów.

W przedstawiony sposób protokół ten umożliwia jakościowe zrozumienie populacji uwodnionego glicyny w atmosferze. Przy założeniu początkowego stężenia izolowanej glicyny wynoszącego 2,9 mln cząsteczek na centymetr sześcienny, zauważamy, że nieuwodniona glicyna (n=0) jest najliczniejszym gatunkiem w większości warunków, z wyjątkiem T=216.65K i RH=100%. Dihydrat (n=2), który wykazuje najniższą sekwencyjną energię swobodną Gibbsa uwodnienia w trzech badanych temperaturach, jest najliczniejszym hydratem w rozważanych tutaj warunkach. Przewiduje się, że monohydrat (n=1) i większe hydraty (n≥3) występują w ilościach pomijalnych. Po analizie Rysunku 3 można stwierdzić, że liczebność klastrów n = 1–4 jest powiązana ze stabilnością i naprężeniem w sieci wiązań wodorowych klastrów. W klastrach tych cząsteczki wody są związane wiązaniami wodorowymi z grupą kwasu karboksylowego glicyny w geometrii ściśle przypominającej różne pierścieniowe struktury wiązań wodorowych, co czyni je szczególnie stabilnymi.

figure-results-7
Rysunek 1: Schematyczny opis aktualnej procedury. Duża pula struktur szacunkowych generowanych przez algorytm genetyczny (GA) jest udoskonalana poprzez serię optymalizacji geometrii PW91, aż do uzyskania zestawu struktur zbieżnych. Obliczane są częstotliwości drgań tych struktur, które następnie służą do wyznaczenia energii swobodnej Gibbsa tworzenia, wykorzystywanej z kolei do obliczenia stężeń równowagowych klastrów w warunkach otoczenia. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

figure-results-8
Rycina 2: Reprezentatywna struktura katalogów dla każdego klastra. Skrypty wewnętrzne zawarte w niniejszym protokole wymagają struktury katalogów przedstawionej powyżej, gdzie n oznacza 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 dla optymalizacji i końcowych obliczeń wibracyjnych na siatkach integracyjnych typu ultrafine. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

figure-results-9
Rysunek 3: Reprezentatywne struktury o niskiej energii dla Gly(H2O)n=0-5. Klastry te stanowią globalne minima energii elektronowej zoptymalizowane na poziomie teorii PW91/6-311++G**. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

figure-results-10
Rycina 4: Zależność od temperatury Gly(H2O)n=0-5 przy 100% wilgotności względnej i ciśnieniu 1 atm. Stężenie hydratów podano w jednostkach molecules cm-3. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

figure-results-11
Rysunek 5: Zależność Gly(H2O)n=0-5 od wilgotności względnej przy temperaturze 298.15 K i ciśnieniu 1 atm. Stężenie hydratów podano w jednostkach molecules cm-3. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

E[PW91/6-311++G**]216,65 K273,15 K298,15 K
LB-UFZPVE∆HS∆G∆HS∆G∆HS∆G
woda-76.43050013.041.7242.595.542.1744.443.082.3745.141.96
glicyna-284.43483848.552.6569.5336.143.7073.8132.094.2275.6130.22

Tabela 1: Energie monomerów. Energie elektronowe podano w jednostkach Hartree, natomiast wszystkie pozostałe wielkości w jednostkach kcal mol-1. Woda i glicyna zostały zoptymalizowane na poziomie teorii PW91/6-311++G**, a następnie obliczono ich częstotliwości wibracyjne. Korekty termodynamiczne dla ciśnienia 1 atm i temperatury 298.15 K obliczono przy użyciu skryptu thermo.pl.

E[PW91/6-311++G**]0 K216.65 K273.15 K298.15 K
nnazwaLB-UFZPVE∆HS∆G∆HS∆G∆HS∆G
1gly-h2o-1-360.8848163.963.6180.1250.225.1286.2745.525.8588.8343.33
2gly-h2o-2-437.3376379.334.5390.8664.176.4698.7858.817.40102.0656.30
3gly-h2o-3-513.7862094.525.67105.0877.428.08114.9471.199.23119.0068.27
4gly-h2o-4-590.23667109.806.03104.9891.308.78116.2184.4010.11120.8781.14
5gly-h2o-5-666.68845125.807.26121.70106.6910.47134.8399.4412.01140.2496.00

Tabela 2: Energie klastrów. Energie struktur Gly(H2O)n=1-5 o najniższej energii, znalezione przy użyciu procedury opisanej na Rysunku 1. Energie elektronowe podano w jednostkach Hartree, natomiast wszystkie pozostałe wartości w jednostkach kcal mol-1.

Całkowita hydratacja: Gly + nH2O <-> Gly(H2O)nHydratacja sekwencyjna: Gly(H2O)n-1 + H2O <-> Gly(H2O)n
E[PW91/6-311++G**]216.65273.15298.15216.65273.15298.15
nnazwa układuLB-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)
1gly-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
2gly-h2o-2-26.22-21.53-23.10-9.27-23.11-5.66-23.09-4.06-14.00-11.68-12.49-5.59-12.50-3.79-12.50-2.99
3gly-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
4gly-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
5gly-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 hydratacji. Całkowita energia hydratacji oraz energia hydratacji sekwencyjnej dla Gly(H2O)n=1-5 w jednostkach kcal mol-1. Tutaj E[PW91/6-311++G**] oznacza zmianę energii elektronowej, ∆E(0) to zmiana energii skorygowana o energię drgań w punkcie zerowym (ZPVE), ∆H(T) to zmiana entalpii w temperaturze T, a ∆G(T) to zmiana energii swobodnej Gibbsa podczas hydratacji każdego klastra Gly(H2O)n=1-5.

Równowagowy rozkład hydratów w funkcji temperatury i wilgotności względnej 
T=298,15KT=273,15 KT=216,65K
Gly(H2O)nRH=100%RH=50%wilgotność względna = 20%RH=100%RH=50%RH=20%RH=100%RH=50%RH=20%
01,3E+062,2E+062,7E+061,1E+062,0E+062,7E+066,1E+051,5E+062,5E+06
12,3E+051,9E+059,5E+042,0E+051,9E+059,9E+041,2E+051,5E+059,5E+04
21,0E+064,3E+058,4E+041,3E+066,1E+051,3E+051,8E+061,1E+063,0E+05
32,8E+055,8E+044,5E+033,2E+057,4E+046,3E+033,1E+059,6E+041,0E+04
41,1E+041,1E+033,4E+011,3E+041,5E+035,0E+011,1E+041,8E+037,5E+01
57,5E+033,9E+024,9E+001,2E+047,2E+029,7E+002,4E+041,9E+033,1E+01

Tabela 4: Równowagowe stężenia hydratów Gly(H2O)n=0-5 w funkcji temperatury (T=298,15K, 273,15K, 216,65K) oraz wilgotności względnej (RH=100%, 50%, 20%). Stężenie hydratów podano w jednostkach cząsteczek cm-3, przyjmując wartości eksperymentalne56,57,58, dla [Gly]0 = 2,9 x 106 cm-3 oraz [H2O] = 7,7 x 1017 cm-3, 1,6 x 1017 cm-3 i 9,9 x 1014 cm-3 przy 100% wilgotności względnej i T = odpowiednio 298,15 K, 273,15 K oraz 216,65 K59.

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)62 i model efektywnego potencjału fragmentów (EFP)63 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-2X64 i wB97X-V65, lub metodę mechaniki kwantowej (QM), taką jak teorie perturbacji Møllera-Plesseta 66,67,68 (MPn) i metody sprzężonego klastra 69 (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 podstawowych70 i ekstrapolacji do limitu pełnego zbioru podstawowego (CBS)71 . 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/CBS32. 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 chemicznej26, 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 wodnych18. Z powodzeniem zastosowaliśmy skakanie po basenie (BH) do badania złożonego PES (H2O)1174, 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, (H2O)n=10-20 wykazało, że GA konsekwentnie znajdowało globalne minimum szybciej niż BH75. 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 aminy32. 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.

Podziękowania

Ten projekt był wspierany przez granty CHE-1229354, CHE-1662030, CHE-1721511 i CHE-1903871 z National Science Foundation (GCS), Arnold and Mabel Beckman Foundation Beckman Scholar Award (AGG) oraz Barry M. Goldwater Scholarship (AGG). Wykorzystano zasoby obliczeniowe o wysokiej wydajności należące do konsorcjum MERCURY (http://www.mercuryconsortium.org).

Materiały

Lista materiałów użytych w tym artykule
NazwaFirmaNumer katalogowyKomentarze
Avogadrohttps://avogadro.ccProgram do wizualizacji molekularnej o otwartym kodzie źródłowym
Gaussian [09/16] Oprogramowaniehttp://www.gaussian.com/Komercyjny program do struktury elektronicznej ab initio
MOPAC 2016http://openmopac.net/MOPAC2016.htmlPółempiryczny program o otwartym kodzie źródłowym
OGOLEM Softwarehttps://www.ogolem.orgGlobalny program optymalizacyjny oparty na algorytmach genetycznych
OpenBabelhttp://openbabel.org/wiki/Main_PageBiblioteka chemiczno-komputerowa typu open source
calcRotConsts.pyShields Group, Wydział Chemii, UniwersytetSkrypt Python do obliczania stałych obrotowych
calcSymmetry.cshGrupa Shields, Wydział Chemii, UniwersytetFurmanaSkrypt powłoki do obliczania liczby symetrii cząsteczki o podanych współrzędnych kartezjańskich
combine-GA.cshGrupa Shields, Wydział Chemii, UniwersytetSkrypt powłoki do łączenia stałych energetycznych i rotacyjnych z różnych katalogów GA
combine-QM.cshGrupa Shields, Wydział Chemii, UniwersytetSkrypt powłoki do łączenia stałych energetycznych i rotacyjnych z różnych katalogów
gaussianE.cshShields Group, Wydział Chemii, UniwersytetSkrypt powłoki do ekstrakcji energii Gaussian 09
gaussianFreqs.cshShields Group, Wydział Chemii, UniwersytetSkrypt powłoki do wyodrębniania częstotliwości drgań Gaussa 09
getrotconstsGrupa Shields, Wydział Chemii, UniwersytetFurmanaPlik wykonywalny do obliczania stałych rotacyjnych przy danych współrzędnych kartezjańskich cząsteczki
getRotConsts-dft-lb.cshShields Group, Wydział Chemii, UniwersytetSkrypt powłoki do obliczania stałych rotacyjnych dla partii struktur zoptymalizowanych pod kątem DFT o dużej podstawie
getRotConsts-dft-lb-ultrafine.cshGrupa Shields, Wydział Chemii, UniwersytetFurmanaSkrypt powłoki do obliczania stałych obrotowych dla partii ultradrobnych struktur zoptymalizowanych pod kątem DFT
getRotConsts-dft-sb.cshGrupa Shields, Wydział Chemii, UniwersytetSkrypt powłoki do obliczania stałych obrotowych dla partii struktur zoptymalizowanych pod kątem DFT o małej podstawie
getRotConsts-GA.cshGrupa Shields, Wydział Chemii, UniwersytetFurmanaSkrypt powłoki do obliczania stałych rotacyjnych dla partii struktur zoptymalizowanych pod kątem algorytmu genetycznego
global-minimum-coords.xyzGrupa Shields, Wydział Chemii, UniwersytetWspółrzędne kartezjańskie globalnych struktur minimalnych gly-(h2o)n, gdzie n=0-5
make-thermo-gaussian.cshGrupa Shields, Wydział Chemii, UniwersytetFurmanaSkrypt 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.ogoGrupa Shields, Wydział Chemii, UniwersytetOgolem przykładowy plik wejściowy
ogolem-submit-script.pbsShields Group, Wydział Chemii, Furman UniversityPBS plik wsadowy do obliczeń Ogolem
README.docxShields Group, Wydział Chemii, UniwersytetWyjaśnienia pomagające czytelnikom efektywnie korzystać ze skryptów
runogolem.cshShields Group, Wydział Chemii, UniwersytetSkrypt powłoki do uruchamiania OGOLEM
run-pw91-lb.cshShields Group, Wydział Chemii, UniwersytetSkrypt powłoki do uruchamiania partii obliczeń optymalizacyjnych DFT o
dużej podstawierun-pw91-lb-ultrafine.cshGrupa Shields, Wydział Chemii, UniwersytetSkrypt powłoki do uruchamiania partii ultradrobnych obliczeń optymalizacyjnych DFT
run-pw91-sb.cshGrupa Shields, Wydział Chemii, UniwersytetSkrypt powłoki do uruchamiania partii obliczeń optymalizacyjnych DFT o małej podstawie
run-thermo-pw91.cshGrupa Shields, Wydział Chemii, UniwersytetFurmanaSkrypt powłoki do obliczania poprawek termodynamicznych dla partii struktur zoptymalizowanych pod kątem DFT
similarityAnalysis.pyShields Group, Wydział Chemii, UniwersytetSkrypt Pythona do określania unikalnych struktur na podstawie stałych rotacyjnych i
symetriiGrupa Shields, Wydział Chemii, UniwersytetPlik wykonywalny do obliczania symetrii molekularnej o współrzędnych
symmetry.c(C) 1996, 2003 S. Patchkovskii, Serguei.Patchkovskii@sympatico.caKod C do określania symmstrii molekularnej cząsteczki o danych współrzędnych kartezjańskich
template-marcy.pbsGrupa Shields, Wydział Chemii, UniwersytetSzablon skryptu przesyłania PBS, który wykorzystuje OGOLEM
template-pw91.comShields Group, Wydział Chemii, UniwersytetFurmanaSzablon wejściowy Gaussian 09
template-pw91-HL.comShields Group, Wydział Chemii, UniwersytetSzablon 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-gasSkrypt open source Perl do obliczania poprawek
termodynamicznych gazu doskonałegogly-h2o-n.xlsxShields Group, Wydział Chemii, UniwersytetArkusz kalkulacyjny Excel dla pełnego protokołu
table-1.xlsxGrupa Shields, Wydział Chemii, UniwersytetArkusz kalkulacyjny Excel
table-2.xlsxGrupa Shields, Wydział Chemii, Uniwersytet FurmanaArkusz kalkulacyjny Excel
table-3.xlsxGrupa Shields, Wydział Chemia, Arkuszkalkulacyjny Excel
table-4.xlsxGrupa Shields, Wydział Chemii,Arkusz kalkulacyjny Excel
water.xyzGrupa Shields, Wydział Chemii, UniwersytetWspółrzędne kartezjańskie wody
glycine.xyzGrupa Shields, Wydział Chemii, UniwersytetWspółrzędne kartezjańskie glicyna
Furmana Furmana Furmana QM Furmana Furmana Furmana Furmana Furmana Furmana Furmana Furmana Furmana Furmana Furmana Furmana energii Furmana kartezjańskich Furmana Furmana Furmana Furmana Uniwersytetu FurmanaUniwersytetu Furmana Furmana Furmana

Bibliografia

  1. Foster, P., Ramaswamy, V. Climate Change 2007 The Scientific Basis. Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., Averyt, K. B., Tignor, M., Miller, H. L. , Cambridge University Press. Cambridge, U.K. (2007).
  2. Kulmala, M., et al. Toward direct measurement of atmospheric nucleation. Science. 318 (5847), 89-92 (2007).
  3. Sipila, M., et al. The role of sulfuric acid in atmospheric nucleation. Science. 327 (5970), 1243-1246 (2010).
  4. Jiang, J., et al. First measurement of neutral atmospheric cluster and 1 - 2 nm particle number size distributions during nucleation events. Aerosol Science and Technology. 45 (4), (2011).
  5. Dunn, M. E., Pokon, E. K., Shields, G. C. Thermodynamics of forming water clusters at various Temperatures and Pressures by Gaussian-2, Gaussian-3, Complete Basis Set-QB3, and Complete Basis Set-APNO model chemistries; implications for atmospheric chemistry. Journal of the American Chemical Society. 126 (8), 2647-2653 (2004).
  6. Pickard, F. C., Pokon, E. K., Liptak, M. D., Shields, G. C. Comparison of CBSQB3, CBSAPNO, G2, and G3 thermochemical predictions with experiment for formation of ionic clusters of hydronium and hydroxide ions complexed with water. Journal of Chemical Physics. 122, 024302(2005).
  7. Pickard, F. C., Dunn, M. E., Shields, G. C. Comparison of Model Chemistry and Density Functional Theory Thermochemical Predictions with Experiment for Formation of Ionic Clusters of the Ammonium Cation Complexed with Water and Ammonia; Atmospheric Implications. Journal of Physical Chemistry A. 109 (22), 4905-4910 (2005).
  8. Alongi, K. S., Dibble, T. S., Shields, G. C., Kirschner, K. N. Exploration of the Potential Energy Surfaces, Prediction of Atmospheric Concentrations, and Vibrational Spectra of the HO2•••(H2O)n (n=1-2) Hydrogen Bonded Complexes. Journal of Physical Chemistry A. 110 (10), 3686-3691 (2006).
  9. Allodi, M. A., Dunn, M. E., Livada, J., Kirschner, K. N. Do Hydroxyl Radical-Water Clusters, OH(H2O)n, n=1-5, Exist in the Atmosphere. Journal of Physical Chemistry A. 110 (49), 13283-13289 (2006).
  10. Kirschner, K. N., Hartt, G. M., Evans, T. M., Shields, G. C. In Search of CS2(H2O)n=1-4 Clusters. Journal of Chemical Physics. 126, 154320(2007).
  11. Hartt, G. M., Kirschner, K. N., Shields, G. C. Hydration of OCS with One to Four Water Molecules in Atmospheric and Laboratory Conditions. Journal of Physical Chemistry A. 112 (19), 4490-4495 (2008).
  12. Morrell, T. E., Shields, G. C. Atmospheric Implications for Formation of Clusters of Ammonium and 110 Water Molecules. Journal of Physical Chemistry A. 114 (12), 4266-4271 (2010).
  13. Temelso, B., et al. Quantum Mechanical Study of Sulfuric Acid Hydration: Atmospheric Implications. Journal of Physical Chemistry A. 116 (9), 2209(2012).
  14. Husar, D. E., Temelso, B., Ashworth, A. L., Shields, G. C. Hydration of the Bisulfate Ion: Atmospheric Implications. Journal of Physical Chemistry A. 116 (21), 5151-5163 (2012).
  15. Bustos, D. J., Temelso, B., Shields, G. C. Hydration of the Sulfuric Acid – Methylamine Complex and Implications for Aerosol Formation. Journal of Physical Chemistry A. 118 (35), 7430-7441 (2014).
  16. Wales, D. J., Scheraga, H. A. Global optimization of clusters, crystals, and biomolecules. Science. 27 (5432), 1368-1372 (1999).
  17. Day, M. B., Kirschner, K. N., Shields, G. C. Global search for minimum energy (H2O)n clusters, n = 3 - 5. The Journal of Physical Chemistry A. 109 (30), 6773-6778 (2005).
  18. Shields, R. M., Temelso, B., Archer, K. A., Morrell, T. E., Shields, G. C. Accurate predictions of water cluster formation, (H2O)n=2-10. The Journal of Physical Chemistry A. 114 (43), 11725-11737 (2010).
  19. Temelso, B., Archer, K. A., Shields, G. C. Benchmark structures and binding energies of small water clusters with anharmonicity corrections. The Journal of Physical Chemistry A. 115 (43), 12034-12046 (2011).
  20. Temelso, B., Shields, G. C. The role of anharmonicity in hydrogen-bonded systems: The case of water clusters. The Journal of Chemical Theory and Computation. 7 (9), 2804-2817 (2011).
  21. Von Freyberg, B., Braun, W. Efficient search for all low energy conformations of polypeptides by Monte Carlo methods. The Journal of Computational Chemistry. 12 (9), 1065-1076 (1991).
  22. Rakshit, A., Yamaguchi, T., Asada, T., Bandyopadhyay, P. Understanding the structure and hydrogen bonding network of (H2O)32 and (H2O)33: An improved Monte Carlo temperature basin paving (MCTBP) method of quantum theory of atoms in molecules (QTAIM) analysis. RSC Advances. 7 (30), 18401-18417 (2017).
  23. Deaven, D. M., Ho, K. M. Molecular geometry optimization with a genetic algorithm. Physical Review Letters. 75, 288-291 (1995).
  24. Hartke, B. Application of evolutionary algorithms to global cluster geometry optimization. Applications of Evolutionary Computation in Chemistry. , Springer. Berlin. (2004).
  25. Dieterich, J. M., Hartke, B. OGOLEM: Global cluster structure optimization for arbitrary mixtures of flexible molecules. A multiscaling, object-oriented approach. Molecular Physics. 108 (3-4), 279-291 (2010).
  26. Herb, J., Nadykto, A. B., Yu, F. Large ternary hydrogen-bonded pre-nucleation clusters in the Earth's atmosphere. Chemical Physics Letters. 518, 7-14 (2011).
  27. Ortega, I. K., et al. From quantum chemical formation free energies to evaporation rates. Atmospheric Chemistry and Physics. 12 (1), 225-235 (2012).
  28. Elm, J., Bilde, M., Mikkelsen, K. V. Influence of Nucleation Precursors on the Reaction Kinetics of Methanol with the OH Radical. Journal of Physical Chemistry A. 117 (30), 6695-6701 (2013).
  29. Loukonen, V., et al. Enhancing effect of dimethylamine in sulfuric acid nucleation in the presence of water - a computational study. Atmospheric Chemistry and Physics. 10 (10), 4961-4974 (2010).
  30. Temelso, B., Phan, T. N., Shields, G. C. Computational study of the hydration of sulfuric acid dimers: implications for acid dissociation and aerosol formation. Journal of Physical Chemistry A. 116 (39), 9745-9758 (2012).
  31. Jiang, S., et al. Study of Cl-(H2O)n (n = 1-4) using basin-hopping method coupled with density functional theory. Journal of Computational Chemistry. 35 (2), 159-165 (2014).
  32. Temelso, B., et al. Effect of mixing ammonia and alkylamines on sulfate aerosol formation. Journal of Physical Chemistry A. 122 (6), 1612-1622 (2018).
  33. Perdew, J. P., Ruzsinszky, A., Tao, J. Prescription for the design and selection of density functional approximations: More constraint satisfaction with fewer fits. Journal of Chemical Physics. 123, 062201(2005).
  34. Riplinger, C., Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. Journal of Chemical Physics. 138, 034106(2013).
  35. Riplinger, C., Pinski, P., Becker, U., Valeev, E. F., Neese, F. Sparse maps--A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory. Journal of Chemical Physics. 144 (2), 024109(2016).
  36. Kildgaard, J. V., Mikkelsen, K. V., Bilde, M., Elm, J. Hydration of atmospheric molecular clusters: a new method for systematic configurational sampling. Journal of Physical Chemistry A. 122 (22), 5026-5036 (2018).
  37. González, Á Measurement of areas on a sphere Using Fibonacci and latitude-longitude lattices. Mathematical Geosciences. 42, 49-64 (2010).
  38. Karaboga, D., Basturk, B. On the performance of artificial bee colony (ABC) algorithm. Applied Soft Computing. 8 (1), 687-697 (2008).
  39. Zhang, J., Doig, M. Global optimization of rigid molecules using the artificial bee colony algorithm. Physical Chemistry Chemical Physics. 18 (4), 3003-3010 (2016).
  40. Kubecka, J., Besel, V., Kurten, T., Myllys, N., Vehkamaki, H. Configurational sampling of noncovalent (atmospheric) molecular clusters: sulfuric acid and guanidine. Journal of Physical Chemistry A. 123 (28), 6022-6033 (2019).
  41. Grimme, S., Bannwarth, C., Shushkov, P. A Robust and accurate tight-binding quantum chemical method for structures, vibrational frequencies, and noncovalent Interactions of large molecular systems parametrized for all spd-block elements (Z = 1-86). Journal of Chemical Theory and Computation. 13 (5), 1989-2009 (2017).
  42. Buck, U., Pradzynski, C. C., Zeuch, T., Dieterich, J. M., Hartke, B. A size resolved investigation of large water clusters. Physical Chemistry Chemical Physics. 16 (15), 6859(2014).
  43. Forck, R. M., et al. Structural diversity in sodium doped water trimers. Physical Chemistry Chemical Physics. 14 (25), 9054-9057 (2012).
  44. Witt, C., Dieterich, J. M., Hartke, B. Cluster structures influenced by interaction with a surface. Physical Chemistry Chemical Physics. 20 (23), 15661-15670 (2018).
  45. Freitbert, A., Dieterich, J. M., Hartke, B. Exploring self-organization of molecular tether molecules on a gold surface by global structure optimization. The Journal of Computational Chemistry. 40 (22), 1978-1989 (2019).
  46. Stewart, J. J. P. Optimization of parameters for semiempirical methods VI: More modifications to the NDDO approximations and re-optimization of parameters. The Journal of Molecular Modeling. 19 (1), 1-32 (2013).
  47. Stewart, J. J. P. MOPAC2012 Computational Chemistry. , Available from: http://openmopac.net (2012).
  48. Burke, K., Perdew, J. P., Wang, Y. Derivation of a generalized gradient approximation: The PW91 density functional. Electronic Density Functional Theory. , Springer. Boston, MA. 81-111 (1998).
  49. Frisch, M. J., et al. Gaussian 09, Revision A.02. , Gaussian, Inc. Wallingford, CT. (2016).
  50. Ditchfield, R., Hehre, W. J., Pople, J. A. Self-consistent molecular-orbital methods. IX. An extended Gaussian-type basis for molecular-orbital studies of organic molecules. The Journal of Chemical Physics. 54 (2), 724(1971).
  51. Elm, J., Bilde, M., Mikkelsen, K. V. Assessment of density functional theory in predicting structures and free energies of reaction of atmospheric prenucleation clusters. The Journal of Chemical Theory and Computation. 8 (6), 2071-2077 (2012).
  52. Elm, J., Mikkelsen, K. V. Computational approaches for efficiently modelling of small atmospheric clusters. Chemical Physics Letters. 615, 26-29 (2014).
  53. Bayucan, A., et al. PBS Portable Batch System. , MRJ Technology Solutions. Mountain View, CA. (1999).
  54. O'Boyle, N. M., et al. Open Babel: An open chemical toolbox. Journal of Cheminformatics. 3, 33(2011).
  55. Csaszar, A. G. Conformers of gaseous glycine. Journal of the American Chemical Society. 114 (24), 9568-9575 (1992).
  56. Zhang, Q., Anastasio, C. Free and combined amino compounds in atmospheric fine particles (PM2.5) and fog waters from Northern California. Atmospheric Environment. 37 (16), 2247-2258 (2003).
  57. Matsumoto, K., Uematsu, M. Free amino acids in marine aerosols over the western North Pacific Ocean. Atmospheric Environment. 39 (11), 2163-2170 (2005).
  58. Mandalakis, M., Apostolaki, M., Stephanou, E. G. Trace analysis of free and combined amino acids in atmospheric aerosols by gas chromatography-mass spectrometry. Journal of Chromatography A. 1217 (1), 143-150 (2010).
  59. Seinfeld, J. H., Pandis, S. N. Atmospheric Chemistry and Physics, 3rd Ed. , John Wiley & Sons. Hoboken, N.J. (2016).
  60. Myllys, N., Elm, J., Halonen, R., Kurten, T., Vehkamaki, H. Coupled cluster evaluation of atmospheric acid-base clusters with up to 10 molecules. The Journal of Physical Chemistry A. 120 (4), 621-630 (2016).
  61. Elm, J., Bilde, M., Mikkelsen, K. V. Assessment of binding energies of atmospherically relevant clusters. Physical Chemistry Chemical Physics. 15 (39), (2013).
  62. Elstner, M. The SCC-DFTB method and its application to biological systems. Theoretical Chemistry Accounts. 116 (1-3), 316-325 (2006).
  63. Kaliman, I. A., Slipchenko, L. V. LIBEFP: A new parallel implementation of the effective fragment potential method as a portable software library. The Journal of Computational Chemistry. 34 (26), 2284-2292 (2013).
  64. Zhao, Y., Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and trasition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theoretical Chemistry Accounts. 120 (1-3), 215-241 (2008).
  65. Mardirossian, N., Head-Gordon, M. wB97X-V: A 10-parameter, range-separated hybrid, generalized gradient approximation density functional with nonlocal correlation, designed by a survival-of-the-fittest strategy. Physical Chemistry Chemical Physics. 16 (21), 9904-9924 (2014).
  66. Head-Gordon, M., Pople, J. A., Frisch, M. J. MP2 energy evaluation by direct methods. Chemical Physics Letters. 153 (6), 503-506 (1988).
  67. Pople, J. A., Seeger, R., Krishnan, R. Variational configuration interaction methods and comparison with perturbation theory. The International Journal of Quantum Chemistry. 12, 149-163 (1977).
  68. Pople, J. A., Binkley, J. S., Seeger, R. Theoretical models incorporating electron correlation. The International Journal of Quantum Chemistry. 10 (10), 1-19 (1976).
  69. Monkhorst, H. J. Calculation of properties with the coupled-cluster method. The International Journal of Quantum Chemistry. 12 (11), 421-432 (1977).
  70. Klopper, W., Manby, F. R., Ten-No, S., Valeev, E. F. R12 methods in explicitly correlated molecular electronic structure theory. International Reviews in Physical Chemistry. 25, 427-468 (2006).
  71. Hattig, C. Optimization of auxiliary basis sets for RI-MP2 and RI-CC2 calculations: Core-valence and quintuple-z basis sets for H to Ar and QZVPP basis sets for Li to Kr. Physical Chemistry Chemical Physics. 7 (1), 59-66 (2005).
  72. Barone, V. Anharmonic vibrational properties by a fully automated second-order perturbative approach. The Journal of Chemical Physics. 122, 014108(2005).
  73. Barone, V. Vibrational zero-point energies and thermodynamic functions beyond the harmonic approximation. The Journal of Chemical Physics. 120 (7), 3059-3065 (2004).
  74. Temelso, B., et al. Exploring the Rich Potential Energy Surface of (H2O)11 and Its Physical Implications. Journal of Chemical Theory and Computation. 14 (2), 1141-1153 (2018).
  75. Kabrede, H., Hentschke, R. Global minima of water clusters (H2O)N, N≤25, described by three empirical potentials. Journal of Physical Chemistry B. 107 (16), (2003).
  76. Steber, A. L., et al. Capturing the Elusive Water Trimer from the Stepwise Growth of Water on the Surface of a Polycyclic Aromatic Hydrocarbon Acenaphthene. Journal of Physical Chemistry Letters. 8 (23), 5744-5750 (2017).
  77. Perez, C., et al. Corrannulene and its complex with water: A tiny cup of water. Physical Chemistry Chemical Physics. 19 (22), 14214-14223 (2017).

Przedruki i uprawnienia

Tagi

Algorytm genetycznypr bkowanie konfiguracyjnemetody p empirycznemetody mechaniki kwantowejpoprawki termodynamiczneenergia swobodna Gibbsauwodnione klastry glicyny
Obliczanie stężeń atmosferycznych klastrów molekularnych na podstawie termochemii ab initio
Skontaktuj się z nami
Poleć bibliotece

Prywatność

Warunki korzystania z usługi

Zasady

Prawa autorskie © 2026 MyJoVE Corporation. Wszelkie prawa zastrzeżone.