Finite-temperature density functional calculation with polarizable continuum model in electrochemical environment

Available online at www.sciencedirect.com

Chemical Physics Letters 451 (2008) 158–162 www.elsevier.com/locate/cplett

Finite-temperature density functional calculation with polarizable continuum model in electrochemical environment Kazuya Shiratori a, Katsuyuki Nobusada a,b,* b

a Department of Structural Molecular Science, The Graduate University for Advanced Studies, Okazaki 444-8585, Japan Department of Theoretical and Computational Molecular Science, Institute for Molecular Science, Okazaki 444-8585, Japan

Received 9 October 2007; in final form 27 November 2007 Available online 4 December 2007

Abstract We present a numerical methodology to calculate electronic structures of a molecule in electrochemical environment. The methodology is based on the finite-temperature density functional theory (FTDFT) and allows us to study electronic properties of a molecule at a fixed chemical potential l. The approach is applied to a reaction of NOþ þ e ¢ NO in chemical equilibrium. The solvent effect is taken into account by a conductor-like polarizable continuum model (C-PCM). We demonstrate that the method combined with C-PCM (FTDFT/C-PCM) successfully describes the electronic structures of the molecule in electrochemical environment. Ó 2007 Elsevier B.V. All rights reserved.

1. Introduction Electrochemical processes have historically been investigated in a wide range of interests in electrochemical cell, corrosion, membrane potential, and very recently in photochemical cells [1]. To analyze those electrochemical processes in detail, we are required to clarify electronic structures of a system in electrochemical environment. The problems arising from modeling the electrochemical processes can be reduced to two parts. The first problem is difficulty in carrying out electronic structure calculations at a constant l, and the other one is how to appropriately describe the reactant–solvent and the reactant–electrode interactions. The conventional ab initio calculations are directed toward obtaining electronic structures at a constant number of electrons, N. Although several studies have been devoted to development of the methods calculating electronic structures at a chemical potential l [2–4], their meth* Corresponding author. Address: Department of Theoretical and Computational Molecular Science, Institute for Molecular Science, Okazaki 444-8585, Japan. Fax: +81 564 53 4660. E-mail address: [email protected] (K. Nobusada).

0009-2614/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2007.11.081

ods are still substantially based on the constant N calculations. Therefore, it is desirable to develop an alternative method to directly calculate electronic structures at a constant l [5–8]. Finite-temperature density functional theory (FTDFT) [9] treats a system in a grand canonical ensemble average and thus one can propose a numerical method based on FTDFT to describe electrochemical processes. In the present study, we develop a method based on the standard ab initio quantum chemistry approach by using the FTDFT calculation to discuss the electronic properties of a molecule in electrochemical environment. The standard ab initio supercell calculation incorporating with FTDFT has been carried out for an electrochemical system [5]. Nevertheless, we stress here that the quantum chemistry approach has two advantages over the supercell one: no need to consider the charge neutrality and an intuitive approach to understand chemical processes. We here mainly focus on developing a numerical method to calculate electronic structures at a constant chemical potential l. Specifically, the solvent effect is approximated by a conductorlike polarizable continuum model (C-PCM) [10,11], and we assume outer-sphere processes in which the reactant–electrode interaction is negligible. The FTDFT/C-PCM

K. Shiratori, K. Nobusada / Chemical Physics Letters 451 (2008) 158–162

calculation is applied to the reaction of NOþ þ e ¢ NO in chemical equilibrium. 2. Methods

NO+(g,298.15K,1mol/l) + e-(g,298.15K,1bar) ΔGbar mol/l

FTDFT is a generalization of DFT to a grand canonical ensemble first developed by Mermin [9] extending the Hohenberg–Kohn theorem [12]. In this theory, the grand potential X can be given by a functional form of electron density n as follows:

¼ i wi ðxÞ; 1 X 2 nðxÞ ¼ fj jwj ðxÞj ;

ð2Þ ð3Þ

j

1 þ exp



;

i l kB h

