Summary
Smälter och vätskor är allestädes närvarande vektorer för masstransport i naturliga system. Vi har utvecklat ett paket med öppen källkod för att analysera ab initio molekylär-dynamics simuleringar av sådana system. Vi beräknar strukturella (bindning, klusterisering, kemisk speciation), transport (diffusion, viskositet) och termodynamiska egenskaper (vibrationsspektrum).
Abstract
Vi har utvecklat ett Python-baserat open source-paket för att analysera resultaten som härrör från ab initio molekylär-dynamik simuleringar av vätskor. Paketet är bäst lämpat för tillämpningar på naturliga system, som silikat och oxid smälter, vattenbaserade vätskor och olika superkritiska vätskor. Paketet är en samling Python-skript som innehåller två stora bibliotek som behandlar filformat och kristallografi. Alla skript körs på kommandoraden. Vi föreslår ett förenklat format för att lagra atombanorna och relevant termodynamisk information om simuleringarna, som sparas i UMD-filer, stående för Universal Molecular Dynamics. UMD-paketet möjliggör beräkning av en rad strukturella, transport- och termodynamiska egenskaper. Från och med parfördelningsfunktionen definierar den bindningslängder, bygger en interatomär anslutningsmatris och bestämmer så småningom den kemiska speciationen. Att bestämma den kemiska artens livslängd gör det möjligt att göra en fullständig statistisk analys. Sedan beräknar dedikerade skript medelvärdesförskjutningarna för atomerna såväl som för den kemiska arten. Den implementerade självkorrelationsanalysen av atomhastigheterna ger diffusionskoefficienterna och vibrationsspektrumet. Samma analys som tillämpas på påfrestningarna ger viskositeten. Paketet är tillgängligt via GitHubs webbplats och via sin egen dedikerade sida i ERC IMPACT-projektet som open access-paket.
Introduction
Vätskor och smälter är aktiva kemiska och fysiska transportvektorer i naturliga miljöer. De förhöjda frekvenserna av atomisk diffusion gynnar kemiska utbyten och reaktioner, den låga viskositeten i kombination med varierande flytkraft gynnar stor massöverföring, och kristallsmälttäthetsrelationer gynnar skiktning inuti planetariska kroppar. Frånvaron av ett periodiskt gitter, typiska höga temperaturer som krävs för att nå det smälta tillståndet och svårigheten för släckning gör den experimentella bestämningen av en serie uppenbara egenskaper, som densitet, diffusion och viskositet, extremt utmanande. Dessa svårigheter gör alternativa beräkningsmetoder starka och användbara verktyg för att undersöka denna klass av material.
Med tillkomsten av datorkraft och tillgången till superdatorer används för närvarande två stora numeriska atomistiska simuleringstekniker för att studera det dynamiska tillståndet hos ett icke-kristallint atomistiskt system, Monte Carlo1 och molekylär dynamik (MD)1,2. I Monte Carlo-simuleringar provtas det konfigurationsutrymmet slumpmässigt. Monte Carlo-metoder visar linjär skalning i parallellisering om alla provtagningsobservationer är oberoende av varandra. Resultatens kvalitet beror på kvaliteten på slumptalsgeneratorn och provtagningens representativitet. Monte Carlo-metoder visar linjär skalning i parallellisering om provtagningen är oberoende av varandra. I molekylär dynamik (MD) provtas det konfigurationsmässiga utrymmet av tidsberoende atombanor. Från en given konfiguration beräknas atombanorna genom att integrera de newtonska rörelseekvationerna. De interatomära krafterna kan beräknas med hjälp av modellinteratomiska potentialer (i klassisk MD) eller med hjälp av första principer metoder (i ab initio, eller första principer, MD). Kvaliteten på resultaten beror på banans längd och dess förmåga att inte lockas till lokala minima.
Molekyldynamiksimuleringar innehåller en mängd information, allt relaterat till systemets dynamiska beteende. Termodynamiska genomsnittliga egenskaper, som intern energi, temperatur och tryck, är ganska standard att beräkna. De kan extraheras från utdatafilerna för simuleringarna och medelvärdet, medan kvantiteter som är direkt relaterade till atomernas rörelse samt deras ömsesidiga förhållande måste beräknas efter extraktion av atompositioner och hastigheter.
Följaktligen har mycket arbete ägnats åt att visualisera resultaten, och olika paket finns tillgängliga idag, på olika plattformar, öppen källkod eller inte [Ovito3, VMD4, Vesta5, Travis6, etc.]. Alla dessa visualiseringsverktyg hanterar effektivt interatomära avstånd, och som sådan möjliggör de effektiv beräkning av parfördelningsfunktioner och diffusionskoefficienter. Olika grupper som utför storskaliga molekylära dynamiksimuleringar har proprietär programvara för att analysera olika andra egenskaper som härrör från simuleringarna, ibland i shareware eller andra former av begränsad tillgång till samhället, och ibland begränsad i omfattning och användning till vissa specifika paket. Sofistikerade algoritmer för att extrahera information om interatomär bindning, geometriska mönster och termodynamik utvecklas och implementeras i några av dessa paket3,4,5,6,7, etc.
Här föreslår vi UMD-paketet - ett paket med öppen källkod skrivet i Python för att analysera utdata från molekylära dynamiksimuleringar. UMD-paketet möjliggör beräkning av ett brett spektrum av strukturella, dynamiska och termodynamiska egenskaper (figur 1). Paketet är tillgängligt via GitHubs webbplats (https://github.com/rcaracas/UMD_package) och via en dedikerad sida (http://moonimpact.eu/umd-package/) av ERC IMPACT-projektet som ett open access-paket.
För att göra det universellt och lättare att hantera, är vårt tillvägagångssätt att först extrahera all information relaterad till det termodynamiska tillståndet och atombanorna från utdatafilen för den faktiska molekylärdynamikkörningen. Den här informationen lagras i en dedikerad fil, vars format är oberoende av det ursprungliga MD-paketet där simuleringen kördes. Vi namnger dessa filer "umd" filer, som står för Universal Molecular Dynamics. På detta sätt kan vårt UMD-paket enkelt användas av alla ab initio-grupper med vilken programvara som helst, allt med en minimal anpassningsinsats. Det enda kravet för att använda det nuvarande paketet är att skriva lämplig toser från utdata från den specifika MD-programvaran till umd-filformatet, om detta ännu inte finns. För närvarande tillhandahåller vi sådana tolkare för VASP8- och QBox9-paketen.
Bild 1: Flödesschema för UMD-biblioteket.
Fysiska egenskaper är i blått och större Python-skript och deras alternativ är i rött. Klicka här för att se en större version av den här figuren.
Umd-filerna är ASCII-filer; typisk förlängning är "umd.dat" men inte obligatorisk. Alla analyskomponenter kan läsa ASCII-filer i umdformatet, oavsett det faktiska namntillägget. Några av de automatiska skript som är utformade för att utföra snabb storskalig statistik över flera simuleringar letar dock specifikt efter filer med umd.dat-tillägget. Varje fysisk egenskap uttrycks på en rad. Varje rad börjar med ett nyckelord. På detta sätt är formatet mycket anpassningsbart och gör det möjligt att läggas till nya egenskaper i umd-filen, samtidigt som dess läsbarhet bevaras i alla versioner. De första 30 raderna i navelfilen för simuleringen av pyrolit vid 4,6 GPa och 3000 K, som används nedan i diskussionen, visas i figur 2.
Figur 2: Början av umdfilen som beskriver simuleringen av flytande pyrolit vid 4,6 GPa och 3000 K.
Rubriken följs av beskrivningen av varje ögonblicksbild. Varje egenskap är skriven på en rad som innehåller namnet på den fysiska egenskapen, värdet eller enheterna och enheterna, alla åtskilda av blanksteg. Klicka här för att se en större version av den här figuren.
Alla umdfiler innehåller en rubrik som beskriver innehållet i simuleringscellen: antalet atomer, elektroner och atomtyper, samt detaljer för varje atom, såsom dess typ, kemiska symbol, antal valenselektroner och dess massa. En tom rad markerar slutet på huvudet och separerar det från huvuddelen av umdfilen.
Sedan beskrivs varje steg i simuleringen. För det första anges de momentana termodynamiska parametrarna, var och en på en annan linje, som anger i) parameterns namn, som energi, påfrestningar, motsvarande hydrostatiskt tryck, densitet, volym, gitterparametrar etc., ii) dess värde och iii) dess enheter. En tabell som beskriver atomerna kommer härnäst. En rubrikrad ger de olika måtten, som kartesiska positioner, hastigheter, avgifter etc., och deras enheter. Sedan beskrivs varje atom på en rad. Av grupper om tre, motsvarande de tre x- y-z-axlarna , är posterna: de reducerade positionerna, de kartesiska positionerna vikta i simuleringscellen, de kartesiska positionerna (som korrekt tar hänsyn till det faktum att atomer kan korsa flera enhetsceller under en simulering), atomhastigheterna och atomkrafterna. De två sista posterna är skalärer: laddning och magnetiskt ögonblick.
Två stora bibliotek säkerställer att hela paketet fungerar korrekt. Det umd_process.py biblioteket behandlar umdfilerna, som att läsa och skriva ut. det crystallography.py biblioteket behandlar all information relaterad till den faktiska atomstrukturen. Den underliggande filosofin i crystallography.py bibliotek är att behandla gitteret som ett vektoriskt utrymme. Enhetscellparametrarna tillsammans med deras orientering representerar basvektorerna. "Rymden" har en serie skalärattribut (specifik volym, densitet, temperatur och specifikt antal atomer), termodynamiska egenskaper (inre energi, tryck, värmekapacitet etc.) och en serie tensoriella egenskaper (stress och elasticitet). Atomer fyller det här utrymmet. Klassen "Lattice" definierar denna ensemble, tillsammans med olika få korta beräkningar, som specifik volym, densitet, erhållande av det ömsesidiga gitteret från det direkta etc. Klassen "Atomer" definierar atomerna. De kännetecknas av en serie skaläregenskaper (namn, symbol, massa, antal elektroner etc.) och en serie vektoregenskaper (positionen i rymden, antingen i förhållande till den vektoriska grunden som beskrivs i Gitterklassen, eller i förhållande till universella kartesiska koordinater, hastigheter, krafter etc.). Bortsett från dessa två klasser innehåller crystallography.py bibliotek en serie funktioner för att utföra en mängd olika tester och beräkningar, till exempel atomavstånd eller cellförökning. Det periodiska systemet för elementen ingår också som en ordlista.
De olika komponenterna i umdpaketet skriver flera utdatafiler. Som en allmän regel är de alla ASCII-filer, alla deras poster separeras av flikar och de görs så självförklarande som möjligt. Till exempel anger de alltid tydligt den fysiska egenskapen och dess enheter. Umd.dat-filerna följer helt denna regel.
Protocol
1. Analys av molekyldynamiken körs
Paketet är tillgängligt via GitHubs webbplats (https://github.com/rcaracas/UMD_package) och via en särskild sida (http://moonimpact.eu/umd-package/) av ERC IMPACT-projektet som ett open access-paket.
- Extrahera varje specifik uppsättning fysiska egenskaper med ett eller flera dedikerade Python-skript från paketet. Kör alla skript på kommandoraden. De använder alla en serie flaggor, som är så konsekventa som möjligt från ett skript till ett annat. Flaggorna, deras betydelse och standardvärden sammanfattas alla i tabell 1.
| Flagga | Betydelse | Skript som använder det | Standardvärde |
| -h | Kort hjälp | alla | |
| -f | UMD-filnamn | alla | |
| -I | Termiseringssteg som ska kasseras | alla | 0 |
| -I | Indatafil som innehåller de interatomära bindningarna | Artbildning | bonds.input |
| -s | Provtagning av frekvensen | msd, speciation | 1 (varje steg övervägs) |
| -a | Lista över atomer eller anjoner | Artbildning | |
| -c | Lista över katjoner | Artbildning | |
| -Jag | Bindningslängd | Artbildning | 2 |
| -t | Temperatur | vibrationer, reologi | |
| -v | Diskretisering av bredden på provtagningsfönstret i banan för analysen av medelvärdet av kvadratförskjutningen | Msd | 20 |
| -z | Diskretisering av provtagningsfönstrets start för banan för analysen av medelvärdet av kvadratförskjutningen | Msd | 20 |
Tabell 1: De vanligaste flaggorna som används i UMD-paketet och deras vanligaste betydelse.
- Börja med att omvandla utdata från MD-simuleringen som utförs i en första princip kod, som VASP8 eller QBox9, till en UMD-fil.
- Om MD-simuleringarna gjordes i VASP, då vid kommandoradstypen:
VaspParser.py -f-i
där –f-flaggan definierar namnet på VASP OUTCAR-filen och termiseringslängden –i.
OBS: Det första steget, definierat av –i gör det möjligt att kassera de första stegen i simuleringarna, som representerar termualiseringen. I en typisk molekylär-dynamics körning representerar den första delen av beräkningen termualiseringen, det vill säga den tid det tar systemet för alla atomer att beskriva en Gaussisk-liknande fördelning av temperaturen och för hela systemet att uppvisa fluktuationer i temperatur, tryck, energi etc. runt jämviktsvärden. Denna termiseringsdel av simuleringen bör inte beaktas vid analys av vätskans statistiska egenskaper.
- Om MD-simuleringarna gjordes i VASP, då vid kommandoradstypen:
- Förvandla . umd filer till . xyz filer för att underlätta visualisering på olika andra paket, som VMD4 eller Vesta5. Vid kommandoradstypen:
umd2xyz.py -f-i -s
där –f definierar namnet på . umd fil, –i definierar termiseringsperioden som ska kasseras, och – s frekvensen för provtagningen av banan som lagras i . umd fil. Standardvärden är –i 0 –s 1, dvs. med tanke på alla steg i simuleringen, utan att några kasseras. - Återför umdfilen till POSCAR-filer av VASP-typ med det umd2poscar.py-skriptet. ögonblicksbilder av simuleringarna kan väljas med en fördefinierad frekvens. Vid kommandoradstypen:
umd2poscar.py -f-i -l -s
där –l representerar det sista steget som ska omvandlas till POSCAR-fil. Standardvärden är -i 0 -l 10000000 -s 1. Detta värde på –l är tillräckligt stort för att täcka en typisk hel bana.
2. Utför den strukturella analysen
- Kör skriptet gofrs_umd.py för att beräkna pdf-funktionen (pair distribution function) gʙ(r) för alla par av atomtyperna A och B (bild 3). Utdata skrivs i en ASCII-fil, tabbavgränsad, med tillägget gofrs.dat. Vid kommandoradstypen:
gofrs_umd.py -f-s < Sampling_Frequency > -d -i
OBS: Standardvärdena är Sampling_Frequency (frekvensen för provtagning av banan) = 1 steg; DiskretizationInterval (för plottning av g(r)) = 0,01 Å; InitialStep (antal steg i början av banan som ignoreras) = 0. Radiell PDF, gᴀ(r) är det genomsnittliga antalet atomer av typ B på ett avstånd d_ᴀ inom ett sfäriskt skal av radie r och tjocklek dr centrerat på atomerna av typ A (figur 3):
med ρ atomtätheten, NA och NB antalet atomer av typ A och B, och δ(r−rᴀʙ) deltafunktionen som är lika med 1 om atomerna A och B ligger på ett avstånd mellan r och r+dr. Abscissa av det första maximala av g(r) ger den högsta sannolikhetsbindningslängden mellan atomerna av typ A och B, vilket är närmast ett genomsnittligt bindningsavstånd som vi kan bestämma. Det första minimibeloppet begränsar omfattningen av den första samordningssfären. Därför ger integralen över PDF-filen upp till det första minimum det genomsnittliga samordningsnumret. Summan av Fouriertransformerna av g'r(r) för alla par av atomtyperna A och B ger vätskans diffraktionsmönster, som erhållits experimentellt med en diffractometer. Men i verkligheten, som ofta saknas de höga ordningssamordningssfärerna från g(r), kan diffraktionsmönstret inte erhållas i sin helhet.
Bild 3: Bestämning av parfördelningsfunktionen.
a) För varje atom av en art (t.ex. röd) räknas alla atomer av den samordnande arten (t.ex. grå och/eller röd) som en funktion av avstånd. b) Den resulterande avståndsfördelningsgrafen för varje ögonblicksbild, som i detta skede endast är en samling deltafunktioner, beräknas sedan över alla atomer och alla ögonblicksbilder och viktas av den ideala gasfördelningen för att generera (c) den parfördelningsfunktion som är kontinuerlig. Det första minimumet av g(r) är radien för den första samordningssfären, som senare används i speciationsanalysen. Klicka här för att se en större version av den här figuren.
- Extrahera de genomsnittliga interatomiska bindningsavstånden som radierna i de första samordningssfärerna. För detta identifierar du positionen för det första maximala av funktionerna gʙ(r): rita gofrs.dat-filen i ett kalkylbladsprogram och sök efter maxima och minima för varje atompar.
- Identifiera radien för den första samordningssfären, som det första minimumet av PDF, gᴀᷱ(r), med hjälp av kalkylbladsprogram. Detta är grunden för hela den strukturella analysen av vätskan; PDF-filen ger den genomsnittliga bindningsstatusen för atomerna i vätskan.
- Extrahera avstånden från den första minima, dvs abscissa, och skriv dem i en separat fil, kallad till exempel bonds.input. Alternativt kan du köra ett av analyze_gofr skripten i UMD-paketet för att identifiera maxima- och minima-funktionerna för funktionerna gʙ(r). Vid kommandoradstypen:
analyze_gofr_semi_automatic.py - Klicka på positionen för den maximala och minsta av funktionen g(r) som visas i diagrammet som öppnas av programmet. Skriptet skannar automatiskt den aktuella mappen, identifierar alla gofrs.dat-filer och utför analysen för var och en av dem. Klicka igen på det maximala och lägsta i fönstret varje gång skriptet behöver en utbildad första gissning.
- Öppna och titta på den automatiskt genererade filen bonds.input som innehåller de interatomära bindningsavstånden.
3. Utför speciationsanalysen
- Beräkna topologin för bindning mellan atomerna, med hjälp av begreppet anslutning inom grafteori: atomerna är noderna och de interatomära bindningarna är vägarna. Det speciation_umd.py skriptet behöver de interatomära bindningsavstånd som definieras i filen bonds.input .
OBS: Anslutningsmatrisen är konstruerad vid varje tidssteg: två atomer som ligger på ett avstånd som är mindre än radien för motsvarande första samordningssfär anses vara bundna, dvs. Olika atomnätverk byggs genom att atomerna behandlas som noder i ett diagram vars anslutningar definieras av detta geometriska kriterium. Dessa nätverk är atomarter, och deras ensemble definierar atomspekationen i just den vätskan (figur 4).
Figur 4: Identifiering av atomklustren.
Koordinationen polyhedra definieras med hjälp av interatomära avstånd. Alla atomer på ett avstånd som är mindre än en angiven radie anses vara bundna. Här motsvarar tröskelvärdet den första samordningssfären (de ljusröda cirklarna), definierad i figur 1. Polymerisation och därmed kemiska arter erhålls från nätverken av de bundna atomerna. Notera det centrala Red1Grey2-klustret, som är isolerat från de andra atomerna, som bildar en oändlig polymer. Klicka här för att se en större version av den här figuren.
- Kör speciationsskriptet för att erhålla anslutningsmatrisen och få koordinationspolyhedra eller polymerisationen. Vid kommandoradstypen:
speciation_umd.py -f-s -i -l -c -a -m -r
där flaggan -i ger filen med de interatomära bindningsavstånden, som producerades till exempel i föregående steg. Du kan också köra skriptet med en enda längd för alla bindningar som definieras av flaggan -l.
OBS: Flaggan -c anger de centrala atomerna och flaggan -a liganderna. Både centrala atomer och ligands kan vara av olika slag; I detta fall måste de separeras med kommatecken. Flaggan -m anger den minsta tid som en art måste leva för att beaktas i analysen. Som standard är den här minimitiden noll, alla förekomster räknas i den slutliga analysen.- Kör skriptet speciation_umd.py med flaggan –r 0, som tar prov på anslutnings diagrammet på den första nivån för att identifiera samordnings polyhedra. Till exempel kan en central atom, betecknad som en katjon , vara omgiven av en eller flera anjoner (figur 4). Speciation-skriptet identifierar varenda en av koordinationspolyhedra. Det viktade genomsnittet av all samordning polyhedra ger samordningsnumret, identiskt med det som erhålls från integrationen av PDF. Vid kommandoradstypen:
speciation_umd.py -f-i -c -a -r 0
OBS: Genomsnittliga samordningsnummer i vätskor är fraktionerade tal. Denna fraktionalitet kommer från samordningens genomsnittliga kännetecken. Definitionen baserad på speciation ger en mer intuitiv och informativ representation av vätskans struktur, där de relativa proportionerna av de olika arterna, dvs. - Kör skriptet speciation_umd.py med flaggan –r 1, som tar prov på anslutnings diagrammet på alla djupnivåer för att erhålla polymerisationen. Nätverket genom atomgrafen har ett visst djup, eftersom atomer binds längre bort till andra bindningar (t.ex. i sekvenser av alternerande katjoner och anjoner) (figur 4).
- Kör skriptet speciation_umd.py med flaggan –r 0, som tar prov på anslutnings diagrammet på den första nivån för att identifiera samordnings polyhedra. Till exempel kan en central atom, betecknad som en katjon , vara omgiven av en eller flera anjoner (figur 4). Speciation-skriptet identifierar varenda en av koordinationspolyhedra. Det viktade genomsnittet av all samordning polyhedra ger samordningsnumret, identiskt med det som erhålls från integrationen av PDF. Vid kommandoradstypen:
- Öppna de två filerna . folk.dat och . stat.dat i följd. dessa utgör utdata från speciationsskriptet. Varje kluster är skrivet på en rad och anger dess kemiska formel, den tid då den bildades, den tidpunkt då den dog, dess livstid, en matris med listan över atomer som bildar detta kluster. Plotta livslängden för varje atomkluster av alla kemiska arter som finns i simuleringen som finns i .popul.dat-filen (figur 5).
- Plotta populationsanalysen med överflöd av varje art, som finns i . stat.dat fil. Denna analys, både absolut och relativ, motsvarar den faktiska statistiken över samordningspolyhedra för ärendet -r 0. För polymerisation, med -r 1 måste detta behandlas noggrant eftersom viss normalisering över det relativa antalet atomer kan behöva appliceras. Överflöd motsvarar integralen över livstiden. Den. stat.dat fil också listar storleken på varje kluster, dvs. hur många atomer som bildar det.
4. Beräkningskoefficienter för diffusion
- Extrahera atomernas genomsnittliga kvadratförskjutningar (MSD) som en funktion av tiden för att erhålla självdrepartiviteten. Standardformeln för msd är:
där prefaktorerna är renormaliseringar. Med MSD-verktyget finns det olika sätt att analysera de dynamiska aspekterna av vätskorna.
OBS: T är den totala tiden för simuleringen och N α är antalet atomer av typ α. Den första tiden t0 är godtycklig och sträcker sig över den första halvan av simuleringen. Ninit är antalet inledande tider. τ är bredden på det tidsintervall som msd beräknas över. Dess maximala värde är halva simuleringens tidslängd. I typiska msd-implementeringar startar varje fönster i slutet av föregående. Men en glesare provtagning kan påskynda beräkningen av msd, utan att ändra msd-systemets lutning. För detta startar i-th-fönstret vid punkt t0(i), men (i+1)-fönstret börjar vid tidpunkt t0(i) + τ + v, där värdet på v är användardefinierat. På samma sätt ökas fönstrets bredd i diskreta steg som definieras av användaren, som sådan: τ(i) = τ(i-1) + z. Värdena z ("horisontellt steg") och v ("vertikalt steg") är positiva eller noll; Standardvärdet för båda är 20. - Beräkna msd-systemet med hjälp av serien med msd_umd-skript . Deras utdata skrivs ut i en . msd.dat fil, där msd för varje atomtyp, atom eller kluster skrivs ut på en kolumn som en funktion av tiden.
- Beräkna medelvärdet av msd för varje atomtyp. Msd beräknas för varje atom och beräknas sedan i genomsnitt för varje atomtyp. Utdatafilen innehåller en kolumn för varje atomtyp. Vid kommandoradstypen:
msd_umd.py -f-z -v -b - Beräkna msd för varje atom. Msd beräknas för varje atom och beräknas sedan i genomsnitt för varje atomtyp. Utdatafilen innehåller en kolumn för varje atom i simuleringen och sedan en kolumn för varje atomtyp. Denna funktion gör det möjligt att identifiera atomer som sprider sig i två olika miljöer, som vätska och gas, eller två vätskor. Vid kommandoradstypen:
msd_all_umd.py -f-z -v -b - Beräkna msd för den kemiska arten. Använd populationen av kluster som identifieras med speciationsskriptet och skrivs ut i . .dat filen. Msd beräknas för varje enskilt kluster. Utdatafilen innehåller en kolumn för varje kluster. För att undvika att överväga storskaliga polymerer, sätt en gräns för klustrets storlek. dess standard är 20 atomer. Vid kommandoradstypen:
msd_cluster_umd.py -f-p -s -b -c
Standardvärden är: –b 100 –s 1 –c 20.
- Beräkna medelvärdet av msd för varje atomtyp. Msd beräknas för varje atom och beräknas sedan i genomsnitt för varje atomtyp. Utdatafilen innehåller en kolumn för varje atomtyp. Vid kommandoradstypen:
- Rita msd med hjälp av en kalkylbladsbaserad programvara (bild 6). Identifiera lutningsändringen i en loggloggrepresentation av MSD kontra tid. Separera den första delen, vanligtvis kort, som representerar den ballistiska regimen, dvs. bevarandet av atomernas hastighet efter kollisioner. Den andra längre delen representerar diffusiv regim, dvs. spridning av atomernas hastighet efter kollisioner.
- Beräkna diffusionskoefficienterna från msd-systemets lutning som:
där Z är antalet frihetsgrader (Z = 2 för diffusion i plan, Z = 3 för diffusion i rymden), och t är tidssteget.
5. Tidskorrelationsfunktioner
- Beräkna tids korrelationsfunktionerna som ett mått på systemets tröghet med hjälp av den allmänna formeln:
A kan vara en mängd tidsberoende variabler, såsom atompositioner, atomhastigheter, påfrestningar, polarisering etc., var och en ger – via Green-Kubo-relationerna12,13 – olika fysiska egenskaper, ibland efter en ytterligare omvandling. - Analysera atomhastigheterna för att erhålla vibrationsspektrumet hos vätskan och alternativt uttryck för atomiska självdiffusionskoefficienter.
- Kör skriptet vibr_spectrum_umd.py för att beräkna funktionen för automatisk korrelation (VAC) för varje atomtyp och utför dess snabb-Fourier-transformering. Vid kommandoradstypen:
vibr_spectrum_umd.py -f-t
där –t är den temperatur som måste definieras av användaren. Skriptet skriver ut två filer: . vels.scf.dat fil med funktionen VAC för varje atomtyp och . vibr.dat fil med vibrationsspektrumet som förmultnats på varje atomart och det totala värdet. - Öppna och läs vels.scf.dat. Rita VAC-funktionen från filen vels.scf.dat med hjälp av kalkylbladsliknande programvara.
- Behåll den verkliga delen av Fourier VAC. Detta är vad som ger vibrationsspektrumet, som en funktion av frekvens:
där m är atommassorna. - Rita vibrationsspektrumet från vibr.dat-filen med hjälp av kalkylbladsliknande programvara (bild 7). Identifiera det finita värdet vid ω=0 som motsvarar vätskans diffusiva karaktär och spektrumets olika toppar vid ändlig frekvens. Identifiera varje atomtyps deltagande i vibrationsspektrumet.
OBS: Nedbrytningen på atomtyper visar att olika atomer har olika ω=0-bidrag, vilket motsvarar deras diffusionskoefficienter. Spektrumets allmänna form är mycket jämnare med färre funktioner än för en motsvarande fast ämne. - Vid skalet, läs integralen över vibrationsspektrumet, vilket ger diffusionskoefficienterna för varje atomart.
OBS: Termodynamiska egenskaper kan erhållas genom integration från vibrationsspektrumet, men resultaten bör användas med försiktighet på grund av två approximationer: integrationen är giltig inom den kvasiharmoniska approximationen, som inte nödvändigtvis håller vid höga temperaturer; och den gasliknande delen av spektrumet som motsvarar diffusionen måste kasseras. Integrationen bör då endast göras över den gitterliknande delen av spektrumet. Men denna separation kräver vanligtvis flera ytterligare efterbehandlingssteg och beräkningar14, som inte omfattas av det nuvarande UMD-paketet.
- Kör skriptet vibr_spectrum_umd.py för att beräkna funktionen för automatisk korrelation (VAC) för varje atomtyp och utför dess snabb-Fourier-transformering. Vid kommandoradstypen:
- Kör viscosity_umd.py skript för att analysera själva korrelationen av komponenternas stress tensor för att uppskatta smältningens viskositet. Vid kommandoradstypen:
viscosity_umd.py -f-i -s -o -l
OBS: Den här funktionen är utforskande och eventuella resultat måste tas med försiktighet. För det första noggrant kontrollera viskositetens konvergens med avseende på simuleringens längd.- Härled vätskans viskositet från stressens självkorrelation tensor15 som:
Där V respektive T är volymen respektive temperaturen är κB Boltzmann-konstanten och σ ij den ij-diagonala komponenten i stress-tensor, uttryckt i kartesiska koordinater. - Använd en mer adekvat passform för att få en mer robust uppskattning av viskositeten15,16 och undvik bruset från den automatiska korrelationsfunktionen för stress-tensor som kan uppstå från simuleringarnas begränsade storlek och ändliga varaktighet. För stressvändarens automatiska korrelationsfunktion använder du följande funktionella form15,16 som ger goda resultat:
där A, B, τ1, τ2 och ω är passformsparametrar. Efter integrering blir uttrycket för viskositeten:
- Härled vätskans viskositet från stressens självkorrelation tensor15 som:
6. Termodynamiska parametrar som härrör från simuleringarna.
- Kör averages.py för att extrahera medelvärdena och spridningen (som standardavvikelse) för tryck, temperatur, densitet och intern energi från umdfilerna . Vid kommandoradstypen:
averages.py -f-s
med –s 0 som standard. - Beräkna medelvärdets statistiska fel med hjälp av blockeringsmetoder.
OBS: Det finns olika smaker av denna metod. Efter Allens och Tildesley2s arbete är det vanligt att genomsnitt över sekvenser av tidsblock, av allt längre längd, uppskatta standardavvikelsen med avseende på det aritmetiska medelvärdet17. Konvergens kan uppnås i gräns för många och tillräckligt långa blockstorlekar, när provtagningen inte är korrelerad. Även om det faktiska tröskelvärdet för konvergensen vanligtvis måste väljas manuellt.- Använd halveringsmetoden18: från och med det ursprungliga dataprovet, i varje steg κ, halvera antalet prov genom att i genomsnitt halvera antalet prover genom att i genomsnitt över vartannat motsvarande på varandra följande prover från föregående steg κ−1:
- Kör skriptet fullaverages.py för att utföra den fullständiga statistiska analysen, inklusive felet för medelvärdet. Vid kommandoradstypen:
fullaverages.py -s-u
Skriptet automatiseras till den grad att du söker efter alla UMD.dat-filer i den aktuella katalogen och utför analysen för dem alla. Standardvärden är –s 0 –u 0. För -u 0 är utgången minimal, och för -u 1 är utdata i sin helhet, med flera alternativa enheter utskrivna. Det här skriptet kräver grafiskt stöd, eftersom det skapar en grafisk bild för att kontrollera konvergensen för att uppskatta felet på medelvärdet.
- Använd halveringsmetoden18: från och med det ursprungliga dataprovet, i varje steg κ, halvera antalet prov genom att i genomsnitt halvera antalet prover genom att i genomsnitt över vartannat motsvarande på varandra följande prover från föregående steg κ−1:
Representative Results
Pyrolit är en modell av flerkomponentssilikatsmälta (0,5Na2O 2CaO 1.5Al2O3 4FeO 30MgO 24SiO2) som bäst approximerar sammansättningen av bulksilikatjorden - det geokemiska genomsnittet eller hela vår planet, förutom dess järnbaserade kärna19. Den tidiga jorden dominerades av en serie storskaliga smältande händelser20, den sista kan ha uppslukat hela planeten, efter dess kondensation för protolunarskivan21. Pyrolit representerar den bästa approximanten till sammansättningen av sådana planetariska magmahav. Följaktligen studerade vi i stor utsträckning de fysiska egenskaperna hos pyrolitsmältning i temperaturområdet 3 000\u20125 000 K och 0\u2012150 GPa-tryckområdet från ab initio molekylär-dynamiska simuleringar i VASP-implementeringen. Dessa termodynamiska förhållanden kännetecknar helt jordens mest extrema magmahavsförhållanden. Vår studie är ett utmärkt exempel på en framgångsrik användning av UMD-paketet för hela den djupgående analysen av smältorna22. Vi beräknade fördelningen och de genomsnittliga bindningslängderna, vi spårade förändringarna i katjon-syre samordning och jämförde våra resultat med tidigare experimentella och beräkningsstudier på amorfa silikater av olika kompositioner. Vår djupgående analys hjälpte till att sönderdela standardkoordineringsnummer i sina grundläggande beståndsdelar, beskriva förekomsten av exotisk samordning polyhedra i smältan och extrahera livstider för all samordning polyhedra. Den beskrev också vikten av provtagning i simuleringar både när det gäller banans längd och antalet atomer som finns i det system som modelleras. När det gäller efterbehandlingen är UMD-analysen oberoende av dessa faktorer, men de bör beaktas vid tolkningen av resultaten från UMD-paketet. Här visar vi några exempel på hur UMD-paketet kan användas för att extrahera flera karakteristiska egenskaper hos smältningarna, med en applikation för att smälta pyrolit.
Si-O-parfördelningsfunktionen som erhålls från gofrs_umd.py-skriptet visar att radien för den första samordningssfären, som är det första minimumet av g(r)-funktionen, ligger cirka 2,5 angstroms vid T = 3000 K och P = 4,6 GPa. Maximalt g(r) ligger på 1.635 Å - detta är den bästa approximationen till böjlängden. Den långa svansen beror på temperaturen. Med denna gräns som Si-O-bindningsavstånd visar speciationsanalysen att SiO4-enheter, som kan pågå i upp till några picosekunder, dominerar smältan (figur 5). Det finns en viktig del av smältan som visar partiell polymerisation, vilket återspeglas av närvaron av dimers som Si2O7 och trimers som Si3Ox-enheter. Deras motsvarande livstid är i picosekundens ordning. Polymerer med högre ordning har alla betydligt kortare livslängd.
Figur 5: Livslängden för de kemiska arterna Si-O.
Speciationen identifieras i en flerkomponentsmälta vid 4,6 GPa och 3000 K. Etiketterna markerar siO3-, SiO4- och SiO5-monomererna och de olika SixOy-polymererna. Klicka här för att se en större version av den här figuren.
De olika värdena för de vertikala och horisontella stegen, som definieras av flaggorna –z och –v ovan, ger olika provtagningar av msd-systemet (figur 6). Även stora värden av z och v är tillräckliga för att definiera sluttningarna och därmed diffusionskoefficienterna för de olika atomerna. Vinsten i tid för efterbehandlingen är anmärkningsvärd när man går till stora värden av z och v. Msd erbjuder ett mycket starkt valideringskriterium för simuleringarnas kvalitet. Om diffusionsdelen av msd inte är tillräckligt lång är det ett tecken på att simuleringen är för kort och inte når vätsketillståndet i statistisk mening. Minimikravet för den diffusiva delen av msd beror i hög grad på systemet. Man kan kräva att alla atomer ändrar sin plats minst en gång i smältans struktur för att den ska betraktas som en fluid10. Ett utmärkt exempel med tillämpningar inom planetvetenskap är komplexa silikat smälter vid höga tryck nära eller till och med under deras liquiduslinje11. Si-atomerna, de stora nätverksbildande katjonerna, byter plats efter mer än två dussin picosekunder. Simuleringar som är kortare än detta tröskelvärde skulle vara betydligt underprovtagning av det möjliga konfigurationsutrymmet. Men eftersom de koordinerande anjonerna, nämligen O-atomerna, rör sig snabbare än de centrala Si-atomerna, kan de kompensera för en del av Si långsamma rörlighet. Som sådan skulle hela systemet faktiskt kunna omfatta en bättre provtagning av det konfigurationsutrymme som endast antas från Si-förskjutningarna.
Figur 6: Medel kvadratförskjutningar (MSD).
MSD illustreras för några atomtyper av en multikomponents silikatsmälta. Provtagningen med olika horisontella och vertikala steg, z och v, ger konsekventa resultat. Fasta cirklar: -z 50 –v 50. Öppna cirklar: -z 250 –v 500. Klicka här för att se en större version av den här figuren.
Slutligen ger de atomiska VAC-funktionerna smältens vibrationsspektrum. Figur 7 visar spektrumet vid samma tryck- och temperaturförhållanden som ovan. Vi representerar bidragen från Mg-, Si- och O-atomer, liksom det totala värdet. Vid noll frekvens finns ett begränsat värde av spektrumet, vilket motsvarar smältens diffusionstecken. Extraktion av de termodynamiska egenskaperna från vibrationsspektrumet måste ta bort denna gasliknande diffusiva karaktär från noll men också för att korrekt ta hänsyn till dess förfall vid högre frekvenser.
Bild 7: Vibrationsspektrumet hos pyrolitsmältning.
Den verkliga delen av Fourier-transformeringen av den atomiska hastighets-självkorrelationsfunktionen ger vibrationsspektrumet. Här beräknas spektrumet för en multikomponents silikatsmältning. Vätskor har en icke-noll gasliknande diffusiv karaktär med noll frekvens. Klicka här för att se en större version av den här figuren.
Discussion
UMD-paketet har utformats för att fungera bättre med ab initio-simuleringar, där antalet ögonblicksbilder vanligtvis är begränsat till tiotusentals till hundratusentals ögonblicksbilder, med några hundra atomer per enhetscell. Större simuleringar är också dragbara förutsatt att den dator där efterbearbetningskörningarna har tillräckligt med aktiva minnesresurser. Koden utmärker sig genom de olika egenskaper som den kan beräkna och genom sin licens med öppen källkod.
Umd.dat-filerna är lämpliga för de ensembler som bevarar antalet partiklar oförändrade under simuleringen. UMD-paketet kan läsa filer som härrör från beräkningar där simuleringsrutans form och volym varierar. Dessa täcker de vanligaste beräkningarna, som NVT och NPT, där antalet partiklar, N, temperatur T, volym, V och/ eller tryck, P, hålls konstant.
För tiden börjar parfördelningsfunktionen samt alla skript som behöver uppskatta de interatomära avstånden, som speciationsskripten, fungerar endast för ortogonala enhetsceller, vilket betyder för kubiska, tetragonala och ortopediska celler, där vinklarna mellan axlarna är 90°.
De viktigaste utvecklingslinjerna för version 2.0 är avlägsnande av ortogonalitetsbegränsningen för avstånd och lägga till fler funktioner för speciationsskripten: att analysera enskilda kemiska bindningar, att analysera de interatomära vinklarna och att genomföra den andra samordningssfären. Med hjälp av externt samarbete arbetar vi med att portera koden till en GPU för snabbare analys i större system.
Disclosures
Författarna har inget att avslöja.
Acknowledgments
Detta arbete stöddes av Europeiska forskningsrådet (ERC) inom ramen för Europeiska unionens forsknings- och innovationsprogram Horisont 2020 (bidragsavtalsnummer 681818 IMPACT till RC), av direktoratet för extrem fysik och kemi vid observatoriet för djup kol och av Norges forskningråd genom sitt finansieringsprogram Centers of Excellence, projektnummer 223272. Vi bekräftar åtkomst till GENCI-superdatorerna genom stl2816-serien av eDARI-datorbidrag, till Irene AMD-superdatorn genom PRACE RA4947-projektet och Fram superdator genom UNINETT Sigma2 NN9697K. FS stöddes av ett Marie Skłodowska-Curie-projekt (bidragsavtal 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. 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. , Oxford University Press. (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).
