Grand canonical electronic density-functional theory: Algorithms and applications to electrochemistry

Density-functional theory (DFT) enables theoretical elucidation of reaction mechanisms at complex catalyst surfaces, making it now possible to design efficient heterogeneous catalysts for various industrial applications from first principles, for example, for high-temperature gas-phase transformation of hydrocarbons to a variety of valuable chemical products.^1,2 The extension of this predictive power to electrocatalysis would be highly valuable for an even broader class of technological problems, including a cornerstone of future technology for renewable energy: converting solar energy to chemical fuels by electrochemical water splitting and carbon dioxide reduction.³ Accurately describing electrochemical phenomena, however, presents two additional challenges.

First, the electrolyte, typically consisting of ions in a liquid solvent, strongly affects the energetics of structures and reactions at the interface. Treating liquids directly in DFT requires expensive molecular dynamics to sample the thermodynamic phase space of atomic configurations. Historically, a number of continuum solvation models that empirically capture liquid effects have enabled theoretical design of liquid-phase catalysts.^4,5 More recently, empirical solvation models suitable for solid-liquid interfaces,^6–8 joint density-functional theory (JDFT) for efficiently treating liquids with an atomic-scale structure,⁹ and minimally empirical solvation models derived from JDFT^10,11 have made great strides towards reliable yet efficient treatment of electrochemical systems.

Second, electrons can flow in and out of the electrode as electrochemical reactions proceed. Changes in electronic charge of electrode surfaces and adsorbates can be especially important because the electrolyte stabilizes charged configurations with a counter charge from the ionic response. For example, reduction of formic acid on platinum at experimentally relevant potentials is dominated by formate ions rather than neutral molecules at the surface.¹² Proton adsorption on stepped and polycrystalline surfaces involves displacing oxidatively adsorbed water at relevant potentials, resulting in non-integer charge transfers and an anomalous pH dependence deviating from the Nernst equation.¹³ 

Accounting for the electrolyte response using our solvation models,^8,11 and adjusting the electron number to match experimentally relevant electrode potentials, realistic predictions of electrochemical reaction mechanisms have now become possible.¹⁴ In particular, application of this methodology to the reduction of CO on Cu(111) predicts onset potentials for methane and ethene formation with 0.05 V accuracy in comparison to experiment, for a wide range of pH varying from 1 to 12.¹⁵ However, conventional DFT software and algorithms are optimized for solving the quantum-mechanical problem at a fixed electron number, requiring extra work (both manual and computational) to calculate properties for a specified electrode potential.

Electric potentials and fields play an important role in fields besides electrochemistry. Density-functional theory approaches accounting for electric potential have been developed in special cases for field emission from metal surfaces using a jellium model¹⁶ and for calculating capacitance in metal-insulator-metal^17,18 and carbon nanotube systems.¹⁹ Calculating non-equilibrium transport of electrons in nanoscale systems also requires accounting for potential difference between reservoirs in a DFT calculation.²⁰ First-principles molecular dynamics approaches have been developed to emulate fixed potential using fluctuating numbers of electrons between time steps.²¹ However, in all these cases, each involved self-consistent DFT calculation contains a fixed number of electrons and is carried out using a conventional canonical-ensemble algorithm.

This paper introduces algorithms for grand canonical DFT, where electron number adjusts automatically to target a specified electron chemical potential (related to electrode potential), thereby enabling efficient and intuitive first-principles treatment of electrochemical phenomena. Section II summarizes the theoretical background of first-principles electrochemistry using JDFT and continuum solvation models and sets up the fundamental basis of grand-canonical DFT. Then, Section III introduces the modifications necessary to make two distinct classes of DFT algorithms, the self-consistent field (GC-SCF) method and the variational free energy minimization using auxiliary Hamiltonians (GC-AuxH), directly converge the grand free energy of electrons at fixed potential. Sections IV B and IV C establish the algorithm parameter(s) that optimize the iterative convergence of the GC-SCF and GC-AuxH methods, respectively, while Section IV D compares the performance of these algorithms for a number of prototypical electrochemical systems. Finally, Section IV E demonstrates the utility of grand canonical DFT by solving an electrochemical mystery: the identity of the second voltammetric peak in the under-potential deposition (UPD) of copper on platinum in chloride-containing electrolytes.

The exact Helmholtz free energy A of a system of interacting electrons in an external potential V(r) at a finite temperature T satisfies the Hohenberg-Kohn-Mermin variational theorem,^22,23

where $AHKMT$ is a universal functional that depends only on the electron density n(r) (and temperature) and not on the external potential. However, constructing approximations for this unknown universal functional that accurately capture the energies and geometries of chemical bonds in terms of the density alone is extremely challenging, partly because the quantum mechanics of the electrons is completely implicit in A_HKM[n] (dropping the T labels here onward for notational simplicity; all the functionals below depend on temperature).

Most practical approximations in electronic density-functional theory follow the Kohn-Sham approach²⁴ that includes the exact free energy of a non-interacting system of electrons with the same density n(r). The universal functional is typically split as

where A_ni[n] is the non-interacting free energy (which we describe in detail below), the second “Hartree” term E_H[n] is the mean-field Coulomb interaction between electrons (using atomic units $e,me,ℏ,kB=1$ throughout), and the final “exchange-correlation” term E_XC[n] captures the remainder which is not known and must hence be approximated. The exact free energy of non-interacting electrons is

which includes the kinetic energy and entropy contributions of a set of orthonormal single-particle orbitals $ψi(𝐫)$ with occupation factors $fi∈[0,1]$⁠. Above, the single-particle entropy function is $S(f)=−flogf−(1−f)log(1−f)$⁠. Note that we let the orbital index i include spin degrees of freedom as well and therefore do not introduce factors of 2 for spin degeneracy. For notational convenience, we also let i include Bloch wave-vectors in the Brillouin zone (in addition to spin and band indices) for periodic systems.

The minimization in (3) is constrained such that the density of the non-interacting system

the density of the real interacting system. Performing this minimization with Lagrange multipliers V_KS(r) (the Kohn-Sham potential) to enforce the density constraint and $εi$ (Kohn-Sham eigenvalues) to enforce orbital normalization constraints results in a set of single-particle Schrödinger-like equations

for stationarity with respect to $ψi(𝐫)$⁠, and the Fermi occupation condition

for stationarity with respect to f_i. Here, the electron chemical potential μ appears as a Lagrange multiplier to enforce the electron number constraint and is chosen so that $∑fi=N$⁠, the number of electrons in the system.

