The *ab initio* computational treatment of electrochemical systems requires an appropriate treatment of the solid/liquid interfaces. A fully quantum mechanical treatment of the interface is computationally demanding due to the large number of degrees of freedom involved. In this work, we develop a computationally efficient model where the electrode part of the interface is described at the density-functional theory (DFT) level, and the electrolyte part is represented through an implicit solvation model based on the Poisson-Boltzmann equation. We describe the implementation of the linearized Poisson-Boltzmann equation into the Vienna *Ab initio* Simulation Package, a widely used DFT code, followed by validation and benchmarking of the method. To demonstrate the utility of the implicit electrolyte model, we apply it to study the surface energy of Cu crystal facets in an aqueous electrolyte as a function of applied electric potential. We show that the applied potential enables the control of the shape of nanocrystals from an octahedral to a truncated octahedral morphology with increasing potential.

## I. INTRODUCTION

Interfaces between dissimilar systems are ubiquitous in nature and frequently encountered and employed in scientific and engineering applications. Of particular interest are solid/liquid interfaces that form the crux of many technologically essential systems such as electrochemical interfaces,^{1–3} corrosion, lubrication, and nanoparticle synthesis.^{4–6} A comprehensive understanding of the behavior and properties of such systems requires a detailed microscopic description of the solid/liquid interfaces.

Given the heterogeneous nature and the complexities of solid/liquid interfaces,^{7} experimental studies often need to be supplemented by theoretical undertakings. A complete theoretical investigation of such an interface based on an explicit description of the electrolyte and solute involves solving systems of nonlinear equations with a large number of degrees of freedom, which quickly becomes prohibitively expensive even when state-of-the-art computational resources are employed.^{8} This calls for the development of efficient and accurate implicit computational frameworks for the study of solid/liquid interfaces.^{8–14}

We have previously developed a self-consistent implicit solvation model^{15,16} that describes the dielectric screening of a solute embedded in an implicit solvent, where the solute is described by density-functional theory (DFT). We implemented this solvation model into the widely used DFT code VASP.^{17} This software package, VASPsol, is freely available under the open-source Apache License, Version 2.0 and has been used by us and others to study the effect of the presence of solvent on a number of materials and processes, such as the catalysis of ketone hydrogenation,^{18} ionization of sodium superoxide in sodium batteries,^{19} the chemistry at electrochemical interfaces,^{20,21} hydrogen adsorption on platinum surfaces,^{22} the surface stability and phase diagram of solvated mica,^{23} etc.^{24,25}

In this work, we describe the extension of our solvation model VASPsol^{15,16} to include the effects of mobile ions in the electrolyte through the solution of the linearized Poisson-Boltzmann equation and then apply the method to electrochemical systems to demonstrate the power of this approach. Figure 1 illustrates how the solid/electrolyte system is divided into spatial regions that are described by the explicit DFT and the implicit solvent. This solvation model allows for the description of charged systems and the study of the electrochemical interfacial systems under applied external voltage. Section II outlines the formalism and presents the derivation of the energy expression. Section III describes the details of the implementation, and Sec. IV validates the method against previous computational studies of electrochemical systems. Finally, in Sec. V, we apply the approach to study how an applied electric potential can control the equilibrium shape of Cu nanocrystals.

## II. THEORETICAL FRAMEWORK

Our solvation model couples an implicit description of the electrolyte to an explicit quantum-mechanical description of the solute. The implicit solvation model incorporates the dielectric screening due to the permittivity of the solvent and the electrostatic shielding due to the mobile ions in the electrolyte. The solute is described explicitly using DFT.

To obtain the energy of the combined solute/electrolyte system, we spatially divide the material system into three regions: (i) the solute that we describe explicitly using DFT, (ii) the electrolyte that we describe using the linearized Poisson-Boltzmann equation, and (iii) the interface region that electrostatically couples the DFT and Poisson-Boltzmann equation. Figure 1 illustrates the regions for a Pt(111) surface embedded in a 1M aqueous electrolyte. The interface between the electrolyte and the solute system is formed self-consistently through the electronic density, and the interaction between the electrolyte and the solute is described by the electrostatic potential as explained in Sec. II A.

