Simulating electrochemical systems by combining the finite field method with a constant potential electrode

TL;DR

This paper introduces a method combining finite electric field application with constant potential electrodes in molecular dynamics to better simulate electrochemical interfaces, capturing ion adsorption and charge fluctuations.

Contribution

It presents a novel simulation approach that enables fixed-potential electrochemical interface modeling within classical molecular dynamics frameworks.

Findings

Successfully simulates electrode charge fluctuations and ion adsorption.

Creates paired electric double layers with aligned dipoles.

Facilitates efficient electrochemical interface simulations.

Abstract

A better understanding of interfacial mechanisms is needed to improve the performances of electrochemical devices. Yet, simulating an electrode surface at fixed electrolyte composition remains a challenge. Here we apply a finite electric field to a single electrode held at constant potential and in contact with an aqueous ionic solution, using classical molecular dynamics. The polarization yields two electrochemical interfaces on opposite sides of the same metal slab. While the net charge on one electrode surface is the opposite of the net charge on the other, maintaining overall charge neutrality of the metal. The electrode surface charges fluctuations are compensated by the adsorption of ions from the electrolyte, forming a pair of electric double layers with aligned dipoles. This opens the way towards the efficient simulation of electrochemical interfaces using any flavor of molecular dynamics, from classical to first principles-based methods.

Equations

$$\rho (\mathbf{r})=\sum _{i=1}^N{q}_i\delta (\mathbf{r}-{\mathbf{r}}_i)+\sum _{j=1}^M{q}_j{\eta }^3{\pi }^{3/2}\exp \left[-{\eta }^2(\mathbf{r}-{\mathbf{r}}_j{)}^2\right]$$

$$H^{PBC}=K({\mathbf{p}}^N)+U({\mathbf{r}}^N,{\mathbf{q}}^M),$$

$${\Psi }_j^{PBC}=\frac{\partial U_C^{PBC}}{\partial q_j}={\Psi }_{\mathrm{J}}$$

$$U_C^{PBC}=\frac{1}{2}\iint \frac{\rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime })}{4\pi {ϵ}_0|\mathbf{r}-{\mathbf{r}}^{\prime }|}\mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }.$$

$$H_E=H^{PBC}-\Omega \mathbf{P}\cdot \mathbf{E},$$

$${\mathbf{P}}_{\mathrm{itinerant}}(t)={\mathbf{P}}_{\mathrm{itinerant}}(0)+\frac{1}{\Omega }\sum _{i=1}^N{q}_i\Delta {\mathbf{r}}_i(t),$$

$${\mathbf{P}}_{\mathrm{conductor}}(t)=\frac{1}{\Omega }\sum _{j=1}^M{q}_j(t){\mathbf{r}}_j.$$

$${\Psi }_j={\Psi }_j^{PBC}-{\mathbf{r}}_j\cdot \mathbf{E}={\Psi }_{\mathrm{J}}$$

$$\frac{\mathrm{d}{E}_z}{\mathrm{d}z}=4\pi \rho (z),$$

$$C=\frac{1}{S}\frac{⟨Q^2⟩-⟨Q{⟩}^2}{k_{B}T},$$

Full text

Despite many advances over the past decades Calle-Vallejo and Koper (2012); Magnussen and Groß (2019), the efficient simulation of full electrochemical cells at the molecular scale, using electronic structure based calculations, remains a daunting task. This is due to their slab structure, since the minimal experimental setup consists of an electrolyte between two electrodes. The system is generally simplified by simulating one interface only, but the main conceptual difficulty is to find a way to charge the electrode surface at fixed composition of the electrolyte. Several methods have recently emerged Bonnet et al. (2012); Melander et al. (2019); Nishihara and Otani (2017), where the system is allowed to exchange electrons with a reservoir at fixed voltage (grand-canonical approach), but they all rely on the use of continuum descriptions for the electrolyte. These models, which are generally based on a Poisson-Boltzmann theory Nattino et al. (2019), remain mostly qualitative and an atomistic description would be preferable (this is also true because the solvent may actively participate to electrochemical reactions Le et al. (2017); Bouzid and Pasquarello (2018)). This is almost impossible to do since it would be necessary to remove/insert ions to counterbalance the electrode charge fluctuations.

