Assessmentofdensityfunctionalapproximationsforthehemibondedstructureofthewaterdimerradicalcation PAPER PCCP

(1)

Cite this: Phys. Chem. Chem. Phys ., 2012, 14, 10705–10712

Assessment of density functional approximations for the hemibonded structure of the water dimer radical cation

Piin-Ruey Pan,

^a

You-Sheng Lin,

^a

Ming-Kang Tsai,

^b

Jer-Lai Kuo*

^c

and Jeng-Da Chai*

^ad

Received 6th April 2012, Accepted 6th June 2012 DOI: 10.1039/c2cp41116d

Due to the severe self-interaction errors associated with some density functional approximations, conventional density functionals often fail to dissociate the hemibonded structure of the water dimer radical cation (H₂O)₂⁺into the correct fragments: H₂O and H₂O⁺. Consequently, the binding energy of the hemibonded structure (H₂O)₂⁺is not well-defined. For a comprehensive comparison of different functionals for this system, we propose three criteria: (i) the binding energies, (ii) the relative energies between the conformers of the water dimer radical cation, and (iii) the dissociation curves predicted by different functionals. The long-range corrected (LC) double-hybrid functional,oB97X-2(LP) [J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2009, 131, 174105], is shown to perform reasonably well based on these three criteria. Reasons that LC hybrid functionals generally work better than conventional density functionals for hemibonded systems are also explained in this work.

I. Introduction

Water can be decomposed when it is exposed to high-energy flux. The products of water radiolysis may contain various radical species, e.g. hydrogen atoms (H), hydroxide radicals (OH), oxygen anions (O^!), and water cations (H₂O⁺), depending on the radiation infrastructure setup. For example the overall decomposition scheme activated byb particles has been outlined by Garrett et al. in 2005¹ where three main channels of decomposition were listed. The cationic channel leads to the formation of ionized water living for about several tens of femtoseconds and hydrated electrons, followed by the generation of hydronium (H₃O⁺) and OH radicals through proton transfer process.^2,3 The energized-neutral and anionic channels could result in the cleavage of the oxygen–hydrogen chemical bonds to produce hydrogen and oxygen derivatives, i.e. H, H^!, H₂, O, O^!, OH^! etc. Subsequent chemical reactions can progress further up to the desorption of stable gas molecules H₂and O₂ being driven by those reactive radical species.¹ The cationic channel is therefore particularly interesting due to its dominant products—OH radicals and solvated electrons.

The smallest system to understand the chemical dynamics of ionized water is the water dimer radical cation (H₂O)₂⁺, and it

has been approached by several experimental studies in the past. Angel and Stace reported the predominant H₃O⁺–OH central core from a collision-induced fragmentation experiment⁴ against the earlier theoretical assignment of a charge-resonance hydrazine structure.⁵Dong et al. observed a weak signal corres- ponding the formation of (H₂O)₂⁺ near the low-mass side of (H₂O)₂H⁺using a 26.5 eV soft X-ray laser.⁶Gardenier, Johnson, and McCoy reported the argon-tagged predissociation infrared spectra of (H₂O)₂⁺ and assigned its structural pattern as a charge-localized H3O⁺–OH complex.⁷Recently, Fujii’s group reported the infrared spectroscopic observations of larger (H₂O)_n⁺clusters, n = 3–11,⁸where the OH radical vibrational signal was clearly identified for n % 5 clusters, but the vibrational signature of the OH radical becomes inseparable due to the overlap with the H-bonded OH stretch in n4 6. As is evidenced in the earlier studies,^7,8theoretical investigations such as ab initio electronic structure theory and density functional theory (DFT) play an important role in understanding the infrared spectroscopic features of the ionized water clusters.

Because high-level ab initio calculations are computationally prohibited for larger ionized water clusters, e.g. fully solvated cationic moieties, a reliable DFT method is necessary.

In earlier theoretical reports, two minimum structures of the water dimer radical cation were identified: the proton transferred structure and the hemibonded structure, as shown in Fig. 1.^9–12The previous DFT calculations have shown that many exchange–

correlation (XC) functionals fail to predict reasonable results^9–11 giving rise to the presence of the hemibonding interaction. The hemibonding interaction, which could be theoretically located in (H₂O)_n⁺systems, is notorious for the serious self-interaction

aDepartment of Physics, National Taiwan University, Taipei 10617, Taiwan. E-mail: [email protected]

bDepartment of Chemistry, National Taiwan Normal University, Taipei 11677, Taiwan

cInstitute of Atomic and Molecular Sciences, Academia Sinica, Taipei 10617, Taiwan. E-mail: [email protected]

dCenter for Theoretical Sciences and Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan

PCCP Dynamic Article Links

www.rsc.org/pccp PAPER

Downloaded by National Taiwan University on 12 July 2012 Published on 07 June 2012 on http://pubs.rsc.org | doi:10.1039/C2CP41116D

