# Ions at electrochemical interfaces: From explicit to implicit molecular solvent descriptions

Swetha Nair, Guillaume Jeanmairet, Benjamin Rotenberg

<figure><img src="image_1.png" /></figure>

## To cite this version:

Swetha Nair, Guillaume Jeanmairet, Benjamin Rotenberg. Ions at electrochemical interfaces: From explicit to implicit molecular solvent descriptions. The Journal of Chemical Physics, 2025, 163 (1), 10.1063/5.0270981. hal-05142589

HAL Id: hal-05142589

https://hal.sorbonne-universite.fr/hal-05142589v1

Submitted on 3 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ions at electrochemical interfaces: from explicit to implicit molecular solvent descriptions

Swetha Nair,¹ Guillaume Jeanmairet,¹,² and Benjamin Rotenberg¹,², a)

¹) Sorbonne Université, CNRS, Physicochimie des Electrolytes et Nanosystèmes Interfaciaux, F-75005 Paris,
France

²) Réseau sur le Stockage Electrochimique de l'Energie (RS2E), FR CNRS 3459, 80039 Amiens Cedex,
France

(Dated: 19 May 2025)

We investigate the interplay between electronic screening inside a metal and screening by a polar molecular solvent, focusing on their impact on the charge induced by an ion and the solvent structure at the interface. To that end, we consider atomistically resolved electrodes within the Thomas-Fermi model of screening and describe the molecular solvent either explicitly via classical molecular dynamics or implicitly using Molecular Density Functional Theory (MDFT). Specifically, we examine the effect of screening by tuning the Thomas-Fermi screening length $l_{TF}$, the ion charge by considering $Na^+$ and $Cl^-$ and the solvent nature by studying water and acetonitrile. Consistently with our previous findings without solvent, $l_{TF}$ significantly affects the charge distribution inside the metal. However, $l_{TF}$ has no significant impact on the interfacial solvent structure, suggesting that its effect on the charge distribution induced inside the metal by the ion is essentially due to how the metal responds to the (same) external charge distribution, including the solvent, even though the coupling between both sides of the interface may play a secondary role. Furthermore, MDFT accurately reproduces fine details of the interfacial solvent structure around the ion at a fraction of the computational cost of MD simulations. These results highlight the relevance of MDFT as a powerful tool to model electrochemical systems at the molecular level.

## I. INTRODUCTION

The behavior of electrode-electrolyte interfaces is crucial for various electrochemical applications including energy storage and electrocatalysis. Their properties arise from the interplay between the metal's electronic response to charges from ions and solvent molecules, and the subsequent reorganization of these species based on the charge distribution at the metal surface. The Thomas-Fermi model¹,² provides a simple description of electronic screening within the metal over a characteristic length, $l_{TF}$, determined by the electronic density of states at the Fermi level³. On the electrolyte side, the screening of electrostatic interactions between ions in polar solvents is characterized by the Bjerrum length, while ionic screening in dilute electrolytes is governed by the Debye screening length. The interplay between these screening mechanisms plays a significant role in shaping the interfacial properties, including induced charge, ion-ion interactions, electric double layer structure, capacitance, and the overall thermodynamics and kinetics of electrochemical processes.

The mutual influence of metal, solvent and ions have been extensively studied within continuum electrostatics, notably though not only by Kornyshev and coworkers⁴⁻¹⁴. The Thomas-Fermi model describes the electronic response of the metal by a dielectric function of the form $\epsilon(k) = \epsilon_{met} \times (1 + k_{TF}^2/k^2)$, with $k_{TF} = l_{TF}^{-1}$ the reciprocal of the Thomas-Fermi screening length and

$k$ the wave vector, which ensures perfect screening for $k \ll k_{TF}$. The perfect metal limit then corresponds to $l_{TF} \to 0$. An important finding was the establishment of a characteristic length for the lateral decay of the charge/potential induced in the metal by an ion scaling as $l_{TF}\epsilon_{solv}/\epsilon_{met}$, i.e. an increase of the screening length by a factor $\epsilon_{solv}$ with respect to the absence of solvent. In parallel, the charge distribution of the metal in response to an external charge was studied using electronic Density Functional Theory (DFT) with both simplified (1D) geometries¹⁵⁻¹⁷ and with more advanced atomically resolved surfaces¹⁸,¹⁹. Recent advances in simulation techniques have further improved the modeling of electrolyte interfaces, incorporating polarization effects between media of different permittivities and explicitly accounting for ions²⁰⁻²⁷. Both ab initio²⁸⁻³⁴ and classical molecular dynamics³⁵⁻⁴⁷ have provided key insights into electrode-electrolyte interfaces and the impact of electrolyte-induced charge within the metal on interfacial properties⁴⁸⁻⁵⁷.

Recent approaches have also extended the description of Thomas-Fermi metals to account for electronic screening and its effects on interfacial thermodynamics in classical molecular simulations⁵⁸⁻⁶⁰. The metallicity of an electrode, which indicates its electronic structure, also influences electron transfer kinetics, as illustrated by recent experiments demonstrating that the kinetics of heterogeneous electron transfer reactions at well-defined twisted bilayer graphene (TBG) surfaces can be adjusted by varying the interlayer Moiré twist angle between the two layers⁶¹. Furthermore, recent molecular simulations combined with continuum dielectric theory have provided deeper insights into these observations, establish-

<footnote>a) benjamin.rotenberg@sorbonne-universite.fr</footnote>
ing a link between the twist angle and the Thomas-Fermi screening length through the density of states (DOS) of TBG$^{62}$. This relationship between DOS and Thomas-Fermi screening length has also been investigated in classical molecular simulations$^{63}$, specifically tuning the description of the charge density around electrode atoms in fluctuating charge models. Recently, we investigated the charge induced in a perfect metal ($l_{TF} = 0$) by an ion near an atomically resolved gold electrode$^{64}$ and extended this study to include the screening effect within the metal, examining the charge induced by an ion in vacuum in a Thomas-Fermi metal using a constant potential fluctuating charge model$^{65}$. We described the influence of the Thomas-Fermi screening length ($l_{TF}$) and the ion's position relative to a graphite electrode on the induced charge distribution, comparing atomistic results with predictions from continuum electrostatics.

In the present work, we further explore the interplay between electronic screening inside the metal and screening by a polar molecular solvent and their consequences on the charge induced by an ion and the solvent structure at the interface. To that end, we employ the Thomas-Fermi model for an atomistically resolved electrode, and describe the molecular solvent either explicitly using classical molecular dynamics (MD) or implicitly within classical Density Functional Theory (DFT). While classical MD simulations provide a more detailed description, sampling the equilibrium properties entails a high computational cost. Classical DFT, whereby the solvent is described implicitly via its local density instead of explicit molecules, offer a less demanding alternative, but rely on approximations typically built on Fundamental Measure Theory of hard-spheres$^{66}$ that prevent a realistic description of molecular solvents$^{47}$. Molecular Density Functional Theory (MDFT) provides an alternative route$^{47,67-79}$, based on a Taylor expansion of the functional, which allows to describe molecular systems. Building on our previous study of a model a water-graphene capacitor with polarizable electrode without considering the role of the Thomas-Fermi screening length$^{75}$, here we compare MDFT with MD to study the charge induced by an ion in a solvent within a Thomas-Fermi metal. To that end, we consider ions with different charges (Na$^+$ and Cl$^-$) in two polar solvents (water and acetonitrile) with different molecular structures and permittivities, near a graphite electrode with varying Thomas-Fermi screening length. The various methods and considered systems are presented in Section II, while the results pertaining to the effects of the screening length on the induced charge distribution and on the solvent structure are provided in Sections III and IV, respectively.

## II. METHODS AND SYSTEM

We first describe the Thomas-Fermi model used in constant-potential classical MD simulations in Sec-

tion II A and molecular DFT in Section II B. We then present the system considered in Section II C, before introducing the observables of interest in Section II D.

### A. Thomas Fermi model within constant-potential classical MD simulations

In recent years, fluctuating-charge models have been widely used for simulating electrochemical systems within classical MD simulations involving electrodes maintained at a fixed potential. First proposed by Siepmann and Sprik to study water/electrode interfaces with classical MD$^{35}$, these models were adapted to electrochemical cells$^{80}$ and more recently extended to capture Thomas-Fermi screening$^{58,60}$. Within this framework, each electrode atom $i$, positioned at $\mathbf{r}_i$ and bearing charge $q_i$, contributes to the total charge density as:

$$ \rho_{\text{elec}}(\mathbf{r}) = \sum_{i \in \text{elec}} \frac{q_i}{(2\pi w^2)^{3/2}} e^{-|\mathbf{r} - \mathbf{r}_i|^2 / 2w^2}, \quad (1) $$

where $w$ is a fixed width of the Gaussian distribution around each atom. The instantaneous charges of all electrode atoms depend on the external charge distribution, i.e. the positions $\mathbf{r}^N$ of atoms in the electrolyte. Specifically, for each configuration of the electrolyte they are determined by the constraints of constant potential imposed on each atom and of global electroneutrality of the system$^{46,81}$. The computation of the set of electrode charges $\mathbf{q} = \{q_i\}_{i \in \text{elec}}$ involves the potential energy of the system, which owing to the form of electrostatic interactions can be written as:

$$ U(\mathbf{r}^N, \mathbf{q}) = \frac{\mathbf{q}^T \mathbf{A} \mathbf{q}}{2} - \mathbf{q}^T \mathbf{B}(\mathbf{r}^N) + C(\mathbf{r}^N), \quad (2) $$

where the first term accounts for electrostatic interactions between electrode atoms, the second for those between electrode and electrolyte, and the third for all other interactions. The function $C$ is independent of electrode charges $\mathbf{q}$ and includes electrostatic interactions within the electrolyte, as well as all non-electrostatic (e.g. Lennard-Jones) interactions of the classical force field. Both the matrix $\mathbf{A}$ and the vector $\mathbf{B}$ account for the Gaussian nature of the charge distribution of electrode atoms and the periodic boundary conditions (PBC) employed in molecular simulations of condensed matter systems. While this description was initially introduced for perfect metals, Scalfi et al. extended it to the Thomas-Fermi model, by including the kinetic energy of the free electron gas as in the original Thomas-Fermi approach. By making further approximations to the description of the local density in the fluctuating charge approach, the Thomas-Fermi model was recast in the same form (Eq. 2) simply by introducing an additional term in the matrix describing electrode-electrode interactions, namely$^{58}$:

$$ \mathbf{A}(l_{\text{TF}}) = \mathbf{A}_0 + \frac{l_{\text{TF}}^2 \bar{\rho}_{at}}{\epsilon_0} \mathbf{I}, \quad (3) $$
where $A_0$ is the matrix for $l_{TF} = 0$, $\bar{\rho}_{at}$ the atomic density of the electrode, $\epsilon_0$ the vacuum permittivity and $I$ is identity matrix. This description, already implemented in several simulation packages$^{82–84}$ allows to tune the metallic character of the metal without additional computational cost.

## B. Molecular Density Functional Theory (MDFT)

One of the objectives of the present work is to assess the quality of Molecular Density Functional Theory (MDFT) to model the solvent in which the electrodes are immersed. Two particular properties should be accurately reproduced: the solvent structure and the induced charges. MDFT is based on the ansatz of classical DFT, a theory that is most naturally formulated in the grand-canonical ensemble. It states that for a given external potential $V$, there exists a unique functional $\Omega_V$ of the classical particle density $\rho$. At its minimum, which is reached for the equilibrium particle density $\rho_{eq}$, this functional is equal to the grand potential$^{85,86}$, where the subscript highlights the dependence on the external potential. In MDFT$^{47,70,74,76}$, which was originally designed to study solvation problems, the working functional is actually the difference between the grand potential functional and the grand potential of the homogeneous system $\Omega_0$, which is obtained in the absence of an external perturbation

$$F[\rho] = \Omega_V[\rho] - \Omega_0 \quad (4)$$

In the present work, the classical particle density represents the solvent *i.e* either water or acetonitrile. Since we consider rigid molecular models (see Section II C), knowledge of the position of the center of mass $r$ and the absolute orientation $\Omega = (\theta, \phi, \psi)$ of a molecule is sufficient to fully characterize the atomic coordinates. Consequently, the molecular solvent density $\rho(r, \Omega)$ also depends on the position and orientation. The external potential, which causes the solvent density to depart from homogeneity, is generated by the electrodes and by, if present, a cation or an anion. The functional of Eq. 4 is usually decomposed into the sum of several terms:

$$F[\rho(r, \Omega)] = F_{id}[\rho(r, \Omega)] + F_{exc}[\rho(r, \Omega)] + F_{ext}[\rho(r, \Omega)]. \quad (5)$$

The ideal term $F_{id}$ is the functional associated with the corresponding non-interacting fluid, while the excess term $F_{exc}$ originates from the interaction between the solvent molecules. As opposed to the ideal term, which is exactly known, the excess term requires some approximation. In this work, the excess functional is expressed as the sum of the so-called hypernetted chain (HNC) functional and a bridge functional. The HNC functional, corresponds to a Taylor expansion of the excess functional around a homogeneous reference solvent density, truncated at second order. Even at this quadratic order, the computation of the HNC functional is not a trivial task. Indeed, it requires the knowledge of the two