Optimizing the total free energy functional (1) with A_HKM[n] implemented by (2, 3) then yields the stationarity condition with respect to electron density,

Conventional density-functional theory calculations then amount to self-consistently solving the Kohn-Sham equation (5) along with (7), coupled via the electron density constraint (4).

The Kohn-Sham potential is arbitrary up to an overall additive constant: changing this constant introduces a rigid shift in the eigenvalues $ϵi$ and the electron chemical potential μ, but does not affect the occupations f_i, electron density, or free energy. For finite systems (of any charge) and for neutral systems that are infinitely periodic in one or two directions, this arbitrariness can be eliminated by requiring that the potential vanishes infinitely far away from the system, giving meaning to the absolute values of μ and $ϵi$ as being referenced to “zero at infinity.” However, for materials that are infinite in all three dimensions, such as periodic crystalline solids, there is no analogous natural choice for the zero of potential, making the absolute reference for μ and $ϵi$ meaningless.

More importantly, for systems that are periodic in one, two, or all three directions, the net charge per unit cell must be zero, otherwise the energy per unit cell becomes infinite. In terms of a finite system size L and then taking the limit $L→∞$⁠, the energy per unit cell of a system with net charge per unit cell diverges $∝lnL$ for one periodic direction, $∝L$ for two periodic directions, and $∝L2$ for three periodic directions. Therefore, the number of electrons per unit cell is physically constrained to keep the unit cell neutral in systems with any periodicity, and only finite systems (like molecules and ions) have number of electrons as a degree of freedom.

For describing electrochemical systems and electrocatalytic mechanisms, we are typically interested in adsorbed species in a solid-electrolyte interface which exchange electrons with the solid (electrode). In these systems, the solid surface, which we would describe in a density-functional calculation as a slab periodic in two directions, does have a net charge per unit cell that depends on the electrode potential. In contrast to the discussion at the end of Sec. II A, this is physically possible (i.e., has a finite energy) because the electrolyte contains mobile ions that respond by locally increasing the concentration of ions of opposite charge near the surface, thereby neutralizing the unit cell. (The charge per unit areas of the electrode and electrolyte are equal and opposite.)

Next to an electrolyte, the charge of a partially periodic system (one-dimensional “wires” or two-dimensional “slabs”) is no longer constrained, allowing the number of electrons to vary. Now, the absolute reference for the electron chemical potential μ does become physically meaningful, and μ now corresponds to the electrode potential that controls the number of electrons in the electrode.

However, treating electrochemical systems using electronic density-functional theory alone is extremely challenging for a variety of reasons. First, treating liquids requires a statistical average over a large number of atomic configurations to integrate over thermodynamic phase space. For this, techniques such as molecular dynamics typically require calculations of at least 10⁴–10⁵ atomic configurations (instead of just one for a solid). Second, such calculations require a large number of liquid molecules to minimize finite size errors in the molecular dynamics. For example, for electrolytes with a realistic ionic concentration of 0.1 M, there is on average one ion for a few hundred solvent molecules. Making statistical errors in the ion number manageable in such calculations therefore requires $>103$ solvent molecules with $>104$ electrons, contrasted with typical 10–100 atoms with 100–10³ electrons in the electrode slab + adsorbate of interest. Combined, these factors make density-functional molecular dynamics simulations of electrochemical systems prohibitively expensive computationally, in additional to being difficult to set up and analyze.

A viable alternative to the above direct approach is to employ joint density-functional theory⁹ (JDFT), a variational theorem akin to the Hohenberg-Kohn theorem that makes it possible to describe the free energy of a solvated system in terms of the electron density n(r) for the solute and in terms of a set of nuclear densities ${Nα(𝐫)}$ (where $α$ indexes nuclear species) of the solvent (or electrolyte). Specifically, the equilibrium free energy of the combined solute and solvent systems minimizes

where V(r) is the external electron potential, $Vα(𝐫)$ is the external potential on the liquid nuclei, and A_JDFT is a universal functional independent of these external potentials. Separating out the Hohenberg-Kohn electronic density functional A_HK[n] for the solute, the total free energy is

In practice, the functional A_diel, much like A_HKM, is unknown and needs to be approximated. Importantly, the liquid is now described directly in terms of its average density rather than individual configurations, and the expensive quantum-mechanical theory of the electrons is restricted to the solute alone, thereby addressing both the sampling and system-size problems that make density-functional molecular dynamics prohibitively expensive.

Typically, we are interested in a situation where the external potentials on the liquid are zero, and the liquid only interacts with itself, and with the electrons and nuclei of the solute. In this situation, we can perform the optimization over liquid densities and define the implicit functional A_JDFT[n] = A_HKM[n] + A_diel[n], where we define $AX[n]≡min{Nα(𝐫)}A∼X[n(𝐫),{Nα(𝐫)}]$ for X = JDFT and X = diel. We will work with these reduced functionals below for simplicity.

The framework of joint density-functional theory encompasses an entire hierarchy of solvation theories. Further separating the solvation term A_diel into a classical density functional for the liquid^25–27 and an electron-liquid interaction functional unlocks the full potential of JDFT to describe an atomic-scale structure in the liquid without statistical sampling. Starting from “full-JDFT” and performing perturbation theory for linear-response of the liquid result in the non-empirical SaLSA solvation model¹⁰ which continues to capture the atomic-scale nonlocality in the liquid response and introduces no fit parameters for the electric response of the solvent.

At the simplest end of the JDFT hierarchy, there are continuum solvation models that neglect the nonlocality of the liquid response and replace it by that of an empirically determined dielectric cavity (optionally with Debye screening due to electrolytes). This includes our recent CANDLE solvation model¹¹ that builds on the stability of SaLSA for highly polar systems, and earlier solvation models suitable for molecules and less-polar systems such as our GLSSA13 model⁷ (or its equivalent, VASPsol⁸) and the comparable Self-Consistent Continuum Solvation (SCCS) model.^6,28 Even the traditional quantum-chemistry finite-system solvation models such as the PCM series⁵ and the SMx series⁴ can be mapped on to this class of solvation models.

For simplicity, we will work here with this simplest class of continuum solvation models. The theoretical considerations and algorithms in this work focus primarily on the electronic density-functional theory component, are largely agnostic to the internals of the solvation model, and therefore straightforwardly generalize up the JDFT hierarchy to the more detailed and complex solvation theories. Essentially, all the simple electron-density-based continuum solvation models^6–8,11 can be summarized abstractly as^7,10