ΔGs (NO+)

NO+(aq,298.15K,1mol/l) + e-(g,298.15K,1bar)

ΔGox(NO) NO+(aq,298.15K,1mol/l) + e-(SHE,298.15K)

-eEabs,r(H+/H2) - e χ aq

Fig. 1. Born–Haber cycle of NO.

ð1Þ

where E is internal energy, h is temperature, and S is entropy. Similarly to the conventional Kohn–Sham (KS) DFT approach, the following finite temperature KS equation is derived from the minimum principle for the grand potential given at fixed temperature and chemical potential [13]   Z 1 2 nðx0 Þ 0  r þ dx þ vxc ½nðxÞ þ vðxÞ wi ðxÞ 2 jx  x0 j

1

ΔGs (NO+) NO(g,298.15K,1bar)

2.1. FTDFT calculation

fi ¼

NO+(g,298.15K,1bar) + e-(g,298.15K,1bar)

ΔGion(NO)

X½n ¼ E½n  hS½n  lN ½n;

159

ð4Þ

where kB is the Boltzmann’s constant, vxc ðxÞ is the exchange-correlation potential, v is the external potential, and wi , i , and fi are respectively the ith KS orbitals, KS orbital energies, and occupation numbers. Eq. (2) is formally equivalent to the conventional KS equation at zero temperature but the electron density nðxÞ depends on temperature through the fractional occupation numbers fi . Once one solves the finite temperature KS equation (Eq. (2)) satisfying Eqs. (3) and (4) in the same way as in the conventional KSDFT approach, the equilibrium electron density at l, h, and v is obtained. 2.2. Absolute electrode potential In the present calculation, the value of the reduced absolute electrode potential [14] for the standard hydrogen electrode (SHE), Eabs,r(H+/H2), is set equal to be 4.24 V vs vacuum. Then the electrode potential relative to SHE of an arbitrary electrode at electrochemical equilibrium EðOx=RedÞ is calculated by l EðOx=RedÞ ¼   4:24 V vs SHE; ð5Þ e where e is the elementary charge. The value, 4.24, is the same as the one by Camaioni et al. [15], but different from those by other authors, 4.34 [14] or 4.31 [16]. This difference arises from the variety for the thermochemical treatment of an electron and an uncertainty of the solvation

energy of the hydrogen ion. To calculate this value, we used a value of 361.7 kcal/mol for the formation free energy of hydrogen ion following the electron convention based on the Fermi–Dirac statistics [17,18], and a value of 264.0 kcal/mol for the solvation energy of the hydrogen ion not including the energy to across the solvent surface [19]. 2.3. Procedure to define the cavities The solvent effect is treated at the level of C-PCM in which the size of the cavity is a key ingredient to attain an appropriate description of solvation. The cavity is formed as an interlocked superposition of atomic spheres. The atomic radii of the neutral molecule NO were given in terms of the van der Waals values scaled by a factor of 1.2. For the charged molecule NOþ , the atomic radii were given so that the corresponding experimental solvation energy was reproduced. The solvation energy calculated by C-PCM does not include the energy to across the solvent surface, and corresponds to the one at the equal standard states both in gas phase and in solute phase [20]. Therefore, we calculated the solvation energy of NOþ ; DGs ðNOþ Þ  evaq ¼ 84:4 kcal=mol, by using the Born–Harber cycle shown in Fig. 1, where vaq is the surface potential of the aqueous solution. In this calculation, we used the values of DG ion ðNOÞ ¼ 213:69 kcal=mol [21,22], DG ox ðNOÞ ¼ 33:44 kcal=mol [23], and DGbar mol=l ¼ 1:90 kcal=mol. The radii were determined so as to reproduce the solvation energy, and they were interpolated as a linear function of the molecular charge rN ðqÞ ¼ 0:1904q þ 1:8600;

ð6Þ

rO ðqÞ ¼ 0:1867q þ 1:8240;

ð7Þ

