User-Defined Electrostatic Potentials in DFT Supercell Calculations: Implementation and Application to Electrified Interfaces

TL;DR

This paper introduces a flexible method to apply custom electric fields in DFT calculations, enabling better study of electrified interfaces and electrochemical processes.

Contribution

A new implementation for applying arbitrary electrostatic potentials in DFT supercell calculations using the VASP–Python interface.

Findings

The implementation allows for energy and force corrections in electrified DFT simulations.

The method was successfully applied to various systems like electrified surfaces and solid-water interfaces.

The approach provides flexibility and control within standard VASP software.

Abstract

Introducing electric fields into density functional theory
(DFT)
calculations is essential for understanding electrochemical processes,
interfacial phenomena, and the behavior of materials under applied
bias. However, applying user-defined electrostatic potentials in DFT
is nontrivial and often requires direct modification to the specific
DFT code. In this work, we present an implementation for supercell
DFT calculations under arbitrary electric fields and discuss the required
corrections to the energies and forces. The implementation is realized
through the recently released VASP–Python interface, enabling
the application of user-defined fields directly within the standard
VASP software and providing great flexibility and control. We demonstrate
the application of this approach with diverse case studies, including
molecular adsorption on electrified surfaces, field ion microscopy,
electrochemical solid–water interfaces, and implicit solvent
models.

Figures

Flowchart of the VASP-Python plugin setup during a Molecular Dynamics-DFT simulation. ΔF I , ΔE I , n e–, V ext and Z I refer to the forces, total energies, number of electrons, external potential and nuclear charges, respectively.

(a) Energy of a hydrogen atom in a constant potential V ext before and after the core correction (see text). The black dashed line represents the energy equal to 1e · V ext. (b) The force in the z direction of a hydrogen atom in a linear potential corresponding to a constant field Eext before and after the correction. The black dashed line represents the force equal to 1e·Eext .

Schematics of simulating charged surface using the (a) Ne computational counter electrode (CCE) setup and the (b) charge density counter electrode (CDCE) setup. The working electrode has a net charge of n electrode, which is compensated by the counter electrode, thus creating an electric field E across the simulation cell. The black lines represent the electrostatic potential Φ.

(a) Electrostatic potential of the Au(111) surface for several electrode charges; the red curve corresponds to the Gaussian charge-density counter electrode (CDCE). (b) Adsorption energies of a single H atom on Au(111) as a function of the applied electric field. (c) Charge-density difference for n electrode = 0, computed with eq and an isosurface value of 1.5 × 10–6 e/Å3. (d) Charge-density difference for n electrode = – 0.5 e –, computed with eq and an isosurface value of 1.5 × 10–7 e/Å3. Yellow denotes excess electronic charge (negative charge density), while cyan denotes a deficit of electronic charge (positive charge density).

AIMD simulations of an acetaldehyde molecule adsorbed on the Au(111) surface under applied field, showing the dependence of Au–O distance with n electrode.

(a) Initial and (b) final snapshot of a 5 ps ab initio molecular-dynamics (AIMD) trajectory of the (952) kinked Li surface at 300 K under an applied electric field of 1.25 V Å–1. The atom colored white is the kink atom, which is emitted to an adatom position after ≈ 4.5 ps of simulation time. Lighter atoms correspond to surface atoms.

(a) Au(111)/H2O slab with the CDCE setup. (b) Time evolution of the electrode charge and the corresponding electrode potential. The target potential Φ0 – Φzc (horizontal line) is set to – 2.0 V. (c) Oxygen density profiles obtained from AIMD simulations of Au(111)/H2O under different applied voltages and at the potential of zero charge.

(a) Contour plot of the solvent charge density ρsolvent MM of a Mg2+ ion solvated in water. (b) Comparison between the electrostatic potential of a single Mg2+ in vacuum and in the solvated QM/MM method.

Full text

Introduction

1

