Artykuł metodologiczny

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

6.2K wyświetleń

⸱

DOI:

10.3791/61534

⸱

17 września 2021

 ,  ,  ,  ,  , 
1Ecole Normale Supérieure de Lyon, Laboratory of Geology of Lyon UMR5276, CNRS, 2Centre for Earth Evolution and Dynamics (CEED), University of Oslo, 3CEA, DAM, DIF

W tym artykule

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 stopy są aktywnymi wektorami transportu chemicznego i fizycznego w środowiskach naturalnych. Podwyższone tempo dyfuzji atomowej sprzyja wymianom chemicznym i reakcjom, niska lepkość w połączeniu ze zmienną wypornością sprzyja dużemu transportowi masy, a relacje gęstości między kryształem a stopem sprzyjają warstwowaniu wewnątrz ciał planetarnych. Brak sieci okresowej, typowo wysokie temperatury wymagane do osiągnięcia stanu stopionego oraz trudności z hartowaniem sprawiają, że eksperymentalne wyznaczenie szeregu oczywistych właściwości, takich jak gęstość, dyfuzja i lepkość, jest niezwykle wymagające. Trudności te sprawiają, że alternatywne metody obliczeniowe stanowią potężne i użyteczne narzędzia do badania tej klasy materiałów.

Wraz z rozwojem mocy obliczeniowej i dostępnością superkomputerów, do badania stanu dynamicznego niekrystalicznego układu atomistycznego stosuje się obecnie dwie główne techniki numerycznych symulacji atomistycznych: metodę Monte Carlo1 oraz dynamikę molekularną (MD)1,2. W symulacjach Monte Carlo przestrzeń konfiguracyjna jest próbkowana losowo; metody Monte Carlo wykazują liniową skalowalność w obliczeniach równoległych, jeśli wszystkie obserwacje z próbek są od siebie niezależne. Jakość wyników zależy od jakości generatora liczb pseudolosowych oraz reprezentatywności próbkowania. Metody Monte Carlo wykazują liniową skalowalność w paralelizacji, jeśli próbkowanie jest niezależne. W dynamice molekularnej (MD) przestrzeń konfiguracyjna jest próbkowana za pomocą zależnych od czasu trajektorii atomowych. Punktem wyjścia jest określona konfiguracja, a trajektorie atomowe są obliczane poprzez całkowanie newtonowskich równań ruchu. Siły międzyatomowe mogą być obliczane przy użyciu modelowych potencjałów międzyatomowych (w klasycznym MD) lub metod pierwszych zasad (w MD ab initio, czyli z pierwszych zasad). Jakość wyników zależy od długości trajektorii oraz jej zdolności do unikania uwięzienia w minimach lokalnych.

Symulacje dynamiki molekularnej zawierają ogromną ilość informacji, wszystkie związane z zachowaniem dynamicznym układu. Właściwości uśrednione termodynamicznie, takie jak energia wewnętrzna, temperatura i ciśnienie, są stosunkowo standardowe do obliczenia. Można je wyodrębnić z plików wyjściowych symulacji i uśrednić, natomiast wielkości bezpośrednio związane z ruchem atomów, jak i z ich wzajemnymi relacjami, muszą zostać obliczone po wyekstrahowaniu położeń i prędkości atomów.

W związku z tym ogromny wysiłek poświęcono wizualizacji wyników; obecnie na różnych platformach dostępnych jest wiele pakietów, zarówno otwartoźródłowych, jak i komercyjnych [Ovito3, VMD4, Vesta5, Travis6, itd.]. Wszystkie te narzędzia do wizualizacji efektywnie przetwarzają odległości międzyatomowe, a co za tym idzie, umożliwiają sprawne obliczanie funkcji rozkładu radialnego oraz współczynników dyfuzji. Różne grupy przeprowadzające wielkoskalowe symulacje dynamiki molekularnej dysponują autorskim oprogramowaniem do analizy innych właściwości wynikających z symulacji, które czasem jest udostępniane jako shareware lub w innych formach ograniczonego dostępu dla społeczności, a niekiedy ma ograniczony zakres i zastosowanie w ramach konkretnych pakietów. W niektórych z tych pakietów opracowano i zaimplementowano zaawansowane algorytmy służące do wyodrębniania informacji o wiązaniach międzyatomowych, wzorcach geometrycznych i termodynamice3,4,5,6,7, itd.