### A. Total energy

We define the interface region between the solute and the electrolyte using the electronic charge density of the solute, $n(r\u2192)$. The electrolyte occupies the volume of the simulation cell where the solute’s electronic charge density is essentially zero. In the interface region where the electronic charge density increases rapidly from zero toward the change from their values inside of the solute, the properties of the implicit solvation model, e.g., the permittivity, *ϵ*, and the Debye screening length, *λ*_{D}, change from their bulk values of the electrolyte to the vacuum values in the explicit solute region. The scaling of the electrolyte properties is given by a shape function,^{10,15,26–28}

Following our previous work^{15} and the work by Arias *et al.*,^{10,26–28} the total free energy, *A*, of the system consisting of the solute and electrolyte is expressed as

where *A*_{TXC} is the kinetic and exchange-correlation contribution from DFT, *ϕ* is the net electrostatic potential of the system, and *ρ*_{s} and *ρ*_{ion} are the total charge density of the solute and the ion charge density of the electrolyte, respectively.

The total solute charge density, *ρ*_{s}, is the sum of the solute electronic and nuclear charge densities, $n(r\u2192)$ and $N(r\u2192)$, respectively,

The ion charge density of the electrolyte, *ρ*_{ion}, is given by

where $ci$ is the concentration of ionic species *i*, *z*_{i} denotes the formal charge, and *q* is the elementary charge. The concentration of ionic species, *c*_{i}, is given by a Boltzmann factor of the electrostatic energy^{29} that is modulated in the interface region by the shape function, $\zeta (n(r\u2192))$,

where $ci0$ is the bulk concentration of ionic species *i*, *k*_{B} is the Boltzmann constant, and *T* is the temperature.

The relative permittivity of the electrolyte, $\u03f5(r\u2192)$, is assumed to be a local functional of the electronic charge density of the solute and modulated by the shape function,^{10,15,26–28}

where *ϵ*_{b} is the bulk relative permittivity of the solvent.

The cavitation energy, *A*_{cav}, describes the energy required to form the solute cavity inside the electrolyte and is given by

where the effective surface tension parameter *τ* describes the cavitation, the dispersion, and the repulsion interaction between the solute and the solvent that are not captured by the electrostatic terms.^{13}

In comparison with our previous solvation model,^{15} the main changes to the free energy expression are the inclusion of the electrostatic interaction of the ion charge density in the electrolyte with the system’s electrostatic potential, *ϕ*, and the nonelectrostatic contribution to the free energy from the mobile ions in the electrolyte, *A*_{ion}. In our work, we assume this nonelectrostatic contribution from the ions to consist of just the entropy term,

where *S*_{ion} is the entropy of mixing of the ions in the electrolyte, which for small electrostatic potentials can be approximated to first order as^{29}

### B. Minimization

To obtain the stationary point of the free energy, *A*, given by Eq. (2), we set the first order variations of the free energy with respect to the system potential, *ϕ*, and the electronic charge density, *n*(*r*), to zero.^{10,15,26–28} Taking the variation of $A[n(r\u2192),\varphi (r\u2192)]$ with respect to the electronic charge density, $n(r\u2192)$, yields the typical Kohn-Sham Hamiltonian^{30} with the following additional term in the local part of the potential:

Taking the variation of $A[n(r\u2192),\varphi (r\u2192)]$ with respect to $\varphi (r\u2192)$ yields the generalized Poisson-Boltzmann equation,

where *ϵ*(*n*(*r*)) is the relative permittivity of the solvent as a local functional of the electronic charge density.

We further simplify the system by considering the case of electrolytes with only two types of ions present, whose charges are equal and opposite, i.e., $c10$ = $c20$ = *c*^{0} and *z*_{1} = −*z*_{2} = *z*. Then, the ion charge density of the electrolyte becomes^{29}

