Analiza stopów i płynów z symulacji dynamiki molekularnej ab initio za pomocą pakietu UMD

Podsumowanie

Topnienia i płyny są wszechobecnymi wektorami transportu masowego w systemach naturalnych. Opracowaliśmy pakiet open-source do analizy ab initio symulacji dynamiki molekularnej takich systemów. Obliczamy właściwości strukturalne (wiązanie, klasteryzacja, specjacja chemiczna), transportowe (dyfuzja, lepkość) i termodynamiczne (widmo drgań).

Streszczenie

Opracowaliśmy oparty na Pythonie pakiet open-source do analizy wyników wynikających z symulacji ab initio dynamiki molekularnej płynów. Pakiet najlepiej nadaje się do zastosowań w systemach naturalnych, takich jak stopione krzemiany i tlenki, płyny na bazie wody i różne płyny nadkrytyczne. Pakiet jest zbiorem skryptów Pythona, które zawierają dwie główne biblioteki zajmujące się formatami plików i krystalografią. Wszystkie skrypty są uruchamiane w wierszu poleceń. Proponujemy uproszczony format do przechowywania trajektorii atomowych i istotnych informacji termodynamicznych symulacji, który jest zapisywany w plikach UMD, co oznacza Uniwersalną Dynamikę Molekularną. Pakiet UMD pozwala na obliczenie szeregu właściwości strukturalnych, transportowych i termodynamicznych. Zaczynając od funkcji rozkładu par, definiuje długości wiązań, buduje macierz połączeń międzyatomowych i ostatecznie określa specjację chemiczną. Określenie czasu życia gatunku chemicznego pozwala na przeprowadzenie pełnej analizy statystycznej. Następnie dedykowane skrypty obliczają przemieszczenia średniokwadratowe dla atomów, a także dla związków chemicznych. Zaimplementowana analiza samokorelacji prędkości atomowych pozwala uzyskać współczynniki dyfuzji i widmo drgań. Ta sama analiza zastosowana do naprężeń daje lepkość. Pakiet jest dostępny za pośrednictwem strony internetowej GitHub oraz na własnej stronie poświęconej projektowi ERC IMPACT w formie pakietu otwartego dostępu.

Wprowadzenie

Płyny i stopione substancje są aktywnymi chemicznymi i fizycznymi wektorami transportu w naturalnym środowisku. Podwyższone tempo dyfuzji atomowej sprzyja wymianom chemicznym i reakcjom, niska lepkość w połączeniu ze zmiennym wyporem sprzyja dużemu transferowi masy, a relacje gęstości kryształów i stopów sprzyjają warstwowaniu się wewnątrz ciał planetarnych. Brak sieci okresowej, typowe wysokie temperatury wymagane do osiągnięcia stanu stopionego oraz trudność w hartowaniu sprawiają, że eksperymentalne określenie szeregu oczywistych właściwości, takich jak gęstość, dyfuzja i lepkość, jest niezwykle trudne. Trudności te sprawiają, że alternatywne metody obliczeniowe są silnymi i użytecznymi narzędziami do badania tej klasy materiałów.

Wraz z pojawieniem się mocy obliczeniowej i dostępnością superkomputerów, dwie główne techniki numerycznych symulacji atomistycznych są obecnie wykorzystywane do badania stanu dynamicznego niekrystalicznego systemu atomistycznego, Monte Carlo¹ i dynamika molekularna (MD)^1,2. W symulacjach Monte Carlo przestrzeń konfiguracyjna jest próbkowana losowo; Metody Monte Carlo pokazują skalowanie liniowe w równoległości, jeśli wszystkie obserwacje próbkowania są od siebie niezależne. Jakość wyników zależy od jakości generatora liczb losowych i reprezentatywności próbkowania. Metody Monte Carlo pokazują skalowanie liniowe w równoległości, jeśli próbkowanie jest od siebie niezależne. W dynamice molekularnej (MD) przestrzeń konfiguracyjna jest próbkowana przez zależne od czasu trajektorie atomowe. Zaczynając od danej konfiguracji, trajektorie atomowe są obliczane przez całkowanie newtonowskich równań ruchu. Siły międzyatomowe można obliczyć za pomocą modelowych potencjałów międzyatomowych (w klasycznym MD) lub przy użyciu metod pierwszych zasad (in ab initio, czyli pierwszych zasad, MD). Jakość wyników zależy od długości trajektorii i jej zdolności do unikania przyciągania do lokalnych minimów.

Symulacje dynamiki molekularnej zawierają mnóstwo informacji, wszystkie związane z dynamicznym zachowaniem systemu. Średnie właściwości termodynamiczne, takie jak energia wewnętrzna, temperatura i ciśnienie, są raczej standardowe do obliczenia. Można je wyodrębnić z pliku wyjściowego symulacji i uśrednić, podczas gdy wielkości związane bezpośrednio z ruchem atomów, jak również z ich wzajemnym stosunkiem, muszą być obliczone po ekstrakcji pozycji i prędkości atomów.

W związku z tym, wiele wysiłku poświęcono wizualizacji wyników, a różne pakiety są dostępne dzisiaj, na różnych platformach, open source lub nie [Ovito³, VMD⁴, Vesta⁵, Travis⁶, itd.]. Wszystkie te narzędzia wizualizacyjne skutecznie radzą sobie z odległościami międzyatomowymi, a tym samym umożliwiają efektywne obliczanie funkcji rozkładu par i współczynników dyfuzji. Różne grupy przeprowadzające symulacje dynamiki molekularnej na dużą skalę mają zastrzeżone oprogramowanie do analizy różnych innych właściwości wynikających z symulacji, czasami w shareware lub innych formach ograniczonego dostępu do społeczności, a czasami ograniczone pod względem zakresu i zastosowania do niektórych konkretnych pakietów. W niektórych z tych pakietów opracowywane i implementowane są zaawansowane algorytmy do wydobywania informacji o wiązaniach międzyatomowych, wzorcach geometrycznych i termodynamice^3,4,5,6,7, itd.