Simulating the response of molecules and materials to external electrostatic potentials and fields is of fundamental importance to chemistry, physics, biology and materials science.? External electric fields may originate from various environmental sources, including solvents, ?,? electrodes, ?−? ? charged tips from measuring instruments such as in atom probe tomography,? from biological membranes and ion channels. ?,? The role these fields play is critical, as they can drive chemical reactions, ?−? ? induce modifications to surfaces,? influence the folding of proteins? and determine the selectivity of reaction pathways. ?,?

Density Functional Theory (DFT) is commonly used to capture the quantum mechanical behavior of molecules and materials under the effect of an external field, due to its favorable balance between computational cost and accuracy. ?,? DFT calculations can be performed with a variety of electronic structure codes, such as the Vienna Ab Initio Simulation Package (VASP), ?,? ORCA,? QUANTUM ESPRESSO,? Molpro,? GPAW,? S/PHI/nX ?,? and many others. These electronic structure codes support the use of external fields and potentials with varying flexibility, from the use of point charges, external plugin files, in-code commands or relying on direct modification of the base code from the user. Well-established, built-in uses of external potentials include the dipole correction? or implicit solvent methods. ?−? ? Beyond these, ongoing efforts seek other external field capabilities, such as for modeling solvated electrified surfaces in electrochemistry ?−? ? or solvation effects when referring to quantum mechanics/molecular mechanics (QM/MM) approaches. ?,?

Implementing these external electrostatic potentials in practice has often required direct modification of the source codes of electronic structure packages by users or research groups. Although such efforts have led to progress on specific scientific applications, they demand intimate knowledge of large and complex code bases and are additionally prone to the introduction of errors, difficult to reproduce or maintain across software versions. Custom patches hinder portability and subsequent extensions by others, and ad hoc edits to individual subroutines can inadvertently cause bugs and regressions in other parts of the code. These challenges highlight the need for modular, well-documented interfaces between the source code and the user, that allow for the interior modification of the complex electronic structure codes, while insulating their core functionality.

The freedom to modify selected components and quantities in the source code, gives rise to a high flexibility of computational setups and systems that can be explored. Nevertheless, careful validation of the computational setup is essential to ensure correct electrostatic boundary conditions, in addition to the physical consistency of the forces and energies under the effect of an external potential. In practical terms, this entails ensuring charge neutrality and avoiding simulation artifacts due to long-range Couloumb interactions, thermodynamic inconsistencies and so forth. Additionally, correction terms must be added to the forces and energies, due to their electrostatic interactions with the external field. Lastly, it is important to note that upon modification of the software package, the responsibility of ensuring numerical and physical correctness shifts from the software developer to the user.

In this context, we discuss the theoretical and computational details for the setup of arbitrary electrostatic potentials in DFT codes while ensuring correct forces, energies and electrostatic boundary conditions. For this purpose, we employ the recently released VASP 6.5.0 version,? which allows for the straightforward modification of the local potential and other internal quantities within the source code through a Python scripting interface. We make use of these features to apply electric fields across interfaces, to control the voltage in electrochemical simulations, and to impose an external solvation potential. These methodologies enable a wide range of applications, such as the modeling of solid–liquid interfaces, electrocatalysis, solvation effects and materials. Beyond the physical and numerical consistency, the proposed computational setup between the Python plugin and the electronic structure source code ensures robustness, modularity and extensibility.

Theoretical and Computational Details

2

Workflow of the VASP-Python Interface

2.1

The Python interface introduced in VASP provides direct access to the local potential, nuclear charges and forces, as well as total energies. These quantities can be modified by the Python script during initialization and dynamically during the ionic loop for an energy minimization or molecular dynamics run. A schematic overview of the workflow between the electronic structure code and the Python plugin is shown in Figure. For the purpose of applying an external electrostatic potential, V ext is added to the total potential during the electronic loop. As such, the charge density n′(r, V ext) under the applied bias is calculated by the VASP code. However, by design of the VASP-Python interface, the interaction between V ext and the nuclear cores are not included automatically. In order to obtain the correct energies and forces, these contributions (force correction ΔF _ I _ and energy correction ΔE _ I _, where the subscript I represents the Ith atom) should also be added through the Python plugin, which we discuss in Sec. 2.2.

*Flowchart of the VASP-Python plugin setup during a Molecular Dynamics-DFT simulation. ΔF