Here, the first term is the electrostatic solute-solvent interaction energy given by the difference between the solvent-screened Coulomb interaction and the bare Coulomb interaction $K^$ of the total solute charge density, $ρel(𝐫)=n(𝐫)+ρnuc(𝐫)$⁠, where $ρnuc(𝐫)$ is the solute nuclear charge density. The screened Coulomb interaction term above is expressed in terms of the bare Coulomb interaction $K^$ and the solvent susceptibility operator $χ^$⁠, defined for the local-response models by

Here, $ϵb$ is the bulk dielectric constant of the solvent and $κ=4π∑iNiZi2/T$ is the inverse Debye screening length in vacuum of a set of ionic species of charge Z_i with concentrations N_i (the net inverse Debye screening length in the presence of the background dielectric is $κ/ϵb$⁠). The cavity shape function s(r) modulates the dielectric and ion response and varies from zero in the solute region (no solvent present) to unity in the solvent at full bulk density. Finally, the second term of (10), A_cav, empirically captures effects beyond mean-field electrostatics, such as the free energy of cavity formation in the liquid and dispersion interactions between the solute and the solvent.²⁹ Different solvation models at this level of the hierarchy only differ in the details of how s(r) is determined from the electron density (or even from atomic positions and fit radii for the PCM⁵ and SMx⁴ solvation models) and in the details of A_cav.

In practice, evaluating the first term of (10) requires the calculation of $ϕ(𝐫)=(K^−1−χ^)−1ρel(𝐫)$⁠, which is the net electrostatic potential including screening by the solvent (and electrolyte). The second part of the first term in (10) represents the electrostatic interactions between electrons and the nuclei in vacuum which cancels corresponding terms in the vacuum DFT functional (specifically the Hartree, Ewald, and long-range part of the local pseudopotential terms), so that all long-range terms that contribute to the net free energy are present in the first term of (10). Finally, substituting (11) and $K^−1=−∇2/(4π)$ into the definition of $ϕ(𝐫)$ results in the linearized Poisson-Boltzmann (modified Helmholtz) equation,

Without ionic screening (second term of (12)), the absolute reference for $ϕ(𝐫)$ is undetermined and does not contribute to the bound charge induced in the liquid, $χ^ϕ(𝐫)$ (given by the first term of (11) alone), because $∇(constant)=0$⁠. However, with ionic screening, a constant shift of $ϕ(𝐫)$ does affect the second terms of (11) and (12), producing a charge response in the liquid and making the absolute reference of the electrostatic potential $ϕ(𝐫)$ meaningful. In particular, integrating (12) over space yields

From (11), the left hand side is the negative of the total bound charge in the liquid, while the right hand side is the total charge of the solute system. Therefore, the continuum electrolyte automatically compensates for any net charge in the solute and makes the complete system neutral. As a side effect, the absolute reference of the electrostatic potential is meaningful and automatically corresponds to “zero at infinity.” This happens because the Green’s function of (12) is $exp(−κr/ϵb)/r$ in the bulk liquid with finite κ (instead of 1/r), causing $ϕ(𝐫)$ to exponentially decay to zero outside the solute region (where $ρel=0$⁠).

Reference 30 gives a detailed version of the above discussion including rigorous proofs of the absoluteness of the reference for $ϕ(𝐫)$⁠, and numerical details and algorithms for efficiently solving (12) with periodic boundary conditions in plane-wave basis DFT calculations. The result that the electrolyte neutralizes charges in the solute is true more generally, even for solvation models with nonlinear⁷ and/or nonlocal response.¹⁰ However, as Ref. 7 describes in detail, for the general nonlinear case, unit cell neutralization is not automatic and must be imposed using a Lagrange multiplier constraint; this Lagrange multiplier then fixes the absolute value of $ϕ(𝐫)$ such that it is zero at infinity.

The key conclusion of this discussion is that continuum electrolytes present two advantages. First, automatic charge neutralization implies that we can perform well-defined calculations with net charge per unit cell in the solute. Second, meaningful absolute reference values for potentials implies that the electron chemical potential μ (which determines the electron occupations in (6)) is also referenced to zero at infinity. Consequently, μ is related directly to the potential U of the electrode providing the electron reservoir in experiments. This potential is typically referenced to the standard hydrogen electrode (SHE), so that

where $μSHE$ is the absolute position of the standard hydrogen electrode relative to the vacuum level (i.e., the zero at infinity reference). Calculations can employ either the experimental estimate $μSHE=−4.44$ eV to relate the absolute and relative potential scales^31,32 or a theoretical calibration based on the calculated and measured potentials of zero charge of solvated metal surfaces.³⁰ The latter approach has the advantage of minimizing systematic errors in the solvation model since they cancel between the calculations used for calibration and prediction. For example, with the CANDLE solvation model we use below, the calibrated $μSHE=−4.66$ eV.¹¹ 

Using joint density-functional theory of continuum solvation models to treat electrolytes as described above, we can calculate the Helmholtz free energy and electrochemical potential (μ, and hence U) of specific microscopic configurations of adsorbates on electrode surfaces at a fixed number of electrons N in the solute subsystem (electrode + adsorbates). However, in electrochemical systems, N is an artificial constraint because electrons can freely exchange between the electrode and an external circuit. Instead, experiments set the electrode potential U, and hence the electron chemical potential μ, and N adjusts accordingly as a dependent variable. Thermodynamically, this corresponds to switching the electrons from the finite-temperature, fixed-number canonical ensemble to the finite-temperature, fixed-potential grand-canonical ensemble. Correspondingly, the relevant free energy minimized at equilibrium is the grand free energy $Φ=A−μN$⁠, instead of the Helmholtz free energy A.

The most straightforward approach to fixed-potential DFT for electrochemistry is to repeat conventional fixed-charge DFT calculations at various electron numbers N to reach a target chemical potential μ. Using a steepest descent/secant method to optimize N results in convergence of μ typically within 10 iterations.¹⁴ However this approach is inefficient since it requires multiple DFT calculations to calculate the grand free energy and charge of a single configuration, typically taking 3 times the time of a single fixed-charge calculation (not 10 times because subsequent calculations have a better starting point and converge quicker).¹⁴ 