Tutaj proponujemy pakiet UMD - pakiet open-source napisany w Pythonie do analizy wyników symulacji dynamiki molekularnej. Pakiet UMD pozwala na obliczanie szerokiego zakresu właściwości strukturalnych, dynamicznych i termodynamicznych (Rysunek 1). Pakiet jest dostępny za pośrednictwem strony internetowej GitHub (https://github.com/rcaracas/UMD_package) oraz na dedykowanej stronie (http://moonimpact.eu/umd-package/) projektu ERC IMPACT w formie pakietu otwartego dostępu.

Aby uczynić to uniwersalnym i łatwiejszym w obsłudze, nasze podejście polega na tym, aby najpierw wyodrębnić wszystkie informacje związane ze stanem termodynamicznym i trajektoriami atomowymi z pliku wyjściowego rzeczywistego przebiegu dynamiki molekularnej. Informacje te są przechowywane w dedykowanym pliku, którego format jest niezależny od oryginalnego pakietu MD, w którym przeprowadzono symulację. Nazywamy te pliki "umd", co oznacza Uniwersalną Dynamikę Molekularną. W ten sposób nasz pakiet UMD może być łatwo używany przez dowolną grupę ab initio z dowolnym oprogramowaniem, a wszystko to przy minimalnym wysiłku adaptacyjnym. Jedynym wymogiem do korzystania z obecnego pakietu jest zapisanie odpowiedniego parsera z danych wyjściowych określonego oprogramowania MD do formatu pliku umd, jeśli ten jeszcze nie istnieje. Na razie udostępniamy takie parsery dla pakietów VASP⁸ oraz QBox⁹.

Rysunek 1: Schemat blokowy biblioteki UMD.
Właściwości fizyczne są zaznaczone na niebiesko, a główne skrypty Pythona i ich opcje są na czerwono. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Pliki umd to pliki ASCII; typowe rozszerzenie to "umd.dat", ale nie jest to obowiązkowe. Wszystkie komponenty analizy mogą odczytywać pliki ASCII w formacie umd, niezależnie od rzeczywistego rozszerzenia nazwy. Jednak niektóre z automatycznych skryptów zaprojektowanych do szybkiego wykonywania statystyk na dużą skalę w kilku symulacjach specjalnie szukają plików z rozszerzeniem umd.dat. Każda właściwość fizyczna jest wyrażona w jednym wierszu. Każdy wiersz zaczyna się od słowa kluczowego. W ten sposób format jest wysoce elastyczny i pozwala na dodawanie nowych właściwości do pliku umd, przy jednoczesnym zachowaniu jego czytelności we wszystkich wersjach. Pierwsze 30 wierszy pliku umd z symulacji pirolitu przy 4,6 GPa i 3000 K, użytego poniżej w dyskusji, są pokazane w Rysunek 2.

Rysunek 2: Początek pliku umd opisującego symulację ciekłego pirolitu w temperaturze 4,6 GPa i temperaturze 3000 K.
Po nagłówku następuje opis każdej migawki. Każda właściwość jest zapisana w jednym wierszu, zawierającym nazwę właściwości fizycznej, wartość (wartości) i jednostki, wszystkie oddzielone spacjami. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Wszystkie pliki umd zawierają nagłówek opisujący zawartość komórki symulacyjnej: liczbę atomów, elektronów i typów atomów, a także szczegóły dotyczące każdego atomu, takie jak jego typ, symbol chemiczny, liczba elektronów walencyjnych i jego masa. Pusta linia oznacza koniec nagłówka i oddziela go od głównej części pliku umd.

Następnie każdy krok symulacji jest szczegółowo opisany. Po pierwsze, podane są chwilowe parametry termodynamiczne, każdy na innej linii, określając (i) nazwę parametru, taką jak energia, naprężenia, równoważne ciśnienie hydrostatyczne, gęstość, objętość, parametry sieci itp., (ii) jego wartość (wartości) oraz (iii) jego jednostki. Następna jest tabela opisująca atomy. Wiersz nagłówka zawiera różne miary, takie jak pozycje kartezjańskie, prędkości, ładunki itp., oraz ich jednostki. Następnie każdy atom jest szczegółowo opisany w jednym wierszu. W grupach po trzy, odpowiadających trzem osiom x, y, z, wpisy to: pozycje zredukowane, pozycje kartezjańskie złożone w komórce symulacji, pozycje kartezjańskie (które prawidłowo uwzględniają fakt, że atomy mogą przechodzić przez kilka komórek elementarnych podczas symulacji), prędkości atomowe i siły atomowe. Ostatnie dwa wpisy to skalary: ładunek i moment magnetyczny.

Dwie główne biblioteki zapewniają prawidłowe funkcjonowanie całego pakietu. Biblioteka umd_process.py zajmuje się plikami umd, takimi jak odczytywanie i drukowanie. Biblioteka crystallography.py zajmuje się wszystkimi informacjami związanymi z rzeczywistą strukturą atomową. Filozofia leżąca u podstaw biblioteki crystallography.py polega na traktowaniu kraty jako przestrzeni wektorowej. Parametry komórki elementarnej wraz z ich orientacją reprezentują wektory bazowe. "Przestrzeń" ma szereg atrybutów skalarnych (właściwa objętość, gęstość, temperatura i określona liczba atomów), właściwości termodynamiczne (energia wewnętrzna, ciśnienie, pojemność cieplna itp.) oraz szereg właściwości tensoralnych (naprężenie i sprężystość). Atomy zasiedlają tę przestrzeń. Klasa "Lattice" definiuje ten zespół, wraz z kilkoma krótkimi obliczeniami, takimi jak objętość właściwa, gęstość, otrzymanie odwrotności sieci z sieci bezpośredniej itp. Klasa "Atoms" definiuje atomy. Charakteryzują się one szeregiem właściwości skalarnych (nazwa, symbol, masa, liczba elektronów itp.) oraz szeregiem właściwości wektorowych (położenie w przestrzeni, albo względem bazy wektorowej opisanej w klasie Lattice, albo względem uniwersalnych współrzędnych kartezjańskich, prędkości, sił itp.). Oprócz tych dwóch klas biblioteka crystallography.py zawiera szereg funkcji do wykonywania różnych testów i obliczeń, takich jak odległości atomowe lub mnożenie komórek. Układ okresowy pierwiastków jest również dołączony jako słownik.

Różne komponenty pakietu umd zapisują kilka plików wyjściowych. Z reguły wszystkie są plikami ASCII, wszystkie ich wpisy są oddzielone zakładkami i są tak oczywiste, jak to tylko możliwe. Na przykład, zawsze wyraźnie wskazują one na własność fizyczną i jej jednostki. Pliki umd.dat są w pełni zgodne z tą regułą.

Protokół

1. Analiza przebiegów molekularno-dynamicznych

UWAGA: Pakiet jest dostępny na stronie internetowej GitHub (https://github.com/rcaracas/UMD_package) oraz na dedykowanej stronie (http://moonimpact.eu/umd-package/) projektu ERC IMPACT jako pakiet o otwartym dostępie.

1. Wyodrębnij każdy określony zestaw właściwości fizycznych przy użyciu co najmniej jednego dedykowanego skryptu języka Python z pakietu. Uruchom wszystkie skrypty w wierszu poleceń; Wszystkie wykorzystują serię flag, które są tak spójne, jak to tylko możliwe, od jednego skryptu do drugiego. Flagi, ich znaczenie i wartości domyślne podsumowano w tabeli 1.

flaga	znaczenie	Skrypt korzystający z niego	Wartość domyślna
-h	Krótka pomoc	cały
-f	Nazwa pliku UMD	cały
-ja	Etapy termalizacji, które należy wyrzucić	cały	0
-ja	Plik wejściowy zawierający wiązania międzyatomowe	Specjacji	obligacje.dane wejściowe
-s	Pobieranie próbek częstotliwości	MSD, specjacja	1 (każdy krok jest brany pod uwagę)
-a	Lista atomów lub anionów	Specjacji
-c	Wykaz kationów	Specjacji
-l	Długość wiązania	Specjacji	cyfra arabska
-t	temperatura	drgania, reologia
-v	Dyskretyzacja szerokości okna pobierania próbek trajektorii dla analizy przemieszczenia średniokwadratowego	Msd	20
-z	Dyskretyzacja początku okna pobierania próbek trajektorii dla analizy przemieszczenia średniokwadratowego	Msd	20

Tabela 1: Najczęściej używane flagi w pakiecie UMD i ich najczęstsze znaczenie.

2. Zacznij od przekształcenia danych wyjściowych symulacji MD wykonanej w kodzie opartych na pierwszych zasadach, takim jak VASP⁸ lub QBox⁹, do pliku UMD.
   1. Jeśli symulacje MD zostały wykonane w VASP, to w wierszu poleceń wpisz:
      VaspParser.py -f -i
      gdzie –f określa nazwę pliku VASP OUTCAR, a –i długość termalizacji.
      UWAGA: Krok początkowy, zdefiniowany przez –i, pozwala na pominięcie pierwszych kroków symulacji, które reprezentują termalizację. W typowym przebiegu dynamiki molekularnej pierwsza część obliczeń reprezentuje termalizację, tj. czas, w którym układ potrzebuje do opisania przez wszystkie atomy rozkładu temperatury podobnego do Gaussa i dla całego układu do wykazania wahań temperatury, ciśnienia, energii itd. wokół wartości równowagi. Ta termalizacja części symulacji nie powinna być brana pod uwagę podczas analizy właściwości statystycznych płynu.
3. Przekształć plik . umd do pliku . xyz, aby ułatwić wizualizację na różnych innych pakietach, takich jak VMD⁴ lub Vesta⁵. W wierszu poleceń wpisz:
   umd2xyz.py -f -i -s
   gdzie –f definiuje nazwę elementu . UMD, –I określa okres termalizacji, który ma zostać odrzucony, a –S częstotliwość próbkowania trajektorii przechowywanej w pliku . Plik UMD. Wartości domyślne to –i 0 –s 1, tj. uwzględniające wszystkie etapy symulacji, bez ich odrzucenia.
4. Za pomocą skryptu umd2poscar.py można odwrócić plik umd do plików POSCAR typu VASP; migawki symulacji można wybierać z predefiniowaną częstotliwością. W wierszu poleceń wpisz:
   umd2poscar.py -f -i -l -s
   gdzie –l oznacza ostatni krok, który ma zostać przekształcony w plik POSCAR. Wartości domyślne to -i 0 -l 10000000 -s 1. Ta wartość –l jest wystarczająco duża, aby pokryć typową całą trajektorię.

2. Wykonaj analizę statyczno-wytrzymałościową

1. Uruchom skrypt gofrs_umd.py, aby obliczyć funkcję rozkładu par (PDF) g_ᴀʙ(r) dla wszystkich par typów atomowych A i B (Rysunek 3). Dane wyjściowe są zapisywane w jednym pliku ASCII, rozdzielonym tabulatorami, z rozszerzeniem gofrs.dat. W wierszu poleceń wpisz:
   gofrs_umd.py -f -s < Sampling_Frequency > -d -i
   UWAGA: Wartości domyślne to Sampling_Frequency (częstotliwość próbkowania trajektorii) = 1 krok; Interwał dyskretyzacji (do wykreślania g(r)) = 0,01 A; InitialStep (liczba kroków na początku trajektorii, które są odrzucane) = 0. Promieniowy PDF, g_ᴀʙ(r) to średnia liczba atomów typu B w odległości d_ᴀʙ w sferycznej powłoce o promieniu r i grubości dr wyśrodkowanej na atomach typu A (Rysunek 3):
   
   
   przy ρ gęstości atomowej, N_A i N_B liczbie atomów typu A i B oraz δ(r−r_ᴀʙ) funkcji delta, która jest równa 1, jeśli atomy A i B leżą w odległości między r i r+dr. Odcięta pierwszego maksimum g_ᴀʙ(r) daje najwyższe prawdopodobieństwo długości wiązania między atomami typu A i B, która jest najbliższa średniej odległości wiązania, którą możemy określić. Pierwsze minimum wyznacza zakres pierwszej sfery koordynacyjnej. Stąd całka nad plikiem PDF do pierwszego minimum daje średnią liczbę koordynacyjną. Suma transformacji Fouriera g_ᴀʙ(r) dla wszystkich par typów atomowych A i B daje wzór dyfrakcyjny płynu, uzyskany eksperymentalnie za pomocą dyfraktometru. Jednak w rzeczywistości, ponieważ często w g_ᴀʙ(r) brakuje sfer koordynacyjnych wyższego rzędu, nie można uzyskać wzoru dyfrakcyjnego w całości.

Rysunek 3: Wyznaczanie funkcji rozkładu par.
a) Dla każdego atomu jednego gatunku (na przykład czerwonego) wszystkie atomy gatunków koordynujących (na przykład szary i/lub czerwony) są liczone w funkcji odległości. (b) Wynikowy wykres rozkładu odległości dla każdej migawki, która na tym etapie jest tylko zbiorem funkcji delta, jest następnie uśredniana dla wszystkich atomów i wszystkich migawek i ważona przez idealny rozkład gazu w celu wygenerowania (c) funkcji rozkładu par, która jest ciągła. Pierwsze minimum g(r) to promień pierwszej sfery koordynacyjnej, wykorzystany w dalszej części analizy specjacji. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

2. Wyodrębnij średnie odległości wiązań międzyatomowych jako promienie pierwszych sfer koordynacyjnych. W tym celu określ położenie pierwszego maksimum funkcji g_ᴀʙ(r): narysuj plik gofrs.dat w aplikacji arkusza kalkulacyjnego i wyszukaj maksima i minima dla każdej pary atomów.
3. Zidentyfikuj promień pierwszej sfery koordynacyjnej jako pierwsze minimum pliku PDF, g_ᴀʙ(r), za pomocą oprogramowania do obsługi arkuszy kalkulacyjnych. Jest to podstawa całej analizy statyczno-wytrzymałościowej płynu; plik PDF daje średni stan wiązania atomów w płynie.
4. Wyodrębnij odległości pierwszych minimów, czyli odciętych, i zapisz je w osobnym pliku, zwanym na przykład bonds.input. Alternatywnie uruchom jeden z analyze_gofr skryptów pakietu UMD, aby zidentyfikować maksima i minima funkcji g_ᴀʙ(r). W wierszu poleceń wpisz:
   analyze_gofr_semi_automatic.py
5. Kliknij pozycję maksimum i minimum funkcji g_ᴀʙ(r) wyświetlaną na wykresie otwieranym przez program. Skrypt automatycznie skanuje bieżący folder, identyfikuje wszystkie pliki gofrs.dat i przeprowadza analizę dla każdego z nich. Kliknij ponownie maksimum i minimum w oknie za każdym razem, gdy skrypt wymaga wykształconego początkowego zgadnięcia.
6. Otwórz i spójrz na automatycznie wygenerowany plik o nazwie bonds.input, który zawiera odległości wiązań międzyatomowych.

3. Wykonaj analizę specjacji

1. Oblicz topologię wiązań między atomami, używając koncepcji łączności w teorii grafów: atomy są węzłami, a wiązania międzyatomowe są ścieżkami. Skrypt speciation_umd.py wymaga odległości wiązań międzyatomowych zdefiniowanych w pliku bonds.input.
   UWAGA: Macierz połączeń jest konstruowana w każdym kroku czasowym: dwa atomy, które leżą w odległości mniejszej niż promień odpowiadającej im pierwszej sfery koordynacyjnej, są uważane za połączone, tj. połączone. Różne sieci atomowe są budowane poprzez traktowanie atomów jako węzłów na grafie, których połączenia są definiowane przez to kryterium geometryczne. Sieci te są gatunkami atomowymi, a ich zespół określa specjację atomową w tym konkretnym płynie (Rysunek 4).

Rysunek 4: Identyfikacja klastrów atomowych.
Wielościany koordynacyjne są definiowane za pomocą odległości międzyatomowych. Wszystkie atomy znajdujące się w odległości mniejszej niż określony promień są uważane za związane. Tutaj próg odpowiada pierwszej sferze koordynacyjnej (jasnoczerwone kółka), zdefiniowanej w Rysunek 1. Polimeryzacja, a co za tym idzie związki chemiczne, są otrzymywane z sieci związanych atomów. Zwróć uwagę na centralną gromadę Red₁Grey₂, która jest odizolowana od innych atomów, które tworzą nieskończony polimer. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

2. Uruchom skrypt specjacyjny, aby uzyskać macierz połączeń i uzyskać wielościany koordynacyjne lub polimeryzację. W wierszu poleceń wpisz:
   speciation_umd.py -f -s -i -l -c -a -m -r
   gdzie flaga -i daje plik z odległościami wiązań międzyatomowych, który został utworzony na przykład w poprzednim kroku. Alternatywnie, uruchom skrypt o jednej długości dla wszystkich wiązań zdefiniowanych przez -l flag.
   UWAGA: Flaga -c określa atomy centralne, a flaga -a ligandy. Zarówno atomy centralne, jak i ligandy mogą być różnych typów; W takim przypadku muszą być oddzielone przecinkami. Flaga -m określa minimalny czas, przez jaki gatunek musi żyć, aby mógł zostać uwzględniony w analizie. Domyślnie ten minimalny czas wynosi zero, a wszystkie wystąpienia są zliczane w końcowej analizie.
   1. Uruchom skrypt speciation_umd.py z flagą –r 0, która próbkuje graf połączeń na pierwszym poziomie w celu zidentyfikowania wielościanów koordynacyjnych. Na przykład centralny atom, oznaczony jako kation, może być otoczony przez jeden lub więcej anionów (Rysunek 4). Skrypt specjacyjny identyfikuje każdy z wielościanów koordynacyjnych. Średnia ważona wszystkich wielościanów koordynacyjnych daje liczbę koordynacyjną, identyczną z liczbą uzyskaną w wyniku integracji pliku PDF. W wierszu poleceń wpisz:
      speciation_umd.py -f -i -c -a -r 0
      UWAGA: Średnie liczby koordynacyjne w płynach są liczbami ułamkowymi. Ta frakcyjność wynika ze średniej charakterystyki koordynacji. Definicja oparta na specjacji daje bardziej intuicyjną i pouczającą reprezentację struktury płynu, w której określa się ilościowo względne proporcje różnych gatunków, tj. koordynacje.
   2. Uruchom skrypt speciation_umd.py z flagą –r 1, który próbkuje wykres połączeń na wszystkich poziomach głębokości w celu uzyskania polimeryzacji. Sieć przechodząca przez wykres atomowy ma pewną głębokość, ponieważ atomy są związane dalej z innymi wiązaniami (np. w sekwencje naprzemiennych kationów i anionów) (Rysunek 4).
3. Otwórz dwa pliki . popul.dat oraz . stat.dat kolejno; Stanowią one dane wyjściowe skryptu specjacyjnego. Każdy klaster jest zapisany w jednym wierszu, określając jego wzór chemiczny, czas powstania, czas śmierci, czas życia, matrycę z listą atomów tworzących ten klaster. Wykreśl czas życia każdego klastra atomowego wszystkich związków chemicznych znalezionych w symulacji, jak znaleziono w pliku .popul.dat (Rysunek 5).
4. Wykreśl analizę populacji z liczebnością każdego gatunku, jak znaleziono w . stat.dat pliku. Analiza ta, zarówno bezwzględna, jak i względna, odpowiada rzeczywistym statystykom wielościanów koordynacyjnych dla przypadku -r 0; W przypadku polimeryzacji, przy -R 1 należy to traktować ostrożnie, ponieważ może być konieczne zastosowanie pewnej normalizacji względem względnej liczby atomów. Obfitość odpowiada całki w ciągu wszystkich wcieleń. Ten. stat.dat plik zawiera również rozmiar każdego klastra, tj. ile atomów go tworzy.

4. Oblicz współczynniki dyfuzji

1. Wyodrębnij średnie przemieszczenia kwadratowe (MSD) atomów w funkcji czasu, aby uzyskać samodyfuzyjność. Standardowa formuła MSD to:
   
   gdzie prefaktorami są renormalizacje. Dzięki narzędziu MSD istnieją różne sposoby analizy dynamicznych aspektów płynów.
   UWAGA: T to całkowity czas symulacji, a N_α to liczba atomów typu α. Początkowy czas t₀ jest dowolny i obejmuje pierwszą połowę symulacji. N_init to liczba początkowych razy. τ oznacza szerokość przedziału czasu, w którym oblicza się zaburzenia układu mięśniowo-szkieletowego; Jego maksymalna wartość to połowa czasu trwania symulacji. W typowych implementacjach MSD każde okno zaczyna się na końcu poprzedniego. Jednak rzadsze próbkowanie może przyspieszyć obliczenia MSD, bez zmiany nachylenia MSD. W tym celu i-te okno rozpoczyna się o czasie t₀(i), ale (i+1)-te okno zaczyna się o czasie t₀(i) + τ + v, gdzie wartość v jest zdefiniowana przez użytkownika. Podobnie szerokość okna jest zwiększana w dyskretnych krokach zdefiniowanych przez użytkownika, w następujący sposób: τ(i) = τ(i-1) + z. Wartości z ("krok poziomy") i v ("krok pionowy") są dodatnie lub zerowe; Wartość domyślna dla obu tych opcji to 20.
2. Obliczanie MSD przy użyciu serii skryptów msd_umd. Ich dane wyjściowe są drukowane w formacie . msd.dat plik, w którym MSD każdego typu atomowego, atomu lub klastra jest drukowany w jednej kolumnie w funkcji czasu.
   1. Oblicz średnią wartość MSD każdego typu niepodzielnego. MSD są obliczane dla każdego atomu, a następnie uśredniane dla każdego typu atomu. Plik wyjściowy zawiera jedną kolumnę dla każdego typu niepodzielnego. W wierszu poleceń wpisz:
      msd_umd.py -f -z -v -b
   2. Oblicz MSD każdego atomu. MSD są obliczane dla każdego atomu, a następnie uśredniane dla każdego typu atomu. Plik wyjściowy zawiera jedną kolumnę dla każdego atomu w symulacji, a następnie jedną kolumnę dla każdego typu atomowego. Ta funkcja pozwala na identyfikację atomów, które dyfundują w dwóch różnych środowiskach, takich jak ciecz i gaz lub dwie ciecze. W wierszu poleceń wpisz:
      msd_all_umd.py -f -z -v -b
   3. Oblicz MSD związków chemicznych. Użyj populacji klastrów zidentyfikowanych za pomocą skryptu specjacyjnego i wydrukowanych w pliku . popul.dat pliku. Wartości MSD są obliczane dla każdego pojedynczego klastra. Plik wyjściowy zawiera jedną kolumnę dla każdego klastra. Aby uniknąć rozważania polimerów na dużą skalę, należy ograniczyć wielkość klastra; Jego wartość domyślna to 20 atomów. W wierszu poleceń wpisz:
      msd_cluster_umd.py -f -p -s -b -c
      UWAGA: Wartości domyślne to: –b 100 –s 1 –c 20.
3. Wykreśl MSD za pomocą oprogramowania opartego na arkuszu kalkulacyjnym (Rysunek 6). W logarytmicznej reprezentacji MSD w funkcji czasu zidentyfikuj zmianę nachylenia. Oddziel pierwszą część, zwykle krótką, która reprezentuje reżim balistyczny, czyli zasadę zachowania prędkości atomów po zderzeniach. Druga, dłuższa część reprezentuje reżim dyfuzyjny, czyli rozpraszanie prędkości atomów po zderzeniach.
4. Oblicz współczynniki dyfuzji na podstawie nachylenia MSD jako:
   
   gdzie Z jest liczbą stopni swobody (Z = 2 dla dyfuzji w płaszczyźnie, Z = 3 dla dyfuzji w przestrzeni), a t jest krokiem czasowym.

5. Funkcje korelacji czasowej

1. Oblicz funkcje korelacji czasowej jako miarę bezwładności układu, korzystając ze wzoru ogólnego:
   
   A może być różnymi zmiennymi zależnymi od czasu, takimi jak położenie atomów, prędkości atomowe, naprężenia, polaryzacja itp., z których każda daje - poprzez relacje Green-Kubo ^12,13 - różne właściwości fizyczne, czasami po dalszym przekształceniu.
2. Przeanalizuj prędkości atomowe, aby uzyskać widmo drgań cieczy i alternatywne wyrażenie atomowych współczynników samodyfuzji.
   1. Uruchom skrypt vibr_spectrum_umd.py, aby obliczyć funkcję automatycznej korelacji prędkości-prędkości atomowej (VAC) dla każdego typu atomowego i wykonać jej szybką transformację Fouriera. W wierszu poleceń wpisz:
      vibr_spectrum_umd.py -f -t
      gdzie –t oznacza temperaturę, która musi zostać zdefiniowana przez użytkownika. Skrypt drukuje dwa pliki: . vels.scf.dat plik z funkcją VAC dla każdego typu niepodzielnego, a plik . vibr.dat plik z widmem drgań rozłożonym na każdy gatunek atomowy i całkowitą wartość.
   2. Otwórz i przeczytaj vels.scf.dat. Wykreśl funkcję VAC z pliku vels.scf.dat za pomocą oprogramowania przypominającego arkusz kalkulacyjny.
   3. Zachowaj rzeczywistą część Fourier VAC. To właśnie daje widmo drgań w funkcji częstotliwości:
      
      gdzie m są masami atomowymi.
   4. Wykreśl widmo drgań z pliku vibr.dat za pomocą oprogramowania podobnego do arkusza kalkulacyjnego (Rysunek 7). Zidentyfikuj skończoną wartość przy ω = 0, która odpowiada dyfuzyjnemu charakterowi płynu i różnym pikom widma o skończonej częstotliwości. Zidentyfikuj udział każdego typu atomu w widmie oscylacyjnym.
      UWAGA: Rozkład typów atomowych pokazuje, że różne atomy mają różne wkłady ω=0, odpowiadające ich współczynnikom dyfuzji. Ogólny kształt widma jest znacznie gładszy i ma mniej cech niż w przypadku odpowiadającej mu bryły.
   5. Na powłoce odczytaj całkę w widmie drgań, co daje współczynniki dyfuzji dla każdego rodzaju atomowego.
      UWAGA: Właściwości termodynamiczne można uzyskać przez całkowanie ze spektrum oscylacyjnego, ale wyniki należy stosować ostrożnie ze względu na dwa przybliżenia: całkowanie jest ważne w przybliżeniu quasi-harmonicznym, które niekoniecznie zachodzi w wysokich temperaturach; a gazopodobna część widma odpowiadająca dyfuzji musi zostać odrzucona. Integracja powinna być wtedy wykonywana tylko na siatkowej części widma. Ale ta separacja zwykle wymaga kilku dalszych kroków przetwarzania końcowego i obliczeń¹⁴, które nie są objęte obecnym pakietem UMD.
3. Uruchom skrypt viscosity_umd.py, aby przeanalizować samokorelację tensora naprężeń składników w celu oszacowania lepkości stopu. W wierszu poleceń wpisz:
   viscosity_umd.py -f -i -s -o -l
   UWAGA: Ta funkcja ma charakter eksploracyjny i wszelkie wyniki należy traktować z ostrożnością. Po pierwsze, dokładnie sprawdź zbieżność lepkości w odniesieniu do długości symulacji.
   1. Wyprowadź lepkość płynu z własnej korelacji tensora naprężenia¹⁵ jako:
      
      gdzie V i T są odpowiednio objętością i temperaturą, κ_B jest stałą Boltzmanna, a σ_ij składową pozadiagonalną ij tensora naprężenia, wyrażoną we współrzędnych kartezjańskich.
   2. Użyj bardziej odpowiedniego dopasowania, aby uzyskać bardziej wiarygodne oszacowanie lepkości^15,16 i uniknąć szumu funkcji automatycznej korelacji tensora naprężenia, który może wynikać ze skończonego rozmiaru i skończonego czasu trwania symulacji. Dla funkcji automatycznej korelacji tensora naprężenia użyj następującej postaci funkcyjnej^15,16, która daje dobre wyniki:
      
      gdzie A, B, τ₁, τ₂ i ω są parametrami dopasowania. Po całkowaniu wyrażenie lepkości staje się następujące:
      

6. Parametry termodynamiczne wynikające z symulacji.

1. Uruchom averages.py, aby wyodrębnić średnie wartości i rozpiętość (jako odchylenie standardowe) dla ciśnienia, temperatury, gęstości i energii wewnętrznej z plików umd. W wierszu poleceń wpisz:
   averages.py -f -s
   z -s 0 jako domyślnym.
2. Oblicz błąd statystyczny średniej, korzystając z metod blokowania.
   UWAGA: Istnieją różne smaki tej metody. Podążając za pracami Allena i Tildesleya², powszechne jest uśrednianie sekwencji bloków czasowych o coraz większej długości i szacowanie odchylenia standardowego w odniesieniu do średniej arytmetycznej¹⁷. Konwergencja może być osiągnięta w granicy wielu i wystarczająco długich rozmiarów bloków, gdy próbkowanie nie jest skorelowane. Chociaż rzeczywista wartość progowa dla zbieżności zwykle musi być wybrana ręcznie.
   1. Użyj metody halvingu¹⁸: zaczynając od początkowej próbki danych, w każdym kroku κ, zmniejsz o połowę liczbę próbek, uśredniając przez każde dwie odpowiadające im kolejne próbki z poprzedniego kroku κ−1:
      
   2. Uruchom skrypt fullaverages.py, aby przeprowadzić pełną analizę statystyczną, w tym błąd średniej. W wierszu poleceń wpisz:
      fullaverages.py -s -u
      UWAGA: Skrypt jest zautomatyzowany do tego stopnia, że wyszukuje wszystkie pliki .umd.dat w bieżącym katalogu i przeprowadza analizę dla wszystkich z nich. Wartości domyślne to –s 0 –u 0. Dla -u 0 wyjście jest minimalne, a dla -u 1 wyjście jest pełne, z wypisanymi kilkoma alternatywnymi jednostkami. Ten skrypt wymaga wsparcia graficznego, ponieważ tworzy obraz graficzny do sprawdzania zbieżności w celu oszacowania błędu średniej.

Wyniki

Pirolit to modelowy wieloskładnikowy stopiony krzemian (0,5Na_{, 2O}, 2CaO, 1,5Al,_{2O, 3}, 4FeO, 30MgO, 24SiO₂), który najlepiej przybliża skład Ziemi krzemianowej - średnią geochemiczną lub całą naszą planetę, z wyjątkiem jej jądra opartego na żelazie¹⁹. Wczesna Ziemia była zdominowana przez serię wielkoskalowych topnień ²⁰, ostatnie z nich mogło pochłonąć całą planetę, po jej zagęszczeniu na dysku protoksiężycowym²¹. Pirolit stanowi najlepsze przybliżenie składu chemicznego takich oceanów magmy w skali planetarnej. W związku z tym szeroko zbadaliśmy właściwości fizyczne stopionego pirolitu w zakresie temperatur 3 000–20125 000 K i zakresu ciśnienia 0–150 GPa z symulacji dynamiki molekularnej ab initio w implementacji VASP. Te warunki termodynamiczne całkowicie charakteryzują najbardziej ekstremalne warunki oceanu magmy na Ziemi. Nasze badanie jest doskonałym przykładem udanego wykorzystania pakietu UMD do całej dogłębnej analizy melts²². Obliczyliśmy rozkład i średnie długości wiązań, prześledziliśmy zmiany w koordynacji kationowo-tlenowej i porównaliśmy nasze wyniki z wcześniejszymi badaniami eksperymentalnymi i obliczeniowymi na amorficznych krzemianach o różnym składzie. Nasza dogłębna analiza pomogła zdekomponować standardowe liczby koordynacyjne na ich podstawowe składniki, określić obecność egzotycznych wielościanów koordynacyjnych w stopie i wyodrębnić czasy życia wszystkich wielościanów koordynacyjnych. Podkreślono również znaczenie próbkowania w symulacjach zarówno pod względem długości trajektorii, jak i liczby atomów obecnych w modelowanym układzie. Jeśli chodzi o post-processing, analiza UMD jest niezależna od tych czynników, jednak powinny one być brane pod uwagę przy interpretacji wyników dostarczanych przez pakiet UMD. W tym miejscu pokazujemy kilka przykładów, w jaki sposób pakiet UMD może być wykorzystany do wydobycia kilku charakterystycznych cech stopów, z zastosowaniem do stopionego pirolitu.

Funkcja rozkładu par Si-O uzyskana ze skryptu gofrs_umd.py pokazuje, że promień pierwszej sfery koordynacyjnej, która jest pierwszym minimum funkcji g(r), leży około 2,5 angstrema przy T = 3000 K i P = 4,6 GPa. Maksimum g(r) wynosi 1,635 A - jest to najlepsze przybliżenie długości gięcia. Długi ogon wynika z temperatury. Wykorzystując ten limit jako odległość wiązania Si-O, analiza specjacji pokazuje, że jednostki SiO₄, które mogą trwać do kilku pikosekund, dominują w topnieniu (Rysunek 5). Istnieje ważna część stopu, która wykazuje częściową polimeryzację, co odzwierciedla obecność dimerów, takich jak Si₂O₇ i trimerów, takich jak jednostki Si₃O_x. Odpowiadający im czas życia jest rzędu pikosekundy. Wszystkie polimery wyższego rzędu mają znacznie krótszą żywotność.

Rysunek 5: Czas życia związków chemicznych Si-O.
Specjacja jest identyfikowana w wieloskładnikowym stopie w temperaturze 4,6 GPa i 3000 K. Etykiety oznaczają monomery SiO₃, SiO₄ i SiO₅ oraz różne polimery Si_xO_y. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Różne wartości kroków pionowych i poziomych, zdefiniowane przez flagi –z i –v powyżej, dają różne próbki MSD (Rysunek 6). Nawet duże wartości z i v wystarczą, aby określić nachylenia, a tym samym współczynniki dyfuzji różnych atomów. Oszczędność czasu na przetwarzanie końcowe jest niezwykła w przypadku przechodzenia do dużych wartości z i v. MSD oferuje bardzo silne kryterium walidacji jakości symulacji. Jeśli część dyfuzyjna MSD nie jest wystarczająco długa, oznacza to, że symulacja jest zbyt krótka i nie osiąga stanu płynnego w sensie statystycznym. Minimalne wymagania dotyczące części dyfuzyjnej układu mięśniowo-szkieletowego w dużym stopniu zależą od systemu. Można wymagać, aby wszystkie atomy przynajmniej raz zmieniły swoje miejsce w strukturze stopu, aby można go było uznać za płyn¹⁰. Doskonałym przykładem mającym zastosowanie w naukach planetarnych są złożone topnienia krzemianów pod wysokim ciśnieniem bliskim lub nawet poniżej ich linii likwidusu¹¹. Atomy Si, główne kationy tworzące sieć, zamieniają się miejscami po ponad dwudziestu pikosekundach. Symulacje krótsze niż ten próg oznaczałyby znaczne zaniżenie próbkowania możliwej przestrzeni konfiguracyjnej. Ponieważ jednak aniony koordynujące, a mianowicie atomy O, poruszają się szybciej niż centralne atomy Si, mogą one kompensować część powolnej ruchliwości Si. W związku z tym cały system może rzeczywiście pokryć lepsze próbkowanie przestrzeni konfiguracyjnej, niż zakłada się tylko na podstawie przemieszczeń Si.

Rysunek 6: Przemieszczenia średniokwadratowe (MSD).
MSD zilustrowano dla kilku typów atomowych wieloskładnikowego stopionego krzemianu. Pobieranie próbek z różnymi krokami poziomymi i pionowymi, z i v, daje spójne wyniki. Okręgi pełne: -z 50 –v 50. Otwarte okręgi: -z 250 –v 500. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Na koniec, atomowe funkcje VAC dają widmo drgań stopu. Rysunek 7 pokazuje widmo w tych samych warunkach ciśnienia i temperatury jak powyżej. Reprezentujemy wkład atomów Mg, Si i O, a także wartość całkowitą. Przy zerowej częstotliwości występuje skończona wartość widma, która odpowiada dyfuzyjnemu charakterowi stopu. Ekstrakcja właściwości termodynamicznych ze spektrum drgań musi usunąć ten gazowy charakter dyfuzyjny od zera, ale także odpowiednio uwzględnić jego rozpad przy wyższych częstotliwościach.

Rysunek 7: Widmo drgań stopionego pirolitu.
Rzeczywista część transformaty Fouriera funkcji samokorelacji prędkość atomowa daje widmo drgań. Tutaj widmo jest obliczane dla wieloskładnikowego stopionego krzemianu. Płyny mają niezerowy gazowy charakter dyfuzyjny przy zerowej częstotliwości. Kliknij tutaj, aby zobaczyć większą wersję tego rysunku.

Dyskusja

Pakiet UMD został zaprojektowany tak, aby lepiej działał z symulacjami ab initio, w których liczba migawek jest zwykle ograniczona do dziesiątek do setek tysięcy migawek, z kilkuset atomami na komórkę elementarną. Większe symulacje są również wykonalne, pod warunkiem, że maszyna, na której działa przetwarzanie końcowe, ma wystarczającą ilość zasobów pamięci aktywnej. Kod wyróżnia się różnorodnością właściwości, które może obliczyć, oraz licencją open source.

Pliki umd.dat są odpowiednie dla zespołów, które zachowują niezmienioną liczbę cząstek w całej symulacji. Pakiet UMD może odczytywać pliki pochodzące z obliczeń, w których kształt i objętość pola symulacji są różne. Obejmują one najczęstsze obliczenia, takie jak NVT i NPT, w których liczba cząstek, N, temperatura T, objętość, V i/lub ciśnienie P, są utrzymywane na stałym poziomie.

Na początku czasu funkcja rozkładu par, jak również wszystkie skrypty potrzebne do oszacowania odległości międzyatomowych, takie jak skrypty specjacyjne, działają tylko dla ortogonalnych komórek elementarnych, czyli dla komórek sześciennych, czworokątnych i rombowych, gdzie kąty między osiami wynoszą 90°.

Głównymi kierunkami rozwoju wersji 2.0 są usunięcie ograniczenia ortogonalności dla odległości i dodanie większej liczby funkcji do skryptów specjacyjnych: analiza poszczególnych wiązań chemicznych, analiza kątów międzyatomowych i implementacja drugiej sfery koordynacyjnej. Z pomocą współpracy zewnętrznej pracujemy nad przeniesieniem kodu na procesor graficzny w celu szybszej analizy w większych systemach.

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
getopt biblioteka
			kodzie
źródłowym biblioteka matplotlib biblioteka
		systemu operacyjnego	otwartym kodzie źródłowym
	Python Software Foundation	w wersji 2 i 3	Biblioteka
		losowa	typu open source Biblioteka
			typu open source Biblioteka
			kodzie źródłowym Biblioteka
			source Biblioteka
			sys