I , ΔE

I , n

e–, V ext and Z

I refer to the forces, total energies, number of electrons, external potential and nuclear charges, respectively.*

Additionally, the VASP-Python plugin enables changes to the number of electrons (n _ e _−), as well as the nuclear core charge (Z _ I _), for each ionic step. This enables the application of an external electric field that changes with time during an ab initio molecular dynamics (AIMD) run. An important example where this functionality is useful is potential-controlled simulations of electrochemical interfaces, which we discuss in more detail in Sec. 2.3.

In this work, we refer only to quasi-static fields that can be modeled within the framework of the Born–Oppenheimer approximation. This contrasts with time-dependent oscillating fields, such as those due to electromagnetic waves, which are relevant for simulating spectroscopic phenomena.

Energy and Force Corrections

2.2

In this section, we derive the modifications of the total energy and forces under the external electrostatic field and identify the terms that need to be included through the VASP-Python plugin. The derivation is adapted from applying the dipole correction to DFT calculations.? The energy functional without V ext can be written as

$$[eqn]$$

where T[n], E ^ e–e ^[n], and E ^ion–ion^ represent the kinetic, electron–electron, and ion–ion interaction energy. V ^ion^ is the ionic pseudopotential. Upon the application of V ext, the energy functional changes to

$$[eqn]$$

Here n ^ion^(r) represents the density of the nuclear core charges. They are represented as point charges in VASP.

The two terms added in eq represent the interaction between the electron charge with the external field and the nucleus charge with the external field, respectively. When V ext is added to the total potential in VASP, n(r) relaxes self-consistently to n′(r, V ext). As such, VASP also computes the energy functional as E[n′(r)], including the effect of V ext to all terms in eq. However, for the two additional terms in eq, by default the VASP output energy includes only the interaction between V ext and the electron density but not the nucleus. Therefore, an energy correction needs to be added for the Ith nucleus as

$$[eqn]$$

Similarly, in the VASP force output, the interaction between V ext and the nucleus charge is also lacking, and a force correction should be added to each nucleus

$$[eqn]$$

where $[eqn]$ represents the externally applied electric field corresponding to V ext. We note that the derivations here for incorporating the external potential into the Hamiltonian is general to all DFT codes, not limited to VASP.

To benchmark the proposed energy and force correction scheme, we perform tests on the model system of a neutral hydrogen atom in vacuum (Figure). In Figurea, a constant potential is added through the VASP-Python interface. Since the system is charge neutral, it is expected that the energy of the system does not change. However, the uncorrected energy (blue line) shows a linear dependence on V ^ ext ^, which scales with 1e · V ^ ext ^ (black dashed line). This behavior matches eq, as the Z _ I _ value for hydrogen is 1. After the energy correction, the energy of the hydrogen atom is invariant with respect to V ^ ext ^.

(a) Energy of a hydrogen atom in a constant potential V ext before and after the core correction (see text). The black dashed line represents the energy equal to 1e · V ext. (b) The force in the z direction of a hydrogen atom in a linear potential corresponding to a constant field Eext before and after the correction. The black dashed line represents the force equal to 1e·Eext .

Similarly, in Figureb, we test the force correction by putting the hydrogen atom in a constant electric field $[eqn]$ . As the hydrogen atom is charge-neutral and its dipole moment due to field polarization is negligible, the force should be close to zero. However, again we observe that the uncorrected force (blue line) has a strong linear dependence on $[eqn]$ . The dependence follows the line $[eqn]$ (black dashed line), as derived in eq. These tests demonstrate that the interaction between the nuclear core charges and the external potential is indeed not included in the VASP-calculated energies and forces, and that the corrections in eqs and ? need to be applied.

Simulating Electrified Surfaces with Potential

Control

2.3

