Analyzing Melts and Fluids from Ab Initio Molecular Dynamics Simulations with the UMD Package

Razvan Caracas; Anais Kobsch; Natalia  V. Solomatova; Zhi Li; Francois Soubiran; Jean-Alexis Hernandez

doi:10.3791/61534

JoVE Journal > Chemistry

Please note that all translations are automatically generated. Click here for the English version.

Chimica

Analyse van smelten en vloeistoffen van Ab Initio moleculaire dynamica simulaties met het UMD-pakket

Published: September 17, 2021

doi:

10.3791/61534

Razvan Caracas^1,2, Anais Kobsch¹, Natalia V. Solomatova¹, Zhi Li¹, Francois Soubiran^1,3, Jean-Alexis Hernandez^1,2

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

Summary

Smelten en vloeistoffen zijn alomtegenwoordige vectoren van massatransport in natuurlijke systemen. We hebben een open-source pakket ontwikkeld om ab initio moleculair-dynamische simulaties van dergelijke systemen te analyseren. We berekenen structurele (binding, clusterisatie, chemische speciatie), transport (diffusie, viscositeit) en thermodynamische eigenschappen (trillingsspectrum).

Abstract

We hebben een op Python gebaseerd open-source pakket ontwikkeld om de resultaten te analyseren die voortkomen uit ab initio moleculair-dynamische simulaties van vloeistoffen. Het pakket is het meest geschikt voor toepassingen op natuurlijke systemen, zoals silicaat- en oxidesmelt, vloeistoffen op waterbasis en verschillende superkritische vloeistoffen. Het pakket is een verzameling Python-scripts die twee belangrijke bibliotheken bevatten die zich bezighouden met bestandsindelingen en met kristallografie. Alle scripts worden uitgevoerd op de opdrachtregel. We stellen een vereenvoudigd formaat voor om de atomaire trajecten en relevante thermodynamische informatie van de simulaties op te slaan, die wordt opgeslagen in UMD-bestanden, wat staat voor Universal Molecular Dynamics. Het UMD-pakket maakt de berekening van een reeks structurele, transport- en thermodynamische eigenschappen mogelijk. Beginnend met de paarverdelingsfunctie definieert het bindingslengtes, bouwt het een interatomaire connectiviteitsmatrix en bepaalt het uiteindelijk de chemische soortvorming. Het bepalen van de levensduur van de chemische soort maakt het mogelijk om een volledige statistische analyse uit te voeren. Vervolgens berekenen speciale scripts de gemiddelde-kwadratische verplaatsingen voor de atomen en voor de chemische soorten. De geïmplementeerde zelfcorrelatieanalyse van de atoomsnelheden levert de diffusiecoëfficiënten en het trillingsspectrum op. Dezelfde analyse toegepast op de spanningen levert de viscositeit op. Het pakket is beschikbaar via de GitHub-website en via een eigen speciale pagina van het ERC IMPACT-project als open-access pakket.

Introduction

Vloeistoffen en smelten zijn actieve chemische en fysische transportvectoren in natuurlijke omgevingen. De verhoogde snelheden van atomaire diffusie bevorderen chemische uitwisselingen en reacties, de lage viscositeit in combinatie met variërend drijfvermogen bevordert grote massaoverdracht en kristalsmeltdichtheidsrelaties bevorderen gelaagdheid in planetaire lichamen. De afwezigheid van een periodiek rooster, typische hoge temperaturen die nodig zijn om de gesmolten toestand te bereiken en de moeilijkheid om te doven, maken de experimentele bepaling van een reeks voor de hand liggende eigenschappen, zoals dichtheid, diffusie en viscositeit, uiterst uitdagend. Deze moeilijkheden maken alternatieve computationele methoden sterk en nuttige hulpmiddelen voor het onderzoeken van deze klasse van materialen.

Met de komst van rekenkracht en de beschikbaarheid van supercomputers worden momenteel twee belangrijke numerieke atomistische simulatietechnieken gebruikt om de dynamische toestand van een niet-kristallijn atomistisch systeem, Monte ^Carlo1 en moleculaire dynamica (MD)^1,2, te ^bestuderen. In Monte Carlo-simulaties wordt de configuratieruimte willekeurig bemonsterd; Monte Carlo-methoden tonen lineaire schaalvergroting in parallellisatie als alle bemonsteringswaarnemingen onafhankelijk van elkaar zijn. De kwaliteit van de resultaten hangt af van de kwaliteit van de random number generator en de representativiteit van de bemonstering. Monte Carlo-methoden tonen lineaire schaalvergroting in parallellisatie als de bemonstering onafhankelijk van elkaar is. In de moleculaire dynamica (MD) wordt de configuratieruimte bemonsterd door tijdsafhankelijke atomaire trajecten. Uitgaande van een bepaalde configuratie worden de atomaire trajecten berekend door de Newtoniaanse bewegingsvergelijkingen te integreren. De interatomaire krachten kunnen worden berekend met behulp van modelinteratomaire potentialen (in klassieke MD) of met behulp van first-principles-methoden (in ab initio, of first-principles, MD). De kwaliteit van de resultaten hangt af van de lengte van het traject en het vermogen om niet aangetrokken te worden door lokale minima.

Moleculaire dynamica simulaties bevatten een overvloed aan informatie, allemaal gerelateerd aan het dynamische gedrag van het systeem. Thermodynamische gemiddelde eigenschappen, zoals interne energie, temperatuur en druk, zijn vrij standaard om te berekenen. Ze kunnen worden geëxtraheerd uit de uitvoerbestanden van de simulaties en worden gemiddeld, terwijl grootheden die rechtstreeks verband houden met de beweging van de atomen en met hun onderlinge relatie moeten worden berekend na extractie van de atomaire posities en snelheden.

