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

Density functional theory (DFT) combined with non-equilibrium Green’s functions (NEGF) is a powerful approach to model quantum transport under external bias potentials at reasonable computational cost. In this work, we present a new interface between the popular mixed Gaussian/plane waves electronic structure package, CP2K, and the NEGF, code SMEAGOL, the most feature-rich implementation of DFT-NEGF available for CP2K to date. The CP2K+SMEAGOL interface includes the implementation of current induced forces. We verify this implementation for a variety of systems: an infinite 1D Au wire, a parallel-plate capacitor, and a Au–H2–Au junction. We find good agreement with SMEAGOL calculations performed with SIESTA for the same systems and with the example of a solvated Au wire demonstrating for the first time that DFT-NEGF can be used to perform molecular dynamics simulations under bias of large-scale condensed phase systems under realistic operating conditions.


Introduction
Improving the performance of electro-catalytic reactions such as in batteries, solar cells and capacitors is essential to address growing energy demands and costs.While much progress has been made, there remains a significant gap between the theoretical understanding of microscopic phenomena and the macroscopic outcomes of experiments. 1 Simulation methods such density functional theory (DFT) can provide valuable information about the structure and dynamical properties of electrochemical (EC) systems, 2,3 however conventional methodologies do not allow for the study of systems under realistic operating conditions such as under applied potential and current flow.Indeed, these methodologies are implemented within a canonical framework, but EC transformation is intrinsically grand canonical.To describe the natural environment of EC transformation, we need to perform molecular dynamics under bias and have an explicit open-boundary description of the electrons, which must be free to enter and leave the computational cell.CP2K, 10,23 CP2K+SMEAGOL is the most feature-rich and complete to-date.The support for DFT-NEGF natively in CP2K is limited to Γ-point only transport calculations, without forces required for geometry optimisation or molecular dynamics. 10e remainder of the paper is organised as follows.Section 2 summarises briefly the theory of DFT-NEGF and the implementation in CP2K+SMEAGOL, followed by results at both zero bias (Section 3.1) and finite bias (Section 3.2-3.3).We also show that CP2K+SMEAGOL can be used to perform molecular dynamics of large solvated systems (Section 3.4), and discuss the performance and as well as methodologies that may be employed to accelerate the dynamics in future work (Section 3.5).Concluding remarks are made in Section 4.
A representation of the typical arrangement used in a DFT-NEGF calculation is shown in Fig. 1, composed of a central region termed the 'extended molecule' (EM) attached to two semi-infinite electrodes or 'leads'.The two leads are kept at different chemical potential by a battery and are able to exchange electrons with the extended region.
The retarded Green's function G(ϵ), referred to as the Green's function in the remainder of this work, is obtained by performing an inversion of the DFT Hamiltonian matrix H combined with the overlap matrix S and the leads self-energies Σ L and Σ R , 5,24 where δ + is an infinitesimal positive number.The self energies are in general non-hermitian matrices that contain information about the electronic structure of the leads and possibly the bias.In principle the leads self-energies Σ L and Σ R can differ, however in practice this is challenging to converge and therefore for all calculations performed in this work the left and right leads are identical.
Similarly to SIESTA+SMEAGOL, the Hamiltonian and the non-equilibrium density of the central region are calculated self-consistently using the NEGF scheme, with boundary conditions set by self-energies that are static and independent of the charge density in the central region.This assumption holds if the central region is sufficiently large, ensuring that changes in charge density are screened before reaching the boundaries.For the metallic leads typically used, which have a short screening length, this requirement is met within a few atomic layers.In calculations, these layers need to be included at the boundaries of the central region to extend the leads.With the NEGF method, the currents entering and leaving the EM are balanced, allowing to evaluate the electron density in the EM, the electrostatic potential and the current.
The electron density is split into left and right contributions D = D L + D R , calculated as an integral of the Green's function where f is the Fermi-Dirac distribution of the left lead described by the chemical potential µ L , k B the Boltzmann constant, T L is the electronic temperature of the left lead and ρ L (ϵ) is the spectral density matrix calculated from the Green's function and the broadening function Γ L of the left electrode, The calculation of the electron density as an integral of the Green's function across the entire energy space in Eq. 2 would be very demanding, and therefore the integral is split further into an equilibrium contribution D L eq , which can be integrated on a complex contour, and a non-equilibrium contribution ∆ R neq , which needs to be integrated along the real axis.Notably, the non-equilibrium contribution to the density only needs to be evaluated for energies within the bias window.The total density is therefore written as where and We could re-write Eq. 5 as D = D R eq + ∆ L neq , where L and R are exchanged.The original version of SMEAGOL calculates the density matrix using either approach, or as an average of both. 8This approach however does not correctly occupy bound states within the bias window, discussed further in Section 3.2.As such we have updated SMEAGOL to more closely follow the work of Brandbyge et al., 24 calculating the density matrix as a weighted sum of the two integrals where Following the self-consistent calculation of the non-equilibrium density a number of transport quantities can be calculated using the DFT-NEGF approach, such as energy dependent total transmission coefficient 8,25 and the current, calculated using the Landauer-Buttiker formalism 26 An electric current can significantly change the forces experienced by a nano-system, and therefore its dynamics.This is at the origin of electromigration phenomena and of the rupture of nanowires under bias conditions.These forces are non-conservative, 27 thus the way the electric current affects the dynamics of an atom can vary along the trajectory of that atom, causing instabilities.The forces are calculated via the time derivative of the atomic momenta, resulting from the average of the gradients (with respect to atomic positions) of the full Hamiltonian, according to the expression: The Hamiltonian in this expression is calculated using the non-equilibrium density obtained by solving the NEGF problem, which embed instantaneously the information about the stationary charge redistribution induced by current and bias.We use these non-conservative forces to perform Born-Oppenheimer molecular dynamics under bias and currents.We would like to highlight that this expression for the force acting on ions is generally different from the Hellmann-Feynman expression, F i = −∇ R i ⟨Ψ| Ĥ|Ψ⟩.However, according to Rungger et al. 28 and Di Ventra et al. 29 if the state is stationary (with the only time dependence being through the exponential term) and is square integrable, the two expressions become formally identical.
Another consideration for any DFT-NEGF implementation is to correctly express in term of the Green functions, the force term arising from the derivative of the density matrix, where, I is the ion index, K ij is the Kohn-Sham matrix, Ω ij is the energy weighted density matrix, and ϕ i/j are localised orbital functions.Specifically, the Ω ij matrix is calculated in the same manner as the density matrix (Eq.2), with an additional energy term ϵ in the integral. 28tably, this term is often neglected in existing implementations of current-induced forces, 32 however without it the forces are incorrect (Section 3.1).This implementation enables the modelling of the 'electron wind' component of the current-induced force, which models the transfer of momentum from electrons to ions under fixed boundary conditions, as discussed by Dundas et al. 33 and Zhang et al. 28 We neglect any further electron-ion coupling, including coupling with electrode phonons.We assume that we are in the limit where the fluctuating forces leading to Joule heating ?and the phase change of the electronic wave functions due to atomic motion ?are not significant.