(2)

errors (SIEs) associated with some density functional approximations. Both local density approximation (LDA) and generalized gradient approximations (GGAs) were reported to contain non-negligible amounts of SIEs for describing the hemibonded structure.^9,10It has been suggested to adopt hybrid functionals with larger fractions of the exact Hartree–Fock (HF) exchange for more accurate results for the hemibonded structure.^9,10 However, as the SIEs of functionals become larger at the dissociation limit, these suggested functionals can yield spurious barriers on their dissociation curves,¹⁰ which can lead to unphysical results in molecular dynamics simulations.

Clearly, more stringent criteria for choosing suitable functionals are needed. In this work, we propose three different criteria for a comprehensive comparison of different functionals for this system.

II. Computational methods

Calculations are performed on the optimized geometries of the two structures of the water dimer cation and the transition state between them, optimized with the ab initio MP2 theory¹³ and various XC functionals involving BLYP,^14,15PBE,¹⁶and M06L,¹⁷which are pure density functionals (i.e. the fraction of HF exchange aHF = 0.00), B97¹⁸ with aHF = 0.19, B3LYP^14,15,19 with aHF = 0.20, PBE0²⁰ with aHF = 0.25, M06²¹withaHF= 0.27, M05²²withaHF= 0.28, BH&HLYP²³

withaHF= 0.50, M06-2X²¹withaHF= 0.54, M05-2X²⁴with aHF = 0.56, M06HF²⁵ with aHF = 1.00, the oB97 series (oB97,²⁶oB97X,²⁶oB97X-D,²⁷andoB97X-2(LP)²⁸), which are long-range corrected (LC) hybrid functionals (i.e. the fraction of HF exchange depends on the interelectronic distance²⁹), and the double-hybrid functional B2PLYP³⁰with aHF= 0.53. The DFT and MP2 calculations are performed with the 6-311++G(3df,3pd) basis set, where the reference values of binding energy are obtained from ref. 12. For efficiency, the resolution-of-identity (RI) approximation³¹is used for calculations with the MP2 correlation (using sufficiently large auxiliary basis sets).

The CCSD(T) dissociation curves are calculated on the fixed monomer geometries (using the CCSD(T) optimized geometry of ref. 11), with the aug-cc-pVTZ basis set. Note that although the ZPE corrected energy of the proton transferred structure (or, referred to as the Ion structure in ref. 11) is inconsistent with ref. 12. However, adopting the geometries obtained from ref. 11 yields results that are consistent with ref. 12.

All of the calculations are performed with the development version of Q-Chem 3.2.³²As the basis set superposition error (BSSE) for the ionized water dimer has been shown to be insignificant (if the diffuse basis functions are adopted),^9,11we do not perform BSSE correction throughout this paper.

III. Results and discussion

The ZPE corrected binding energies and relative energies of the water dimer cation calculated by various XC functionals are shown in Tables 1 and 2, respectively. The calculated dissociation curves for the hemibonded structure are shown in Fig. 2. A summary of the results based on these three different criteria is shown in Table 4. The notation used for characterizing statistical errors is as follows: mean signed errors (MSEs), mean absolute errors (MAEs) and root-mean-square (RMS) errors.

Fig. 1 (a) The proton transferred structure, (b) the hemibonded structure, and (c) the transition state between the structures of (a) and (b).

Table 1 Binding energies (in kcal mol^!1) of the ionized water dimer

Method aHF

Proton transferred structure Transition state Hemibonded structure

MSE MAE RMS

E Error E Error E Error

BLYP 0.00 !45.62 2.10 — — !52.89 18.17 — — —

PBE 0.00 !47.37 3.85 — — !53.93 19.21 — — —

M06L 0.00 !45.00 1.48 !40.85 12.45 !48.24 13.52 9.15 9.15 10.65

B97 0.19 !45.33 1.81 !38.39 9.99 !45.88 11.16 7.65 7.65 8.71

B3LYP 0.20 -45.94 2.42 !38.43 10.02 !45.67 10.95 7.80 7.80 8.68

PBE0 0.25 !46.61 3.09 !36.75 8.35 !43.95 9.23 6.89 6.89 7.40

M06 0.27 !45.83 2.31 !36.13 7.73 !42.43 7.71 5.92 5.92 6.44

M05 0.28 !45.35 1.83 !35.46 7.06 !41.31 6.59 5.16 5.16 5.68

BH&HLYP 0.50 !45.97 2.45 !29.55 1.15 !35.11 0.39 1.33 1.33 1.58

B2PLYP 0.53 !45.04 1.52 !32.58 4.18 !40.96 6.24 3.98 3.98 4.42

M06-2X 0.54 !47.05 3.53 !32.21 3.81 !40.13 5.41 4.25 4.25 4.33

M05-2X 0.56 !46.76 3.24 !31.99 3.59 !39.35 4.63 3.82 3.82 3.87

