Modelling electrochemical systems with finite field molecular dynamics

Physical chemistry of electric double layers and ionic solutions plays a fundamental role in energy related applications such as electrocatalysis, super-capacitors, fuel cells, lithium/sodium ion batteries. A realistic representation of these electrochemical systems requires treating electronic, structural and dynamic properties on an equal footing. Density functional theory based molecular dynamics (DFTMD) is perhaps the only approach that can provide a consistent atomistic description. However, one of the challenges in DFTMD modelling of electrochemical systems is the slow convergence of the polarization P, where P is the central quantity connecting to all electrical properties (the dielectric constant, the Helmholtz capacitance and the ionic conductivity). Here, we summarize our recent progress on developing finite field MD for computing electrical properties. Discussions on notable extensions and outlook for future works are also given.

. The schematic picture of EDLs at oxide-NaCl electrolyte interfaces. Bridging oxygen atoms of the oxide surface are protonated at one interface and adsorbed water molecules are deprotonated at the other interface. Compensating electrolyte Na + and Cl − ions would form two parallel EDLs eventually. The Helmholtz layer of EDL is highlighted by its molecular structure.
using density functional theory (DFT) methods. Here, the electrode is represented by metal clusters and slabs. Often, the double layer is merely treated by a water bilayer at the electrode surface [10] which, although useful as a first approximation, has been shown to be unstable with respect to thermal fluctuations. This suggests that a sufficiently large region of bulk needs to be included [11]. A number of implicit solvation models have been developed toward the needs of including explicit EDL descriptions at metal-electrolyte interfaces and reducing computational overheads [12][13][14][15]. In the meanwhile, recent progress on the full atomistic DFT modelling of metal-electrolyte interfaces can be found in reference [16][17][18].
In comparison with clean metal electrode interfaces [19,20], new complexities are added in the case of oxide electrodes [21,22]. The surface of oxides is known to undergo hydroxylation in the presence of adsorbed water molecules, thereby forming either acid-type or base-type adsorbed hydroxyl groups (figure 1). The condition in which the surface is proton charge neutral is the so-called point of zero proton charge (PZPC) [23][24][25]. When deprotonated at high pH, acidic sites acquire a negative charge which is compensated by hydrated cations in solution. Similarly, protonation of basic sites at low pH creates a positively charged surface with hydrated anions as the neutralizing counter charge. The surface charge density in these compact "protonic" EDLs at the aqueous interface of oxides can exceed the electronic charge density at metal/electrolyte interfaces by an order of magnitude. Formation of these protonic EDLs is a general phenomenon observed for all oxides, insulating main group oxides, such as silicates and aluminates as well as the conducting oxides formed by transition metal oxides. The study of protonic EDL is therefore an interdisciplinary area where colloid science and electrochemistry overlap, and the focus in of this review.

Ionic solutions
In ionic solutions, most solvents are made of polar molecules. Naturally, water occupies the leading position in the development of electrochemistry. Simple electrolytes which are completely dissociated in dilute solution were first used as prototype systems for the development of analytical theories such as the well-known Debye-Hückel theory [26]. This was followed by Kirkwood's theory to link the macroscopic dielectrics of polar liquids to the local orientation of molecules (i.e. the g-factor) [27] and the theoretical developments building on it by Deutch [28] and Wertheim [29].
Theories of ionic solutions were mostly done in dilute electrolyte solutions [30]. As the electrolyte concentration increases, ion-ion interactions become very important and determine the thermodynamic and dynamical properties of electrolyte solutions. At higher concentrations, several types of non-ideal phenomena show up, like ion-pairing [31] in which cations and anions associate and form metastable neutral pairs and even ionic aggregates. When the concentration further increases (e.g. 21 m LiTFSI in water), the system enters the regime of so-called 'water in salt' electrolyte (WiSE) [32]. This leads to remarkable electrochemical stability and opens up the potential window from 1.23 V to 3 V for lithium-ion battery applications. Moreover, the cross-correlated movements of different ions (of equal or opposite charge) can notably reduce the ionic conductivity and lead to deviation from the Nernst-Einstein relation [33]. Concentration effects can also alter the ion transport mechanism from the vehicular mechanism to the structural mechanism. In the former case, the ion travels with its solvation shell and diffuses slower than the solvent; in the latter case, the ion hops among the solvent molecules and diffuses faster than the solvent. In aqueous electrolyte solutions, the classic example for the vehicular mechanism is Li + diffusion and the counterpart for the structural mechanism is H + /OH − diffusion via the Grotthuss mechanism [34]. In non-aqueous organic electrolytes, the common ion transport mechanism is the vehicular mechanism. However, it has been shown that in highly concentrated sulfolane-based liquid electrolytes Li + conducts via hopping instead [35].
Experimental techniques such as x-ray adsorption spectroscopy, neutron scattering, vibrational spectroscopy, pulsed field gradient NMR and impedance spectroscopy provide valuable information regarding structural and transport properties of ionic solutions [36][37][38][39][40]. Complementary to experiments, simulation techniques-in particular molecular dynamics (MD) simulations-can be rather useful to provide insight and prediction regarding the solvation structure of ions and the mechanism of ion transport. Interested readers should look up the recent progress on the development of more accurate (reactive) force fields [41] as well as on analytical theories of liquids and electrolytes [42], which is out of the scope of this review.

