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

Technisches Resümee: Ermöglichung von Ab-initio-Molekulardynamik unter Bias via der CP2K+SMEAGOL-Schnittstelle

Problemstellung
Obwohl die Dichtefunktionaltheorie (DFT) in Kombination mit Nichtgleichgewichts-Green-Funktionen (NEGF) ein Standardansatz zur Modellierung des Quantentransports unter externen Bias-Potenzialen ist, sind konventionelle Methodiken auf statische Systeme beschränkt. Bestehende Implementierungen operieren typischerweise innerhalb eines kanonischen Rahmens und berechnen Transporteigenschaften bei festen Atompositionen. Dies erzeugt eine signifikante Lücke zwischen dem theoretischen Verständnis und den makroskopischen Ergebnissen von Experimenten, insbesondere für elektrochemische (EC) Systeme, deren Transformationen intrinsisch grand-kanonisch sind. Um die natürliche Umgebung von EC-Transformationen genau zu beschreiben, erfordern Simulationen Molekulardynamik (MD) unter Bias sowie eine explizite Open-Boundary-Beschreibung von Elektronen, die in die Rechenzelle eintreten und sie verlassen können. Darüber hinaus fehlt es bestehenden DFT-NEGF-Implementierungen oft an der Fähigkeit, die für die Geometrieoptimierung und Dynamik notwendigen strominduzierten Kräfte zu berechnen, oder sie sind auf spezifische DFT-Codes (z. B. SIESTA) beschränkt, denen Funktionen moderner Pakete, die für kondensierte Phasen optimiert sind, fehlen.

Methodik
Die Autoren präsentieren eine neue Schnittstelle zwischen dem populären Mixed-Gaussian/Plane-Wave-Elektronikstrukturpaket CP2K und dem NEGF-Code SMEAGOL. Diese Schnittstelle, bezeichnet als CP2K+SMEAGOL, ist als eigenständige Bibliothek konzipiert, die mit jeder DFT-Software verknüpft werden kann, wobei der ursprüngliche SMEAGOL-Code so weit wie möglich wiederverwendet wird.

Wesentliche technische Implementierungen umfassen:

• Strominduzierte Kräfte: Die Schnittstelle implementiert die Berechnung nicht-konservativer Kräfte, die unter Bias auf Ionen wirken. Diese Kräfte werden aus der zeitlichen Ableitung der Imp momenta der Atome abgeleitet, unter Verwendung der Gradienten des vollständigen Hamiltonians bezüglich der Atompositionen. Entscheidend ist, dass die Implementierung den Kraftterm enthält, der aus der Ableitung der Dichtematrix ($\Omega_{ij}$) resultiert, welcher die Komponente des „Elektronenwinds“ (Impulsübertrag von Elektronen auf Ionen) berücksichtigt. Dieser Term wird in anderen Implementierungen oft vernachlässigt, ist jedoch für die qualitative Korrektheit essenziell.
• Konstruktion der Dichtematrix: Um gebundene Zustände innerhalb des Bias-Fensters korrekt zu besetzen, aktualisiert die Schnittstelle die Berechnung der Dichtematrix. Anstatt eine einzelne Kontur-Evaluierung oder einen einfachen Durchschnitt von Links-/Rechts-Beiträgen zu verwenden, nutzt sie ein gewichtetes Doppel-Kontur-Schema (folgend nach Brandbyge et al.), das Integrale von beiden Leads mit Gewichten zusammenfasst, die durch die Nichtgleichgewichts-Beiträge bestimmt werden.
• Simulationsaufbau: Die Methode ermöglicht Born-Oppenheimer-Molekulardynamik unter Bias. Das System wird mit einer zentralen „erweiterten Molekül“-Region modelliert, die an semi-infinite Elektroden angeschlossen ist. Die Selbstenergien der Leads werden als statisch und unabhängig von der Ladendichte der zentralen Region behandelt, eine Annahme, die gültig ist, wenn die zentrale Region groß genug ist, um Ladungsänderungen abzuschirmen.

Wesentliche Ergebnisse und Validierung
Die Autoren validierten die CP2K+SMEAGOL-Schnittstelle gegen Referenzberechnungen, die mit SIESTA+SMEAGOL an mehreren Systemen durchgeführt wurden:

1. Null-Bias-Kräfte (Unendlicher Au-Draht): Die Schnittstelle reproduziert erfolgreich die Atomkräfte bei Null-Bias und stimmt mit CP2K-eigenen Berechnungen überein. Sie zeigt zudem auf, dass das Vernachlässigen des Energie-Dichtematrix-Terms ($\Omega$) zu qualitativ falschen Kräften führt.
2. Parallelplattenkondensator: Die Schnittstelle reproduziert den erwarteten Potenzialabfall über einem Kondensator und zeigt ein konstantes Potenzial in den Leads sowie einen linearen Abfall im Vakuum. Sie bestätigt, dass das gewichtete Doppel-Kontur-Schema notwendig ist, um symmetrische Ladungsdichten und korrekte Hartree-Potenziale bei hohem Bias zu erzeugen, während Single-Kontur-Evaluierungen asymmetrische, inkorrekte Ergebnisse liefern.
3. Au-H$_2$-Au-Verbindung Geometrieoptimierung: Unter angelegtem Bias (bis zu 1,5 V) sagt die Schnittstelle korrekt die Elongation der H-H-Bindung und atomare Verschiebungen voraus, konsistent mit früheren SIESTA+SMEAGOL-Ergebnissen. Diese Änderungen werden dem elektrischen Strom und nicht dem elektrischen Feld zugeschrieben.
4. Molekulardynamik eines solvatisierten Au-Drahtes: Die Autoren führten die ersten DFT-NEGF-Molekulardynamik-Simulationen an einem großskaligen System der kondensierten Phase (ein monoatomarer Au-Draht, solvatisiert durch 166 Wassermoleküle) durch.
   • Stabilität: Die Dynamik war stabil mit minimalem langfristigem Energie-Drift, trotz der Open-System-Natur.
   • Elektronenmigration: In einer 200 fs Trajektorie bei 1 V Bias beobachteten die Simulationen Effekte der Elektronenmigration, bei denen sich das lokale Minimum der Elektronendichte vom zentralen Drahtatom zum rechtsseitigen Atom verschob, was konsistent mit der Abschirmung der positiven Ladung der rechten Elektrode ist.
   • Transmission: Die Transmissionsspektren des solvatisierten Drahtes zeigten eine qualitative Ähnlichkeit zu einem idealen unendlichen Draht, wobei Unterschiede auf zusätzliche Tunnelkanäle zurückzuführen sind, die durch Wassermoleküle geöffnet werden.

Performance
Der primäre Rechenengpass ist die Evaluierung der Elektronendichte als Integral der Green-Funktion, was 96 Energiepunkt-Berechnungen (64 reelle Achse, 32 komplexe Achse) pro SCF-Schritt erfordert, anstatt einer einzelnen Hamiltonian-Diagonalisierung.

• Parallelisierung: Der Code nutzt MPI für Realraum-Berechnungen und k-Punkte sowie OpenMP für DO-Schleifen und LAPACK-Routinen. Die Performance skaliert gut bis zu 32 OpenMP-Threads.
• Kosten: Für ein solvatisiertes Au-Draht-System war die CP2K+SMEAGOL-Berechnung etwa 3,2-mal langsamer pro SCF-Schritt als eine Standard-CP2K-Berechnung, obwohl 8-mal mehr Kerne (2048 vs. 256) verwendet wurden. Die Autoren merken an, dass dies eine Steigerung der Kosten um eine Größenordnung darstellt, was lange MD-Trajektorien ohne weitere Beschleunigung (z. B. via ScaLAPACK oder Machine-Learning-Kraftfeldern) herausfordernd macht.

Bedeutung und Ansprüche
Das Paper behauptet, dass CP2K+SMEAGOL die am reichsten ausgestattete und vollständigste DFT-NEGF-Implementierung für CP2K ist, die die native CP2K-Implementierung, welche auf $\Gamma$-Punkt-Transport ohne Kräfte beschränkt ist, übertrifft.

Die primäre Bedeutung dieser Arbeit liegt in der Demonstration, dass DFT-NEGF verwendet werden kann, um Molekulardynamik-Simulationen unter Bias für großskalige Systeme der kondensierten Phase unter realistischen Betriebsbedingungen durchzuführen. Nach Kenntnis der Autoren stellt dies die erste solche Berechnung dar. Die Schnittstelle ermöglicht die Untersuchung elektrokatalytischer Reaktionen und elektrochemischer Systeme, bei denen das Zusammenspiel von angelegtem Potenzial, Stromfluss und atomischer Dynamik entscheidend ist, und schließt damit die Lücke zwischen mikroskopischen theoretischen Modellen und makroskopischen experimentellen Bedingungen.