oB97 0.00–1.00 !45.92 2.40 !33.63 5.23 !41.66 6.94 4.86 4.86 5.20

oB97X 0.16–1.00 !46.13 2.61 !34.34 5.94 !42.20 7.48 5.34 5.34 5.72

oB97X-D 0.22–1.00 !45.71 2.19 !35.33 6.93 !43.14 8.42 5.85 5.85 6.42

oB97X-2(LP) 0.68–1.00 !45.42 1.90 !29.17 0.77 !37.83 3.11 1.93 1.93 2.15

M06HF 1.00 !48.36 4.84 !28.23 !0.17 !35.75 1.03 1.90 2.01 2.86

MP2 1.00 !43.95 0.43 !25.16 !3.24 !30.03 !4.69 !2.50 2.79 3.30

CCSD(T)^a 1.00 !43.52^b 0.00 !28.39 0.00 !34.72 0.00 0.00 0.00 0.00

aThe CCSD(T) results, taken from ref. 12, are adopted as the reference.^bThe ZPE corrected binding energy of the proton transferred structure calculated by CCSD(T) in ref. 11 is inconsistent with ref. 12. However, adopting the geometry of ref. 12 will yield results that are consistent with ref. 12.

(3)

A. Criterion I: binding energies

Table 1 shows the binding energies of the two structures of the water dimer radical cation and the transition state between them. The reference data obtained from ref. 12 are based on the CCSD(T) calculations with the aug-cc-pVQZ basis set.

Since the errors of XC functionals for the hemibonded structure are much larger than those for the proton transferred structure, we focus our discussion on the hemibonded structure. From Table 1, the functionals with MAE less than 2.5 kcal mol^!1are oB97X-2(LP), M06HF, and BH&HLYP. Table 1 also confirms the trend that has been mentioned previously: functionals with larger fractions of HF exchange give more accurate results for the hemibonded structure. To give reasonable results for the hemibonded structure, the aHFof a global hybrid functional should be at least larger than 0.43, as observed in ref. 11 for the MPW1K functional. Although M06HF, containing a full HF exchange, gives a small error of the hemibonded structure, it yields a large error for the proton transferred structure (due to the incomplete cancelation of errors between the exact exchange and semilocal correlation), as shown in Table 1. Although functionals with aHF larger than 0.43 have been suggested, functionals with aHF Z 0.5 are not always reliable, which can be observed from the errors of the hemibonded structure

calculated by B2PLYP, M05-2X and M06-2X. But this trend still holds: the results of M05-2X and M06-2X are much better than those of M05 and M06. This means that although the energy of the hemibonded structure is sensitive to the aHF

values in XC functionals, it may also be affected by the associated density functional approximations (DFAs). Also note that some GGA functionals, such as BLYP and PBE, cannot predict that the transition state between the two structures of the water dimer radical cation.

The HF exchange included in theoB97 series is given by EoB97 series

HF exchange= E^LR!HFx + CxE^SR!HFx , (1) where

E^LR_x ^!HF ¼ !1 2

X

s

X

occu:

i;j

Z Z

c^#isðr¹Þc^#jsðr²Þ

&erfðor12Þ

r12 cjsðr¹Þcisðr²Þdr¹dr2;

ð2Þ

and

E_x^SR^!HF ¼ !1 2

X

s

X

occu:

i;j

Z Z

c^#_isðr¹Þc^#jsðr²Þ

&erfcðor¹²Þ

r₁₂ c_jsðr¹Þcisðr²Þdr¹dr2; ð3Þ

Here r12' |r¹²| = |r1! r²| (atomic units are used throughout this paper). The parametero defines the range of the splitting operators. The coefficients for the oB97 series are listed in Table 3. Since the fraction of HF exchange in theoB97 series depends on the interelectronic distance r₁₂, the trend mentioned previously is not as obvious as the global hybrid functionals.

But it is clear that the oB97X-2(LP), a LC double-hybrid functional, gives the most accurate results compare to the other functionals in theoB97 series.

As mentioned previously, functionals with large fractions of HF exchange may perform well for the hemibonded structure where the serious SIE takes place, they may perform unsatisfactorily for the other structures. Therefore, we also consider another criterion: the relative energies between the three structures, as proposed by Cheng et al.¹²

B. Criterion II: relative energies

In this criterion, the ground-state energy of the proton transferred structure is set to the zero point, i.e. both of the ground state energies of the hemibonded structure and that of the transition state are relative to the proton transferred structure, as shown in Table 2. The previously recommended functional, M06HF, performs poorly here.

In this criterion, it is obvious that functionals without the exact HF exchange, such as BLYP and PBE, overstabilize the hemibonded structure and wrongly predict the hemibonded structure to be more stable than the proton transferred one.

In addition to the previously recommended functionals based on Criterion I, the M05-2X and the M06-2X functionals also give accurate relative energies here. Although they give results that are not accurate enough for the binding energies, they yield good relative energies between those three structures of the water dimer radical cation.