A more efficient approach would be to directly optimize the Kohn-Sham functional in the grand-canonical ensemble. For a solvated system, this corresponds to a small modification of the minimization problem (8) with A_JDFT given by (9) and with A_HKM evaluated using the Kohn-Sham approach ((2) and (3)). The only differences are that the Lagrange multiplier term $−μ(∑ifi−N)$ that enforced the electron number constraint is replaced by $−μ∑ifi$ which implements the Legendre transform of the Helmholtz free energy to the grand free energy, and that the fixed-N constraint is removed. The electron occupation factors are still Fermi functions (6), but μ is specified as an input (instead of being adjusted to match a specified N). Operationally, this amounts again to solving the Kohn-Sham eigenvalue problem (5) self-consistently with

where there is now a single extra contribution to the potential due to the solvent (electrolyte) and with the aforementioned changes to the Lagrange multipliers and constraints.

This conceptually simple modification, however, presents numerical challenges to the algorithms commonly used for solving the Kohn-Sham problem, such as the self-consistent field (SCF) method, which solves the Kohn-Sham eigenvalue problem from an input density (or Kohn-Sham potential), adjusting this density (or potential) iteratively until self-consistency is achieved. A common instability for metallic systems in the SCF method is “charge sloshing”: the electron density oscillates spatially between iterations instead of converging. Switching to the grand-canonical ensemble can significantly exacerbate this problem: electrons can now additionally slosh between the system and the electron reservoir. Here, we present modifications to the SCF method and an alternate algorithm that allow reliable and efficient convergence for grand-canonical Kohn-Sham DFT.

Given electron density $nin(i)$ (𝐫), the Kohn-Sham equation (5) with potential V_KS given by (15) defines orbitals ${ψj}$ and eigenvalues ${εj}$⁠. At a given electron chemical potential μ, these in turn define the occupations { f_j} given by (6) and a new electron density $nout(i)$(𝐫) given by (4). The Self-Consistent Field (SCF) method attempts to find n⁽ⁱ⁾(r) such that $nout(i)(𝐫)=nin(i)(𝐫)$⁠.

There are two algorithmic ingredients to this method: solution of the Kohn-Sham eigenvalue equations and optimization of the electron density. The eigenvalue equations remain unchanged between conventional and fixed-potential calculations. To solve these, we use the standard Davidson algorithm.³³ The difficulty in fixed-potential calculations arises in the electron density optimization, which we discuss below.

A robust and frequently used algorithm for charge-density optimization is Pulay mixing³⁴ with Kerker preconditioning.³⁵ Briefly, Pulay mixing assumes that the residual $R[nin(𝐫)]≡nout(𝐫)−nin(𝐫)$ is approximately linear in the input electron density n_in(r) and calculates the optimum input electron density as a linear combination of previous iterations,

Minimization of the norm of the corresponding residual

with the constraint $∑iαi=1$ yields a set of linear equations determining the coefficients $αi$⁠, where $M^$ is the metric for defining the norm of the residual. Finally, the next input density is obtained by mixing the optimum input density with its corresponding output density,

where $K^$ is the Kerker preconditioning operator.

The metric $M^$ and preconditioner $K^$ serve to balance the influence of different components of the electron density to the optimization procedure. These are usually defined in reciprocal space, expanding $n(𝐫)=∑𝐆n∼(𝐆)ei𝐆⋅𝐫$⁠, where G are reciprocal lattice vectors. In reciprocal space, the Hartree potential in (15) takes the form $V∼H(𝐆)=n∼(𝐆)4π/G2$⁠, causing small G (long-wavelength) variations of the electron density to produce larger changes in the Kohn-Sham potential than large G ones. Uncompensated, this makes the electron density optimization unstable against long-wavelength perturbations (the charge-sloshing instability). The Kerker preconditioner

mitigates this problem by suppressing the contribution of the problematic small G components in determining the next input density. The metric

enhances the contribution of small G components in the residual norm, thereby prioritizing the optimization of these components when determining $ninopt$⁠. The wavevectors q_K and q_M, which control the balance between short and long-wavelength components, and the prefactor A, which sets the maximum fraction of n_out that contributes to the next n_in, can be adjusted to optimize the convergence of the SCF method. See Ref. 36 for a detailed discussion of the Pulay-Kerker SCF approach and its performance for conventional fixed electron number (canonical) DFT calculations.

The G = 0 component of the electron density, which equals $N/Ω$⁠, where Ω is the unit cell volume, remains fixed in canonical DFT calculations and is therefore excluded from the metric and the preconditioner. This is no longer true in fixed-potential DFT, where the number of electrons changes. With the above prescriptions, the Kerker preconditioner (19) $→0$ as $G→0$ which will prevent the electron number from changing between SCF iterations. Likewise, the $G→0$ divergence of the Pulay metric (20) causes the residual norm to become undefined when the electron number changes between iterations.

In order to generalize the Pulay-Kerker SCF approach for fixed-potential DFT, we therefore need to fix the $G→0$ behavior of both the preconditioner and the metric. The need for the preconditioner and the metric arose from the reciprocal space Coulomb operator $4π/G2$ in the Hartree potential. In a uniform electrolyte with Debye screening, the reciprocal space Coulomb operator instead takes the form $4π/(ϵbG2+κ2)$ (from (12) in reciprocal space with s = 1). Since the electrolyte is responsible for fixing the indeterminacy of the G = 0 component of the potential, a reasonable ansatz for extending Pulay-Kerker to fixed potential calculations is replacing G² with $G2+qκ2$⁠, where

In Section IV B, we show that setting

and

indeed makes the grand canonical self-consistent field (GC-SCF) method function efficiently, with optimum convergence for $qκ$ given by (21).

An alternate approach to solving the Kohn-Sham equations is to directly minimize the total (free-)energy functional (1), with A_HKM given by (2, 3), in terms of the Kohn-Sham orbitals as independent variables. For joint density-functional theory, this amounts to

where n(r) is now derived from ${ψi}$ and {f_i} as given by (4). In the above minimization, the orbitals ${ψi}$ must be orthonormal, the occupation factors must satisfy $0≤fi≤1$⁠, and optionally, $∑ifi=N$ for the canonical fixed electron number case.

For insulators at T = 0, the occupations f_i are known in advance. The constrained optimization over orthonormal ${ψi}$ is most efficiently carried out using a preconditioned conjugate-gradients (CG) algorithm on unconstrained orbitals using the analytically continued energy functional approach.³⁷ Briefly, this approach evaluates the energy functional (24) on a set of orthonormal orbitals, which are a functional of the unconstrained orbitals used for minimization (see Ref. 37 for further details).