where the ratio of the radii of N and O was fixed to be their van der Waals radii, 1.55:1.52 [24]. 3. Computational details We assume that vxc is represented by the Becke threeparameter hybrid exchange functional with the Lee– Yang–Parr correlation functional (B3LYP) [25,26]

160

K. Shiratori, K. Nobusada / Chemical Physics Letters 451 (2008) 158–162

Grand potential (eV)

2.5 2.0

E = -0.84 E = 1.16

1.5

E = 3.16 (V vs SHE)

1.0 0.5 0.0 1.0

Charge

0.8 0.6 0.4 0.2 0.0 1.00 Fig. 2. Flow chart of the computational procedure based on the present FTDFT/C-PCM calculation. l0 is the desired chemical potential and r is a set of effective radii for C-PCM. The dashed boxes indicate the additional numerical processes that are not required in the conventional DFT calculation for isolated molecules with a fixed N and in the vacuum.

although vxc in the FTDFT approach should depend on temperature in a narrow sense. The KS orbitals are expanded in Dunning’s augmented correlation-consistent basis set (aug-cc-pVDZ) [27,28]. The dielectric constant of water in C-PCM is set equal to be 78.39. The KSFTDFT approach mentioned in the previous section is similar to the conventional KSDFT one except for an extra procedure to specify the fractional occupation numbers fi . Therefore, we carried out the calculations with the GAMESS package of quantum chemistry programs [29] in which the numerical methodology of FTDFT is implemented. We preliminarily carried out the present FTDFT calculations by solving Eq. (2) self-consistently allowing the fractional occupation numbers given by Eq. (4). This numerical procedure showed, however, very slow convergence. This is because the total number of electrons N changes at each step of the self-consistent calculations. In other words, the strength to screen the nuclear charge by the electrons fluctuates, so that the KS orbital energies change heavily during the calculation. To avoid this problem, we solve Eq. (2) at a given l through finding the appropriate value of N with the bisection method. The flow chart of the computation is shown in Fig. 2. 4. Results and discussion 4.1. FTDFT/C-PCM calculation of NO in electrochemical environment Fig. 3a shows the grand potential curves for l ¼ 3:40; 5:40; and 7:40 eV. These values are equal to the electrode potentials E = 0.84, 1.16, and 3.16 V vs SHE,

1.04

1.08

1.12

1.16

Internuclear distance (Å) Fig. 3. (a) Grand potential curves of NO with different electrode potentials E. Each potential curve is shifted so that the energy of the bottom is equal to be zero and (b) charge of the system with different electrode potentials.

