Dielectric properties of nano-confined water: a canonical thermopotentiostat approach

We introduce a novel approach to sample the canonical ensemble at constant temperature and applied electric potential. Our approach can be straightforwardly implemented into any densityfunctional theory code. Using thermopotentiostat molecular dynamics simulations allows us to compute the dielectric constant of nano-confined water without any assumptions for the dielectric volume. Compared to the commonly used approach of calculating dielectric properties from polarization fluctuations, our thermopotentiostat technique reduces the required computational time by two orders of magnitude.

We introduce a novel approach to sample the canonical ensemble at constant temperature and applied electric potential. Our approach can be straightforwardly implemented into any densityfunctional theory code. Using thermopotentiostat molecular dynamics simulations allows us to compute the dielectric constant of nano-confined water without any assumptions for the dielectric volume. Compared to the commonly used approach of calculating dielectric properties from polarization fluctuations, our thermopotentiostat technique reduces the required computational time by two orders of magnitude.
Molecular dynamics (MD) has become an indispensable tool to efficiently simulate the behaviour of a wide range of systems. Experiments, however, are often performed using basic thermodynamic variables that are different from the ones easily accessible in simulations. The need for constant temperature as opposed to the much simpler constant energy simulations is widely recognized as one of the most important examples of this kind. Consequently, significant effort has been directed at developing thermostats [1][2][3][4][5][6][7][8][9][10][11] with the dual purpose of (i) efficiently sampling the canonical ensemble and (ii) enabling direct control of the temperature. With the advent of robust techniques to apply electric fields in densityfunctional calculations [12][13][14][15][16][17][18][19][20][21][22][23][24][25], there has been continuous interest to use MD simulations to study electrically triggered processes involving electron transfer reactions, such as electrochemical reactions, field desorption, and quantum transport. For this purpose, it is necessary to conceive a potentiostat, in analogy to the theory of thermostats, in order to incorporate the electric potential as a thermodynamic degree of freedom.
The first approach of this kind was pioneered by Bonnet et al. [17]. They chose a computational setup that is grand-canonical with respect to the electronic charge, cf. Fig. 1a. The system is then coupled to an external potentiostat, where the charge is treated as an extended spatial coordinate. Analogous to how the mechanical momenta are obtained from the derivatives of the Hamiltonian with respect to the spatial coordinates, a fictitious momentum of the charge is obtained from the derivative of the total energy with respect to the charge. This fictitious momentum is then coupled to standard Nose-Hoover dynamics. The grand-canonical nature of their setup led Bonnet et al. to recognize that "[...] in its present formulation, a requirement for implementing the potentiostat scheme is the existence of an energy function E(r i , n e ) that is differentiable with respect to the total electronic charge. This implies the ability to treat noninteger numbers of electrons and, in general, non-neutral systems." [17]. Unfortunately, in the context of densityfunctional calculations the total energy as a function of * wippermann@mpie.de the number of electrons is a notoriously difficult quantity to compute. Furthermore, the electronic charge is a single degree of freedom. Yet, thermostating single degrees of freedom by the Nose-Hoover method requires a chain of Nose-Hoover thermostats [6]. Both requirements are serious obstacles in implementing this potentiostat concept in existing density-functional theory (DFT) codes. In fact, none of the commonly used DFT codes [26] has a potentiostat scheme that provides the opportunity to study electrochemical systems with molecular dynamics and explicit water.
In the present study we show that the origin of these difficulties lies in the formally equivalent construction of the potentiostat to a thermostat: existing approaches consider a transfer of energy (thermostat) or charge (potentiostat) from the DFT cell to an external reservoir. While such a grand-canonical coupling can be straightforwardly implemented for an energy exchange this task is severely harder for a charge transfer since non-neutral DFT systems have to be considered. Moreover, there is also a fundamental difference between thermostats and potentiostats. For a thermostat, the system is considered to be embedded in an external temperature bath. In the physical realization of such a setup, heat transfer occurs only at the boundaries between the system and the temperature bath. Therefore, in simulations, it is common to model the simulation cell, the thermostat, and the energy exchange between them as separate entities. For a potentiostat, however, the energy exchange with the voltage source is mediated by the electric field. Since the electric field permeates the system, energy exchange occurs throughout the whole system and not just at the boundaries.
Thus, the electric field plays the same conceptual role for a potentiostat as the temperature bath for a thermostat. Yet, the electric field is an integral part of the system. In contrast to a thermostat where the bath is external, only the control mechanism of a potentiostat is external, but not the bath. We therefore propose to construct a potentiostat using the electric field as the control parameter. As will be derived in the following, this allows us to remove the need for treating charged systems. The actual implementation requires only quantities that are readily accessible in standard DFT codes, which makes it easy to integrate this potentiostat into existing electronic structure codes.
We start our derivation from the Hamiltonian given by Bonnet et al. [17] for the grand-canonical potentiostat scheme sketched in Fig. 1a: The explicit variables to describe atoms exposed to an external electric field are the atoms' positions q i , their momenta p i and the kinetic energy of the system K(p i ).
The single electrode present in the grand-canonical setup carries a charge n and is included in the simulation cell. A direct consequence of having a single charged electrode is that the system is not charge neutral but has an excess charge n ex = n. The potential energy V (q i , n ex ), which describes the interatomic interaction, thus explicitly depends on the (non-zero) excess charge. We note that the potential energy V (q i , n ex = 0) is simply the Born-Oppenheimer surface as computed in any DFT code. The term −nΦ 0 accounts for the potential of the external charge reservoir. To derive an equation of motion for the potentiostat from this Hamiltonian requires an explicit expression for the derivative of the charged system's energy V (q i , n ex ) with respect to the excess charge n ex [17].
A key concept of our proposed approach is to remove the need to compute V (q i , n ex ) explicitly. As an intermediate step to achieve this aim, we consider the system shown in Fig. 1b, which is based on concepts of the modern theory of polarization [12]. Here, the simulation cell contains the dielectric displacement field D created by moving a charge n from the left to the right boundary. The charge itself is outside the simulation cell. The total of the left and the right charge is zero, i.e., the inner region of the cell remains charge neutral. Without loss of generality, V can then be partitioned into the regular interatomic potential energy at zero excess charge V (q i , n ex = 0) and the electric energy E(q i , n) due to the presence of the D-field. The electric energy E is given by [14], where P is the polarization density. Thus, the difficult task of calculating V (q i , n ex = n) for the charged system is substituted by the much easier calculation of V (q i , n ex = 0) + E(q i , n) for a corresponding neutral system.
In a final step we extend this idea to a dielectric placed between two electrodes. The electrodes are connected to a voltage source with potential difference Φ 0 and an internal resistance R, cf. Fig. 1c. The corresponding Hamiltonian is [42]: Here, we explicitly include the two electrodes with charge n and −n that create the displacement field D. C 0 is the bare capacitance of the electrodes in vacuum without the dielectric and n p = Φ/C 0 − n is the polarization bound charge at the left electrode surface due to the polarization of the dielectric [27], where Φ is the instantaneous voltage Various computational setups used to include the effect of an applied electric field. a) The simulation cell contains a single charged electrode. It is grand-canonical and has an excess charge nex = n [15,17]. b) In the modern theory of polarization (MTP), either the electric field E or the displacement field D is used as the basic variable [14]. The simulation system is always charge-neutral. c) The simulation cell contains two electrodes with charge n and −n. In analogy to the MTP, it is canonical and charge-neutral. measured directly across the electrodes. Note that the bound charge n p is an implicit function of q i and n. The Hamiltonian of the canonical potentiostat (Eq. 2) allows us to obtain ∂H/∂n analytically: The reason why we are now able to obtain an analytical expression is that in the canonical system, n is no longer the net charge n ex of the system exchanged with an external bath as in the grand-canonical potentiostat approach. Rather, it is based on a charge transfer from one electrode to the other which leaves the total charge of the system unchanged. Hence, ∂H/∂n is determined by the infinitesimal amount of energy ∂H required to transfer an infinitesimal amount of charge ∂n between the electrodes. Thereby, Eq. 3 connects the derivative ∂H/∂n, that has to be computed for the microscopic quantum-mechanical system, to an external macroscopic quantity: the instantaneous voltage Φ. Since our total system is charge neutral, Φ is determined directly by the total dipole moment of the charge distribution within our simulation cell. This quantity is readily accessible in DFT codes. In the following, we will use these properties of Eq. 3 to construct an equation of motion for the voltage Φ and, subsequently, for the electrode charge n, the control parameter of our potentiostat.
A thermostated system evolves under the influence of its internal, energy-conserving Hamiltonian and extra forces that drive the exchange of energy with the thermal bath. Similarly, a potentiostat scheme requires a force-like term that drives the necessary changes of the electrode charge to keep the average potential constant.
To obtain this force we recast Eq. 3 in differential form: The first term in Eq. 4 expresses how the potential Φ will evolve under the system's internal, energy-conserving Hamiltonian. The second term, f dt, is the extra forcelike term that couples to the electrode charge. This extra-force term controls on the one hand the constancy of the (average) potential. On the other hand, it has to mimic the statistical fluctuations that occur in a finite system. These two aspects are balanced by the fluctuation-dissipation theorem [28] and ensure that the system stays in the NTΦ ensemble.
Changing the potential of a capacitive system is always a dissipative process, i.e., only adapting the electrode charge to control the voltage would drain energy from the system. To avoid this energy drain, any dissipation must be accompanied by a corresponding fluctuation that returns the removed energy. In other words, the applied electric field itself must have a finite temperature: the energy dissipated by the potential control mechanism must equal on average the energy gained from thermal potential fluctuations. For the electrical circuit shown in Fig. 1c, Johnson [29] and Nyquist [30] derived already in 1928 the relation between fluctuation and dissipation. This relation is now known as a specific case of the fluctuation-dissipation theorem (FDT), and determines the variance of the potential fluctuations as well as its distribution in frequency space. Using Ohm's Law and Kirchhoff's 2nd Law, the current through the circuit shown in Fig. 1c In conjunction with the FDT, we can then express the potentiostat force f dt as a stochastic differential equation (SDE) [42]: τ Φ := (RC 0 ) −1 and dW t are the relaxation time constant of the potentiostat and a stochastic noise term given by a Wiener process, respectively. The deterministic dissipation term in Eq. 6 is equivalent to Ohm's Law. The stochastic term provides the thermal fluctuations. Together, they satisfy the FDT exactly, also in frequency space [42]. Far from equilibrium, the deterministic part in Eq. 6 dominates and drives the system towards the target potential with a relaxation time of τ Φ . Close to the target potential, the deterministic term becomes small and the stochastic term takes over to ensure that in thermodynamic equilibrium the canonical ensemble is correctly sampled. The SDE Eq. 6 can be solved by employing Itō calculus [31]. Using this calculus and Eq. 5, Eq. 6 can be integrated analytically to a finite time step [42] where N is a Gaussian random number with N = 0 and N 2 = 1. Eq. 7 is the central result of our derivation: All quantities entering the right side can be either directly obtained from the DFT calculation or are known from the specific setup. Thus, the electrode charge n can be directly computed at each MD time step. This central result allows us to include the potentiostating process in any standard discrete time step MD scheme since the only quantity that needs to be extracted from the MD is the potential Φ(t). The scheme can be used equally well to perform ab initio or empirical potential MD.
To validate our canonical thermopotentiostat scheme and to demonstrate its performance we consider a topic that had recently gained a lot of attention. Recent experiments [32] showed that a water film confined to a few nm thickness changes its dielectric behaviour from the bulk dielectric constant of 80 down to 2. Thus, the presence of solid-water interfaces appears to modify the dielectric response of water from a highly polarizable medium, which is considered to be the origin of the unique solvation behaviour of water, down to a response that is close to the vacuum dielectric constant. Understanding and being able to qualitatively describe this mechanism is crucial since interfacial water is omnipresent.
Due to the relevance of this question in fields as diverse as electrochemistry, corrosion and electrocatalysis, several computational studies addressed dielectric properties of nano-confined water [33][34][35]. These studies used either the variance of the total dipole moment fluctuations per volume, using Kirkwood-Fröhlich theory [36], or the theory of polarization fluctuations [37]. They require as an input quantity the dielectric volume to compute the dielectric constant. However, the exact location of the boundary between the electrode and the dielectric is ill-defined in the presence of adsorbates, thermal motion of the electrode surface or in the context of explicit electronic structure calculations. For this reason, past studies often reported only dipole fluctuations perpendicular to the electrode surface, but not the dielectric constant ⊥ itself [33].
Our thermopotentiostat MD allows us to address this issue directly, since our setup shown in Fig. 1c exactly reproduces the experimental situation. In analogy to experiment, we compute the static dielectric constant as Before studying in detail the thickness dependent properties of nano-confined water let us use this system as a model to test our thermopotentiostat scheme. We therefore perform classical MD simulations [38] of liquid TIP3P water confined between two parallel electrodes. Numerical details are given in the supplementary material [42]. We first investigated the effect of our potential control scheme on the temperature of an NVE ensemble. Applying the potential control mechanism alone without the fluctuation term, i.e., the upper part of Eq. 7, dissipates thermally induced fluctuations of the total dipole moment. Thus, in absence of an explicit thermostat, the NVE ensemble shows a constant loss of kinetic energy and the system cools down, cf. Fig. 2. Without a thermostat, but using the full fluctuation-dissipation relation Eq. 7, this spurious energy transfer disappears. The temperature correctly fluctuates around the target temperature T = 350 K set in Eq. 7. This result clearly demonstrates that our potentiostat alone is not only able to control the potential, but also the temperature, even in the absence of an explicit thermostat.
We now discuss the thickness-dependence of the dielectric properties. In Fig. 3a, we compare our calculated static dielectric constant ⊥ as a function of the water layer thickness to the experimental data from Ref. [32]. Red triangles and blue squares denote data points obtained from our thermopotentiostat MD at Φ 0 = 1 V and Φ 0 = 4 V, respectively. Experimental data points measured by Fumagalli et al. [32] are shown as grey circles. Consistent with the measurements, our results display a pronounced decrease of ⊥ compared to the static dielectric constant of bulk liquid water bulk that persists for electrode separations exceeding 100 nm. Based on this level of agreement, we therefore expect our TIP3P water model to correctly capture the impact of an interface on the dielectric properties of water.
In order to understand the origin of the decreasing ⊥ with decreasing d, we computed the local static dielectric constant ⊥ (z) as a function of the normal distance z to the electrode surface [42]. Fig. 3b shows the inverse dielectric profile −1 ⊥ (z) for an electrode separation of d = 8 nm. At the position of the electrode surface z = 0, −1 ⊥ drops sharply and intersects the water bulk value at ∼3Å above the surface. With further increasing z, −1 ⊥ assumes negative values for interfacial water and then approaches the bulk water value in an oscillatory fashion. At a normal distance of ∼9Å, the dielectric constant of bulk liquid water bulk is recovered. The behaviour in the ∼9Å thick interface layer reflects the well-known density modulations of water close to interfaces [39,40], cf. Fig.  3b lower part.
We will now test whether the presence of the relatively thin layer of interfacial water with modified dielectric properties explains the observed huge decrease of the total static dielectric constant. Guided by the dielectric profile shown in Fig. 3b, we partition the dielectric profile into three regions: (i) a hydrophobic gap between electrode and surface with a thickness of d gap = 2Å, (ii) an interfacial water region consisting of the first two water layers with a thickness of d if = 5.5Å and (iii) the remaining approximately bulk-like region (Fig. 3b). The existence of a hydrophobic gap, which is clearly visible in our data, was also suggested by Niu et al. [41] based on spectroscopic data. The effective dielectric constant of each region is obtained by integrating over the dielectric profile, yielding gap = 1.2 and if = 17.3 for the hydrophobic gap and interfacial water regions, respectively. We further assume that each region is approximately independent of d. The surrogate electrostatic model becomes then a simple plate capacitor with multiple dielectrics (inset in Fig. 3a). The total dielectric constant is given by the analytical expression The solid blue line in Fig. 3a denotes ⊥ (d) obtained by the surrogate model. Although here ⊥ (d) was obtained from a single explicit calculation for an electrode separation of d = 8 nm, the solid blue line accurately reproduces all other data points obtained from explicit thermopotentiostat MD simulations at different separations d. These findings confirm that the local dielectric properties of wa-ter close to the interface are indeed responsible for the observed reduction of ⊥ compared to bulk .
The calculation of dielectric profiles from polarization fluctuations [37] requires hundreds of nanoseconds of statistical sampling, in practice enforcing the use of classical MD. In contrast, the stochastic canonical sampling of our thermopotentiostat MD in conjunction with finite electric field techniques turns out to be extremely efficient: our expression for ⊥ [42] allows us to rely purely on thermodynamic averages rather than variances. Thus, the dielectric profiles shown in Fig. 3b converged within less than 4 ns, reducing the required computational time by more than two orders of magnitude.
In conclusion, we devised a novel thermopotentiostat approach to sample the canonical ensemble at constant temperature and applied electric potential. Our approach (i) avoids the need to treat non-neutral simulation cells, (ii) requires only quantities that can be directly obtained from density-functional theory simulations and (iii) is straightforward to implement in any standard ab initio molecular dynamics package. To demonstrate the performance of our approach we computed the thicknessdependent dielectric properties of nano-confined water.
We showed that the presence of interfaces strongly modifies the dielectric constant of an interfacial water region. This region is spatially highly confined with a thickness of only ∼1 nm (roughly two water-layers). This thin region with modified dielectric properties is shown to fully capture the experimentally observed anisotropic dielectric response, persisting for distances exceeding 100 nm. The computational efficiency of our approach is improved by more than two orders of magnitude compared to previous ones. In conjunction with ab initio MD, we expect our thermopotentiostat to open the door towards accurate and efficient simulations of equilibrium properties, such as interfacial dielectric constants, as well as processes, such as electron transfer and electrochemical reactions.
We thank L. Fumagalli for providing the raw experimental data from Ref. [ In the modern theory of polarization, the Hamiltonian at constant dielectric displacement D is given by [1] where E KS and P are the Kohn-Sham energy and polarization density, respectively. For a plate capacitor, according to Gauss's Law, the value ||D|| of the dielectric displacement field between the plates is equal to the charge density n A on the plates. Analogously, we introduce an effective bound charge density np A on the plates that is equal to the value −||P|| of the polarization density between the plates. With the definition of the capacitance C 0 = 0 A d and Ω = Ad we rewrite the electrical energy E as: We note that this step is optional. Instead of deriving an equation of motion for n or Φ as done in our letter, it is also possible to use Eq. S1 in order to obtain an equivalent equation of motion for D or E in an analogous way. In the following, we use the charge n as the basic variable. This allows us to introduce the Hamiltonian: (2) The potential energy V (q i , n ex = 0) at zero excess charge n ex is simply the Born-Oppenheimer total energy surface as computed in any DFT code. If an explicit counter electrode is used in the DFT calculation, e.g., as in Ref. [2], or the generalized dipole correction [3], the Kohn-Sham energy already includes the contribution from the electric field, so that E KS = V (q i , n ex = 0) + Alternatively, the contribution from the electric field can be evaluated explicitly using the modern theory of polarization [1].

EQUATIONS OF MOTION
The complete, self-contained set of equations of motion that samples the canonical ensemble at constant potential Φ 0 is given by: Eqs. S3 and S4 are the standard Hamiltonian equations of motion, integrated by any molecular dynamics code. The term g i is used to couple to an external thermostat, cf. section Thermopotentiostat MD with an explicit thermostat. In the following, we outline how to derive the single additional equation of motion for the charge dynamics Eq. S5.
In the main text, we took the derivative d dt ∂H ∂n in order to obtain Eqs. 4 and 5 as: We then equate the current flowing through the capacitance dn/dt to the current −(Φ − Φ 0 )/R flowing in reverse through the voltage source and its internal resistivity R (Ohm's Law and Kirchhoff's 2nd Law): Here, we add a fluctuation termΦ and a stochastic noise term given by a Wiener process dW t , respectively. This term must be introduced in order to satisfy the fluctuation-dissipation theorem (FDT) [4]. Note that dW t plays an analogous conceptual role as dt. It is an integration variable and represents a stochastic time step for integration. Integration is performed by Ito calculus [5], see the Practical Implementation section. The FDT [4] determines both the variance and the frequency spectrum of the fluctuations. Classically, per frequency interval dν, the variance of the fluctuating voltage at a resistance R is given by σ 2 U dν = 4k B T Rdν [6,7]. In conjunction with the capacitance C 0 , the resistance R in our system forms an RC low pass. The variance of Φ is obtained by integrating the noise spectral density over the bandwidth of the RC low pass: Eq. S7 allows us to construct a suitable fluctuation term Φ. The energy loss in Eq. S6 due to the dissipation by the resistance R must be exactly equal, on average, to the energy gained through the fluctuations inΦ, while the variance, and also the frequency spectrum of the fluctuations, are determined by the FDT. Thereby, the applied electric field itself is fluctuating and has a finite temperature. In analogy to the Langevin [8] and BDP [9,10] thermostats, we recognize that Eq. S6 formally represents a stochastic differential equation known as the Ornstein-Uhlenbeck process The expectation value and variance of the Ornstein-Uhlenbeck process are given by [5]: With τ Φ := RC 0 , k := τ −1 Φ and D/2k := σ 2 Φ we write: Eq. 6 provides a correction potential to the Hamiltonian evolution of the potential at constant dielectric displacement (constant charge). The corresponding charge dynamics given by Eq. S5 is obtained trivially by multiplying Eq. 6 with C 0 .

THE NOSÉ-HOOVER POTENTIOSTAT AS A LIMITING CASE OF OUR THERMOPOTENTIOSTAT
Bonnet et al. proposed an empirically motivated potentiostat, introducing a fictitious momentum P n for the potential and and electronic mass M n . These are used to perform Nosé-Hoover dynamics on the charge n, witḣ n = P n M n (S11) (S12) P 2 n /(2M n ) enters the Nosé-Hoover equation of motion foṙ ξ as a fictitious kinetic energy. From our scheme, it now becomes apparent that these quantities permit a direct physical interpretation. With P n := −Rn (S13) Eq. S11 is identical to Ohm's law. The fictitious kinetic energy P 2 n /(2M n ) then becomes the electrostatic potential energy n(Φ − Φ 0 )/2. Inserting P n and M n into Eq. S12 yields: Comparing Eqs. S15 and S5, we see that the deterministic dξ in Bonnet et al.'s potentiostat is substituted in our approach by a fully stochastic term in Eq. S5.

SPECTRAL CONSIDERATIONS
The FDT determines both the variance and also the frequency spectrum of the fluctuations. In absence of a dielectric, the variance of the voltage fluctuations given by Eq. S7 is distributed in frequency space according to the spectral density function [11] of the ideal fully relaxed Ornstein-Uhlenbeck process with τ = τ Φ and c = 2k B T /(τ Φ C 0 ). In Fig. S1 we show the spectral density of the variance of the potential fluctuations for the d = 8 nm cell. Below the cutoff frequency (2πτ Φ ) −1 , the slope of the spectral density is zero (white noise region), whereas it is -2 above (1/f 2 noise region), cf., e. g., Ref. [11] for a detailed discussion. Note that the fluctuations are always independent of the dielectric properties of the system under consideration. Instead, the dielectric properties determine how the system responds to the excitation provided by the fluctuations.
Comparing the expressions for the Nose-Hoover potentiostat and our thermopotentiostat it becomes clear that Eq. S17 approximates the spectrum required by the FDT with a single discrete frequency. For poorly ergodic systems, as is the case here, it is generally recognized that a single frequency is insufficient. Using a Nose-Hoover chain as proposed by Bonnet et al. [12], the spectrum is then approximated by a set of discrete frequencies. For water in particular, which has a rather discrete spectrum itself with large phonon gaps, it is advantageous to utilize a continuous spectrum. By design, Eq. S18 satisfies the required spectrum exactly. Compared to the reversible integrator proposed by Martyna et al. [13,14] or a self-consistent solution of the Nose-Hoover equations of motion, Eq. S18 is straightforward to integrate since it has an analytical propagator, cf. next section.

PRACTICAL IMPLEMENTATION
The practical implementation of our scheme is extraordinarily easy: compared to any standard molecular dynamics code, there is only a single additional equation of motion, namely Eq. S5. The Eqs. S3 and S4 can be solved by any standard discrete time step MD scheme, e.g., using the velocity verlet algorithm. Eq. S5 is a stochastic differential equation of a specific type, the Ornstein-Uhlenbeck process: The Ornstein-Uhlenbeck process can be integrated analytically to a finite time step [15]: N is a Gaussian random number with N = 0 and N 2 = 1. With x := Φ − Φ 0 , τ := RC 0 , D := k B T /C 0 , we use Eqs. 6 and S20 to integrate Eq. S5 to: As mentioned above, the implementation of our approach is simple. For example, we performed a thermopotentiostat simulation using only the internal scripting language of the LAMMPS [16] molecular dynamics package. We provide a complete thermopotentiostat example simulation for TIP3P water sandwiched between two electrodes as part of the supplementary information. The provided input file "control.inp" can be executed by the standard LAMMPS distribution without any modifications to the LAMMPS code at all.

CHOICE OF C0
C 0 is the geometric capacitance of the two electrode system in absence of the dielectric. It can be estimated using the classical electrostatic expression C 0 = 0 A/d.
Alternatively, it can be measured by an explicit MD simulation for the bare electrodes without the dielectric using C 0 = n / Φ . It then enters Eq. 7 as a parameter and is held constant during the simulation.

THERMOPOTENTIOSTAT MD WITH AN EXPLICIT THERMOSTAT
The approach presented here does not require an additional thermostat, since the coupling to the electric field fluctuations is sufficient to thermostat the total kinetic energy of the system. Still, for computational efficiency an explicit thermostat may be desirable, e.g. to achieve fast equilibration. Introducing an explicit thermostat amounts to choosing a suitable expression for g i in Eq. S3. In standard Langevin dynamics g i is given by [10] Note the conceptual similarity to Eq. S5, where the capacitance C 0 plays the role of the mass. For the BDP thermostat, a global g is used which is the same for every kinetic degree of freedom i [10]: (S22) Here, N f is the number of kinetic degrees of freedom and K (K) is the instantaneous (mean) kinetic energy, respectively. Furthermore, we expect also the Lowe-Andersen thermostat [17] to work well in conjunction with our thermopotentiostat approach.

NUMERICAL DETAILS
Water is described by the TIP3P [18] potential. The electrodes are modeled as single sheets of fixed atoms, where interactions with water are described by the TIP3P O-O pair potential. We do not include any pair potential for the electrode atoms and the hydrogen atoms of the water molecules. The electrode thereby affects the orientation of interfacial water only by the reduced dimensionality at the interface and by the electrode charge. Periodic boundary conditions are applied in x and y direction; non-periodic boundary conditions are used along the direction of the electrode surface normal. Simulations were carried out either in the NVE ensemble or the NVT [9] ensemble at T = 350 K, using a time step of 40 atomic units (≈ 0.97 fs). All ensembles were first equilibrated without a field for 1 ns using the Langevin thermostat and subsequently for another 2 ns using the BDP thermostat. Sampling time was 4 ns. Coulomb long-range interactions were treated by a particle-particle particle mesh solver [19] with a precision of 10 −5 . All calculations were performed for rigid water in order to enable the use of a longer time step. Bond lengths and angles were constrained using the SHAKE algorithm [20,21].
The shape of the spectral density shown in Fig. S1 motivates the choice of the thermopotentiostat's relaxation time constant τ Φ . A smaller τ Φ implies a wider white noise region. Decreasing the relaxation time thereby leads to a faster sampling of the phase space, but at the same time disturbs the dynamics of the trajectory more. In practice, we chose a relaxation time of τ Φ = 100 fs. This value is large enough to ensure that the vibrational density of states obtained from our TIP3P thermopotentiostat MD simulations is virtually indistinguishable from a pure NVE simulation. Simultaneously, the required sampling time to obtain converged dielectric properties is still in range of first principles simulations.

CALCULATION OF DIELECTRIC PROFILES
In analogy to the constant charge simulations in Ref. [22], we calculate first the electric field profile ∆E ⊥ (z) induced by our thermopotentiostat's displacement field D ⊥ : Subscript 0 denotes zero external electric field. The local polarization density m ⊥ (z) is obtained from the dielectric's charge density ρ with The local dielectric profile is then extracted from the linear response relation ∆E ⊥ ≈ [ 0 ⊥ ] −1 D ⊥ . It is instructive to use our theoretical values for the dielectric constant gap and bulk to estimate if for interfacial water from Fumagalli et al.'s measurements. Assuming a hydrophobic gap of constant size with d gap = 2Å, the optimum fit is obtained for if = 2.85, cf. pink solid line in Fig. 3a, which is significantly lower than what we obtain from our thermopotentiostat simulations. Alternatively, assuming if = 17.3 for interfacial water to be independent of the specific electrode surface, a fit with identical accuracy is obtained for d gap = 4Å.