Results
We present validation of CP2K+SMEAGOL at both zero bias (Section 3.1) and finite bias (Section 3.2-3.3),for a variety of systems.We also show that CP2K+SMEAGOL can be used to perform molecular dynamics of large solvated systems (Section 3.4), and discuss the performance and as well as methodologies that may be employed to accelerate the dynamics in future work (Section 3.5).

Zero bias forces
A simple validation of our CP2K+SMEAGOL interface is to ensure that the forces calculated at zero bias with CP2K+SMEAGOL are equal to the forces calculated only with CP2K.We use the same system setup as Zhang et al., 28 an infinite Au chain composed of 9 atoms.The central Au atom is displaced by 1 Å in the +x direction, such that there is a restoring force acting towards the equilibrium position.The left lead, extended molecule and right lead are each composed of 3 Au atoms.The LDA functional is used, with a single-ζ basis set for the Au(6s, 5d) electrons.The same system setup is used to perform reference calculations in SIESTA+SMEAGOL, using an explicit cutoff radius for Au(6s) of 6.5 Bohr and Au(5d) of 5.85 Bohr.Unless otherwise stated, the same computational setup is used for calculations in this work.
Fig. 2 shows the x component of the force calculated with both CP2K and CP2K+SMEAGOL at zero bias, as well as reference calculations with SIESTA and SIESTA+SMEAGOL.While the exact value of the force differs between CP2K and SIESTA, they are both consistent with the forces calculated with SMEAGOL at zero bias.We also demonstrate the effect of neglecting the forces calculated from the energy density matrix term Ω, as performed in some DFT-NEGF codes, 32 resulting in qualitatively incorrect forces.