An important application of field-dependent DFT calculations is to simulate electrified surfaces and interfaces.? This is achieved by setting up a counter electrode which compensates the charge on the electrified model surface, thus including an electric field across the slab system. For a discussion of different approaches to set up a counter electrode charge we refer to ref ? Here we use a neon (Ne) computational counter electrode (CCE), which has been developed and used in previous works. ?,? It allows a precise adjustment and control of the counter electrode charge by changing the nuclear core charge of Ne. The CCE model setup is shown schematically in Figurea. We denote the net charge on the working electrode as n electrode. In the CCE setup, the Ne counter electrode carries a net charge of – n electrode and therefore the cell is charge neutral. This configuration allows us to induce a tunable electric field $[eqn]$ . The magnitude of the field can be controlled by the position of the CCE and/or the charge n electrode on the electrode.

Schematics of simulating charged surface using the (a) Ne computational counter electrode (CCE) setup and the (b) charge density counter electrode (CDCE) setup. The working electrode has a net charge of n electrode, which is compensated by the counter electrode, thus creating an electric field E across the simulation cell. The black lines represent the electrostatic potential Φ.

In this work, we introduce a computational setup beyond CCE, which we term as the charge density counter electrode (CDCE). ?,? In the CDCE scheme, instead of having explicit atoms for the counter electrode, we introduce an external counter charge density, ρ_ext_, with total charge equal to – n electrode. ρ_ext_ is constructed as a sheet charge with a Gaussian distribution in the z-direction. In practice, we calculate the corresponding V ext of ρ_ext_ by solving the Poisson equation with fast Fourier transformation (FFT)

$$[eqn]$$

V ext is then added to the total potential through the VASP-Python plugin. As such, the – n electrode charge is not accounted for in the DFT self-consistency loop. The total number of electrons in the system is then set to

$$[eqn]$$

Here Σ_ I _ Z _ I _ sums over the nuclear core charges of all atoms in the system, and is equal to the total number of electrons in a charge-neutral cell. This setup establishes the same $[eqn]$ across the cell as CCE (Figureb).

The major advantage of going from CCE to CDCE is to increase the range of electric field strength that can be applied. In the CCE method, the maximum field strength applicable is constrained by the Ne band gap. When a sufficiently large field is applied across the cell, the Fermi level of the system goes below the valence band maximum or above the conduction band minimum of Ne, resulting in dielectric breakthrough. This constraint is eliminated when the Ne atoms are replaced by a Gaussian charge density. However, we note that other constraints exist for the maximum field strength, for example the band gap of the electrolyte and the alignment of the Fermi level with respect to the vacuum level. A more detailed discussion can be found in ref ?.

Electrochemical experiments are routinely carried out under constant electrode potential. To correctly represent the experimental setup in AIMD simulations, it is essential to define realistic electrostatic boundary conditions. For a potentiostat, this includes the thermal fluctuations in the electrode charge. One such approach is the thermopotentiostat technique proposed by Deißenbeck. ?,? The key ingredient of the method is to introduce the fluctuation term in the electrode charge n electrode, to mimic the macroscopic constant-potential dynamics within the small DFT cell. Given an electrode charge at a certain time step n electrode(t), the electrode charge at the next time step, n electrode(t + Δt) is given by

$$[eqn]$$

where C 0 is the bare capacitance of the electrodes in vacuum, Φ­(t) is the instantaneous potential measured at the dipole correction, Φ_0_ is the target potential, τ_Φ_ is the relaxation time constant, T is the temperature and k B is Boltzmann’s constant. N is a random number with zero mean and variance one. Equation essentially sets n electrode at every time step to achieve Φ_0_ on average, with the consistent fluctuations according to statistical physics.

In ref ?., the thermopotentiostat calculations were performed by directly modifying the VASP source code. In this work, we reimplement the method using the VASP-Python interface to enable such simulations with the standard VASP version. Additionally, we have extended the thermopotentiostat method based on the Ne CCE to the CDCE as well.

Implementation in VASP

2.4

In this section, we provide a practical guide on how to use the VASP Python plugin to apply an arbitrary electric field in a simulation. This is done by creating an additional input file for a VASP calculation named vasp_plugin.py. This Python file provides callback functions that can access certain quantities during the simulation and allow the user to modify some of them. For the purpose of this work, we use three of these functions.