respectively. In what follows we will abbreviate ‘vs SHE’. The potential curves are shifted so that the grand potential at each equilibrium internuclear distance Req , (i.e., the bottom of the curve) is equal to be zero because their relative values have no meaning. The curves for E ¼ 0:84 and 3:16 V are approximately quadratic functions with Req ¼ ˚ , respectively. On the other hand, the curve 1:15 and 1.06 A for E ¼ 1:16 V does not show such quadratic dependence ˚ ) is equal to that for on R although its Req (=1.15 A E ¼ 0:84 V. To explain this E-dependence of the grand potential, the charge of the system is shown as a function of R in Fig. 3b. It is clearly seen from the figure that the charges for E = 0.84 and 3.16 V are 0.0 and 1.0, respectively in the whole range of R. This means that the grand potential curves for E = 0.84 and 3.16 V are identical to the potential energy curves of the neutral and cationic states, respectively, except for the term of lN in Eq. (1). The term of hS in Eq. (1) has a finite value. Nevertheless, in the present calculation, it makes no significant contribution to the grand potential because the energy gap is very large compared with the value of hS at h = 298.15 K. As distinct from the above two cases, the charge for E = 1.16 V depends on R. The charge is equivalent to either 1.0 (the cationic state) ˚ ) or longer or 0.0 (the neutral state) in the shorter (< 1.02 A ˚ (>1.12 A) internuclear distance, whereas N is a fractional number and the charge monotonically decreases as a func˚ < R < 1.12 A ˚ . This is tion of R in the range of 1.02 A ˚ ˚ because this range (1.02 A < R < 1.12 A) corresponds to a transition region between the grand potential curves of the neutral and cationic states. In other words, the system

K. Shiratori, K. Nobusada / Chemical Physics Letters 451 (2008) 158–162

is in a mixed state. Therefore, this figure indicates that the present FTDFT/C-PCM method succeeds in directly calculating the electronic states as a function of l. 4.2. Applicability of FTDFT/C-PCM calculation In this section, we further discuss the applicability of the FTDFT/C-PCM approach in comparison with the statistically averaged DFT one. In the statistically averaged DFT, b is given by the expectation value of any operator A XX b ¼ Tr½ AC b ¼ b N i i; h Ai pN i hWN i j AjW ð8Þ b¼ C

XX N

N

i

pN i jWN i ihWN i j;

ð9Þ

i

exp½ðEN i  lN Þ=ðk B hÞ pN i ¼ P P ; exp½ðEN i  lN Þ=ðk B hÞ N

ð10Þ

i

where EN i and jWN i i are respectively ith eigenvalue and eigenfunction in a pure state with N electrons. Fig. 4 shows the charge of NO in electrochemical environment as a function of the electrode potential E. The FTDFT result indicates that the charge transfer occurs in the range from E ’ 0:3 to 1.4 V while the statistically averaged DFT calculation shows discontinuous dependence on E. In the extension of the Hohenberg–Kohn theorem by Perdew et al. to the system with a fractional number of electrons N, an isolated open system was described by a statistical mixture, and they demonstrated that the energy calculated by using DFT should show the derivative discontinuity with respect to N [30]. Fig. 4 shows that the present FTDFT calculation fails to describe the derivative discontinuity. This is because the FTDFT approach directly addresses systems with a fractional number of electrons but the B3LYP functional is unsuitable for such systems. The functionals that represent the derivative discontinuity have been extensively developed [31,32]. Although the statistical averaged method gives the proper derivative discontinuity, Eq. (8) illustrates that the aver-

161

aged method becomes computationally demanding and might be infeasible when a number of states with different Ni contribute to the average. At this time, the computational time of the FTDFT calculation is more than four times longer than that of the statistically averaged DFT calculation. This is partly due to the fact that the numerical algorithm of the FTDFT calculation at a constant l is still at an early stage of development. 5. Concluding remarks We have developed a numerical method to study electronic structures of a molecule in electrochemical environment. The present method is based on the Mermin’s finitetemperature density functional theory (FTDFT), and the solvent effect was treated at the level of a conductor-like polarizable continuum model (C-PCM). The numerical methodology was successfully applied to a system of NOþ þ e ¢ NO in chemical equilibrium. We have presented that the FTDFT/C-PCM calculation gives a reasonable prediction for the electronic structure at a constant l of a molecule interacting with solvent molecules. In the present study, we focused on calculating electronic structure at a constant l. The solvent effect was simply modeled by C-PCM and the reactant–electrode interaction was not treated explicitly. However, the present FTDFT method has nothing to do with theoretical levels of description of these properties. Therefore, the method is straightforwardly extended to more real electrochemical processes by combining with sophisticated procedures describing the solvent effect and the reactant–electrode interaction. Furthermore, the present methodology can be applied to nonequilibrium electrochemical processes and thus it allows us to calculate reaction coordinate of charge transfer processes between electrodes and molecules as studied by Marcus [33]. This Letter is our starting point to establish a methodology to calculate molecular properties in electrochemical environment. Acknowledgements

Chemical potential (eV) -4.0

-5.0

-5.0

-7.0

-8.0

1.0

Charge

0.8

The present work was supported by a Grant-in-Aid (No. 18066019) and by the Next Generation Super Computing Project, Nanoscience Program, from MEXT, Japan. References

0.6 0.4 0.2 0.0 0.0

1.0

2.0

3.0

Electrode potential (V vs SHE) Fig. 4. Charge of NO in electrochemical environment as a function of the ˚. electrode potential E. The internuclear distance is fixed to be R = 1.04 A The solid and dash-dotted lines denote the results obtained by the FTDFT and statistically averaged DFT calculations, respectively.

[1] J.O’M. Bockris, S.U.M. Khan, Surface Electrochemistry, Plenum, New York, 1993. [2] H. Nakatsuji, J. Chem. Phys. 87 (1987) 4995. [3] I. Tavernelli, R. Vuilleumier, M. Sprik, Phys. Rev. Lett. 88 (2002) 213002. [4] Y. Cai, A.B. Anderson, J. Phys. Chem. B 108 (2004) 9829. [5] A.Y. Lozovoi, A. Alavi, J. Kohanoff, R.M. Lynden-Bell, J. Chem. Phys. 115 (2001) 1661. [6] S. Jacobi, R. Baer, J. Chem. Phys. 123 (2005) 044112. [7] M. Otani, O. Sugino, Phys. Rev. B 73 (2006) 115407. [8] C.D. Taylor, S.A. Wasileski, J.-S. Filhol, M. Neurock, Phys. Rev. B 73 (2006) 165402.

162 [9] [10] [11] [12] [13] [14] [15] [16] [17]

[18] [19] [20] [21]

K. Shiratori, K. Nobusada / Chemical Physics Letters 451 (2008) 158–162 N.D. Mermin, Phys. Rev. 137 (1965) A1441. V. Barone, M. Cossi, J. Phys. Chem. A 102 (1998) 1995. J. Tomasi, B. Mennucci, R. Cammi, Chem. Rev. 105 (2005) 2999. P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864. W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. J. Llano, L.A. Eriksson, J. Chem. Phys. 117 (2002) 10193. D.M. Camaioni, M. Dupuis, J. Bentley, J. Phys. Chem. A 107 (2003) 5778. S. Trasatti, J. Chem. Soc., Faraday Trans. 1 70 (1974) 1752. S.G. Lias, J.E. Bartmess, Gas-phase ion thermochemistry in NIST chemistry WebBook, NIST standard reference database number 69, in: P.J. Linstrom, W.G. Mallard, June 2005, National Institute of Standards and Technology, Gaithersburg MD, 20899 . J.E. Bartmess, J. Phys. Chem. 98 (1994) 6420. M.D. Tissandier et al., J. Phys. Chem. A 102 (1998) 7787. A. Ben-Naim, J. Solution Chem. 30 (2001) 475. P.J. Linstrom, W.G. Mallard (Eds.), NIST Chemistry WebBook, NIST Standard Reference Database Number 69, National Institute of

[22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

Standards and Technology, Gaithersburg MD, 20899, June 2005, . D.R. Lide (Ed.), CRC Handbook of Chemistry and Physics, 86th ed., CRC, Boca Raton, 2005. A.J. Bard (Ed.), Encyclopedia of Electrochemistry of the Elements, vol. 8, Marcel Dekker, New York, 1978, p. 325. A. Bondi, J. Phys. Chem. 68 (1964) 441. A.D. Becke, J. Chem. Phys. 98 (1993) 5648. C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. T.H. Dunning Jr., J. Chem. Phys. 90 (1989) 1007. R.A. Kendall, T.H. Dunning Jr., R.J. Harrison, J. Chem. Phys. 96 (1992) 6796. M.W. Schmidt et al., J. Comput. Chem. 14 (1993) 1347. J.P. Perdew, R.G. Parr, M. Levy, J.L. Balduz Jr., Phys. Rev. Lett. 49 (1982) 1691. J.B. Krieger, Y. Li, G.J. Iafrate, Phys. Rev. A 45 (1992) 101. O.A. Vydrov, G.E. Scuseria, J.P. Perdew, J. Chem. Phys. 126 (2007) 154109. R.A. Marcus, J. Chem. Phys. 24 (1956) 966.