Bias-dependent diffusion of H 2 O molecules on an Al(111) surface
Satoshi Hagiwara, Chunping Hu, Satomichi Nishihara, Minoru Otani
aa r X i v : . [ c ond - m a t . o t h e r] D ec APS/123-QED
Bias-dependent diffusion of H O molecules on an Al(111) surface
Satoshi Hagiwara , ∗ Chunping Hu , Satomichi Nishihara , and Minoru Otani † National Institute of Advanced Industrial Science and Technology (AIST),1-1-1, Umezono, Tsukuba,Ibaraki 305-8568, Japan AdvanceSoft Corporation, 4-3,Kanda Suruga-dai, Chiyoda-ku,Tokyo 101-0062, Japan (Dated: December 21, 2020)
Abstract
We investigate the process by which a water molecule diffuses on the surface of an Al(111)electrode under constant bias voltage by first-principles density functional theory. To understandthe diffusion path of the water on the Al(111), we calculated the minimum energy path (MEP)determined by the nudged elastic band method in combination with constant electron chemicalpotential (constant- µ e ) methods. The simulation shows that the MEP of the water molecule,its adsorption site, and the activation barrier strongly depend on the applied bias voltage. Thisstrong dependence of the water diffusion process on the bias voltage is in good agreement with theresult of a previous scanning tunneling microscopy (STM) experiment. The agreement betweenthe theoretical and experimental results implies that accurate treatment of bias voltage plays asignificant role in understanding the interaction between the electric field and the surface of thematerial. Comparative studies of the diffusion process with the constant total number of electrons(constant- N e ) scheme show that the absence of strong interaction between the molecular dipole andthe electric field leads to a different understanding of how water diffuses on a metal surface. Theproposed constant- µ e scheme is a realistic tool for the simulation of reactions under bias voltagenot only using STM but also at the electrochemical interface. ∗ Electronic address: [email protected] † Electronic address: [email protected] . INTRODUCTION Fundamental studies on metal/water interfaces [1–4] have attracted much attention be-cause understanding the process of water reactions at the interface plays a central role ina wide variety of applications such as catalysis[5] and fuel cells[6]. A recent experimentto investigate the electrochemical interface reportedly showed that the molecular structureof water strongly depends on the electrode potential [7]. Furthermore, an experiment inwhich scanning tunneling microscopy (STM) was used found that the applied bias voltageaffected the activation barrier of water diffusion on a Pt surface [8]. To further understandthe water diffusion processes and reactions, it is necessary to clarify the extent to whichwater adsorption and diffusion on the surface depends on the bias voltage.First-principles density functional theory (DFT) [9, 10] is a powerful tool for investigat-ing the molecular adsorption and diffusion at the surface. However, previous theoreticalstudies on water diffusion processes were mainly carried out without the bias-voltage [11–14]. Generally, the electrode potential is related to the electron chemical potential ( µ e ) ofthe electrode, and the difference in µ e between the two subjects defines the bias voltage.Thus, for a system to which a fixed bias voltage is applied, the value of µ e of the electrode(or sample) surface must be constant. Therefore, to control the bias voltage in the simu-lation, we need to control µ e [15–18] beyond the constraint of a fixed number of electrons(constant- N e ).The quantum mechanical theory formulated under a grand-canonical ensemble is indis-pensable for a simulation with the fixed µ e (constant- µ e )[15–18]. Two flexible simulationmethods within the constant- µ e scheme have been proposed: The one is the fictitious chargeparticle (FCP) method developed by Bonnet et al. [19], and the other is the grand-canonicalself-consistent field (GCSCF) method introduced by Sundararaman et al [20]. These grand-canonical methods based on DFT can be applied not only to the electrochemical interface[20–24] but also to a system subjected to a bias voltage such as in the STM experiment.In this study, we carried out a first-principles study of the adsorption and diffusionprocesses a water molecule undergoes on the surface of Al(111) by varying the bias voltageusing the FCP and GCSCF methods. To understand the dependence of water diffusion onthe bias voltage, we computed the minimum energy path using the nudged elastic band(NEB) method [25–27] under constant bias voltage. Specifically, in the case of aluminum, it2 IG. 1: (Color online). (a) Schematic view of the model under a bias voltage. The cyan, red,and pink spheres represent Al, O, and H atoms, respectively. The coordinate z =0 is defined asthe center of the Al/H O system. Two ESM regions with dielectric constants ǫ =1 and ǫ = ∞ ,respectively, are attached to the left and right ends of the supercell to represent vacuum and thecounter electrode (metal). (b) Adsorption sites on the Al(111) surface. The cyan, green, and bluespheres represent the first-, second-, and third-layer Al atoms, respectively, counting from the sideof the counter electrode. is essential to understand the fundamental processes at the metal/water interface [12, 28, 29];thus, we used an aluminum electrode in our study.The paper is organized as follows: Section II provides a brief description of the basicideas of the constant- µ e schemes, their computational procedures, and the computationaldetails for the DFT and NEB calculations. In Sec. III, we discuss the results obtained bythe NEB combined with the constant- µ e scheme. Finally, the conclusions of our study arepresented in Sec. IV. II. METHODS AND COMPUTATIONAL DETAILS
Here, we describe the computational methods and details. First, we provide an overviewof the constant- µ e method under the boundary condition of the effective screening medium(ESM) technique [30]. The next subsection presents a discussion of the computationalprocedures of the FCP and GCSCF methods. Finally, we provide the computational detailsof the proposed method. 3 . Constant- µ e plus ESM method First, we briefly describe the basics of the constant- µ e scheme combined with the ESMmethod [30]. The ESM technique, which was developed by Otani and Sugino, is a powerfultool for studying various material surfaces under repeated slab approximation. Figure 1(a)shows a model of a constant- µ e calculation combined with the ESM used in this study. TwoESM regions with dielectric constants ǫ = 1 and ǫ = ∞ , respectively, are attached to thetwo ends of the supercell to represent the vacuum and the counter electrode (Metal). Theelectrostatic potential at the counter electrode was set to zero as the reference potential.Thus, the ESM enables us to always compare the energies measured from the same referencelevel.To achieve the constant- µ e condition, we use the FCP [19] and GCSCF methods [20].Both of these methods can be used compute the electronic structure and atomic geometryunder the given target chemical potential µ t . µ t is imposed by a potentiostat at a potentialΦ t , as shown in Fig. 1(a), i.e. µ t = − e Φ t , where e is the charge of an electron. We candefine the grand-potential Ω, instead of the total energy functional E tot , as follows:Ω = E tot − ( N e − N ) µ t = E tot − ∆ N e µ t , (1)where ∆ N e is a fictitious charge particle, which is the difference between the total number ofelectrons in the system N e and that in the neutral system N . Then, we minimize the valueof Ω for the atomic positions and ∆ N e , where ∆ N e is not a constant but a dynamic variableduring the entire minimization procedure. Because both the FCP and GCSCF methodsconverge to the same physical state, we expect these two methods to yield the same resultsunder the same computational conditions. B. Computational procedure of constant- µ e Here, we discuss the practical minimization procedure for Ω in the FCP and GCSCFmethods. Figure 2(a) and (b), respectively, show the flow charts for the calculations withFCP and GCSCF, where we show a series of flows for the DFT calculation with geometryoptimization. The SCF and geometry optimization loops show the blue and green shadedareas, respectively. In the following two subsections, we discuss the FCP and GCSCF4
IG. 2: (Color online). Calculation flows for (a) FCP and (b) GCSCF methods. The greenand blue shaded areas, respectively, indicate the calculation loops for the geometry optimizationand self-consistent field. Here, H KS , ε KS i , and ψ i denote the Kohn-Sham (KS) Hamiltonian, KS-eigenvalue, and KS-wavefunctions, respectively. The electron charge density is obtained by n ( r ) = P i f i | ψ i ( r ) | , where f i is the occupation number for each state i . F α denotes the forces acting oneach atom labeled by α . methods using these flows.
1. FCP method
The FCP employs a grand-canonical ensemble by the system connecting to the fictitiouspotentiostat, as shown in Fig. 1(a), and minimizes Ω under the constraint of constant- µ e in the loop in which the geometry is optimized, which is shown as the green shaded area.In the SCF loop shown as the blue shaded area, the Kohn-Sham equation is solved withthe fixed number of electrons. Therefore, the system reaches constant- µ e via simultaneouslyoptimizing not only the atomic positions but also the total number of electrons N e . To5ptimize N e , we define a fictitious force for N e as F e = − ∂ Ω ∂N e = − µ + µ t , (2)where, µ = ∂E tot /∂N e implies the instantaneous µ e , and yields the electrode potential as µ = − e Φ. To obtain N e for µ t , we minimize F e using the quasi-Newton algorithm byusing the Broyden-Fletcher-Goldfarb-Shanno minimization method [31–34]. This methodis commonly used for geometry optimization in conjunction with DFT calculations. In thegeometry optimization with the quasi-Newton algorithm, the following equation updates alldegrees of freedom for the atomic positions at the k -th iteration ( x k ), until all forces actingon the atoms become zero, as follows: x k +1 = x k + h k f k . (3)Here, f k and h k denote the forces acting on the atoms and Hessian, respectively. In con-ventional geometry optimization, the motion of N atoms in three dimensions produces 3 N degrees of freedom. In the FCP, both F e and N e are included in f k and x k , respectively.Thus, we explore the solution of eq. (3) within the space of the 3 N + 1-th dimension. How-ever, we cannot directly treat F e and N e on equal footing with f k and x k because the units of N e and F e are different from the atomic positions and forces. To address this, we introducethe effective charge position ( N ′ e = αN e ), where α is a scaling factor unit in bohr/ e . F e isalso scaled by α as follows, F ′ e = ( − µ + µ t ) /α. (4)Here, the definition of α is α = L max /V max C . L max and V max are the upper-bound of thechange in the length and voltage at each step, respectively. C is the capacitance determinedby the formula of the parallel-plate capacitor: C = 14 π SL . (5)Here, S denotes the surface area. Although the original definition of L is the distancebetween the parallel plates of the capacitor, we approximately use the half-length of theunit cell in the z-direction for convenience. In the FCP optimization, we add N ′ e and F ′ e to x k and f k , respectively.In eq. (3), h k plays a role in determining the step width of not only the new atomicpositions but also the new charge N e . For the Hessian component of N e , we use the first6erivative of µ with respect to the excess charge ( − ∂µ/∂N e ). In the present implementation,we use the approximate inverse of the density of states (DOS) at µ for − ∂µ/∂N e (1 /ρ ( µ )).Generally, for a large DOS system near µ t , µ gradually approaches µ t because of the smallHessian for N e evaluated by 1 /ρ ( µ ). Thus, the convergence behavior of the FCP methoddepends on the DOS near µ t .
2. GCSCF method
Here, we briefly review the GCSCF method discussed in Ref. 20. The GCSCF reachesconstant- µ e during the SCF loop, is shown as the blue shaded area in Fig. 2. Because theformulation of GCSCF is simple, its calculation flow is essentially the same as that of theconventional DFT with geometry optimization. However, in the GCSCF, N e is a variableat each SCF step. Generally, we can evaluate N e by summing the occupied Kohn–Sham(KS) orbitals with the given µ t . However, such a simple method for evaluating N e violatesthe numerical stability of the SCF [15]. Therefore, it is necessary to modify the numericalalgorithms to determine the i -th transient Fermi energy ( µ i ) and update the electron densityduring the SCF loop.Now, we explain the algorithm for determining the µ i . In the GCSCF, we graduallyapproach µ i to µ t during the SCF as follows: First, we evaluate the Fermi energy ε i F atthe i -th SCF step by the total number of electrons at the i -th SCF step, using an ordinalmethod. Second, µ i is determined by simply mixing ε i F and µ t as follows: µ i = βµ t + (1 − β ) ε i F , (6)where β is the mixing factor for µ e (0 < β < µ i is determined, and we finally update the occupation number and theelectron density using µ i and the total number of electrons.Next, we discuss the charge mixing scheme using the direct inversion of the iterativesubspace (DIIS) method [35] within the GCSCF framework. At the i -th SCF step, weupdate the electron density by solving the KS eq. with an input electron density n i in , andthen the updated electron density is used as a n i +1in . However, to obtain a more appropriatevalue of n i +1in , the updated electron density is mixed with the electron density of the previousSCF steps. To accelerate the convergence of the SCF, we usually use the optimized electron7ensity obtained by the solution of the DIIS method as n i +1in . Usually, to stabilize the DIISacceleration, the metric and Kerker preconditioning operators [36] ˆ M and ˆ K are introduced,and the DIIS method updates the electron density without altering the term of n ( G = ),where G is the reciprocal lattice vector. In contrast, the GCSCF requires the total numberof electrons to be updated in the SCF loop. Therefore, to update N e , we introduce thedumping factor of Q K ( >
0) to ˆ K and ˆ M as follows: h G | ˆ K | G i = | G | + Q | G | + q + Q , (7) h G | ˆ M | G i = 4 π | G | + Q . (8)Here, q K is the original damping factor of the Kerker preconditioning operator. When Q K is set to zero, ˆ K and ˆ M revert to their original values. By introducing Q K , ˆ K , and ˆ M at G = becomes a finite value, and we can always update the total number of electronsin the SCF loop. Once both Ω and µ i have converged, the ordinal geometry optimizationprocedure provides the new atomic positions.
3. NEB method combined with constant- µ e Here, we briefly discuss the NEB method in combination with the constant- µ e methods.The NEB method [25–27] describes the minimum energy path (MEP) for a chemical reactionby combining the images of the first and final states. The MEP is determined by solvingEq. (3) until the forces acting on each image become zero. Thus, in the conventional NEBmethod, we explore the solution of eq. (3) within the space of the 3 N × N im dimension,which corresponds to the motion of N atoms in three dimensions with N im images of theMEP. In the NEB combined with the FCP, eq. (3) is solved within the (3 N + 1) × N im dimension because we extended the optimization space, as discussed in the previous section.In contrast, when used in combination with the GCSCF, the NEB does not require specialmodification of the optimization procedure for the conventional NEB framework becauseof its straightforward formulation. In this work, we employ the Broyden method [37] as aquasi-Newton algorithm for determining the MEPs.8 . Computational details Here, we provide the computational details of this study. All calculations were performedusing the
Quantum Espresso package [38, 39], which is a DFT code within the plane-wavebasis sets and the ultrasoft-pseudopotential [40] framework. We implemented the FCP andGCSCF routines in combination with the ESM method in the package. We used the five-layered slab model in the p(2 ×
2) supercell to represent the Al(111) surface with a singlewater molecule, which corresponds to a 0.25 monolayer coverage, as shown in Fig. 1(a)].The two ESM regions with ǫ = 1 and ǫ = ∞ are located at a distance of ∼ . × . × .
05 ˚A ) [41] was used to construct the surface slab. The cut-off ener-gies for the wavefunctions and charge density were 40 Ry and 320 Ry, respectively. Theexchange-correlation functional was the spin unpolarized version of Perdew–Wang 91 withinthe generalized-gradient approximation [42]. k -point sampling used a 5 × × F α < . × − Ry/Bohr with the bottom of three Al layersfixed at bulk truncated positions. In the FCP calculation, the convergence thresholds for Ωand F e , respectively, are set to 1 . × − Ry and 1 . × − eV. In the GCSCF calculation,we decrease the threshold of the convergence criteria for Ω to 1 . × − Ry, and used theThomas–Fermi charge-mixing scheme[43]. The NEB [26, 27] calculation under the constant- µ e condition, which enables the determination of the MEP between two stable endpoints,was carried out to determine the activation barriers and diffusion paths of H O moleculesat the surface. The NEB calculation was conducted with ten discrete images for the MEPwith a path threshold of 0 .
03 eV / ˚A. The activation barriers and diffusion paths convergedwell for the number of images and path threshold. III. RESULTS AND DISCUSSIONS
Here, we discuss the results of water adsorption and diffusion on the Al(111) surface ascalculated using the NEB method with constant- µ e .9 . Adsorption site of water-molecule on Al(111) First, we briefly discuss the adsorption energy and µ e at the neutral Al(111) surface witha H O molecule for adsorption sites, as shown in Fig 1(b). The calculation shows that theon-top sites are the most stable for H O adsorption, and the adsorption energy we obtainedfor these sites is 0.23 eV. These results are consistent with the previous DFT result [12]. Forthe MEP, we considered H O diffusion from one stable on-top site to the next on-top sitevia a bridge site [44]. To apply the bias voltage U , we measure µ t from µ e at the potentialof zero charges (PZC), which was − . U is defined as U = ( µ t − µ PZC ) /e [45], and the three values of U are examined: U = 0V, − µ PZC is µ e at the PZC.For comparison purposes, we also carried out the constant- N e calculations by evaluatingthe excess charges (∆ N e ) induced by the bias voltage. The results we obtained for ∆ N e under U = −
1V and U = +1V are +0 . e and − . e , respectively. The difference inthe absolute values of ∆ N e between U = 1V and −
1V implies that the response of surfaceelectrons to the applied bias potential deviates from the linear response regime. In theconstant- N e calculation, we imposed the obtained ∆ N e during the diffusion processes. B. Activation barrier of water-diffusion on Al(111)
Table I presents the results of the activation barriers E a obtained by the NEB with theconstant- µ e and - N e schemes. First, we briefly compare the E a obtained by the FCP andGCSCF methods ( E FCPa and E GCSCFa ). The results of E FCPa are almost the same as thoseof E GCSCFa . Because the FCP and GCSCF methods converge to the same physical state,this agreement between the methods is reasonable. Hereafter, unless otherwise specified,we discuss the results of E a using those derived with FCP as being representative of theconstant- µ e scheme.The result of E a for the constant- N e with ∆ N e = 0 . e is very close to that of theprevious climbing image NEB[12]. In contrast, constant- µ e with U = 0 .
0V produces a muchlower E a compared to E a by the constant- N e with ∆ N e = 0 . e . Under U = − . E a increases to 0.161eV, which is also smaller than the counterpart value of 0.184 eV obtainedwith the constant- N e scheme with N e = +0 . e . By switching the value of U from 0 .
0V to10
ABLE I: The activation energies E a (in eV) of H O diffusion on the Al(111) under bias voltages of U = 0 . − . .
0V using the constant- µ e scheme. The results of E a for the constant- N e scheme with excess charges ∆ N e are also listed. The superscripts of E a by the constant- µ e scheme,respectively, denote the results obtained by the FCP and GCSCF methods.constant- µ e constant- N e U E
FCPa E GCSCFa ∆ N e E a . . e − . . e . − . e +1 . E a decreases from 0.092eV to 0.024eV. E a with U = +1 .
0V is still lower than thatdetermined by using the constant- N e method with ∆ N e = − . e . The above comparisonbetween the constant- µ e and - N e methods shows that although they yield similar trends for E a either as a result of applying bias voltages or introducing excess charges, the results arequantitatively quite different. Compared to the previous experimental data of H O-diffusion[8], our results of E a by the constant- µ e are in good agreement with the observation thatthe E a of water diffusion changes significantly with respect to the bias voltage. Thus,the simulation with the constant- µ e method plays an important role in reproducing theexperimental conditions for a constant bias-potential. C. MEPs for water-diffusion on Al(111)
Here, we discuss the results of the water diffusion path obtained by the NEB with theconstant- µ e and - N e schemes. Figure 3 (a)–(c), respectively, show consecutive images of thewater diffusion along the MEPs at applied bias voltage of U = 0V, +1 . − . O, we regard the first and last images as identical, and the water moleculesare bonded to the surface Al atom via the O atom. Along the diffusion path, the appliedbias voltage drastically alters the dipole direction of the water molecule. We can evaluatethe change in the direction of the H O dipole by obtaining the tilt angle between the dipolenormal and the surface ( θ ). For the first H O image, the values of θ under U = 0 . IG. 3: (Color online). Consecutive views of the MEP (ten images) for water diffusion on Al(111):results of constant- µ e scheme with bias-voltages of (a) 0.0V, (b) +1.0V, (c) -1.0V, and constant- N e scheme with the excess charge of (d) 0 . e , (e) − . e , and (f) +0 . e . The cyan, pink, andred spheres represent Al, H, and O atoms, respectively. The upper and lower panels show top andside views of the calculation cell, respectively. +1 . − .
0V are 76 ◦ , 85 ◦ , and 71 ◦ , respectively. This result indicates that the changesin the adsorption structure of H O resulting from the bias voltage are small for the firstimages of MEPs. However, the value of θ changes drastically in the intermediate imagesas a consequence of changes in the bias voltage. The dipole direction near the bridge siteat U = 0 .
0V and +1 .
0V becomes nearly perpendicular to the Al surface, and the H atomsorientate themselves downward. In contrast, the dipole directions for all images tend to beparallel to the surface at U = − . N e scheme.Here, the first and last images of H O in (d), (e), and (f) are the same as those in (a), (b),and (c), respectively. However, the changes in θ in the intermediate images are much smallerthan those in the constant- µ e . Thus, we interpret the difference in the results of E a betweenthe constant- µ e and - N e , as presented in Table I, as the difference in the MEPs. Therefore,this dependence of the H O geometry along the MEP on the applied bias voltage indicatesthe importance of the interaction between the water dipole and external electric fields.12
IG. 4: (Color online). Changes in the (a) grand-potential ∆Ω, (b) total energies ∆ E tot , (c)introduced excess electrons ∆ N e , and (d) chemical potential ∆ µ e as a function of the H O diffusionpath. The results of (a) and (c) are obtained under U = 0 . − . . N e = 0 . e , +0 . e , and − . e . The open circles and crosssymbols denote the results obtained by the FCP and GCSCF methods, respectively. The resultsof constant- N e are represented by diamond symbols. The solid and dashed lines are intended toguide the eyes. D. Details of H O diffusion along the MEPs
Figure 4 shows the results of the analysis with NEB. Before discussing the details of thediffusion properties, we briefly discuss the results between the FCP and GCSCF methodsshown in Fig. 4(a) and (c). Overall, the results of changes in the grand-potential ∆Ω andexcess charge ∆ N e along the MEPs obtained by the FCP (represented by open circles) andGCSCF (represented by cross symbols) are the same within the computational accuracy.These results indicate that we successfully implemented the constant- µ e methods. The maindifference between the FCP and GCSCF methods is the computational procedure discussedin Sec. II. Because, as discussed above, we used the extended Hessian in the FCP method, theconvergence behavior depends on the inverse of the DOS near µ t . In contrast, the GCSCFmethod directly optimizes µ e during a single SCF calculation. Because these constant- µ e µ e condition.Figure 4(a) presents the results of the ∆Ω, where the black, red, and blue circles, respec-tively, denote the values under U = 0 . − . . E a of the diffusion ofwater on Al(111) depends on the external electric field, as listed in Table I. For U = +1 . O molecules on Al(111) change from theon-top sites to sites in the vicinity of the bridge sites. In contrast, the results of ∆ E tot shown in Fig. 4(b) obtained by the constant- N e scheme do not alter the sign of ∆ E tot forany values of N e . Therefore, a comparison of the results of the constant- µ e and - N e schemeswould necessitate careful adjustment of the bias voltage to determine the stable adsorptionsite of H O on Al(111) in the STM experiments and at the electrochemical interface.Figure 4(c) shows the results of ∆ N e along the diffusion pathway. In the first images,the values of ∆ N e are the same as those used in the constant- N e calculations. In the nextfew images, the value of ∆ N e increases, and then it decreases in the intermediate images ofthe MEPs. This result is the consequence of introducing excess charges from an externalpotentiostat to Al(111) to maintain a constant bias voltage. The results obtained for ∆ µ e by using the constant- N e scheme highly depend on the MEPs shown in Fig. 4(d). Here, wedefine ∆ µ e as the difference between the µ e s in the MEPs and that in the first image of theMEPs of the PZC. Because the total number of electrons is fixed, this result originates fromthe charge transfer between the H O adsorbate and the Al electrode. This charge transferalters the height of the dipole barrier for the substrate along with the MEPs. Because thechange in the surface dipole barrier alters the work function that is directly related to theelectrode potential, ∆ µ e highly depends on the H O diffusion path. Thus, these differencesin the control mechanism of the surface charge between the constant- µ e and - N e schemesprovide the different MEPs for H O diffusion at Al(111), as shown in Fig. 3.14
V. SUMMARY
In summary, we demonstrated the realistic simulation of the bias-dependent diffusion ofH O on the Al(111) surface using NEB calculations within the constant- µ e scheme. Ourresults showed that significant differences exist in the activation barrier energies and MEPs,and that this depends on whether the applied bias voltage is positive or negative. A com-parison of the constant- µ e and - N e schemes also showed that the conventional constant- N e scheme does not provide a good description of molecular diffusion under a bias voltage owingto the absence of strong interaction between the molecular dipole and the electric field. Incomparison, the FCP and GCSCF methods produced the same results within the compu-tational accuracy. We expect the proposed scheme to find a wide variety of applications inthe simulation of STM experiments and electrochemical reactions under constant electrodepotentials. Acknowledgments
C.H. and M.O. thank Prof. Osamu Sugino for valuable discussions. This work was sup-ported by MEXT as the “Program for Promoting Research on the Supercomputer Fugaku”(Fugaku Battery & Fuel Cell Project), Grant Number JPMXP1020200301. The compu-tations were performed using the supercomputers of Research Center for ComputationalScience (Okazaki, Japan), and Institute for Solid State Physics and Information TechnologyCenter at the University of Tokyo. [1] M. A. Henderson, Surf. Sci. Rep. , 1 (2002).[2] A. Michaelides, V. Ranea, P. De Andres, and D. King, Phys. Rev. Lett. , 216102 (2003).[3] S. Meng, E. Wang, and S. Gao, Phys. Rev. B , 195404 (2004).[4] S. Schnur and A. Groß, New J. Phys. , 125003 (2009).[5] J. Carrasco, A. Hodgson, and A. Michaelides, Nat. Mater. , 667 (2012).[6] H. Ogasawara, B. Brena, D. Nordlund, M. Nyberg, A. Pelmenschikov, L. Pettersson, andA. Nilsson, Phys. Rev. Lett. , 276102 (2002).[7] T. Utsunomiya, Y. Yokota, T. Enoki, and K.-i. Fukui, Chem. Commun. , 15537 (2014).
8] K. Motobayashi, L. ´Arnad´ottir, C. Matsumoto, E. M. Stuve, H. J´onsson, Y. Kim, andM. Kawai, ACS Nano , 11583 (2014).[9] P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964).[10] W. Kohn and L. J. Sham, Phy. Rev. , A1133 (1965).[11] V. A. Ranea, A. Michaelides, R. Ram´ırez, P. L. de Andres, J. A. Verg´es, and D. A. King,Phys. Rev. Lett. , 136104 (2004).[12] V. A. Ranea, J. Chem. Phys. , 204702 (2012).[13] S. Karkare, L. Boulet, A. Singh, R. Hennig, and I. Bazarov, Phys. Rev. B , 035408 (2015).[14] T. B. Rawal, S. Hong, A. Pulkkinen, M. Alatalo, and T. S. Rahman, Phys. Rev. B , 035444(2015).[15] A. Lozovoi, A. Alavi, J. Kohanoff, and R. Lynden-Bell, J. Chem. Phys. , 1661 (2001).[16] I. Tavernelli, R. Vuilleumier, and M. Sprik, Phys. Rev. Lett. , 213002 (2002).[17] W. B. Schneider and A. A. Auer, Beilstein J. Nanotech. , 668 (2014).[18] U. Benedikt, W. B. Schneider, and A. A. Auer, Phys. Chem. Chem. Phys. , 2712 (2013).[19] N. Bonnet, T. Morishita, O. Sugino, and M. Otani, Phys. Rev. Lett. , 266101 (2012).[20] R. Sundararaman, W. A. Goddard III, and T. A. Arias, J. Chem. Phys. , 114104 (2017).[21] T. Ikeshoji and M. Otani, Phys. Chem. Chem. Phys. , 4447 (2017).[22] R. Sundararaman, K. Letchworth-Weaver, and K. A. Schwarz, The Journal of chemical physics , 144105 (2018).[23] J. Haruyama, T. Ikeshoji, and M. Otani, J. Phys. Chem. C , 9804 (2018).[24] S. E. Weitzner, S. A. Akhade, J. B. Varley, B. C. Wood, E. B. Duoss, S. E. Baker, andM. Otani, J. Phys. Chem. Lett. (2020).[25] G. Mills, H. J´onsson, and G. K. Schenter, Surface Science , 305 (1995).[26] G. Henkelman, B. P. Uberuaga, and H. J´onsson, J. Chem. Phys. , 9901 (2000).[27] G. Henkelman and H. J´onsson, J. Chem. Phys. , 9978 (2000).[28] A. Michaelides, V. Ranea, P. De Andres, and D. King, Phys. Rev. B , 075409 (2004).[29] J. Li, Y. Li, S. Zhu, and F. Wang, Phys. Rev. B , 153415 (2006).[30] M. Otani and O. Sugino, Phys. Rev. B , 115407 (2006).[31] C. G. Broyden, J. Inst. Appl. Math. , 222 (1970).[32] R. Fletcher, Comput. J. , 317 (1970).[33] D. Goldfarb, Math. Comput. , 23 (1970).
34] D. F. Shanno, Math. Comput. , 647 (1970).[35] P. Pulay, J. Comput. Chem. , 556 (1982).[36] G. P. Kerker, Phys. Rev. B , 3082 (1981).[37] C. G. Broyden, Math. Comput. , 577 (1965).[38] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L.Chiarotti, M. Cococcioni, I. Dabo, et al., J. Phys.: Condens. Matter , 395502 (2009).[39] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, M. Cococcioni, et al., J. Phys.: Condens. Matter , 465901 (2017).[40] D. Vanderbilt, Phys. Rev. B , 7892 (1990).[41] S. Popovi´c, B. Grˇzeta, V. Ilakovac, R. Kroggel, G. Wendrock, and H. L¨offler, Phys. StatusSolidi A , 273 (1992).[42] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, andC. Fiolhais, Phys. Rev. B , 6671 (1992).[43] D. Raczkowski, A. Canning, and L. W. Wang, Phys. Rev. B , 121101 (2001).[44] The first and the last images of the NEB calculation were prepared by the most stable adsorp-tion structure at a bias U . Then the intermediate images were initialized using the interpolationscheme both for the geometry and charge.[45] The sign of U is set as in the STM: A positive sample bias drives the electron injected intothe sample and increases the Fermi energy with respect to the µ pzc ..