body correlation function of the homogeneous solvent $c^{(2)}(r, \Omega, \Omega')$ that depends on a distance and two orientations. Moreover, its calculation is computationally challenging since it involves a spatial and an orientational convolution. The expression of the HNC functional and details about its practical implementation can be found in our previous work$^{74}$. The bridge functional contains in principle all terms beyond the quadratic order and cannot be computed in practice. It is approximated here using a simple, angular independent, weighted density approximation (WDA) functional$^{77}$.

The last term of Eq. 5 is the external functional, which represents here the interaction between the solvent particles and the external potential $V$ generated by the electrodes and, if present, the cation or anion. It reads

$$F_{ext}[\rho(r, \Omega)] = \iint \rho(r, \Omega)V(r, \Omega)d\mathbf{r}d\Omega. \quad (6)$$

In the present case, $V$ is decomposed as the sum of a Lennard-Jones term $V_{LJ}$, an electrostatic term due to fixed charges $V_{ES, fixed}$ and an electrostatic term due to the fluctuating electrode charges $V_{ES, fluct}$

$$V(r, \Omega) = V_{LJ}(r, \Omega) + V_{ES, fixed}(r, \Omega) + V_{ES, fluct}(r, \Omega; q). \quad (7)$$

The external potentials due to the Lennard-Jones and fixed charges (cation or anion) are computed only once, at the beginning of the computation. For a given set of electrode charges, $q$, it is possible to compute $V_{ES, fluct}$. Minimizing the functional of equation 5, we obtain the equilibrium solvent density $\rho_{eq}(r, \Omega)$. Since the solvent density is inhomogeneous, the corresponding charge density, $\rho_c$, does not vanish

$$\rho_c(r) = \iint \rho_{eq}(r', \Omega)\sigma(r - r', \Omega)d\mathbf{r}'d\Omega. \quad (8)$$

In Eq. 8, $\sigma(r, \Omega)$ is the charge density of a water molecule with orientation $\Omega$ located at the origin, $r = 0$. The solvent charge density, in turn, impacts the electrostatic potential felt by the electrode atoms. In practice, the electrostatic potential generated by the solvent charges ($\rho_c$ discretized on a grid) enters the vector $\mathbf{B}$ in Eq. 2. Following the procedure described in previous work, we perform successive functional minimization and electrode charge computation. At convergence, we obtain the equilibrium solvent density and the electrode charges fulfilling both the electroneutrality and fixed potential conditions.

## C. Computational details

In order to study the charge induced in a Thomas-Fermi metal by an ion in a molecular solvent, we have considered the two systems depicted in Fig. 1. Both contain a metallic graphite electrode, held at a fixed potential and an insulating graphite electrode. The metallic electrode consists of 9 graphite planes with a total of
<figure><img src="image_3.png" /><figcaption>FIG. 1. Snapshots of typical molecular configurations considered in this work, with water (a) or acetonitrile (b) as solvent. Both panels include a graphite electrode with explicit atoms, an ion (here a $Na^+$ cation positioned at $z_{ion} = 5.1$ Å from the first atomic plane of the electrode, shown in cyan), and an insulating wall with the same graphite structure (shown in dark grey). The color of the electrode atoms represents the charge induced by the ion and solvent molecules, with red and blue indicating negative and positive values, respectively. For a Thomas-Fermi screening length of $l_{TF} = 0$ Å, as shown here, almost no charge is found beyond the first atomic plane in contact with the liquid. Periodic boundary conditions are applied in the *x* and *y* directions.</figcaption></figure>

4320 atoms (to allow for the possibility of spreading the induced charged within the electrode when the screening length is larger than the interplane distance), while the insulating wall is made of 3 atomic layers only, with a total of 1440 carbon atoms. The interatomic distance within each plane is $d_{CC} = 1.42$ Å, and the spacing between planes is $a = 3.354$ Å. The width of the Gaussian atomic charge entering Eq. 8 is $w = 0.4$ Å, and the Thomas-Fermi length $l_{TF}$ ranges from 0 Å to 5 Å. Both systems contain a $Na^+$ or $Cl^-$ ion either in water (2160 molecules) or in acetonitrile (789 molecules).

Periodic boundary conditions are used in the directions parallel to the walls, with lateral box dimensions of $L_x = 34.10$ Å, $L_y = 36.92$ Å. The system is not periodic in the *z* direction. The distance between the electrode and insulating wall planes in contact with the pure liquid was equilibrated for 0.5 ns in the *NPT* ensemble at 1 atm by applying a constant force to the electrodes (pistons) with $l_{TF} = 0$ Å, yielding distances of $L_z = 55.088$ Å and $L_z = 54.070$ Å for the water and acetonitrile systems, respectively. The ion is fixed at a distance $z_{ion} = 5.1$ Å from the electrode plane, above a C atom.

This distance was chosen based on equilibrium density profiles of $Na^+$ ion near the graphite electrode surface in a 1 M $NaCl$ aqueous solution (see Fig. A1 in the appendix of our previous work$^{65}$). In order to investigate the effect of the ion charge and of the solvent, we used the same positions for the $Cl^-$ anion, both in water and in acetonitrile. Finally, we also performed simulations of the same systems in the absence of ion.

Electrostatic interactions are computed using 2D Ewald summation method taking into account the Gaussian distribution of electrode atoms and the ion's point charge$^{83}$. Non-electrostatic interactions are introduced through truncated and shifted Lennard-Jones (LJ) potentials, with parameters from Ref. 87 for $Na^+$ and $Cl^-$ and from Ref. 88 for carbon (graphite electrode and opposite wall), with Lorentz-Berthelot mixing rules. Water molecules are modeled using the rigid SPC/E force field$^{89}$, while acetonitrile (MeCN) molecules are described by a rigid coarse-grained three-site model involving a unified description of the methyl group$^{90,91}$. The partial charges of the SPC/E water model are -0.8476e for oxygen and +0.4238e for hydrogen, with *e* the elementary charge; for the three-site MeCN model, they are -0.398e for nitrogen, +0.129e for the central carbon and +0.269e for the methyl group. Simulations are performed in the *NVT* ensemble for 8 ns with a time step of 2 fs, using a Nosé-Hoover chain thermostat with a time constant of 1 ps to maintain a temperature of 298 K. All simulations were performed with the molecular dynamics code Metallwalls$^{82}$, and electrode charges were calculated at each time step via matrix inversion$^{81}$, taking into account the constraint of overall electroneutrality (which implies that the total charge of the electrode compensates that of the ion) in addition to the constant-potential one. Electrode charges and solvent configurations were sampled every 100 steps (0.2 ps) and 1000 steps (2 ps), respectively, for subsequent analysis.

The same systems were studied using MDFT by removing the explicit solvent molecules and employing the functional introduced in Eq 5 to implicitly account for the solvent. The functional was minimized within a $34.10 \times 36.92 \times 120$ Å$^3$ simulation cell, discretized on a $102 \times 111 \times 120$ grid, corresponding to a resolution of approximately 0.33 Å in each direction. To address the 3D PBC used in MDFT, a hard wall perpendicular to the electrode surface was introduced to avoid self interaction. This wall also ensures that the interlayer space between graphite sheets remained free of water molecules. The three Euler angles ($\theta, \phi, \psi$) were discretized into 196 orientations, corresponding to a truncation of the expansion of the HNC functional to $m_{max} = 3$ (see Ref. 74). For the aqueous system, an angle-independent Weighted Density Approximation (WDA) bridge functional was added to the HNC functional$^{77,92}$. However, as the WDA bridge functional has not been parameterized for acetonitrile, the excess functional in this case was limited to the HNC functional.
### D. Observable properties

In order to analyze the charge distribution inside the metal and the structure of the solvent at the interface, in the following sections we will consider several observables. The equilibrium charge density is obtained by averaging Eq. 1 over the configurations in the canonical ensemble, which simply amounts to replacing the instantaneous charges $q_i$ (which depend on the microscopic configuration of the liquid) by the ensemble average $\langle q_i \rangle$. It is more conveniently analyzed by integrating over the depth of the electrode to obtain the 2D charge density

$$ \sigma_{\text{ind}}(x, y) = \int_{-\infty}^{0} \rho_{\text{elec}}(x, y, z) \, dz. \qquad (9) $$

In practice, it is more efficient to perform the integration over $z$ analytically and directly reconstruct the 2D distribution as

$$ \sigma_{\text{ind}}(x, y) = \sum_{i \in \text{elec}} \frac{\langle q_i \rangle}{2\pi w^2} e^{-\frac{[(x-x_i)^2+(y-y_i)^2]}{2w^2}}. \qquad (10) $$

This quantity is computed on a 2D grid with a bin width of 0.1 Å, where the charge in each bin is calculated using the analytical integral of the 2D Gaussian distributions (see Ref. 65). In order to analyze the radial decay away from the ion located at $(x, y, z) = (0, 0, z_{\text{ion}})$, we also introduce the radial charge density:

$$ \sigma_{\text{ind}}(r) = \frac{1}{2\pi r} \iint \sigma_{\text{ind}}(x, y) \delta\left(r - \sqrt{x^2 + y^2}\right) dx dy, \qquad (11) $$

with $\delta$ Dirac's delta function as well as its running integral:

$$ Q_{\text{ind}}(r) = \int_{0}^{r} \sigma_{\text{ind}}(r') 2\pi r' dr'. \qquad (12) $$

## III. CHARGE INDUCED INSIDE THE METAL

We first focus on the influence of both the metallicity and the nature of the solvent on the charge induced in the metal by an interfacial ion. In addition, we assess the ability of MDFT to accurately predict the induced charge by comparing the results obtained through functional minimization with those of MD simulations. The structure of the solvent in all cases will be discussed in detail in Section IV.

### A. Two-dimensional distribution of the induced charge

#### 1. Sodium ion at the water/graphite interface

Starting with the aqueous system depicted in Fig. 1a, we show in Fig. 2 the surface charge density $\sigma_{\text{ind}}(x, y)$ (see Eq. 10) induced by a $Na^+$ ion located at $(x, y) =$

(0, 0) and a distance $z_{\text{ion}} = 5.1$ Å from the first atomic plane of the electrode, for three values of the Thomas-Fermi screening: $l_{\text{TF}} = 0, 3$ and $5$ Å. The MDFT predictions (panels a-c) are in excellent agreement with the MD results (panels d-f). In all cases, the atomic structure of graphite is clearly visible in the modulation of the induced charge density. While the metal is overall negatively charged (its total charge compensates that of the $Na^+$ ion), a pattern of slightly positive local charges (in blue) is observed for $l_{\text{TF}} = 0$ Å. As already discussed in the case of gold$^{64}$, these regions correspond to atoms in the second atomic layer, with a slightly positive charge induced by the neighboring negative induced charge in the first atomic plane in contact with the liquid. The charge in this first plane is more negative below the ion and decays away from it, as discussed in more detail with the radial charge distribution in Section III B. As the screening length increases, the lateral decay becomes slower and leads to a homogeneous distribution. This behavior parallels observations made in previous studies, where increasing $l_{\text{TF}}$ led to broader charge distributions$^{65}$ and reduced charge localization, attributed to the energetic cost of high local charge densities$^{60}$.

#### 2. Sodium ion at the acetonitrile/graphite interface

We now consider a $Na^+$ ion dissolved in acetonitrile, at the same position as in the water case. Fig 3 shows the surface charge density for $l_{\text{TF}} = 0, 3$ and $5$ Å. In contrast to the water case (note in particular the different range of $\sigma_{\text{ind}}$ compared to Fig. 2), for $l_{\text{TF}} = 0$ Å, the acetonitrile system exhibits a positive charge in the first electrode plane below the ion, as shown by the blue color in panels (a) and (d), surrounded by a localized negative charge, followed by a smooth decay away from the ion. As will be discussed in Section IV, this reflects the presence of acetonitrile molecules with electronegative nitrogen atoms solvating the cation. This subtle effect is well reproduced by MDFT, even though the local charge is more positive below the ion than in the MD case. MDFT also captures well the effect of screening inside the metal, which leads to a more delocalized charge distribution (in particular there is no sign change, unlike for $l_{\text{TF}} = 0$ Å), even though the agreement with MD is not quantitative. In particular, the ability of MDFT to accurately capture the complex pattern of the induced charge density indicates a proper description of the solvation structure (see Section IV B).

### B. Radial distribution of the induced charge

#### 1. Sodium ion in water and acetonitrile

In order to facilitate the discussion of the effect of the Thomas-Fermi screening length and the comparison between MDFT and MD, we now turn to the radial induced
<figure><img src="image_4.png" /><figcaption>FIG. 2. Surface charge density $\sigma_{ind}(x, y)$ (see Eq. 9) induced by a Na$^+$ cation in liquid water located at $(x, y) = (0, 0)$ and a distance $z_{ion} = 5.1$ Å from the first atomic plane of the electrode, for $l_{TF} = 0, 3$ and $5$ Å. Panels (a)-(c) indicate the results from MDFT, panels (d)-(f) the ones from MD. The colormap is common for all panels.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. Surface charge density $\sigma_{ind}(x, y)$ (see Eq. 9) induced by a Na$^+$ cation in liquid acetonitrile located at $(x, y) = (0, 0)$ and a distance $z_{ion} = 5.1$ Å from the first atomic plane of the electrode, for $l_{TF} = 0, 3$ and $5$ Å. Panels (a)-(c) indicate the results from MDFT, panels (d)-(f) the ones from MD. The colormap is common for all panels.</figcaption></figure>

charge density $\sigma_{ind}(r)$ and its running integral $Q_{ind}(r)$ introduced in Eqs. 11 and 12 respectively. The results are shown in Fig. 4, for Na$^+$ in water in panels 4a and 4b, and in acetonitrile in panels 4c and 4d. In each panel, the color indicates the value of $l_{TF}$ and MD results (symbols) are compared with the MDFT ones (dashed lines). As discussed previously in the case of vacuum$^{65}$, the oscillations observed in $\sigma_{ind}(r)$ with both solvents (panels 4a and 4c) arise from the atomic structure of the electrode. Their magnitude reflects the influence of the solvent. The ability of MDFT to capture these oscillations further establishes its suitability to describe the solvent implicitly

in realistic atomistic models (ion and explicit electrode atoms). These oscillations are smoothed out in the running integral, $Q_{ind}(r)$, which does not converge to $-q_{ion}$ (which corresponds to the whole electrode) but rather to a value close to $-\frac{\pi}{4}q_{ion}$ because the corners of the electrode are not included when integrating radially up to half of the simulation box. Overall, MDFT is able to reproduce the MD ones, including both the effect of $l_{TF}$ and that of the solvent (water vs acetonitrile).

A closer examination of the water case (panels 4a and 4b) reveals a small discrepancy for $l_{TF} = 0$ Å. This is probably due to subtle differences in the solvent struc-
<figure><img src="image_6.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0270981

<figure><img src="image_7.png" /><figcaption>FIG. 4. (a, c) Radially averaged induced surface charge density $\sigma_{ind}(r)$ (see Eq. 11), and (b, d) the corresponding running integrals (see Eq. 12) for the charge density induced in a graphite electrode by a $Na^+$ cation located at $(x, y) = (0, 0)$ and positioned $z_{ion} = 5.1$ Å above the electrode's first atomic plane. Colors refer to Thomas-Fermi lengths $l_{TF} = 0, 3$ and $5$ Å as indicated in panel (b). Panels (a) and (b) refer to $Na^+$ in water, while (c) and (d) refer to $Na^+$ in acetonitrile.</figcaption></figure>

ture predicted by MD and MDFT, mainly in the first adsorbed layer at the graphite interface and in the first solvation shell around the sodium ion. This is a known limitation of the excess functional, which underestimates the tetrahedral order around charged solutes. As $l_{TF}$ increases, the local charge below the ion ($r = 0$) decreases, while the lateral extent of the charge distribution broadens and becomes more homogeneous due to the combined effects of the screening inside the metal and the presence of the solvent.

The results for acetonitrile (panels 4c and 4d) are similar. As already observed in the 2D maps of Fig. 3a and 3, for $l_{TF} = 0$ Å $\sigma_{ind}(r)$ displays a positive region near $r = 0$. In that case, MDFT overestimates $\sigma_{ind}(r)$ compared to MD, in line with previous findings that MDFT tends to overestimate solvation structures relative to MD simulations⁶⁸. As $l_{TF}$ increases, the local charge below the ion ($r = 0$) decreases, while the lateral extent of the charge distribution broadens, improving the agreement between MDFT and MD results. Interestingly, for larger $r$ values, $l_{TF}$ exhibits a reversed ordering compared to the water solvent. For the radial integral $Q_{ind}(r)$, the agreement between MDFT and MD remains quantitative.

As mentioned in the introduction, Kornyshev and coworkers⁶ predicted using continuum electrostatics that the characteristic length for the long-range lateral decay

<figure><img src="image_8.png" /><figcaption>FIG. 5. (a) Radially averaged induced surface charge density $\sigma_{ind}(r)$ (see Eq. 11), and (b) the corresponding running integrals (see Eq. 12) for the charge density induced in a graphite electrode by a $Na^+$ cation located at $(x, y) = (0, 0)$ and positioned $z_{ion} = 5.1$ Å above the electrode's first atomic plane. Symbols are the results in the absence of solvent, while dashed and solid lines correspond to acetonitrile and water, respectively. Colors refer to Thomas-Fermi lengths $l_{TF} = 0, 3$ and $5$ Å as indicated in panel (b).</figcaption></figure>

of the charge/potential induced in the metal by an ion scales as $l_{TF}\epsilon_{solv}/\epsilon_{met}$. The MD results in the presence of water and acetonitrile (already shown in Fig. 4) are compared to the vacuum case, discussed in more detail in Ref. 65, in Fig. 5. The influence of $l_{TF}$ is much more visible (in particular in the running integral of panel 5b) in vacuum than in the presence of a solvent. For both solvents, the broadening of the distribution with respect to the vacuum case, consistent with the increase predicted from continuum electrostatics, is in fact such that the induced charge density is almost homogeneous beyond a few Å from the ion, regardless of $l_{TF}$. Such a homogeneity (also visible in Figs. 2 and 3) is due to the large decay length compared to the lateral dimensions of the box (see Ref. 65 for more discussion of the effect of periodic boundary conditions in the vacuum case) and prevents us from a quantitative test of the scaling as $l_{TF}\epsilon_{solv}/\epsilon_{met}$. We note that this predicts no effect of the solvent perturbation in the $l_{TF} = 0$ case, unlike what is observed in simulations. However this prediction is only valid for a single ion in an infinite system, in the limit of large distances from the ion, so that it is not in contradiction with the MD result that are sensitive to the periodic boundary conditions in the lateral directions.

## 2. Chloride ion in water and acetonitrile

To further investigate MDFT's predictive capabilities, we replace $Na^+$ with $Cl^-$ at the same position in both solvents. While the 2D distributions of the induced charge are similar for MDFT and MD (see Section S1 of the supplementary material), the radial distributions $\sigma_{ind}(r)$ indicate some quantitative deviations, most notably for water (see Figs. 6a and 6b). This is another consequence
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0270981

<figure><img src="image_9.png" /><figcaption>FIG. 6. (a, c) Radially averaged induced surface charge density $\sigma_{ind}(r)$ (see Eq. 11), and (b, d) the corresponding running integrals (see Eq. 12) for the charge density induced in a graphite electrode by a $Cl^-$ anion located at $(x, y) = (0, 0)$ and positioned $z_{ion} = 5.1$ Å above the electrode's first atomic plane. Colors refer to Thomas-Fermi lengths $l_{TF} = 0, 3$ and $5$ Å as indicated in panel (b). Panels (a) and (b) refer to $Cl^-$ in water, while (c) and (d) refer to $Cl^-$ in acetonitrile.</figcaption></figure>

of the limitations of the excess functional, which is known to perform better for positively charged solutes than for negatively charged ones, due to its shortcomings for the description of hydrogen bonding and to the orientation of OH bonds toward the anion. In acetonitrile, where no such complex bonding interactions occur, the functional is more accurate, and we observe excellent agreement between MD and MDFT for both the radial distribution and the its running integral (see Figs. 6c and 6d).

## IV. SOLVENT STRUCTURE

Overall, the above results confirm both the impact of metallicity on the charge induced inside the metal and the consistency between MDFT and MD for this observable. In this section, we explore how the Thomas-Fermi screening length influences the structure of the solvent at the interface with graphite, first in the absence of ion, then investigating the solvation of cations and anions at the interface.

### A. Interface between graphite and pure solvents

#### 1. Solvent layering

Fig. 7 presents a comparison of water's oxygen ($\rho_O$) and hydrogen ($\rho_H$) density profiles computed using both MD and MDFT. The two methods yield qualitatively similar results, with two distinct molecular layers, beyond which the homogeneous densities are rapidly recovered. However, a closer examination reveals some differences. For the oxygen density in Fig. 7a, the two main peaks of the density profile are smaller in MDFT compared to MD. In Fig. 7b, the maximum of the first hydrogen density peak is located at 10.5 Å in MDFT, slightly closer to the electrode than in MD (10.8 Å). In addition, this first peak is sharper in MD than in MDFT.

<figure><img src="image_10.png" /><figcaption>FIG. 7. Oxygen (a) and hydrogen (b) density profiles at the graphite/water interface for $l_{TF}=0$ Å. The MD results are shown in red, while the MDFT results are represented by black dashed lines. Panels (c) and (d) compare the oxygen density profiles for different Thomas-Fermi lengths $l_{TF}$ (indicated by the color), for MD and MDFT, respectively, while hydrogen density profiles are shown in panels (e) and (f) for MD and MDFT, respectively. The profiles are zoomed in for clarity.</figcaption></figure>

Furthermore, we compare the oxygen (Figs. 7c and 7d) and hydrogen (Figs. 7e and 7f) density profiles for various values of $l_{TF}$ ranging from 0 to 5 Å as indicated in the color bar. The results are shown for both MD and MDFT. We observe that the solvent density remains largely unchanged with varying $l_{TF}$ values. This aligns closely with the observations of Scalfi et al.⁵⁸, where water density near the electrode showed minimal variation
<figure><img src="image_11.png" /></figure>

with different $l_{TF}$ values in the absence of applied voltage across the interface. This suggest that, under these conditions, the metallicity of the electrode has a limited impact on the structure of the solvent at the interface. This is further supported by the fact that the O and H densities on the neutral wall (which corresponds to the insulating limit $l_{TF} \to \infty$), visible on the right sides of panels Figs. 7a and 7b, are comparable to that on the metallic one (only symmetric with respect to the center of the pore). This conclusion might not hold for non-zero applied voltage$^{58}$.

<figure><img src="image_12.png" /><figcaption>FIG. 8. Nitrogen (a), carbon (b), and methyl group (c) density profiles at the acetonitrile/graphite interface for $l_{TF} = 0$ Å, as obtained by MDFT (black dashed lines) and MD (blue solid lines).</figcaption></figure>

We also consider the predictions by MDFT and MD of the nitrogen (Fig. 8a), carbon (Fig. 8b), and methyl group (Fig. 8c) density profiles at the acetonitrile/graphite interface for $l_{TF} = 0$ Å. Both MDFT and MD predict an interfacial structure with 3-4 distinct layers, with density maxima of decreasing intensity away from the interface. However, closer examination reveals small differences between the predictions of the two methods. For all profiles, the first solvation peak predicted by MDFT is positioned closer to the electrode than with MD. Notably, the carbon and methyl group profiles reveal significant differences between MDFT and MD, with peak maxima at 4.06 Å for MDFT and 4.36 Å for MD. For nitrogen, the difference is smaller, with peak maxima at 4.06 Å for MDFT and 4.16 Å for MD. Much like water, the solvent structure remains largely unaffected by changes in the Thomas-Fermi length (not shown).

## 2. Solvent structure in the first adsorbed layer on graphite

<figure><img src="image_13.png" /><figcaption>FIG. 9. Lateral correlations between solvent molecules in the first adsorbed layer on the graphite electrode, for (a) water, with the 2D radial distribution for oxygen atoms, and (b) acetonitrile, with the 2D radial distribution for the central carbon atoms, for a range of $l_{TF}$ indicated by the color bar.</figcaption></figure>

Since the electrode metallicity has a clear impact on the surface charge density within the electrode, as discussed in Section III, whereas it virtually has no effect on the density profiles in the direction perpendicular to the surface, we now examine its possible influence on the 2D structure of the solvent within the first adsorbed layer (defined from the density profiles of Figs. 7 and 8). The 2D radial distribution functions between oxygens of water molecules and between central carbon atoms of acetonitrile molecules within that layer are shown presented in Figs. 9a and Figs. 9b, respectively, for Thomas-Fermi screening lengths varying from 0 to 5 Å. There is a clear in-plane structure, with a first peak in the radial distribution functions located at 2.7 Å for water and 4.8 Å for acetonitrile. These distances are not due to a templating effect due to the graphite lattice (the 2D densities
along the plane do not display any particular feature, not shown) and rather reflect the lateral correlations between molecules sliding over the surface. While the presence of the wall does not change the typical distance between water molecules compared to the bulk, it results in an increase of the linear acetonitrile molecules (≈ 4 Å in the bulk) due to packing constraints. Importantly, we observe no significant effect of the Thomas-Fermi length on this in-plane solvent structure, in contrast with that of the induced charge. This is consistent with previous findings on a nanocapacitor consisting of a 1 M aqueous NaCl solution between graphite electrodes⁶⁰.

## B. Ion solvation at the solvent/graphite interface

We now turn to the solvation of ions by water and acetonitrile at the interface with graphite, by considering the prototypical cases sodium cation and chloride anion. The effect of the ionic charge is twofold: Firstly, it structures the solvent in its vicinity and secondly, it induces a non-zero charge on the metallic electrode, which may in turn modify the structure of the interfacial solvent. The induced charge was described in Section III and we consider here only the solvent structure. To that end, we introduce the radially average solvent charge density in planes perpendicular to the electrode containing the ion:

$$ \rho_{\text{solv}}(r, z) = \left\langle \sum_{k} q_{k} \delta(r_{k} - r) \delta(z_{k} - z) \right\rangle, \quad (13) $$

where the sum runs over solvent atoms with partial charge $q_k$ and position $(r_k, z_k)$ expressed in cylindrical coordinates. We first examine the effect of the ion charge for $l_{TF} = 0$ Å, before turning to the effect of $l_{TF}$, and compare MDFT with MD in all cases.

### 1. Effect of ion charge

Fig. 10 shows the average solvent charge density in the $(r, z)$ plane for $Na^+$ and $Cl^-$ ions in water and acetonitrile. The alternating signs of the charge density in panels 10a, 10b, 10e and 10f reflect the characteristic orientation of water molecules around the ions and at the graphite surface. Far from the ions, the charge density profiles depend only on $z$ (indicating limited lateral ordering) and are consistent with the atomic density profiles of Fig. 7, with some dangling OH bonds extending out the first molecular layer. Close to the ion, water molecules are oriented with their oxygen atoms pointing toward the sodium cation, while the preferred orientation is flipped in the vicinity of the chloride atom, as expected. MDFT is able to reproduce faithfully the rather complex solvation behaviour for both ions, although minor differences arise in the magnitude of the charge densities.

Similar conclusions can be drawn in the acetonitrile case, shown in panels 10c, 10d, 10g and 10h. Far from

the ions, we do not observe the alternating positively and negatively charged planes perpendicular to the electrode. Instead, we observe an oscillating charge within the first adsorbed layer and a more diffuse charge pattern beyond, suggesting the packing of the linear molecules parallel to the surface, with the alternating charge pattern further reflecting some dipolar ordering of the molecules. The nitrogen atom is positioned closer to the $Na^+$ cation, while it is the methyl groups which points towards $Cl^-$. This can be attributed to acetonitrile's linear structure and the higher electronegativity of nitrogen (as described within this classical force field by its negative partial charge) compared to the carbon atom and methyl group. MDFT effectively captures these asymmetric solvation behaviors resulting with an overall good agreement with MD.

Finally, It is worth noticing that the noisy patterns far from the ion and the electrode in MD is absent in MDFT, which predicts a well defined charge density in the whole $(r, z)$ plane. Indeed, the solvent density predicted by MDFT is obtained by minimizing the functional, while the one predicted by MD suffer from statistical error (which can be reduced using force-based estimators⁹³⁻⁹⁶). This is particularly true in region where there is no preferred solvent orientation, which is the case far from the solute and the electrode.

### 2. Effect of Thomas-Fermi screening length

Finally, we examine the effect of the Thomas-Fermi screening length on the solvent structure near the ion. To that end, we extend the previous analysis (Fig. 10) to three different $l_{TF}$ values: 0, 3, and 5 Å. Fig. 11 shows the average charge density of water molecules in the $(r, z)$ plane for a $Na^+$ ion at these $l_{TF}$ values. We observe no significant change in the solvent charge density when varying $l_{TF}$ in both MDFT and MD. This is consistent with the pure solvent case discussed in Section IVA, as well as with the results of Loche et al.⁵⁴ showing a very limited effect of the metallicity of graphite on the potential of mean force for the adsorption of an ion in water. These observations on the solvent structure also hold in acetonitrile, as shown in Fig. 12) and or the $Cl^-$ anion in both solvents (see Supplementary Figs. S2 and S3). In all cases, MDFT consistently captures the asymmetric solvation structure and solvent patterns, demonstrating its robustness in reproducing solvation behavior accurately.

Overall, the results from Section IV indicate that there is no significant effect of $l_{TF}$ on the interfacial structure (at least in the absence of applied voltage). This suggests that the effect of $l_{TF}$ on the charge distribution induced inside the metal by the ion (Section III) is essentially due to how the metal responds to the (same) external charge distribution, even though the coupling between both sides of the interface may play a secondary role.
<figure><img src="image_14.png" /><figcaption>FIG. 10. Average charge density distribution (see Eq. 13) of solvent molecules around Na⁺ and Cl⁻ ions in water (a,b,e,f) and acetonitrile (c,d,g,h) for *l*<sub>TF</sub> = 0 Å obtained using MDFT (a-d) and MD (e-h). Insets illustrate typical configurations of solvent molecules in the first solvation shell, with water represented using red for oxygen and white for hydrogen, and acetonitrile using yellow for nitrogen, orange for carbon, and purple for the methyl group.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 11. Average charge density of water molecules in the (*r*, *z*) plane (see Eq. 13) for the Na⁺ ion at three different Thomas-Fermi screening lengths (*l*<sub>TF</sub> = 0, 3, and 5 Å). Panels (a), (c), and (e) show the MDFT results, while panels (b), (d), and (f) present the MD results.</figcaption></figure>

V. CONCLUSION AND PERSPECTIVES

We investigated the charge induced in a metallic elec-
trode by an ion in solution. Using both classical molecu-
lar simulations and molecular Density Functional Theory,
we explored the impact of several factors. The effect of
the screening inside the metal was studied by tuning the
Thomas-Fermi length. The role of the ion, as well as the

impact of the solvent, were assessed by considering both
a Na⁺ cation and a Cl⁻ anion dissolved in either water or
acetonitrile. Consistently with previous molecular-scale
studies on the interface between a Thomas-Fermi metal
and vacuum or water, we find that at *l*<sub>TF</sub> = 0 Å, the
induced surface charge is highly localized below the ion
and spreads as *l*<sub>TF</sub> increases. The main conclusions of
the present work are twofold.
<figure><img src="image_16.png" /><figcaption>FIG. 12. Average charge density of acetonitrile molecules in the (r, z) plane (see Eq. 13) for the Na⁺ ion at three different Thomas-Fermi screening lengths (l<sub>TF</sub> = 0, 3, and 5 Å). Panels (a), (c), and (e) show the MDFT results, while panels (b), (d), and (f) present the MD results.</figcaption></figure>

Firstly, there is no significant effect of the Thomas-Fermi screening length on the interfacial structure (at least in the absence of applied voltage). This suggests that the effect of $l_{TF}$ on the charge distribution induced inside the metal by the ion is essentially due to how the metal responds to the external charge distribution. The coupling between both sides of the interface may play a secondary role. Since the interfacial structure also depends on the identity of the metal, it would be interesting to further explore the effect of the electrode lattice structure and of the description of metal-solvent interactions beyond the ones considered here – using the present classical constant-potential approach (e.g. tuning the LJ interactions⁹⁷, or introducing three-body potentials as was done for Pt-water interactions³⁵), or more elaborate QM/MM (quantum mechanics / molecular mechanics) descriptions²⁹,³²,⁹⁸,⁹⁹.

Secondly, MDFT is able to capture accurately the solvation of ions in water and acetonitrile at an atomistically resolved electrode. It reproduces subtle features of the solvent structure at a fraction of the computational cost required to converge the 3D structure with MD simulations. MDFT typically allows for a reduction of the computational cost by three orders of magnitude. As an example, in the pure water case it took roughly 7400 CPUh to obtain 8 ns of simulation, while the functional optimization requested 6 CPUh. While we have focused here on the induced charge and the solvent structure around an interfacial ion, we are currently also exploring the ability of MDFT to reproduce the effect of metallicity on interfacial thermodynamics, which can be

investigated with MD as described in Ref 60.

As discussed in Ref. 65 for the charge induced in an atomistically resolved electrode by an ion in vacuum, the predictions of continuum electrostatics capture most of the features observed with the atomistic description (except the oscillations due to the atomic sites of the graphite lattice), provided that the ion-surface distance and the Thomas-Fermi screening length are larger than the inter-atomic distances within the electrode. This requires however a careful definition of the position of the effective interface between the metal and vacuum. The present work provides a molecular basis (both with MD and MDFT) for a systematic quantitative comparison with analytical predictions from continuum electrostatics in the presence of an implicit solvent, either described only by its permittivity⁶ or by more elaborate models interpolating between the long- and short-range screening properties of the solvent¹⁰⁰⁻¹⁰⁴.

Since the induced charges and solvent structure are well reproduced, the next step is to consider thermodynamic and kinetic properties. The potential of mean force as a function of the distance from the surface, for which semi-analytical models have been proposed¹⁰⁵, is particularly suited for MDFT since the solvation free energy is a direct output of the functional minimization. Moreover we established a framework to compute Marcus free energy curves with MDFT, allowing to access the activation free energy. MDFT provides a computationally efficient method to access these quantities, while maintaining a molecular description of the solvent⁷⁶,⁷⁸. The present work shows that it should be possible to also in-

<footnote>This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0270981</footnote>
clude the impact of metallicity by tuning of the Thomas-Fermi length into such such studies, at a more affordable cost than MD simulations⁶². This would open the way to a systematic study of the influence of the nature of the ion (size and charge), of the solvent and of the surface (atomic structure and screening length) on adsorption ion and electron transfer kinetics.

## SUPPLEMENTARY MATERIAL

Additional results (charged induced inside the metal and solvent structure) in the case of a Cl⁻ anion are provided as Supplementary Material.

## ACKNOWLEDGMENTS

We thank Luc Belloni for useful discussions and for providing us with the MDFT functional for acetonitrile. This project received funding from the European Research Council under the European Union's Horizon 2020 research and innovation program (grant agreement no. 863473). GJ acknowledges financial support from the French Agence Nationale de la Recherche (project ANR BRIDGE AAP CE29). This work was performed using HPC resources from GENCI-IDRIS (Grant 2024-AD010912966R2).

## AUTHOR DECLARATIONS

### Conflict of interest

There is no conflict of interest to declare

### Author contributions

**Swetha Nair:** Conceptualization (supporting); Formal analysis (equal); Investigation (lead); Methodology (supporting); Software (supporting); Writing/Original Draft Preparation (lead); Writing - review & editing (equal). **Guillaume Jeanmairet:** Conceptualization (equal); Formal analysis (equal); Investigation (supporting); Methodology (equal); Software (lead); Supervision (supporting); Writing - review & editing (equal). **Benjamin Rotenberg:** Conceptualization (equal); Formal analysis (equal); Funding Acquisition (lead); Investigation (supporting); Methodology (equal); Supervision (lead); Writing - review & editing (equal).

## DATA AVAILABILITY

**To be completed consistently with the final version of the manuscript:** The original data presented in this study are openly available in

Zenodo at [DOI/URL] or [reference/accession number].
¹L. Thomas, Proc. Cambridge Phil. Roy. Soc. **23**, 542 (1927).

²E. Fermi, Rend. Accad. Naz. Lincei **6**, 602 (1927).

³N. W. Ashcroft and N. D. Mermin, *Solid state physics* (Holt, Rinehart and Winston, New York, NY, 1976).

⁴A. A. Kornyshev, A. I. Rubinshtein, and M. A. Vorotyntsev, *physica status solidi* (b) **84**, 125 (1977).

⁵A. A. Kornyshev and M. A. Vorotyntsev, J. Phys. C: Solid State Phys. **11**, L691 (1978).

⁶M. A. Vorotyntsev and A. A. Kornyshev, Soviet Journal of Experimental and Theoretical Physics **51**, 509 (1980).

⁷A. A. Kornyshev, W. Schmickler, and M. A. Vorotyntsev, Phys. Rev. B **25**, 5244 (1982).

⁸P. G. Dzhavakhidze, A. A. Kornyshev, and L. I. Krishtalik, J. Electroanal. Chem. **228**, 329 (1987).

⁹D. K. Phelps, A. A. Kornyshev, and M. J. Weaver, J. Phys. Chem. **94**, 1454 (1990).

¹⁰N. B. Luque and W. Schmickler, *Electrochim. Acta* **71**, 82 (2012).

¹¹A. A. Kornyshev, N. B. Luque, and W. Schmickler, J. Solid State Electrochem. **18**, 1345 (2014).

¹²C. C. Rochester, A. A. Lee, G. Pruessner, and A. A. Kornyshev, ChemPhysChem **14**, 4121 (2013).

¹³J. Comtet, A. Niguès, V. Kaiser, B. Coasne, L. Bocquet, and A. Siria, Nat. Mater. **16**, 634 (2017).

¹⁴V. Kaiser, J. Comtet, A. Niguès, A. Siria, B. Coasne, and L. Bocquet, Faraday Discuss. **199**, 129 (2017).

¹⁵N. D. Lang and W. Kohn, Phys. Rev. B **7**, 3541 (1973).

¹⁶W. Schmickler and D. Henderson, The Journal of Chemical Physics **80**, 3381 (1984).

¹⁷N. V. Smith, C. T. Chen, and M. Weinert, Phys. Rev. B **40**, 7565 (1989).

¹⁸J. Jung, J. E. Alvarellos, E. Chacón, and P. García-González, *Journal of Physics: Condensed Matter* **19**, 266008 (2007).

¹⁹E. K. Yu, D. A. Stewart, and S. Tiwari, Phys. Rev. B **77**, 195406 (2008).

²⁰S. Tyagi, A. Arnold, and C. Holm, J. Chem. Phys. **127**, 154723 (2007).

²¹S. Tyagi, A. Arnold, and C. Holm, J. Chem. Phys. **129**, 204102 (2008).

²²A. Arnold, K. Breitsprecher, F. Fahrenberger, S. Kesselheim, O. Lenz, and C. Holm, Entropy **15**, 4569 (2015).

²³M. Girotto, A. P. dos Santos, and Y. Levin, The Journal of Chemical Physics **147**, 074109 (2017).

²⁴T. D. Nguyen, H. Li, D. Bagchi, F. J. Solis, and M. Olvera de la Cruz, Computer Physics Communications **241**, 80 (2019).

²⁵O. Maxian, R. P. Peláez, L. Greengard, and A. Donev, The Journal of Chemical Physics **154**, 204107 (2021).

²⁶J. Yuan, H. S. Antila, and E. Luijten, The Journal of Chemical Physics **154**, 094115 (2021).

²⁷A. P. dos Santos, F. Jiménez-Ángeles, A. Ehlen, and M. Olvera de la Cruz, Physical Review Research **5**, 043174 (2023).

²⁸S. Sakong and A. Groß, J. Chem. Phys. **149**, 084705 (2018).

²⁹J. D. Elliott, A. Troisi, and P. Carbone, Journal of Chemical Theory and Computation **16**, 5253 (2020).

³⁰J.-B. Le, A. Chen, L. Li, J.-F. Xiong, J. Lan, Y.-P. Liu, M. Ian-nuzzi, and J. Cheng, JACS Au **1**, 569 (2021).

³¹S. Sakong, J. Huang, M. Eikerling, and A. Groß, Current Opinion in Electrochemistry **33**, 100953 (2022).

³²K. Takahashi, H. Nakano, and H. Sato, The Journal of Chemical Physics **157**, 234107 (2022).

³³R. Sundararaman, D. Vigil-Fowler, and K. Schwarz, Chemical Reviews **122**, 10651 (2022).

³⁴A. Grisafi, A. Bussy, M. Salanne, and R. Vuilleumier, Physical Review Materials **7**, 125403 (2023).

³⁵J. I. Siepmann and M. Sprik, The Journal of Chemical Physics **102**, 511 (1995).

³⁶E. Spohr, *Electrochimica Acta* **44**, 1697 (1999).

³⁷A. A. Kornyshev, E. Spohr, and M. A. Vorotyntsev, “Electrochemical interfaces: At the border line,” in *Encyclopedia of Electrochemistry* (John Wiley & Sons, Ltd, 2007) pp. 33–132.

³⁸S. K. Reed, P. A. Madden, and A. Papadopoulos, The Journal of Chemical Physics **128**, 124701 (2008).

³⁹M. V. Fedorov and A. A. Kornyshev, Chem. Rev. **114**, 2978 (2014).

⁴⁰N. Onofrio and A. Strachan, The Journal of Chemical Physics **143**, 054109 (2015).

⁴¹K. Breitsprecher, K. Szuttor, and C. Holm, J. Phys. Chem. C **119**, 22445 (2015).

⁴²T. Liang, A. C. Antony, S. A. Akhade, M. J. Janik, and S. B. Sinnott, J. Phys. Chem. A **122**, 631 (2018).

⁴³H. Nakano and H. Sato, The Journal of Chemical Physics **151**, 164123 (2019).

⁴⁴K. Takahashi, H. Nakano, and H. Sato, The Journal of Chemical Physics **157**, 014111 (2022).

⁴⁵R. S. Sitlapersad, A. R. Thornton, and W. K. den Otter, The Journal of Chemical Physics **160**, 034107 (2024).

⁴⁶L. Scalfi, M. Salanne, and B. Rotenberg, Annual Review of Physical Chemistry **72**, 189 (2021).

⁴⁷G. Jeanmairet, B. Rotenberg, and M. Salanne, Chemical Reviews **122**, 10860 (2022).

⁴⁸A. P. Willard, S. K. Reed, P. A. Madden, and D. Chandler, Faraday Discuss. **141**, 423 (2009).

⁴⁹D. T. Limmer, A. P. Willard, P. Madden, and D. Chandler, Proc. Natl. Acad. Sci. U.S.A. **110**, 4200 (2013).

⁵⁰A. P. Willard, D. T. Limmer, P. A. Madden, and D. Chandler, J. Chem. Phys. **138**, 184702 (2013).

⁵¹E. Paek, A. J. Pak, and G. S. Hwang, J. Chem. Phys. **142**, 024701 (2015).

⁵²M. Buraschi, S. Sansotta, and D. Zahn, J. Phys. Chem. C **124**, 2002 (2020).

⁵³A. Serva, L. Scalfi, B. Rotenberg, and M. Salanne, The Journal of Chemical Physics **155**, 044703 (2021).

⁵⁴P. Loche, L. Scalfi, M. Ali Amu, O. Schullian, D. J. Bonthuis, B. Rotenberg, and R. R. Netz, The Journal of Chemical Physics **157**, 094707 (2022).

⁵⁵S. Ntim and M. Sulpizi, Physical Chemistry Chemical Physics **22**, 10786 (2020).

⁵⁶S. Ntim and M. Sulpizi, Physical Chemistry Chemical Physics **25**, 22619 (2023).

⁵⁷S. Ntim and M. Sulpizi, The Journal of Physical Chemistry B **128**, 1936 (2024).

⁵⁸L. Scalfi, T. Dufils, K. G. Reeves, B. Rotenberg, and M. Salanne, The Journal of Chemical Physics **153**, 174704 (2020).

⁵⁹A. Schlaich, D. Jin, L. Bocquet, and B. Coasne, Nature Materials **21**, 237 (2021).

⁶⁰L. Scalfi and B. Rotenberg, Proceedings of the National Academy of Sciences **118**, e2108769118 (2021).

⁶¹Y. Yu, K. Zhang, H. Parks, M. Babar, S. Carr, I. M. Craig, M. Van Winkle, A. Lyssenko, T. Taniguchi, K. Watanabe, V. Viswanathan, and D. K. Bediako, Nature Chemistry **14**, 267 (2022).

⁶²L. Coello Escalante and D. T. Limmer, Nano Letters **24**, 14868 (2024).

⁶³K. Goloviznina, J. Fleischhaker, T. Binninger, B. Rotenberg, H. Ers, V. Ivanistsev, R. Meissner, A. Serva, and M. Salanne, Adv. Mater. **36**, 2405230 (2024).

⁶⁴G. Pireddu, L. Scalfi, and B. Rotenberg, The Journal of Chemical Physics **155**, 204705 (2021).

⁶⁵S. Nair, G. Pireddu, and B. Rotenberg, Molecular Physics, e2365990 (2024).

⁶⁶Y. Rosenfeld and N. W. Ashcroft, Physical Review A **20**, 1208 (1979).

⁶⁷R. Ramirez, R. Gebauer, M. Mareschal, and D. Borgis, Phys. Rev. E **66**, 031206 (2002).

⁶⁸S. Zhao, R. Ramirez, R. Vuilleumier, and D. Borgis, The Journal of Chemical Physics **134**, 194102 (2011).

⁶⁹D. Borgis, L. Gendre, and R. Ramirez, The Journal of Physical Chemistry B **116**, 2504 (2012).
70G. Jeanmairet, M. Levesque, R. Vuilleumier, and D. Borgis, The Journal of Physical Chemistry Letters **4**, 619 (2013).

71G. Jeanmairet, M. Levesque, and D. Borgis, The Journal of Chemical Physics **139**, 154101 (2013).

72G. Jeanmairet, M. Levesque, V. Sergiievskyi, and D. Borgis, The Journal of Chemical Physics **142**, 154112 (2015).

73G. Jeanmairet, N. Levy, M. Levesque, and D. Borgis, Journal of Physics: Condensed Matter **28**, 244005 (2016).

74L. Ding, M. Levesque, D. Borgis, and L. Belloni, The Journal of Chemical Physics **147**, 094107 (2017).

75G. Jeanmairet, B. Rotenberg, D. Borgis, and M. Salanne, The Journal of Chemical Physics **151**, 124111 (2019).

76G. Jeanmairet, B. Rotenberg, M. Levesque, D. Borgis, and M. Salanne, Chem. Sci. **10**, 2130 (2019).

77D. Borgis, S. Luukkonen, L. Belloni, and G. Jeanmairet, The Journal of Chemical Physics **155**, 024117 (2021).

78T.-Y. Hsu, R. Berthin, A. Serva, K. Reeves, M. Salanne, and G. Jeanmairet, The Journal of Chemical Physics **157**, 094103 (2022).

79D. Borgis, D. Laage, L. Belloni, and G. Jeanmairet, Chem. Sci. **14**, 11141 (2023).

80S. K. Reed, O. J. Lanning, and P. A. Madden, J. Chem. Phys. **126**, 084704 (2007).

81L. Scalfi, D. T. Limmer, A. Coretti, S. Bonella, P. A. Madden, M. Salanne, and B. Rotenberg, Physical Chemistry Chemical Physics **22**, 10480 (2020).

82A. Marin-Laflèche, M. Haefele, L. Scalfi, A. Coretti, T. Dufils, G. Jeanmairet, S. K. Reed, A. Serva, R. Berthin, C. Bacon, S. Bonella, B. Rotenberg, P. A. Madden, and M. Salanne, Journal of Open Source Software **5**, 2373 (2020).

83A. Coretti, C. Bacon, R. Berthin, A. Serva, L. Scalfi, I. Chubak, K. Goloviznina, M. Haefele, A. Marin-Laflèche, B. Rotenberg, S. Bonella, and M. Salanne, The Journal of Chemical Physics **157**, 184801 (2022).

84L. J. V. Ahrens-Iwers, M. Janssen, S. R. Tee, and R. H. Meißner, The Journal of Chemical Physics **157**, 084801 (2022).

85N. D. Mermin, Phys. Rev. **137**, A1441 (1965).

86R. Evans, Adv. Phys. **28**, 143 (1979).

87L. X. Dang, Journal of the American Chemical Society **117**, 6954 (1995).

88T. Werder, J. H. Walther, R. L. Jaffe, T. Halicioglu, and P. Koumoutsakos, The Journal of Physical Chemistry B **107**, 1345 (2003).

89H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, The Journal of Physical Chemistry **91**, 6269 (1987).

90P. A. M. Duncan M.F. Edwards and I. R. McDonald, Molecular Physics **51**, 1141 (1984).

91C. Merlet, M. Salanne, B. Rotenberg, and P. A. Madden, Electrochimica Acta **101**, 262 (2013).

92D. Borgis, S. Luukkonen, L. Belloni, and G. Jeanmairet, The Journal of Physical Chemistry B **124**, 6885 (2020).

93D. Borgis, R. Assaraf, B. Rotenberg, and R. Vuilleumier, Molecular Physics **111**, 3486 (2013).

94S. W. Coles, D. Borgis, R. Vuilleumier, and B. Rotenberg, The Journal of Chemical Physics **151**, 064124 (2019).

95B. Rotenberg, The Journal of Chemical Physics **153**, 150902 (2020), publisher: American Institute of Physics.

96S. W. Coles, E. Mangaud, D. Frenkel, and B. Rotenberg, The Journal of Chemical Physics **154**, 191101 (2021), publisher: American Institute of Physics.

97C. Bacon, A. Serva, C. Merlet, P. Simon, and M. Salanne, Electrochimica Acta **455**, 142380 (2023).

98J. D. Elliott, M. Chiricotto, A. Troisi, and P. Carbone, Carbon **207**, 292 (2023).

99A. Grisafi and M. Salanne, The Journal of Chemical Physics **161**, 024109 (2024), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0218379/20037849/024109_1.5.0218379.pdf.

100H. Berthoumieux and F. Paillusson, The Journal of Chemical Physics **150**, 094507 (2019).

101M. Vatin, A. Porro, N. Sator, J.-F. Dufrêche, and H. Berthoumieux, Molecular Physics **119**, e1825849 (2021), publisher: Taylor & Francis eprint: https://doi.org/10.1080/00268976.2020.1825849.

102G. Monet, F. Bresme, A. Kornyshev, and H. Berthoumieux, Physical Review Letters **126**, 216001 (2021), publisher: American Physical Society.

103J. G. Hedley, H. Berthoumieux, and A. A. Kornyshev, The Journal of Physical Chemistry C **127**, 8429 (2023), publisher: American Chemical Society.

104J. G. Hedley, K. Coshic, A. Aksimentiev, and A. A. Kornyshev, Physical Review X **14**, 031042 (2024).

105A. Robert, H. Berthoumieux, and M.-L. Bocquet, Physical Review Letters **130**, 076201 (2023), publisher: American Physical Society.