Summary

Analyse von Schmelzen und Flüssigkeiten aus Ab Initio Molekulardynamiksimulationen mit dem UMD-Paket

Published: September 17, 2021
doi:

Summary

Schmelzen und Flüssigkeiten sind allgegenwärtige Vektoren des Massentransports in natürlichen Systemen. Wir haben ein Open-Source-Paket entwickelt, um ab initio molekulardynamische Simulationen solcher Systeme zu analysieren. Wir berechnen strukturelle (Bindung, Clusterisierung, chemische Speziation), Transport (Diffusion, Viskosität) und thermodynamische Eigenschaften (Schwingungsspektrum).

Abstract

Wir haben ein Python-basiertes Open-Source-Paket entwickelt, um die Ergebnisse von ab initio molekulardynamischen Simulationen von Flüssigkeiten zu analysieren. Das Paket eignet sich am besten für Anwendungen an natürlichen Systemen wie Silikat- und Oxidschmelzen, Flüssigkeiten auf Wasserbasis und verschiedenen überkritischen Flüssigkeiten. Das Paket ist eine Sammlung von Python-Skripten, die zwei Hauptbibliotheken enthalten, die sich mit Dateiformaten und mit Kristallographie befassen. Alle Skripts werden über die Befehlszeile ausgeführt. Wir schlagen ein vereinfachtes Format vor, um die atomaren Trajektorien und relevanten thermodynamischen Informationen der Simulationen zu speichern, die in UMD-Dateien gespeichert sind, was für Universal Molecular Dynamics steht. Das UMD-Paket ermöglicht die Berechnung einer Reihe von strukturellen, Transport- und thermodynamischen Eigenschaften. Ausgehend von der Paarverteilungsfunktion definiert es Bindungslängen, baut eine interatomare Konnektivitätsmatrix auf und bestimmt schließlich die chemische Speziation. Die Bestimmung der Lebensdauer der chemischen Spezies ermöglicht die Durchführung einer vollständigen statistischen Analyse. Dann berechnen dedizierte Skripte die mittleren quadratischen Verschiebungen für die Atome sowie für die chemische Spezies. Die implementierte Selbstkorrelationsanalyse der Atomgeschwindigkeiten ergibt die Diffusionskoeffizienten und das Schwingungsspektrum. Die gleiche Analyse, die auf die Spannungen angewendet wird, ergibt die Viskosität. Das Paket ist über die GitHub-Website und über eine eigene Seite des ERC IMPACT-Projekts als Open-Access-Paket verfügbar.

Introduction

Flüssigkeiten und Schmelzen sind aktive chemische und physikalische Transportvektoren in natürlichen Umgebungen. Die erhöhten Raten der atomaren Diffusion begünstigen den chemischen Austausch und die Reaktionen, die niedrige Viskosität in Verbindung mit unterschiedlichem Auftrieb begünstigt einen großen Stoffaustausch und die Kristallschmelzedichtebeziehungen begünstigen die Schichtung innerhalb planetarer Körper. Das Fehlen eines periodischen Gitters, typische hohe Temperaturen, die erforderlich sind, um den geschmolzenen Zustand zu erreichen, und die Schwierigkeit beim Abschrecken machen die experimentelle Bestimmung einer Reihe offensichtlicher Eigenschaften wie Dichte, Diffusion und Viskosität extrem schwierig. Diese Schwierigkeiten machen alternative Berechnungsmethoden zu starken und nützlichen Werkzeugen für die Untersuchung dieser Materialklasse.

Mit dem Aufkommen der Rechenleistung und der Verfügbarkeit von Supercomputern werden derzeit zwei wichtige numerische atomistische Simulationstechniken eingesetzt, um den dynamischen Zustand eines nichtkristallinen atomistischen Systems zu untersuchen, Monte Carlo1 und Molekulardynamik (MD)1,2. In Monte-Carlo-Simulationen wird der Konfigurationsraum zufällig abgetastet; Monte-Carlo-Methoden zeigen eine lineare Skalierung in Parallelisierung, wenn alle Stichprobenbeobachtungen unabhängig voneinander sind. Die Qualität der Ergebnisse hängt von der Qualität des Zufallszahlengenerators und der Repräsentativität der Stichprobe ab. Monte-Carlo-Methoden zeigen lineare Skalierung in Parallelisierung, wenn die Stichprobe unabhängig voneinander ist. In der Molekulardynamik (MD) wird der Konfigurationsraum durch zeitabhängige atomare Trajektorien abgetastet. Ausgehend von einer gegebenen Konfiguration werden die atomaren Trajektorien durch Integration der Newtonschen Bewegungsgleichungen berechnet. Die interatomaren Kräfte können mit Hilfe von modellierten interatomaren Potentialen (in klassischer MD) oder mit Methoden der ersten Prinzipien (in ab initio oder first-principles, MD) berechnet werden. Die Qualität der Ergebnisse hängt von der Länge der Trajektorie und ihrer Fähigkeit ab, nicht von lokalen Minima angezogen zu werden.

Molekulardynamik-Simulationen enthalten eine Fülle von Informationen, die alle mit dem dynamischen Verhalten des Systems zusammenhängen. Thermodynamische Durchschnittseigenschaften wie innere Energie, Temperatur und Druck sind eher Standard zu berechnen. Sie können aus den Ausgabedateien der Simulationen extrahiert und gemittelt werden, während Größen, die sich direkt auf die Bewegung der Atome sowie auf ihre gegenseitige Beziehung beziehen, nach Extraktion der atomaren Positionen und Geschwindigkeiten berechnet werden müssen.

Folglich wurde viel Mühe darauf verwendet, die Ergebnisse zu visualisieren, und verschiedene Pakete sind heute auf verschiedenen Plattformen verfügbar, Open Source oder nicht [Ovito3, VMD4, Vesta5, Travis6 usw.]. Alle diese Visualisierungswerkzeuge arbeiten effizient mit interatomaren Abständen und ermöglichen so die effiziente Berechnung von Paarverteilungsfunktionen und Diffusionskoeffizienten. Verschiedene Gruppen, die groß angelegte molekulardynamische Simulationen durchführen, verfügen über proprietäre Software, um verschiedene andere Eigenschaften zu analysieren, die sich aus den Simulationen ergeben, manchmal in Shareware oder anderen Formen des eingeschränkten Zugangs zur Community und manchmal begrenzt in Umfang und Verwendung auf einige bestimmte Pakete. Ausgeklügelte Algorithmen zur Extraktion von Informationen über interatomare Bindungen, geometrische Muster und Thermodynamik werden entwickelt und in einigen dieser Pakete implementiert3,4,5,6,7 usw.