In each of these functions, constants provides a dataclass that contains given calculated quantities during the VASP run. These quantities are accessible to the user, but cannot be changed. additions, on the other hand, contains quantities that can be changed by the user. The list of quantities included in each dataclass can be found in the VASP documentation. For example, in the local_potential function, we have constants

So in the local_potential function in vasp_plugin.py, we can access the energy cutoff (ENCUT) with constants.ENCUT (similarly for all other accessible quantities), and we can change the total potential with

To call these functions defined in vasp_plugin.py, one also needs to set the corresponding flags to true in the INCAR file:

In the following sections, we briefly explain the role of each function and discuss related technical details.

local_potential

2.4.1

The task of the local_potential routine is to modify the electrostatic potential that VASP uses internally. The modification is performed by adding an external contribution V ext to the existing potential:

Here delta_potential is a three-dimensional array that must have the same mesh as VASP’s electrostatic potential and charge density (dimensions NGX × NGY × NGZ). The local_potential function has access to the charge density constants.charge_density and the dimensions can be directly obtained.

In addition to updating the potential, two further corrections are required. First, VASP does not include the interaction of the external potential with the nuclear charges, so the corresponding energy term has to be added explicitly:

The quantity delta_energy is defined as

$$[eqn]$$

where the sum runs over all nuclei I with charge Z _ I _ located at R _ I _.

Second, special care is needed when a dipole correction is required. Consider the CDCE configuration illustrated in Figureb. In this setup V ext represents the field generated by the counter-electrode charge – n electrode. The total number of electrons that VASP should treat is

$$[eqn]$$

with n neutral being the electron count of a neutral cell.

Because VASP does not contain the charge density of the counter electrode, it interprets the system as non-neutral. If the dipole correction is activated in VASP (via the tag LDIPOL = .TRUE.), the program aborts, since the dipole correction is only implemented for charge-neutral cells. To circumvent this limitation we disable VASP’s built-in dipole correction and implement it manually inside the local_potential routine. The contribution of the counter-electrode charge is computed in the helper function calc_dipole and added to delta_potential before the potential update:

Here rho is the charge density in VASP (excluding the counter electrode charge density), rho_external is the charge density of the externally added counter electrode. The calc_dipole function then computes the total dipole of the system and outputs the dipole charge density that should be added in the vacuum region rho_dipole. The electrostatic_potential function then computes the corresponding potential profile of the charge density, which is added to the total potential. By doing so the dipole correction is applied consistently, even for simulations that involve a net charge introduced by the external electrode.

Force_and_stress

2.4.2

The force_and_stress routine implements the force correction that appears in eq, i.e., the force exerted on each nucleus by the externally applied field. For the CDCE configuration a short-range repulsive “wall” is also added near the counter-electrode to keep the water molecules confined; this wall force is obtained from the helper function wall. Both contributions are combined and added to the forces that VASP returns:

In this work, we consider only calculations with fixed cell size, so the stress contributions are not considered.

Occupancies

2.4.3

The occupancies routine is called at the end of every ionic step. It enables a time-dependent bias during a molecular-dynamics run, as required by the thermopotentiostat scheme. In the charge-density counter-electrode setup (Figureb), the routine performs two tasks.

• (i) Update of the electrode charge The change in the electrode charge, Δn electrode ≡ dq, for the next step is evaluated from eq: Here C0 is the system capacitance, phi0 the target bias, and phi the instantaneous bias computed at the current step. The value of dq is stored so that the local_potential function can read it and update the external potential V ext.
• (ii) Adjustment of the total electron count Because the electrode charge changes, the total number of electrons in the simulation must be modified in accordance with eq. This is done by updating the NELECT variable:

Additional Development Features

2.4.4

Two capabilities that were employed in this work are available only in the development branch of VASP and have not yet been released in the official distribution:

• 1. Ne computational counter electrode. Instead of modifying V ext and NELECT, the number of core electrons of neon (ZVAL) is changed. In VASP 6.5.0 the ZVAL field is not part of the additions dataclass used by occupancies, but this will be added in a forthcoming release.

Force-drift correction. By default VASP applies a drift-correction that forces the net force on all atoms to vanish. When an external potential is present this constraint is no longer appropriate. In the development version the drift correction can be disabled with the flag This option will also become available in the next official VASP release.