and the Poisson-Boltzmann equation becomes

For small arguments $x=zq\varphi kBT\u226a1$, sinh(*x*) → *x* and we obtain the linearized Poisson-Boltzmann equation

with

where *λ*_{D} is the Debye length that characterizes the dimension of the electrochemical double layer.

The linearized Poisson-Boltzmann equation given by Eqs. (14) and (15) and the additional local potential, *V*_{solv}, of Eq. (10) define our implicit electrolyte model. This model has four key approximations: (i) The ionic entropy term of Eq. (9) treats the ionic solution as an ideal system and assumes that any interactions between the ions are small compared to *k*_{B}*T*. (ii) The cations and anions of the electrolyte have equal and opposite charges, simplifying the Poisson-Boltzmann equation. (iii) The electrostatic potential in the electrolyte region is small such that *zq*Φ ≪ *k*_{B}*T*, which leads to the linear approximation of the Poisson-Boltzmann Eq. (14). To be quantitative, arguments to the sinh function in Eq. (13) of less than 0.25 lead to errors below 1%. Hence, in the region where the shape function is unity (and thus $\u2207\u2192\u22c5\u03f5\u22480$), charge excesses of up to 0.5 *c*^{0} are still faithfully reproduced. (iv) The ions are treated as point charges. The finite size of the ions in solution would limit their maximum concentration at the interface. These four approximations could be overcome and will be considered in future extensions of the VASPsol model.

## III. IMPLEMENTATION

We implement the implicit electrolyte model described above into the widely used DFT software Vienna *Ab initio* Software Package (VASP).^{17} VASP is a parallel plane-wave DFT code that supports both ultrasoft pseudopotentials^{31,32} and the projector-augmented wave (PAW)^{33} formulation of pseudopotentials. Our software module, VASPsol, is freely distributed as an open-source package and hosted on GitHub at https://github.com/henniggroup/VASPsol.^{16}

The main modifications in the code are the evaluation of the additional contributions to the total energy and the local potential, given by Eqs. (2) and (10), respectively. Corrections to the local potential require the solution of the linearized Poisson-Boltzmann equation given by Eq. (13) in each self-consistent iteration. Our implementation solves the equation in reciprocal space and makes efficient use of fast Fourier transformations (FFTs).^{15} This enables our implementation to be compatible with the Message Passing Interface (MPI) and to take advantage of the memory layout of VASP. We use a preconditioned conjugate gradient algorithm to solve the linearized Poisson-Boltzmann equation with the preconditioner $G2+\kappa 2\u22121$, where *G* is the magnitude of the reciprocal lattice vector and *κ*^{2} is the inverse of the square of the Debye length, *λ*_{D}.

The solution of the linearized Poisson-Boltzmann Eq. (13) provides a natural reference electrostatic potential by setting the potential to zero in the bulk of the electrolyte. However, plane-wave DFT codes, such as VASP, implicitly set the average electrostatic potential in the simulation cell to zero, not the potential in the electrolyte region. For simplicity, we do not modify the VASP reference potential but provide the shift in reference potential that needs to be added to the Kohn Sham eigenvalues and the Fermi level. Furthermore, the shift of the electrostatic potential to align the potential in the electrolyte region to zero, or to any other value, modifies the energy of the system. This energy change is given by Δ*E*_{ref} = *n*_{electrode} Δ*U*_{ref}, where *n*_{electrode} is the net charge of the electrode slab and Δ*U*_{ref} is the change in reference potential, i.e., the shift to align the potential in the electrolyte region to zero.

To validate the change in reference potential, Fig. 2 shows the grand-canonical electronic energy, *F*(*U*), as a function of the applied potential, *U*, of the Pt (111) electrode slab. The grand-canonical electronic energy, *F*, is the Legendre transformation of the free energy, A, of the system, *F*(*U*) = *A*(*n*) − *n*_{electrode}*U*, and is always lowered when charge is added or removed from the neutral slab. The grand canonical energy, *F*, exhibits the expected quadratic behavior^{34} and the maximum at the neutral slab when the energy change Δ*E*_{ref} due to the change in reference potentials is included.

