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

Technical Summary: Enabling Ab-Initio Molecular Dynamics under Bias via the CP2K+SMEAGOL Interface

Problem Statement
While Density Functional Theory (DFT) combined with Non-Equilibrium Green's Functions (NEGF) is a standard approach for modeling quantum transport under external bias potentials, conventional methodologies are limited to static systems. Existing implementations typically operate within a canonical framework, calculating transport properties at fixed atomic positions. This creates a significant gap between theoretical understanding and the macroscopic outcomes of experiments, particularly for electrochemical (EC) systems where transformations are intrinsically grand canonical. To accurately describe the natural environment of EC transformations, simulations require molecular dynamics (MD) under bias and an explicit open-boundary description of electrons that can enter and leave the computational cell. Furthermore, existing DFT-NEGF implementations often lack the ability to calculate current-induced forces necessary for geometry optimization and dynamics, or are restricted to specific DFT codes (e.g., SIESTA) that lack features present in state-of-the-art packages optimized for condensed phase simulations.

Methodology
The authors present a new interface between the popular mixed Gaussian/plane-wave electronic structure package CP2K and the NEGF code SMEAGOL. This interface, termed CP2K+SMEAGOL, is designed to be a standalone library linkable with any DFT software, reusing the original SMEAGOL code as much as possible.

Key technical implementations include:

• Current-Induced Forces: The interface implements the calculation of non-conservative forces acting on ions under bias. These forces are derived from the time derivative of atomic momenta, utilizing the gradients of the full Hamiltonian with respect to atomic positions. Crucially, the implementation includes the force term arising from the derivative of the density matrix ($\Omega_{ij}$), which accounts for the "electron wind" component (momentum transfer from electrons to ions). This term is often neglected in other implementations but is shown to be essential for qualitative correctness.
• Density Matrix Construction: To correctly occupy bound states within the bias window, the interface updates the density matrix calculation. Instead of using a single contour evaluation or a simple average of left/right contributions, it employs a weighted double-contour scheme (following Brandbyge et al.) that sums integrals from both leads with weights determined by the non-equilibrium contributions.
• Simulation Setup: The method allows for Born-Oppenheimer molecular dynamics under bias. The system is modeled with a central "extended molecule" region attached to semi-infinite electrodes. The self-energies of the leads are treated as static and independent of the central region's charge density, an assumption valid when the central region is sufficiently large to screen charge changes.

Key Results and Validation
The authors validated the CP2K+SMEAGOL interface against reference calculations performed with SIESTA+SMEAGOL across several systems:

1. Zero-Bias Forces (Infinite Au Wire): The interface successfully reproduces atomic forces at zero bias, matching CP2K-only calculations. It also demonstrates that neglecting the energy density matrix term ($\Omega$) leads to qualitatively incorrect forces.
2. Parallel-Plate Capacitor: The interface reproduces the expected potential drop across a capacitor, showing a constant potential in the leads and a linear drop in the vacuum. It confirms that the weighted double-contour scheme is necessary to produce symmetric charge densities and correct Hartree potentials at high bias, whereas single-contour evaluations yield asymmetric, incorrect results.
3. Au-H$_2$-Au Junction Geometry Optimization: Under applied bias (up to 1.5 V), the interface correctly predicts the elongation of the H-H bond and atomic displacements, consistent with previous SIESTA+SMEAGOL results. The changes are attributed to electric current rather than the electric field.
4. Molecular Dynamics of a Solvated Au Wire: The authors performed the first DFT-NEGF molecular dynamics simulations on a large-scale condensed phase system (a monoatomic Au wire solvated by 166 water molecules).
   • Stability: The dynamics were stable with minimal long-term energy drift, despite the open-system nature.
   • Electron Migration: In a 200 fs trajectory at 1 V bias, the simulation observed electron migration effects, where the local electron density minimum shifted from the central wire atom to the rightmost atom, consistent with screening the positive charge of the right electrode.
   • Transmission: Transmission spectra for the solvated wire showed qualitative similarity to an ideal infinite wire, with differences attributed to additional tunneling channels opened by water molecules.

Performance
The primary computational bottleneck is the evaluation of the electron density as an integral of the Green's function, requiring 96 energy point computations (64 real axis, 32 complex axis) per SCF step, replacing a single Hamiltonian diagonalization.

• Parallelization: The code utilizes MPI for real-space computations and k-points, and OpenMP for DO loops and LAPACK routines. Performance scales well up to 32 OpenMP threads.
• Cost: For a solvated Au wire system, the CP2K+SMEAGOL calculation was approximately 3.2 times slower per SCF step than a standard CP2K calculation, despite using 8 times more cores (2048 vs. 256). The authors note this represents an order-of-magnitude increase in cost, making long MD trajectories challenging without further acceleration (e.g., via ScaLAPACK or machine learning force fields).

Significance and Claims
The paper claims that CP2K+SMEAGOL is the most feature-rich and complete DFT-NEGF implementation available for CP2K to date, surpassing the native CP2K implementation which is limited to $\Gamma$-point transport without forces.

The primary significance of this work is the demonstration that DFT-NEGF can be used to perform molecular dynamics simulations under bias for large-scale condensed phase systems under realistic operating conditions. To the authors' knowledge, this represents the first such calculations performed. The interface enables the study of electro-catalytic reactions and electrochemical systems where the interplay between applied potential, current flow, and atomic dynamics is critical, bridging the gap between microscopic theoretical models and macroscopic experimental conditions.