Enabling Ab-Initio Molecular Dynamics under Bias: The CP2K+SMEAGOL Interface for Integrating Density Functional Theory and Non-Equilibrium Green Functions

Technische Samenvatting: Het mogelijk maken van Ab-Initio Moleculaire Dynamica onder Bias via de CP2K+SMEAGOL Interface

Probleemstelling
Hoewel Density Functional Theory (DFT) gecombineerd met Non-Equilibrium Green's Functions (NEGF) een standaardmethode is voor het modelleren van kwantumtransport onder externe biaspotentialen, zijn conventionele methodologieën beperkt tot statische systemen. Bestaande implementaties opereren doorgaans binnen een canoniek kader, waarbij transporteigenschappen worden berekend bij vaste atomaire posities. Dit creëert een aanzienlijke kloof tussen het theoretisch begrip en de macroscopische uitkomsten van experimenten, met name voor elektrochemische (EC) systemen waar transformaties inherent grand-canonisch zijn. Om de natuurlijke omgeving van EC-transformaties accuraat te beschrijven, vereisen simulaties moleculaire dynamica (MD) onder bias en een expliciete open-boundary beschrijving van elektronen die de computationele cel kunnen binnenkomen en verlaten. Bovendien missen bestaande DFT-NEGF implementaties vaak het vermogen om stroom-geïnduceerde krachten te berekenen die noodzakelijk zijn voor geometrie-optimalisatie en dynamica, of zijn ze beperkt tot specifieke DFT-codes (bijv. SIESTA) die niet beschikken over de functies die aanwezig zijn in state-of-the-art pakketten geoptimaliseerd voor gecondenseerde fasen.

Methodologie
De auteurs presenteren een nieuwe interface tussen de populaire gemengde Gaussian/plane-wave elektronische structuurpackage CP2K en de NEGF-code SMEAGOL. Deze interface, getiteld CP2K+SMEAGOL, is ontworpen als een standalone library die gelinkt kan worden met elke DFT-software, waarbij de originele SMEAGOL-code zo veel mogelijk wordt hergebruikt.

Belangrijke technische implementaties zijn onder meer:

• Stroom-geïnduceerde krachten: De interface implementeert de berekening van niet-conservatieve krachten die op ionen inwerken onder bias. Deze krachten worden afgeleid van de tijd afgeleide van de atomaire impulsen, gebruikmakend van de gradiënten van de volledige Hamiltoniaan met betrekking tot de atomaire posities. Cruciaal is dat de implementatie de krachtterm bevat die voortvloeit uit de afgeleide van de dichtheidsmatrix ($\Omega_{ij}$), wat de "elektronenwind"-component (momentumoverdracht van elektronen naar ionen) account houdt. Deze term wordt in andere implementaties vaak genegeerd, maar wordt getoond essentieel te zijn voor kwalitatieve correctheid.
• Constructie van de Dichtheidsmatrix: Om gebonden toestanden binnen de bias-window correct te bezetten, werkt de interface de berekening van de dichtheidsmatrix bij. In plaats van een enkele contour-evaluatie of een eenvoudige gemiddelde van links/rechts bijdragen te gebruiken, hanteert het een gewogen dubbele-contour methode (volgens Brandbyge et al.) die integralen van beide leads optelt met gewichten bepaald door de niet-evenwicht bijdragen.
• Simulatie Setup: De methode maakt Born-Oppenheimer moleculaire dynamica onder bias mogelijk. Het systeem wordt gemodelleerd met een centrale "uitgebreide molecuul"-regio die is gekoppeld aan semi-oneindige elektroden. De zelf-energieën van de leads worden behandeld als statisch en onafhankelijk van de ladingdichtheid van de centrale regio, een aanname die geldig is wanneer de centrale regio voldoende groot is om ladingveranderingen af te schermen.

Belangrijkste Resultaten en Validatie
De auteurs hebben de CP2K+SMEAGOL interface gevalideerd tegen referentieberekeningen uitgevoerd met SIESTA+SMEAGOL over verschillende systemen:

1. Nul-bias krachten (Oneindige Au-draad): De interface reproduceert succesvol de atomaire krachten bij nul bias, overeenkomend met CP2K-berekeningen. Het demonstreert ook dat het negeren van de energie-dichtheidsmatrix term ($\Omega$) leidt tot kwalitatief incorrecte krachten.
2. Parallelle-plaat Condensator: De interface reproduceert de verwachte potentiaalval over een condensator, waarbij een constante potentiaal in de leads wordt getoond en een lineaire daling in het vacuüm. Het bevestigt dat de gewogen dubbele-contour methode noodzakelijk is om symmetrische ladingdichtheden en correcte Hartree-potentialen bij hoge bias te produceren, terwijl enkelvoudige contour-evaluaties asymmetrische, incorrecte resultaten opleveren.
3. Au-H$_2$-Au Junction Geometrie-optimalisatie: Onder toegepaste bias (tot 1,5 V) voorspelt de interface correct de verlenging van de H-H binding en atomaire verplaatsingen, consistent met eerdere SIESTA+SMEAGOL resultaten. De veranderingen worden toegeschreven aan elektrische stroom in plaats van het elektrische veld.
4. Moleculaire Dynamica van een Gesolveerde Au-draad: De auteurs hebben de eerste DFT-NEGF moleculaire dynamica simulaties op een grootschalig gecondenseerd fasesysteem uitgevoerd (een monoatomair Au-draad gesolveerd door 166 watermoleculen).
   • Stabiliteit: De dynamica was stabiel met minimale langetermijn energie-drift, ondanks de open-systeem natuur.
   • Elektronenmigratie: In een traject van 200 fs bij 1 V bias observeerde de simulatie effecten van elektronenmigratie, waarbij het lokale elektronendichtheid-minimum verschoof van het centrale draadatoom naar het rechtermost atoom, consistent met het afschermen van de positieve lading van de rechter elektrode.
   • Transmissie: Transmissiespectra voor de gesolveerde draad vertoonden een kwalitatieve gelijkenis met een ideale oneindige draad, waarbij verschillen werden toegeschreven aan aanvullende tunnelkanalen geopend door watermoleculen.

Prestaties
De primaire computationele bottleneck is de evaluatie van de elektronendichtheid als een integraal van de Green-functie, wat 96 energiepunt-berekeningen vereist (64 reële as, 32 complexe as) per SCF-stap, ter vervanging van een enkele Hamiltoniaan diagonalisatie.

• Parallellisatie: De code maakt gebruik van MPI voor real-space berekeningen en k-punten, en OpenMP voor DO-loops en LAPACK-routines. De prestaties schalen goed tot 3ate 32 OpenMP threads.
• Kosten: Voor een gesolveerd Au-draad systeem was de CP2K+SMEAGOL berekening ongeveer 3,2 keer trager per SCF-stap dan een standaard CP2K berekening, ondanks het gebruik van 8 keer meer cores (2048 vs. 256). De auteurs merken op dat dit een orde van grootte toename in kosten vertegenwoordigt, wat lange MD-trajecten uitdagend maakt zonder verdere versnelling (bijv. via ScaLAPACK of machine learning krachtvelden).

Betekenis en Claims
Het artikel claimt dat CP2K+SMEAGOL de meest feature-rijke en complete DFT-NEGF implementatie is die tot nu toe beschikbaar is voor CP2K, waarmee het de native CP2K implementatie overtreft die beperkt is tot $\Gamma$-punt transport zonder krachten.

De primaire betekenis van dit werk is de demonstratie dat DFT-NEGF gebruikt kan worden voor het uitvoeren van moleculaire dynamica simulaties onder bias voor grootschalige gecondenseerde fasesystemen onder realistische operationele condities. Naar weten van de auteurs betreft dit de eerste dergelijke berekeningen die zijn uitgevoerd. De interface maakt het mogelijk om elektro-katalytische reacties en elektrochemische systemen te bestuderen waarbij de interactie tussen toegepaste potentiaal, stroomstroom en atomaire dynamica cruciaal is, waardoor de kloof tussen microscopische theoretische modellen en macroscopische experimentele condities wordt overbrugd.