## IV. VALIDATION

We further validate the model by comparison with existing experimental and computational data. First, we compute the potential of zero charge (PZC) of various metallic slabs and compare them against the experimental values. Second, we calculate the surface charge density of a Pt(111) slab as a function of the applied external potential and compare the resulting value of the double-layer capacitance with previous computations and experiments.

The DFT calculations for these two benchmarks and the following application are performed with VASP and the VASPsol module using the PAW formalism describing the electron-ion interactions and the PBE approximation for the exchange-correlation functional.^{17,33,36} The Brillouin-zone integration employs an automatic mesh with 50 *k*-points per inverse Ångstrom with only one *k*-point in the direction perpendicular to the slab.

We consider an electrolyte that consists of an aqueous solution of monovalent anions and cations of 1M concentration. At room temperature, this electrolyte has a relative permittivity of *ϵ*_{b} = 78.4 and a Debye length of *λ*_{D} = 3.04 Å. The parameters for the shape function, $\zeta (n(r\u2192))$, are taken from Refs. 15 and 28 to be *n*_{c} = 0.0025 Å^{−1} and *σ* = 0.6.

To properly resolve the interfacial region between the solute and the implicit electrolyte region and to obtain accurate values for the cavitation energy requires a high plane-wave basis set cutoff energy of 1000 eV. However, in our validation calculations and applications, we observe that the contribution of the cavitation energy to the solvation energies is negligibly small. We, therefore, set the effective surface tension parameter *τ* = 0 for all following calculations. This also removes the requirement for such a high cutoff energy, and we find that a cutoff energy of 600 eV results in converged surface energies and PZC.

To change the applied potential, we adjust the electron count and then determine the corresponding potential from the shift in electrostatic potential, as reflected in the shift in the Fermi level. This takes advantage of the fact that the electrostatic potential goes to zero in the electrolyte region for the solution of the Poisson-Boltzmann equation, providing a reference for the electrostatic potential.

Figure 3 compares the PZC for Cu, Ag, and Au surface facets calculated with the implicit electrolyte model as implemented in VASPsol with experimental data. The electrode PZC is defined as the electrostatic potential of a neutral metal electrode and is given by the Fermi energy measured relative to the reference potential. A natural choice of reference for the implicit electrolyte model is the electrostatic potential in the bulk of the electrolyte, which is zero in any solution of the Poisson-Boltzmann equation. In experiments, the reference is frequently chosen as the standard-hydrogen electrode (SHE). The dashed line in Fig. 3 presents the fit of $Upred=Uexp+\Delta USHEpred$ to obtain the shift, $\Delta USHEpred$, between the experimental and the computed PZC, i.e., the potential of the SHE for our electrolyte model. The resulting $\Delta USHEpred=4.6$ V compares well with the previously reported computed value of 4.7 V for a similar solvation model.^{28}

Figure 4 shows the calculated surface charge density, *σ*, of a Pt(111) surface as a function of applied electrostatic potential, *U*. We find that the PZC for the Pt(111) electrode is 0.85 V, in good agreement with previous computational studies.^{26,28} However, this value is at the high end of the experimentally reported results.^{35,37} This discrepancy is likely due to adsorption on Pt(111). The work by Sekong and Groß shows that Pt(111) is covered by hydrogen at potentials below 0.5 V, followed by the so-called double layer region between 0.5 and 0.75 V with no adsorbates present and OH adsorption at more positive potentials.^{25} The presence of adsorbates in the experiments alters the PZC and precludes a direct comparison with our calculations.

The linear slope of *σ*(*U*) is the double-layer capacitance. Fitting a quadratic polynomial to the data and evaluating the slope at the PZC yields a double-layer capacitance for the Pt(111) surface at 1M concentration of 14 *μ*F/cm^{2}. This agrees perfectly with the previously reported computed result^{28} and is close to the experimental value of 20 *μ*F/cm^{2}.^{37}