Case Studies

3

By utilizing Python scripting within VASP, we are able to apply an arbitrary external electric field in a DFT calculation. In this section, we present four case studies utilizing the implemented method: (1) adsorption on electrified surfaces; (2) field ion microscopy; (3) electrochemical interfaces; and (4) implicit solvation model. These examples highlight the range of field-induced physical and chemical processes that can be modeled with DFT and the flexibility and robustness of our implementation.

Adsorption on Electrified Surfaces

3.1

To demonstrate the capability of our workflow to treat electrified interfaces, we use the Au(111) surface as a benchmark system and impose a series of electrode charges n electrode. Each chosen n electrode generates a corresponding electrostatic bias (or voltage) Φ, which can be obtained directly from the macroscopic dipole moment of the slab

$$[eqn]$$

where μ_dip_ is the dipole moment of the simulation cell, ϵ_0_ the vacuum permittivity, and A the surface area of the cell. The laterally averaged electrostatic potentials resulting from the Au(111) calculations are shown in Figurea; the red Gaussian-shaped curve represents the CDCE. For n electrode = −0.5 e ^–^ a voltage of Φ = 5.8 V is induced in the Au(111) slab. This approach for applying a bias to a DFT calculation is readily transferable to any slab system: the position, shape, and magnitude of the external charge density ρ_ext_ (and thus n electrode) can be tuned, provided the resulting field does not exceed the dielectric breakdown limit. ?,?

(a) Electrostatic potential of the Au(111) surface for several electrode charges; the red curve corresponds to the Gaussian charge-density counter electrode (CDCE). (b) Adsorption energies of a single H atom on Au(111) as a function of the applied electric field. (c) Charge-density difference for n electrode = 0, computed with eq and an isosurface value of 1.5 × 10–6 e/Å3. (d) Charge-density difference for n electrode = – 0.5 e –, computed with eq and an isosurface value of 1.5 × 10–7 e/Å3. Yellow denotes excess electronic charge (negative charge density), while cyan denotes a deficit of electronic charge (positive charge density).

Identifying active sites and reaction mechanisms is central to the development of electrocatalysts. Surface adsorption under an applied potential reveals the thermodynamics of the interface at a given electrochemical condition and is therefore key to uncovering reaction pathways.? As a showcase, we examine hydrogen adsorption on Au(111) using the CDCE approach. With Au(111) as the reference surface and a single H atom as the adsorbate, the adsorption energy under a bias Φ is evaluated as

$$[eqn]$$

where E tot ^Au/H^(Φ) and E tot ^Au^(Φ) are the total energies of the slab with and without the adsorbate, respectively, and E tot ^H_2_ ^ is the total energy of an isolated H_2_ molecule. The resulting hydrogen adsorption energies as a function of the applied field are displayed in Figureb. Three adsorption sites are considered: on-top, hcp-hollow, and fcc-hollow. At zero field the fcc-hollow site is the most stable, but the hcp-hollow site becomes energetically favored over the fcc-hollow site at an electric field of approximately 0.05 V Å^–1^.

The strength of these adsorbate–surface bonds is determined by the amount of electron transfer, which controls the kinetics of the adsorption and desorption of intermediates during the catalytic cycle.? Electron redistribution during the catalytic cycle is therefore critical to understanding catalytic performance.? To visualize and analyze these charge transfer processes, one can calculate the charge density difference between the adsorbate on the surface and the bare surface and isolated species. With the CDCE methodology, this is also possible under applied potential, where the electronic charge density difference is given by

$$[eqn]$$

with zc representing zero electrode charge, i.e., no external applied potential.

Figurec illustrates this charge density difference, which is dominated by the charge transfer of the surface to the H atom irrespective of Φ, meaning it is difficult to visually differentiate between the charge density difference calculated by eq when n electrode = −0.5 e ^–^ or n electrode = 0. Therefore, to isolate the charge transfer occurring upon applying an external bias from the charge transfer of the Au–H bond formation, we calculate the double charge density difference, given by

$$[eqn]$$