The general case of metallic systems and/or finite T additionally requires the optimization of { f_i}, which is challenging for non-linear optimization algorithms because of the inequality constraints $0≤fi≤1$⁠. One possibility is to update the fillings from the Kohn-Sham eigenvalues using (6) after every few steps of the CG algorithm,³⁸ but this hinders the convergence of CG because the functional effectively changes each time the fillings are altered. The “Ensemble DFT” approach³⁹ rectifies this convergence issue by optimizing the occupation factors at fixed orbitals in an inner loop and performing the optimization of orbitals in an outer loop using the CG method, but this increases the computational cost compared to the case of the insulators. The SCF approach is typically much more computational efficient than these variants of the direct variational minimization algorithm with variable occupations.³⁶ 

An alternate strategy for direct variational minimization with variable occupations is to introduce an auxiliary subspace Hamiltonian matrix H_aux as an independent variable⁴⁰ and setting the occupations $fi=f(ηi)$ in terms of the eigenvalues ${ηi}$ of H_aux instead of the Kohn-Sham eigenvalues ${εi}$⁠. (Here, $f(η)$ is the Fermi function given by (6), and the electron chemical potential μ is chosen so as to satisfy the electron number constraint $∑if(ηi)=N$⁠.) This eliminates the problematic inequality constraints on the occupations, and upon minimization, the auxiliary subspace Hamiltonian $Hauxij$ approaches the true subspace Hamiltonian, $Hsubij=⟨ψi|H^KS|ψj⟩$⁠, where $H^KS=−∇i2/2+VKS(𝐫)$ is the Kohn-Sham Hamiltonian. By choosing the undetermined unitary rotations of the orbitals ${ψi}$ to diagonalize H_aux, Ref. 40 further shows that the gradient of the free energy with respect to the independent variables simplifies to

and

These gradients are used to perform line minimization along a search direction in the space of independent variables. With every update of the independent variables, the orbitals are re-orthonormalized and the unitary rotations of the orbitals are updated to keep the auxiliary Hamiltonian diagonal.

In the preconditioned CG algorithm, the next search direction is obtained as a linear combination of the current search direction and the preconditioned gradients given by

and

where $T^inv$ and K are preconditioners. The role of preconditioning is to balance the weight of different directions in the minimization space in the explored search directions, and ideally the preconditioner equals the inverse of the Hessian (which is difficult to compute exactly). For the orbital directions, the standard preconditioner $T^inv$ resembles the inverse of the dominant kinetic energy operator in the Kohn-Sham Hamiltonian.³⁷ For the auxiliary Hamiltonian direction, the preconditioner removes the Fermi function derivatives and finite difference factors from (26) in order to more equitably weight all components of H_aux. The preconditioning factor K controls the overall contribution of the H_aux components relative to that of the orbital components and is adjusted to achieve optimum convergence.

In this auxiliary Hamiltonian (AuxH) approach, the orbitals and occupations are continuously and simultaneously optimized to minimize the total free energy, resulting in better convergence and computational efficiency comparable to the fixed-occupations insulator case, and competitive with the SCF method even for metallic systems. See Ref. 40 for further details on the algorithms and performance comparisons for conventional fixed electron number calculations.

Now, for fixed potential calculations, we set the occupations to Fermi functions of the auxiliary Hamiltonian eigenvalues at a specified μ, instead of selecting μ based on the electron number constraint. Correspondingly, the second term of the auxiliary Hamiltonian gradient (26), which arises from this constraint, drops out, and the algorithm requires no further modification.

The convergence rate of this algorithm, however, is sensitive to the preconditioning factor K and we propose a modified heuristic to update K automatically and continuously. At the end of each line minimization, the derivative of the optimized free energy A_min with respect to K can be evaluated from the overlap between the auxiliary Hamiltonian gradient and search direction.⁴⁰ The line-minimized energy is optimum for $∂Amin/∂K=0$⁠. Therefore, if variation of A_min with respect to K is convex, we should increase K if we find $∂Amin/∂K<0$ and vice versa. However convexity is often lost if K is initialized at too high a value. Therefore, our heuristic tries to zero $∂Amin/∂K$ while limiting the contribution of the auxiliary Hamiltonian gradient. In particular, we update

at the end of each line minimization, where $fsat(x)≡x/1+x2$ to saturate the factor by which K can change in one iteration, g_aux is the overlap of the H_aux components of the gradient and preconditioned gradient, and g_tot is the total overlap of the gradient and preconditioned gradients (orbital + H_aux). Finally, we reset the conjugate-gradient algorithm (i.e., set the search direction to the negative of the preconditioned gradient direction) after K has increased or decreased by a factor greater than e². We do this because dynamically changing the preconditioner technically invalidates the strict orthogonality of the CG search direction with previous directions. Section IV C shows that this heuristic exceeds the convergence obtained with fixed K, while Section IV D shows that the grand-canonical auxiliary Hamiltonian (GC-AuxH) algorithm consistently outperforms GC-SCF for fixed-potential calculations.

We implement all algorithms and perform all calculations using the open-source plane-wave density-functional theory software, JDFTx.⁴¹ Below, we specify computational and convergence parameters in atomic units (distances in bohrs $a0≈0.529$ Å  and energies in hartrees $Eh≈27.2$ eV), but present any physically relevant properties in conventional units (Å, eV). All calculations in this work employ the PBE⁴² exchange-correlation functional with GBRV ultrasoft pseudopotentials⁴³ at a kinetic energy cutoff of 20 E_h for Kohn-Sham orbitals and 100 E_h for the charge density. The metal surface calculations use inversion-symmetric slabs of at least five layers, with at least 15 Å  vacuum separation and truncated Coulomb potentials⁴⁴ to minimize interactions with periodic images. For Brillouin zone integration, we use a Fermi smearing of 0.01 E_h and a Monkhorst-Pack k-point mesh along the periodic directions with the number of k-points chosen such that the effective supercell is larger than 30 Å  in each direction. We use the CANDLE solvation model to describe the effect of liquid water and Debye screening due to 1M electrolyte, which we showed recently to most accurately capture the solvation of highly charged negative and positive solutes.¹¹ We emphasize that the methods and algorithms described above do not rely on specific choices for the pseudopotential, exchange-correlation functional, k-mesh, or solvation model; we keep these computational parameters constant here for consistency.