Here we propose an alternative route to simulate electrochemical cells. Our approach is based on the coupling of a finite field with a system consisting in an electrolyte and a single electrode. Finite fields methods, developed in the framework of the modern theory of polarization, consist in imposing a macroscopic field (electric field Souza et al. (2002), polarization Diéguez and Vanderbilt (2006) or electric displacement Stengel et al. (2009)) via an extended Hamiltonian accounting for the interaction between the system and the fixed field. They have recently been adapted and applied to the study of electrical double layers at solid/liquid interfaces, and more precisely charged Zhang and Sprik (2016) or polar Sayer et al. (2017) insulators/electrolyte interfaces. Electrochemical systems are by nature more complex since they involve metallic electrodes whose charge distribution is not fixed but depends on the surrounding medium and the applied potential. Due to the long simulation times related to the relaxation of the electrical double layer (which is typically longer than the nanosecond), we establish here a proof of concept by using a classical molecular dynamics (MD) setup. Indeed, even if metals can only be accurately described in the framework of quantum mechanics, models have been developed to reproduce the electrostatics in classical or mixed quantum/classical simulations Siepmann and Sprik (1995); Golze et al. (2013). In this contribution we extend the finite electric field method to an electrolyte interacting with such a model metallic electrode.

In the following, we focus on the description of the electrostatics, the Van der Waals interactions being represented by the conventional Lennard-Jones model. The charge density in any point of space is given by

[TABLE]

where the first term is the contribution of the electrolyte, represented by a distribution of point charges with $q_{i}$ the partial charge of the atom $i\in[1,N]$ and ${\bf r}_{i}$ its position; $\delta$ is the Dirac distribution. The second term represents the atoms of the metallic electrode, in which each site $j\in[1,M]$ is immobile (with position ${\bf r}_{j}$) and carries a charge $q_{j}$ which is spatially distributed following a gaussian charge distribution of width $\eta^{-1}$. In order to represent the metallic character of the electrodes, the latter charges are allowed to fluctuate in response to the electrolyte fluctuations and thus are part of the microscopic degrees of freedom Siepmann and Sprik (1995); Limmer et al. (2013). The Hamiltonian of the system is written as

[TABLE]

where $K(\mathbf{p}^{N})$ is the kinetic energy which depends on the ion momenta $\mathbf{p}^{N}=\{\mathbf{p}_{1}...\mathbf{p}_{N}\}$ and $U(\mathbf{r}^{N},\mathbf{q}^{M})$ the potential energy, which depends on the ion positions $\mathbf{r}^{N}=\{\mathbf{r}_{1}...\mathbf{r}_{N}\}$ and charges of the electrode $\mathbf{q}^{M}=\{q_{1}...q_{M}\}$. We have appended a superscript PBC (periodic boundary conditions) to indicate that the electrostatic energies and forces are computed using standard Ewald summation to account. All the atoms within an electrode are held at constant potential by enforcing the following condition on each atom

[TABLE]

where $\Psi_{\rm J}$ is the prescribed potential of electrode J to which the atom $j$ belongs and $U^{PBC}_{C}$ is the Coulombic contribution to the energy, given by:

[TABLE]

Formally, solving the set of self-consistent equations given by Equation 3 is equivalent to minimizing $U^{PBC}_{C}-\sum_{j}\Psi_{\rm J}q_{j}$ with respect to the charges. Since this function is quadratic in the fluctuating charges $q_{j}$ (see Eq 4), the minimization can be efficiently performed with conjugate gradients. Note that we add an additional constraint by forcing the sum of the electrode charges to be null Gingrich and Wilson (2010).

The minimal components of an electrochemical cell are two electrodes and an electrolyte between them. In classical MD, the computational cost is not prohibitive so it is relatively easy to simulate complete systems instead of a single electrode surface. The conventional setup to simulate such systems is illustrated in the top panel of Figure 1. In the following we will refer to this system as the electrolyte-centered supercell (ECS). It is simulated using the 2D Ewald summation Kawata and Mikami (2001); Reed et al. (2007) since the two electrodes are held at different potentials ($\Psi_{\rm J}=\Psi_{\rm left}$ or $\Psi_{\rm right}$). We note $x$ and $y$ the two directions along which PBCs are used in this setup.

Finite field ($\mathbf{E}$) simulations can be performed using the extended Hamiltonian introduced by Stengel and Vanderbilt in Stengel et al. (2009) and is written as:

[TABLE]

where $H^{PBC}$ is the Hamiltonian defined by Eq 2, $\Omega$ is the volume of the supercell and $\mathbf{P}$ the polarization per unit volume. In the modern theory of polarization, the dipole moment of a unit cell involving infinite periodic systems is viewed as a multivalued quantity, since it depends on the choice of the position of the periodic boundaries. Nevertheless, this is not a significant issue since only differences in polarization matter in the dynamics and in the calculation of physical properties, in practice via the itinerant polarization Caillol (1994):

[TABLE]