W niniejszej pracy proponujemy pakiet UMD – otwartoźródłowy pakiet napisany w języku Python do analizy wyników symulacji dynamiki molekularnej. Pakiet UMD umożliwia obliczanie szerokiego zakresu właściwości strukturalnych, dynamicznych i termodynamicznych (Rysunek 1). Pakiet jest dostępny za pośrednictwem strony GitHub (https://github.com/rcaracas/UMD_package) oraz dedykowanej strony (http://moonimpact.eu/umd-package/) projektu ERC IMPACT jako pakiet w otwartym dostępie.

Aby uczynić to rozwiązanie uniwersalnym i łatwiejszym w obsłudze, nasze podejście polega na tym, aby najpierw wyodrębnić z pliku wyjściowego właściwego przebiegu dynamiki molekularnej wszystkie informacje dotyczące stanu termodynamicznego oraz trajektorii atomowych. Informacje te są przechowywane w dedykowanym pliku, którego format jest niezależny od oryginalnego pakietu MD, w którym przeprowadzono symulację. Pliki te nazywamy plikami „umd”, co oznacza Universal Molecular Dynamics. Dzięki temu nasz pakiet UMD może być z łatwością wykorzystany przez każdą grupę zajmującą się metodami ab initio z dowolnym oprogramowaniem, przy minimalnym nakładzie pracy związanej z adaptacją. Jedynym wymogiem korzystania z obecnego pakietu jest napisanie odpowiedniego parsera, który przetworzy dane wyjściowe konkretnego oprogramowania MD na format pliku umd, jeśli taki parser jeszcze nie istnieje. Obecnie udostępniamy takie parsery dla pakietów VASP8 oraz QBox9.

Schemat analizy danych dynamiki molekularnej; wiązania międzyatomowe, dyfuzja atomowa, termodynamika.
Rycina 1: Schemat blokowy biblioteki UMD.
Właściwości fizyczne zaznaczono kolorem niebieskim, a główne skrypty Python wraz z ich opcjami kolorem czerwonym. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Pliki umd są plikami ASCII; typowe rozszerzenie to „umd.dat”, choć nie jest ono obowiązkowe. Wszystkie komponenty analityczne potrafią odczytywać pliki ASCII w formacie umd, niezależnie od rzeczywistego rozszerzenia nazwy. Jednakże niektóre skrypty automatyczne, zaprojektowane do szybkiego przeprowadzania statystyk na wielką skalę w wielu symulacjach, poszukują konkretnie plików z rozszerzeniem umd.dat. Każda właściwość fizyczna jest zapisana w jednej linii. Każda linia zaczyna się od słowa kluczowego. Dzięki temu format jest wysoce adaptacyjny i pozwala na dodawanie nowych właściwości do pliku umd, zachowując jednocześnie jego czytelność pomiędzy wersjami. Pierwsze 30 linii pliku umd z symulacji pyrolitu przy 4,6 GPa i 3000 K, wykorzystanej poniżej w dyskusji, przedstawiono na Ryc. 2.

Dane wyjściowe z symulacji dynamiki molekularnej; energia, tensor naprężeń, parametry strukturalne.
Rysunek 2: Początek pliku umd opisującego symulację ciekłego pyrolitu przy 4.6 GPa i 3000 K.
Za nagłówkiem następuje opis każdego migawki (snapshot). Każda właściwość jest zapisana w jednej linii i zawiera nazwę właściwości fizycznej, wartość(y) oraz jednostki, oddzielone od siebie spacjami. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

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

Następnie szczegółowo opisany jest każdy krok symulacji. Najpierw podawane są chwilowe parametry termodynamiczne, z których każdy znajduje się w osobnej linii, określając (i) nazwę parametru, np. energię, naprężenia, równoważne ciśnienie hydrostatyczne, gęstość, objętość, parametry sieci itp., (ii) jego wartość(y) oraz (iii) jednostki. Następnie następuje tabela opisująca atomy. Linia nagłówka zawiera różne miary, takie jak pozycje kartezjańskie, prędkości, ładunki itp. oraz ich jednostki. Potem każdy atom jest szczegółowo opisany w jednej linii. W grupach po trzy, odpowiadających trzem osiom x, y, z, wpisy stanowią: pozycje zredukowane, pozycje kartezjańskie złożone do komórki symulacyjnej, pozycje kartezjańskie (uwzględniające fakt, że atomy mogą przejść przez kilka komórek elementarnych podczas symulacji), prędkości atomowe oraz siły atomowe. Dwie ostatnie pozycje to wartości skalarne: ładunek i moment magnetyczny.

Za prawidłowe funkcjonowanie całego pakietu odpowiadają dwie główne biblioteki. Biblioteka umd_process.py zajmuje się plikami umd, na przykład ich odczytywaniem i drukowaniem. Biblioteka crystallography.py przetwarza wszystkie informacje związane z właściwą strukturą atomową. Podstawową filozofią biblioteki crystallography.py jest traktowanie sieci jako przestrzeni wektorowej. Parametry komórki elementarnej wraz z ich orientacją reprezentują wektory bazy. „Przestrzeń” posiada szereg atrybutów skalarnych (objętość właściwa, gęstość, temperatura i konkretna liczba atomów), właściwości termodynamiczne (energia wewnętrzna, ciśnienie, ciepło właściwe itp.) oraz szereg właściwości tensorowych (naprężenie i sprężystość). Przestrzeń ta jest wypełniona atomami. Klasa „Lattice” definiuje ten zespół, a także zawiera kilka krótkich obliczeń, takich jak objętość właściwa, gęstość, wyznaczanie sieci odwrotnej z bezpośredniej itd. Klasa „Atoms” definiuje atomy. Są one charakteryzowane przez szereg właściwości skalarnych (nazwa, symbol, masa, liczba elektronów itp.) oraz szereg 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ły itp.). Poza tymi dwiema klasami biblioteka crystallography.py zawiera szereg funkcji służących do przeprowadzania różnych testów i obliczeń, takich jak odległości atomowe czy mnożenie komórek. Dołączony jest również układ okresowy pierwiastków w formie słownika.

Różne komponenty pakietu umd generują kilka plików wyjściowych. Co do zasady są to pliki ASCII, w których wszystkie wpisy są oddzielone tabulatorami, a ich struktura jest maksymalnie intuicyjna. Na przykład zawsze jasno wskazują one właściwość fizyczną oraz jednostki. Pliki umd.dat w pełni spełniają tę zasadę.

Protokół

1. Analiza przebiegów dynamiki molekularnej