The self-consistent field (GC-SCF) algorithm summarized in Section III A depends on several parameters that control its iterative convergence. All these parameters are common to the conventional fixed-charge version (SCF) and the fixed-potential variant (GC-SCF) introduced here, except for the low-frequency cutoff wavevector $qκ$ that is necessary to allow the net electron number to change in the fixed-potential case. Figure 1 compares the dependence of GC-SCF convergence on $qκ$ for a prototypical calculation of an electrochemical system: a Cu(111) surface treated using a five-layer inversion-symmetric slab surrounded by 1M aqueous non-adsorbing electrolyte treated using the CANDLE solvation model,¹¹ at a fixed potential of $μ=−0.208Eh$ (1 V SHE). (The remaining GC-SCF parameters are set to their default values which we discuss below.) The best convergence is obtained with $qκ=0.17$ which corresponds to the Debye screening length of the electrolyte (21). For small $qκ$⁠, including the conventional case of $qκ=0$⁠, the number of electrons does not respond sufficiently quickly, stalling at about 0.1 electrons from the converged value, correspondingly with the free energy stalling at about 0.1 E_h (⁠$≈2.7$ eV) away from the converged value. Larger $qκ$ causes the electron number to change too rapidly, hindering convergence and eventually leading to a divergence as seen for the case of $qκ=0.7$⁠. An issue remains in the convergence independent of $qκ$⁠: after initial convergence, the free energy oscillates at the $10−6Eh$ level, while the electron number oscillates at the 10⁻³ level.

Keeping $qκ$ at this optimum value given by (21), we next examine the dependence of GC-SCF convergence on the remaining algorithm parameters for the same example system. Figure 2 shows the dependence on the Kerker mixing wavevector q_K, which helps stabilize the GC-SCF algorithm against long-wavelength charge oscillations. Optimal convergence is obtained for the typical recommended value³⁶ of $0.8a0−1$ (⁠$≈1.5$ Å⁻¹). As expected, convergence is relatively insensitive to the exact choice of q_K, as long as q_K does not become so small that convergence is ruined by charge sloshing. Notice that the final convergence beyond the $10−6Eh$ and 10⁻³ electron level remains an issue that is not resolved for any choice of q_K.

Next, Figure 3 shows the variation of GC-SCF convergence with the wavevector q_M controlling the reciprocal-space metric used by the Pulay algorithm. The convergence is entirely insensitive to this choice, and we henceforth set $qM=qK=0.8a0−1$ (the recommended value³⁶). Again, the final convergence issue remains unaffected by the choice of q_M.

Finally, Figure 4 compares the dependence of GC-SCF convergence on the maximum Kerker mixing fraction A, which effectively controls what fraction of the new electron density is mixed into the current value. We find nominally the best convergence for A = 0.5, but the performance of other values is not much worse. Smaller values of A lead to greater stability initially, but marginally slower convergence later on, while larger values of A lead to greater oscillations initially, but faster convergence later on. Regardless, as before, good convergence is obtained until the free energy reaches the $10−6Eh$ level, but continues to oscillate at that level beyond that point.

This final convergence issue may not affect practical calculations where relevant energy differences are at the $10−3Eh$ level or higher. However, smooth exponential convergence to the final answer is desirable as this makes it easier to determine when the target accuracy has been reached. Unfortunately, no combination of GC-SCF parameters achieves uniformly smooth convergence for the fixed-potential case. On the other hand, with our $qκ$ modification, the GC-SCF algorithm at least converges to the $10−6Eh$ level independent of system size (Figure 5), with the number of cycles required for convergence remaining mostly unchanged with increasing number of Cu(111) layers, and with increasing thickness of the solvent region.

Note that the standard fixed-charge SCF method converges the energy smoothly and exponentially in most cases, including in our implementation in JDFTx. Our implementation of the GC-SCF method in JDFTx uses exactly the same code, except for the preconditioner and metric modifications (due to $qκ$⁠) presented here. Therefore we believe that the final convergence difficulty in GC-SCF is a property of the algorithm itself, rather than an implementation issue.

The grand-canonical auxiliary Hamiltonian (GC-AuxH) approach discussed in Section III B directly minimizes the total free energy of the system without assuming any models for physical properties of the system (such as dielectric response models that are built into the SCF mixing schemes). This algorithm contains a single parameter K, which weights the relative contributions of the Kohn-Sham orbital and subspace Hamiltonian degrees of freedom in the conjugate-gradients search direction for free energy minimization.

Figure 6 shows the dependence of iterative convergence of the GC-AuxH algorithm on this preconditioning parameter K for the same Cu(111) test problem considered above. If the preconditioning factor K is held fixed, the rate of convergence is sensitive to the choice of K, with the optimum choice being $K≈0.3$ for this system. With the preconditioner auto-adjusted using the heuristic given by (29), we find that indeed the convergence picks up from that of the sub-optimal K = 1 towards that of the optimal value. More importantly, we observe smooth exponential convergence of both the free energy and the electron number, in contrast to our experiences with the GC-SCF method.

Figure 7 further shows that this smooth convergence sustains with changing system size. In particular, the convergence is virtually unchanged with the thickness of the solvent regions, but slows down slightly with increasing number of copper layers in the surface slabs.

Having analyzed and optimized the convergence of the GC-SCF and GC-AuxH methods, we now compare the performance of these algorithms for a few different cases. In this comparison, we also include the present state of the art: the “Loop” method which uses a secant method to adjust the number of electrons in an outer loop to match the specified electron chemical potential.¹⁴ For a fair comparison, we use the fixed-charge SCF method in the inner loops, because it achieves the fastest convergence; this algorithm works equally well with the AuxH method in the inner loop, but is then marginally slower. In each test case, we start all three algorithms from the same starting point: the converged state of the corresponding neutral (fixed-charge) calculation. Different test cases effectively perturb the potential by different amounts, thereby testing the algorithms for a range of proximities between initial and final states. Also, we now compare the wall time between algorithms, because there is no straightforward correspondence between GC-SCF cycles and conjugate-gradient iterations of the GC-AuxH method. The relative wall-time performance of these fairly distinct algorithms will depend to an extent on details of code optimization for each. However, there is no perfect metric for comparing these algorithms and wall time suffices for a rough qualitative comparison.