Bijgevolg is er veel aandacht besteed aan het visualiseren van de resultaten en zijn er vandaag verschillende pakketten beschikbaar, op verschillende platforms, open source of niet [^Ovito3, ^VMD4, ^Vesta5, ^Travis6, enz.]. Al deze visualisatietools gaan efficiënt om met interatomaire afstanden en maken als zodanig de efficiënte berekening van paarverdelingsfuncties en diffusiecoëfficiënten mogelijk. Verschillende groepen die grootschalige moleculaire dynamica-simulaties uitvoeren, hebben eigen software om verschillende andere eigenschappen te analyseren die voortvloeien uit de simulaties, soms in shareware of andere vormen van beperkte toegang tot de gemeenschap, en soms beperkt in omvang en gebruik tot sommige specifieke pakketten. Geavanceerde algoritmen om informatie te extraheren over interatomaire binding, geometrische patronen en thermodynamica worden ontwikkeld en geïmplementeerd in sommige van deze pakketten3,4,5,6,7, enz.

Hier stellen we het UMD-pakket voor – een open-source pakket geschreven in Python om de output van moleculaire dynamicasimulaties te analyseren. Het UMD-pakket maakt de berekening van een breed scala aan structurele, dynamische en thermodynamische eigenschappen mogelijk (figuur 1). Het pakket is beschikbaar via de GitHub-website (https://github.com/rcaracas/UMD_package) en via een speciale pagina (http://moonimpact.eu/umd-package/) van het ERC IMPACT-project als een open-access pakket.

Om het universeel en gemakkelijker te hanteren te maken, is onze aanpak om eerst alle informatie met betrekking tot de thermodynamische toestand en de atomaire trajecten uit het uitvoerbestand van de werkelijke moleculair-dynamische run te extraheren. Deze informatie wordt opgeslagen in een speciaal bestand, waarvan de indeling onafhankelijk is van het oorspronkelijke MD-pakket waarin de simulatie werd uitgevoerd. We noemen deze bestanden “umd” -bestanden, wat staat voor Universal Molecular Dynamics. Op deze manier kan ons UMD-pakket gemakkelijk worden gebruikt door elke ab initio-groep met elke software, allemaal met een minimale inspanning van aanpassing. De enige vereiste om het huidige pakket te gebruiken, is om de juiste parser van de uitvoer van de specifieke MD-software naar het umd-bestandsformaat te schrijven, als dit nog niet bestaat. Voorlopig leveren wij dergelijke parsers voor de ^VASP8 en de ^QBox9 pakketten.

Figuur 1: Stroomdiagram van de UMD-bibliotheek.
Fysieke eigenschappen zijn in blauw en de belangrijkste Python-scripts en hun opties zijn in het rood. Klik hier om een grotere versie van deze figuur te bekijken.

De umd-bestanden zijn ASCII-bestanden; typische uitbreiding is “umd.dat” maar niet verplicht. Alle analysecomponenten kunnen ASCII-bestanden van het umd-formaat lezen, ongeacht de werkelijke naamextensie. Sommige van de automatische scripts die zijn ontworpen om snelle grootschalige statistieken uit te voeren over verschillende simulaties, zoeken echter specifiek naar bestanden met de extensie umd.dat. Elke fysieke eigenschap wordt uitgedrukt op één regel. Elke regel begint met een trefwoord. Op deze manier is het formaat zeer aanpasbaar en maakt het mogelijk om nieuwe eigenschappen toe te voegen aan het umd-bestand, terwijl de leesbaarheid in alle versies behouden blijft. De eerste 30 regels van het umd-bestand van de simulatie van pyroliet bij 4,6 GPa en 3000 K, die hieronder in de discussie worden gebruikt, zijn weergegeven in figuur 2.

Figuur 2: Het begin van het umd-bestand waarin de simulatie van vloeibaar pyroliet bij 4,6 GPa en 3000 K wordt beschreven.
De koptekst wordt gevolgd door de beschrijving van elke momentopname. Elke eigenschap wordt op één regel geschreven, met de naam van de fysieke eigenschap, de waarde(n) en de eenheden, allemaal gescheiden door spaties. Klik hier om een grotere versie van deze figuur te bekijken.

Alle umd-bestanden bevatten een header die de inhoud van de simulatiecel beschrijft: het aantal atomen, elektronen en atoomtypen, evenals details voor elk atoom, zoals het type, het chemische symbool, het aantal valentie-elektronen en de massa. Een lege regel markeert het einde van de koptekst en scheidt deze van het hoofdgedeelte van het umd-bestand.

Vervolgens wordt elke stap van de simulatie gedetailleerd beschreven. Eerst worden de momentane thermodynamische parameters gegeven, elk op een andere lijn, waarbij (i) de naam van de parameter wordt gespecificeerd, zoals energie, spanningen, equivalente hydrostatische druk, dichtheid, volume, roosterparameters, enz., (ii) de waarde (en) en (iii) de eenheden. Een tabel die de atomen beschrijft, komt daarna. Een koplijn geeft de verschillende maten weer, zoals cartesiaanse posities, snelheden, ladingen, enz., En hun eenheden. Vervolgens wordt elk atoom op één lijn gedetailleerd. Door groepen van drie, overeenkomend met de drie x-, y-, z-assen , zijn de ingangen: de gereduceerde posities, de cartesische posities gevouwen in de simulatiecel, de Cartesiaanse posities (die goed rekening houden met het feit dat atomen tijdens een simulatie verschillende eenheidscellen kunnen doorkruisen), de atoomsnelheden en de atoomkrachten. De laatste twee ingangen zijn scalairen: lading en magnetisch moment.

Twee grote bibliotheken zorgen voor de goede werking van het hele pakket. De umd_process.py bibliotheek behandelt de umd-bestanden, zoals lezen en afdrukken. De crystallography.py bibliotheek behandelt alle informatie met betrekking tot de eigenlijke atomaire structuur. De onderliggende filosofie van de crystallography.py bibliotheek is om het rooster te behandelen als een vectoriële ruimte. De eenheidscelparameters vertegenwoordigen samen met hun oriëntatie de basisvectoren. De “Ruimte” heeft een reeks scalaire attributen (specifiek volume, dichtheid, temperatuur en specifiek aantal atomen), thermodynamische eigenschappen (interne energie, druk, warmtecapaciteit, enz.) en een reeks tensoriale eigenschappen (spanning en elasticiteit). Atomen bevolken deze ruimte. De klasse “Lattice” definieert dit ensemble, naast verschillende korte berekeningen, zoals specifiek volume, dichtheid, het verkrijgen van het reciproke rooster uit het directe, enz. De klasse “Atomen” definieert de atomen. Ze worden gekenmerkt door een reeks scalaire eigenschappen (naam, symbool, massa, aantal elektronen, enz.) en een reeks vectoriële eigenschappen (de positie in de ruimte, hetzij ten opzichte van de vectoriële basis beschreven in de klasse Lattice, of ten opzichte van universele Cartesische coördinaten, snelheden, krachten, enz.). Afgezien van deze twee klassen bevat de crystallography.py bibliotheek een reeks functies om een verscheidenheid aan tests en berekeningen uit te voeren, zoals atomaire afstanden of celvermenigvuldiging. Het periodiek systeem van de elementen is ook opgenomen als woordenboek.

De verschillende componenten van het umd-pakket schrijven verschillende uitvoerbestanden. Als algemene regel zijn het allemaal ASCII-bestanden, al hun vermeldingen worden gescheiden door tabbladen en ze zijn zo vanzelfsprekend mogelijk gemaakt. Ze geven bijvoorbeeld altijd duidelijk de fysieke eigenschap en zijn eenheden aan. De umd.dat bestanden voldoen volledig aan deze regel.

Protocol

1. Analyse van de moleculair-dynamische runs OPMERKING: Het pakket is beschikbaar via de GitHub-website (https://github.com/rcaracas/UMD_package) en via een speciale pagina (http://moonimpact.eu/umd-package/) van het ERC IMPACT-project als een open-access pakket. Extraheer elke specifieke set fysieke eigenschappen met behulp van een of meer speciale Python-scripts uit het pakket. Voer alle scripts uit op de opdrachtregel; ze gebruiken allemaal een reeks vlaggen, die zo consistent mogelijk zijn van het ene script naar het andere. De vlaggen, hun betekenis en standaardwaarden zijn allemaal samengevat in tabel 1. Vlag Betekenis Script dat het gebruikt Standaardwaarde -h Korte hulp alle -f UMD-bestandsnaam alle -i Thermalisatiestappen die moeten worden weggegooid alle 0 -i Invoerbestand met de interatomaire bindingen Soortvorming bonds.invoer -s Bemonstering van de frequentie msd, soortvorming 1 (elke stap wordt overwogen) -een Lijst van atomen of anionen Soortvorming -c Lijst van kationen Soortvorming -l Bond lengte Soortvorming 2 -t Temperatuur trillingen, reologie -v Discretisatie van de breedte van het bemonsteringsvenster van het traject voor de gemiddelde-kwadratische verplaatsingsanalyse Msd 20 -z Discretisatie van het begin van het bemonsteringsvenster van het traject voor de gemiddelde-kwadratische verplaatsingsanalyse Msd 20 Tabel 1: Meest voorkomende vlaggen die in het UMD-pakket worden gebruikt en hun meest voorkomende betekenis. Begin met het transformeren van de uitvoer van de MD-simulatie uitgevoerd in een first-principles code, zoals VASP8 of QBox9, in een UMD-bestand. Als de MD-simulaties in VASP zijn uitgevoerd, typt u op de opdrachtregel:VaspParser.py -f -i waarbij de vlag –f de naam van het VASP OUTCAR-bestand definieert en de –i de thermalisatielengte.OPMERKING: De eerste stap, gedefinieerd door -i, maakt het mogelijk om de eerste stappen van de simulaties, die de thermalisatie vertegenwoordigen, weg te gooien. In een typische moleculair-dynamische run vertegenwoordigt het eerste deel van de berekening de thermalisatie, d.w.z. de tijd die het systeem nodig heeft om alle atomen een Gauss-achtige verdeling van de temperatuur te beschrijven, en voor het hele systeem om fluctuaties van de temperatuur, druk, energie, enz. rond evenwichtswaarden te vertonen. Dit thermalisatiegedeelte van de simulatie mag niet in aanmerking worden genomen bij het analyseren van de statistische eigenschappen van de vloeistof. Transformeer de . umd bestanden in . xyz-bestanden om visualisatie op verschillende andere pakketten, zoals VMD4 of Vesta5, te vergemakkelijken. Typ op de opdrachtregel:umd2xyz.py -f -i -s waarbij –f de naam van de . umd-bestand , -i definieert de thermalisatieperiode die moet worden weggegooid, en -s de frequentie van de bemonstering van het traject dat is opgeslagen in de . umd-bestand . Standaardwaarden zijn –i 0 –s 1, d.w.z. rekening houdend met alle stappen van de simulatie, zonder dat deze worden weggegooid. Draai het umd-bestand om in VASP-type POSCAR-bestanden met behulp van het umd2poscar.py script; snapshots van de simulaties kunnen worden geselecteerd met een vooraf gedefinieerde frequentie. Typ op de opdrachtregel:umd2poscar.py -f -i -l -s waarbij –l de laatste stap vertegenwoordigt die moet worden omgezet in POSCAR-bestand. Standaardwaarden zijn -i 0 -l 10000000 -s 1. Deze waarde van -l is groot genoeg om een typisch heel traject te bestrijken. 2. Voer de structurele analyse uit Voer het gofrs_umd.py script uit om de paarverdelingsfunctie (PDF) gᴀʙ(r) te berekenen voor alle paren van atoomtypen A en B (afbeelding 3). De uitvoer is geschreven in één ASCII-bestand, gescheiden door tabs, met de extensie gofrs.dat. Typ op de opdrachtregel:gofrs_umd.py -f -s -d -i OPMERKING: De standaardwaarden zijn Sampling_Frequency (de frequentie voor het bemonsteren van het traject) = 1 stap; DiscretizationInterval (voor het plotten van de g(r)) = 0,01 Å; InitialStep (aantal stappen in het begin van het traject dat wordt weggegooid) = 0. De radiale PDF, gᴀʙ(r) is het gemiddelde aantal atomen van type B op een afstand d_ᴀʙ binnen een bolvormige schil met straal r en dikte gecentreerd op de atomen van type A (figuur 3):met ρ de atoomdichtheid, NA en NB het aantal atomen van type A en B, en δ(r−rᴀʙ) de deltafunctie die gelijk is aan 1 als de atomen A en B op een afstand liggen tussen r en r+dr. De abscis van het eerste maximum van gᴀʙ(r) geeft de hoogste waarschijnlijkheid bindingslengte tussen de atomen van type A en B, die het dichtst bij een gemiddelde bindingsafstand ligt die we kunnen bepalen. Het eerste minimum bakent de omvang van de eerste coördinatiesfeer af. Vandaar dat de integraal over de PDF tot het eerste minimum het gemiddelde coördinatienummer geeft. De som van de Fouriertransformaties van de gᴀʙ(r) voor alle paren van atoomtypen A en B levert het diffractiepatroon van de vloeistof op, zoals experimenteel verkregen met een diffractometer. In werkelijkheid echter, omdat vaak de hoge orde coördinatiesferen ontbreken in de gᴀʙ(r), kan het diffractiepatroon niet in zijn geheel worden verkregen. Figuur 3: Bepaling van de paarverdelingsfunctie.(a) Voor elk atoom van één soort (bijvoorbeeld rood) worden alle atomen van de coördinerende soort (bijvoorbeeld grijs en/of rood) geteld als functie van de afstand. (b) De resulterende afstandsverdelingsgrafiek voor elke momentopname, die in dit stadium slechts een verzameling deltafuncties is, wordt vervolgens gemiddeld over alle atomen en alle momentopnamen en gewogen met de ideale gasverdeling om (c) de paarverdelingsfunctie te genereren die continu is. Het eerste minimum van de g(r) is de straal van de eerste coördinatiebol, die later in de soortvormingsanalyse wordt gebruikt. Klik hier om een grotere versie van deze figuur te bekijken. Haal de gemiddelde interatomaire bindingsafstanden eruit als de radii van de eerste coördinatiesferen. Identificeer hiervoor de positie van het eerste maximum van de gᴀʙ(r)-functies: plot het gofrs.dat bestand in een spreadsheettoepassing en zoek naar de maxima en minima voor elk paar atomen. Identificeer de straal van de eerste coördinatiesfeer, als het eerste minimum van de PDF, gᴀʙ(r), met behulp van spreadsheetsoftware. Dit is de basis voor de gehele structurele analyse van de vloeistof; de PDF geeft de gemiddelde bindingsstatus van de atomen in de vloeistof. Extraheer de afstanden van de eerste minima, d.w.z. de abscis, en schrijf ze in een apart bestand, bijvoorbeeld bonds.input genoemd. U kunt ook een van de analyze_gofr scripts van het UMD-pakket uitvoeren om de maxima en de minima van de gᴀʙ(r)-functies te identificeren. Typ op de opdrachtregel:analyze_gofr_semi_automatic.py Klik op de positie van het maximum en het minimum van de gᴀʙ(r)-functie die wordt weergegeven in de grafiek die door het programma wordt geopend. Het script scant automatisch de huidige map, identificeert alle gofrs.dat bestanden en voert de analyse voor elk van hen uit. Klik nogmaals op het maximum en het minimum in het venster telkens wanneer het script een weloverwogen eerste gok nodig heeft. Open en bekijk het automatisch gegenereerde bestand met de naam bonds.input dat de interatomaire bindingsafstanden bevat. 3. Voer de soortvormingsanalyse uit Bereken de topologie van binding tussen de atomen, met behulp van het concept van connectiviteit binnen de grafentheorie: de atomen zijn de knooppunten en de interatomaire bindingen zijn de paden. Het speciation_umd.py script heeft de interatomaire bindingsafstanden nodig die zijn gedefinieerd in het bestand bonds.input .OPMERKING: De connectiviteitsmatrix wordt bij elke tijdstap geconstrueerd: twee atomen die op een afstand liggen die kleiner is dan de straal van hun overeenkomstige eerste coördinatiebol, worden beschouwd als gebonden, d.w.z. verbonden. Verschillende atomaire netwerken worden gebouwd door de atomen te behandelen als knooppunten in een grafiek waarvan de verbindingen worden gedefinieerd door dit geometrische criterium. Deze netwerken zijn de atomaire soort en hun ensemble definieert de atomaire soortvorming in die specifieke vloeistof (figuur 4). Figuur 4: Identificatie van de atoomclusters.De coördinatiepolyhedra worden gedefinieerd met behulp van interatomaire afstanden. Alle atomen op een afstand kleiner dan een bepaalde straal worden geacht gebonden te zijn. Hier komt de drempel overeen met de eerste coördinatiebol (de lichtrode cirkels), gedefinieerd in figuur 1. Polymerisatie en dus chemische soorten worden verkregen uit de netwerken van de gebonden atomen. Let op de centrale Red1Grey2-cluster, die geïsoleerd is van de andere atomen, die een oneindig polymeer vormen. Klik hier om een grotere versie van deze figuur te bekijken. Voer het soortvormingsscript uit om de connectiviteitsmatrix te verkrijgen en de coördinatiepolyhedra of de polymerisatie te verkrijgen. Typ op de opdrachtregel:speciation_umd.py -f -s -i -l -c -a -m -r waarbij de vlag -i het bestand geeft met de interatomaire bindingsafstanden, die bijvoorbeeld in de vorige stap zijn geproduceerd. U kunt het script ook uitvoeren met één enkele lengte voor alle bindingen die zijn gedefinieerd door de vlag -l.OPMERKING: De vlag -c geeft de centrale atomen aan en de vlag -a de liganden. Zowel centrale atomen als liganden kunnen van verschillende typen zijn; in dit geval moeten ze worden gescheiden door komma’s. De -m-vlag geeft de minimale tijd aan die een soort moet leven om in de analyse te worden meegenomen. Standaard is deze minimale tijd nul, waarbij alle voorvallen in de uiteindelijke analyse worden meegeteld. Voer het speciation_umd.py script uit met de vlag -r 0, die de connectiviteitsgrafiek op het eerste niveau bemonstert om de coördinatiepolyhedra te identificeren. Een centraal atoom, aangeduid als een kation , kan bijvoorbeeld worden omgeven door een of meer anionen (figuur 4). Het soortvormingsscript identificeert elk van de coördinatiepolydra. Het gewogen gemiddelde van alle coördinatiepolydra geeft het coördinatienummer, identiek aan het nummer dat is verkregen uit de integratie van de PDF. Typ op de opdrachtregel:speciation_umd.py -f -i -c -a -r 0OPMERKING: Gemiddelde coördinatienummers in vloeistoffen zijn fractionele getallen. Deze fractionaliteit komt voort uit het gemiddelde kenmerk van de coördinatie. De definitie op basis van soortvorming levert een meer intuïtieve en informatieve weergave van de structuur van de vloeistof op, waarbij de relatieve verhoudingen van de verschillende soorten, d.w.z. coördinaties, worden gekwantificeerd. Voer het speciation_umd.py script uit met de vlag -r 1, die de connectiviteitsgrafiek op alle diepteniveaus bemonstert om de polymerisatie te verkrijgen. Het netwerk door de atoomgrafiek heeft een bepaalde diepte, omdat atomen verder weg gebonden zijn aan andere bindingen (bijvoorbeeld in reeksen van afwisselend kationen en anionen) (figuur 4). Open de twee bestanden . popul.dat en . stat.dat achtereenvolgens; deze vormen de output van het soortschrift. Elke cluster is op één regel geschreven, met vermelding van de chemische formule, het tijdstip waarop het is gevormd, het tijdstip waarop het stierf, zijn levensduur, een matrix met de lijst van de atomen die deze cluster vormen. Plot de levensduur van elke atoomcluster van alle chemische soorten die in de simulatie worden aangetroffen, zoals gevonden in het bestand .popul.dat (figuur 5). Zet de populatieanalyse uit met de abundantie van elke soort, zoals te vinden in de . stat.dat bestand. Deze analyse, zowel absoluut als relatief, komt overeen met de feitelijke statistieken van de coördinatiepolyt voor de zaak -r 0; voor het geval van polymerisatie, met -r 1, moet dit zorgvuldig worden behandeld, omdat enige normalisatie over het relatieve aantal atomen mogelijk moet worden toegepast. De overvloed komt overeen met de integraal over de levens. De. stat.dat bestand vermeldt ook de grootte van elk cluster, d.w.z. hoeveel atomen het vormen. 4. Bereken diffusiecoëfficiënten Extraheer de gemiddelde vierkante verplaatsingen (MSD) van de atomen als functie van de tijd om de zelfdiffusie te verkrijgen. De standaardformule van de MSD is:waarbij de prefactoren renormalisaties zijn. Met de MSD-tool zijn er verschillende manieren om de dynamische aspecten van de vloeistoffen te analyseren.OPMERKING: T is de totale tijd van de simulatie en N α is het aantal atomen van het type α. De begintijd t0 is arbitrair en beslaat de eerste helft van de simulatie. Ninit is het aantal begintijden. τ = de breedte van het tijdsinterval waarover de msd worden berekend; de maximale waarde is de helft van de tijdsduur van de simulatie. Bij typische MSD-implementaties begint elk venster aan het einde van het vorige venster. Maar een sparser-bemonstering kan de berekening van de MSD versnellen, zonder de helling van de MSD te veranderen. Hiervoor begint het i-de venster op tijdstip t0(i), maar het (i+1)-de venster begint op tijdstip t0(i) + τ + v, waarbij de waarde van v door de gebruiker is gedefinieerd. Evenzo wordt de breedte van het venster vergroot in discrete stappen die door de gebruiker worden gedefinieerd, als zodanig: τ(i) = τ(i-1) + z. De waarden van z (“horizontale stap”) en v (“verticale stap”) zijn positief of nul; de standaardwaarde voor beide is 20. Bereken de MSD met behulp van de reeks msd_umd scripts. Hun uitvoer wordt afgedrukt in een . msd.dat bestand, waarbij de MSD van elk atoomtype, atoom of cluster in één kolom wordt afgedrukt als functie van de tijd. Bereken de gemiddelde MSD van elk atoomtype. De MSD worden berekend voor elk atoom en vervolgens gemiddeld voor elk atoomtype. Het uitvoerbestand bevat één kolom voor elk atoomtype. Typ op de opdrachtregel:msd_umd.py -f -z -v -b Bereken de MSD van elk atoom. De MSD worden berekend voor elk atoom en vervolgens gemiddeld voor elk atoomtype. Het uitvoerbestand bevat één kolom voor elk atoom in de simulatie en vervolgens één kolom voor elk atoomtype. Deze functie maakt het mogelijk om atomen te identificeren die diffunderen in twee verschillende omgevingen, zoals vloeistof en gas, of twee vloeistoffen. Typ op de opdrachtregel:msd_all_umd.py -f -z -v -b Bereken de MSD van de chemische soort. Gebruik de populatie van clusters die zijn geïdentificeerd met het soortschrift en afgedrukt in de . popul.dat bestand. De MSD wordt berekend voor elk afzonderlijk cluster. Het uitvoerbestand bevat één kolom voor elk cluster. Om te voorkomen dat grootschalige polymeren worden overwogen, stelt u een limiet aan de grootte van het cluster; de standaardwaarde is 20 atomen. Typ op de opdrachtregel:msd_cluster_umd.py -f -p -s -b -c OPMERKING: Standaardwaarden zijn: –b 100 –s 1 –c 20. Plot de MSD met behulp van een spreadsheet-gebaseerde software (Figuur 6). Identificeer in een log-logweergave van de MSD versus de tijd de hellingsverandering. Scheid het eerste deel, meestal kort, dat het ballistische regime vertegenwoordigt, d.w.z. het behoud van de snelheid van atomen na botsingen. Het tweede langere deel vertegenwoordigt het diffusieregime , d.w.z. verstrooiing van de snelheid van atomen na botsingen. Bereken de diffusiecoëfficiënten van de helling van het MSD als volgt:waarbij Z het aantal vrijheidsgraden is (Z = 2 voor diffusie in vlak, Z = 3 voor diffusie in de ruimte), en t de tijdstap is. 5. Tijd correlatie functies Bereken de tijdcorrelatiefuncties als een maat voor de traagheid van het systeem met behulp van de algemene formule:A kan een verscheidenheid aan tijdsafhankelijke variabelen zijn, zoals de atomaire posities, atoomsnelheden, spanningen, polarisatie, enz., Die elk – via de Green-Kubo-relaties12,13 – verschillende fysische eigenschappen opleveren, soms na een verdere transformatie. Analyseer de atoomsnelheden om het trillingsspectrum van de vloeistof en alternatieve expressie van de atomaire zelfdiffusiecoëfficiënten te verkrijgen. Voer het vibr_spectrum_umd.py script uit om de atomaire snelheid-snelheid auto-correlatie (VAC) functie voor elk atoomtype te berekenen en voer de snelle Fouriertransformatie uit. Typ op de opdrachtregel:vibr_spectrum_umd.py -f -t waarbij –t de temperatuur is die door de gebruiker moet worden bepaald. Het script drukt twee bestanden af: de . vels.scf.dat bestand met de functie VAC voor elk atoomtype en de . vibr.dat bestand met het trillingsspectrum dat op elke atoomsoort is afgebroken en de totale waarde. Open en lees de vels.scf.dat. Plot de VAC-functie uit het vels.scf.dat bestand met behulp van spreadsheet-achtige software. Houd het echte deel van de Fourier VAC. Dit is wat het trillingsspectrum oplevert, als functie van de frequentie:waarbij m de atoommassa’s zijn. Plot het trillingsspectrum van het vibr.dat bestand met behulp van spreadsheetachtige software (figuur 7). Identificeer de eindige waarde bij ω=0 die overeenkomt met het diffusieve karakter van de vloeistof en de verschillende pieken van het spectrum bij eindige frequentie. Identificeer de deelname van elk atoomtype aan het trillingsspectrum.OPMERKING: De ontleding op atoomtypen laat zien dat verschillende atomen verschillende ω=0-bijdragen hebben, overeenkomend met hun diffusiecoëfficiënten. De algemene vorm van het spectrum is veel vloeiender met minder functies dan voor een overeenkomstige vaste stof. Lees bij de schaal de integraal over het trillingsspectrum af, wat de diffusiecoëfficiënten voor elke atoomsoort oplevert.OPMERKING: Thermodynamische eigenschappen kunnen worden verkregen door integratie uit het trillingsspectrum, maar de resultaten moeten met voorzichtigheid worden gebruikt vanwege twee benaderingen: de integratie is geldig binnen de quasi-harmonische benadering, die niet noodzakelijkerwijs geldt bij hoge temperaturen; en het gasachtige deel van het spectrum dat overeenkomt met de diffusie moet worden weggegooid. De integratie moet dan alleen over het roosterachtige deel van het spectrum gebeuren. Maar deze scheiding vereist meestal verschillende verdere nabewerkingsstappen en berekeningen14, die niet onder het huidige UMD-pakket vallen. Voer het viscosity_umd.py-script uit om de zelfcorrelatie van de componentenspanningstensor te analyseren om de viscositeit van de smelt te schatten. Typ op de opdrachtregel:viscosity_umd.py -f -i -s -o -l OPMERKING: Deze functie is verkennend en eventuele resultaten moeten met de nodige voorzichtigheid worden genomen. Controleer eerst grondig de convergentie van de viscositeit met betrekking tot de lengte van de simulatie. Leid de viscositeit van de vloeistof af uit de zelfcorrelatie van de spanningstensor15 als:waarbij V en T respectievelijk het volume en de temperatuur zijn, is κB de constante van Boltzmann en σ ij de ij off-diagonale component van de spannings-tensor, uitgedrukt in Cartesische coördinaten. Gebruik een adequatere pasvorm om een robuustere schatting van de viscositeit te verkrijgen15,16 en vermijd de ruis van de spanning-tensor autocorrelatiefunctie die zou kunnen voortvloeien uit de eindige grootte en de eindige duur van de simulaties. Gebruik voor de autocorrelatiefunctie van de spanningstensor de volgende functionele vorm15,16 die goede resultaten oplevert:waarbij A, B, τ1, τ2 en ω geschikte parameters zijn. Na integratie wordt de uitdrukking voor de viscositeit: 6. Thermodynamische parameters die voortvloeien uit de simulaties. Voer averages.py uit om de gemiddelde waarden en de spread (als standaarddeviatie) voor druk, temperatuur, dichtheid en interne energie uit de umd-bestanden te extraheren. Typ op de opdrachtregel:averages.py -f -s met –s 0 als standaard. Bereken de statistische fout van het gemiddelde met behulp van blokkeringsmethoden.OPMERKING: Er zijn verschillende smaken van deze methode. Volgens het werk van Allen en Tildesley2 is het gebruikelijk om over reeksen tijdsblokken, van steeds langere lengte, te gemiddelden en de standaarddeviatie ten opzichte van het rekenkundig gemiddelde te schatten17. Convergentie kan worden bereikt in de limiet van blokgroottes die veel en lang genoeg zijn, wanneer de steekproef niet-gecorreleerd is. Hoewel de werkelijke drempelwaarde voor de convergentie meestal handmatig moet worden gekozen. Gebruik de halveringsmethode18: te beginnen met het eerste gegevensmonster, halveer bij elke stap κ het aantal monsters door het gemiddelde te nemen over elke twee overeenkomstige opeenvolgende monsters van de vorige stap κ−1: Voer het fullaverages.py script uit om de volledige statistische analyse uit te voeren, inclusief de fout van het gemiddelde. Typ op de opdrachtregel:fullaverages.py -s -u OPMERKING: Het script is geautomatiseerd tot het punt van het zoeken naar alle .umd.dat bestanden in de huidige map en het uitvoeren van de analyse voor alle bestanden. Standaardinstellingen zijn –s 0 –u 0. Voor -u 0 is de uitvoer minimaal en voor -u 1 is de uitvoer volledig, met verschillende alternatieve eenheden afgedrukt. Dit script vereist grafische ondersteuning, omdat het een grafische afbeelding creëert voor het controleren van de convergentie voor het schatten van de fout op het gemiddelde.

Representative Results

Pyroliet is een model multi-component silicaat smelt (0.5Na2O 2CaO 1.5Al2O3 4FeO 30MgO 24SiO2) dat de samenstelling van het bulksilicaat aarde het beste benadert – het geochemische gemiddelde of onze hele planeet, behalve de op ijzer gebaseerde kern19. De vroege aarde werd gedomineerd door een reeks grootschalige smeltgebeurtenissen20, de laatste zou de hele planeet hebben overspoeld, na zijn condensatie voor de protolunar-schijf21. Pyroliet vertegenwoordigt de beste benadering van de samenstelling van dergelijke magma-oceanen op planetaire schaal. Daarom hebben we uitgebreid de fysische eigenschappen van pyrolietensmelt bestudeerd in het temperatuurbereik van 3.000\u20125.000 K en het drukbereik van 0\u2012150 GPa van ab initio moleculair-dynamische simulaties in de VASP-implementatie. Deze thermodynamische omstandigheden kenmerken volledig de meest extreme magma-oceaanomstandigheden op aarde. Onze studie is een uitstekend voorbeeld van een succesvol gebruik van het UMD-pakket voor de gehele diepgaande analyse van de melts22. We berekenden de verdeling en de gemiddelde bindingslengtes, we traceerden de veranderingen in kation-zuurstofcoördinatie en vergeleken onze resultaten met eerdere experimentele en computationele studies op amorfe silicaten van verschillende samenstellingen. Onze diepgaande analyse hielp bij het ontleden van standaardcoördinatienummers in hun basisbestanddelen, schetste de aanwezigheid van exotische coördinatiepolydra in de smelt en extractlevensduur voor alle coördinatiepolyhedra. Het schetste ook het belang van bemonstering in simulaties in termen van zowel de lengte van het traject als het aantal atomen dat aanwezig is in het systeem dat is gemodelleerd. Wat de nabewerking betreft, is de UMD-analyse onafhankelijk van deze factoren, maar ze moeten in aanmerking worden genomen bij de interpretatie van de resultaten van het UMD-pakket. Hier laten we een paar voorbeelden zien van hoe het UMD-pakket kan worden gebruikt om verschillende karakteristieke kenmerken van de smelten te extraheren, met een toepassing om pyroliet te smolten. De Si-O-paarverdelingsfunctie verkregen uit het gofrs_umd.py-script laat zien dat de straal van de eerste coördinatiebol, die het eerste minimum van de g(r)-functie is, rond 2,5 angstroms ligt bij T = 3000 K en P = 4,6 GPa. Het maximum van de g(r) ligt op 1,635 Å — dit is de beste benadering van de buiglengte. De lange staart is te wijten aan de temperatuur. Met behulp van deze limiet als de Si-O-bindingsafstand, laat de soortvormingsanalyse zien dat SiO4-eenheden, die tot enkele picoseconden kunnen duren, de smelt domineren (figuur 5). Er is een belangrijk deel van de smelt dat gedeeltelijke polymerisatie vertoont, zoals weerspiegeld door de aanwezigheid van dimeren zoals Si2O7 en trimers zoals Si3Ox-eenheden. Hun overeenkomstige levensduur is in de orde van de picoseconde. Polymeren van hogere orde hebben allemaal een aanzienlijk kortere levensduur. Figuur 5: Levensduur van de Si-O chemische soort.De soortvorming wordt geïdentificeerd in een meercomponentensmelt bij 4,6 GPa en 3000 K. De labels markeren de SiO3-, SiO4- en SiO5-monomeren en de verschillende SixOy-polymeren. Klik hier om een grotere versie van deze figuur te bekijken. De verschillende waarden van de verticale en horizontale stappen, gedefinieerd door de markeringen –z en –v hierboven, leveren verschillende steekproeven van de MSD op (figuur 6). Zelfs grote waarden van z en v zijn voldoende om de hellingen en dus de diffusiecoëfficiënten van de verschillende atomen te definiëren. De tijdwinst voor de nabewerking is opmerkelijk bij grote waarden van z en v. Het VIB biedt een zeer sterk validatiecriterium voor de kwaliteit van de simulaties. Als het diffusiegedeelte van het MSD niet lang genoeg is, is dat een teken dat de simulatie te kort is en de vloeistoftoestand in statistische zin niet bereikt. De minimumeis voor het diffusieve deel van het VIB is sterk afhankelijk van het systeem. Men kan eisen dat alle atomen minstens één keer hun plaats in de structuur van de smelt veranderen om het als een vloeistof te beschouwen10. Een uitstekend voorbeeld met toepassingen in de planetaire wetenschappen is complexe silicaatsmelt bij hoge druk dicht bij of zelfs onder hun liquiduslijn11. De Si-atomen, de belangrijkste netwerkvormende kationen, wisselen na meer dan twee dozijn picoseconden van plaats. Simulaties korter dan deze drempel zouden de mogelijke configuratieruimte aanzienlijk onderbemonsteren. Omdat de coördinerende anionen, namelijk de O-atomen, echter sneller bewegen dan de centrale Si-atomen, kunnen ze een deel van de langzame mobiliteit van Si compenseren. Als zodanig zou het hele systeem inderdaad een betere bemonstering van de configuratieruimte kunnen omvatten dan alleen van de Si-verplaatsingen wordt aangenomen. Figuur 6: Mean-square displacements (MSD).De MSD zijn geïllustreerd voor een paar atomaire typen van een meercomponenten silicaatsmelt. De bemonstering met verschillende horizontale en verticale stappen, z en v, levert consistente resultaten op. Massieve cirkels: -z 50 –v 50. Open cirkels: -z 250 –v 500. Klik hier om een grotere versie van deze figuur te bekijken. Ten slotte leveren de atomaire VAC-functies het trillingsspectrum van de smelt op. Figuur 7 toont het spectrum bij dezelfde druk- en temperatuuromstandigheden als hierboven. We vertegenwoordigen de bijdragen van Mg-, Si- en O-atomen, evenals de totale waarde. Bij nulfrequentie is er een eindige waarde van het spectrum, die overeenkomt met het diffusiekarakter van de smelt. Extractie van de thermodynamische eigenschappen uit het trillingsspectrum moet dit gasachtige diffusiekarakter van nul verwijderen, maar ook goed rekening houden met het verval bij hogere frequenties. Figuur 7: Het trillingsspectrum van pyrolietsmelt.Het reële deel van de Fouriertransformatie van de atomaire snelheid-snelheids-zelfcorrelatiefunctie levert het trillingsspectrum op. Hier wordt het spectrum berekend voor een meercomponenten silicaatsmelt. Vloeistoffen hebben een niet-nul gasachtig diffuus karakter bij nulfrequentie. Klik hier om een grotere versie van deze figuur te bekijken.

Discussion

Het UMD-pakket is ontworpen om beter te werken met ab initio-simulaties, waarbij het aantal snapshots meestal beperkt is tot tientallen tot honderdduizenden snapshots, met een paar honderd atomen per celeenheid. Grotere simulaties zijn ook handelbaar op voorwaarde dat de machine waarop de nabewerking draait over voldoende actieve geheugenbronnen beschikt. De code onderscheidt zich door de verscheidenheid aan eigenschappen die het kan berekenen en door zijn open-source licentie.

De umd.dat bestanden zijn geschikt voor de ensembles die het aantal deeltjes ongewijzigd gedurende de simulatie behouden. Het UMD-pakket kan bestanden lezen die afkomstig zijn van berekeningen waarbij de vorm en het volume van de simulatiebox variëren. Deze omvatten de meest voorkomende berekeningen, zoals NVT en NPT, waarbij het aantal deeltjes, N, temperatuur T, volume, V en / of druk, P, constant wordt gehouden.

Voor de tijd beginnen de paarverdelingsfunctie en alle scripts die de interatomaire afstanden moeten schatten, zoals de soortvormingsscripts, werken alleen voor orthogonale eenheidscellen, wat betekent voor kubische, tetragonale en orthorhombe cellen, waarbij de hoeken tussen de assen 90 ° zijn.

De belangrijkste ontwikkelingslijnen voor versie 2.0 zijn het verwijderen van de orthogonaliteitsbeperking voor afstanden en het toevoegen van meer functies voor de soortvormingsscripts: om individuele chemische bindingen te analyseren, om de interatomaire hoeken te analyseren en om de tweede coördinatiesfeer te implementeren. Met behulp van externe samenwerking werken we aan het porten van de code naar een GPU voor snellere analyse in grotere systemen.

Divulgazioni

The authors have nothing to disclose.

Acknowledgements

Dit werk werd ondersteund door de European Research Council (ERC) in het kader van het Horizon 2020 onderzoeks- en innovatieprogramma van de Europese Unie (subsidieovereenkomst nummer 681818 IMPACT to RC), door het Directoraat Extreme Physics and Chemistry van het Deep Carbon Observatory en door de Onderzoeksraad van Noorwegen via zijn Centers of Excellence-financieringsregeling, projectnummer 223272. We erkennen de toegang tot de GENCI-supercomputers via de stl2816-serie eDARI-computersubsidies, tot de Irene AMD-supercomputer via het PRACE RA4947-project en de Fram-supercomputer via de UNINETT Sigma2 NN9697K. FS werd ondersteund door een Marie Skłodowska-Curie-project (subsidieovereenkomst 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

Riferimenti

Frenkel, D., Smit, B. . Understanding Molecular Simulation. From Algorithms to Applications. , (2001).
Allen, M. P., Tildesley, D. J., Allen, T. . Computer Simulation of Liquids. , (1989).
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).
Humphrey, W., Dalke, A., Schulten, K. VMD: Visual molecular dynamics. Journal of Molecular Graphics & Modeling. 14 (1), 33-38 (1996).
Momma, K., Izumi, F. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. Journal of Applied Crystallography. 44 (6), 1272-1276 (2011).
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).
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).
Kresse, G., Hafner, J. Ab initio Molecular-Dynamics for Liquid-Metals. Physical Review B. 47 (1), 558-561 (1993).
Gygi, F. Architecture of Qbox: A scalable first-principles molecular dynamics code. IBM Journal of Research and Development. 52 (1-2), 137-144 (2008).
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).
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).
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).
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).
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).
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).
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).
Flyvbjerg, H., Petersen, H. G. Error estimates on averages of correlated data. The Journal of Chemical Physics. 91 (1), 461-466 (1989).
Tuckerman, M. E. . Statistical mechanics: theory and molecular simulation. , (2010).
McDonough, W. F., Sun, S. S. The composition of the Earth. Chemical Geology. 120, 223-253 (1995).
Elkins-Taton, L. T. Magma oceans in the inner solar system. Annual Review of Earth and Planetary Sciences. 40, 113-139 (2012).
Lock, S. J., et al. The origin of the Moon within a terrestrial synestia. J. Geophysical Research: Planets. 123, 910-951 (2018).
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

PDF

DOI

DOWNLOAD MATERIALS LIST

Citazione di questo articolo

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).