Implicit self-consistent description of electrolyte in plane-wave density-functional theory

0.97 MB·

by  Kiran Mathew

Implicit self-consistent description of electrolyte in plane-wave density-functional
theory
Kiran Mathew1,2,∗ and Richard G. Hennig2,†
1Department of Materials Science, Cornell University, Ithaca, New York 14853, USA
2Department of Materials Science, University of Florida, Gainesville, Florida 32608, USA
(Dated: January 14, 2016)
Ab-intiocomputationaltreatmentofelectrochemicalsystemsrequiresanappropriatetreatmentof
thesolid/liquidinterfaces. Afullyquantummechanicaltreatmentoftheinterfaceiscomputationally
unfeasible due to the large number of degrees of freedom involved. In this work we describe a
computationally efficient model where the electrode part of the interface is modeled at the density
6
functionaltheory(DFT)levelandtheelectrolytepartisrepresentedthroughanimplicitmodelbased
1
onthePoisson-Boltzmannequation. WedescribetheimplementationofthemodelViennaAb-intio
0
Simulation Package (VASP), a widely used DFT code, followed by validation and benchmarking of
2
theimplementation. Todemonstratetheutilityoftheimplicitelectrolytemodelweapplythemodel

n to study the effect of electrolyte and external voltage on the surface diffusion of sodium atoms on
a the electrode.
J