Computing electrical properties from finite field molecular dynamics
When it comes to the modelling of EDLs and ionic solutions, a realistic representation requires including electronic, structural and dynamic properties. The density functional theory-based molecular dynamics (DFTMD) method is perhaps the only approach that can provide a consistent atomistic description.
Despite that, electrical properties, such as the dielectric constant, the Helmholtz capacitance and the ionic conductivity, were normally considered to be out of reach for DFTMD. One of the challenges therein is the slow convergence of the polarization P, where P is a central quantity which connects to all electrical properties (figure 2).
Our contribution in this area is to develop finite field MD simulation techniques for computing electrical properties [59][60][61][62][63][64][65][66][67]. The constant electric displacement D Hamiltonian, originally designed for treating spontaneous polarization in groundstate ferroelectric systems [76], is a new ensemble in statistical mechanics. We showed that the advantage of constant D simulations in computational electrochemistry is two-fold: a) it significantly speeds up simulations and makes the dielectric constant calculations of polar liquids converged on the tens of picoseconds timescale; b) it eliminates the finite-size error in computation of the Helmholtz capacitance due to periodic boundary conditions (PBCs). This methodology was further extended to study EDLs formed by polar surfaces and the electrolyte solution.
In the following, we will outline briefly the theoretical basis of finite field Hamiltonians and summarize our recent works in computing the dielectric constant of polar liquids and the Helmholtz capacitance at charged insulator-electrolyte interfaces.

Modern theory of polarization and itinerant polarization
The dipole moment in crystalline materials of infinite extension is not uniquely defined but depends on the choice of unit cell. According to the modern theory of polarization [69][70][71], the macroscopic polarization P is a multi-valued quantity and only changes of polarization are physically meaningful. Uniform polarization P in the modern theory of polarization, as well as in the Maxwell-Lorentz continuum theory [72,73], is not based on a multipole expansion of the charge density but is defined in terms of the internal current j: Accordingly, a change in polarization ∆P is related to the transient current density j(r, t) by where Ω is the cell volume. Equation (2) led to the Berry phase (γ) definition of the electronic polarization (of a cubic lattice with cell length L) [69][70][71].
where Ψ is the N-particle wave function and e i(2π/L) q i and r i are the charge and the position of particle i respectively. Finally, the time integral leads to Note that the value of polarization depends on the choice of t = 0 supercell while the change of polarization should not. In the literature of MD studies of electrolyte solutions [74][75][76], polarization corresponding to a time integral of current was referred to as itinerant polarization. P in equations (3) and (5) is the polarization generated by the (transient) flow of all charges in the system, including possible mobile ions. In the light of these, it is clear that the Berry phase polarization developed by the solid state physics community and the itinerant polarization known in the physical chemistry community are indeed the same.
Instead, the so-called cell polarization is defined as, where L stands for the cell matrix and nint(x) is a rounding function to the nearest-integer. The difference between itinerant polarization and cell polarization is that the former preserves the continuity. This means that the particles need to be tracked in this case when leaving the cell. In contrast, the particles would be wrapped into the cell at each time for calculating cell polarization. In Caillol's 1994 paper, he showed that that cell polarization as defined in equation (6) violates the Stillinger-Lovett sum rule for ionic solutions [76]. Instead, itinerant polarization satisfies all key conditions in statistical mechanics.