UWAGA: Pakiet jest dostępny za pośrednictwem strony GitHub (https://github.com/rcaracas/UMD_package) oraz na dedykowanej stronie (http://moonimpact.eu/umd-package/) projektu ERC IMPACT jako pakiet z otwartym dostępem.

  1. Wyodrębnij każdy konkretny zestaw właściwości fizycznych, korzystając z jednego lub kilku dedykowanych skryptów Python z pakietu. Uruchom wszystkie skrypty w wierszu poleceń; wszystkie one wykorzystują zestaw flag, które są w miarę możliwości spójne między poszczególnymi skryptami. Flagi, ich znaczenie oraz wartości domyślne zostały podsumowane w Tabeli 1.
OznaczenieZnaczenieSkrypt z jego użyciemWartość domyślna
-hKrótka pomocwszystko
-fnazwa pliku UMDWszystko
-iEtapy termalizacji do odrzuceniaWszystko0
-iPlik wejściowy zawierający wiązania międzyatomowespecjacjaWięzi.input
-sPróbkowanie częstotliwościMSD, specjacja1 (każdy krok jest brany pod uwagę)
-aLista atomów lub anionówspecjacja
-cLista kationówspecjacja
-lDługość wiązaniaspecjacja2
-tTemperaturawibracje, reologia
-vDyskretyzacja szerokości okna próbkowania trajektorii dla analizy średnio-kwadratowego przemieszczeniamsd20
-zDyskretyzacja początku okna próbkowania trajektorii dla analizy średniokwadratowego przemieszczeniamsd20

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

  1. Zacznij od przekształcenia wyników symulacji MD przeprowadzonej w kodzie opartym na zasadach pierwszych, takim jak VASP8 lub QBox9, do pliku UMD.
    1. Jeśli symulacje MD zostały wykonane w VASP, wpisz w wierszu poleceń:
      VaspParser.py -f -i
      gdzie flaga –f definiuje nazwę pliku OUTCAR z programu VASP, a –i długość termalizacji.
      UWAGA: Krok wstępny, definiowany przez –i, pozwala na odrzucenie pierwszych kroków symulacji, które reprezentują termalizację. W typowym przebiegu dynamiki molekularnej pierwsza część obliczeń reprezentuje termalizację, t. j. czas potrzebny systemowi na to, aby wszystkie atomy opisały rozkład temperatury zbliżony do rozkładu Gaussa, a cały system wykazywał fluktuacje temperatury, ciśnienia, energii itp. wokół wartości równowagowych. Tę część termalizacji symulacji nie należy brać pod uwagę podczas analizy właściwości statystycznych płynu.
  2. Przekształć pliki .umd w pliki .xyz, aby ułatwić wizualizację w różnych innych pakietach, takich jak VMD4 lub Vesta5. W wierszu poleceń wpisz:
    umd2xyz.py -f -i -s
    gdzie –f definiuje nazwę pliku .umd, –i definiuje okres termalizacji do odrzucenia, a –s częstotliwość próbkowania trajektorii zapisanej w pliku .umd. Wartości domyślne to –i 0 –s 1, co oznacza uwzględnienie wszystkich kroków symulacji, bez odrzucania któregokolwiek z nich.
  3. Przekonwertuj plik umd z powrotem do plików POSCAR typu VASP za pomocą skryptu umd2poscar.py; migawki z symulacji można wybrać z predefined częstotliwością. W wierszu poleceń wpisz:
    umd2poscar.py -f -i -l -s
    gdzie –l reprezentuje ostatni krok do przekształcenia w plik POSCAR. Wartości domyślne to -i 0 -l 10000000 -s 1. Ta wartość –l jest wystarczająco duża, aby objąć typową pełną trajektorię.

2. Przeprowadzenie analizy strukturalnej

  1. Uruchom skrypt gofrs_umd.py, aby obliczyć funkcję rozkładu radialnego (PDF) gᴀʙ(r) dla wszystkich par typów atomowych A i B (Rycina 3). Wynik jest zapisywany w jednym pliku ASCII z rozdzielaczem tabulatorem i 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; DiscretizationInterval (do wykreślania g(r)) = 0.01 Å; InitialStep (liczba kroków na początku trajektorii, które zostają pominięte) = 0. Radialna PDF, gᴀʙ(r), to średnia liczba atomów typu B w odległości d_ᴀʙ w obrębie sferycznej powłoki o promieniu r i grubości dr z centrum w atomach typu A (Rycina 3):

    Równanie funkcji rozkładu radialnego, \( g_{AB}(r) \), do analizy strukturalnej na schemacie naukowym.
    gdzie ρ to gęstość atomowa, NA i NB to liczba atomów typu A i B, a δ(r−rᴀʙ) to funkcja delta, która przyjmuje wartość 1, jeśli atomy A i B znajdują się w odległości między r a r+dr. Odscissa pierwszego maksimum gᴀʙ(r) wskazuje najbardziej prawdopodobną długość wiązania między atomami typu A i B, co stanowi najbliższą możliwą do wyznaczenia średnią odległość wiązania. Pierwsze minimum wyznacza zasięg pierwszej sfery koordynacyjnej. Zatem całka z PDF do pierwszego minimum daje średnią liczbę koordynacyjną. Suma transformat Fouriera gᴀʙ(r) dla wszystkich par typów atomowych A i B daje obraz dyfrakcyjny cieczy, taki jaki uzyskuje się eksperymentalnie za pomocą dyfraktometru. Jednak w rzeczywistości, ponieważ często w gᴀʙ(r) brakuje sfer koordynacyjnych wyższego rzędu, obraz dyfrakcyjny nie może zostać uzyskany w całości.

Diagram radialnego rozkładu atomów; wykresy liczby atomów i funkcji rozkładu par.
Rycina 3: Wyznaczanie funkcji rozkładu par.
(a) Dla każdego atomu jednego gatunku (na przykład czerwonego) wszystkie atomy gatunku koordynującego (na przykład szarego i/lub czerwonego) są liczone jako funkcja odległości. (b) Wynikowy wykres rozkładu odległości dla każdego migawki, który na tym etapie jest jedynie zbiorem funkcji delta, jest następnie uśredniany dla wszystkich atomów i wszystkich migawek oraz ważony rozkładem gazu doskonałego, aby wygenerować (c) ciągłą funkcję rozkładu par. Pierwsze minimum g(r) stanowi promień pierwszej sfery koordynacyjnej, wykorzystywany później w analizie specjacji. Kliknij tutaj, aby zobaczyć powiększoną wersję tej ryciny.

  1. Wyznacz średnie odległości wiązań międzyatomowych jako promienie pierwszych sfer koordynacyjnych. W tym celu zidentyfikuj położenie pierwszego maksimum funkcji gᴀʙ(r): wykreśl dane z pliku gofrs.dat w aplikacji arkusza kalkulacyjnego i wyszukaj maksima oraz minima dla każdej pary atomów.
  2. Zidentyfikuj promień pierwszej sfery koordynacyjnej jako pierwsze minimum PDF, gᴀʙ(r), korzystając z oprogramowania arkusza kalkulacyjnego. Stanowi to podstawę całej analizy strukturalnej cieczy; PDF dostarcza informacji o średnim stanie wiązań atomów w cieczy.
  3. Wyodrębnij wartości odległości pierwszych minimów, tj. wartości odciętych, i zapisz je w osobnym pliku, nazwanym na przykład bonds.input. Alternatywnie uruchom jeden ze skryptów analyze_gofr z pakietu UMD, aby zidentyfikować maksima i minima funkcji gᴀʙ(r). W wierszu poleceń wpisz:
    analyze_gofr_semi_automatic.py
  4. Kliknij w położenie maksimum i minimum funkcji gᴀʙ(r) wyświetlonej na wykresie otwieranym przez program. Skrypt automatycznie skanuje bieżący folder, identyfikuje wszystkie pliki gofrs.dat i przeprowadza analizę dla każdego z nich. Klikaj ponownie w maksimum i minimum w oknie za każdym razem, gdy skrypt będzie wymagał wstępnego przybliżenia.
  5. Otwórz i przejrzyj automatycznie wygenerowany plik o nazwie bonds.input, który zawiera odległości wiązań międzyatomowych.

3. Przeprowadzenie analizy specjacyjnej

  1. Oblicz topologię wiązań między atomami, wykorzystując koncepcję spójności w ramach teorii grafów: atomy stanowią węzły, a wiązania międzyatomowe są ścieżkami. Skrypt speciation_umd.py wymaga zdefiniowania odległości wiązań międzyatomowych w pliku bonds.input.
    UWAGA: Macierz spójności jest konstruowana w każdym kroku czasowym: dwa atomy, których odległość jest mniejsza niż promień ich odpowiedniej pierwszej sfery koordynacyjnej, są uznawane za związane, tj. połączone. Różne sieci atomowe są budowane poprzez traktowanie atomów jako węzłów w grafie, którego połączenia są definiowane przez to kryterium geometryczne. Sieci te stanowią gatunki atomowe, a ich zespół definiuje specjację atomową w danej cieczy (Rycina 4).

Schemat struktury krystalicznej; wizualizacja układu cząsteczkowego; koncepcja fizyki ciała stałego.
Rysunek 4: Identyfikacja klastrów atomowych.
Wielościany koordynacyjne są definiowane przy użyciu odległości międzyatomowych. Wszystkie atomy w odległości mniejszej niż określony promień są uznawane za związane. W tym przypadku próg odpowiada pierwszej sferze koordynacyjnej (jasnoczerwone okręgi), zdefiniowanej na Rysunku 1. Polimeryzację, a tym samym gatunki chemiczne, wyznacza się z sieci związanych atomów. Należy zwrócić uwagę na centralny klaster Red1Grey2, który jest odizolowany od pozostałych atomów tworzących nieskończony polimer. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

  1. Uruchom skrypt specjacji w celu uzyskania macierzy łączności oraz wielościanów koordynacyjnych lub stopnia polimeryzacji. W wierszu poleceń wpisz:
    speciation_umd.py -f -s -i -l -c -a -m -r
    gdzie flaga -i wskazuje plik z odległościami wiązań międzyatomowych, który został wygenerowany np. w poprzednim kroku. Alternatywnie uruchom skrypt z jedną, stałą długością dla wszystkich wiązań, zdefiniowaną za pomocą flagi -l.
    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 należy oddzielić je przecinkami. Flaga -m określa minimalny czas życia gatunku, aby został on uwzględniony w analizie. Domyślnie czas ten wynosi zero, co oznacza, że wszystkie wystąpienia są liczone w końcowej analizie.
    1. Uruchom skrypt speciation_umd.py z flagą –r 0, która analizuje graf łączności na pierwszym poziomie w celu zidentyfikowania wielościanów koordynacyjnych. Na przykład atom centralny, określany jako kation, może być otoczony przez jeden lub więcej anionów (Rysunek 4). Skrypt specjacji identyfikuje każdy z wielościanów koordynacyjnych. Średnia ważona wszystkich wielościanów koordynacyjnych daje liczbę koordynacyjną, identyczną z tą uzyskaną z całkowania PDF. W wierszu poleceń wpisz:
      speciation_umd.py -f -i -c -a -r 0
      UWAGA: Średnie liczby koordynacyjne w płynach są wartościami ułamkowymi. Wynika to z uśrednionego charakteru koordynacji. Definicja oparta na specjacji daje bardziej intuicyjną i informacyjną reprezentację struktury płynu, w której kwantyfikowane są względne proporcje różnych gatunków, tj. koordynacji.
    2. Uruchom skrypt speciation_umd.py z flagą –r 1, która analizuje graf łączności na wszystkich poziomach głębokości w celu wyznaczenia polimeryzacji. Sieć w grafie atomowym ma określoną głębokość, ponieważ atomy są wiązane dalej z innymi wiązaniami (np. w sekwencjach naprzemiennych kationów i anionów) (Rysunek 4).
  2. Otwórz kolejno dwa pliki: .popul.dat oraz .stat.dat; stanowią one wynik działania skryptu specjacji. Każdy klaster jest zapisany w jednej linii, z określeniem jego wzoru chemicznego, czasu powstania, czasu rozpadu, czasu życia oraz macierzy z listą atomów tworzących ten klaster. Wykreśl czas życia każdego klastra atomowego dla wszystkich gatunków chemicznych znalezionych w symulacji, zgodnie z danymi w pliku .popul.dat (Rysunek 5).
  3. Wykreśl analizę populacji z obfitością każdego gatunku, zgodnie z danymi w pliku .stat.dat. Analiza ta, zarówno w wartościach bezwzględnych, jak i względnych, odpowiada rzeczywistej statystyce wielościanów koordynacyjnych dla przypadku -r 0; w przypadku polimeryzacji (-r 1) należy do niej podejść ostrożnie, gdyż może być konieczne zastosowanie normalizacji względem względnej liczby atomów. Obfitość odpowiada całce z czasów życia. Plik .stat.dat zawiera również rozmiar każdego klastra, czyli liczbę atomów, z których się składa.

4. Obliczanie współczynników dyfuzji

  1. Wyznacz średnie przesunięcie kwadratowe (MSD) atomów jako funkcję czasu, aby uzyskać współczynnik samodyfuzji. Standardowy wzór na MSD to:
    Równanie średniego przesunięcia kwadratowego, MSD(τ), opisujące analizę statystyczną ruchu cząstek.
    gdzie przedmnożniki są renormalizacjami. Za pomocą narzędzia MSD można analizować dynamiczne aspekty płynów na różne sposoby.
    UWAGA: T to całkowity czas symulacji, a Nα to liczba atomów typu α. Czas początkowy t0 jest dowolny i obejmuje pierwszą połowę symulacji. Ninit to liczba czasów początkowych. τ to szerokość przedziału czasowego, dla którego oblicza się MSD; jego maksymalna wartość wynosi połowę 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 wykresu MSD. W tym przypadku i-te okno zaczyna się w czasie t0(i), ale (i+1)-te okno zaczyna się w czasie t0(i) + τ + v, gdzie wartość v jest definiowana przez użytkownika. Podobnie szerokość okna jest zwiększana w dyskretnych krokach zdefiniowanych przez użytkownika, zgodnie z wzorem: τ(i) = τ(i-1) + z. Wartości z („krok poziomy”) i v („krok pionowy”) są dodatnie lub równe zero; wartość domyślna dla obu wynosi 20.
  2. Oblicz MSD, korzystając z serii skryptów msd_umd. Ich wynik jest zapisywany w pliku .msd.dat, w którym MSD każdego typu atomowego, atomu lub klastra jest zapisane w jednej kolumnie jako funkcja czasu.
    1. Oblicz średnie MSD dla każdego typu atomowego. MSD są obliczane dla każdego atomu, a następnie uśredniane dla każdego typu atomowego. Plik wynikowy zawiera jedną kolumnę dla każdego typu atomowego. W wierszu poleceń wpisz:
      msd_umd.py -f -z -v -b
    2. Oblicz MSD dla każdego atomu. MSD są obliczane dla każdego atomu, a następnie uśredniane dla każdego typu atomowego. Plik wynikowy zawiera jedną kolumnę dla każdego atomu w symulacji, a następnie jedną kolumnę dla każdego typu atomowego. Funkcja ta pozwala na identyfikację atomów dyfundujących w dwóch różnych środowiskach, np. w cieczy i gazie lub w dwóch cieczach. W wierszu poleceń wpisz:
      msd_all_umd.py -f -z -v -b
    3. Oblicz MSD dla gatunków chemicznych. Użyj populacji klastrów zidentyfikowanych skryptem do specjacji i zapisanych w pliku .popul.dat. MSD są obliczane dla każdego pojedynczego klastra. Plik wynikowy zawiera jedną kolumnę dla każdego klastra. Aby uniknąć uwzględniania polimerów wielkoskalowych, wprowadź limit wielkości klastra; domyślnie wynosi on 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 do arkuszy kalkulacyjnych (Rycina 6). W reprezentacji log-log MSD w funkcji czasu zidentyfikuj zmianę nachylenia. Oddziel pierwszą część, zazwyczaj krótką, która reprezentuje reżim balistyczny, tj. zachowanie prędkości atomów po zderzeniach. Druga, dłuższa część reprezentuje reżim dyfuzyjny, tj. rozproszenie prędkości atomów po zderzeniach.
  4. Oblicz współczynniki dyfuzji z nachylenia MSD według wzoru:
    Równanie współczynnika dyfuzji D=MSD/2Zt; analiza teoretyczna; badanie dynamiki ruchu.
    gdzie Z to liczba stopni swobody (Z = 2 dla dyfuzji w płaszczyźnie, Z = 3 dla dyfuzji w przestrzeni), a t to krok czasowy.

5. Funkcje korelacji czasowej

  1. Oblicz funkcje korelacji czasowej jako miarę bezwładności układu, korzystając z ogólnego wzoru:
    Wzór na funkcję autokorelacji C(τ)=1/τΣA(t+τ)A(t), równania matematyczne.
    A może reprezentować różne zmienne zależne od czasu, takie jak pozycje atomów, prędkości atomowe, naprężenia, polaryzację itp., z których każda — poprzez relacje Greena-Kubo12,13 — pozwala wyznaczyć różne właściwości fizyczne, czasem po dodatkowej transformacji.
  2. Przeanalizuj prędkości atomowe, aby uzyskać widmo wibracyjne cieczy oraz alternatywny sposób wyznaczenia współczynników samodyfuzji atomowej.
    1. Uruchom skrypt vibr_spectrum_umd.py, aby obliczyć funkcję autokorelacji prędkości-prędkości (VAC) dla każdego typu atomu i wykonać jej szybką transformatę Fouriera. W wierszu poleceń wpisz:
      vibr_spectrum_umd.py -f -t
      gdzie –t to temperatura, która musi zostać zdefiniowana przez użytkownika. Skrypt generuje dwa pliki: plik .vels.scf.dat z funkcją VAC dla każdego typu atomu oraz plik .vibr.dat z widmem wibracyjnym rozbitym na poszczególne gatunki atomowe oraz wartością całkowitą.
    2. Otwórz i odczytaj plik vels.scf.dat. Wykreśl funkcję VAC z pliku vels.scf.dat, korzystając z oprogramowania typu arkusz kalkulacyjny.
    3. Zatrzymaj rzeczywistą część transformaty Fouriera VAC. To ona pozwala uzyskać widmo wibracyjne jako funkcję częstotliwości:
      Wzór na gęstość widmową, powiązany z badaniami wzbudzenia optycznego i absorpcji przejściowej.
      gdzie m to masy atomowe.
    4. Wykreśl widmo wibracyjne z pliku vibr.dat, korzystając z oprogramowania typu arkusz kalkulacyjny (Rysunek 7). Zidentyfikuj skończoną wartość przy ω=0, która odpowiada charakterowi dyfuzyjnemu cieczy, oraz poszczególne piki widma przy skończonej częstotliwości. Określ udział każdego typu atomu w widmie wibracyjnym.
      UWAGA: Rozkład na typy atomowe pokazuje, że różne atomy mają różne wkłady przy ω=0, co odpowiada ich współczynnikom dyfuzji. Ogólny kształt widma jest znacznie gładszy i posiada mniej cech charakterystycznych niż w przypadku odpowiadającego mu ciała stałego.
    5. W powłoce systemowej odczytaj całkę z widma wibracyjnego, która pozwala wyznaczyć współczynniki dyfuzji dla każdego gatunku atomowego.
      UWAGA: Właściwości termodynamiczne można uzyskać poprzez całkowanie widma wibracyjnego, jednak wyniki należy interpretować z ostrożnością ze względu na dwa przybliżenia: całkowanie jest poprawne w ramach przybliżenia kwazi-harmonicznego, które niekoniecznie jest spełnione w wysokich temperaturach; ponadto należy odrzucić część widma o charakterze gazowym, odpowiadającą dyfuzji. Całkowanie powinno być zatem wykonywane tylko po części widma o charakterze sieciowym. Taki podział wymaga jednak zazwyczaj kilku kolejnych etapów post-processingu i obliczeń14, które nie są objęte obecnym pakietem UMD.
  3. Uruchom skrypt viscosity_umd.py, aby przeanalizować autokorelację składowych tensora naprężeń w celu oszacowania lepkości stopu. W wierszu poleceń wpisz:
    viscosity_umd.py -f -i -s -o -l
    UWAGA: Ta funkcja ma charakter badawczy i wszelkie wyniki należy traktować z ostrożnością. Przede wszystkim dokładnie sprawdź zbieżność lepkości w stosunku do długości symulacji.
    1. Wyznacz lepkość cieczy z autokorelacji tensora naprężeń15 zgodnie z wzorem:
      Wzór na lepkość η=V/3kBT∑∫<σij(t+τ)σij(t)>dτ, równanie mechaniki statystycznej.
      gdzie V i T oznaczają odpowiednio objętość i temperaturę, κB jest stałą Boltzmanna, a σij to poza-diagonalna składowa ij tensora naprężeń wyrażona w współrzędnych kartezjańskich.
    2. Zastosuj bardziej odpowiedni dopasowanie, aby uzyskać bardziej wiarygodną estymację lepkości15,16 i uniknąć szumu funkcji autokorelacji tensora naprężeń, który może wynikać z ograniczonego rozmiaru i czasu trwania symulacji. Dla funkcji autokorelacji tensora naprężeń zastosuj następującą postać funkcjonalną15,16, która daje dobre wyniki:
      Mechanika statystyczna, równanie: funkcja korelacji czasowej, wyrazy zaniku, analiza diagramu.
      gdzie A, B, τ1, τ2 oraz ω są parametrami dopasowania. Po całkowaniu wyrażenie na lepkość przyjmuje postać:
      Równanie dynamiki płynów ηkBT/V=Aτ₁+Bτ₂/(1+ω²τ₂²) w analizie wzorów naukowych.

6. Parametry termodynamiczne wynikające z symulacji.

  1. Uruchom program averages.py, aby wyodrębnić wartości średnie oraz rozrzut (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
    gdzie domyślną wartością jest –s 0.
  2. Oblicz błąd statystyczny średniej, stosując metody blokowe.
    UWAGA: Istnieje wiele wariantów tej metody. Zgodnie z pracą Allena i Tildesley2, powszechnym podejściem jest wyciąganie średniej z sekwencji bloków czasowych o coraz większej długości i szacowanie odchylenia standardowego względem średniej arytmetycznej17. Konwergencję można osiągnąć w granicy dużej liczby i dostatecznie długich rozmiarów bloków, gdy próbkowanie jest nieskorelowane. Choć rzeczywista wartość progowa konwergencji zazwyczaj musi zostać wybrana ręcznie.
    1. Zastosuj metodę połowienia18: zaczynając od początkowej próbki danych, w każdym kroku κ zmniejszaj liczbę próbek o połowę, obliczając średnią z każdych dwóch kolejnych próbek odpowiadających poprzedniemu krokowi κ−1:
      Wzór równowagi statycznej; równanie \(S^k_i = \frac{S^{k-1}_{2i} + S^{k-1}_{2i+1}}{2}\).
    2. Uruchom skrypt fullaverages.py, aby przeprowadzić pełną analizę statystyczną, w tym wyznaczyć błąd średniej. W wierszu poleceń wpisz:
      fullaverages.py -s -u
      UWAGA: Skrypt jest zautomatyzowany w taki sposób, że wyszukuje wszystkie pliki .umd.dat w bieżącym katalogu i wykonuje analizę dla każdego z nich. Wartości domyślne to –s 0 –u 0. Dla -u 0 wynik jest minimalny, a dla -u 1 wynik jest pełny, z wydrukowaniem kilku alternatywnych jednostek. Skrypt ten wymaga wsparcia graficznego, ponieważ generuje obraz graficzny służący do sprawdzania konwergencji przy szacowaniu błędu średniej.

Wyniki

Pyrolite jest modelowym wieloskładnikowym stopem krzemianowym (0,5Na2O 2CaO 1,5Al2O3 4FeO 30MgO 24SiO2), który najlepiej przybliża skład całkowity krzemianowej Ziemi — średnią geochemiczną naszej planety z pominięciem jej żelaznego jądra19. Wczesna Ziemia była zdominowana przez serię wielkoskalowych zdarzeń topnienia20; ostatnie z nich mogło objąć całą planetę po jej kondensacji dla protolunarno-ziemskiego dysku21. Pyrolite stanowi najlepsze przybliżenie składu takich oceanów magmy w skali planetarnej. W związku z tym szczegółowo zbadaliśmy właściwości fizyczne stopu pyrolitu w zakresie temperatur 3 000‒5 000 K i w zakresie ciśnień 0‒150 GPa za pomocą symulacji dynamiki molekularnej ab initio w implementacji VASP. Te warunki termodynamiczne w pełni charakteryzują najbardziej ekstremalne warunki panujące w ziemskich oceanach magmy. Nasza praca jest doskonałym przykładem skutecznego wykorzystania pakietu UMD do kompleksowej, pogłębionej analizy stopów22. Obliczyliśmy rozkłady i średnie długości wiązań, prześledziliśmy zmiany w koordynacji kation-tlen oraz porównaliśmy nasze wyniki z wcześniejszymi badaniami eksperymentalnymi i obliczeniowymi dotyczącymi krzemianów amorficznych o różnych składach. Nasza szczegółowa analiza pomogła rozłożyć standardowe liczby koordynacyjne na ich podstawowe składniki, wskazać obecność egzotycznych wielościanów koordynacyjnych w stopie oraz wyznaczyć czasy życia dla wszystkich wielościanów koordynacyjnych. Wskazano w niej również znaczenie próbkowania w symulacjach, zarówno pod kątem długości trajektorii, jak i liczby atomów obecnych w modelowanym systemie. Jeśli chodzi o postprocessing, analiza w pakiecie UMD jest niezależna od tych czynników, jednak należy je brać pod uwagę podczas interpretacji wyników dostarczonych przez ten pakiet. Poniżej przedstawiamy kilka przykładów wykorzystania pakietu UMD do wyłonienia szeregu charakterystycznych cech stopów w odniesieniu do stopionego pyrolitu.

Funkcja rozkładu par Si-O uzyskana za pomocą skryptu gofrs_umd.py pokazuje, że promień pierwszej sfery koordynacyjnej, który stanowi pierwsze minimum funkcji g(r), wynosi około 2.5 angstromów przy T = 3000 K i P = 4.6 GPa. Maksimum g(r) znajduje się przy 1.635 Å — stanowi to najlepsze przybliżenie długości wiązania. Długi ogon wynika z temperatury. Przyjmując tę granicę jako odległość wiązania Si-O, analiza specjacji wykazuje, że w stopie dominują jednostki SiO4, które mogą utrzymywać się przez maksymalnie kilka pikosekund (Rysunek 5). Istotna część stopu wykazuje częściową polimeryzację, o czym świadczy obecność dimerów, takich jak Si2O7, oraz trimerów, takich jak jednostki Si3Ox. Ich czas życia jest rzędu pikosekund. Polimery wyższego rzędu mają znacznie krótszy czas życia.

Wykres rozkładu czasu życia klastrów tlenku krzemu, analiza składu SiOx, wyniki spektroskopii zaniku.
Rycina 5: Czas życia chemicznych form Si-O.
Speciacja została zidentyfikowana w stopie wieloskładnikowym przy 4,6 GPa i 3000 K. Etykiety oznaczają monomery SiO3, SiO4 i SiO5 oraz różne polimery SixOy. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Różne wartości kroków pionowych i poziomych, definiowane przez powyższe flagi –z i –v, dają różne próbkowania MSD (Rysunek 6). Nawet duże wartości z i v są wystarczające do zdefiniowania nachyleń, a tym samym współczynników dyfuzji różnych atomów. Zysk czasu podczas przetwarzania danych jest znaczący przy przejściu na duże wartości z i v. MSD stanowi bardzo silne kryterium walidacji jakości symulacji. Jeśli część dyfuzyjna MSD nie jest wystarczająco długa, jest to sygnał, że symulacja jest zbyt krótka i w sensie statystycznym nie osiągnięto stanu płynnego. Minimalne wymagania dla części dyfuzyjnej MSD w dużym stopniu zależą od systemu. Można wymagać, aby wszystkie atomy zmieniły swoje miejsce w strukturze stopu co najmniej raz, aby można go było uznać za płyn10. Doskonałym przykładem z zastosowaniami w naukach planetarnych są złożone stopy krzemianów pod wysokim ciśnieniem, blisko lub nawet poniżej ich linii liquidus11. Atomy Si, będące głównymi kationami tworzącymi sieć, zmieniają miejsca po więcej niż dwóch tuzinach pikosekund. Symulacje krótsze niż ten próg wiązałyby się ze znacznym niedopróbkowaniem możliwej przestrzeni konfiguracyjnej. Jednakże, ponieważ aniony koordynujące, a mianowicie atomy O, poruszają się szybciej niż centralne atomy Si, mogą one zrekompensować część powolnej mobilności Si. W ten sposób cały system mógłby w rzeczywistości objąć lepsze próbkowanie przestrzeni konfiguracyjnej, niż można by założyć, opierając się wyłącznie na przemieszczeniach Si.

Wykres średniego przesunięcia kwadratowego, czas vs MSD dla Mg, O, Si; analiza obliczeniowa, badanie dynamiczne.
Rycina 6: Średnie przesunięcia kwadratowe (MSD).
Przedstawiono MSD dla kilku typów atomów w wieloskładnikowym stopie krzemianowym. Próbkowanie z różnymi krokami poziomymi i pionowymi, z i v, daje spójne wyniki. Pełne koła: -z 50 –v 50. Puste koła: -z 250 –v 500. Kliknij tutaj, aby zobaczyć powiększoną wersję tej ryciny.

Ostatecznie funkcje atomowe VAC dają widmo wibracyjne stopu. Rysunek 7 przedstawia widmo przy tych samych warunkach ciśnienia i temperatury co powyżej. Przedstawiono wkłady atomów Mg, Si i O, a także wartość całkowitą. Przy częstotliwości zero występuje skończona wartość widma, która odpowiada charakterowi dyfuzyjnemu stopu. Wyznaczenie właściwości termodynamicznych z widma wibracyjnego wymaga usunięcia tego gazopodobnego charakteru dyfuzyjnego przy zerze, a także odpowiedniego uwzględnienia jego zaniku przy wyższych częstotliwościach.

Wykres spektroskopii Ramanowskiej przedstawiający intensywność w funkcji częstotliwości dla Mg, Si, O oraz wszystkich pierwiastków.
Rycina 7: Widmo wibracyjne stopu pyrolitu.
Część rzeczywista transformaty Fouriera funkcji autokorelacji prędkości atomów pozwala na wyznaczenie widma wibracyjnego. W tym przypadku widmo obliczono dla wieloskładnikowego stopu krzemianowego. Ciecze charakteryzują się niezerowym, gazopodobnym charakterem dyfuzyjnym przy częstotliwości zero. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

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.

Oświadczenia

Autorzy nie mają nic do ujawnienia.

Podziękowania

Ta praca była wspierana przez Europejską Radę ds. Badań Naukowych (ERC) w ramach programu Unii Europejskiej Horyzont 2020 badań i innowacji (numer umowy o grant 681818 IMPACT dla RC), przez Dyrekcję ds. Fizyki Ekstremalnej i Chemii Obserwatorium Głębokiego Węgla oraz przez Norweską Radę ds. Badań poprzez program finansowania Centrów Doskonałości, numer projektu 223272. Potwierdzamy dostęp do superkomputerów GENCI w ramach serii grantów obliczeniowych eDARI stl2816, do superkomputera Irene AMD w ramach projektu PRACE RA4947 oraz superkomputera Fram za pośrednictwem UNINETT Sigma2 NN9697K. FS otrzymał wsparcie w ramach działania "Maria Skłodowska-Curie" (umowa o grant ABISSE nr 750901).

Materiały

Lista materiałów użytych w tym artykule
NazwaFirmaNumer katalogowyKomentarze
getopt biblioteka
kodzie
źródłowym biblioteka matplotlib biblioteka
systemu operacyjnegootwartym kodzie źródłowym
Python Software Foundationw wersji 2 i 3Biblioteka
losowatypu open source Biblioteka
typu open source Biblioteka
kodzie źródłowym Biblioteka
source Biblioteka
sys
glob o otwartym numpy o otwartym kodzie źródłowym biblioteka o otwartym kodzie źródłowym oprogramowanie Python o RE Scipy o otwartym podprocesów typu open typu open source

Bibliografia

  1. Frenkel, D., Smit, B. Understanding Molecular Simulation. From Algorithms to Applications. , Elsevier. (2001).
  2. Allen, M. P., Tildesley, D. J., Allen, T. Computer Simulation of Liquids. , Oxford University Press. (1989).
  3. Zepeda-Ruiz, L. A., Stukowski, A., Oppelstrup, T., Bulatov, V. V. Probing the limits of metal plasticity with molecular-dynamics simulations. Nature Publishing Group. 550 (7677), 492-495 (2017).
  4. Humphrey, W., Dalke, A., Schulten, K. VMD: Visual molecular dynamics. Journal of Molecular Graphics & Modeling. 14 (1), 33-38 (1996).
  5. Momma, K., Izumi, F. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. Journal of Applied Crystallography. 44 (6), 1272-1276 (2011).
  6. Brehm, M., Kirchner, B. TRAVIS - A free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories. Journal of Chemical Information and Modeling. 51 (8), 2007-2023 (2011).
  7. Stixrude, L. Visualization-based analysis of structural and dynamical properties of simulated hydrous silicate melt. Physics and Chemistry of Minerals. 37 (2), 103-117 (2009).
  8. Kresse, G., Hafner, J. Ab initio Molecular-Dynamics for Liquid-Metals. Physical Review B. 47 (1), 558-561 (1993).
  9. Gygi, F. Architecture of Qbox: A scalable first-principles molecular dynamics code. IBM Journal of Research and Development. 52 (1-2), 137-144 (2008).
  10. Harvey, J. P., Asimow, P. D. Current limitations of molecular dynamic simulations as probes of thermo-physical behavior of silicate melts. American Mineralogist. 100 (8-9), 1866-1882 (2015).
  11. Caracas, R., Hirose, K., Nomura, R., Ballmer, M. D. Melt-crystal density crossover in a deep magma ocean. Earth and Planetary Science Letters. 516, 202-211 (2019).
  12. Green, M. S. Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids. The Journal of Chemical Physics. 22 (3), 398-413 (1954).
  13. Kubo, R. Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems. Journal of the Physical Society of Japan. 12 (6), 570-586 (1957).
  14. Lin, S. T., Blanco, M., Goddard, W. A. The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids. The Journal of Chemical Physics. 119 (22), 11792-11805 (2003).
  15. Meyer, E. R., Kress, J. D., Collins, L. A., Ticknor, C. Effect of correlation on viscosity and diffusion in molecular-dynamics simulations. Physical Review E. 90 (4), 1198-1212 (2014).
  16. Soubiran, F., Militzer, B., Driver, K. P., Zhang, S. Properties of hydrogen, helium, and silicon dioxide mixtures in giant planet interiors. Physics of Plasmas. 24 (4), 041401-041407 (2017).
  17. Flyvbjerg, H., Petersen, H. G. Error estimates on averages of correlated data. The Journal of Chemical Physics. 91 (1), 461-466 (1989).
  18. Tuckerman, M. E. Statistical mechanics: theory and molecular simulation. , Oxford University Press. (2010).
  19. McDonough, W. F., Sun, S. S. The composition of the Earth. Chemical Geology. 120, 223-253 (1995).
  20. Elkins-Taton, L. T. Magma oceans in the inner solar system. Annual Review of Earth and Planetary Sciences. 40, 113-139 (2012).
  21. Lock, S. J., et al. The origin of the Moon within a terrestrial synestia. J. Geophysical Research: Planets. 123, 910-951 (2018).
  22. Solomatova, N. V., Caracas, R. Pressure-induced coordination changes in a pyrolitic silicate melt from ab initio molecular dynamics simulations. Journal of Geophysical Research: Solid Earth. 124, 11232-11250 (2019).

Przedruki i uprawnienia

Tagi

Funkcja rozkładu radialnegospecjacja chemicznaśredni kwadrat przemieszczeniawspółczynniki dyfuzjiwidmo wibracyjneobliczanie lepkościanaliza strukturalnawłaściwości transportowe
Analiza stopów i płynów z symulacji dynamiki molekularnej ab initio za pomocą pakietu UMD
Skontaktuj się z nami
Poleć bibliotece

Prywatność

Warunki korzystania z usługi

Zasady

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