Our implementation of the electrolyte model neglects the finite volume of the ions in the electrolyte. Hence, at larger applied potentials, the model may predict unphysically large ion concentrations to screen the surface charge. To understand how this approximation limits the applied potential, we perform calculations on a Cu(111) slab in an electrolyte with 2M concentration for a range of applied potentials from −3 to +3 V.

Figure 5 shows the extrema of the net ion concentration in the double-layer region as a function of the applied potential (written by VASPsol into the file RHOION). The negative and positive concentrations correspond to excess cations and anions near the surface, respectively. For a large range of potentials from −2 to +3 V, the net ion concentration depends nearly linear on the applied potential. At negative potentials below −2 V, we observe that electrons leak from the slab into the electrolyte region, rendering the results of the electrolyte model unphysical. Hence, care should be taken when utilizing the model at large negative potentials.

The inset of Fig. 5 illustrates the average net ion concentration along the direction perpendicular to the Cu(111) surface for an applied potential of −1 V. The ion concentration decays exponentially away from the solid-liquid interface and approaches zero in the bulk electrolyte. We selected an electrolyte region of 30 Å to illustrate the need to converge the results with respect to the size of the electrolyte domain. Here, the size is not quite sufficient to reach a negligible ion concentration. While this does not affect the considered extrema of the ion concentration, care must be taken to ensure convergence of the results with the size of the electrolyte region. As general guidance, an electrolyte region >10*λ* should suffice.^{38}

At a large applied potential of +3 V, we find a net ion concentration of 0.2M, which is 10% of the electrolyte concentration, far from the 50% excess that would introduce a 1% error, as discussed following Eq. (14). Here, the concentrations of anions and cations are 2.1 and 1.9M, respectively, and the error introduced by the linear approximation to the sinh function is less than 0.1%. Thus, even for a considerable potential of +3 V, the linear electrolyte model does not result in unphysically large ion concentrations.

Finally, we validate the accuracy of the implemented analytic force evaluation. As an example, we have chosen the adsorption of a Na atom on the Au(111) surface. The surface charge was set to 0.4 *e*^{−} for the symmetric p(2 × 2) surface, corresponding to an electrochemical potential of about −1.9 V vs the SHE, while the neutral Na@Au(111) surface corresponds to a potential of −2.6 V. Figure 6 shows that the agreement between the analytical gradient and the minimum for the energy is excellent, provided a sufficiently high cutoff energy of 600 eV is used.

## V. CRYSTAL SHAPE CONTROL BY APPLIED POTENTIAL

The implicit electrolyte model, VASPsol, enables the study of electrode surfaces in realistic environments and under various conditions, such as of an electrode immersed in an electrolyte with an applied external potential. To illustrate the utility of the solvation model, we determine how an applied potential changes the equilibrium shape of a Cu crystal in a 1M electrolyte. The crystal shape is controlled by the surface energies. We show that the surface energies are sensitive to the applied potential and that each facet follows different trends, leading to opportunities for the practical control of the shape of nanocrystals.

Using the MPInterfaces framework,^{40} we construct slabs of minimum 10 Å thickness for the (111), (100), (110), (210), (221), (311), and (331) facets of Cu. The different facets for Cu are chosen based on the crystal facets included in the MaterialsProject database.^{39} We calculate the surface energy of each slab in vacuum and electrolyte by relaxing the top and bottom layers and keeping the middle layers of atoms fixed to their bulk positions. We employ a vacuum spacing of 30 Å to ensure that the electrostatic interactions between periodic images of the slabs are negligibly small. From the surface energies, we determine the resulting shape using the Wulff construction as implemented in the pymatgen framework.^{39}