Hier schlagen wir das UMD-Paket vor – ein in Python geschriebenes Open-Source-Paket, um die Ausgabe molekulardynamischer Simulationen zu analysieren. Das UMD-Paket ermöglicht die Berechnung eines breiten Spektrums struktureller, dynamischer und thermodynamischer Eigenschaften (Abbildung 1). Das Paket ist über die GitHub-Website (https://github.com/rcaracas/UMD_package) und über eine eigene Seite (http://moonimpact.eu/umd-package/) des ERC IMPACT-Projekts als Open-Access-Paket verfügbar.

Um es universell und einfacher zu handhaben zu machen, besteht unser Ansatz darin, zunächst alle Informationen über den thermodynamischen Zustand und die atomaren Trajektorien aus der Ausgabedatei des eigentlichen Molekulardynamiklaufs zu extrahieren. Diese Informationen werden in einer dedizierten Datei gespeichert, deren Format unabhängig vom ursprünglichen MD-Paket ist, in dem die Simulation ausgeführt wurde. Wir nennen diese Dateien “umd”-Dateien, was für Universal Molecular Dynamics steht. Auf diese Weise kann unser UMD-Paket von jeder ab initio-Gruppe mit jeder Software einfach und mit minimalem Anpassungsaufwand verwendet werden. Die einzige Voraussetzung für die Verwendung des vorliegenden Pakets besteht darin, den entsprechenden Parser aus der Ausgabe der jeweiligen MD-Software in das umd-Dateiformat zu schreiben, falls dieses noch nicht vorhanden ist. Vorerst stellen wir solche Parser für die Pakete VASP8 und QBox9 zur Verfügung.

Figure 1
Abbildung 1: Flussdiagramm der UMD-Bibliothek.
Physische Eigenschaften sind blau und die wichtigsten Python-Skripte und ihre Optionen sind rot. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Die umd-Dateien sind ASCII-Dateien; Typische Erweiterung ist “umd.dat”, aber nicht obligatorisch. Alle Analysekomponenten können ASCII-Dateien im umd-Format lesen, unabhängig von der tatsächlichen Namenserweiterung. Einige der automatischen Skripte, die entwickelt wurden, um schnelle umfangreiche Statistiken über mehrere Simulationen hinweg durchzuführen, suchen jedoch speziell nach Dateien mit der Erweiterung umd.dat. Jede physische Eigenschaft wird in einer Zeile ausgedrückt. Jede Zeile beginnt mit einem Schlüsselwort. Auf diese Weise ist das Format sehr anpassungsfähig und ermöglicht das Hinzufügen neuer Eigenschaften zur umd-Datei, während gleichzeitig die Lesbarkeit in allen Versionen erhalten bleibt. Die ersten 30 Zeilen der umd-Datei der Simulation von Pyrolit bei 4,6 GPa und 3000 K, die unten in der Diskussion verwendet werden, sind in Abbildung 2 dargestellt.

Figure 2
Abbildung 2: Der Anfang der umd-Datei, die die Simulation von flüssigem Pyrolit bei 4,6 GPa und 3000 K beschreibt.
Auf den Header folgt die Beschreibung jedes Snapshots. Jede Eigenschaft wird in eine Zeile geschrieben, die den Namen der physischen Eigenschaft, die Werte und die Einheiten enthält, die alle durch Leerzeichen getrennt sind. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Alle umd-Dateien enthalten einen Header, der den Inhalt der Simulationszelle beschreibt: die Anzahl der Atome, Elektronen und Atomtypen sowie Details für jedes Atom, wie seinen Typ, sein chemisches Symbol, die Anzahl der Valenzelektronen und seine Masse. Eine leere Zeile markiert das Ende des Headers und trennt ihn vom Hauptteil der umd-Datei.

Dann wird jeder Schritt der Simulation detailliert beschrieben. Zuerst werden die momentanen thermodynamischen Parameter jeweils auf einer anderen Linie angegeben, wobei (i) der Name des Parameters wie Energie, Spannungen, äquivalenter hydrostatischer Druck, Dichte, Volumen, Gitterparameter usw., (ii) seine Werte und (iii) seine Einheiten angegeben werden. Eine Tabelle, die die Atome beschreibt, kommt als nächstes. Eine Kopfzeile gibt die verschiedenen Maße wie kartesische Positionen, Geschwindigkeiten, Ladungen usw. und ihre Einheiten an. Dann wird jedes Atom auf einer Zeile detailliert beschrieben. Durch Dreiergruppen, die den drei x-, y-, z-Achsen entsprechen, sind die Einträge: die reduzierten Positionen, die kartesischen Positionen, die in die Simulationszelle gefaltet sind, die kartesischen Positionen (die die Tatsache berücksichtigen, dass Atome während einer Simulation mehrere Einheitszellen durchqueren können), die Atomgeschwindigkeiten und die Atomkräfte. Die letzten beiden Einträge sind Skalare: Ladung und magnetisches Moment.

Zwei Hauptbibliotheken gewährleisten das ordnungsgemäße Funktionieren des gesamten Pakets. Die umd_process.py Bibliothek befasst sich mit den umd-Dateien, wie Lesen und Drucken. Die crystallography.py Bibliothek befasst sich mit allen Informationen, die sich auf die tatsächliche atomare Struktur beziehen. Die zugrunde liegende Philosophie der crystallography.py Bibliothek besteht darin, das Gitter als Vektorraum zu behandeln. Die Elementarzellenparameter stellen zusammen mit ihrer Orientierung die Basisvektoren dar. Der “Raum” hat eine Reihe von skalaren Attributen (spezifisches Volumen, Dichte, Temperatur und spezifische Anzahl von Atomen), thermodynamische Eigenschaften (innere Energie, Druck, Wärmekapazität usw.) und eine Reihe von Tensorialeigenschaften (Spannung und Elastizität). Atome bevölkern diesen Raum. Die Klasse “Lattice” definiert dieses Ensemble, zusammen mit verschiedenen kurzen Berechnungen, wie spezifisches Volumen, Dichte, Erhalten des reziproken Gitters aus dem direkten usw. Die Klasse “Atoms” definiert die Atome. Sie zeichnen sich durch eine Reihe skalarer Eigenschaften (Name, Symbol, Masse, Anzahl der Elektronen usw.) und eine Reihe von vektoriellen Eigenschaften (die Position im Raum, entweder relativ zu der in der Lattice-Klasse beschriebenen vektoriellen Basis oder relativ zu universellen kartesischen Koordinaten, Geschwindigkeiten, Kräften usw.) aus. Abgesehen von diesen beiden Klassen enthält die crystallography.py Bibliothek eine Reihe von Funktionen zum Ausführen einer Vielzahl von Tests und Berechnungen, z. B. atomare Abstände oder Zellmultiplikation. Das Periodensystem der Elemente ist auch als Wörterbuch enthalten.

Die verschiedenen Komponenten des umd-Pakets schreiben mehrere Ausgabedateien. In der Regel handelt es sich bei allen um ASCII-Dateien, alle ihre Einträge sind durch Tabulatoren getrennt und sie sind so selbsterklärend wie möglich gestaltet. Zum Beispiel zeigen sie immer klar die physische Eigenschaft und ihre Einheiten an. Die umd.dat-Dateien entsprechen vollständig dieser Regel.

Protocol

1. Analyse der Molekulardynamik verläuft HINWEIS: Das Paket ist über die GitHub-Website (https://github.com/rcaracas/UMD_package) und über eine eigene Seite (http://moonimpact.eu/umd-package/) des ERC IMPACT-Projekts als Open-Access-Paket verfügbar. Extrahieren Sie jeden spezifischen Satz physischer Eigenschaften mit einem oder mehreren dedizierten Python-Skripts aus dem Paket. Führen Sie alle Skripts über die Befehlszeile aus. Sie alle verwenden eine Reihe von Flags, die von einem Skript zum anderen so konsistent wie möglich sind. Die Flags, ihre Bedeutung und die Standardwerte sind alle in Tabelle 1 zusammengefasst. Flagge Bedeutung Skript, das es verwendet Standardwert -h Kurze Hilfe alle -f UMD-Dateiname alle -i Zu verwerfende Thermalisierungsschritte alle 0 -i Eingabedatei mit den interatomaren Bindungen Speziation Bindungen.Eingabe -s Abtastung der Frequenz msd, Artbildung 1 (jeder Schritt wird berücksichtigt) -a Liste der Atome oder Anionen Speziation -c Liste der Kationen Speziation -l Haftlänge Speziation 2 -t Temperatur Vibrationen, Rheologie -v Diskretisierung der Breite des Abtastfensters der Trajektorie für die Mittelwert-Quadrat-Verschiebungsanalyse Msd 20 -z Diskretisierung des Beginns des Abtastfensters der Trajektorie für die Mittelwert-Quadrat-Verschiebungsanalyse Msd 20 Tabelle 1: Die häufigsten im UMD-Paket verwendeten Flags und ihre häufigste Bedeutung. Beginnen Sie mit der Umwandlung der Ausgabe der MD-Simulation, die in einem First-Principles-Code wie VASP8 oder QBox9 ausgeführt wurde, in eine UMD-Datei. Wenn die MD-Simulationen in VASP durchgeführt wurden, geben Sie in der Befehlszeile Folgendes ein:VaspParser.py -f -i wobei das Flag –f den Namen der VASP OUTCAR-Datei und das –i die Thermalisierungslänge definiert.HINWEIS: Der erste Schritt, definiert durch –i, ermöglicht es, die ersten Schritte der Simulationen, die die Thermalisierung darstellen, zu verwerfen. In einem typischen molekulardynamischen Lauf stellt der erste Teil der Berechnung die Thermalisierung dar, d.h. die Zeit, die das System benötigt, um eine Gaußsche Verteilung der Temperatur zu beschreiben, und für das gesamte System, um Schwankungen der Temperatur, des Drucks, der Energie usw. um Gleichgewichtswerte zu zeigen. Dieser Thermalisierungsteil der Simulation sollte bei der Analyse der statistischen Eigenschaften des Fluids nicht berücksichtigt werden. Transformieren Sie die . umd-Dateien in . xyz-Dateien , um die Visualisierung auf verschiedenen anderen Paketen wie VMD4 oder Vesta5 zu erleichtern. Geben Sie in der Befehlszeile Folgendes ein:umd2xyz.py -f -i -s wobei –f den Namen von definiert. umd-Datei definiert –i den zu verwerfenden Thermalisierungszeitraum und –s die Häufigkeit der Abtastung der trajektorie, die in der gespeichert ist. umd-Datei . Standardwerte sind –i 0 –s 1, d. h. unter Berücksichtigung aller Schritte der Simulation, ohne dass welche verworfen werden. Wandeln Sie die umd-Datei mithilfe des umd2poscar.py-Skripts in POSCAR-Dateien vom Typ VASP um. Schnappschüsse der Simulationen können mit einer vordefinierten Häufigkeit ausgewählt werden. Geben Sie in der Befehlszeile Folgendes ein:umd2poscar.py -f -i -l -s wobei –l den letzten Schritt darstellt, der in eine POSCAR-Datei umgewandelt werden soll. Die Standardwerte sind -i 0 -l 10000000 -s 1. Dieser Wert von -l ist groß genug, um eine typische gesamte Flugbahn abzudecken. 2. Führen Sie die Strukturanalyse durch Führen Sie das skript gofrs_umd.py aus, um die Paarverteilungsfunktion (PDF) gᴀʙ(r) für alle Paare der Atomtypen A und B zu berechnen (Abbildung 3). Die Ausgabe wird in einer ASCII-Datei geschrieben, tabulatorgetrennt, mit der Erweiterung gofrs.dat. Geben Sie in der Befehlszeile Folgendes ein:gofrs_umd.py -f -s -d -i HINWEIS: Die Standardwerte sind Sampling_Frequency (die Häufigkeit für die Abtastung der Leitkurve) = 1 Schritt; DiskretisierungIntervall (zur Darstellung des g(r)) = 0,01 Å; InitialStep (Anzahl der Schritte am Anfang der Leitkurve, die verworfen werden) = 0. Das radiale PDF gᴀʙ(r) ist die durchschnittliche Anzahl von Atomen vom Typ B in einem Abstand d_ᴀʙ innerhalb einer kugelförmigen Schale mit Radius r und Dicke dr, zentriert auf den Atomen des Typs A (Abbildung 3):mit ρ die Atomdichte, NA und NB die Anzahl der Atome vom Typ A und B und δ(r−rᴀʙ) die Deltafunktion, die gleich 1 ist, wenn die Atome A und B in einem Abstand zwischen r und r+dr liegen. Die Abszisse des ersten Maximums von gᴀʙ(r) ergibt die höchste Wahrscheinlichkeitsbindungslänge zwischen den Atomen vom Typ A und B, die einem durchschnittlichen Bindungsabstand am nächsten kommt, den wir bestimmen können. Das erste Minimum grenzt die Ausdehnung des ersten Koordinationsbereichs ab. Daher ergibt das Integral über dem PDF bis zum ersten Minimum die durchschnittliche Koordinationsnummer. Die Summe der Fourier-Transformationen des gᴀʙ(r) für alle Paare der Atomtypen A und B ergibt das Beugungsmuster der Flüssigkeit, wie es experimentell mit einem Diffraktometer erhalten wurde. Da jedoch in der Realität die Koordinationskugeln hoher Ordnung oft im gᴀʙ(r) fehlen, kann das Beugungsmuster nicht in seiner Gesamtheit erhalten werden. Abbildung 3: Bestimmung der Paarverteilungsfunktion.(a) Für jedes Atom einer Spezies (z. B. rot) werden alle Atome der koordinierenden Spezies (z. B. grau und/oder rot) als Funktion der Entfernung gezählt. (b) Der resultierende Entfernungsverteilungsgraph für jeden Schnappschuss, der in diesem Stadium nur eine Sammlung von Deltafunktionen ist, wird dann über alle Atome und alle Schnappschüsse gemittelt und mit der idealen Gasverteilung gewichtet, um (c) die paarweise Verteilungsfunktion zu erzeugen, die kontinuierlich ist. Das erste Minimum des g(r) ist der Radius der ersten Koordinationskugel, der später in der Speziationsanalyse verwendet wird. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen. Extrahieren Sie die mittleren interatomaren Bindungsabstände als Radien der ersten Koordinationskugeln. Identifizieren Sie dazu die Position des ersten Maximums der gᴀʙ(r)-Funktionen: Plotten Sie die gofrs.dat-Datei in einer Tabellenkalkulationsanwendung und suchen Sie nach den Maxima und Minima für jedes Atompaar. Identifizieren Sie den Radius der ersten Koordinationskugel als erstes Minimum der PDF-Datei gᴀʙ(r) mithilfe einer Tabellenkalkulationssoftware. Dies ist die Grundlage für die gesamte Strukturanalyse des Fluids; Das PDF liefert den durchschnittlichen Bindungsstatus der Atome in der Flüssigkeit. Extrahieren Sie die Abstände der ersten Minima, d.h. der Abszisse, und schreiben Sie sie in eine separate Datei, die beispielsweise bonds.input heißt. Alternativ können Sie eines der analyze_gofr Skripte des UMD-Pakets ausführen, um die Maxima und die Minima der Funktionen gᴀʙ(r) zu identifizieren. Geben Sie in der Befehlszeile Folgendes ein:analyze_gofr_semi_automatic.py Klicken Sie auf die Position des Maximums und des Minimums der Funktion gᴀʙ(r), die in der vom Programm geöffneten Grafik angezeigt wird. Das Skript scannt automatisch den aktuellen Ordner, identifiziert alle gofrs.dat Dateien und führt die Analyse für jede von ihnen durch. Klicken Sie jedes Mal erneut auf das Maximum und das Minimum im Fenster, wenn das Skript eine fundierte Ersteinschätzung benötigt. Öffnen und sehen Sie sich die automatisch generierte Datei bonds.input an, die die interatomaren Bindungsabstände enthält. 3. Führen Sie die Speziationsanalyse durch Berechnen Sie die Topologie der Bindung zwischen den Atomen unter Verwendung des Konzepts der Konnektivität innerhalb der Graphentheorie: Die Atome sind die Knoten und die interatomaren Bindungen sind die Pfade. Das speciation_umd.py Skript benötigt die interatomaren Bindungsabstände, die in der Datei bonds.input definiert sind .HINWEIS: Die Konnektivitätsmatrix wird in jedem Zeitschritt aufgebaut: Zwei Atome, die in einem Abstand liegen, der kleiner ist als der Radius ihrer entsprechenden ersten Koordinationskugel, gilt als gebunden, d.h. verbunden. Verschiedene atomare Netzwerke werden aufgebaut, indem die Atome als Knoten in einem Graphen behandelt werden, deren Verbindungen durch dieses geometrische Kriterium definiert werden. Diese Netzwerke sind die atomaren Spezies, und ihr Ensemble definiert die atomare Artbildung in dieser bestimmten Flüssigkeit (Abbildung 4). Abbildung 4: Identifizierung der atomaren Cluster.Die Koordinationspolyeder werden über interatomare Abstände definiert. Alle Atome in einem Abstand, der kleiner als ein bestimmter Radius ist, gelten als gebunden. Hier entspricht der Schwellenwert der ersten Koordinationskugel (den hellroten Kreisen), die in Abbildung 1 definiert ist. Polymerisation und damit chemische Spezies werden aus den Netzwerken der gebundenen Atome gewonnen. Beachten Sie den zentralen Red1Grey2-Cluster, der von den anderen Atomen isoliert ist, die ein unendliches Polymer bilden. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen. Führen Sie das Speziationsskript aus, um die Konnektivitätsmatrix und die Koordinationspolyeder oder die Polymerisation zu erhalten. Geben Sie in der Befehlszeile Folgendes ein:speciation_umd.py -f -s -i -l -c -a -m -r wobei das -i-Flag die Datei mit den interatomaren Bindungsabständen angibt, die beispielsweise im vorherigen Schritt erzeugt wurde. Alternativ können Sie das Skript mit einer einzigen Länge für alle Bindungen ausführen, die durch das Flag -l definiert sind.HINWEIS: Das Flag -c gibt die zentralen Atome und das Flag -a die Liganden an. Sowohl Zentralatome als auch Liganden können von unterschiedlichem Typ sein; in diesem Fall müssen sie durch Kommas getrennt werden. Die -m-Flagge gibt die Mindestdauer an, die eine Art leben muss, um bei der Analyse berücksichtigt zu werden. Standardmäßig ist diese Mindestzeit Null, wobei alle Vorkommen letztendlich gezählt werden. Führen Sie das skript speciation_umd.py mit dem Flag –r 0 aus, das den Konnektivitätsgraphen auf der ersten Ebene abtastet, um die Koordinationspolyeder zu identifizieren. Zum Beispiel kann ein zentrales Atom, das als Kation bezeichnet wird, von einem oder mehreren Anionen umgeben sein (Abbildung 4). Das Artbildungsskript identifiziert jedes einzelne der Koordinationspolyeder. Der gewichtete Durchschnitt aller Koordinationspolyeder ergibt die Koordinationsnummer, die mit der aus der Integration des PDF erhaltenen identisch ist. Geben Sie in der Befehlszeile Folgendes ein:speciation_umd.py -f -i -c -a -r 0HINWEIS: Durchschnittliche Koordinationszahlen in Fluiden sind Bruchzahlen. Diese Fraktionalität ergibt sich aus der durchschnittlichen Charakteristik der Koordination. Die auf der Artbildung basierende Definition liefert eine intuitivere und informativere Darstellung der Struktur der Flüssigkeit, wobei die relativen Anteile der verschiedenen Arten, d.h. Koordinationen, quantifiziert werden. Führen Sie das speciation_umd.py Skript mit dem Flag -r 1 aus, das den Konnektivitätsgraphen auf allen Tiefenebenen abtastet, um die Polymerisation zu erhalten. Das Netzwerk durch den Atomgraphen hat eine gewisse Tiefe, da Atome weiter entfernt an andere Bindungen gebunden sind (z. B. in Sequenzen von abwechselnden Kationen und Anionen) (Abbildung 4). Öffnen Sie die beiden Dateien . popul.dat und . stat.dat nacheinander; diese stellen die Ausgabe des Speziationsskripts dar. Jeder Cluster wird auf einer Zeile geschrieben und gibt seine chemische Formel, die Zeit, zu der er sich gebildet hat, die Zeit, zu der er starb, seine Lebensdauer, eine Matrix mit der Liste der Atome, die diesen Cluster bilden, an. Zeichnen Sie die Lebensdauer jedes Atomclusters aller in der Simulation gefundenen chemischen Spezies auf, wie in der .popul.dat-Datei (Abbildung 5) zu finden. Zeichnen Sie die Populationsanalyse mit der Häufigkeit jeder Art, wie in der . stat.dat Datei. Diese Analyse, sowohl absolut als auch relativ, entspricht der tatsächlichen Statistik der Koordinationspolyeder für den Fall -r 0; Im Falle der Polymerisation muss dies bei -r1 sorgfältig behandelt werden, da möglicherweise eine Gewisse Normalisierung über die relative Anzahl von Atomen angewendet werden muss. Die Häufigkeit entspricht dem Integral über die Lebenszeiten. Das. stat.dat Datei listet auch die Größe jedes Clusters auf, d.h. wie viele Atome ihn bilden. 4. Diffusionskoeffizienten berechnen Extrahieren Sie die mittleren quadratischen Verschiebungen (MSD) der Atome als Funktion der Zeit, um die Selbstdiffusivität zu erhalten. Die Standardformel des MSD lautet:wobei die Präfaktoren Renormalisierungen sind. Mit dem MSD-Tool gibt es verschiedene Möglichkeiten, die dynamischen Aspekte der Fluide zu analysieren.HINWEIS: T ist die Gesamtzeit der Simulation und N α ist die Anzahl der Atome vom Typ α. Die anfängliche Zeit t0 ist beliebig und erstreckt sich über die erste Hälfte der Simulation. Ninit ist die Anzahl der Anfangszeiten. τ ist die Breite des Zeitintervalls, über das die MSD berechnet werden; sein Maximalwert ist die Hälfte der Zeitdauer der Simulation. In typischen MSD-Implementierungen beginnt jedes Fenster am Ende des vorherigen. Eine spärlichere Stichprobe kann jedoch die Berechnung der MSD beschleunigen, ohne die Steigung der MSD zu verändern. Dazu beginnt das i-te Fenster zum Zeitpunkt t0(i), aber das (i+1)-te Fenster beginnt zum Zeitpunkt t0(i) + τ + v, wobei der Wert von v benutzerdefiniert ist. In ähnlicher Weise wird die Breite des Fensters in diskreten Schritten erhöht, die vom Benutzer als solche definiert werden: τ(i) = τ(i-1) + z. Die Werte von z (“horizontaler Schritt”) und v (“vertikaler Schritt”) sind positiv oder Null; Der Standardwert für beide ist 20. Berechnen Sie den MSD mithilfe der Reihe msd_umd Skripts. Ihre Ausgabe wird in einem gedruckt. msd.dat Datei, in der der MSD jedes atomaren Typs, Atoms oder Clusters als Funktion der Zeit auf einer Spalte gedruckt wird. Berechnen Sie die durchschnittliche MSD jedes atomaren Typs. Die MSD werden für jedes Atom berechnet und dann für jeden Atomtyp gemittelt. Die Ausgabedatei enthält eine Spalte für jeden atomaren Typ. Geben Sie in der Befehlszeile Folgendes ein:msd_umd.py -f -z -v -b Berechnen Sie den MSD jedes Atoms. Die MSD werden für jedes Atom berechnet und dann für jeden Atomtyp gemittelt. Die Ausgabedatei enthält eine Spalte für jedes Atom in der Simulation und dann eine Spalte für jeden Atomtyp. Diese Funktion ermöglicht die Identifizierung von Atomen, die in zwei verschiedenen Umgebungen wie Flüssigkeit und Gas oder zwei Flüssigkeiten diffundieren. Geben Sie in der Befehlszeile Folgendes ein:msd_all_umd.py -f -z -v -b Berechnen Sie die MSD der chemischen Spezies. Verwenden Sie die Population von Clustern, die mit dem Speziationsskript identifiziert und in der . popul.dat-Datei . Die MSD werden für jeden einzelnen Cluster berechnet. Die Ausgabedatei enthält eine Spalte für jeden Cluster. Um zu vermeiden, dass großflächige Polymere in Betracht gezogen werden, begrenzen Sie die Größe des Clusters. Der Standardwert ist 20 Atome. Geben Sie in der Befehlszeile Folgendes ein:msd_cluster_umd.py -f -p -s -b -c HINWEIS: Die Standardwerte sind: –b 100 –s 1 –c 20. Zeichnen Sie den MSD mithilfe einer tabellenbasierten Software (Abbildung 6). Identifizieren Sie in einer Log-Log-Darstellung des MSD im Vergleich zur Zeit die Steigungsänderung. Trennen Sie den ersten Teil, normalerweise kurz, der das ballistische Regime darstellt, d.h. die Erhaltung der Geschwindigkeit von Atomen nach Kollisionen. Der zweite längere Teil stellt das diffusive Regime dar, d.h. die Streuung der Geschwindigkeit von Atomen nach Kollisionen. Berechnen Sie die Diffusionskoeffizienten aus der Steigung des MSD wie folgt:wobei Z die Anzahl der Freiheitsgrade ist (Z = 2 für Diffusion in der Ebene, Z = 3 für die Diffusion im Raum) und t der Zeitschritt ist. 5. Zeitkorrelationsfunktionen Berechnen Sie die Zeitkorrelationsfunktionen als Maß für die Trägheit des Systems unter Verwendung der allgemeinen Formel:A kann eine Vielzahl von zeitabhängigen Variablen sein, wie die atomaren Positionen, Atomgeschwindigkeiten, Spannungen, Polarisation usw., die jeweils über die Green-Kubo-Beziehungen12,13 unterschiedliche physikalische Eigenschaften ergeben, manchmal nach einer weiteren Transformation. Analysieren Sie die Atomgeschwindigkeiten, um das Schwingungsspektrum der Flüssigkeit und den alternativen Ausdruck der atomaren Selbstdiffusionskoeffizienten zu erhalten. Führen Sie das skript vibr_spectrum_umd.py aus, um die VAC-Funktion (Atomic Velocity-Velocity Auto-Correlation) für jeden atomaren Typ zu berechnen und die schnelle Fourier-Transformation durchzuführen. Geben Sie in der Befehlszeile Folgendes ein:vibr_spectrum_umd.py -f -t wobei –t die Temperatur ist, die vom Benutzer definiert werden muss. Das Skript druckt zwei Dateien: die . vels.scf.dat Datei mit der VAC-Funktion für jeden atomaren Typ und . vibr.dat Datei mit dem auf jeder Atomspezies zerlegten Schwingungsspektrum und dem Gesamtwert. Öffnen und lesen Sie die vels.scf.dat. Zeichnen Sie die VAC-Funktion aus der Datei vels.scf.dat mit tabellenähnlicher Software. Behalten Sie den echten Teil des Fourier VAC. Dies ergibt das Schwingungsspektrum als Funktion der Frequenz:wobei m die Atommassen sind. Zeichnen Sie das Schwingungsspektrum aus der datei vibr.dat mit tabellenähnlicher Software (Abbildung 7). Identifizieren Sie den endlichen Wert bei ω=0, der dem diffusiven Charakter des Fluids und den verschiedenen Spitzen des Spektrums bei endlicher Frequenz entspricht. Identifizieren Sie die Beteiligung jedes atomaren Typs am Schwingungsspektrum.HINWEIS: Die Zerlegung auf Atomtypen zeigt, dass verschiedene Atome unterschiedliche ω=0-Beiträge haben, entsprechend ihren Diffusionskoeffizienten. Die allgemeine Form des Spektrums ist viel glatter mit weniger Merkmalen als bei einem entsprechenden Volumenkörper. Lesen Sie an der Schale das Integral über das Schwingungsspektrum, das die Diffusionskoeffizienten für jede Atomspezies ergibt.HINWEIS: Thermodynamische Eigenschaften können durch Integration aus dem Schwingungsspektrum erhalten werden, aber die Ergebnisse sollten wegen zweier Näherungen mit Vorsicht verwendet werden: Die Integration ist innerhalb der quasi-harmonischen Näherung gültig, die bei hohen Temperaturen nicht unbedingt gilt; und der gasähnliche Teil des Spektrums, der der Diffusion entspricht, muss verworfen werden. Die Integration sollte dann nur über den gitterartigen Teil des Spektrums erfolgen. Diese Trennung erfordert jedoch in der Regel mehrere weitere Nachbearbeitungsschritte und Berechnungen14, die vom vorliegenden UMD-Paket nicht abgedeckt werden. Führen Sie das skript viscosity_umd.py aus, um die Selbstkorrelation des Spannungstensors der Komponenten zu analysieren und die Viskosität der Schmelze zu schätzen. Geben Sie in der Befehlszeile Folgendes ein:viscosity_umd.py -f -i -s -o -l HINWEIS: Diese Funktion ist explorativ und alle Ergebnisse müssen mit Vorsicht betrachtet werden. Überprüfen Sie zunächst gründlich die Konvergenz der Viskosität in Bezug auf die Länge der Simulation. Leiten Sie die Viskosität der Flüssigkeit aus der Selbstkorrelation des Spannungstensors15 wie folgt ab:wobei V und T das Volumen bzw. die Temperatur sind, κB die Boltzmann-Konstante und σ ij die ij-off-diagonale Komponente des Spannungstensors, ausgedrückt in kartesischen Koordinaten. Verwenden Sie eine angemessenere Anpassung, um eine robustere Schätzung der Viskosität15,16 zu erhalten, und vermeiden Sie das Rauschen der Spannungstensor-Autokorrelationsfunktion, das sich aus der endlichen Größe und der endlichen Dauer der Simulationen ergeben könnte. Verwenden Sie für die Autokorrelationsfunktion des Spannungstensors die folgende Funktionsform15,16, die gute Ergebnisse liefert:wobei A, B, τ1, τ2 und ω Anpassungsparameter sind. Nach der Integration wird der Ausdruck für die Viskosität: 6. Thermodynamische Parameter, die sich aus den Simulationen ergeben. Führen Sie averages.py aus, um die Durchschnittswerte und die Spreizung (als Standardabweichung) für Druck, Temperatur, Dichte und interne Energie aus den umd-Dateien zu extrahieren. Geben Sie in der Befehlszeile Folgendes ein:averages.py -f -s mit –s 0 als Standard. Berechnen Sie den statistischen Fehler des Durchschnitts mit blockierenden Methoden.HINWEIS: Es gibt verschiedene Geschmacksrichtungen dieser Methode. In Anlehnung an die Arbeit von Allen und Tildesley2 ist es üblich, über Sequenzen von Zeitblöcken von immer längerer Länge zu mittelen und die Standardabweichung in Bezug auf das arithmetische Mittel zu schätzen17. Konvergenz kann im Grenzbereich von blockgrößen mit vielen und langen Blöcken erreicht werden, wenn die Stichprobe nicht korreliert ist. Der tatsächliche Schwellenwert für die Konvergenz muss jedoch in der Regel manuell ausgewählt werden. Verwenden Sie die Halbierungsmethode18: Beginnend mit der anfänglichen Datenstichprobe halbieren Sie bei jedem Schritt κ die Anzahl der Stichproben, indem Sie über jeweils zwei entsprechende aufeinanderfolgende Stichproben aus dem vorherigen Schritt κ−1 mittelen: Führen Sie das fullaverages.py-Skript aus, um die vollständige statistische Analyse einschließlich des Fehlers des Mittelwerts durchzuführen. Geben Sie in der Befehlszeile Folgendes ein:fullaverages.py -s -u HINWEIS: Das Skript wird so weit automatisiert, dass nach allen .umd.dat-Dateien im aktuellen Verzeichnis gesucht und die Analyse für alle durchgeführt wird. Die Standardwerte sind –s 0 –u 0. Für -u 0 ist die Ausgabe minimal, und für -u 1 ist die Ausgabe vollständig, wobei mehrere alternative Einheiten ausgedruckt werden. Dieses Skript erfordert grafische Unterstützung, da es ein grafisches Bild zur Überprüfung der Konvergenz zur Schätzung des Fehlers auf dem Mittelwert erstellt.

Representative Results

Pyrolit ist eine Modell-Mehrkomponenten-Silikatschmelze (0,5Na2O 2CaO 1,5Al2O3 4FeO 30MgO 24SiO2), die sich am besten der Zusammensetzung des Massensilikats Erde annähert – dem geochemischen Durchschnitt oder unserem gesamten Planeten, mit Ausnahme seines eisenbasierten Kerns19. Die frühe Erde wurde von einer Reihe groß angelegter Schmelzereignisse dominiert20, das letzte könnte den gesamten Planeten verschlungen haben, nachdem es für die Protolunarscheibe kondensiert hatte21. Pyrolit stellt den besten Näherungswert für die Zusammensetzung solcher planetarischen Magmaozeane dar. Folglich haben wir die physikalischen Eigenschaften der Pyrolitschmelze im Temperaturbereich von 3.000 bis 20125.000 K und im Druckbereich von 0 bis 0 bis 2012150 GPa aus den molekulardynamischen Simulationen von ab initio in der VASP-Implementierung umfassend untersucht. Diese thermodynamischen Bedingungen charakterisieren vollständig die extremsten Magma-Ozeanbedingungen der Erde. Unsere Studie ist ein hervorragendes Beispiel für den erfolgreichen Einsatz des UMD-Pakets für die gesamte tiefgehende Analyse der Schmelzen22. Wir berechneten die Verteilung und die durchschnittlichen Bindungslängen, verfolgten die Veränderungen in der Kationen-Sauerstoff-Koordination und verglichen unsere Ergebnisse mit früheren experimentellen und computergestützten Studien an amorphen Silikaten verschiedener Zusammensetzungen. Unsere eingehende Analyse half dabei, Standard-Koordinationszahlen in ihre Grundbestandteile zu zerlegen, das Vorhandensein exotischer Koordinationspolyeder in der Schmelze zu skizzieren und Lebensdauern für alle Koordinationspolyeder zu extrahieren. Es skizzierte auch die Bedeutung der Probenahme in Simulationen sowohl in Bezug auf die Länge der Flugbahn als auch auf die Anzahl der Atome, die in dem modellierten System vorhanden sind. Was die Nachbearbeitung betrifft, so ist die UMD-Analyse unabhängig von diesen Faktoren, sie sollten jedoch bei der Interpretation der vom UMD-Paket bereitgestellten Ergebnisse berücksichtigt werden. Hier zeigen wir einige Beispiele, wie das UMD-Paket verwendet werden kann, um mehrere charakteristische Merkmale der Schmelzen zu extrahieren, mit einer Anwendung auf geschmolzenen Pyrolit. Die si-O-Paarverteilungsfunktion, die aus dem gofrs_umd.py Skript erhalten wurde, zeigt, dass der Radius der ersten Koordinationskugel, die das erste Minimum der g(r)-Funktion ist, bei T = 3000 K und P = 4,6 GPa bei etwa 2,5 Angström liegt. Das Maximum des g(r) liegt bei 1,635 Å – dies ist die beste Annäherung an die Biegelänge. Der lange Schwanz ist auf die Temperatur zurückzuführen. Unter Verwendung dieser Grenze als Si-O-Bindungsabstand zeigt die Speziationsanalyse, dass SiO4-Einheiten, die bis zu einigen Pikosekunden dauern können, die Schmelze dominieren (Abbildung 5). Es gibt einen wichtigen Teil der Schmelze, der eine partielle Polymerisation zeigt, wie sie sich in der Anwesenheit von Dimeren wie Si2O7 und Trimern wie Si3Ox-Einheiten widerspiegelt. Ihre entsprechende Lebensdauer liegt in der Größenordnung der Pikosekunde. Polymere höherer Ordnung haben alle eine wesentlich kürzere Lebensdauer. Abbildung 5: Lebensdauer der chemischen Si-O-Spezies.Die Artbildung wird in einer Mehrkomponentenschmelze bei 4,6 GPa und 3000 K identifiziert. Die Etiketten kennzeichnen die Monomere SiO3, SiO4 und SiO5 sowie die verschiedenen SixOy-Polymere. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen. Die unterschiedlichen Werte der vertikalen und horizontalen Schritte, die durch die Flags –z und –v oben definiert sind, ergeben verschiedene Stichproben des MSD (Abbildung 6). Schon große Werte von z und v reichen aus, um die Steigungen und damit die Diffusionskoeffizienten der verschiedenen Atome zu definieren. Der Zeitgewinn für die Nachbearbeitung ist bemerkenswert, wenn man zu großen Werten von z und v geht. Der MSD bietet ein sehr starkes Validierungskriterium für die Qualität der Simulationen. Wenn der Diffusionsteil des MSD nicht lang genug ist, ist dies ein Zeichen dafür, dass die Simulation zu kurz ist und den fluiden Zustand im statistischen Sinne nicht erreicht. Die Mindestanforderung an den diffusiven Teil des MSD hängt stark vom System ab. Man kann verlangen, dass alle Atome ihren Standort mindestens einmal in der Struktur der Schmelze ändern, damit sie als Flüssigkeit betrachtet werden kann10. Ein hervorragendes Beispiel für Anwendungen in den Planetenwissenschaften sind komplexe Silikatschmelzen bei hohen Drücken nahe oder sogar unterhalb ihrer Liquiduslinie11. Die Si-Atome, die wichtigsten netzwerkbildenden Kationen, wechseln nach mehr als zwei Dutzend Pikosekunden ihre Standorte. Simulationen, die kürzer als dieser Schwellenwert sind, würden den möglichen Konfigurationsraum erheblich unterschätzen. Da sich die koordinierenden Anionen, nämlich die O-Atome, jedoch schneller bewegen als die zentralen Si-Atome, können sie einen Teil der langsamen Beweglichkeit von Si kompensieren. Als solches könnte das gesamte System in der Tat eine bessere Abtastung des Konfigurationsraums abdecken, als nur aus den Si-Verschiebungen angenommen. Abbildung 6: Mittlere quadratische Verschiebungen (MSD).Die MSD sind für einige atomare Typen einer Mehrkomponenten-Silikatschmelze dargestellt. Die Probenahme mit verschiedenen horizontalen und vertikalen Schritten, z und v, liefert konsistente Ergebnisse. Feste Kreise: -z 50 –v 50. Offene Kreise: -z 250 –v 500. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen. Schließlich liefern die atomaren VAC-Funktionen das Schwingungsspektrum der Schmelze. Abbildung 7 zeigt das Spektrum bei den gleichen Druck- und Temperaturbedingungen wie oben. Wir stellen die Beiträge der Mg-, Si- und O-Atome sowie den Gesamtwert dar. Bei der Nullfrequenz gibt es einen endlichen Wert des Spektrums, der dem Diffusionscharakter der Schmelze entspricht. Die Extraktion der thermodynamischen Eigenschaften aus dem Schwingungsspektrum muss diesen gasartigen diffusiven Charakter von Null entfernen, aber auch seinen Zerfall bei höheren Frequenzen richtig berücksichtigen. Abbildung 7: Das Schwingungsspektrum der Pyrolitschmelze.Der reelle Teil der Fourier-Transformation der atomaren Geschwindigkeits-Geschwindigkeits-Selbstkorrelationsfunktion ergibt das Schwingungsspektrum. Hier wird das Spektrum für eine Mehrkomponenten-Silikatschmelze berechnet. Flüssigkeiten haben bei Nullfrequenz einen gasartigen diffusiven Charakter ungleich Null. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Discussion

Das UMD-Paket wurde entwickelt, um besser mit Ab-initio-Simulationen zu arbeiten, bei denen die Anzahl der Snapshots typischerweise auf Zehn- bis Hunderttausende von Snapshots mit einigen hundert Atomen pro Einheitszelle begrenzt ist. Größere Simulationen sind ebenfalls handhabbar, sofern die Maschine, auf der die Nachbearbeitung läuft, über genügend aktive Speicherressourcen verfügt. Der Code zeichnet sich durch die Vielfalt der Eigenschaften aus, die er berechnen kann, und durch seine Open-Source-Lizenz.

Die umd.dat-Dateien eignen sich für die Ensembles, die die Anzahl der Partikel während der gesamten Simulation unverändert beibehalten. Das UMD-Paket kann Dateien lesen, die aus Berechnungen stammen, bei denen Form und Volumen der Simulationsbox variieren. Diese decken die gängigsten Berechnungen wie NVT und NPT ab, bei denen die Anzahl der Partikel, N, Temperatur T, Volumen, V und / oder Druck P, konstant gehalten werden.

Für den Zeitbeginn funktionieren die Paarverteilungsfunktion sowie alle Skripte, die die interatomaren Abstände schätzen müssen, wie die Speziationsskripte, nur für orthogonale Einheitszellen, dh für kubische, tetragonale und orthorhombische Zellen, bei denen die Winkel zwischen den Achsen 90 ° betragen.

Die Hauptentwicklungslinien für Version 2.0 sind die Aufhebung der Orthogonalitätsbeschränkung für Entfernungen und das Hinzufügen weiterer Funktionen für die Speziationsskripte: die Analyse einzelner chemischer Bindungen, die Analyse der interatomaren Winkel und die Implementierung der zweiten Koordinationssphäre. Mit Hilfe der externen Zusammenarbeit arbeiten wir daran, den Code auf eine GPU zu portieren, um eine schnellere Analyse in größeren Systemen zu ermöglichen.

Disclosures

The authors have nothing to disclose.

Acknowledgements

Diese Arbeit wurde vom Europäischen Forschungsrat (ERC) im Rahmen des Forschungs- und Innovationsprogramms Horizon 2020 der Europäischen Union (Förderkennzeichnungsnummer 681818 IMPACT to RC), vom Direktion extreme Physik und Chemie des Deep Carbon Observatory und vom norwegischen Forschungsrat über sein Förderprogramm für Exzellenzzentren, Projektnummer 223272, unterstützt. Wir bestätigen den Zugang zu den GENCI-Supercomputern durch die eDARI-Computing-Grants der stl2816-Serie, zum Irene AMD-Supercomputer durch das PRACE RA4947-Projekt und zum Fram-Supercomputer durch den UNINETT Sigma2 NN9697K. FS wurde durch ein Marie-Skłodowska-Curie-Projekt unterstützt (Finanzhilfevereinbarung ABISSE Nr. 750901).

Materials

getopt library open-source
glob library open-source
matplotlib library open-source
numpy library open-source
os library open-source
Python software The Python Software Foundation Version 2 and 3 open-source
random library open-source
re library open-source
scipy library open-source
subprocess library open-source
sys library open-source

References

  1. Frenkel, D., Smit, B. . Understanding Molecular Simulation. From Algorithms to Applications. , (2001).
  2. Allen, M. P., Tildesley, D. J., Allen, T. . Computer Simulation of Liquids. , (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. , (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).

Play Video

Cite This Article
Caracas, R., Kobsch, A., Solomatova, N. V., Li, Z., Soubiran, F., Hernandez, J. Analyzing Melts and Fluids from Ab Initio Molecular Dynamics Simulations with the UMD Package. J. Vis. Exp. (175), e61534, doi:10.3791/61534 (2021).

View Video