3 PACSnumbers:
1

] I. INTRODUCTION 100
i ε
c n diel
s Interfacesbetweendissimilarsystemsareubiquitousin el
- 50
Solute Electrolyte
l nature and frequently encountered and employed in sci-
r
t entific and engineering applications. Of particular inter-
m
estaresolid/liquidinterfacesthatformthecruxofmany 0
. technologically important systems such as electrochem-
t n
a ical interfaces,1–3 nanoparticle synthesis,4–6 and protein ion
m –50
folding and membrane formation.7 A comprehensive un-
Inter-
- derstanding of the behavior and properties of such sys-
nd tems requires a detailed microscopic description of the –100 VHartree face
o solid/liquid interfaces. 25 30 35 40
c Given the heterogeneous nature and the complexities Position (Å)
[ ofsolid/liquidinterfaces,8 itisoftenthecasethatexper-

1 imental studies need to be supplemented by theoretical FIG. 1: Spatial decomposition of the solid/electrolyte sys-
v undertakings. A complete theoretical study of such in- tem into the solute, interface, and electrolyte regions for a
6 terface based onan explicit descriptionof the electrolyte bccNa(110)surfaceembeddedinanelectrolytewithrelative
4 andsoluteinvolvessolvingsystemsofnonlinearequations permittivity (cid:15) = 78.4 and monovalent cations and anions
r
3 withenormousamountofdegreesoffreedomandquickly with concentration 1M. The solute is described by density-
3 become prohibitively expensive even when state of the functional theory, the electrolyte by an implicit solvation
0 art computational resources are employed.9 This calls model. The interface regions is formed self-consistently as
1. for the development of efficient implicit computational a functional of the electronic density of the solute. All prop-
0 frameworks for the study of solid/liquid interfaces.9–14 erties are shown as a percentage of their maximum absolute
6 value.
We have previously developed a self-consistent solva-
1
tion model15,16 that describes the dielectric screening of
:
v a solute embedded in an implicit solvation model for
i
X the solvent, where the solute is described by density- in the solvent through the solution of the linear Poisson-
functional theory (DFT). We implemented this solva- Boltzmann equation and then apply the method to elec-
r
a tionmodelintothewidelyusedDFTcodeVASP.17 This trochemical systems to illustrate the power of this ap-
software package, VASPsol, is freely available under the proach. Figure1illustrateshowthesolid/electrolytesys-
GPL v3 license and has been used by us and others to tem is divided into spatial regions that are described by
study the effect of the presence of solvent on a number the explicit DFT and the implicit solvation model. This
of materials and processes, such as the catalysis of ke- newsolvationmodelallowsforthedescriptionofcharged
tonehydrogenation,18theinterfacialchemistryinelectro- systems and the study of the electrochemical interfacial
chemical interfaces,19 hydrogen adsorption on platinum systems under applied external voltage. Section II de-
surfaces,20 etc. scribes the formalism and outlines the derivation of the
Inthisworkwedescribetheextensionofoursolvation energy expression. In Sec. III we describe the details of
model VASPsol15,16 to include the effects of mobile ions theimplementationandvalidatethemethodagainstpre-
2
vious computational studies of electrochemical systems. where A is the kinetic and exchange-correlation con-
TXC
Finally, in Sec. IV we apply the method to study the ef- tributionfromDFT,φisthenetelectrostaticpotentialof
fect of an applied external potential on the adsorption thesystem,andρ andρ arethetotalchargedensityof
s ion
and surface diffusion for metal electrodes. thesoluteandtheionicchargedensityoftheelectrolyte,
respectively. Thenetsolutechargedensity,ρ ,isthesum
s
ofthesoluteelectronicandnuclearchargedensities,n((cid:126)r)
II. THEORETICAL FRAMEWORK and N((cid:126)r), respectively,
ρ ((cid:126)r)=n((cid:126)r)+N((cid:126)r). (3)
Our solvation model couples an implicit description of s
theelectrolytetoaquantum-mechanicalexplicitdescrip- The net ionic charge density, ρ , is given by
ion
tion of the solute. The implicit solvation model incor-
(cid:88)
porates the dielectric screening due to the permittivity ρ ((cid:126)r)= qz c ((cid:126)r), (4)
ion i i
of the solvent and the electrostatic shielding due to the
i
mobile ions in the electrolyte. The solute is described
explicitly using DFT. where ci((cid:126)r) is the concentration of ionic species i with
To obtain the energy of the combined so- charge ziq, and q is the elementary charge. The con-
lute/electrolytesystem,wespatiallydividethematerials centration of ionic species, ci, is given by a Boltzmann
system into threeregions: (i)the solute thatwe describe factoroftheelectrostaticenergy24 andmodulatedinthe
explicitlyusingDFT,(ii)theelectrolytethatwedescribe interface region by the shape function, S(n((cid:126)r)),
usingthelinearPoisson-Boltzmannequationand(iii)the (cid:18) (cid:19)
z qφ((cid:126)r)
interface region, that electrostatically couples the DFT c ((cid:126)r)=S[n((cid:126)r)]c0exp − i , (5)
i i kT
and Poisson-Boltzmann equation. Fig. 1 illustrates the
regionsforaNa(110)surfaceembeddedina1Maqueous
where c0 is the bulk concentration of ionic species
electrolyte. The interface between the electrolyte and i
i, k is the Boltzmann constant and T the tempera-
the solute system is formed self-consistently through
ture. The relative permittivity of the electrolyte, (cid:15)((cid:126)r),
the electronic density and the interaction between the
is assumed to be a local functional of the electronic
electrolyte and solute is described by the electrostatic
charge density of the solute and modulated by the shape
potential as explained in the following sections.
function,11,15,21–23
(cid:15)(n((cid:126)r))=1+((cid:15) 1)S(n((cid:126)r)), (6)
b
A. Total energy −
where (cid:15) is the bulk relative permittivity of the solvent.
b
We define the interface region between the solute and The cavitation energy, Acav, describes the energy re-
the electrolyte using the electronic charge density of the quired to form the solute cavity inside the electrolyte
solute, n((cid:126)r). The electrolyte occupies the volume of the and is given by14
simulation cell where the solute charge density goes to (cid:90)
zero. In the interface region where the electronic charge A =τ S d3r. (7)
cav
|∇ |
density increases rapidly from zero towards the average
bulk value of the solute, the properties of the implicit In comparison to our previous solvation model,15 the
solvation model, e.g. the permittivity (cid:15) and the Debye
main changes to the free energy expression are the in-
screeninglengthλ ,changefromtheirbulkvaluesofthe
D clusionoftheelectrostaticinteractionoftheioniccharge
electrolyte to the vacuum values in the explicit region.
density in the electrolyte with the system’s electrostatic
The scaling of the electrolyte properties is given by a
potential, φ, and the non-electrostatic contribution to
shape function11,15,21–23
the free energy from the mobile ions in the electrolyte,
A . Inthisworkweassumethattobejusttheentropic
1 (cid:26)log(n/n )(cid:27) ion
S(n((cid:126)r))= erfc c . (1) contribution,24 i.e.
2 σ√2
A =kTS , (8)
ion ion
Following our previous work15 and the work by Arias
et al.,11,21–23 the total free energy, A, of the system con- where Sion is the entropy of mixing of the ions in the
sisting of the solute and electrolyte is expressed as electrolyte, which for small electrostatic potentials can
be approximated to first order as
(cid:90)
A[n((cid:126)r),φ((cid:126)r)]=ATXC[n((cid:126)r)]+ φ((cid:126)r)ρs((cid:126)r)d3r S =(cid:90) (cid:88)c ln(cid:18)ci(cid:19)d3r
(cid:90) φ2 (cid:90) 1 ion i i c0i
− (cid:15)((cid:126)r)|∇8π| d3r+ 2φ((cid:126)r)ρion((cid:126)r)d3r (cid:90) (cid:88)c ziqφd3r. (9)
+A +A , (2) ≈− i kT
cav ion i
3
B. Minimization 5.5
Cu
To obtain the stationary point of the free energy, A, )
V5.0 Au
given by Eq. (2), we set the first order variations of the (
C
free energy with respect to the system potential, φ, and Z Ag
P
theelectronicchargedensity,n(r),tozero.11,15,21–23Tak- al4.5
ingthevariationofA[n((cid:126)r),φ((cid:126)r)]withrespecttotheelec- tic
e
tronicchargedensity,n((cid:126)r),yieldsthetypicalKohn-Sham r
o
Hamiltonian25 with the following additional term in the he4.0
T
local part of the potential
Vsolv =δ(cid:15)δ(nn)|∇8φπ|2 +φδρδnion 3−.50.8 −0.6 E−xp0e.r4ime−nt0a.l2PZC0(.0V,vsS0.H2E) 0.4 0.6
δ S δS
+τ |∇ | +kT ion. (10)
δn δn FIG. 2: Comparison between the computed and the exper-
imental Potential of Zero Charge(PZC) with respect to the
TakingthevariationofA[n((cid:126)r),φ((cid:126)r)]withrespecttoφ((cid:126)r) standard hydrogen electrode (SHE). The dashed line is a fit
yields the generalized Poisson-Boltzmann equation ofVpred =Vexp+Vpredtodeterminethetheoreticalpotential
SHE
of the SHE.
(cid:126) (cid:15)(cid:126)φ= ρ ρ , (11)
s ion
∇· ∇ − −
where (cid:15)(n(r)) is the relative permittivity of the solvent pseudopotentials26,27 as well as the PAW28 formulation
as a functional of the electronic charge density. of pseudopotentials.
Wefurthersimplifythesystembyconsideringthecase The main modifications in the code are the evalua-
ofelectrolyteswithonlytwotypesofionspresent,whose tion of the additional contributions to the total energy
chargesareequalandopposite,i.e. c01=c02=c0andz1= and the local potential given by Eqs. (2) and (10), re-
z2 =z. Thentheionicchargedensityoftheelectrolyte spectively. Corrections to the local potential require the
−becomes24 solution of the generalized Poisson-Boltzmann equation
(cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21) (13) in each self-consistent iteration. Our implementa-
zqφ zqφ
ρ =qzc0 exp − exp tion solves the equation in reciprocal space and makes
ion kT − kT efficient use of fast Fourier transformations (FFT).15
(cid:18) (cid:19)
zqφ This enables our implementation to be compatible with
= 2qzc0sinh (12)
− kT the Message Passing Protocol (MPI) and to take ad-
vantage of the memory layout of VASP. We use a pre-
and the Poisson-Boltzmann equation becomes conditioned conjugate gradient algorithm to solve the
generalized Poisson-Boltzmann equation with the pre-
(cid:18) (cid:19)
(cid:126) (cid:15)(cid:126)φ= ρ +2qzc0sinh zqφ . (13) conditioner (cid:0)G2+κ2(cid:1)−1, where G is the magnitude of
∇· ∇ − s kT
the reciprocal lattice vector and κ2 is the inverse of the
square of the Debye length, λ . Our software mod-
For small arguments x = zqφ 1, sinh(x) x and we D
kT (cid:28) → ule, VASPsol, is freely distributed as an open-source
obtain the linearized Poisson-Boltzmann equation
packageandhostedonGitHubathttps://github.com/
henniggroup/VASPsol.16
(cid:126) (cid:15)(cid:126)φ κ2φ= ρ (14)
∇· ∇ − − s We benchmark the model against existing experimen-
tal and computational data.21 First, we compute the po-
with
tentialofzerocharge(PZC)ofvariousmetallicslabsand
(cid:18)2c0z2q2(cid:19) 1 compare them against the experimental values. Second,
κ2 = = , (15)
kT λ2 we compute the surface charge density of a Pt(111) slab
D
asfunctionoftheappliedexternalpotentialandcompare
where λ is the Debye length that characterizes the di- the resulting value of the double layer capacitance with
D
mension of the electrochemical double layer. previous computation and experiment.
The DFT calculations for these two benchmarks are
performed with VASP and the VASPsol module using
III. IMPLEMENTATION AND VALIDATION projected augmented wave (PAW) formalism describ-
ing the electron-ion interactions and the PBE approxi-
The implicit electrolyte model described above is im- mation for the exchange-correlation functional.17,29 The
plemented in the widely used DFT software Vienna Brillouin-zone integration uses an 8 8 1 k-points. A
Ab-intio Software Package (VASP).17 VASP is a paral- high plane-wave basis set cutoff ener×gy o×f 800 eV is used
lel plane-wave DFT code that supports both ultra-soft to ensure proper resolution of the interfacial region be-
4
30
)
2
m
c
/ 20
C
µ
(
σ 10
e, FIG. 4: Pt adatom(black) adsorbed on three different sites
g
har 0 on the Pt 111 surface. Red, green and blue indicate top,
c second and third layers of Pt 111 respectively
e
c
urfa−10
S
5 A
0.0 0.5 1.0 1.5 2.0 2.5
B
Electrodepotential(V,vsSHE)
) 0 C
V
FIG. 3: Surface charge density of the Pt(111) surface as a (e
functionofappliedelectrostaticpotentialcomputedwiththe ds
Ea
implicit electrolyte module, VASPsol. The slope of the curve 5
determines the capacitance of the dielectric double layer. −
10
−
tweenthesoluteandimplicit. Weconsideranelectrolyte 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.5
− − − −
thatconsistsofanaqueoussolutionofmonovalentanions ElectrodepotentialvsSHE(V)
and cations of 1M concentration. At room temperature
this electrolyte has a relative permittivity of (cid:15) = 78.4 FIG. 5: Plot of voltage vs the adsorption energy
b
andaDebyelengthof3˚A.Tochangetheappliedpoten-
tial,weadjusttheelectroncountandthendeterminethe
correspondingpotentialfromtheshiftinelectrostaticpo- capacitanceandthecomputedvalueintheliteraturecan
tential as reflected in the shift in Fermi level. This takes be attributed to the difference in pseudo-potentials.
advantageofthefactthattheelectrostaticpotentialgoes
to zero in the electrolyte region for the solution of the
Poisson-Boltzmann equation.
IV. APPLICATIONS
Figure 2 compares the electrode potential of zero
charge (PZC) calculated with the implicit electrolyte
The implicit electrolyte model, VASPsol, enables the
model as implemented in VASPsol with experimental
studytheelectrodesurfacesinrealisticenvironmentsand
data. The electrode potential of zero charge is defined
under various conditions, such as of an electrode im-
astheelectrostaticpotentialofaneutralmetalelectrode
mersed in an electrolyte with an external potential ap-
and is given by the Fermi energy with respect to a ref-
plied. To illustrate the utility of the solvation model, we
erence potential. A natural choice of reference for the
applyittostudytheadsorptionenergyanddiffusionen-
implicitelectrolytemodelistheelectrostaticpotentialin
ergy barrier of an adatom on a metallic electrode that is
the bulk of the electrolyte, which is zero in any solution
immersed in a 1M aqueous electrolyte as a function of
of the Poisson-Boltzmann equation. In experiments, the
applied potential.
reference is frequently chosen as the standard-hydrogen
Stability of fuel cell materials plays a significant
electrode (SHE). The dashed line in Fig. 2 presents the
role in the the development of polymer electrolyte fuel
fit of Vpred = Vexp +Vpred to obtain shift between the
SHE cells(PEFCs). PEFC performance loss under steady-
experimentalPZCandthecomputedPZC,. i.e. thepo-
state and cycling conditions has been attributed in part
tential of the SHE for our electrolyte model. The result-
to a loss of electrochemically active surface area of the
ingVpred =4.6Vcomparesverywellwiththepreviously platinum cathode electro-catalyst31. In this work we
SHE
reported computed value of 4.7 V for a similar solvation study the dissolution of Pt electrode by looking at the
model.23 adsorption energies of Pt adatom on Pt 111 surfaces.
Figure 3 shows the calculated surface charge den- We apply the implicit electrolyte model to simulate the
sity, σ, of the Pt(111) surface as a function of applied Pt electrode under external potential. For the adsorp-
electrostatic potential, V. The linear slope of σ(V) is tion, three adsorption sites are considered, labeled A, B
the double-layer capacitance. The resulting double-layer and C, as shown in Figure 4. The plot of adsorption
capacitance for the Pt(111) surface at 1M concentra- energy as a function of the externally applied potential
tion is 15 µF/cm2, which closely agrees with the pre- is shown in the Figure 5. From the plot it can be ob-
viously reported computed23 and experimental30 values servedthattheadsorptionenergyincreasesasthepoten-
of 14 µF/cm2 and 20 µF/cm2, respectively. The dis- tial is ramped up. This implies that the Pt adatom is
crepancy between the VASPsol value of the double layer less likely to adsorb on to the Pt 111 surface at higher
5
potentials i.e rate of dissolution is higher at larger po- by comparing the calculated electrode potential of zero
tentials. Theobservationiscorroboratedbypreviousex- chargeforvariousmetalelectrodeswithexperimentalval-
perimentalstudies31,32,thusdemonstratingthepractical ues and by comparing the calculated double-layer capac-
applicability of the model. itance of Pt(111) in a 1M solution and to previously cal-
culations and experiments. To illustrate the usefulness
of the model and the implementation, we compute the
V. CONCLUSIONS adsorption energy of Pt adatoms on a Pt(111) electrode
asfunctionoftheappliedpotential. Confirmingprevious
Solid/liquid interfaces between electrodes and elec- experimental observations, we find that the Pt dissolu-
trolytes in electrochemical cells present a complex sys- tionrateincreasesasafunctionoftheexternallyapplied
tem that requires the insights provided by computa- potential.
tional studies to supplement and explain experimental
observations. We developed and implemented a self-
consistent computational framework that provides an ef- Acknowledgments
ficient implicit description of electrode/electrolyte inter-
faces within density-functional theory through the solu- This work was supported by the National Science
tion of the generalized linear Poisson-Boltzmann equa- Foundation under the CAREER award No. DMR-
tion. Theelectrolytemodelisimplementedinthewidely 1056587 and under award No. ACI-1440547. This re-
used DFT code VASP and the implementation is made search used computational resources of the Texas Ad-
availableasafreeandopen-sourcemodule,VASPsol. We vanced Computing Center under Contract Number TG-
show that this model enables ab-initio studies of elec- DMR050028N and of the University of Florida Research
trochemistry at the DFT level. We validate the model Computing Center.
∗ Electronic address: [email protected] Physics 140, 084106 (2014).
† Electronic address: [email protected]fl.edu 16 K. Mathew and R. G. Hennig, VASPsol - Solvation model
1 D.Wang,H.L.Xin,Y.Yu,H.Wang,E.Rus,D.A.Muller, fortheplanewaveDFTcodeVASP,https://github.com/
and H. D. Abrua, J. Am. Chem. Soc. 132, 17664 (2010). henniggroup/VASPsol (2015).
2 D.-H. Ha, M. A. Islam, and R. D. Robinson, Nano Lett. 17 G. Kresse and J. Furthmu¨ller, Comput. Mater. Sci. 6, 15
12, 5122 (2012). (1996).
3 G. Liu, E. Luais, and J. J. Gooding, Langmuir 27, 4176 18 C.Michel,J.Zaffran,A.M.Ruppert,J.Matras-Michalska,
(2011). M. Jedrzejczyk, J. Grams, and P. Sautet, Chemical Com-
4 J.J.Choi,C.R.Bealing,K.Bian,K.J.Hughes,W.Zhang, munications 50, 12450 (2014).
D.-M.Smilgies,R.G.Hennig,J.R.Engstrom,andT.Han- 19 N. Kumar, K. Leung, and D. J. Siegel, Journal of The
rath, J. Am. Chem. Soc. 133, 3131 (2011). Electrochemical Society 161, E3059 (2014).
5 D. F. Moyano and V. M. Rotello, Langmuir 27, 10376 20 S. Sakong, M. Naderian, K. Mathew, R. G. Hennig, and
(2011). A. Groß, The Journal of chemical physics 142, 234107
6 T.-F. Liu, Y.-P. Chen, A. A. Yakovenko, and H.-C. Zhou, (2015).
Journal of the American Chemical Society 134, 17358 21 K. Letchworth-Weaver and T. A. Arias, Phys. Rev. B 86,
(2012), pMID: 23043458. 075140 (2012).
7 K. A. Dill, S. Bromberg, K. Yue, K. M. Fiebig, D. P. Yee, 22 R. Sundararaman, D. Gunceler, K. Letchworth-
P. D. Thomas, and H. S. Chan, Protein Science: A Publi- Weaver, and T. A. Arias. (2012), JDFTx.
cation of the Protein Society 4, 561602 (1995). http://jdftx.sourceforge.net.
8 V.A.AgubraandJ.W.Fergus,JournalofPowerSources 23 D. Gunceler, K. Letchworth-Weaver, R. Sundararaman,
268, 153 (2014), ISSN 0378-7753. K. A. Schwarz, and T. A. Arias, Modelling Simul. Mater.
9 S.Ismail-BeigiandT.A.Arias,Comp.Phys.Comm.128, Sci. Eng. 21, 074005 (2013).
1 (2000). 24 K. A. Sharp and B. Honig, pp. 7684–7692 (1990).
10 J. A. White, E. Schwegler, G. Galli, and F. Gygi, Journal 25 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
of Chemical Physics 113 (2000). 26 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
11 S. A. Petrosyan, A. A. Rigos, and T. A. Arias, J. Phys. 27 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
Chem. B 109, 15436 (2005). 28 P. E. Bl¨ochl, Phys. Rev. B 50, 17953 (1994).
12 C. J. Cramer and D. G. Truhlar, Accounts of Chemical 29 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.
Research 41, 760 (2008). Lett. 77, 3865 (1996).
13 J.LischnerandT.A.Arias,J.Chem.Phys.B114, 1946 30 A. Tymosiak-Zieliska and Z. Borkowska, Electrochimica
(2010). Acta 46, 3063 (2001), ISSN 00134686.
14 O. Andreussi, I. Dabo, and N. Marzari, J. Chem. Phys 31 X.Wang,R.Kumar,andD.J.Myers,Electrochemicaland
136, 064102 (2012). Solid-State Letters 9, A225 (2006), ISSN 10990062.
15 K. Mathew, R. Sundararaman, K. Letchworth-Weaver, 32 P. J. Ferreira, G. J. la O, Y. Shao-Horn, D. Morgan,
T. A. Arias, and R. G. Hennig, The Journal of Chemical R. Makharia, S. Kocha, and H. A. Gasteiger, Journal
6
of The Electrochemical Society 152, A2256 (2005), ISSN 00134651.