Parallel-plate capacitor
Another key validation of our CP2K+SMEAGOL interface is to reproduce the potential drop across a parallel-plate capacitor, a popular benchmark system for DFT-NEGF implementations. 8,24The structure used is shown in Fig. 1 Fig. 3 also demonstrates that the Hartree potential calculated using a single contour evaluation of the Green's function (Eq.5) is qualitatively incorrect at high bias, with asymmetric charge density at the left and right electrodes.Using our newly implemented weighted double contour (Eq.8) scheme in SMEAGOL we produce a qualitatively correct Hartree potential, with symmetric charge density on the left and right electrode.Brandbyge et al. 24 attributed this difference to the presence of bound states within the bias window that only couple to one electrode, and therefore a single contour evaluation does not correctly populate these states.As the use of the weighted double contour is only necessary at high applied bias, generally above 2 V in our experience, we only use a single contour evaluation of the Green's function for all other calculations in this work.

Geometry optimisation of an Au-H 2 -Au junction
38][39]39 Having demonstrated that CP2K+SMEAGOL can reproduce forces and charge density consistent with SIESTA+SMEAGOL, we perform calculations for the Au-H 2 -Au junction based on the work of Bai et al. 38  Fig. 4 shows the mean displacement of the unconstrained atoms against the current flow, and the elongation of the H-H bond for both CP2K+SMEAGOL and SIESTA+SMEAGOL.
While the exact H-H bond lengths are different, the elongation as a function of an applied bias up to 1.5 V bias is consistent.The minor differences in the absolute bond lengths can be attributed to the difference in the basis sets used for the calculations.By comparing the results to a hydrogen molecule in vacuum in the presence of an effective electric field, Bai et al. 38 confirmed in their work that the observed geometrical changes are are dominated by the electric current rather than by the electric field.
We have verified that there is no significant difference in the geometry when the Green's function is evaluated using a weighted double contour (Eq.8) instead of the single contour (Eq.5) evaluation used in the original work of Bai et al. 38