and shown in Figured. Here we visualize the polarization the Au–H bond under the effect of Φ. To quantify the charge transfer of eqs and ?, the charge density differences can be integrated as

$$[eqn]$$ $$[eqn]$$

Equation gives a charge of 1.29 and 1.28 e ^–^ for n electrode equal to zero and −0.5 e ^–^, respectively, while eq gives a charge 0.09 e ^–^ for n electrode = −0.5 e ^–^. This shows how the formation of the interaction of H with the Au dominates the charge transfer process; an additionally potential bias then only slightly further polarizes the surface–adsorbate. The positively charged gold surface polarizes the electron density along the Au–H bond toward the surface.

The CDCE methodology can be additionally used for AIMD, enabling a wide range of investigations. As a test system, we perform AIMD simulations to probe the interaction of an acetaldehyde molecule with the Au(111) surface under a range of electrode charges (Figure). The finite temperature simulation allows us to explore the diverse configurations of the organic molecule under applied field. The positively polarized surface, i.e., n electrode < 0, interacts strongly with the electronegative oxygen of the aldehyde, resulting in a small Au–O distance. As n electrode increases and the surface becomes negatively polarized, the oxygen is repelled. This results in the rotation of the acetaldehyde molecule with changing electrode charge. These simulations illustrate how our methodology can be used to explore dynamic systems under applied field, leading to possible mechanistic insights on reactions or interface processes.

AIMD simulations of an acetaldehyde molecule adsorbed on the Au(111) surface under applied field, showing the dependence of Au–O distance with n electrode.

Atom Probe tomography

3.2

Modeling electrified surfaces is of central importance for field ion microscopy (FIM) and atom probe tomography (APT), where surface atoms undergo field-induced desorption, ionization, and evaporation. Recent work on the Li(110) surface demonstrated that a strong electric field can alter the preferred adsorption site of Li adatoms and even enable barrier-less surface diffusion.? The same phenomenon can be captured directly in AIMD simulations that employ the CDCE setup.

Here we simulated a kinked Li(952) surface at 300 K. The top-view of the initial structure is shown in Figurea, with the kink atom highlighted in white. In the absence of an external field the kink atom remains bound to the step edge and no diffusion events occur. When a uniform field of 1.25 V Å^–1^ is applied, the kink atom spontaneously detaches and migrates across the surface within ∼4.5 ps of AIMD time, as illustrated in Figureb.

(a) Initial and (b) final snapshot of a 5 ps ab initio molecular-dynamics (AIMD) trajectory of the (952) kinked Li surface at 300 K under an applied electric field of 1.25 V Å–1. The atom colored white is the kink atom, which is emitted to an adatom position after ≈ 4.5 ps of simulation time. Lighter atoms correspond to surface atoms.

These calculations provide a direct, atomistic confirmation that strong electric fields can destabilize kink/step atoms and promote their release as adatoms at finite temperature. The observed behavior corroborates thermodynamic analyses predicting the destabilization of surface adatom states under high fields,? and offers valuable insight into the mechanisms governing APT experiments.

Electrochemical Interfaces

3.3

Thus, far we have described how to impose an external electric field on a vacuum-terminated slab. Extending the approach to solid–liquid interfaces requires additional care to avoid that the computational counter electrode affects water structure and dynamics. The Ne computational counter electrode (Ne-CCE) offers an inert, wide-band gap electrode that has been successfully employed to study Mg corrosion. ?,? However, its hydrophobic character limits the maximum field that can be sustained across the solid–liquid boundary. This limitation motivates the development of fully electronic counter electrodes, such as a Gaussian charge density analogous to the CDCE. A comprehensive discussion of the available approaches can be found in ref ?.

Figure illustrates the combination of the CDCE with the thermopotentiostat method for the Au(111)/water interface. During an AIMD run the charge on the gold surface, n electrode, is updated at every time step by the VASP–Python plugin according to eq. This feedback loop enforces a constant electrode potential, enabling truly potential-controlled AIMD simulations of solid–liquid interfaces. The resulting dynamics of n electrode and the instantaneous electrode potential Φ are displayed in Figureb. As the simulation progresses, the thermopotentiostat adjusts n electrode so that Φ approaches the prescribed target Φ_0_.