Thermodynamic and mechanical basis of finite field Hamiltonians
The central quantity behind the finite field method is the electric enthalpy functional [68]. The electric enthalpy of a system of volume Ω is written as:where ν = r N , p N are the classical degrees of freedom of the N atom system. The electronic degrees of freedom φ α , α = 1 . . . M, such as the orbital coefficients specifying the one-electron orbitals and determining the induced dipoles, are collectively represented by φ. P(ν, φ) is the macroscopic polarization density for the microscopic state specified by ν and φ. E is the homogeneous Maxwell field which needs be distinguished from the external field E 0 . The Hamiltonian associated with the electric enthalpy F (ν, E) is obtained as the variational minimum in the electronic degrees of freedom: is the value of φ for which the energy is minimal at a given atomic configuration ν and electric field E. The potential energy function in equation (8) can be compared to the Born-Oppenheimer surface of electronic structure calculation. The equivalent of the electric enthalpy of equation (8) at finite temperature is the free energy of the ensemble generated by the same Hamiltonian: Z E is the electric field dependent partition function where β = 1/k B T is the inverse temperature. Note that the combinatorial prefactor 1/(h 3 N N!) was omitted.
Taking the derivative gives the polarization equation (11) can be regarded as the electric equation of state for a uniform insulator.
In 2009, Stengel, Spaldin and Vanderbilt (SSV) had expanded the constant E method to the constant D method [68]. They introduced a new functional, the electric internal energy functional U (ν, φ, D).
Minimizing in the electronic degrees of freedom φ yields the effective Hamiltonian associated with the electric internal energy U (ν, D), i.e. the constant D counterpart of equation (8).
where φ 0 (ν, D) is the value of φ for which the energy is minimal at given atomic configuration ν and displacement field D.
The corresponding electric internal energy can also be obtained from the D dependent partition function Z D as with Note that U(D) is a (Helmholtz) free energy, the same as F(E). Evaluating the derivative with respect to the control variable D, we recover the macroscopic field The second identity follows from which defines the dielectric displacement in the Maxwell theory. When equation (16) is substituted in the Maxwell field expression for electrical work (see Landau and Lifshitz [77]) It is evident that U in equation (14) is indeed the electric free energy. This further justifies E as the macroscopic Maxwell field instead of the external field E 0 . If the electric enthalpy F(E) of equation (9) and internal energy U(D) of equation (14) are Legendre transforms of each other, then the E derivative of F in equation (11) should give −D. However, this is not the case. To fix this, SSV adjusted the definition of F(E) toF(E) following Landau and Lifshitz [77].
Then with equation (17) dF The term added −ΩE 2 /8π is a constant in the ensemble generated by the electric enthalpy Hamiltonian of equation (8) and will therefore not affect any ensemble average. Recently, one of the authors in this review showed that finite uniform E can be accounted for by adding it as a new degree of freedom in an extended Lagrangian [65]. Representing the uniform polarization P as the time integral of the internal current and E as the time derivative of a uniform vector field A, one can define an extended Lagrangian coupling A to the total current j t (internal plus external) and hence derive a Hamiltonian resembling the minimal coupling Hamiltonian of electrodynamics. By transforming the j t · A coupling to P · D, the resulting Hamiltonian is identical to the constant D Hamiltonian given in equation (13). This provides a clarification of the mechanical foundation of the finite field Hamiltonians.
Before closing this section, it should be noted that the polarization P in finite field Hamiltonians includes all charges including mobile ions in the system. Thus, the source of D is just the charge on the 'virtual' electrodes at infinity (the tin foil) acting as the electric boundary conditions for Ewald summation. Consequently, the constant D Hamiltonian corresponds to open circuit conditions where this charge is fixed and the constant E Hamiltonian corresponds to constant voltage conditions ( figure 3). This is a distinct difference with respect to the practice in the physical chemistry of ionic solutions where D is the electric field generated by the mobile ions treated as external charge.

Connection of finite field Hamiltonians to the Yeh-Berkowitz correction and the dipole correction
When it comes to modelling electrochemical interfaces, the hybrid constant electric displacement D Hamiltonian should be used instead. 'Hybrid' means the D field and the corresponding electric boundary condition are only applied in the direction of the surface normal. The hybrid constant D Hamiltonian is written as where P is in the same direction as D and P equals to M/Ω where M is the corresponding dipole moment. Writing equation (21) in terms of the dipole moment M we have Setting D = 0 in equation (22), we recognize the Yeh-Berkowitz dipole screened slab Ewald Hamiltonian [78] which is familiar to the physical chemistry community. Finite D adds the expected coupling between the dipole moment and the applied field (the −DM term) plus a vacuum energy of the applied field (the ΩD 2 /(8π) term). If we take the derivative of the hybrid constant D Hamiltonian over the component R I of nuclear position (in the same direction as D and P) to get the corresponding force F I : where M n (ν) is the polarization due to nuclear charges and M e is the electronic polarization. In the case of a point charge system, M e is zero.
because the electronic state φ 0 (ν, D) is a stationary point and the second term vanishes at the variational minimum. Similarly, = 0 Therefore, the force on nucleus I in the hybrid Hamiltonian turns out to be: where Q I is the nuclear charge of atom I in electronic structure calculations and the partial charge of atom I in classical force fields. The corresponding field generated from the dipole moment E M at D = 0 is: Then we recognize that this is the same expression as the dipole correction known in the solid state community [79,80]. In other words, the hybrid constant D Hamiltonian at D = 0, the Yeh-Berkowitz correction and the dipole correction are simply equivalent.
What is the difference then? In the authors' view, the main difference is that both the Yeh-Berkowitz correction and the dipole correction originate purely from the electrostatics and consider the inclusion of a vacuum slab as a prerequisite for the fix. Instead, finite field Hamiltonians, as elaborated in this section, have a firm thermodynamic foundation and subsequently lift the conceptual constraint that a vacuum slab must be the integral part of the supercell modelling. More importantly, they allow to define, control (constant E), calculate (constant D) a potential difference over a periodic cell seemingly in violation of the PBCs [76,81]. Indeed, these conceptual breakthroughs led to a number of new developments in the modelling of electrochemical systems with finite field MD as summarized below.