Fig. 2 (a) Dissociation curves for the hemibonded structure calculated by various XC functionals. (b) Dissociation curves for the hemibonded structure calculated by MP2 and double-hybrid functionals.

(4)

In fact, functionals that are unable to give reasonable binding energies for the hemibonded structure may be traced back to the predicted dissociation curves of the hemibonded systems. Since many functionals fail to dissociate the hemibonded structure of the water dimer radical cation into the correct fragments, the definition of the binding energy is not well-defined.

Therefore, the entire dissociation curve from the hemibonded structure should be concerned.

C. Criterion III: dissociation behavior

Due to the severe SIEs associated with DFAs, systems with three-electron hemibonds, such as the hemibonded structure of the water dimer radical cation, are especially difficult for conventional density functionals. Many XC functionals cannot dissociate it into the correct fragments, H2O and H2O⁺(ionic state), i.e. they predict that the hemibonded structure should be dissociated into two fragments, each of which carries half of positive charge (covalent state). Fig. 2 shows the dissociation curves calculated by various XC functionals. Note that the discontinuous points in Fig. 2 near R = 2.5 angstrom for the oB97X-2(LP) and 3 angstrom for the M06HF functional are the respective broken-symmetry points.

The BH&HLYP functional, which has been suggested in the equilibrium ground-state energy calculation by the earlier reports^9,11 has a spurious barrier on its dissociation curve. Although we do not present the dissociation curve of the MPW1K functional, we expect it will suffer from the same problem as BH&HLYP.

The spurious barrier can be removed if the 100% exact exchange is adopted in a functional (e.g. M06HF), but its

shortcoming is described in the previous subsection and is thus not recommended. This shortcoming can be greatly reduced by the use of LC hybrid functionals, such as theoB97 series or the other LC hybrid functionals.³³The LC functionals retain the full HF exchange at the long range, while the good cancelation of errors between the semilocal exchange and correlation functionals are retained at the short range.²⁶

In the following, we will explain why LC hybrid functionals do not suffer from a spurious energy barrier as global hybrid functionals do. An estimate of the SIE of symmetric radical cations by global hybrid or pure density functionals has been derived as:^10,34