where $\Delta\mathbf{r}_{i}$ is the displacement of the atom between time $t=0$ and $t$ for the ”unfolded” trajectory, i.e. not taking jumps in position (hence polarization) across the periodic boundaries. The system must be periodic in the direction in which the finite electric field is applied, which implies the use of 3D PBCs. A field $E$ corresponds to a drop of Poisson potential across the cell $\Delta\Psi=-EL_{z}$ where $L_{z}$ is the length of the box in the direction of the field. From the practical point of view, a consequence of the use of 3D PBCs is that now the two electrodes of the ECS necessarily merge, yielding a single electrode at fixed potential $\Psi_{\rm J}=\Psi_{\rm electrode}$. The simulation cell can then be represented with the electrode at its center, yielding the conductor-centered supercell (CCS) shown in the bottom panel of Figure 1. Contrarily to the constant applied potential simulations, it is now the presence of the finite field that induces the polarization of the electrode and potential drop at the two electrode/electrolyte interfaces.

Coupling fluctuating charges to model conductors with finite field simulations requires two important adaptations of both methods. On the one hand, the cell polarization includes a contribution from the fluctuating charges as:

[TABLE]

This additional term does not depend on the position of the electrode inside the supercell since we enforce $\sum_{j}q_{j}=0$. On the other hand, the determination of the partial charges via the constraint of fixed potential includes an additional contribution to the electrostatic energy and potential due to the finite field. From the extended Hamiltonian (Eq 5) and the expression of the electrode contribution to polarization (Eq 7), one obtains the generalization of Eq 3 as:

[TABLE]

As an extension of the constant applied potential case, the self-consistent expressions given by Eq 8 are now solved by minimizing $U_{C}^{PBC}-\Omega~{}\mathbf{P}\cdot\mathbf{E}-\sum_{j}\Psi_{\rm J}q_{j}$ with respect to the charges. Since ${\bf P}$ is a linear function of the charges, the performances of the conjugate gradient minimizer are not affected by this setup.

In order to test this new approach, we simulate two systems, one at constant applied potential between two distinct electrodes and the other with the finite field method and a single fixed-potential electrode as shown on Figure 1. The electrode(s) consist in a model structure made of a cubic crystal with the NaCl lattice constant, with the (111) plane facing the liquid. The intermolecular interactions consist in electrostatic interactions and Lennard-Jones potentials using Lorentz-Berthelot mixing. The electrode sites Lennard-Jones parameters are the ones of Cl-  Joung and Cheatham (2008). The electrolyte consists in an aqueous solution of NaCl composed of 603 SPC/E water molecules Berendsen et al. (1987) and 20 ion pairs Joung and Cheatham (2008). The cross-sectional area is 2.20 nm2 and the length along the $z$ axis $L_{z}=11$ nm.

The simulations were performed with a timestep of 2 fs in the NVT ensemble at 298 K using a Nosé-Hoover thermostat with a coupling constant of 0.4 ps. The systems were equilibrated during 10 ns before a production run of 10 ns.

For the finite field simulations, we use the CCS configuration with 3D PBCs. The single electrode, which is made of 12 planes of atoms, is set at null potential and the field $E$ ranges from 0 to 2 V nm*-1*. For the constant applied potential simulations, we setup the ECS configuration with 2D PBCs. In this setup the two electrodes are of equal dimensions (i.e. 6 planes of atoms each), kept under constant potentials $\Psi_{\rm left}$ and $\Psi_{\rm right}$ such that $\Delta\Psi=\Psi_{\rm left}-\Psi_{\rm right}=-EL_{z}$, using the same values for $E$ as above. In both series of simulations, the value of the $\eta$ parameter for the gaussian charges has been set to 0.5052 Å -1 following ref. Reed et al. (2007). Electroneutrality is enforced during the charge calculation process Gingrich and Wilson (2010).

A first validation is provided by comparing the polarization of the electrodes with the two setups. To do this we split the single electrode in the CCS setup in two parts. This is easily made since one side accumulates positive charge and the other is exactly opposite, while the centre is almost neutral. Firstly, the average accumulated charge $Q$ on the positive side is compared with the ECS setup on Figure 2A for a wide range of applied fields (potentials). The agreement between the two methods is excellent since the points are almost superimposed. This is true not only for voltages up to 6 V where the charge increases linearly with the applied potential (which reflects a constant differential capacitance Limmer et al. (2013)), but also up to 20 V for which the variation is not linear anymore (see the inset). The probability distributions of the instantaneous values of this quantity (Figure 2B) are also identical, which shows that the two methods sample the same configuration space. Finally, Figure 2C illustrates the instantaneous electrode charges for a given electrolyte configuration. The two systems are indistinguishable, showing that even at the local scale the finite field method yields a correct representation of the electrode/electrolyte interface.