Condition
Formula

Four basic formulas
Following Kirkwood-Fröhlich theory [27,82], it should be possible to obtain an estimate of the dielectric constant from polarization fluctuations either under E = 0 or D = 0 dynamics. These results should in principle agree with finite field simulations with sufficiently small E or D fields, as shown before [59]. A summary of four basic formulas for computing dielectric constant of polar liquids is given in table 1 [59,61,[83][84][85]. Further, we found that the ratio of polarization fluctuations at E = 0 and D = 0 leads to [59]: For polarizable polar liquids, i.e. the ones in electronic structure calculations and polarizable force fields, the optical dielectric constant ϵ ∞ can be computed using either of the following two equations.
where P is the total polarization at fixed geometry ν in the direction where the field is applied. ϵ ∞ = 1 for non-polarizable polar liquids. There are two important aspects regarding these four basic formulas given in table 1. First, switching from E = 0 to D = 0 not only leads to the suppression of polarization fluctuations (figure 4a) but also speeds up the relaxation time [59]. The relaxation of polarization in the classical Debye theory is exponential (as seen in figure 4b). The difference between relaxation times at D = 0 and E = 0 may be interpreted as a longitudinal versus transverse (or Debye) relaxation time τ L = τ D /ϵ [87,88]. Second, these formulas differ in practice because of different requirements on the accuracy. At least several nanoseconds are needed to converge computation of the dielectric constant of liquid water at ambient conditions, in the case of E = 0.
Despite a much faster relaxation time at D = 0, the troublesome inverse relation between fluctuations and the dielectric constant requires the accuracy in the second moment to be proportionally higher. Thus, finite E and D versions may be better choices. This also motivated us to explore an alternative solution to this long-standing problem, as summarized in the next section.

Finite field methods and the Kirkwood g-factor
In Kirkwood's 1939 paper, he postulated that it should be possible to express the dielectric constant ε of polar liquids in terms of a short-range orientational correlation function [27]. Building on Onsager's local field approach [89], Kirkwood gives his famous formula as N is the number of polar molecules in a system of volume Ω. µ is the value of the dipole of a molecule in the polar liquid. The orientational correlations are contained in a single number, the Kirkwood g-factor g K . Setting g K = 1 in equation (31) recovers Onsager's mean field approximation [89]. For correlated polar liquids g K is obtained as the asymptotic value of the r-dependent Kirkwood g-factor where µ 1 is the dipole of a reference molecule 1 at the center of a sphere of radius r. M 1 (r) is the sum of total dipoles µ i in that sphere. Because local orientational correlations are short-ranged, Kirkwood argued that G K (r) should approach a constant g K beyond a certain distance r K . This distance is the Kirkwood correlation length and g K = G K (r)| r≥rK [27]. r-dependence of Kirkwood g-factor depends on electric boundary conditions. This was analyzed in detail by Caillol [90] using an embedding dielectric medium model of dielectric constant ε ′ . By identifying ϵ′ = ∞ with E = 0 and ε ′ = 0 with D = 0, we showed in our work [61] that a superposition of the correlation function (equation (32)) evaluated under E = 0 and D = 0 conditions leads to the desired short-ranged Kirkwood g-factor G Kc (r) without the need for prior knowledge of ε.
This formula is general and can be applied to both polarizable and non-polarizable polar liquids. The more interesting feature of the short-ranged Kirkwood g-factor G Kc (r) is that its convergence depends on r. The dielectric constant ε(r) computed from G Kc (r) at larger r values shows the poor convergence familiar from the computation of ε from fluctuations of the volume dipole (table 1). In contrast, for distances smaller than the Kirkwood correlation length, ε(r) is to a good approximation converged on the tens of picoseconds timescale, which made the computation of the dielectric constant of liquid water possible using DFTMD simulations [61].
The formula given in equation (33) requires two separated simulations which may not be so convenient in practice. Since the corresponding linear combination of E = 0 and D = 0 fluctuations is a 2/3 and 1/3, one can in practice obtain the same result using a hybrid constant D Hamiltonian equation (21) (with D = 0) for the isotropic polar liquids instead. In other words, the short-ranged Kirkwood g-factor G Kc (r) can also be generated using the following formula: The reference value of ε for bulk SPC/E water at ambient conditions is about 72, see reference [59]. Adapted from reference [61] with CC-BY licence.
Note that D = 0 means it is the hybrid constant D Hamiltonian, e.g. E x = 0, E y = 0 and D z = 0. The comparison between G Kc (r) obtained from the superposition of E = 0 and D = 0 simulations and G K (r) D=0 obtained directly from the hybrid D = 0 simulation is shown in figure 5.

The Helmholtz capacitance at charged insulator-electrolyte interfaces
The differential capacitance C EDL at oxide-electrolyte interfaces is a key probe of the structure of EDLs.
Generally speaking, C EDL is made of three terms [91,92]: The first term C SC is originated from the electron accumulation or depletion in the space charge (SC) layer, which may be up to 100 nm thick. The second term C H is the Helmholtz capacitance because of the acid-base reactions in response to the pH. The last term C GC is the so-called Gouy-Chapman capacitance, due to the distribution of dissolved electrolyte.
In equation (35), C SC will normally be the leading term because it is the smallest among these three capacitances. However, its contribution has been eliminated at the flatband potential [91]. For photoelectrocatalysis applications happening at the high ionic strength, the capacitance of the diffuse layer becomes higher, which makes the inverse C GC term negligible. Based on these considerations, computing the Helmholtz capacitance C H at both non-polar and polar surfaces is what we focus on.

Computing the Helmholtz capacitance of non-polar surfaces from the supercell polarization
Under PBCs, proton charges at non-polar surfaces can be set up either in a symmetric way or in an asymmetric way. In the former case, two sides of insulator surfaces take up an equal amount and the same kind of surface charges [93][94][95][96][97][98][99]; in the later case, they take up an equal amount but opposite kinds of surface charges [100,101]. The difference is that the composition of the system does not change when charging up in the asymmetric setup, which resembles the experimental setting and makes it convenient for comparing total energies. Therefore, this is the setup we used in our works.
In this setup as shown in figure 6(a), the cell is made of two parallel EDLs which lead to a net polarization. Since the overall potential difference over the cell has to be zero because of PBCs, there must be a finite electric field in the insulator (e.g. a vacuum slab in figure 6(a)). A finite field means the EDL is not fully compensated and bears a net charge Q net according to Gauss's theorem, where Q net is the sum of the proton charge density on the non-polar surfaces and the polarization charge density from solvent molecules and ions. Moreover, this net charge in EDLs is the sign of a finite-size error in the computation of C H .
Following Maxwell's equation, the electric field E z (z) and the laterally averaged charge density ρ(z) along the z axis are related by Since E z (z e ) = 0 in the bulk electrolyte (z = z e ), the net EDL charge is then given as where A is the surface area and E d = E z (−L/2) is the field in the dielectric (oxide slab). For our simple system shown in figure 6a, E d is constant and larger than zero, it means Q net < 0. This further suggests that EDL is fully compensated only in the limit of E d = 0.
To derive the electric equations of state and the conditions satisfying Q net = 0, we used a minimal Stern model shown in figure 6b as the continuum counterpart of figure 6a. As a first approximation, we assume these two EDLs equivalent. Each EDL is characterized by the dielectric constant ϵ H , the thickness l H , and the field E H . The compensating charge density σ from the electrolyte solution is not necessarily the same as the imposed surface charge density σ 0 . In addition, the electrolyte as an ionic conductor is characterized by the thickness l e , the dielectric constant ϵ e = ∞, and the field E e = 0. Finally, the corresponding parameters for the insulator slab are the thickness l d , the dielectric constant ϵ d and the field E d . Note that the supercell size To show the equivalence of the atomistic system and the continuum model, it is instructive to re-derive the expression of Q net from the surface theorem [102]. The total charge σ d on the insulator plane (see figure 6b) is Table 2. Three formulas for computing the Helmholtz capacitance CH at non-polar insulator electrolyte interfaces.

Condition
Formula where P H and P d are polarizations in the Helmholtz layer and in the insulator slab respectively. σ d is related to E H and E d as Note that E H and E d are in opposite directions as shown in figure 6b, so do P H and P d . The total counter charge σ H is Because E e = 0 in the electrolyte, σ H can be expressed as Finally, combining equations (38) and (40), this leads to the expression of Q net as which according to equations (39) and (41) is equal to Therefore, both the atomistic system and the continuum model lead to the same expression of the surface net charge.
Through the analysis of the electric equations of state of this Stern-type model, we have devised two approaches to compute the size-independent C H of charged non-polar surfaces under PBCs [60]. The first approach is using constant E simulations, where E field is in the same direction as the surface normal. By recovering the zero net charge (ZNC) state of EDL, the corresponding E field directly gives C H [60]. The second approach is using constant D simulations. The derivative of P with respect to σ 0 provides an efficient way to estimate C H [60].
In the second approach, C H was obtained without the need to locate the ZNC state [60]. This means the corresponding formula can be derived alternatively instead restoring to the Stern-type model. Therefore, through thermodynamic relations, we came up with the third approach for computing C H using P at E = 0 and the composite dielectric constant ϵ ⊥ of the whole system in the direction of the surface normal [64]. A summary of these three methods is given in table 2. Note that L in the denominator does not mean that C H vanishes when L → ∞, because ∆V = −LE ZNC is the potential over the supercell which converges to the sum of the potentials over the double layers. The same applies to ∆V = 4πL⟨∆P⟩ at the constant D.
One important aspect of applying the finite D formula in table 2 is to realize that the electric displacement inherits the multi-valued nature of polarization [60].
Polarization in two branches are related according to where n is an integer and e/A is the quantum of polarization [71]. As mentioned before, itinerant polarization needs to follow the same branch staring from t = 0. Since the physical state of the system must be the same regardless of the choice of the branch, the dynamics of the atoms should not differ even if we change branches. Hence, when the polarization is changing the branch, so must the electric displacement in order to conserve the electric field In addition, one should be aware that finite field methods may fail if a long-lived resonance state in the presence of a field is not a well-defined solution in the context of electronic structure calculations. This is known as the Zener breakdown and suggests that the limiting field strength E ≪ E g /L where E g is the band gap of the system and L is the box length in the direction of the field [70].

The case of polar surfaces
The terminations of ionic crystals can leave the solid with a surface carrying a net charge. According to the classification of Tasker [103], Type-II terminations acquire a unit cell with a zero dipole (e.g. − -+2 -−) while type-III terminations lead to a unit cell with net dipole (e.g. + -). The latter case is a polar termination which generates an internal electric field and leads to an unsustainable potential difference across the slab of the increasing width. The familiar example of Type-II termination is the (111) surface of CaF 2 (fluorite) and that of Type-III termination is the (111) surface of NaCl (rocksalt). Other examples of Type-III surfaces include (0001) surfaces of the corundum (Al 2 O 3 , Fe 2 O 3 ) and wurtzite (ZnO) structure [104]. Polar surfaces can be stabilized by a surface reconstruction which eliminates the bulk dipole moment. Such a reconstruction is necessarily non-stoichiometric as it must change the net surface charge [104]. For polar interfaces with vacuum, the non-stoichiometric reconstruction is often observed to occur by removal or addition of ions. However, it is also possible for a transfer of electrons to provide the required surface charge, which is termed as an electronic reconstruction. For polar surfaces in contact with an electrolyte, the compensating charge may instead be provided by an exchange of ions with the electrolyte, and this leads to the formation of an EDL at the polar surface.
To make the simulation of the semi-infinite system possible, we explored a Stern-type continuum model similar to the one outlined in the previous section. Once again, we came up with two methods: one using a constant E field, and the other using a constant D field. These are summarized below and corresponding formula for computing C H at the polar surface are given in table 3.
The first method introduces a compensating E field in order to obtain the state of compensating net charge (CNC) [62]. This CNC condition, E = E CNC , was empirically determined by searching for the E which eliminates the potential drop over the polar crystal slab, as required to isolate the EDL capacitance. It is interesting to note that the compensating charge provided by the electrolyte under the CNC condition is σ CNC = (n + 1)σ 0 /2n. In the limit n → ∞, the prefactor tends toward 1/2, which agrees with the famous Tasker rule for rocksalt (111) polar surfaces [103,105].
The second method is to simulate CNC by keeping the electric displacement D fixed instead. By deriving the expression of D CNC and validating it with both force field based MD and DFTMD simulations [66], it was found that D CNC involves only structural parameters and can therefore be known a priori. This saves considerable computational time. Care must be taken, however, as the exact expression for D CNC depends on the gauge. One can, for instance, switch between an insulator centred supercell (ICS) and an electrolyte centered supercell (ECS). Unsurprisingly, the ECS gauge involves placing the electrolyte in the middle, and the solid on the sides. However, depending on whether the supercell boundary is located in an interplanar layer with electric field E 1 or in a layer with electric field E 2 , one may further distinguish between an 'ECS1' Table 3. Four formulas for computing the Helmholtz capacitance CH at polar ionic solid electrolyte interfaces.

Condition
Formula and 'ECS2' . An example of an ECS1 gauge is shown in figure 7. What is important to note is that, given proper accounting, the final C H does not depend on the choice of gauge. Indeed this must be the case, as the capacitance is a physical observable.

DFTMD application to the TiO 2 rutile (110) surface
As the first application of finite field DFTMD to oxide systems, we applied the hybrid constant D simulation to the electrified rutile TiO 2 (110)-NaCl electrolyte interface [67]. The electronic structure of rutile TiO 2 (110)-NaCl electrolyte systems was solved applying DFT in the Perdew-Burke-Ernzerhof (PBE) approximation [106] using the CP2K suit of programs [107,108]. Double-ζ basis sets with one additional polarization functions (DZVP) optimized for molecular systems (MOLOPT) [109] and a charge density cutoff of 320 Ry were used. Core electrons were taken into account using the dual-space Goedecker-Teter-Hutter (GTH) pseudopotentials [110].
The model system consisted of a symmetric periodic slab of five O-Ti-O trilayers with a (4x2) surface cell, 110 water molecules, 5 Na + and 5 Cl − ions in a periodic cell of 11.90 Å × 13.20 Å × 38.28 Å. To keep the composition fixed in all setups, the negatively charged side of the TiO 2 slab is formed by removing H + from absorbed water molecules and the positively charged side is formed by adding removed H + to oxygen sites of the other side of TiO 2 slab. The integration time-step is 0.5 fs and MD trajectories were collected for about 15ps for each step-up at D = 0 condition after initial equilibration at E = 0 condition. The initial configuration for different charge densities was obtained from our previous work of charged insulator-electrolyte systems using a classical point-charge like model [60]. The Bussi-Donadio-Parrinello thermostat [111] was used to keep the temperature at 330 K throughout all simulations.
The macroscopic polarization is computed using both the Resta formula [112] and the maximally-localized Wannier functions (MLWFs) [113], as implemented in CP2K. This is an important technical check to ensure the branch consistency because of the multi-valued nature of the macroscopic polarization and the electric displacement as noted in the previous section.
Since P z is proportional to the surface charge density σ 0 , P z should be zero at the PZPC. This provides a reference for testing the convergence of DFTMD simulations in computing C H . As shown in figure 8, M z relaxes to zero almost instantaneously when switching from E = 0 to D = 0. Within 10 picoseconds which is the typical timescale in DFTMD simulations, the time average of P z reaches the expected value as 0.02 D.
Following the formula at constant D shown in table 2, one obtains C H which is about 76 µF/cm 2 . Despite of very different starting configurations at surface charge q = 2e and q = 4e, the resulting C H are in excellent agreement with each other (81 µF/cm 2 vs. 72 µF/cm 2 ). Note that experimental estimated capacitances for the same system are quite scattered, ranging from 64 to 160 µF/cm 2 [114][115][116].
From this decomposition of the overall Helmholtz capacitance into protonic C + H and deprotonic C − H using the macroscopic averaging technique [117], it is found that C − H is about 50% higher than C + H (101 µF/cm 2 vs. 67 µF/cm 2 at q = 2e; 67 µF/cm 2 vs. 59 µF/cm 2 at q = 4e). The suggests a higher Helmholtz capacitance at high pH, which is in agreement with titration experiments [116,118,119] and has also been seen in ZnO [120]     6. Notable extensions and outlook for future works 6.1. Modelling metal-electrolyte interfaces EDL at metal-electrolyte interfaces is a long-standing topic in electrochemistry. Instead of taking protonic charge due to the acid-base chemistry, charge at metal surfaces is originated from the rearrangement of electronic density in response to the external potential or the electric field.
One widely used method for constant potential simulation at metal-electrolyte interfaces is the Siepmann-Sprik method [122], which was originally developed for mimicking scanning tunneling microscope on probing the adsorption of water molecules at a metal surface and adapted later for modelling electrochemical cells [123]. In the model, Gaussian charge distributions q i η 3 π 3/2 exp [−η 2 (r − r i ) 2 ] were put on each atom of the electrodes and the actual value q i at each MD step is determined by requiring that the potential is strictly constant in an electrode with the difference between electrodes matching the imposed cell potential ∆Φ.
As mentioned previously, the distinction of D in finite field Hamiltonians from the one used in physical chemistry of ionic solutions is that the source of D is just the charge on the 'virtual' electrodes at infinity. This implies that ∆Φ =−EL can be directly interpreted as the voltage over the orthorombic supercell in which E field is along a side with length L. Indeed, as shown recently by Salanne and co-workers [121], finite field simulations with 3D PBC give the same result as the constant potential ∆Φ simulation with 2D PBC (figure 10). This proof-of-concept study opens the door to an efficient simulation of metal-electrolyte electrochemical interfaces. Nevertheless, its DFTMD extension remains to be seen.

Modelling electrolyte solutions
Besides electrified solid-electrolyte interfaces and polar liquids, finite field Hamiltonians have been explored for modelling electrolyte solutions as well. It has been shown that the constant E Hamiltonian allows to derive the linear response relation for the ionic conductivity in a much straightforward manner by avoiding Hamiltonians depending on current densities and vector potentials [124]. More importantly, it is found that applying the constant D Hamiltonian to the electrolyte system satisfies the Stillinger-Lovett condition, i.e. ⟨E⟩ = 0. In other words, 4πP = D for bulk electrolyte solutions.
As seen in figure 11, 4π⟨P x,ion ⟩ = D x and 4π⟨P x,wat ⟩ = 0 [124]. Since ⟨E⟩ = 0 holds only for conductors and purely (non-dissociable) water can not make a conductor without ions, this result, that the response of polarization to the D field comes entirely from the ions, may not come as a total surprise .
In addition, it is worth to mention that simulations of the electrolyte solutions with the constant D Hamiltonian need to be done by coupling a D field to the polarizations of both ions and solvent molecules. Otherwise [125], this likely violates not only the Stillinger-Lovett condition but also disobeys the definition of P in finite field Hamiltonians which should include both solvent molecules' and mobile ions' contributions.

Electromechanical coupling at interfaces
Coupling of electric field and stress at solid-liquid interfaces is a rather interesting but difficult topic. The Maxwell stress tensor σ M is a second rank tensor quantifying the stress generated by long-range electrostatic forces. Adding the stress tensor σ due to all other supposedly short-range forces gives the total stress tensor σ t = σ + σ M controlling the mechanical equilibrium. The difficulty has already arisen at the theoretical level not mentioning the simulations, because there is no unique way of partitioning the stress into electrostatic and short-range contributions.
Recently, an attempt has been made to connect two collective properties-total surface tension σ t and a finite dipolar surface potential-at a liquid water vapour interface using finite field MD [126]. An upside down parabola with a maximum shifted away from zero field was found when determining the response of the tangential component of the surface tension to the application of an electric field normal to the surface using finite field MD simulations ( figure 12). This leads to the first approximation of the electromechanical surface potential χ σ 0 according to equation (46). Note that χ σ 0 is a related but different quantity from the electrochemical surface potential χ S and the surface dipole potential χ d .
where γ T is the tangential component of the surface tension as evaluated in the MD simulation. z v and z l are positions in the vapour phase and in the bulk region (defining as p N = p T ) of the liquid .

Finite field MD simulations with machine learning potential
Simulations of electrical properties necessitate long timescales that are normally beyond the reach of standard DFTMD. One way to tackle this timescale challenge is to explore finite field MD simulations to speed up the convergence of the polarization P as discussed in this review. The other way to solve this problem is to make use of (reactive) force fields to access longer timescales. One promising approach in this direction is to devise Figure 12. Field-dependent tangential component γ T of the surface tension of the liquid water vapour interface with the SPC/E model [86]. Lz = 5.55 nm, 8.33 nm and 11.11 nm are the values for the length in the direction normal to the surface for a small, medium and large MD box respectively. A dividing surface is chosen to extract the value of γ T for the interface, which is about 0.9 nm from the onset of the tangential pressure profile p T (z). Reproduced from reference [126] with permission from the PCCP Owner Societies. high-dimensional neural network potentials (NNPs) with DFT quality as proposed by Behler and Parrinello [127], which enables computing ionic conductivities in alkaline electrolyte solutions [128].
Recently, one of the authors and co-workers have taken the initiative and developed an open-source Python library named PiNN (https://github.com/Teoroo-CMC/PiNN/), allowing researchers to easily develop and train state-of-the-art atomic neural network architectures specifically for making chemical predictions and approximating the potential energy surface. In particular, they have designed and implemented an interpretable and high-performing graph convolutional neural network (GCNN) architecture PiNet [129], and demonstrate how the chemical insight "learned" by such a network can be extracted.
What comes next is to include predictions of electronic properties such as the polarization as well as the coupling to the Maxwell field in periodic systems ( figure 13). This will allow us to fully explore the potential of finite field MD simulations for modelling electrochemical systems.