E^SIE( ð1 ! a^HFÞ 1 2! C

! "

Jþ 1 4R

# $

; ð4Þ

where CE 2^!1/3, (0.5! C) E !0.29, and J is the Coulomb self-interaction energy for the ionic state (in the case that the bond electron is localized at either of the two fragments). For three-electron-bonded radical cations, the bonding between the fragments is accomplished by the delocalizedb electron, which dominates the total SIE.¹⁰

We have derived an estimate of the SIE of symmetric radical cations by LC hybrid functionals, with the details arranged in the appendix, and the result is

E^SIE( ð1 ! C^xÞ 1

2! BðoÞC

! "

J^SRðoÞ þerfcðoRÞ 4R

# $

; ð5Þ where B(o) and J^SR(o) are constants with respect to R. But they depend ono, i.e. for LC hybrid functionals with different o values, their B(o) and J^SR(o) are different. The dependence of J^SR(o) on o is defined by

J^SRðoÞ ¼1 2

Z Z

r_bðr¹Þerfcðor¹²Þ

r12 r_bðr²Þdr¹dr2; ð6Þ and for smallo,

B(o) E 1 ! 0.254o(bohr^!1). (7) Table 2 Relative energies (in kcal mol^!1) between the three structures of the water dimer radical cation

Method aHF

Proton transferred structure Transition state Hemibonded structure

MSE MAE RMS

E Error E Error E Error

BLYP 0.00 0.00 0.00 — — !7.27 16.07 — — —

PBE 0.00 0.00 0.00 — — !6.56 15.36 — — —

M06L 0.00 0.00 0.00 4.14 10.86 !3.24 12.04 11.45 11.45 11.47

B97 0.19 0.00 0.00 6.93 8.07 !0.55 9.35 8.71 8.71 8.74

B3LYP 0.20 0.00 0.00 7.51 7.49 0.27 8.53 8.01 8.01 8.02

PBE0 0.25 0.00 0.00 9.86 5.14 2.66 6.14 5.64 5.64 5.66

M06 0.27 0.00 0.00 9.70 5.30 3.40 5.40 5.35 5.35 5.35

M05 0.28 0.00 0.00 9.89 5.11 4.04 4.76 4.94 4.94 4.94

BH&HLYP 0.50 0.00 0.00 16.41 !1.41 10.86 !2.06 !1.74 1.74 1.77

B2PLYP 0.53 0.00 0.00 12.46 2.54 4.08 4.72 3.63 3.63 3.79

M06-2X 0.54 0.00 0.00 14.85 0.15 6.92 1.88 1.02 1.02 1.33

M05-2X 0.56 0.00 0.00 14.77 0.23 7.41 1.39 0.81 0.81 1.00

oB97 0.00–1.00 0.00 0.00 12.29 2.71 4.27 4.53 3.62 3.62 3.73

oB97X 0.16–1.00 0.00 0.00 11.79 3.21 3.92 4.88 4.05 4.05 4.13

oB97X-D 0.22–1.00 0.00 0.00 10.38 4.62 2.57 6.23 5.42 5.42 5.48

oB97X-2(LP) 0.68–1.00 0.00 0.00 16.25 !1.25 7.59 1.21 !0.02 1.23 1.23

M06HF 1.00 0.00 0.00 20.12 !5.13 12.61 !3.81 !4.47 4.47 4.52

MP2 1.00 0.00 0.00 18.79 !3.79 13.92 !5.12 !4.45 4.45 4.50

CCSD(T)^a 1.00 0.00 0.00 15.14 0.00 8.80 0.00 0.00 0.00 0.00

aThe CCSD(T) results, taken from ref. 12, are adopted as the reference.

Table 3 The coefficients of the SR HF exchange ando for the oB97 series

oB97 oB97X oB97X-D oB97X-2(LP)

o(bohr^!1) 0.4 0.3 0.2 0.3

Cx 0.00 0.16 0.22 0.68

(5)

Note that estimate (4) is the special case with vanishingo of estimate (5). We apply estimates (5) and (7) to the simplest three-electron-bonded system, and the estimated He2+(which is also a three-electron hemibonded system) dissociation curves for pure DFT and LC hybrid functional are shown in Fig. 3(a). For simplicity, C_xis set to zero and the LDA orbital is used for evaluating the Coulomb self-interactions;o is set to 0.4 bohr^!1 for the LC hybrid functional. When the exact dissociation curve approaches zero, the SIE of the LC hybrid functional is already close to a constant, while the SIE of pure DFT is still decreasing. Therefore the dissociation curve by the LC hybrid functional does not display a spurious energy barrier as that of the pure DFT. Another effect of the long-range correction is the reduction of Coulomb self-interaction for the ionic state. In Fig. 3(a), (0.5! C)J E

!92 kcal mol^!1has been modified to [0.5! B(o)C]J^SR(o) E

!40 kcal mol^!1.

The formal feature of a spurious energy barrier is one more turning point (at which the derivative changes sign) on top of the barrier, in addition to the one in the potential well.

Since the ground state energy calculated by a functional is approximately the exact ground state energy plus the SIE produced by that functional,

4R²dE^SIE dR

%%

%% (2oR ffiffiffi

p exp½ ! ðoRÞp ²+ þ erfcðoRÞ: ð9Þ As shown in Fig. 3(b), for a sufficiento value (0.2 bohr^!1or more), this factor can change the decay nature of the derivative magnitude, from power law to exponential. Thus, the SIE derivative magnitude curve of typical LC hybrid functionals can avoid the second intersection with the derivative of the exact curve. Global hybrid functionals simply scale down the SIE derivative curve by a constant, so they cannot avoid the second intersection, unlessaHFapproaches unity. This explains why global hybrid functionals display a spurious energy barrier which LC hybrid functionals avoid.

Note that LC hybrid functionals which do not contain the SR HF exchange, such as oB97, can still be free from the spurious barrier, but will lose the possibility to predict symmetry breaking during the dissociation, i.e. the dissociation curve cannot converge to zero. This means that the SR HF exchange is also important. In fact, the covalent (symmetric) state and the ionic (symmetry-broken) state are nearly degenerate by CCSD(T) calculations.¹⁰But most of the XC functionals cannot predict that these two states are degenerate: due to the serious SIE, they usually overstabilize the covalent state. Therefore, a functional which can predict that the ionic state is more stable than the covalent state (i.e. the hemibonded structure will dissociate into H₂O and H₂O⁺) will give the correct dissociation limit. Very recently, a double-hybrid functional containing a very large fraction of HF exchange (E 79%),³⁶has been shown to be promising for reducing the SIEs in hemibonded systems.

The dissociation curves of the hemibonded structure calculated by double-hybrid functionals and MP2 are shown in Fig. 2(b).

Note that the PT2 calculation should be executed in the stable wave function. For example, although the dissociation curve of the covalent state seems more stable than that of the ionic state, we should choose the dissociation curve of the ionic state as an actual dissociation behavior calculated by MP2. Since the HF theory, which provides the reference orbitals for computing the MP2 correlation energy, predicts the ionic state to be more stable than the covalent one. Thus we choose the dissociation curve of the ionic state rather than that of the covalent state. Fig. 2(b) shows that theoB97X-2(LP) functional and the MP2 theory can predict the correct dissociation limits, while the B2PLYP functional cannot.

There is another way to define the dissociation behavior of the XC functionals: since many of the XC functionals cannot Fig. 3 The spurious energy barrier of the hemibonded systems pre-

dicted by DFT functionals can be illustrated qualitatively. (a) Com- parison of He2+dissociation curves by pure DFT and the LC hybrid functional using estimate (5). Zero level is set to E(He)+E(He⁺) for each method. (b) The first derivative of the exact He2+

dissociation curve and SIE derivative magnitude curves by pure DFT and the LC hybrid functional using estimate (9).

(6)

predict that the hemibonded structure of the ionized water dimer will dissociate into H₂O and H₂O⁺, we set the zeros of the dissociation curves to their respective dissociation limits, as shown in Fig. 4. In this definition, we focus on the potential curve experienced by the two fragments during the dissociation process of the hemibonded structure. Surprisingly, the dissociation curve of theoB97 functional is extremely close to that of the CCSD(T) theory. This means thatoB97 gives the best potential energy curve toward the dissociation process.

The previous suggested functionals,^9,11such as BH&HLYP, yield potential curves that are too shallow. Functionals which predict symmetry-breaking solutions during the dissociation process (e.g. M06HF and oB97X-2(LP)) are found to yield dissociation curves that are narrower than that of the CCSD(T) theory.

Finally, we discuss the dissociation of the proton transferred structure of the water dimer cation. In this structure, the SIEs associated with functionals for this structure are not as large as

those for the hemibonded structure, so all the dissociation curves are very similar, as shown in Fig. 5.

The results of the water dimer radical cation using three criteria we proposed are summarized in Table. 4. We find that the functional which performs well based on these three criteria is theoB97X-2(LP) functional, yielding the accurate binding energies, relative energies, and the correct dissociation limit. However, this functional yields a dissociation curve of the hemibonded structure that is a little narrower than that of CCSD(T). For applications sensitive to the shape of potential of the hemibonded structure, we suggest to use oB97:

although this functional neither gives a dissociation curve of the hemibonded structure that converges to zero nor yields accurate binding energies, this functional gives a dissociation curve which has nearly the same shape as the curve calculated by CCSD(T). Thus, we recommend this functional for researchers who like to perform the molecular dynamics of the water dimer cation.

Fig. 4 Dissociation curves for the hemibonded structure. The zeros of the dissociation curves are set to their respective dissociation limits.

Fig. 5 The dissociation curves for the proton transferred structure of ionized water dimer.

Table 4 Summary of results based on the three criteria

Method aHF

Criteria

Binding Energies Relative Energies

Correct dissociation limit

MSE MAE RMS MSE MAE RMS

BLYP 0.00 — — — — — — No

PBE 0.00 — — — — — — No

M06L 0.00 9.15 9.15 10.65 11.45 11.45 11.47 No

B97 0.19 7.65 7.65 8.71 8.71 8.71 8.74 No

B3LYP 0.20 7.80 7.80 8.68 8.01 8.01 8.02 No

PBE0 0.25 6.89 6.89 7.40 5.64 5.64 5.66 No

M06 0.27 5.92 5.92 6.44 5.35 5.35 5.35 No

M05 0.28 5.16 5.16 5.68 4.94 4.94 4.94 No

BH&HLYP 0.50 1.33 1.33 1.58 !1.74 1.74 1.77 No

B2PLYP 0.53 3.98 3.98 4.42 3.63 3.63 3.79 No

M06-2X 0.54 4.25 4.25 4.33 1.02 1.02 1.33 No

M05-2X 0.56 3.82 3.82 3.87 0.81 0.81 1.00 No

oB97 0.00–1.00 4.86 4.86 5.20 3.62 3.62 3.73 No

oB97X 0.16–1.00 5.34 5.34 5.72 4.05 4.05 4.13 No

oB97X-D 0.22–1.00 5.85 5.85 6.42 5.42 5.42 5.48 No

oB97X-2(LP) 0.68–1.00 1.93 1.93 2.15 !0.02 1.23 1.23 Yes

M06HF 1.00 1.90 2.01 2.86 !4.47 4.47 4.52 Yes

MP2 1.00 !2.50 2.79 3.30 !4.45 4.45 4.50 Yes

CCSD(T)^a 1.00 0.00 0.00 0.00 0.00 0.00 0.00 Yes

aThe CCSD(T) results, taken from ref. 12, are adopted as the reference.

(7)

IV. Conclusions

We have proposed three criteria to examine the performance of density functionals on the water dimer radical cation, and explained why LC hybrid functionals generally work better than conventional density functionals for hemibonded systems.

The previously recommended functional, BH&HLYP, cannot dissociate the hemibonded structure of the water dimer cation into the correct fragments: H₂O and H₂O⁺. Furthermore, the BH&HLYP dissociation curve displays an unphysical repulsive barrier, and is too shallow for molecular dynamics simulations.

Such a spurious barrier could be removed by functionals with very large fractions of HF exchange (e.g. M06HF andoB97X-2(LP)).

LC hybrid functionals, such as theoB97 series, are shown to be accurate for dissociation curves of the hemibonded structure (i.e. no spurious barriers), and thus are suitable for molecular dynamics simulations of larger-size systems. For research which is sensitive to the dissociation curve experienced by the fragments of the hemibonded structure, we recommend the use of theoB97 functional. For researchers who are not sure which criterion is the most important factor during the simulation of the water dimer radical cation, we recommend the use of the oB97X-2(LP) functional, which is overall good throughout the three criteria we proposed.

Appendix

We first reproduce estimate (4) for the SIE of symmetric radical cations by global hybrid or pure density functionals.^10,34Since most XC functionals predict that the covalent state is more stable than the ionic state, the SIE of the covalent state is of interest.

Neglecting the small SIE in the correlation energy, we have E^SIE_cov E J^SIcov+ E^SI_{x, cov}. (A1) Because the SIE is small in the ionic state,

E^SI_{x, ionic}E ! J^SIionic, (A2) it is favorable to express J^SIcov and E^SIx, cov in terms of the Coulomb self-interaction energy for the ionic state J^SI_ionic, or simply J in section III.C. This can be done by substituting the density of the delocalizedb electron (which causes the serious SIE of the hemibonding systems¹⁰) in the covalent state with the density of that electron in the ionic state. Note that there are two situations for the ionic state: one is the electric charge localized in fragment A and the other is the electric charge localized in fragment B. To a good approximation the density of the delocalizedb electron of the covalent state can be expressed as

r^bcovðrÞ (r^b_AðrÞ 2 þr^bBðrÞ

2 : ðA3Þ

The Coulomb self-interaction for the covalent state can be expressed as

J_cov^SI (1

2J_ionic^SI þ 1

4R; ðA4Þ

if one assumes that R is large compared to the spatial extent ofr

bAandr^bB. This expression can be applied to the self-interaction HF exchange energy for the covalent state,

E^{SI, HF}_{x, cov} =! J^SIcov. (A5)

The pure-DFT self-exchange energy for the covalent state is E^{SI, DFT}_{x, cov} E 2Ex(r^bA/2) = CE^SI_x,ionic. (A6) For LDA, C = 2^!1/3E 0.79. Combining

E^SI_{x, cov}=aHFE^{SI, HF}_{x, cov} + (1! aHF)E^{SI, DFT}_{x, cov} (A7) and estimate (A4), one can obtain estimate (4).

The SIE estimate for LC hybrid functionals can be derived in the same manner as the above one for global hybrid functionals. Substituting estimate (A3) into the self-interaction LR HF exchange yields

E_x;cov^SI;LR^!HF( ! 1

2J_ionic^SI;LRðoÞ þerfðoRÞ 4R

# $

; ðA8Þ

where we define J_ionic^SI;LRðoÞ ¼1

2 Z Z

r_bðr¹Þerfðor¹²Þ

r12 r_bðr²Þdr¹dr2: ðA9Þ Since the integration of J^{SI, LR}_ionic (o) is only over one fragment, it is independent of R. Likewise, the SR HF exchange of the covalent state is

E^SI;SR_x;cov^!HF( ! 1

2J_ionic^SI;SRðoÞ þerfcðoRÞ 4R

# $

; ðA10Þ

where we define J_ionic^SI;SRðoÞ ¼1

2 Z Z

rbðr¹Þerfcðor12Þ

r12 rbðr²Þdr¹dr2: ðA11Þ The SR-DFA self-exchange energy for the covalent state is

E^{SI, SR-DFA}x, cov E 2Ex(r^bA/2) = B(o)CE^SIx,ionic. (A12) For SR-LDA,^26,35

BðoÞ ¼ Z

r⁴⁼³_b ðrÞF 2¹⁼³o k_Fb

! "

dr 'Z

r⁴⁼³_b ðrÞF o k_Fb

! "

dr ðA13Þ k_Fb ' (6p²rb(r))^1/3 is the local Fermi wave vector, and the attenuation function is given by

FðlÞ ¼ 1 !2l 3 ½2 ffiffiffi

pperfðl^!1Þ þ , , ,+

( 1 !4lpffiffiffip

3 ; for smalll: ðA14Þ

For small o, and with the density approximated as one electron in a sphere with Bohr radius a₀,

BðoÞ ( 1 ! 3^!5=34 ffiffiffi

32 p ð ffiffiffi

32

p ! 1Þ ffiffiffi p p6

oa0( 1 ! 0:254oðbohr^!1Þ:

ðA15Þ Combining

E^SI_{x, cov}= E^{SI, LR-HF}_{x, cov} + C_xE^{SI, SR!HF}_{x, cov} + (1! Cx)E^{SI, SR!DFA}_{x, cov} (A16) and estimate (A4), one can obtain estimate (5).

Acknowledgements

We thank the support from National Science Council of Taiwan (Grant No. NSC99-2113-M-003-007-MY2, Downloaded by National Taiwan University on 12 July 2012 Published on 07 June 2012 on http://pubs.rsc.org | doi:10.1039/C2CP41116D

(8)

NSC98-2113-M-001-029-MY3 and NSC98-2112-M-002-023- MY3) and NCTS of Taiwan. We are grateful to the support from National Taiwan University (Grant No. 99R70304 and 10R80914-1) and Academia Sinica Research Program on NanoScience and Nano Technology. Computational resources are supported in part by the National Center for High Performance Computing.

References

1 B. C. Garrett, et al., Chem. Rev., 2005, 105, 355.

2 J. A. LaVerne and S. M. Pimblott, J. Phys. Chem. A, 2000, 104, 9820.

3 A. Furuhama, M. Dupuis and K. Hirao, Phys. Chem. Chem. Phys., 2008, 10, 2033–2042.

4 L. Angel and A. J. Stace, Chem. Phys. Lett., 2001, 345, 227.

5 R. N. Barnett and U. Landman, J. Phys. Chem. A, 1997, 101, 164.

6 F. Dong, S. Heinbuch, J. J. Rocca and E. R. Bernstein, J. Chem.

Phys., 2006, 124, 224319.

7 G. H. Gardenier, M. A. Johnson and A. B. McCoy, J. Phys. Chem.

A, 2009, 113, 4772.

8 K. Mizuse, J.-L. Kuo and A. Fujii, Chem. Sci., 2011, 2, 868, DOI:

10.1039/c0sc00604a.

9 M. Sodupe, J. Bertran, L. Rodrı´guez-Santiago and E. J. Baerends, J. Phys. Chem. A, 1999, 103, 166.

10 J. Gra¨fenstein, E. Kraka and D. Cremer, Phys. Chem. Chem.

Phys., 2004, 6, 1096.

11 H. M. Lee and K. S. Kim, J. Chem. Theory Comput., 2009, 5, 976.

12 Q. Cheng, F. A. Evangelista, A. C. Simmonett, Y. Yamaguchi and H. F. Schaefer, J. Phys. Chem. A, 2009, 113, 13779.

13 C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618.

14 A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098.

15 C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785.

16 J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.

17 Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2006, 125, 194101.

18 A. D. Becke, J. Chem. Phys., 1997, 107, 8554.

20 C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158.

21 Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215.

22 Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Phys., 2005, 123, 161103.

24 Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Theory Comput., 2006, 2, 364.

25 Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2006, 110, 13126.

26 J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.

27 J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615.

28 J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2009, 131, 174105.

29 J.-D. Chai and M. Head-Gordon, Chem. Phys. Lett., 2008, 467, 176.

30 S. Grimme, J. Chem. Phys., 2006, 124, 034108.

31 R. A. Kendall and H. A. Fru¨chtl, Theor. Chem. Acc., 1997, 97, 158.

32 Y. Shao, L. Fusti-Molnar, Y. Jung, J. Kussmann, C. Ochsenfeld, S. T. Brown, A. T. B. Gilbert, L. V. Slipchenko, S. V. Levchenko, D. P. O’Neill, R. A. Distasio Jr., R. C. Lochan, T. Wang, G. J. O.

Beran, N. A. Besley, J. M. Herbert, C. Y. Lin, T. Van Voorhis, S. H. Chien, A. Sodt, R. P. Steele, V. A. Rassolov, P. E. Maslen, P. P. Korambath, R. D. Adamson, B. Austin, J. Baker, E. F. C.

Byrd, H. Dachsel, R. J. Doerksen, A. Dreuw, B. D. Dunietz, A. D. Dutoi, T. R. Furlani, S. R. Gwaltney, A. Heyden, S. Hirata, C.-P. Hsu, G. Kedziora, R. Z. Khalliulin, P. Klunzinger, A. M. Lee, M. S. Lee, W. Liang, I. Lotan, N. Nair, B. Peters, E. I. Proynov, P. A. Pieniazek, Y. M. Rhee, J. Ritchie, E. Rosta, C. D. Sherrill, A. C. Simmonett, J. E. Subotnik, H. L. Woodcock III, W. Zhang, A. T. Bell, A. K. Chakraborty, D. M. Chipman, F. J. Keil, A. Warshel, W. J. Hehre, H. F. Schaefer III, J. Kong, A. I. Krylov, P. M. W. Gill and M. Head-Gordon, Phys. Chem.

Chem. Phys., 2006, 8, 3172.

33 Y.-S. Lin, C.-W. Tsai, G.-D. Li and J.-D. Chai, J. Chem. Phys., 2012, 136, 154109.

34 J. Gra¨fenstein, E. Kraka and D. Cremer, J. Chem. Phys., 2004, 120, 524.

35 P. M. W. Gill, R. D. Adamson and J. A. Pople, Mol. Phys., 1996, 88, 1005.

36 J.-D. Chai and S.-P. Mao, Chem. Phys. Lett., 2012, 538, 121.