First, Figure 8 compares the performance of the algorithms for the 5-layer Cu(111) slab used in all the tests so far. The spikes in the free energy seen in the Loop method are the points where the electron number changes in the outer loop and a new SCF convergence at fixed charge begins. Both the GC-SCF and GC-AuxH methods are quite competitive, cutting the time to convergence within $10−6Eh$ in half compared to the Loop method. Given the smooth convergence beyond $10−6Eh$ however, the GC-AuxH method is preferable over GC-SCF for fixed potential calculations. Note that in the fixed-charge case, when SCF converges smoothly, it often outperforms the AuxH method as mentioned above. The advantage of the AuxH and GC-AuxH methods is their stability on account of being variational methods: the free energy is guaranteed to decrease at every step. Consequently, the convergence difficulties of the GC-SCF method make the variationally stable GC-AuxH method relatively more attractive.

Next, we compare these algorithms for more complex text cases. Figure 9 compares the convergence for a five-layer Pt(111) slab fixed to a potential of $μ=−0.171Eh$ (0 V SHE). At this potential, the surface of Pt(111) charges negatively, and the CANDLE solvation model brings the cavity closer to the electrons to capture the more effective solvation of negative charges by liquid water. Additionally, the dielectric response of platinum is more complex than copper due to the partially filled d shell. Both of these factors make this system harder to converge than the previous test case, and therefore the convergence of the GC-AuxH method is no longer clearly a single exponential. Despite this, both the direct grand-canonical methods converge faster than the Loop method, with the GC-AuxH method eking out an advantage in final convergence as before.

Finally, Figure 10 compares the convergence for chloride anions adsorbed at one-third monolayer coverage, in a $3×3$ supercell of a five-layer Pt(111) slab. Despite the increased complexity, the direct grand-canonical methods exhibit the best convergence, edging out the Loop method by a factor of four in wall time now, again with the smoothest convergence for the GC-AuxH method. Due to the systematic convergence advantage of the GC-AuxH method, we use it for all remaining calculations in this work and recommend it as the default general purpose algorithm for converging electrochemical calculations.

In Sec. II and all calculations so far, we used Fermi smearing where the electron occupations are given by the Fermi function (6). Practical k-point meshes typically require the use of a temperature T substantially higher than room temperature (we used 0.01 E_h which is approximately ten times higher), which could result in inaccurate free energies. Such errors can be reduced substantially by changing the functional form of the occupations and the electronic entropy, with the caveat that the smearing width T no longer corresponds to an electron temperature. Common modifications include Gaussian smearing, where the Fermi functions are replaced with error functions, and Cold smearing,⁴⁵ where the functional form is chosen to cancel the lowest order variation of the free energy with T (see Ref. 45 for details).

Figure 11 compares the performance of the preferred GC-AuxH method for various smearing methods. The insets show the variation of the converged free energy and electron number with smearing width T. Gaussian smearing reduces the coefficient of the quadratic T dependence compared to Fermi smearing, while Cold smearing cancels the quadratic dependence altogether, by design. The variation of electron number with smearing width is also reduced by Cold smearing, but to a lesser extent. Notice that the use of Cold smearing marginally slows down the iterative convergence of the GC-AuxH method. However, using Cold smearing at a high width of 0.01 E_h is still faster than using any smearing method with the lower width 0.001 E_h that is close to room temperature (because of the far denser Brillouin zone sampling required for smaller widths). Therefore, it is still advantageous to use Cold smearing at elevated smearing widths, and so we use Cold smearing with a width of 0.01 E_h for the final demonstration below.

The application of sufficiently negative (reductive) potentials on an electrode immersed in a solution containing metal ions reduces those ions and results in bulk electro-deposition of metal on the surface. Additionally, for many pairs of metals, a single monolayer of one metal deposits on a surface of the other at an under potential, that is, at a potential less favorable than for bulk deposition. This phenomenon of under-potential deposition (UPD) has several technological applications since it enables precise synthesis of heterogeneous metal interfaces. It also serves as an archetype for fundamental studies of electrochemical processes (see Ref. 46 for an extensive review), which makes it a perfect example for demonstrating our grand-canonical density-functional theory method.

The basic reason for underpotential deposition is that the heterogeneous binding between the two metals is stronger than the homogeneous binding of the depositing metal to itself. Indeed, metal pairs that exhibit underpotential deposition also display analogous phenomena in vapor adsorption.⁴⁷ However, the process in solution is far more complicated and highly sensitive to the composition of the solution because of competing adsorbates,⁴⁸ as well as to the structure of the electrode surface.⁴⁹ 

The UPD of copper on Pt(111) in the presence of chloride anions is particularly interesting and the subject of considerable debate in the literature. Voltammetry for this system⁵⁰ exhibits two well-separated under-potential peaks, as shown in the background of Figure 12. Certain LEED and in situ X-ray scattering studies of this system⁵¹ find evidence of a $2×2$ bilayer of copper and chloride ions co-adsorbed on the surface at potentials between the two peaks, suggesting that one peak corresponds to the formation of a partial layer, and the second peak to the formation of the full monolayer. In contrast, other studies^50,52,53 do not find this signature and propose that the additional peak arises from adsorption and desorption of chloride ions alone.

To address this debate, we perform grand-canonical density functional theory calculations of various configurations of copper and chlorine adsorbed on a 5-layer Pt(111) slab in $3×3$ and $2×2$ supercells. We determine the most stable configurations at each potential (⁠$μ$⁠) and the potentials at which transitions between configurations occur by comparing their grand free energies. The relevant grand free energy of a configuration α containing $NPtα$ platinum atoms in the slab with $NCuα$ copper and $NClα$ chlorine atoms adsorbed at the surface within the calculation cell is

normalized by the number of surface platinum atoms $NPtsurf$ in order to correctly compare energies of calculations in different supercells. (⁠$NPtsurf=2$ for the unit cell, 6 for the $3×3$ supercell, and 8 for the $2×2$ supercell, accounting for the top and bottom surfaces in the inversion-symmetric setup.) Since no light atoms are present, we safely neglect changes in vibrational contributions to the free energy between adsorbate configurations.