(a) Au(111)/H2O slab with the CDCE setup. (b) Time evolution of the electrode charge and the corresponding electrode potential. The target potential Φ0 – Φzc (horizontal line) is set to – 2.0 V. (c) Oxygen density profiles obtained from AIMD simulations of Au(111)/H2O under different applied voltages and at the potential of zero charge.

The applied potential has a pronounced effect on the interfacial water structure, as evident from the oxygen density profiles in Figurec. At negative potentials the oxygen atoms are repelled from the gold surface, shifting the first density peak away from the interface. Conversely, at positive potentials the surface charge attracts oxygen, pulling the two main peaks of the density profile closer to the metal. These structural changes are accompanied by variations in vibrational spectra, capacitance, and dynamical properties of the water layer, ?,? underscoring the importance of explicit potential control in simulations of electrochemical interfaces.

QM/MM Solvation Models

3.4

Beyond explicit solid–liquid simulations, the ability to modify the external potential V ext makes it straightforward to incorporate a variety of solvation schemes. In particular, one can replace the explicit solvent by an effective external potential V ^solvent^, enabling implicit solvent models or QM/MM (quantum-mechanics/molecular-mechanics) approaches.

Figure illustrates a QM/mean-field-MM implicit solvation workflow implemented using the VASP-Python interface. The system under study is a Mg^2+^ ion dissolved in water. First, a classical molecular-dynamics (MD) simulation of the aqueous environment is performed for several nanoseconds. During this MD run the charge density of the solvent is averaged, yielding the MM solvent charge density ρ_solvent_ ^MM^ shown in Figurea. This charge density is then converted into an electrostatic potential V solvent ^MM^, which is added to the Kohn–Sham Hamiltonian as the external potential V ext. More details about the QM/MM solvation method will be described in a forthcoming paper.

(a) Contour plot of the solvent charge density ρsolvent MM of a Mg2+ ion solvated in water. (b) Comparison between the electrostatic potential of a single Mg2+ in vacuum and in the solvated QM/MM method.

Figureb compares the resulting electrostatic potentials for two cases. In vacuum, the isolated Mg^2+^ ion exhibits a long-range, parabolic electrostatic potential (green line). When the solvent electrostatics are included (blue line), the +2 charge of the magnesium ion is screened by the surrounding water molecules, producing an oscillatory potential profile that reflects the formation of distinct solvation shells around the ion.

Because the solvent contribution is introduced solely through V ext, the QM/MM scheme can be readily extended to more complex geometries, such as surfaces and interfaces. Users can define arbitrary QM regions and couple them to a flexible MM description of the environment, making the approach highly adaptable to a wide range of electrochemical and interfacial problems.

Conclusion

4

In this work we have introduced a flexible scheme for imposing arbitrary electrostatic potentials in density functional theory (DFT) supercell calculations. The approach is implemented via the VASP–Python interface, which allows the user to supply any desired external potential V ext and to have it incorporated directly into the Kohn–Sham Hamiltonian. Because the external field also interacts with the ionic cores, we derived and applied the necessary energy- and force-correction terms to ensure that total energies and atomic forces remain physically meaningful.

The utility of the method is illustrated through several representative case studies:

• Surface adsorption under an applied bias, demonstrating how the adsorption energy can be tuned by the external field.
• Field-ion microscopy simulations, where the impact of external fields on atom diffusion and desorption is captured accurately.
• Electrochemical interfaces, showing the ability to model charged electrodes and to explore potential-dependent reaction pathways.
• Implicit-solvation and QM/MM setups, in which the solvent is represented by an effective potential V ^solvent^ (i.e., the averaged electrostatic potential generated by the MM charge density) that is added to V ext.

By exposing the external potential as a user-controlled input, the implementation provides full control over the electrostatic environment while retaining the standard workflow of VASP calculations. Consequently, it offers a versatile and easily extensible framework for investigating a wide range of field-induced phenomena  from catalysis and corrosion to nanoscale device physics within the well-established DFT paradigm.