Molecular dynamics of a solvated Au wire
Atomic size metallic contacts have been studied extensively, [40][41][42][43][44] with interest both in their fundamental properties as well as for their possible applications in nano-electronics.An important application particularly relevant to this work is in high density memories and  logic applications, e.g. in molecular switches and electrohemical metallisation memories, [45][46][47] where the non equilibrium component of the force is essential to describe the change of status of the system.
Monoatomic Au wires are typically formed through experimental techniques such as scanning tunnelling microscopy (STM) 48 and the mechanically controllable break junction (MCBJ), 49,50 however most studies are performed in vacuum and therefore neglect any solvation effects.In this work we present CP2K and CP2K+SMEAGOL MD calculations for a fully solvated monoatomic Au wire composed of a 4-atom wire, connected via two 3-atom pyramidal tips to two Au(111) slabs.
Starting from an initial structure with a straight Au wire forming Au-Au-Au bond angles of 180°, during DFT cell optimisation performed in vacuum the Au wire relaxes to form a zigzag geometry 42 with an Au-Au-Au bond angle of around 130°.The Au wire was then fully solvated by adding 166 water molecules, equilibrated by performing 30 ns classical MD with a TIP3P water model 51 and Lennard-Jones potentials. 52The Au atoms were then allowed to relax through 9 ps of NVT DFT-MD, where the Au wire straightened and one Au wire atoms was forced into the tip geometry.The PBE functional was used, with triple-ζ plus polarisation functions for the Au(6s, 5d) electrons, O(2p, 2s) electrons and H(1s) electrons.
To decrease the significant computational cost of subsequent CP2K+SMEAGOL calcula- Fig. 6 shows the transmission calculated for the solvated Au wire sampled across 500 fs of DFT-MD, with comparison to the transmission of an ideal infinite Au wire.Qualitatively the transmission is very similar, roughly equal to 1 between -2 eV and 2 eV.The differences in transmission between the solvated wire and the system where water was removed (red and green curves in Fig. 6 ) can be ascribed to the opening of additional tunnelling channels for the electron through the junction when the water is present.More specifically, these channels become available in the same energy range where the water band arises, with the density of states shown in ESI Fig. 5.The peak in transmission for the ideal Au wire between -2 and 0 eV is not present for the solvated Au wire, as this peak is due to Au(5d) states that are neglected in the leads atoms for the solvated Au wire.
We also perform CP2K+SMEAGOL molecular dynamics at an applied bias of 0 V, 0.1 V and 1 V.While the total energy for such a current-carrying open system may not be conserved, 27 in practice we that the dynamics are stable with minimal long-term energy drift.See ESI Section 1.3 for further discussion of energy conservation for molecular dynamics performed for both the solvated Au wire (Section 3.4) and the Au-H 2 -Au junction (Section 3.3).

Performance
The main bottleneck in a CP2K+SMEAGOL calculation is the evaluation of the density as an integral of the Green's function (Eq.2), as many computations of the Green's function be performed under bias, for performing molecular dynamics calculations it becomes problematic.As such, we are currently investigating methods that can be used to accelerate or circumvent the need for these calculations.The use of ScaLAPACK could allow the scaling limit to be improved, using distributed memory parallelism for the matrix inversion and matrix multiplication routines.Alternatively, it could be possible to use CP2K+SMEAGOL to generate training data for a machine learning force field. 57,58Work towards this is ongoing in our group.

Conclusion
In this work we have developed a new interface between the popular DFT package CP2K and the NEGF code SMEAGOL, allowing for CP2K calculations to be performed under an applied potential and current flow.In contrast to the existing DFT-NEGF implementation available natively in CP2K, both k-point sampling and forces are available.
We have verified and benchmarked our interface against systems previously studied with SMEAGOL and the DFT code SIESTA, showing good agreement for both single point calculations and geometry re-optimisation under bias for: an infinite Au wire, a parallelplate capacitor and a Au-H 2 -Au junction.
For the first time we are able to perform large-scale molecular dynamics simulations under realistic operating conditions using the DFT-NEGF approach.For an example of a solvated Au wire, we observe electron migration effects during a short 200 fs MD trajectory at an applied bias of 1 V.We expect that this methodology will become valuable for the emerging field of first principles electrochemistry, for example to model the effect of an applied potential or current flow through an electro-chemical cell.

Figure 1 :
Figure 1: Schematic view of a Au capacitor used for CP2K+SMEAGOL and SIESTA+SMEAGOL calculations.

Figure 2 :
Figure 2: Zero bias tests for an infinite Au wire.(A) Atomic forces calculated with CP2K and CP2K+SMEAGOL.(B) Atomic forces calculated with SIESTA and SIESTA+SMEAGOL.The y and z components of the atomic forces are shown in ESI Fig. 1.

Figure 3 :
Fig.3shows the planar average of the Hartree potential difference across the parallelplate capacitor with and without an applied bias.The CP2K+SMEAGOL and SIESTA+SMEAGOL