Above, $Φα(μ)$ is the free energy of adsorbate configuration α calculated by fixed-potential DFT, which is grand canonical with respect to the electrons at chemical potential μ (related to the electrode potential by (14) as discussed at the end of Section II B). Then, (30) calculates the free energy $Φ∼α(μ)$ which is additionally grand canonical with respect to all relevant atoms with chemical potentials $μPt$⁠, $μCu$ and $μCl$⁠. Several conventions are possible in defining the electron-grand-canonical free energy Φ, depending on what electron number we subtract: change from neutral value, total electron number, or number of valence electrons in pseudopotential DFT calculations. The atom chemical potentials would then respectively correspond to neutral atoms, bare nuclei, or pseudo-nuclei (nuclei + core electrons in pseudopotential). The full grand canonical free energy $Φ∼$ does not depend on this choice. In our JDFTx implementation, we choose the last option above (number of valence electrons and correspondingly atom chemical potentials of the pseudo-nuclei).

The bulk of the platinum electrode sets the Pt chemical potential, $μPt=EPt(s)−μNPt(s)e$⁠, where E_Pt(s) is the DFT energy of a bulk fcc Pt calculation with a single atom in the unit cell, and $NPt(s)e$ is the number of valence electrons in that calculation. The second term here implements the electron counting convention discussed above. Next, copper ions in solution set $μCu$⁠, but directly calculating the free energy of such ions using solvation models is error-prone.¹¹ So, instead, we use the DFT calculated energy, E_Cu(s), of a bulk fcc Cu calculation (containing $NCu(s)e$ valence electrons) and relate it to the free energy of the ion via the experimentally determined standard reduction potential $UCu2+→Cu(s)=0.342$ V SHE.⁵⁴ This yields

where the second term accounts for the change from Cu(s) to Cu²⁺ ions, and the final term accounts for change in ionic concentration from the standard value of 1 mol/l to the current value of [Cu²⁺] (in mol/l). Similarly, chlorine ions in solution set $μCl$⁠, but to minimize DFT errors, we connect to the DFT calculated energy, E_Cl(at), of an isolated chlorine atom (containing $NCl(at)e$ valence electrons), via the experimentally determined atomization energy $ECl2→2Cl(at)=242.6$ kJ/mol,⁵⁵ gas-phase entropy $SCl2(g)=223.1$ J/mol K,⁵⁴ and reduction potential $UCl2(g)→2Cl−=1.358$ V SHE.⁵⁴ Specifically,

where the second term accounts for the change from atomic to gas-phase chlorine, the third term for the change to chloride ions, and the final term for the change in chloride ion concentration to [Cl⁻] (in mol/l).

Figure 12(a) shows the calculated grand free energies as a function of electrode potential for a number of Cu and Cl adsorbate configurations on the surface of Pt(111). At high potentials, the most stable (lowest free energy) configuration is 1/4 Cl coverage (Figure 12(f)), which transitions to a clean Pt surface (Figure 12(e)) at a potential of 0.55 V SHE. Upon further lowering the potential, the stable configuration transitions to a full monolayer of copper with 1/3 Cl coverage (Figure 12(c)) at a potential of 0.46 V SHE. Experimentally, the two voltammogram peaks are at approximately $(0.63±0.04)$ and $(0.51±0.02)$ V SHE, averaging over the forward and reverse direction sweeps. Therefore chlorine desorption and full-copper-monolayer formation are plausible explanations^50,52,53 of the two peaks, with our first-principles predictions reproducing well the peak spacing (0.09 eV versus 0.12 eV in experiment), and placing the absolute locations of the peaks to within 0.07 eV. Similar accuracy has been achieved in comparison to experiment for onset potentials and product selectivity in CO reduction on Cu(111)^15,56 and the oxygen evolution reaction on IrO2(110)⁵⁷ in concurrent work using exactly the same calculation protocol as here: fixed-potential DFT in JDFTx using the PBE exchange-correlation functional and the CANDLE solvation model. The partial 2 × 2 monolayer of copper (Figure 12(d)) proposed by others⁵¹ as the reason for the second peak is not the most stable configuration in our calculations at any potential, lying a significant 0.3 eV above the other phases at relevant potentials.

For comparison, Figure 12(b) shows the analogous results that would be obtained using only conventional vacuum calculations. In the above formalism, this corresponds to assuming $Φα(μ)≈Aα−μNeα$⁠, where $Aα$ is the Helmholtz energy from a neutral vacuum DFT calculation of configuration α containing $Neα$ valence electrons. This approximation results in a single transition directly from the Cl-covered Pt surface to the one with a copper monolayer, predicting a single voltammogram peak in disagreement with experiment. Accurate predictions for electrochemical systems therefore require treating charged configurations stabilized by the electrolyte at relevant electron potentials, now easily accomplished with the methods and algorithms introduced in this work.

This work introduces algorithms for directly converging DFT calculations in the grand-canonical ensemble of electrons, where the number of electrons adjusts to maintain the system at constant electron chemical potential, while ionic response in a continuum solvation model of electrolyte keeps the system neutral. We show that, with appropriate modifications, grand canonical versions of both the self-consistent field (GC-SCF) method and direct free-energy minimization with auxiliary Hamiltonians (GC-AuxH) method are able to rapidly converge the grand free energy of electrons. This substantially improves upon the current state of the art of running an outer loop over conventional fixed-charge DFT calculations. With detailed tests of the convergence of all these algorithms, we show that the GC-AuxH method is the most suitable default choice exhibiting smooth exponential convergence to the minimum.

Grand-canonical DFT directly mimics the experimental condition in electrochemical systems, where electrode potential sets the chemical potential of electrons, and the number of electrons at the electrode surface (including adsorbates in the electrochemical interface) changes continuously in response. Describing this change in charge at the surface plays an important role in accurately modeling several electrochemical phenomena.^12–15 Here, we showcase the new algorithms by analyzing the under-potential deposition (UPD) of copper on platinum in an electrolyte containing chloride ions. We resolve an old debate about the identity of a second under-potential peak, showing that partial copper monolayers are not plausible and that the second peak is due to desorption of chloride ions. We expect the new methods presented here to substantially advance the realistic treatment of electrochemical phenomena in first principles calculations.

R.S. and W.A.G. acknowledge support from the Joint Center for Artificial Photosynthesis (JCAP), a DOE Energy Innovation Hub, supported through the Office of Science of the U.S. Department of Energy under Award No. DE-SC0004993. R.S. and T.A.A. acknowledge support from the Energy Materials Center at Cornell (EMC²), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0001086. Calculations in this work used the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We thank Kendra Letchworth-Weaver, Kathleen Schwarz, Yuan Ping, Hai Xiao, Tao Cheng, Robert J. Nielsen, and Jason Goodpaster for insightful discussions.