Summary
Smelter og væsker er allestedsnærværende vektorer av massetransport i naturlige systemer. Vi har utviklet en åpen kildekode-pakke for å analysere ab initio molekylære dynamikksimuleringer av slike systemer. Vi beregner strukturelle (binding, klynger, kjemisk spesiering), transport (diffusjon, viskositet) og termodynamiske egenskaper (vibrasjonsspektrum).
Abstract
Vi har utviklet en Python-basert åpen kildekode-pakke for å analysere resultatene som stammer fra ab initio molekylære dynamikksimuleringer av væsker. Pakken er best egnet for bruksområder på naturlige systemer, som silikat- og oksidsmelting, vannbaserte væsker og ulike superkritiske væsker. Pakken er en samling Python-skript som inkluderer to store biblioteker som omhandler filformater og krystallografi. Alle skriptene kjøres på kommandolinjen. Vi foreslår et forenklet format for å lagre atombanene og relevant termodynamisk informasjon om simuleringene, som er lagret i UMD-filer, som står for Universal Molecular Dynamics. UMD-pakken tillater beregning av en rekke strukturelle, transport- og termodynamiske egenskaper. Fra og med pardistribusjonsfunksjonen definerer den bindingslengder, bygger en interatomisk tilkoblingsmatrise og bestemmer til slutt den kjemiske spesieringen. Å bestemme levetiden til de kjemiske artene gjør det mulig å kjøre en full statistisk analyse. Deretter beregner dedikerte skript de gjennomsnittlige kvadratiske forskyvningene for atomene så vel som for de kjemiske artene. Den implementerte selvkorrelasjonsanalysen av atomhastighetene gir diffusjonskoeffisientene og vibrasjonsspekteret. Den samme analysen som brukes på påkjenningene, gir viskositeten. Pakken er tilgjengelig via GitHub-nettstedet og via sin egen dedikerte side i ERC IMPACT-prosjektet som åpen tilgangspakke.
Introduction
Væsker og smelter er aktive kjemiske og fysiske transportvektorer i naturlige miljøer. De forhøyede ratene av atomdiffusjon favoriserer kjemiske utvekslinger og reaksjoner, den lave viskositeten kombinert med varierende oppdrift favoriserer stor masseoverføring, og krystall-smeltetetthetsrelasjoner favoriserer lagdeling inne i planetlegemer. Fraværet av et periodisk gitter, typiske høye temperaturer som kreves for å nå smeltet tilstand, og vanskeligheten med slukking gjør den eksperimentelle bestemmelsen av en rekke åpenbare egenskaper, som tetthet, diffusjon og viskositet, ekstremt utfordrende. Disse vanskelighetene gjør alternative beregningsmetoder sterke og nyttige verktøy for å undersøke denne klassen av materialer.
Med ankomsten av datakraft og tilgjengeligheten av superdatamaskiner, brukes to store numeriske atomistiske simuleringsteknikker for tiden for å studere den dynamiske tilstanden til et ikke-krystallinsk atomistisk system, Monte Carlo1 og molekylær dynamikk (MD) 1,2. I Monte Carlo-simuleringer prøves konfigurasjonsrommet tilfeldig; Monte Carlo-metoder viser lineær skalering i parallellisering hvis alle prøvetakingsobservasjoner er uavhengige av hverandre. Kvaliteten på resultatene avhenger av kvaliteten på den tilfeldige tallgeneratoren og representativiteten til prøvetakingen. Monte Carlo-metoder viser lineær skalering parallelt hvis prøvetakingen er uavhengig av hverandre. I molekylær dynamikk (MD) prøves konfigurasjonsrommet av tidsavhengige atombaner. Fra en gitt konfigurasjon beregnes atombanene ved å integrere de newtonske bevegelsesligningene. De interatomære kreftene kan beregnes ved hjelp av modellinteratomiske potensialer (i klassisk MD) eller ved hjelp av første-prinsipper metoder (i ab initio, eller første-prinsipper, MD). Kvaliteten på resultatene avhenger av lengden på banen og dens evne til ikke å bli tiltrukket av lokal minima.
Molekylære dynamikksimuleringer inneholder en mengde informasjon, alt relatert til systemets dynamiske oppførsel. Termodynamiske gjennomsnittsegenskaper, som intern energi, temperatur og trykk, er ganske standard for beregning. De kan trekkes ut fra utgangsfilen(e) av simuleringene og gjennomsnittet, mens mengder relatert direkte til bevegelsen av atomene, samt til deres gjensidige forhold, må beregnes etter utvinning av atomposisjoner og hastigheter.
Følgelig har mye arbeid vært dedikert til å visualisere resultatene, og forskjellige pakker er tilgjengelige i dag, på forskjellige plattformer, åpen kildekode eller ikke [Ovito3, VMD4, Vesta5, Travis6, etc.]. Alle disse visualiseringsverktøyene håndterer effektivt med interatomære avstander, og som sådan tillater de effektiv beregning av parfordelingsfunksjoner og diffusjonskoeffisienter. Ulike grupper som utfører simuleringer av storskala molekylær dynamikk har proprietær programvare for å analysere ulike andre egenskaper som følge av simuleringene, noen ganger i shareware eller andre former for begrenset tilgang til samfunnet, og noen ganger begrenset i omfang og bruk til noen spesifikke pakker. Sofistikerte algoritmer for å trekke ut informasjon om interatomær binding, geometriske mønstre og termodynamikk utvikles og implementeres i noen av disse pakkene3,4,5,6,7, etc.
Her foreslår vi UMD-pakken - en åpen kildekode-pakke skrevet i Python for å analysere utgangen av molekylære dynamikksimuleringer. UMD-pakken muliggjør beregning av et bredt spekter av strukturelle, dynamiske og termodynamiske egenskaper (figur 1). Pakken er tilgjengelig via GitHub-nettstedet (https://github.com/rcaracas/UMD_package) og via en dedikert side (http://moonimpact.eu/umd-package/) av ERC IMPACT-prosjektet som en åpen tilgangspakke.
For å gjøre det universelt og enklere å håndtere, er vår tilnærming først å trekke ut all informasjon relatert til den termodynamiske tilstanden og atombanene fra utdatafilen til selve molekylærdynamikkløpet. Denne informasjonen lagres i en dedikert fil, hvis format er uavhengig av den opprinnelige MD-pakken der simuleringen ble kjørt. Vi navngir disse filene "umd" filer, som står for Universal Molecular Dynamics. På denne måten kan UMD-pakken vår enkelt brukes av enhver ab initio-gruppe med hvilken som helst programvare, alt med en minimal innsats for tilpasning. Det eneste kravet for å bruke den nåværende pakken er å skrive riktig parser fra utgangen av den aktuelle MD-programvaren til UMD-filformatet, hvis dette ennå ikke eksisterer. Foreløpig tilbyr vi slike parsere for VASP8- og QBox9-pakkene.
Figur 1: Flytskjema for UMD-biblioteket.
Fysiske egenskaper er i blått, og store Python-skript og deres alternativer er i rødt. Klikk her for å se en større versjon av denne figuren.
Umd-filene er ASCII-filer; typisk utvidelse er "umd.dat", men ikke obligatorisk. Alle analysekomponentene kan lese ASCII-filer i UMD-formatet, uavhengig av den faktiske navneutvidelsen. Noen av de automatiske skriptene som er utformet for å utføre rask statistikk i stor skala over flere simuleringer, ser imidlertid spesielt etter filer med filtypen UMD.dat. Hver fysiske egenskap uttrykkes på én linje. Hver linje starter med et nøkkelord. På denne måten er formatet svært tilpasningsdyktig og gjør det mulig å legge til nye egenskaper i UMD-filen, samtidig som det bevarer lesbarheten i alle versjoner. De første 30 linjene i umd-filen av simuleringen av pyrolitt ved 4,6 GPa og 3000 K, brukt nedenfor i diskusjonen, vises i figur 2.
Figur 2: Begynnelsen på umd-filen som beskriver simuleringen av flytende pyrolitt ved 4,6 GPa og 3000 K.
Hodet etterfølges av beskrivelsen av hvert øyeblikksbilde. Hver egenskap skrives på én linje, som inneholder navnet på den fysiske egenskapen, verdien(e) og enhetene, alle atskilt med mellomrom. Klikk her for å se en større versjon av denne figuren.
Alle UMD-filer inneholder en overskrift som beskriver innholdet i simuleringscellen: antall atomer, elektroner og atomtyper, samt detaljer for hvert atom, for eksempel type, kjemisk symbol, antall valenselektroner og dens masse. En tom linje markerer slutten på toppteksten og skiller den fra hoveddelen av UMD-filen.
Deretter er hvert trinn i simuleringen detaljert. For det første er de øyeblikkelige termodynamiske parametrene gitt, hver på en annen linje, som spesifiserer (i) navnet på parameteren, som energi, spenninger, tilsvarende hydrostatisk trykk, tetthet, volum, gitterparametere, etc., (ii) dens verdi (er) og (iii) dens enheter. Et bord som beskriver atomene kommer neste gang. En overskriftslinje gir de forskjellige målene, for eksempel kartesiske posisjoner, hastigheter, gebyrer osv., og deres enheter. Deretter er hvert atom detaljert på en linje. Ved grupper på tre, som tilsvarer de tre x-, y-, z-aksene , er oppføringene: de reduserte posisjonene, de kartesiske posisjonene brettet inn i simuleringscellen, de kartesiske posisjonene (som riktig tar hensyn til det faktum at atomer kan krysse flere enhetsceller under en simulering), atomhastighetene og atomkreftene. De to siste oppføringene er skalarer: ladning og magnetisk øyeblikk.
To store biblioteker sikrer at hele pakken fungerer som den skal. Det umd_process.py biblioteket omhandler UMD-filene, som lesing og utskrift. Det crystallography.py biblioteket omhandler all informasjon relatert til den faktiske atomstrukturen. Den underliggende filosofien til det crystallography.py biblioteket er å behandle gitteret som et vektorrom. Enhetscelleparameterne sammen med orienteringen representerer basisvektorene. "Space" har en rekke skalaregenskaper (spesifikt volum, tetthet, temperatur og spesifikt antall atomer), termodynamiske egenskaper (intern energi, trykk, varmekapasitet, etc.), og en rekke tensoriale egenskaper (stress og elastisitet). Atomer fyller dette rommet. "Gitter" -klassen definerer dette ensemblet, sammen med ulike få korte beregninger, som spesifikt volum, tetthet, oppnå det gjensidige gitteret fra den direkte, etc. "Atomer" -klassen definerer atomene. De er preget av en rekke skalaregenskaper (navn, symbol, masse, antall elektroner, etc.) og en rekke vektoregenskaper (posisjonen i rommet, enten i forhold til vektorgrunnlaget beskrevet i Gitterklassen, eller i forhold til universelle kartesiske koordinater, hastigheter, krefter, etc.). Bortsett fra disse to klassene, inneholder det crystallography.py biblioteket en rekke funksjoner for å utføre en rekke tester og beregninger, for eksempel atomavstander eller cellemultiplikasjon. Den periodiske tabellen over elementene er også inkludert som en ordbok.
De ulike komponentene i UMD-pakken skriver flere utdatafiler. Som en generell regel er de alle ASCII-filer, alle oppføringene deres er skilt av faner, og de er laget så selvforklarende som mulig. For eksempel angir de alltid tydelig den fysiske egenskapen og dens enheter. Filene umd.dat er fullstendig i samsvar med denne regelen.
Protocol
1. Analyse av molekylærdynamikken går
MERK: Pakken er tilgjengelig via GitHub-nettstedet (https://github.com/rcaracas/UMD_package) og via en dedikert side (http://moonimpact.eu/umd-package/) av ERC IMPACT-prosjektet som en åpen tilgangspakke.
- Pakk ut hvert enkelt sett med fysiske egenskaper ved hjelp av ett eller flere dedikerte Python-skript fra pakken. Kjør alle skriptene på kommandolinjen. De bruker alle en rekke flagg, som er så konsekvente som mulig fra ett skript til et annet. Flaggene, betydningen og standardverdiene er alle oppsummert i tabell 1.
| Flagg | Betydning | Skript som bruker det | Standardverdi |
| -h | Kort hjelp | alt | |
| -f | UMD filnavn | alt | |
| -Jeg | Termiske trinn som skal kastes | alt | 0 |
| -Jeg | Inndatafil som inneholder de interatomære bindingene | Artsdannelse | obligasjoner.input |
| -s | Prøvetaking av frekvensen | msd, spesiering | 1 (hvert trinn vurderes) |
| -a | Liste over atomer eller anioner | Artsdannelse | |
| -c | Liste over kasjoner | Artsdannelse | |
| -Jeg | Binding lengde | Artsdannelse | 2 |
| -t | Temperatur | vibrasjoner, reologi | |
| -v | Diskretisering av bredden på prøvetakingsvinduet for banen for middelverdiforskyvningsanalysen | Msd | 20 |
| -å | Diskretisering av starten på prøvetakingsvinduet for banen for middelverdiforskyvningsanalysen | Msd | 20 |
Tabell 1: De vanligste flaggene som brukes i UMD-pakken og deres vanligste betydning.
- Begynn med å transformere utgangen av MD-simuleringen som utføres i en førsteprinsipperskode, som VASP8 eller QBox9, til en UMD-fil.
- Hvis MD-simuleringene ble gjort i VASP, skriver du inn følgende på kommandolinjen:
VaspParser.py -f-i
der –f-flagget definerer navnet på VASP OUTCAR-filen, og –i termiseringslengden.
MERK: Det første trinnet, definert av –i, gjør det mulig å kaste de første trinnene i simuleringene, som representerer termaliseringen. I et typisk molekylært dynamikkløp representerer den første delen av beregningen termiseringen, det vil si tiden det tar systemet for alle atomer å beskrive en gaussisk-lignende fordeling av temperaturen, og for hele systemet å vise svingninger i temperatur, trykk, energi, etc. rundt likevektsverdier. Denne termiske delen av simuleringen bør ikke tas i betraktning når du analyserer væskens statistiske egenskaper.
- Hvis MD-simuleringene ble gjort i VASP, skriver du inn følgende på kommandolinjen:
- Transformer . UMD-filer i . xyz-filer for å lette visualisering på forskjellige andre pakker, som VMD4 eller Vesta5. Ved kommandolinjetypen:
umd2xyz.py -f-i -s
der –f definerer navnet på . umd-fil , –i definerer termiseringsperioden som skal kastes, og –er frekvensen av prøvetakingen av banen som er lagret i . UMD-fil . Standardverdiene er –i 0 –s 1, det vil si med tanke på alle trinnene i simuleringen, uten at noen forkastes. - Reverser UMD-filen til POSCAR-filer av TYPEN VASP ved hjelp av umd2poscar.py skriptet. øyeblikksbilder av simuleringene kan velges med en forhåndsdefinert frekvens. Ved kommandolinjetypen:
umd2poscar.py -f-i -l -s
der –l representerer det siste trinnet som skal gjøres om til POSCAR-fil. Standardverdier er -i 0 -l 100000000 -s 1. Denne verdien av –l er stor nok til å dekke en typisk hel bane.
2. Utføre strukturanalysen
- Kjør skriptet gofrs_umd.py for å beregne pardistribusjonsfunksjonen (PDF) g ᴀʙ(r) for alle parene atomtyper A og B (figur 3). Utdataene er skrevet i en ASCII-fil, tab-separert, med utvidelsen gofrs.dat. Ved kommandolinjetypen:
gofrs_umd.py -f-s < Sampling_Frequency > -d -i
MERK: Standardene er Sampling_Frequency (frekvensen for prøvetaking av banen) = 1 trinn; DiskretiseringInterval (for plotting av g(r)) = 0,01 Å; InitialStep (antall trinn i begynnelsen av banen som forkastes) = 0. Radial PDF, g ᴀʙ(r) er gjennomsnittlig antall atomer av type B på avstand d_ᴀʙ innenfor et sfærisk skall av radius r og tykkelse dr sentrert på atomer av type A (figur 3):
med ρ atomtetthet, NA og NB antall atomer av type A og B, og δ(r−r ᴀʙ) deltafunksjonen som er lik 1 hvis atomene A og B ligger i avstand mellom r og r +dr. Abscissa av det første maksimumet av g ᴀʙ(r) gir den høyeste sannsynlighetsbindingslengden mellom atomene av type A og B, som er nærmest en gjennomsnittlig bindingsavstand som vi kan bestemme. Det første minimum avgrenser omfanget av den første koordineringssfæren. Derfor gir integralet over PDF-filen opp til det første minimum det gjennomsnittlige koordineringsnummeret. Summen av Fouriertransformasjonene av g ᴀʙ(r) for alle par atomtyper A og B gir væskens diffraksjonsmønster, som oppnådd eksperimentelt med et diraktometer. Men i virkeligheten, så ofte de høye ordens koordineringssfærene mangler fra g ᴀʙ(r), kan ikke diffraksjonsmønsteret oppnås i sin helhet.
Figur 3: Bestemmelse av parfordelingsfunksjon.
(a) For hvert atom av en art (for eksempel rød), telles alle atomene til koordinerende arter (for eksempel grå og / eller rød) som en funksjon av avstand. (b) Den resulterende avstandsfordelingsgrafen for hvert øyeblikksbilde, som på dette stadiet bare er en samling deltafunksjoner, blir deretter gjennomsnittet over alle atomene og alle øyeblikksbildene og vektet av den ideelle gassfordelingen for å generere (c) parfordelingsfunksjonen som er kontinuerlig. Det første minimumet av g(r) er radiusen til den første koordineringssfæren, som senere brukes i spesieringsanalysen. Klikk her for å se en større versjon av denne figuren.
- Trekk ut de gjennomsnittlige interatomære bindingsavstandene som radier av de første koordineringssfærene. For dette, identifiser plasseringen av det første maksimumet av g ᴀʙ(r)-funksjonene: plott gofrs.dat-filen i et regnearkprogram og søk etter maxima og minima for hvert par atomer.
- Identifiser radiusen til den første koordineringssfæren, som det første minimumet av PDF ᴀʙ(r), ved hjelp av regnearkprogramvare. Dette er grunnlaget for hele strukturanalysen av væsken; PDF-filen gir den gjennomsnittlige bindingsstatusen til atomene i væsken.
- Trekk ut avstandene til den første minima, det vil si abscissa, og skriv dem i en egen fil, kalt for eksempel bonds.input. Alternativt kan du kjøre et av de analyze_gofr skriptene i UMD-pakken for å identifisere maxima- og minima ᴀʙ(r)-funksjonene. Ved kommandolinjetypen:
analyze_gofr_semi_automatic.py - Klikk på posisjonen til maksimum og minimum av g ᴀʙ(r)-funksjonen som vises i grafen som åpnes av programmet. Skriptet skanner automatisk gjeldende mappe, identifiserer alle gofrs.dat filer og utfører analysen for hver enkelt av dem. Klikk igjen på maksimum og minimum i vinduet hver gang skriptet trenger en utdannet innledende gjetning.
- Åpne og se på den automatisk genererte filen kalt bonds.input som inneholder de interatomære bindingsavstandene.
3. Utfør spesieringsanalysen
- Beregn topologien for binding mellom atomene, ved hjelp av begrepet tilkobling innen grafteori: atomene er nodene og de interatomære bindingene er banene. Det speciation_umd.py skriptet trenger de interatomære bindingsavstandene som er definert i filen bonds.input .
MERK: Tilkoblingsmatrisen er konstruert ved hvert trinn: To atomer som ligger på en avstand som er mindre enn radiusen til den tilsvarende første koordineringssfæren, anses å være bundet, det vil si tilkoblet. Ulike atomnettverk er bygget ved å behandle atomene som noder i en graf hvis forbindelser er definert av dette geometriske kriteriet. Disse nettverkene er atomartene, og deres ensemble definerer atomspesifikasjonen i den aktuelle væsken (figur 4).
Figur 4: Identifikasjon av atomklyngene.
Koordineringspolyhedra er definert ved hjelp av interatomære avstander. Alle atomer i en avstand som er mindre enn en bestemt radius, anses å være bundet. Her tilsvarer terskelen den første koordineringssfæren (de lyserøde sirklene), definert i figur 1. Polymerisasjon og dermed kjemiske arter er hentet fra nettverkene til de limte atomene. Legg merke til den sentrale Red1Grey2-klyngen, som er isolert fra de andre atomene, som danner en uendelig polymer. Klikk her for å se en større versjon av denne figuren.
- Kjør spesieringsskriptet for å få tilkoblingsmatrisen og få koordineringspolyhedraen eller polymeriseringen. Ved kommandolinjetypen:
speciation_umd.py -f-s -i -l -c -a -m -r
der -i-flagget gir filen de interatomære bindingsavstandene, som for eksempel ble produsert i forrige trinn. Du kan også kjøre skriptet med én enkelt lengde for alle obligasjoner som er definert av -l -flagget.
MERK: -c-flagget angir de sentrale atomene, og -a flagger ligandene. Både sentrale atomer og ligander kan være av forskjellige typer; I dette tilfellet må de skilles med komma. -m-flagget gir minimumstiden en art må leve for å bli vurdert i analysen. Som standard er denne minimumstiden null, alle forekomster telles i den endelige analysen.- Kjør skriptet speciation_umd.py med flagget –r 0, som tar prøver av tilkoblingsgrafen på første nivå for å identifisere koordineringspolyhedraen. Et sentralt atom, betegnet som en kation , kan for eksempel være omgitt av én eller flere anioner (figur 4). Spesieringsskriptet identifiserer hver eneste av koordineringspolyhedraen. Det vektede gjennomsnittet av all koordineringspolyhedra gir koordineringsnummeret, identisk med det som er oppnådd ved integrasjonen av PDF-filen. Ved kommandolinjetypen:
speciation_umd.py -f-i -c -a -r 0
MERK: Gjennomsnittlige koordinasjonstall i væsker er brøktall. Denne fraksjonaliteten kommer fra den gjennomsnittlige egenskapen til koordineringen. Definisjonen basert på spesiering gir en mer intuitiv og informativ representasjon av væskens struktur, hvor de relative proporsjonene av de forskjellige artene, det vil si koordinasjoner, kvantifiseres. - Kjør speciation_umd.py skript med flagget –r 1, som tar prøver av tilkoblingsgrafen på alle dybdenivåer for å oppnå polymerisasjonen. Nettverket gjennom atomgrafen har en viss dybde, da atomer er bundet lenger bort til andre bindinger (f.eks. i sekvenser av vekslende kasjoner og anioner) (figur 4).
- Kjør skriptet speciation_umd.py med flagget –r 0, som tar prøver av tilkoblingsgrafen på første nivå for å identifisere koordineringspolyhedraen. Et sentralt atom, betegnet som en kation , kan for eksempel være omgitt av én eller flere anioner (figur 4). Spesieringsskriptet identifiserer hver eneste av koordineringspolyhedraen. Det vektede gjennomsnittet av all koordineringspolyhedra gir koordineringsnummeret, identisk med det som er oppnådd ved integrasjonen av PDF-filen. Ved kommandolinjetypen:
- Åpne de to filene . popul.dat og . statistikk.dat fortløpende; Disse utgjør utdataene fra spesieringsskriptet. Hver klynge er skrevet på en linje, og spesifiserer sin kjemiske formel, tidspunktet da den dannet seg, tidspunktet da den døde, dens levetid, en matrise med listen over atomer som danner denne klyngen. Plott levetiden til hver atomklynge av alle de kjemiske artene som finnes i simuleringen som finnes i .popul.dat-filen (figur 5).
- Plott befolkningsanalysen med overflod av hver art, som finnes i . stat.dat fil. Denne analysen, både absolutt og relativ, tilsvarer den faktiske statistikken over koordineringspolyhedraen for saken -r 0; For tilfelle av polymerisasjon, med -r 1 må dette behandles nøye, da noe normalisering over det relative antallet atomer kanskje må påføres. Overfloden tilsvarer integralet i løpet av levetiden. Den. stat.dat filen viser også størrelsen på hver klynge, det vil si hvor mange atomer som danner den.
4. Beregne diffusjonskoeffisienter
- Trekk ut atomenes gjennomsnittlige firkantede forskyvninger (MSD) som en funksjon av tid for å oppnå selvspredning. Standardformelen for MSD er:
der prefaktorene er renormaliseringer. Med MSD-verktøyet er det forskjellige måter å analysere de dynamiske aspektene ved væskene på.
MERK: T er den totale tiden for simuleringen, og N α er antall atomer av type α. Den første gangen t0 er vilkårlig og strekker seg over første halvdel av simuleringen. Ninit er antall innledende tidspunkt. τ er bredden på tidsintervallet som MSD beregnes over. Maksimumsverdien er halvparten av simuleringens tidslengde. I vanlige MSD-implementeringer starter hvert vindu på slutten av det forrige. Men en sparsommelig prøvetaking kan øke beregningen av MSD, uten å endre skråningen av MSD. For dette starter det i-th-vinduet på tidspunktet t0(i), men vinduet (i+1)-th starter på tidspunkt t0(i) + τ + v, der verdien av v er brukerdefinert. På samme måte økes bredden på vinduet i atskilte trinn som er definert av brukeren, for eksempel: τ(i) = τ(i-1) + z. Verdiene for z ("horisontalt trinn") og v ("vertikalt trinn") er positive eller null; Standarden for begge er 20. - Beregn MSD ved hjelp av serien med msd_umd skript. Utdataene skrives ut i en . msd.dat fil, der MSD for hver atomtype, atom eller klynge skrives ut på én kolonne som en tidsfunksjon.
- Beregn gjennomsnittlig MSD for hver atomtype. MSD beregnes for hvert atom og beregnes deretter for hver atomtype. Utdatafilen inneholder én kolonne for hver atomtype. Ved kommandolinjetypen:
msd_umd.py -f-z -v -b - Beregn MSD for hvert atom. MSD beregnes for hvert atom og beregnes deretter for hver atomtype. Utdatafilen inneholder en kolonne for hvert atom i simuleringen, og deretter en kolonne for hver atomtype. Denne funksjonen gjør det mulig å identifisere atomer som sprer seg i to forskjellige miljøer, som væske og gass, eller to væsker. Ved kommandolinjetypen:
msd_all_umd.py -f-z -v -b - Beregn MSD for de kjemiske artene. Bruk populasjonen av klynger som er identifisert med spesieringsskriptet, og skriv ut i . popul.dat fil. MSD beregnes for hver enkelt klynge. Utdatafilen inneholder én kolonne for hver klynge. For å unngå å vurdere store polymerer, plasser en grense på størrelsen på klyngen; Standarden er 20 atomer. Ved kommandolinjetypen:
msd_cluster_umd.py -f-p -s -b -c
MERK: Standardverdier er: –b 100 –s 1 –c 20.
- Beregn gjennomsnittlig MSD for hver atomtype. MSD beregnes for hvert atom og beregnes deretter for hver atomtype. Utdatafilen inneholder én kolonne for hver atomtype. Ved kommandolinjetypen:
- Tegn inn MSD ved hjelp av en regnearkbasert programvare (figur 6). I en loggloggrepresentasjon av MSD i forhold til tid identifiserer du stigningstallet. Skill den første delen, vanligvis kort, som representerer det ballistiske regimet, det vil si bevaring av hastigheten på atomer etter kollisjoner. Den andre lengre delen representerer det diffusive regimet, det vil si spredning av hastigheten på atomer etter kollisjoner.
- Beregn diffusjonskoeffisientene fra stigningstallet til MSD som:
hvor Z er antall frihetsgrader (Z = 2 for diffusjon i plan, Z = 3 for diffusjon i rommet), og t er tidstrinnet.
5. Tidskorrelasjonsfunksjoner
- Beregn tidskorrelasjonsfunksjonene som et mål på tregheten til systemet ved hjelp av den generelle formelen:
En kan være en rekke tidsavhengige variabler, for eksempel atomposisjoner, atomhastigheter, spenninger, polarisering, etc., som hver gir – via Green-Kubo-relasjonene12,13 – forskjellige fysiske egenskaper, noen ganger etter en ytterligere transformasjon. - Analyser atomhastighetene for å oppnå vibrasjonsspekteret til væsken og det alternative uttrykket av de atomiske selvspredningskoeffisientene.
- Kjør vibr_spectrum_umd.py skript for å beregne funksjonen for automatisk korrelasjon med atomhastighet (VAC) for hver atomtype og utføre den raske Fourier-transformeringen. Ved kommandolinjetypen:
vibr_spectrum_umd.py -f-t
der –t er temperaturen som må defineres av brukeren. Skriptet skriver ut to filer: . vels.scf.dat fil med VAC-funksjonen for hver atomtype, og . vibr.dat fil med vibrasjonsspekteret dekomponert på hver atomart og den totale verdien. - Åpne og les vels.scf.dat. Plott VAC-funksjonen fra filen vels.scf.dat ved hjelp av regnearklignende programvare.
- Behold den virkelige delen av Fourier VAC. Dette er det som gir vibrasjonsspekteret, som en funksjon av frekvens:
hvor m er atommassene. - Plott det vibrasjonsmessige spekteret fra vibr.dat-filen ved hjelp av regnearklignende programvare (figur 7). Identifiser den endelige verdien ved ω = 0 som tilsvarer væskens diffusive karakter og de forskjellige toppene i spekteret ved begrenset frekvens. Identifiser deltakelsen av hver atomtype til vibrasjonsspekteret.
MERK: Nedbrytningen på atomtyper viser at forskjellige atomer har forskjellige ω = 0-bidrag, som tilsvarer deres diffusjonskoeffisienter. Den generelle formen på spekteret er mye jevnere med færre funksjoner enn for en tilsvarende solid. - På skallet, les integralet over vibrasjonsspekteret, som gir diffusjonskoeffisientene for hver atomart.
MERK: Termodynamiske egenskaper kan oppnås ved integrasjon fra vibrasjonsspekteret, men resultatene bør brukes med forsiktighet på grunn av to tilnærminger: integrasjonen er gyldig innenfor kvasi-harmonisk tilnærming, som ikke nødvendigvis holder ved høye temperaturer; og den gasslignende delen av spekteret som tilsvarer diffusjonen må kastes. Integrasjonen skal da bare gjøres over den gitterlignende delen av spekteret. Men denne separasjonen krever vanligvis flere ytterligere etterbehandlingstrinn og beregninger14, som ikke dekkes av den nåværende UMD-pakken.
- Kjør vibr_spectrum_umd.py skript for å beregne funksjonen for automatisk korrelasjon med atomhastighet (VAC) for hver atomtype og utføre den raske Fourier-transformeringen. Ved kommandolinjetypen:
- Kjør viscosity_umd.py skript for å analysere selvkorrelasjonen til komponentene stress tensor for å estimere viskositeten til smelten. Ved kommandolinjetypen:
viscosity_umd.py -f-i -s -o -l
MERK: Denne funksjonen er utforskende, og eventuelle resultater må tas med forsiktighet. For det første, kontroller konvergensen av viskositeten grundig med hensyn til lengden på simuleringen.- Utlede viskositeten til væsken fra selvkorrelasjonen av stresset tensor15 som:
der V og T er henholdsvis volumet og temperaturen, κB er Boltzmann konstant og σ ij den ij off-diagonal komponenten av stress-tensor, uttrykt i kartesiske koordinater. - Bruk en mer tilstrekkelig passform for å oppnå et mer robust estimat av viskositeten15,16 og unngå støyen fra stress-tensor auto-korrelasjonsfunksjonen som kan oppstå fra den endelige størrelsen og den begrensede varigheten av simuleringene. For den automatiske korrelasjonsfunksjonen til spennings tensoren, bruk følgende funksjonelle form15,16 som gir gode resultater:
der A, B, τ1, τ2 og ω er tilpasningsparametere. Etter integrering blir uttrykket for viskositeten:
- Utlede viskositeten til væsken fra selvkorrelasjonen av stresset tensor15 som:
6. Termodynamiske parametere som stammer fra simuleringene.
- Kjør averages.py for å trekke ut gjennomsnittsverdiene og spredningen (som standardavvik) for trykk, temperatur, tetthet og intern energi fra UMD-filene . Ved kommandolinjetypen:
averages.py -f-s
med –s 0 som standard. - Beregn den statistiske feilen i gjennomsnittet ved hjelp av blokkeringsmetoder.
MERK: Det er forskjellige smaker av denne metoden. Etter Allens og Tildesley2s arbeid er det vanlig å gjennomsnitte over sekvenser av tidsblokker, med stadig lengre lengde, og estimere standardavviket med hensyn til det aritmetiske gjennomsnittet17. Konvergens kan nås på grensen mellom mange og lange nok blokkstørrelser, når prøvetakingen ikke er korrelert. Selv om den faktiske terskelverdien for konvergens vanligvis må velges manuelt.- Bruk halveringsmetoden18: Fra og med det første dataeksemplet halverer du antallet prøver ved å beregne gjennomsnittet for hvert annet tilsvarende påfølgende utvalg fra forrige trinn κ−1:
- Kjør skriptet fullaverages.py for å utføre hele den statistiske analysen, inkludert feilen i gjennomsnittet. Ved kommandolinjetypen:
fullaverages.py -s-u
MERK: Skriptet automatiseres til det punktet å søke etter alle .umd.dat filer i gjeldende katalog og utføre analysen for dem alle. Standard er –s 0 –u 0. For -u 0 er utdataene minimale, og for -u 1 er utskriften i sin helhet, med flere alternative enheter skrevet ut. Dette skriptet krever grafisk støtte, da det oppretter et grafisk bilde for å kontrollere konvergensen for å estimere feilen på gjennomsnittet.
- Bruk halveringsmetoden18: Fra og med det første dataeksemplet halverer du antallet prøver ved å beregne gjennomsnittet for hvert annet tilsvarende påfølgende utvalg fra forrige trinn κ−1:
Representative Results
Pyrolitt er en modell multikomponent silikatsmelting (0,5Na2O 2CaO 1,5Al2O3 4FeO 30MgO 24SiO2) som best tilnærmer seg sammensetningen av bulksilikatjorden – det geokjemiske gjennomsnittet eller hele planeten vår, bortsett fra den jernbaserte kjernen19. Den tidlige jorden var dominert av en rekke store smeltehendelser20, den siste kunne ha oppslukt hele planeten, etter kondensasjonen for protolunardisken21. Pyrolitt representerer den beste tilnærmingen til sammensetningen av slike planetariske magmahav. Følgelig studerte vi i stor grad de fysiske egenskapene til pyrolittsmelting i 3000 \u20125,000 K temperaturområdet og 0\u2012150 GPa-trykkområdet fra ab initio molekylære dynamikksimuleringer i VASP-implementeringen. Disse termodynamiske forholdene karakteriserer helt jordens mest ekstreme magma havforhold. Vår studie er et utmerket eksempel på en vellykket bruk av UMD-pakken for hele grundig analyse av smeltene22. Vi beregnet fordelingen og den gjennomsnittlige bindingslengden, vi sporet endringene i kation-oksygenkoordinering, og sammenlignet resultatene våre med tidligere eksperimentelle og beregningsstudier på amorfe silikater av ulike komposisjoner. Vår grundige analyse bidro til å dekomponere standard koordinasjonstall i sine grunnleggende bestanddeler, skissere tilstedeværelsen av eksotisk koordineringspolyhedra i smelten, og trekke ut levetider for all koordineringspolyhedra. Den skisserte også viktigheten av prøvetaking i simuleringer både når det gjelder lengden på banen og også antall atomer som finnes i systemet som er modellert. Når det gjelder etterbehandling, er UMD-analysen uavhengig av disse faktorene, men de bør tas i betraktning når de tolker resultatene fra UMD-pakken. Her viser vi noen eksempler på hvordan UMD-pakken kan brukes til å trekke ut flere karakteristiske trekk ved smeltene, med et program for smeltet pyrolitt.
Si-O-parfordelingsfunksjonen hentet fra gofrs_umd.py skript viser at radiusen til den første koordineringssfæren, som er det første minimumet av g(r)-funksjonen, ligger rundt 2,5 angstrom ved T = 3000 K og P = 4,6 GPa. Maksimum av g(r) ligger på 1.635 Å – dette er den beste tilnærmingen til bøyelengden. Den lange halen skyldes temperaturen. Ved å bruke denne grensen som Si-O-bindingsavstand, viser spesieringsanalysen at SiO4-enheter, som kan vare i opptil noen få picoseconds, dominerer smelten (figur 5). Det er en viktig del av smelten som viser delvis polymerisering, som reflektert av tilstedeværelsen av dimmere som Si2O7, og trimers som Si3Ox-enheter. Deres tilsvarende levetid er i picoseconds rekkefølge. Høyere orden polymerer har alle betydelig kortere levetid.
Figur 5: Levetiden til si-O-kjemiske arter.
Spesieringen er identifisert i en flerkomponentsmelting ved 4,6 GPa og 3000 K. Etikettene markerer SiO3-, SiO4- og SiO5-monomerene og de forskjellige SixOy-polymerene. Klikk her for å se en større versjon av denne figuren.
De ulike verdiene i de loddrette og vannrette trinnene, definert av flaggene –z og –v ovenfor, gir ulike utvalg av MSD (figur 6). Selv store verdier av z og v er nok til å definere bakkene og dermed diffusjonskoeffisientene til de forskjellige atomene. Gevinsten i tide for etterbehandling er bemerkelsesverdig når du går til store verdier av z og v. Msd tilbyr et veldig sterkt valideringskriterium for kvaliteten på simuleringene. Hvis diffusjonsdelen av MSD ikke er tilstrekkelig lang, er det et tegn på at simuleringen er for kort, og ikke når væsketilstanden i statistisk forstand. Minimumskravet for den diffusive delen av MSD avhenger sterkt av systemet. Man kan kreve at alle atomene endrer sitt sted minst en gang i smeltestrukturen for at det skal betraktes som en væske10. Et utmerket eksempel med applikasjoner i planetariske vitenskaper er komplekse silikat smelter ved høyt trykk nær eller til og med under deres liquidus line11. Si-atomene, de store nettverksdannende kasjonene, bytter nettsteder etter mer enn to dusin picoseconds. Simuleringer som er kortere enn denne terskelen, vil være betydelig underprøve av det mulige konfigurasjonsområdet. Men som de koordinerende anionene, nemlig O-atomene, beveger seg raskere enn de sentrale Si-atomene, kan de kompensere for en del av siens langsomme mobilitet. Som sådan kunne hele systemet faktisk dekke et bedre utvalg av konfigurasjonsrommet enn antatt bare fra Si-forskyvningene.
Figur 6: Middelverdiforskyvninger (MSD).
Msd er illustrert for noen atomiske typer av en multi-komponent silikat smelte. Prøvetakingen med ulike horisontale og vertikale trinn, z og v, gir konsistente resultater. Heldekkende sirkler: -z 50 –v 50. Åpne sirkler: -z 250 –v 500. Klikk her for å se en større versjon av denne figuren.
Til slutt gir de atomiske VAC-funksjonene det vibrasjonsmessige spekteret av smelten. Figur 7 viser spekteret ved samme trykk- og temperaturforhold som ovenfor. Vi representerer bidragene fra Mg-, Si- og O-atomer, samt den totale verdien. Ved null frekvens er det en begrenset verdi av spekteret, som tilsvarer diffusjonskarakteren til smelten. Ekstraksjon av de termodynamiske egenskapene fra vibrasjonsspekteret må fjerne denne gasslignende diffusive karakteren fra null, men også å ta hensyn til forfallet ved høyere frekvenser.
Figur 7: Vibrasjonsspekteret av pyrolitt smelter.
Den virkelige delen av Fourier-transformasjonen av selvkorrelasjonsfunksjonen for atomhastighet gir vibrasjonsspekteret. Her beregnes spekteret for en flerkomponent silikatsmelting. Væsker har en ikke-null gasslignende diffusiv karakter ved null frekvens. Klikk her for å se en større versjon av denne figuren.
Discussion
UMD-pakken er designet for å fungere bedre med ab initio-simuleringer, der antall øyeblikksbilder vanligvis er begrenset til titusenvis av øyeblikksbilder, med noen få hundre atomer per enhetscelle. Større simuleringer kan også sendes til maskinen som etterbehandlingen kjører på, har nok aktive minneressurser. Koden skiller seg ut med de forskjellige egenskapene den kan beregne og ved åpen kildekode-lisensen.
Umd.dat filene passer til ensemblene som bevarer antall partikler uendret gjennom hele simuleringen. UMD-pakken kan lese filer som stammer fra beregninger der formen og volumet på simuleringsboksen varierer. Disse dekker de vanligste beregningene, som NVT og NPT, der antall partikler, N, temperatur T, volum, V og/eller trykk, P, holdes konstant.
For tiden begynner pardistribusjonsfunksjonen så vel som alle skriptene som trenger å estimere de interatomære avstandene, som spesieringsskriptene, fungerer bare for ortogonale enhetsceller, noe som betyr for kubikkceller, tetragonale og orthorhombic celler, hvor vinklene mellom aksene er 90°.
De viktigste utviklingslinjene for versjon 2.0 er fjerning av ortogonalitetsbegrensningen for avstander og legge til flere funksjoner for spesieringsskriptene: å analysere individuelle kjemiske bindinger, for å analysere de interatomære vinklene, og for å implementere den andre koordineringssfæren. Med hjelp fra eksternt samarbeid jobber vi med å overføre koden til en GPU for raskere analyse i større systemer.
Disclosures
Forfatterne har ingenting å avsløre.
Acknowledgments
Dette arbeidet ble støttet av Det europeiske forskningsrådet (ERC) under Forsknings- og innovasjonsprogrammet European Union Horizon 2020 (tilskuddsavtalenummer 681818 IMPACT to RC), av Extreme Physics and Chemistry Directorate of the Deep Carbon Observatory, og av Norges forskningsråd gjennom sin sentre for fremragende forskning, prosjektnummer 223272. Vi anerkjenner tilgang til GENCI-superdatamaskinene gjennom stl2816-serien med eDARI-databehandlingsstipender, til Irene AMD-superdatamaskinen gjennom PRACE RA4947-prosjektet, og Fram-superdatamaskinen gjennom UNINETT Sigma2 NN9697K. FS ble støttet av et Marie Skłodowska-Curie-prosjekt (tilskuddsavtale ABISSE Nr.750901).
Materials
| Name | Company | Catalog Number | Comments |
| 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
- Frenkel, D., Smit, B. Understanding Molecular Simulation. From Algorithms to Applications. , Elsevier. (2001).
- Allen, M. P., Tildesley, D. J., Allen, T. Computer Simulation of Liquids. , Oxford University Press. (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.
- 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. , Oxford University Press. (2010).
- McDonough, W. F., Sun, S. S.
- 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).