When a finite field is applied to a bulk liquid, due to the PBCs the charges cannot accumulate in a specific region so that there is a net electric field in each point of the supercell. Here the presence of a blocking surface (the electrode) results in the formation of polarized layers on the electrolyte side. The polarization arises from two mechanisms: i) reorientation of the water molecules ii) local charge imbalance by accumulation of one ionic species and depletion of the other. The structure adopted by the liquid may be compared with the case of constant applied potential simulations. As shown on Figure 3, the agreement is again very good for the variation of the charge density and its splitting between water and ionic contributions across the simulation cell. This validates further the use of the single electrode in order to study electrochemical interfaces using 3D PBCs. Additional tests on the electric field and Poisson potential for several applied voltages are provided in the Supplementary Information; they all show the same level of accuracy.

In conclusion, we have demonstrated the possibility to simulate metal slabs with fluctuating surface charge and in the presence of an explicit electrolyte that counterbalances the charge. This is done through the combination of two methods: An applied finite field which polarizes the cell and a constant potential electrode that screens the field in bulk, leading to the formation of two independent electrochemical interfaces. The net charge on one electrode surface is the opposite of the net charge on the other, which maintains the overall charge neutrality of the metal slab. The electrode surface charge fluctuations are compensated by the adsorption of ions from the electrolyte. The method is validated through extensive comparisons with simulations using constant applied potential between two separated electrodes (and 2D PBCs). From the practical point of view, it is much easier to implement in classical MD packages since it avoids the introduction of 2D PBCs. From the computational point of view, the method is also more efficient: The simulation time is reduced by approximately 15 %. Last, but not least, it could much be more easily applied to the case of ab initio molecular dynamics, which would open the door for the first-principles simulation of electrochemical reactions occurring at electrodes, in the presence of an explicit electrolyte.

Supplementary Information

Constant potential and applied electric field

As seen in the main text, fixed field and constant potential agree very well on the charge distribution (Fig 4 A). From here one obtains the electric field profile in the z direction, averaged on the x and y directions, using the Maxwell-Gauss equation

[TABLE]

with the condition $\langle E_{z}\rangle=E$, the applied field. Fig 4 B shows that both methods agree once again. In particular, they both display a null electric field in the bulk electrode and electrolyte as expected for conductors. This is highlighted when we calculate from the electric field profile the Poisson potential profile displayed on fig 4 C, since the potential is flat in the electrode for both methods. We also observe the same potential drop for both methods between the electrode and the electrolyte at each of the interfaces. Thus the capacitance of each half-electrode in the 3D case is the same as the 2D equivalent electrode. The agreement between applied field and constant potential is thus extended to the Poisson potential.

Field profiles for the single elctrode

One feature we need to check is if the single electrode under an applied field does behave like a perfect metal. This implies that the electric field in the bulk of the electrode is really null for every applied electric field. The electric field profile is displayed for the single electrode for multiple values of the applied electric field on fig 5.

Whatever the applied field, the field in the bulk electrode is equal to zero after the second atomic layer. A closer look on this field may be obtained by the calculation of the Poisson potential. This profile is displayed for the single electrode for multiple values of the applied electric field on fig 6.

In the electrode, the Poisson potential profile is flat, corresponding indeed to a zero electric field in the electrode, which does behave like a perfect metal. The plateau observed allows to assign a value of the Poisson potential to the single electrode. We also notice that because of the charge density fluctuations in the electrolyte, the potential difference $\Delta\Psi$ between the two bulk part of the electrolyte on each side of the electrode is not strictly equal to $-EL_{z}$. This slight difference is otherwise not significant and is lower than the fluctuations of $\Delta\Psi$.

Fluctuation-dissipation relation

The capacitance of the system sudied in the main text may be obtained from the $Q(\Delta\Psi)$ relation displayed in the main text, and gives a value of $C=1.94\pm 0.11~{}\mu\mathrm{F.cm}^{-2}$. It may also be calculated using the fluctuation-dissipation theorem:

[TABLE]

where $S$ is the cross-section and $T$ the temperature. The obtained values and a comparison with the fit of the $Q(\Delta\Psi)$ relation are displayed on figure 7, showing that the two approaches are equivalent. This is not surprising since we observed in the main text that the charge distributions for constant electric field and constant applied potential are very similar in terms of both average values and standard deviation. This confirms that the capacitance does not depend on the applied voltage in the range [0:5V], that we could infer from $Q(\Delta\Psi)$ relation.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 771294).