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

1.2 MB·

by  Ravishankar Sundararaman

Grand canonical electronic density-functional theory:
algorithms and applications to electrochemistry
Ravishankar Sundararaman,1,2,a) William A. Goddard III,3,2,b) and Tomas A. Arias4,c)
1)Department of Materials Science and Engineering, Rensselaer Polytechnic Institute, Troy,
NY 12180
2)Joint Center for Artificial Photosynthesis, California Institute of Technology, Pasadena, CA 91125
3)Materials and Process Simulation Center, California Institute of Technology, Pasadena,
CA 91125
4)Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY 14853
(Dated: 1 March 2017)
7
1 First-principles calculations combining density-functional theory and continuum solvation models enable re-
0 alistictheoreticalmodelinganddesignofelectrochemicalsystems. Whenareactionproceedsinsuchsystems,
2 the number of electrons in the portion of the system treated quantum mechanically changes continuously,
withabalancingchargeappearinginthecontinuumelectrolyte. Agrand-canonicalensembleofelectronsata
b
e chemicalpotentialsetbytheelectrodepotentialisthereforetheidealdescriptionofsuchsystemsthatdirectly
F mimics the experimental condition. We present two distinct algorithms, a self-consistent field method (GC-
SCF) and a direct variational free energy minimization method using auxiliary Hamiltonians (GC-AuxH), to
8
2 solvetheKohn-Shamequationsofelectronicdensity-functionaltheorydirectlyinthegrandcanonicalensemble
at fixed potential. Both methods substantially improve performance compared to a sequence of conventional

] fixed-number calculations targeting the desired potential, with the GC-AuxH method additionally exhibiting
h
reliableandsmoothexponentialconvergenceofthegrandfreeenergy. Finally,weapplygrand-canonicalDFT
p
to the under-potential deposition of copper on platinum from chloride-containing electrolytes and show that
-
m chloridedesorption,notpartialcoppermonolayerformation,isresponsibleforthesecondvoltammetricpeak.
e
h
c Density-functional theory (DFT) enables theoretical reliable yet efficient treatment of electrochemical sys-
. elucidation of reaction mechanisms at complex catalyst tems.
s
c surfaces, making it now possible to design efficient het- Second, electrons can flow in and out of the electrode
si erogeneous catalysts for various industrial applications as electrochemical reactions proceed. Changes in elec-
y from first principles, for example for high-temperature tronic charge of electrode surfaces and adsorbates can
h gas-phasetransformationofhydrocarbonstoavarietyof beespeciallyimportantbecausetheelectrolytestabilizes
p valuable chemical products.1,2 The extension of this pre- charged configurations with a counter charge from the
[
dictive power to electrocatalysis would be highly valu- ionic response. For example, reduction of formic acid on

2 able for an even broader class of technological problems, platinum at experimentally relevant potentials is dom-
v including a cornerstone of future technology for renew- inated by formate ions rather than neutral molecules
0 able energy: converting solar energy to chemical fuels at the surface.12 Proton adsorption on stepped and
9
by electrochemical water splitting and carbon dioxide polycrystalline surfaces involves displacing oxidatively-
4
4 reduction.3 Accurately describing electrochemical phe- adsorbed water at relevant potentials, resulting in non-
0 nomena, however, presents two additional challenges. integer charge transfers and an anomalous pH depen-
. First, the electrolyte, typically consisting of ions in a dence deviating from the Nernst equation.13
1
liquid solvent, strongly affects the energetics of struc-
0 Accounting for the electrolyte response using our sol-
tures and reactions at the interface. Treating liquids di-
7 vation models,8,11 and adjusting the electron number to
1 rectly in DFT requires expensive molecular dynamics to match experimentally relevant electrode potentials, real-
: sample the thermodynamic phase space of atomic con-
v istic predictions of electrochemical reaction mechanisms
figurations. Historically, a number of continuum solva-
i have now become possible.14 In particular, application
X tion models that empirically capture liquid effects have
of this methodology to the reduction of CO on Cu(111)
r enabled theoretical design of liquid-phase catalysts.4,5
predicts onset potentials for methane and ethene forma-
a More recently, empirical solvation models suitable for
tion with 0.05 V accuracy in comparison to experiment,
solid-liquid interfaces,6–8 joint density-functional theory
for a wide range of pH varying from 1 to 12.15 How-
(JDFT) for efficiently treating liquids with atomic-scale
ever, conventional DFT software and algorithms are op-
structure,9 andminimally-empiricalsolvationmodelsde-
timized for solving the quantum-mechanical problem at
rived from JDFT,10,11 have made great strides towards
fixedelectronnumber,requiringextrawork(bothmanual
andcomputational)tocalculatepropertiesforaspecified
electrode potential.
a)Electronicmail: [email protected] Electricpotentialsandfieldsplayanimportantrolein
b)Electronicmail: [email protected] fields besides than electrochemistry. Density-functional
c)Electronicmail: [email protected] theory approaches accounting for electric potential have
2
been developed in special cases for field emission from labelshereonwardfornotationalsimplicity; allthefunc-
metal surfaces using a jellium model,16 and for calculat- tionals below depend on temperature).
ing capacitance in metal-insulator-metal17,18 and carbon Most practical approximations in electronic density-
nanotube systems.19 Calculating non-equilibrium trans- functionaltheoryfollowtheKohn-Shamapproach24 that
port of electrons in nanoscale systems also requires ac- includestheexactfreeenergyofanon-interactingsystem
counting for potential difference between reservoirs in a of electrons with the same density n(r). The universal
DFT calculation.20 First-principles molecular dynamics functional is typically split as
approaches have been developed to emulate fixed po-
tential using fluctuating numbers of electrons between A [n]=A [n]+(cid:90) dr(cid:90) dr(cid:48)n(r)n(r(cid:48))+E [n], (2)
time steps.21 However, in all these cases, each involved HKM ni 2|r−r(cid:48)| XC
self-consistent DFT calculation contains a fixed num- (cid:124) (cid:123)(cid:122) (cid:125)
ber of electrons and is carried out using a conventional EH[n]
canonical-ensemble algorithm. where A [n] is the non-interacting free energy (which
ni
This paper introduces algorithms for grand canonical we describe in detail below), the second ‘Hartree’ term
DFT, where electron number adjusts automatically to E [n] is the mean-field Coulomb interaction between
H
target a specified electron chemical potential (related to electrons (using atomic units e,m ,(cid:126),k = 1 through-
e B
electrode potential), thereby enabling efficient and in- out), and the final ‘exchange-correlation’ term E [n]
XC
tuitive first-principles treatment of electrochemical phe- captures the remainder which is not known and must
nomena. Section I summarizes the theoretical back- hence be approximated. The exact free energy of non-
ground of first-principles electrochemistry using JDFT interacting electrons is
and continuum solvation models, and sets up the funda-
mental basis of grand-canonical DFT. Then, section II A [n]= min (cid:88)(cid:18)fi (cid:90) dr|∇ψ (r)|2 −T S(f ) (cid:19),
itnintrcotdculacsessetshoefmDoFdTifiaclagtoiornitshmnesc,etshsaersyeltfo-cmonaskisetetnwtofideilsd- ni {ψ→i(rn)(,rf)i} i (cid:124)2 (cid:123)(cid:122) i (cid:125) E(cid:124)nt(cid:123)r(cid:122)oip(cid:125)y
(GC-SCF) method and the variational free energy min- Kinetic
(3)
imization using auxiliary Hamiltonians (GC-AuxH), di-
which includes the kinetic energy and entropy contri-
rectlyconvergethegrandfreeenergyofelectronsatfixed
butions of a set of orthonormal single-particle orbitals
potential. SectionsIIIBandIIICestablishthealgorithm
ψ (r) with occupation factors f ∈ [0,1]. Above, the
parameter(s) that optimize the iterative convergence of i i
single-particleentropyfunctionisS(f)=−flogf−(1−
the GC-SCF and GC-AuxH methods respectively, while
f)log(1−f). Notethatwelettheorbitalindexiinclude
section IIID compares the performance of these algo-
spin degrees of freedom as well, and therefore do not in-
rithms for a number of prototypical electrochemical sys-
troduce factors of 2 for spin degeneracy. For notational
tems. Finally, section IIIE demonstrates the utility of
convenience, we also let i include Bloch wave-vectors in
grandcanonicalDFTbysolvinganelectrochemicalmys-
the Brillouin zone (in addition to spin and band indices)
tery: theidentityofthesecondvoltammetricpeakinthe
for periodic systems.
under-potentialdeposition(UPD)ofcopperonplatinum
The minimization in (3) is constrained such that the
in chloride-containing electrolytes.
density of the non-interacting system
(cid:88)
f |ψ (r)|2 =n(r), (4)
I. THEORY i i
i
A. Background: electronic density functional theory the density of the real interacting system. Performing
this minimization with Lagrange multipliers V (r) (the
KS
The exact Helmholtz free energy A of a system of Kohn-Sham potential) to enforce the density constraint
interacting electrons in an external potential V(r) at and εi (Kohn-Sham eigenvalues) to enforce orbital nor-
a finite temperature T satisfies the Hohenberg-Kohn- malization constraints, results in a set of single-particle
Mermin variational theorem22,23 Schr¨odinger-like equations
(cid:18) (cid:90) (cid:19)
A=min AT [n(r)]+ drV(r)n(r) , (1) ∇2
n(r) HKM − 2 ψi(r)+VKS(r)ψi(r)=εiψi(r) (5)
where ATHKM is a universal functional that depends only for stationarity with respect to ψi(r), and the Fermi oc-
on the electron density n(r) (and temperature), and not cupation condition
ontheexternalpotential. However,constructingapprox-
imations for this unknown universal functional that ac- 1
f = (6)
curately capture the energies and geometries of chemical i 1+expεi−µ
T
bondsintermsofthedensityaloneisextremelychalleng-
ing, partly because the quantum mechanics of the elec- for stationarity with respect to f . Here, the electron
i
trons is completely implicit in A [n] (dropping the T chemical potential µ appears as a Lagrange multiplier to
HKM
3
enforce the electron number constraint, and is chosen so Next to an electrolyte, the charge of a partially-
(cid:80)
that f =N, the number of electrons in the system. periodic system (one-dimensional ‘wires’ or two-
i
Optimizing the total free energy functional (1) with dimensional‘slabs’)isnolongerconstrained,allowingthe
A [n]implementedby(2,3)thenyieldsthestationar- number of electrons to vary. Now, the absolute reference
HKM
ity condition with respect to electron density, for the electron chemical potential µ does become physi-
callymeaningful,andµnowcorrespondstotheelectrode
δ potential that controls the number of electrons in the
V [n](r)=V(r)+ (E [n]+E [n]). (7)
KS δn(r) H XC electrode.
However, treating electrochemical systems using elec-
Conventionaldensity-functionaltheorycalculationsthen tronic density-functional theory alone is extremely chal-
amounttoself-consistentlysolvingtheKohn-Shamequa- lengingforavarietyofreasons. First,treatingliquidsre-
tions (5) along with (7), coupled via the electron density quiresastatisticalaverageoveralargenumberofatomic
constraint (4). configurations to integrate over thermodynamic phase
TheKohn-Shampotentialisarbitraryuptoanoverall space. For this, techniques such as molecular dynam-
additive constant: changing this constant introduces a ics typically require calculations of at least 104 − 105
rigid shift in the eigenvalues (cid:15) and the electron chemi- atomic configurations (instead of just one for a solid).
i
cal potential µ, but does not affect the occupations f , Second,suchcalculationsrequirealargenumberofliquid
i
electrondensityorfreeenergy. Forfinitesystems(ofany molecules to minimize finite size errors in the molecular
charge) and for neutral systems that are infinitely peri- dynamics. For example, for electrolytes with a realistic
odic in one or two directions, this arbitrariness can be ionic concentration of 0.1 M, there is on average one ion
eliminated by requiring that the potential vanishes in- for a few hundred solvent molecules. Making statistical
finitely far away from the system, giving meaning to the errorsintheionnumbermanageableinsuchcalculations
absolute values of µ and (cid:15) as being referenced to ‘zero therefore requires > 103 solvent molecules with > 104
i
atinfinity’. However,formaterialsthatareinfiniteinall electrons, contrasted with a typical 10−100 atoms with
threedimensions,suchasperiodiccrystallinesolids,there 100−103electronsintheelectrodeslab+adsorbateofin-
is no analogous natural choice for the zero of potential, terest. Combined, these factors make density-functional
making the absolute reference for µ and (cid:15) meaningless. molecular dynamics simulations of electrochemical sys-
i
Moreimportantly,forsystemsthatareperiodicinone, tems prohibitively expensive computationally, in addi-
two or all three directions, the net charge per unit cell tional to being difficult to set up and analyze.
mustbezero, otherwisetheenergyperunitcellbecomes A viable alternative to the above direct approach is to
infinite. In terms of a finite system size L and then tak- employjointdensity-functionaltheory9 (JDFT),avaria-
ingthelimitL→∞,theenergyperunitcellofasystem tionaltheoremakintotheHohenberg-Kohntheoremthat
withnetchargeperunitcelldiverges∝lnLforoneperi- makesitpossibletodescribethefreeenergyofasolvated
odic direction, ∝L for two periodic directions and ∝L2 systemintermsoftheelectrondensityn(r)forthesolute
for three periodic directions. Therefore, the number of andintermsofasetofnucleardensities{N (r)}(where
α
electrons per unit cell is physically constrained to keep α indexes nuclear species) of the solvent (or electrolyte).
theunitcellneutralinsystemswithany periodicity,and Specifically, the equilibrium free energy of the combined
onlyfinitesystems(likemoleculesandions)havenumber solute and solvent systems minimizes
of electrons as a degree of freedom.
(cid:18)
A= min A˜ [n(r),{N (r)}]
JDFT α
n(r),{Nα(r)}
B. Electrochemistry with joint density-functional theory (cid:90) (cid:90) (cid:19)
(cid:88)
+ drV(r)n(r)+ drV (r)N (r) , (8)
α α
Fordescribingelectrochemicalsystemsandelectrocat- α
alytic mechanisms, we are typically interested in ad-
where V(r) is the external electron potential, V (r) is
sorbed species in a solid-electrolyte interface which ex- α
the external potential on the liquid nuclei and A is
change electrons with the solid (electrode). In these sys- JDFT
a universal functional independent of these external po-
tems, the solid surface, which we would describe in a
tentials. Separating out the Hohenberg-Kohn electronic
density-functionalcalculationasaslabperiodicintwodi-
density functional A [n] for the solute, the total free
rections,doeshaveanetchargeperunitcellthatdepends HK
energy is
on the electrode potential. In contrast to the discussion
at the end of the previous section, this is physically pos-
sible (i.e. has a finite energy) because the electrolyte A˜JDFT[n,{Nα}]=AHKM[n]+A˜diel[n,{Nα}]. (9)
(cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125)
contains mobile ions that respond by locally increasing
electronic solvation
theconcentrationofionsofoppositechargenearthesur-
face, thereby neutralizing the unit cell. (The charge per In practice, the functional A , much like A , is un-
diel HKM
unit areas of the electrode and electrolyte are equal and known and needs to be approximated. Importantly, the
opposite.) liquid is now described directly in terms of its average
4
densityratherthanindividualconfigurations,andtheex- ρ (r) = n(r)+ρ (r), where ρ (r) is the solute nu-
el nuc nuc
pensivequantum-mechanicaltheoryoftheelectronsisre- clear charge density. The screened Coulomb interaction
stricted to the solute alone, thereby addressing both the term above is expressed in terms of the bare Coulomb
sampling and system-size problems that make density- interaction Kˆ and the solvent susceptibility operator χˆ,
functional molecular dynamics prohibitively expensive. defined for the local-response models by
Typically, we are interested in a situation where the (cid:18)(cid:15) −1 (cid:19) κ2
external potentials on the liquid are zero, and the liq- χˆ·φ(r)≡∇· b s(r)∇φ(r) − s(r)φ(r). (11)
4π 4π
uid only interacts with itself, and with the electrons and
(cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125)
nuclei of the solute. In this situation, we can perform Dielectric Ionic
the optimization over liquid densities and define the im-
Here, (cid:15) is the bulk dielectric constant of the solvent
plicit functional AJDFT[n] = AHKM[n]+Adiel[n] where and κ =b (cid:112)4π(cid:80) N Z2/T is the inverse Debye screen-
we define A [n] ≡ min A˜ [n(r),{N (r)}] for X i i i
X {Nα(r)} X α ing length in vacuum of a set of ionic species of charge
=JDFTandX =diel. Wewillworkwiththesereduced
Z withconcentrationsN (thenetinverseDebyescreen-
functionals below for simplicity. i i
ing length in the presence of the background dielectric is
The framework of joint density-functional theory en- √
κ/ (cid:15) ). The cavity shape function s(r) modulates the
compasses an entire hierarchy of solvation theories. Fur- b
dielectric and ion response, and varies from zero in the
ther separating the solvation term A into a classical
diel solute region (no solvent present) to unity in the solvent
density functional for the liquid25–27 and an electron-
at full bulk density. Finally, the second term of (10),
liquid interaction functional, unlocks the full potential
A ,empiricallycaptureseffectsbeyondmean-fieldelec-
of JDFT to describe atomic-scale structure in the liquid cav
trostatics, such as the free energy of cavity formation in
without statistical sampling. Starting from ‘full-JDFT’
the liquid and dispersion interactions between the solute
and performing perturbation theory for linear-response
and solvent.29 Different solvation models at this level of
of the liquid results in the non-empirical SaLSA solva-
the hierarchy only differ in the details of how s(r) is de-
tionmodel10 whichcontinuestocapturetheatomic-scale
termined from the electron density (or even from atomic
nonlocality in the liquid response and introduces no fit
positions and fit radii for the PCM5 and SMx4 solvation
parameters for the electric response of the solvent.
models), and in the details of A .
cav
At the simplest end of the JDFT hierarchy, are con-
In practice, evaluating the first term of (10) requires
tinuum solvation models that neglect the nonlocality the calculation of φ(r) = (Kˆ−1 − χˆ)−1ρ (r), which is
of the liquid response and replace it by that of an el
the net electrostatic potential including screening by the
empirically-determined dielectric cavity (optionally with
solvent (and electrolyte). The second part of the first
Debye screening due to electrolytes). This includes our
term in (10) represents the electrostatic interactions be-
recent CANDLE solvation model11 that builds on the
tween electrons and the nuclei in vacuum which can-
stability of SaLSA for highly polar systems, and earlier
cels corresponding terms in the vacuum DFT functional
solvationmodelssuitableformoleculesandless-polarsys-
(specifically the Hartree, Ewald and long-range part of
tems such as our GLSSA13 model7 (or its equivalent,
the local pseudopotential terms), so that all long-range
VASPsol8) and the comparable Self-Consistent Contin-
terms that contribute to the net free energy are present
uum Solvation (SCCS) model.6,28 Even the traditional
in the first term of (10). Finally, substituting (11) and
quantum-chemistry finite-system solvation models such
Kˆ−1 = −∇2/(4π) into the definition of φ(r) results in
as the PCM series5 and the SMx series4 can be mapped
the linearized Poisson-Boltzmann (modified Helmholtz)
on to this class of solvation models.
equation
For simplicity, we will work here with this simplest
classofcontinuumsolvationmodels. Thetheoreticalcon- −∇·((cid:15) s(r)∇φ(r))+κ2s(r)φ(r)=4πρ (r). (12)
b el
siderations and algorithms in this work focus primarily
the electronic density-functional theory component, are Without ionic screening (second term of (12) above),
largely agnostic to the internals of the solvation model the absolute reference for φ(r) is undetermined and does
and therefore straightforwardly generalize up the JDFT not contribute to the bound charge charge induced in
hierarchytothemoredetailedandcomplexsolvationthe- the liquid, χˆφ(r) (given by the first term of (11) alone),
ories. Essentially, all the simple electron-density-based because∇(constant)=0. However,withionicscreening,
continuumsolvationmodels6–8,11 canbesummarizedab- a constant shift of φ(r) does affect the second terms of
stractly as7,10 (11) and (12), producing a charge response in the liquid
and making the absolute reference of the electrostatic
(cid:90) (Kˆ−1−χˆ)−1−Kˆ potentialφ(r)meaningful. Inparticular,integrating(12)
Adiel[n(r)]= drρel(r) 2 ρel(r)+Acav[s]. over space yields
(10) (cid:90) κ2 (cid:90)
Here, the first term is the electrostatic solute-solvent dr4πs(r)φ(r)= ρel(r). (13)
interaction energy given by the difference between
the solvent-screened Coulomb interaction and the bare From (11), the left hand side is the negative of the total
CoulombinteractionKˆ ofthetotalsolutechargedensity, bound charge in the liquid, while the right hand side
5
is the total charge of the solute system. Therefore, the atafixednumberofelectronsN inthesolutesubsystem
continuumelectrolyteautomaticallycompensatesforany (electrode + adsorbates). However, in electrochemical
net charge in the solute and makes the complete system systems, N is an artificial constraint because electrons
neutral. As a side effect, the absolute reference of the can freely exchange between the electrode and an exter-
electrostatic potential is meaningful and automatically nalcircuit. Instead,experimentssettheelectrodepoten-
corresponds to ‘zero at infinity’. This happens because tial U, and hence the electron chemical potential µ, and
√
the Greens function of (12) is exp(−κr/ (cid:15) )/r in the N adjusts accordingly as a dependent variable. Thermo-
b
bulk liquid with finite κ (instead of 1/r), causing φ(r) dynamically, this corresponds to switching the electrons
to exponentially decay to zero outside the solute region from the finite-temperature, fixed-number canonical en-
(where ρ =0). semble to the finite-temperature, fixed-potential grand-
el
Ref.30givesadetailedversionoftheabovediscussion canonical ensemble. Correspondingly, the relevant free
including rigorous proofs of the absoluteness of the ref- energyminimizedatequilibriumisthegrandfreeenergy
erence for φ(r), and numerical details and algorithms for Φ=A−µN, instead of the Helmholtz free energy A.
efficientlysolving(12)withperiodicboundaryconditions The most straightforward approach to fixed-potential
in plane-wave basis DFT calculations. The result that DFTforelectrochemistryistorepeatconventionalfixed-
the electrolyte neutralizes charges in the solute is true charge DFT calculations at various electron numbers N
moregenerally,evenforsolvationmodelswithnonlinear7 to reach a target chemical potential µ. Using a steep-
and/ornonlocalresponse.10 However,asRef.7describes est descent / secant method to optimize N results in
indetail,forthegeneralnonlinearcase,unitcellneutral- convergence of µ typically within 10 iterations.14 How-
ization is not automatic and must be imposed using a ever this approach is inefficient since it requires multiple
Lagrange multiplier constraint; this Lagrange multiplier DFT calculations to calculate the grand free energy and
then fixes the absolute value of φ(r) such that it is zero charge of a single configuration, typically taking 3x the
at infinity. timeofasinglefixed-chargecalculation(not10xbecause
The key conclusion of this discussion is that contin- subsequent calculations have a better starting point and
uum electrolytes present two advantages. First, auto- converge quicker).14
matic chargeneutralization implies wecan performwell- A more efficient approach would be to directly opti-
defined calculations with net charge per unit cell in the mize the Kohn-Sham functional in the grand-canonical
solute. Second, meaningful absolute reference values for ensemble. For a solvated system, this corresponds to
potentials implies that the electron chemical potential a small modification of the minimization problem (8)
µ (which determines the electron occupations in (6)) is with A given by (9) and with A evaluated us-
JDFT HKM
also referenced to zero at infinity. Consequently, µ is re- ing the Kohn-Sham approach (2,3). The only differences
(cid:80)
lated directly to the potential U of the electrode provid- are that the Lagrange multiplier term −µ( f − N)
i i
ing the electron reservoir in experiments. This potential that enforced the electron number constraint is replaced
(cid:80)
istypicallyreferencedtothestandardhydrogenelectrode by −µ f which implements the Legendre transform
i i
(SHE), so that of the Helmholtz free energy to the grand free energy,
and that the fixed-N constraint is removed. The elec-
µ=µ −U, (14)
SHE tron occupation factors are still Fermi functions (6), but
µ is specified as an input (instead of being adjusted to
where µ is the absolute position of the standard hy-
SHE match a specified N). Operationally, this amounts again
drogen electrode relative to the vacuum level (i.e. the
to solving the Kohn-Sham eigenvalue problem (5) self-
zeroatinfinityreference). Calculationscanemployeither
consistently with
theexperimentalestimateµ =−4.44eVtorelatethe
SHE
absoluteandrelativepotentialscales,31,32oratheoretical δ
V [n](r)=V(r)+ (E [n]+E [n]+A [n]),
calibrationbasedon thecalculated andmeasured poten- KS δn(r) H XC diel
tials of zero charge of solvated metal surfaces.30 The lat- (15)
terapproachhastheadvantageofminimizingsystematic where there is now a single extra contribution to the
errors in the solvation model since they cancel between potential due to the solvent (electrolyte), and with the
the calculations used for calibration and prediction. For aforementioned changes to the Lagrange multipliers and
example, with the CANDLE solvation model we use be- constraints.
low, the calibrated µSHE =−4.66 eV.11 This conceptually simple modification, however,
presents numerical challenges to the algorithms com-
monly used for solving the Kohn-Sham problem, such as
C. Grand-canonical density-functional theory the self-consistent field (SCF) method, which solves the
Kohn-Shameigenvalueproblemfromaninputdensity(or
Usingjointdensity-functionaltheoryofcontinuumsol- Kohn-Sham potential), adjusting this density (or poten-
vation models to treat electrolytes as described above, tial)iterativelyuntilself-consistencyisachieved. Acom-
we can calculate the Helmholtz free energy and electro- mon instability for metallic systems in the SCF method
chemical potential (µ, and hence U) of specific micro- is ‘charge sloshing’: the electron density oscillates spa-
scopic configurations of adsorbates on electrode surfaces tiallybetweeniterationsinsteadofconverging. Switching
6
tothegrand-canonicalensemblecansignificantlyexacer- the form V˜ (G) = n˜(G)4π/G2, causing small G (long-
H
bate this problem: electrons can now additionally slosh wavelength)variationsoftheelectrondensitytoproduce
between the system and the electron reservoir. Here, we larger changes in the Kohn-Sham potential than large G
present modifications to the SCF method and an alter- ones. Uncompensated, this makes the electron density
nate algorithm that allow reliable and efficient conver- optimizationunstableagainstlong-wavelengthperturba-
gence for grand-canonical Kohn-Sham DFT. tions (the charge-sloshing instability). The Kerker pre-
conditioner
G2
II. ALGORITHMS K˜(G)=A (19)
G2+q2
K
A. Self-Consistent Field method: Pulay mixing mitigates this problem by suppressing the contribution
of the problematic small G components in determining
Given electron density n(i)(r), the Kohn-Sham equa- the next input density. The metric
in
tions (5) with potential V given by (15) define orbitals
KS G2+q2
{ψj} and eigenvalues {εj}. At a given electron chemi- M˜(G)= M (20)
calpotentialµ, theseinturndefinetheoccupations{f } G2
j
given by (6) and a new electron density n(i)(r) given by enhances the contribution of small G components in the
out
(4). The Self-Consistent Field (SCF) method attempts residual norm, thereby prioritizing the optimization of
to find n(i)(r) such that n(oiu)t(r)=n(ini)(r). tthoresseqcoamnpdoqnen,twshwihchencodnettreorlmtihneinbgalnaoinnpcte.bTehtweeweanvsehvoerct-
There are two algorithmic ingredients to this method: K M
and long-wavelength components, and the prefactor A,
solution of the Kohn-Sham eigenvalue equations and
whichsetsthemaximumfractionofn thatcontributes
optimization of the electron density. The eigenvalue out
to the next n , can be adjusted to optimize the conver-
equations remain unchanged between conventional and in
gence of the SCF method. See Ref. 36 for a detailed dis-
fixed-potential calculations. To solve, these we use the
standard Davidson algorithm.33 The difficulty in fixed- cussionofthePulay-KerkerSCFapproachanditsperfor-
manceforconventionalfixedelectronnumber(canonical)
potential calculations arise in the electron density opti-
DFT calculations.
mization, which we discuss below.
The G = 0 component of the electron density, which
A robust and frequently-used algorithm for charge-
density optimization is Pulay mixing34 with Kerker equals N/Ω, where Ω is the unit cell volume, remains
preconditioning.35Briefly,Pulaymixingassumesthatthe fixed in canonical DFT calculations and is therefore ex-
cluded from the metric and the preconditioner. This is
residualR[n (r)]≡n (r)−n (r)isapproximatelylin-
in out in
no longer true in fixed-potential DFT, where the num-
earintheinputelectrondensityn (r),andcalculatesthe
in
ber of electrons changes. With the above prescriptions,
optimum input electron density as a linear combination
the Kerker preconditioner (19) →0 as G→0 which will
of previous iterations,
preventtheelectronnumberfromchangingbetweenSCF
nopt(r)=(cid:88)α n(i)(r) (16) iterations. Likewise, the G → 0 divergence of the Pulay
in i in
metric(20)causetheresidualnormtobecomeundefined
i
when the electron number changes between iterations.
Minimization of the norm of the corresponding residual In order to generalize the Pulay-Kerker SCF approach
forfixed-potentialDFT,wethereforeneedtofixtheG→
(cid:90)
F({α })=(cid:88)α α drR[n(i)(r)]MˆR[n(j)(r)] (17) 0 behavior of both the preconditioner and the metric.
i i j in in
The need for the preconditioner and the metric arose
ij
fromthereciprocalspaceCoulomboperator4π/G2inthe
(cid:80)
with the constraint iαi =1 yieldsa set of linear equa- Hartree potential. In a uniform electrolyte with Debye
tionsdeterminingthecoefficientsα ,whereMˆ isthemet- screening,thereciprocalspaceCoulomboperatorinstead
i
ricfordefiningthenormoftheresidual. Finally,thenext takes the form 4π/((cid:15) G2 +κ2) (from (12) in reciprocal
b
input density is obtained by mixing the optimum input space with s = 1). Since the electrolyte is responsible
density with its corresponding output density for fixing the indeterminacy of the G = 0 component of
the potential, a reasonable ansatz for extending Pulay-
n(i+1)(r)=nopt(r)+KˆR[nopt(r)], (18) KerkertofixedpotentialcalculationsisreplacingG2with
in in in
G2+q2, where
κ
where Kˆ is the Kerker preconditioning operator.
κ
The metric Mˆ and preconditioner Kˆ serve to bal- q = √ . (21)
κ (cid:15)
ance the influence of different components of the elec- b
tron density to the optimization procedure. These are In Section IIIB below, we show that setting
usually defined in reciprocal space, expanding n(r) =
(cid:80) n˜(G)eiG·r, where G are reciprocal lattice vectors. G2+q2
G K˜(G)=A κ (22)
In reciprocal space, the Hartree potential in (15) takes G2+q2 +q2
κ K
7
and iary subspace Hamiltonian matrix H as an indepen-
aux
dent variable,40 and setting the occupations f = f(η )
G2+q2 +q2 i i
M˜(G)= κ M (23) in terms of the eigenvalues {ηi} of Haux instead of the
G2+qκ2 Kohn-Sham eigenvalues {εi}. (Here, f(η) is the Fermi
function given by (6), and the electron chemical po-
indeed makes the grand canonical self-consistent field
tential µ is chosen so as to satisfy the electron num-
(GC-SCF) method function efficiently, with optimum (cid:80)
ber constraint f(η ) = N.) This eliminates the
convergence for q given by (21). i i
κ problematic inequality constraints on the occupations,
and upon minimization, the auxiliary subspace Hamil-
tonian Hij approaches the true subspace Hamiltonian,
B. Variational minimization: auxiliary Hamiltonian aux
Hij =(cid:104)ψ |Hˆ |ψ (cid:105),whereHˆ =−∇2/2+V (r)isthe
method sub i KS j KS i KS
Kohn-ShamHamiltonian. Bychoosingtheundetermined
unitaryrotationsoftheorbitals{ψ }todiagonalizeH ,
AnalternateapproachtosolvingtheKohn-Shamequa- i aux
Ref.40furthershowsthatthegradientofthefreeenergy
tions is to directly minimize the total (free-)energy func-
with respect to the independent variables simplifies to
tional (1), with A given by (2,3), in terms of the
HKM
Kohn-Sham orbitals as independent variables. For joint (cid:18) (cid:19)
density-functional theory, this amounts to δψδA(r) =f(ηi) HˆKSψi(r)−(cid:88)ψj(r)Hsjuib (25)
i
j
A= min (cid:34)(cid:88)(cid:18)fi (cid:90) dr|∇ψ (r)|2−TS(f )(cid:19) and
{ψi(r),fi} i 2 i i
∂A ∂f(η )
(cid:90) (cid:21) =δ (Hii −η ) i
+EH[n]+EXC[n]+Adiel[n]+ drV(r)n(r) , (24) ∂Haijux ij sub i ∂ηi
−δ ∂µ (cid:88)(Hkk −η )∂f(ηk)
wheren(r)isnowderivedfrom{ψ }and{f }asgivenby ij∂η sub k ∂η
i i i k
k
(4). In the above minimization, the orbitals {ψ } must
i f(η )−f(η )
be orthonormal, the occupation factors must satisfy 0≤ +(1−δ )Hij i j . (26)
(cid:80) ij sub η −η
f ≤1, and optionally, f =N for the canonical fixed i j
i i i
electron number case.
These gradients are used to perform line minimization
For insulators at T =0, the occupations f are known
i alongasearchdirectioninthespaceofindependentvari-
inadvance. Theconstrainedoptimizationoverorthonor-
ables. With every update of the independent variables,
mal {ψ } is most efficiently carried out using a precon-
i the orbitals are re-orthonormalized and the unitary ro-
ditioned conjugate-gradients (CG) algorithm on uncon-
tations of the orbitals are updated to keep the auxiliary
strained orbitals using the analytically-continued energy
Hamiltonian diagonal.
functional approach.37 Briefly, this approach evaluates
In the preconditioned CG algorithm, the next search
the energy functional (24) on a set of orthonormal or-
direction is obtained as a linear combination of the cur-
bitals, which are a functional of the unconstrained or-
rent search direction and the preconditioned gradients
bitals used for minimization. (See Ref. 37 for further
given by
details.)
Thegeneralcaseofmetallicsystemsand/orfiniteT ad- δA
K =Tˆ (27)
ditionallyrequirestheoptimizationof{fi},whichischal- ψi(r) invδψ (r)
i
lengingfornon-linearoptimizationalgorithmsbecauseof
the inequality constraints 0 ≤ f ≤ 1. One possibility is and
i
to update the fillings from the Kohn-Sham eigenvalues
using(6)aftereveryfewstepsoftheCGalgorithm,38but KHaijux =−K(Hsiujb−ηiδij), (28)
thishinderstheconvergenceofCGbecausethefunctional
effectivelychangeseachtimethefillingsarealtered. The where Tˆinv and K are preconditioners. The role of pre-
‘EnsembleDFT’approach39 rectifiesthisconvergenceis- conditioning is to balance the weight of different direc-
suebyoptimizingtheoccupationfactorsatfixedorbitals tions in the minimization space in the explored search
in an inner loop, and performing the optimization of or- directions, and ideally the preconditioner equals the in-
bitals in an outer loop using the CG method, but this verse of the Hessian (which is difficult to compute ex-
increases the computational cost compared to the case actly). For the orbital directions, the standard precondi-
of the insulators. The SCF approach is typically much tioner Tˆ resembles the inverse of the dominant kinetic
inv
more computational efficient than these variants of the energy operator in the Kohn-Sham Hamiltonian.37 For
direct variational minimization algorithm with variable the auxiliary Hamiltonian direction, the preconditioner
occupations.36 removes the Fermi function derivatives and finite differ-
An alternate strategy for direct variational minimiza- ence factors from (26) in order to more equitably weight
tion with variable occupations is to introduce an auxil- all components of H . The preconditioning factor K
aux
8
controlstheoverallcontributionoftheH components III. RESULTS
aux
relativetothatoftheorbitalcomponents,andisadjusted
to achieve optimum convergence. A. Computational details
In this auxiliary Hamiltonian (AuxH) approach, the
We implement all algorithms and perform all cal-
orbitalsandoccupationsarecontinuouslyandsimultane-
culations using the open-source plane-wave density-
ouslyoptimizedtominimizethetotalfreeenergy, result-
functional theory software, JDFTx.41 Below, we spec-
ing in better convergence and computational efficiency
ifycomputationalandconvergenceparametersinatomic
comparable to the fixed-occupations insulator case, and
units (distances in bohrs a ≈ 0.529 ˚A and en-
competitive with the SCF method even for metallic sys- 0
ergies in Hartrees E ≈ 27.2 eV), but present any
tems. See Ref. 40 for further details on the algorithms h
physically relevant properties in conventional units (˚A,
andperformancecomparisonsforconventionalfixedelec-
eV). All calculations in this work employ the PBE42
tron number calculations.
exchange-correlation functional with GBRV ultrasoft
pseudopotentials43 atakineticenergycutoffof20E for
Now, for fixed potential calculations, we set the occu- h
Kohn-Sham orbitals and 100 E for the charge density.
pations to Fermi functions of the auxiliary Hamiltonian h
The metal surface calculations use inversion-symmetric
eigenvalues at a specified µ, instead of selecting µ based
slabs of at least five layers, with at least 15 ˚A vacuum
on the electron number constraint. Correspondingly, the
separation and truncated Coulomb potentials44 to min-
second term of the auxiliary Hamiltonian gradient (26),
imize interactions with periodic images. For Brillouin
which arises from this constraint, drops out, and the al-
zone integration, we use a Fermi smearing of 0.01 E
gorithm requires no further modification. h
and a Monkhorst-Pack k-point mesh along the periodic
directions with the number of k-points chosen such that
The convergence rate of this algorithm, however, is
the effective supercell is larger than 30 ˚A in each direc-
sensitive to the preconditioning factor K and we pro-
tion. We use the CANDLE solvation model to describe
poseamodifiedheuristictoupdateK automaticallyand
the effect of liquid water and Debye screening due to 1M
continuously. At the end of each line minimization, the
electrolyte, whichweshowedrecentlytomostaccurately
derivativeoftheoptimizedfreeenergyA withrespect
min
capturethesolvationofhighly-chargednegativeandpos-
toK canbeevaluatedfromtheoverlapbetweentheaux-
iliary Hamiltonian gradient and search direction.40 The itive solutes.11 We emphasize that the methods and al-
gorithms described above do not rely on specific choices
line-minimized energy is optimum for ∂A /∂K = 0.
min
for the pseudopotential, exchange-correlation functional,
Therefore, if variation of A with respect to K is con-
min
k-mesh or solvation model; we keep these computational
vex, we should increase K if we find ∂A /∂K <0 and
min
parameters constant here for consistency.
viceversa. HoweverconvexityisoftenlostifK isinitial-
ized at too high a value. Therefore, our heuristic tries
to zero ∂A /∂K while limiting the contribution of the
min
B. Convergence of the GC-SCF method
auxiliaryHamiltoniangradient. Inparticular,weupdate
The self-consistent field (GC-SCF) algorithm summa-
rized in section IIA depends on several parameters that
(cid:20) (cid:18) (cid:18) (cid:19)(cid:19) (cid:21)
K ←K×max exp f −∂Amin/∂K , gtot control its iterative convergence. All these parameters
sat g 2g are common to the conventional fixed-charge version
tot aux
(29) (SCF) and the fixed-potential variant (GC-SCF) intro-
at√the end of each line minimization, where fsat(x) ≡ duced here, except for the low-frequency cutoff wavevec-
x/ 1+x2 tosaturatethefactorbywhichK canchange tor q that is necessary to allow the net electron num-
κ
in one iteration, g is the overlap of the H compo- ber to change in the fixed-potential case. Figure 1
aux aux
nents of the gradient and preconditioned gradient, and compares the dependence of GC-SCF convergence on
g is the total overlap of the gradient and precondi- q for a prototypical calculation of an electrochemical
tot κ
tioned gradients (orbital + H ). Finally, we reset the system: a Cu(111) surface treated using a five-layer
aux
conjugate-gradient algorithm (i.e. set the search direc- inversion-symmetricslabsurroundedby1Maqueousnon-
tion to the negative of the preconditioned gradient di- adsorbing electrolyte treated using the CANDLE sol-
rection) after K has increased or decreased by a factor vation model,11 at a fixed potential of µ = −0.208E
h
greater than e2. We do this because dynamically chang- (1V SHE). (The remaining GC-SCF parameters are set
ing the preconditioner technically invalidates the strict to their default values which we discuss below.) The
orthogonality of the CG search direction with previous best convergence is obtained with q = 0.17 which cor-
κ
directions. Section IIIC below shows that this heuristic respondstotheDebyescreeninglengthoftheelectrolyte
exceeds the convergence obtained with fixed K, while (21). For small q , including the conventional case of
κ
section IIID shows that the grand-canonical auxiliary q = 0, the number of electrons does not respond suffi-
κ
Hamiltonian (GC-AuxH) algorithm consistently outper- ciently quickly, stalling at about 0.1 electrons from the
forms GC-SCF for fixed-potential calculations. converged value, correspondingly with the free energy
9
103 qκ = 0 103 qK = 0.0
1 qκ = 0.01 1 qK = 0.2
E| []h0 1100-1 qqκκq κ== = 00 ..001.474 E| []h0 1100-1 qqqKKK === 001...486
ΦΦ - 10-3 qκ = 0.7 ΦΦ - 10-3
| 10-5 | 10-5
-7 -7
10 10
1 1
10 10
|0 |0
N -1 N -1
N - 10 N - 10
| -3 | -3
10 10
0 10 20 30 40 0 10 20 30 40
GC-SCF cycles GC-SCF cycles
FIG. 1. Dependence of GC-SCF convergence on low- FIG. 2. Dependence of GC-SCF convergence on Kerker-
frequency cutoff wavevector q . The upper panel shows the mixing wavevector q . Good convergence is observed near
κ K
convergenceofthegrandfreeenergyΦ(towardsitsfinalvalue the typical recommended value q = 0.8 a−1 (≈ 1.5 ˚A-1).
K 0
Φ )onalogarithmicscale,andthelowerpanelshowsthatof Convergenceisrelativelyinsensitivetoq nearthisvalue,but
0 K
theelectronnumberN (towardsitsfinalvalueN ). Bestcon- becomesunstableforsmallq approachingthelow-frequency
0 K
vergence is obtained with q = 0.17a−1, the inverse Debye- cutoff q . System and remaining details are identical to Fig-
κ 0 κ
screeninglength,whichweuseasthedefaultvaluehenceforth. ure 1.
Results shown here are for a five-layer (5ML) Cu(111) slab
solvated in 1M CANDLE aqueous electrolyte, with potential
fixed to 1V SHE (µ=−0.208Eh), starting from a converged 103 qM = 0.0
neutral calculation of the same slab in vacuum. q = 0.2
1 M
10 q = 0.4
E| []h0 10-1 qqMMM == 01..86
stalling at about 0.1 E (≈ 2.7 eV) away from the con- Φ -3
vcheragnegdetvoaolurea.piLdalyr,gherinqdhκercinagusceosnvtheregeelneccetraonndneuvmenbteuraltlyo Φ| - 1100-5
leading to a divergence as seen for the case of qκ = 0.7. 10-7
An issue remains in the convergence independent of q :
κ 1
10
after initial convergence, the free energy oscillates at the
|0
10−6 E level,whiletheelectronnumberoscillatesatthe N -1
10−3 lehvel. N - 10
| -3
10
Keeping q at this optimum value given by (21), we
κ
next examine the dependence of GC-SCF convergence 0 10 20 30 40
on the remaining algorithm parameters for the same ex-
GC-SCF cycles
ample system. Figure 2 shows the dependence on the
Kerker mixing wavevector qK, which helps stabilize the FIG. 3. Dependence of GC-SCF convergence on Pulay-
GC-SCFalgorithmagainstlong-wavelengthchargeoscil- metric wavevector q . Convergence is relatively insensitive
M
lations. Optimal convergence is obtained for the typical to q , and we set q = q = 0.8 a−1 henceforth. System
M M K 0
recommended value36 of 0.8 a−1 (≈ 1.5 ˚A-1). As ex- and remaining details are identical to Figure 1.
0
pected, convergence is relatively insensitive to the exact
choiceofq ,aslongasq doesnotbecomesosmallthat
K K
convergenceisruinedbychargesloshing. Noticethatthe Finally,Figure4comparesthedependenceofGC-SCF
final convergence beyond the 10−6 Eh and 10−3 electron convergence on the maximum Kerker mixing fraction A,
level remains an issue that is not resolved for any choice which effectively controls what fraction of the new elec-
of qK. tron density is mixed into the current value. We find
Next, Figure3showsthevariationofGC-SCFconver- nominally best convergence for A = 0.5, but the perfor-
gence with the wavevector q controlling the reciprocal- mance of other values is not much worse. Smaller val-
M
space metric used by the Pulay algorithm. The conver- uesofAleadtogreaterstabilityinitially, butmarginally
gence is entirely insensitive to this choice, and we hence- slowerconvergencelateron,whilelargervaluesofAlead
forthsetq =q =0.8a−1 (therecommendedvalue36). to greater oscillations initially, but faster convergence
M K 0
Again, the final convergence issue remains unaffected by later on. Regardless, as before, good convergence is ob-
the choice of q . taineduntilthefreeenergyreachesthe10−6E level,but
M h
10
103 A = 0.1 K = 0.03
A = 0.2 10-3 K = 0.1
1
ΦE| []h0 1100--31 AAA === 001...580 ΦE| []h0 10-5 AutoK-a KK=d j==0u .s133t
Φ - 10 Φ -
| 10-5 |
-7
10
-7
10
1
10 -1
|0 |0 10
N -1 N
N - 10 N - -3
| -3 | 10
10
0 10 20 30 40 0 10 20 30 40 50
GC-SCF cycles Conjugate-Gradient iterations
FIG. 4. Dependence of GC-SCF convergence on maximum FIG. 6. Dependence of the convergence of the GC-
Kerker-mixing fraction A. Convergence is relatively insensi- AuxHvariational-minimizemethodonK,thepreconditioning
tivetoA,exceptthatitslowsdownforsmallAWehenceforth scalefactorforsubspacerotationsgeneratedbytheauxiliary
set A = 0.5 which nominally exhibits the best convergence. Hamiltonian. Near-optimal convergence is obtained when K
System and remaining details are identical to Figure 1. is automatically adjusted using the heuristic given by (29),
which we use by default henceforth. Calculations are for sol-
vated5MLCu(111)at1VSHE,startingfromthecorrespond-
103 5ML + 15a0 ing neutral vacuum calculation, exactly as in Figures 1-4.
101 55MMLL ++ 2355aa0 Note the smooth exponential convergence (without oscilla-
E| []h0 10-1 79MMLL ++ 1155aa000 ttihoensGiCn-eSlCecFtrcoansen.umberandfreeenergy)here,incontrastto
Φ -3
Φ - 10
| 10-5 fixed-potential case. On the other hand, with our q
κ
-7 modification, the GC-SCF algorithm at least converges
10
to the 10−6 E level independent of system size (Fig-
h
1
10 ure5),withthenumberofcyclesrequiredforconvergence
N|0 -1 remaining mostly unchanged with increasing number of
N - 10 Cu(111) layers, and with increasing thickness of the sol-
| 10-3 vent region.
Note that the standard fixed-charge SCF method con-
0 10 20 30 40
verges the energy smoothly and exponentially in most
GC-SCF cycles cases, including in our implementation in JDFTx. Our
implementation of the GC-SCF method in JDFTx uses
FIG. 5. Variation of GC-SCF convergence using the de-
exactlythesamecode, exceptforthepreconditionerand
fault parameters determined above with system size: both
metric modifications (due to q ) presented here. There-
with number of slab layers ranging from 5ML to 9ML with κ
fore we believe that the final convergence difficulty in
fixedvacuumspacing15a ,andwithvacuumspacingranging
0
GC-SCFisapropertyofthealgorithmitself,ratherthan
from 15 a to 35 a with the 5ML slab. Convergence slows
0 0
only marginally with increasing system size, with either vac- an implementation issue.
uum spacing or layer count. Each calculations is for solvated
Cu(111)chargedto1VSHEinCANDLEelectrolyte,starting
fromthestateofthecorrespondingconvergedneutralvacuum
C. Convergence of the GC-AuxH method
calculation.
The grand-canonical auxiliary Hamiltonian (GC-
continues to oscillate at that level beyond that point. AuxH) approach discussed in section IIB, directly min-
This final convergence issue may not affect practical imizes the total free energy of the system without as-
calculations where relevant energy differences are at the suming any models for physical properties of the system
10−3 E level or higher. However, smooth exponential (suchasdielectricresponsemodelsthatarebuiltintothe
h
convergencetothefinalanswerisdesirableasthismakes SCF mixing schemes). This algorithm contains a single
iteasiertodeterminewhenthetargetaccuracyhasbeen parameter K, which weights the relative contributions
reached. Unfortunately, no combination of GC-SCF pa- of the Kohn-Sham orbital and subspace Hamiltonian de-
rameters achieves uniformly smooth convergence for the grees of freedom in the conjugate-gradients search direc-