Figure 4 :
Figure 4: Finite bias geometry optimisation of a Au-H 2 -Au junction.(A) Optimised atomic structure of the Au-H 2 -Au junction, where yellow and white spheres represent Au and H atoms respectively.(B) Mean displacement of the six highlighted atoms as a function of the applied bias.(C) Change in the H-H bond length as a function of the applied bias.
is shown in Fig.4, composed of a hydrogen molecule sandwiched between two 3x3 Au(100) slabs with 7 layers (left) and 6 layers (right).All atoms are constrained in the system except for the 2 H atoms and the 6 Au atoms either side.The PBE functional is used, with a single-ζ basis set for the Au(6s, 5d) electrons and a double-ζ basis set for the H(2s) electrons plus polarisation functions.For the reference SIESTA+SMEAGOL calculations we use an explicit cutoff radius for Au(6s) of 6 Bohr, Au(5d) of 5.5 Bohr, H(2s) of 5.5 Bohr and H(1p) of 1.5 Bohr.The k-point sampling for the semi-infinite leads is 3x3x20 and for the extended molecule is 3x3x1.

Figure 5 :
Figure 5: Molecular dynamics of a solvated Au wire using CP2K+SMEAGOL at an applied bias of 1 V. (A) DFT equilibrated starting geometry of the solvated Au wire, where yellow, white and red spheres represent Au, H and O atoms respectively.The basis set used for the Au atoms is shown above the structure, with each distinct Au region separated by a dotted black line.(B) Bader charges for the Au wire atoms, Au tip atoms and the first layer of the Au slabs calculated on the DFT equilibrated geometry as well as after 100 fs MD and after 200 fs MD.

Figure 6 :
Figure 6: Average transmission calculated using CP2K+SMEAGOL at zero bias.Solvated Au wire (red) and with water molecules removed (green) sampled every 100 fs across 500 fs of DFT-MD, with comparison to the transmission for an ideal infinite Au wire in vacuum (black).The shaded region represents 1 standard deviation of the mean, shown as a solid line.
tions it is necessary to use a minimal single-ζ Au(6s) basis set for the leads atoms, which has been shown in previous work to reproduce key transport properties such as the transmission and current up to an applied bias of around 2 V.53,54 In addition the triple-ζ basis set for the solvated Au wire is reduced to a double-ζ basis set, and to avoid spurious density buildup at the interface between the different basis sets of the extended molecule and the leads additional screening regions are added.See ESI Section 1.3 for further details regarding the basis set choices and size of the leads.All Au atoms are frozen except for the 3-atom Au wire and the 7 atoms belonging to the pyramidal Au tips.With this new configuration 1.5 ps of NVT DFT-MD is performed to re-equilibrate the system, resulting in the structure shown in Fig.5.

Fig. 5
Fig.5shows the Bader charges of the Au wire atoms, as well as the Au tip and first layer of the Au slabs.The Bader charges are calculated as the difference in the reference charge and the atomic contribution to the electron density calculated using the algorithm developed by Henkelman et al,55 in units of the elementary charge |e|.With both CP2K and CP2K+SMEAGOL at zero applied bias the Au atoms belonging to the Au wire are negatively charged (-0.30, -0.32, -0.27), as well as the first atom of the pyramidal Au tips (-0.19, -0.18).The remaining Au atoms of the tips are slightly positively charged: +0.09 to the left of the Au wire, +0.05 to the right of the Au wire, while the first layer of the Au slabs are approximately neutral: -0.03 and -0.01 respectively.With an applied bias of 0.1 V there is no significant change in the Bader charges, while at an applied bias of 1 V the charge on the leftmost Au wire atom increases from -0.19 to -0.16 (+0.03), and the charge on the rightmost Au atom decreases from -0.18 to -0.20 (-0.02).Performing 200 fs CP2K+SMEAGOL MD at an applied bias of 1 V we observe electron migration effects, where the minima in the Bader

Figure 7 :
Figure 7: Speedup of CP2K+SMEAGOL as a function of the number of OpenMP threads for different system sizes.Calculations are performed for a model Au-BDT-Au junction, where the number of atoms is increased from 948 (3940 basis functions) to 3756 (15388 basis functions).