Figure 7 compares the crystal shape of Cu in vacuum obtained in our calculations with that in MaterialsProject. The MaterialsProject calculations used a cutoff energy of 400 eV and a vacuum spacing of 10 Å , both of which are lower than the cutoff energy of 600 eV and the vacuum spacing of 30 Å used in this work. We find that our predicted surface energies agree with those in MaterialsProject within 2% for all the facets. However, we note that even such small deviations in surface energies matter for the details of the crystal shape, where the high index (310) and (221) facets do not show up in our calculation of the equilibrium crystal shape of Cu. These small but observable deviations demonstrate that the crystal shape is sensitive to small variations in the surface energies.

For the calculation of the Cu surface energies in the 1M electrolyte, we vary the number of electrons in the slabs and determine the change in the corresponding applied electrostatic potential, as described in Sec. III. We calculate the surface energies for each of the above facets in vacuum over a range of applied potentials of approximately ±1 eV around the corresponding PZC. We then perform a spline interpolation to obtain the surface energies as a function of the applied potential. This transformation enables the comparison of the surface energies for a given applied potential and is made possible by the absolute reference potential provided by the solvation model. From the results, we obtain the equilibrium crystal shape as a function of applied potentials by the Wulff construction.

For copper, the (111) facets exhibit the lowest surface energy over the whole potential range. Figure 8 shows how the ratios of the surface energies of the copper (hkl) facets relative to the (111) facet vary as a function of applied potential. As expected, the surface energies of all facets systematically increase with potential; however, their ratios relative to the (111) facet decrease. Due to this, the crystal shape changes as a function of applied potential.

Figure 9 illustrates the change in crystal shape with increasing applied potential. At negative potentials, the (111) facets dominate and lead to an octahedral shape. With increasing potential, the ratio of the (100) to (111) surface energies decreases sufficiently to lead to a change in shape to a truncated octahedron. The transition occurs around a potential of −0.56 V for Cu. With a further increase in potential, the area of the (100) facets increases. Even though the ratios of the surface energies of other facets are also lowered, the change is not sufficient for other facets to show up in the equilibrium crystal shape.

The prediction for the changes in the equilibrium shape of the Cu crystal as a function of applied electric potential demonstrates the utility of the VASPsol solvation model for computational electrochemistry simulations using density-functional theory. It furthermore provides the opportunity to design nanocrystals of specific shape through the application of an electrostatic potential during growth at a positive applied potential.

## VI. CONCLUSIONS

Solid/liquid interfaces between electrodes and electrolytes in electrochemical cells present a complex system that benefits from insights provided by computational studies to supplement and explain experimental observations. We developed and implemented a self-consistent computational framework that provides an efficient implicit description of electrode/electrolyte interfaces within density-functional theory through the solution of the linearized Poisson-Boltzmann equation. The electrolyte model is implemented in the widely used DFT code VASP, and the implementation is made available as a free and open-source module, VASPsol, hosted on GitHub at https://github.com/henniggroup/VASPsol. We showed that this model enables *ab initio* electrochemistry studies at the DFT level. We validated the model by comparing the calculated electrode potentials of zero charge for various metal electrodes with experimental values and the calculated double-layer capacitance of Pt(111) in a 1M solution to previous calculations and experiments. We tested the validity of the linear approximation at large applied potentials and the accuracy of the analytic expressions for the forces. To illustrate the usefulness of the model and the implementation, we apply the model to determine the equilibrium shape of the Cu crystal as a function of applied electrostatic potential. We predict a change in shape from an octahedron to a truncated octahedron with increasing potential. These calculations illustrate the utility of the implicit VASPsol model for computational electrochemistry simulations.

## ACKNOWLEDGMENTS

This work was supported by the National Science Foundation under Award Nos. DMR-1542776, OAC-1740251, and CHE-1665305. This research used computational resources of the Texas Advanced Computing Center under Contract No. TG-DMR050028N and of the University of Florida Research Computing Center. Part of this research was performed while the authors were visiting the Institute for Pure and Applied Mathematics (IPAM), which is supported by the National Science Foundation under Grant No. DMS-1440415. S.N.S. thanks the SYSPROD project and AXELERA Pôle de Compétitivité for support (PSMN Data Center).