Accurate Hydrogen Bond Energies within the Density Functional Tight

Article

Accurate Hydrogen Bond Energies within the Density Functional Tight Binding Method A. Domínguez,*,† T. A. Niehaus,‡ and T. Frauenheim† †

Bremen Center for Computational Materials Science, Universität Bremen, Am Fallturm 1, 28359 Bremen, Germany Department of Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany

‡

ABSTRACT: The density-functional-based tight-binding (DFTB) approach has been recently extended by incorporating one-center exchange-like terms in the expansion of the multicenter integrals. This goes beyond the Mulliken approximation and leads to a scheme which treats in a self-consistent way the fluctuations of the whole dual density matrix and not only its diagonal elements (Mulliken charges). To date, only the performance of this new formalism to reproduce excited-state properties has been assessed (Domínguez et al. J. Chem. Theory Comput., 2013, 9, 4901−4914). Here we study the effect of our corrections on the computation of hydrogen bond energies for water clusters and water-containing systems. The limitations of traditional DFTB to reproduce hydrogen bonds has been acknowledged often. We compare our results for a set of 22 small water clusters and watercontaining systems as well as for five water hexadecamers to those obtained with the DFTB3 method. Additionally, we combine our extension with a third-order energy expansion in the charge fluctuations. Our results show that the new formalisms significantly improve upon original DFTB. by the formalism.5 Furthermore, there have been some attempts to automatize the generation of the needed parameters, one recent effort being the development of a semiautomatic parametrization approach for the electronic part of DFTB.6 Other extensions include the development of QM/MM schemes based on DFTB7,8 and the implementation of empirical corrections to account for dispersion forces.9,10 Recently, the level of approximation in DFTB has taken a step forward with the inclusion of third-order terms in the energy expansion, which introduces a higher degree of selfconsistency and accuracy.11−14 This extension becomes particularly important for highly charged molecules, and combined with an empirical correction, it has been shown to improve the parameter transferability to reproduce hydrogen bonding energies and proton affinities.13 Whereas many efforts seem to be now devoted to what appears to be the next generation of DFTB, we believe that there is still plenty of room for improvements in second-order DFTB. Indeed, at this level of theory, two main approximations are applied, namely, the monopole approximation of the density fluctuations and the Mulliken approach for the evaluation of the multicenter integrals. Possible corrections to these approaches have, however, just started to be exploited. For the former approach, a refinement has been proposed, where dipole−monopole interactions are considered,15 whereas we have recently improved on the Mulliken approximation with the so-called onsite correction.16 In this approach, we propose a more accurate evaluation of multicenter integrals within DFTB, which does not imply additional computational effort. This scheme has been shown

1. INTRODUCTION Density functional theory (DFT) has become the method of choice to simulate a wide range of processes in quantum chemistry and computational material science due to its wellbalanced compromise between accuracy and efficiency. To address those problems still escaping from the scope of DFT due to their highly demanding computing resources, several approximate schemes have been developed. The density functional tight-binding (DFTB) approach1,2 has been shown to be a very useful tool, combining, in many cases, the accuracy typical of DFT with the efficiency representative of semiempirical methods. DFTB was originally based on the expansion of the DFT total energy up to the zeroth order in the charge density fluctuations around an input density, which is chosen as the sum of the densities of the isolated atomic constituents of the system in question. This original approach1 had limitations regarding the description of some molecular systems with an electronic density different from the mere superposition of neutral atomic contributions. With the self-consistent-charge (SCC) extension (expansion of the energy up to the secondorder), the method was significantly improved, thus accounting for charge transfer among atoms. Afterward, DFTB has been extended in multiple ways to gradually broaden its applicability. For example, to allow for the study of spin-polarized systems, the method was generalized to a spin-unrestricted formalism.3,4 Broadening the number of chemical elements (and the combination of elements) that can be investigated within the approach has also been a consistent focus. In this sense, relativistic effects have been incorporated in the parametrization process, thus widening the number of atomic species covered © 2015 American Chemical Society

Published: March 12, 2015 3535

DOI: 10.1021/acs.jpca.5b01732 J. Phys. Chem. A 2015, 119, 3535−3544

Article

The Journal of Physical Chemistry A

where f hxc = f H + f xc is the sum of the Coulomb and exchange− correlation (XC) kernels with

to markedly improve the description of singlet-to-triplet and σ → π electronic transitions.16 Despite the outstanding success of DFTB in many fields, serious shortcomings have been disclosed for some applications. The inaccurate description of hydrogen bonded complexes deserves special mention. In particular, DFTB has been shown to poorly reproduce bulk water and water clusters properties partly due to underestimation of the strength of hydrogen bonding interaction.17−20 This issue has been partly overcome with the introduction of a purely empirical shortrange correction to the γ functions involving hydrogen.12 A different empirical correction to hydrogen bond energies within DFTB has also been proposed, which corrects for the underestimated binding energies of neutral water clusters.21 In this contribution, we show that the onsite correction improves upon hydrogen bond energies of neutral, protonated, and hydroxide water clusters as well as methylimidazole−water complexes, and, at the same time, any empiricism in the formalism is eliminated. Apart from the considerable saving of computational time, semiempirical methods such as DFTB possess an important advantage for the description of hydrogen-bonded complexes over DFT approaches employing medium size basis sets, which suffer from nonphysical overbinding due to the basis set superposition error (BSSE). Althought there are several ways to correct for the BSSE, with the application of the counterpoise (CP) correction the most popular one,22 they are rarely employed in practical situations, especially when dealing with large complexes. DFTB has, therefore, a great potential to become an accurate, lightweigth alternative for the study of large H-bonding systems. The paper is organized as follows: In section 2, we describe the method in detail, whereas in section 3, the practical implementation of the method is depicted. Afterward, our results for a chosen set of systems are reported in section 4. Finally, the conclusions of our work are drawn in section 5.

fH (r, r′) =

δ 2E H[ρ] 1 = δρ(r)δρ(r′) |r − r′|

fxc [ρ](r, r′) =

δ 2Exc[ρ] δρ(r)δρ(r′)

(2)

and gxc is the derivative of the XC kernel with respect to the electron density: g xc[ρ](r, r′, r″) =

δ 3Exc[ρ] δρ(r)δρ(r′)δρ(r″)

(3)

2.1. Zeroth-Order Energy. DFTB0. The quantity , in eq 1 is the zeroth-order energy term and is sometimes called the Harris functional. This functional is defined as follows: ,[ρ] =

∑ niεi − EH[ρ] + Exc[ρ] − ∫ Vxc[ρ](r)ρ(r)d r + E NN i

(4)

The sum in eq 4 runs over molecular orbitals, with ni and εi being, respectively, the orbital occupation numbers and energies. Vxc= δExc/ δρ is the XC potential and ENN is the nucleus−nucleus repulsion energy. , [ρ]coincides with the total energy functional E[ρ]at the ground-state density and is stationary about this density.23 The problem of finding the stationary solution of , is equivalent to solving the DFT Kohn−Sham (KS) equations, whereas the error of , [ρ0]with respect to the actual ground-state energy is second order in the error of the guessed density ρ0. The advantage of the Harris scheme over DFT relies on the fact that , does not depend explicitly on the output density ρ, and the eigenvalues εi are the solution of a nonself-consistent Schrödinger equation where the effective potential, Vs, depends on ρ0 only: ⎤ ⎡ 1 0 Ĥ ψi(r) = ⎢− ∇2 + Vs[ρ0 ]⎥ψi(r) = εiψi(r) ⎦ ⎣ 2

2. FORMALISM The DFTB method has been described in detail elsewhere.1−3 Here we provide a brief review of the general formalism, where we focus on the second-order extension and its subsequent refinement. As mentioned in the previous section, DFTB is an approximation to DFT, where the Kohn−Sham total energy is expanded up to the zeroth, second, or third order in the charge density difference, Δρ, with respect to the input density, ρ0. The choice of truncation of the Volterra series defines different levels of approximation of the method. A zeroth-order expansion is the simplest approximation and the first employed formalism. It does not depend on the actual electron density but on an input or reference density, and it does not require a self-consistent treatment. Expansion of the Volterra series up to the second order introduces a dependence on the density fluctuations, thus leading to a self-consistent scheme with the iterative update of the electron density until convergence is reached. It is commonly referred to as self-consistent-charge DFTB (SCC-DFTB) or more recently as DFTB2. The most recent level of approximation entails the series truncation after the third-order term and when combined with the aforementioned empirical correction for hydrogen, it is known as DFTB3. The expression of the DFTB energy up to the third-order reads

(5)

The accuracy of this scheme depends on how well the reference density ρ0 represents the true density of the system. In DFTB, ρ0 is taken as the superposition of spherical atomiclike densities, ρ0 = ∑A ρA, subject to the constraint, ∫ ρA(r) dr = ZA. These atomiclike densities ρA are obtained via a self-consistent DFT atomic calculation where an extra parabolic potential is introduced for confinement purposes. The resulting densities are thus compressed with respect to the free-atom densities and are therefore a more suitable choice for the description of a wide variety of molecules and condensed systems. DFTB also employs a minimal basis set of compressed atomic orbitals (AO). These pseudoatomic wave functions are constructed by the linear combination of Slater-type orbitals (STO). All functions are centered at the corresponding atom. Orbitals belonging to the same atom are therefore orthogonal. By evaluating , [ρ] at ρ0 = ∑A ρA and expanding the XC part of the functional in a cluster series while truncating the expansion after the two-center term, it is possible to write , [ ρ0] in a Tight-Binding form: ,[ρ0 ] =

1 E[ρ] = ,[ρ0 ] + fhxc [ρ0 ](r, r′) Δρ(r)Δρ(r′) d rd r′ 2 1 + g xc[ρ0 ](r, r′, r″) Δρ(r)Δρ(r′) Δρ(r″) d rd r′d r″ + ... 6

∫∫

∑ niεi + Erep i

(6)

where the so-called repulsive energy, Erep, is given by the sum of short-range pairwise potentials. In DFTB, these pairwise potentials are not computed by directly evaluating at ρ0.

∫∫∫

(1) 3536


Article

The Journal of Physical Chemistry A Instead, they are fitted to ab initio results for a chosen reference system(s). In practice, the repulsive potentials are represented by a polynomial or spline for every pair of elements up to a limiting interatomic distance Rc,24 which is usually chosen between 1.5 and 2 times the equilibrium bond length. To evaluate the first term of eq 6, the variational principle is applied with the consequent expansion of the eigenfunctions ψi in the aforementioned set of AO, {ϕα}, thus leading to a secular equation |H 0 − εS| = 0

Thus, second-order corrections due to charge fluctuations lead to a KS Hamiltonian which includes up to first-order corrections. The Hamiltonian is now self-consistent and so the secular equations need to be solved iteratively. Next, two important approximations are applied to eqs 11 and 12, which simplify the evaluation of the required twoelectron integrals considerably. The evaluation of multicenter integrals is one of the main bottlenecks in molecular orbitals methods, and DFTB employs a well-known technique to simplify this process. Such a technique is the Mulliken approach, and it is one of the two main approximations in second-order DFTB. The Mulliken approximation consists of setting every differential overlap of AOs as ϕμϕν ≈ 1/2 Sμν(|ϕμ|2 + |ϕν|2), which automatically discards three- and four-center integrals from the scheme. The multicenter integrals appearing in eqs 11 and 12 can thus be written as a linear combination of one- and two-center integrals:

(7)

where = ⟨ϕα|Ĥ 0|ϕβ⟩ and Sαβ = ⟨ϕα|ϕβ⟩. Hamiltonian matrix elements involving three centers are neglected in DFTB. The rest of the Hamiltonian elements, along with the overlap matrix elements Sαβ, are determined from DFT calculations for every pair of chemical elements with respect to the interatomic distance and are tabulated. Thereby, no integral evaluation has to be performed during geometry optimizations or Molecular Dynamics simulations. 2.2. Second-Order Corrections. DFTB2. Now, we turn to the second order energy term in DFTB. Although second-order DFTB has been described in several occasions, we reintroduce here the necessary quantities and notation employed in the subsequent development of the onsite correction. We will consider a spin-unrestricted scheme.4 This generalization implies to write the second-order contribution to the energy as H0αβ

E(2) =

1 2

στ (μν|f hxc [ρ0 ]|κλ) ≈

E(2) =

(9)

∑∑ στ μνκλ

ΔPσμν

Pσμν

E(2) =

∑μcσμsϕμ,

(11)

−P0,σ μν

στ τ [ρ0 ]|κλ)ΔPκλ ∑ ∑ (μν|f hxc τ

κλ

(14)

(15)

(16)

1 2

∑∑

σ στ ΔqAl Γ Al , Bl ′ΔqBlτ ′

στ AB , ll ′

(17)

where qσAl is the Mulliken population for atom A, angular momentum l, and spin σ. Likewise, by applying the Mulliken (eq 13) and monopole (eq 16) approximations to eq 12, one arrives at

where = denotes the AO density matrix of the difference density Δ ρσ(r). The density matrix elements equal Pσμν= ∑s nsσcσμscσνs for orthogonal KS functions ψsσ, whereas P0,σ μν designates the density matrix elements of the reference system. By applying the variational principle to the energy functional up to the second order, E = , [ρ0] + E(2), we obtain a secular equation as in (7), but with a modified DFTB Hamiltonian, H = H0 + H1, where 1 Hμνσ =

μν

where the functions FA,l depend on the atom and angular momenta related to the AOs, but not on the magnetic number.16 Using the shorthand notation Γ, the second-order energy can be expressed as

∫ ∫ ϕμ(r)ϕν(r)g(r, r′)ϕκ(r′)ϕλ(r′) d rd r′

σ στ τ ΔPμν (μν|f hxc [ρ0 ]|κλ)ΔPκλ

στ

στ στ (μμ|f hxc [ρ0 ]|νν) ≈ (FAl|f hxc [ρ0 ]|FBl ′) = Γ στ Al , Bl ′

(10)

1 2

στ ̃τ [ρ0 ]|νν)ΔPνν ∑ ∑ ΔPμμ̃σ (μμ|f hxc

The second of the two main approximations in SCC-DFTB not only simplifies further the formalism but also ensures that the final total energy expression is rotational invariant. This property was lost when applying the Mulliken approximation. Rotational invariance (RI) is recovered by spherical averaging over AO products or, in other words, by taking the monopolar term in the multipole expansion of the AO product. We then have

As usual, we denote KS orbitals with Roman indices s, t, u, v and AO with Greek letters μ, ν, κ, λ. Furthermore, a general two-electron integral over a kernel g(r,r′) will be abbreviated as

E(2) =

1 2

1 P̃ σ = (Pσ ·S + S·Pσ ) 2

where the sum runs over the spin variables σ = ↑,↓, and Δρσ = ρσ − ρ0,σ. The output spin-dependent density ρσ= ∑s|ψsσ|2, the input density ρ0 = [ρ0, ↑, ρ0, ↓] is chosen spin-unpolarized (ρ0,↑(r)− ρ0,↓(r) = 0, ∀r) and the exchange−correlation component of f hxc now depends on the spin-dependent electron density:

By expanding the KS orbitals into the AO set, ψsσ = eq 8 can be written as

β ∈ {κ , λ}

σ where P̃ μν are the elements of the dual density matrix defined as16,25

(8)

(μν|g |κλ) =

στ (αα|f hxc [ρ0 ]|ββ)

(13)

στ

δ 2Exc[ρ] δρσ (r)δρτ (r′)

∑

By substituting eq 13 into 11, the second order contribution to the energy reads

στ [ρ0 ](r, r′) Δρτ (r′) d rd r′ ∑ ∫ ∫ Δρσ (r)f hxc

f xcστ [ρ](r, r′) =

1 SμνSκλ ∑ 4 α ∈ {μ , ν}

1 Hμνσ =

=

1 στ στ Sμν ∑ ∑ [(μμ|f hxc [ρ0 ]|γγ ) + (νν|f hxc [ρ0 ]|γγ )]ΔPγγ̃τ 2 τ γ 1 στ τ Sμν ∑ ∑ (Γ στ Al , Cl ″ + Γ Bl ′ , Cl ″)ΔqCl ″ 2 τ Cl ″

(18)

From eq 18, it should be noted that the atomic blocks of the Hamiltonian matrix are diagonal due to the orthogonality of the basis functions centered at a common atom. In practice, rather than working with spin-dependent charge fluctuations, it is preferred to evaluate the net charge and spin

(12) 3537


Article

The Journal of Physical Chemistry A populations during the self-consistent process. This is achieved after transformation from the set {ρ↑, ρ↓} to the total density ρ = ρ↑ + ρ↓ and magnetization m = ρ↑ − ρ↓. By this change of variables, the XC kernel can be split, and one arrives at16 Γ στ Al , Bl ′ = γAl , Bl ′ + δσ δτδABWAl , l ′

The derivative (∂γAl,Bl′/∂γAl,Al) is calculated analytically from the interpolation formula 22, whereas the Hubbard derivatives, 0 (∂γAl,Al/∂qAl)|qAl , can be either obtained using DFT selfconsistent-field (SCF) calculations or fitted to reproduce a desired property. In DFTB, third-order energy expansion is commonly accompanied by a purely empirical correction to the γ functions involving hydrogen12 (the modified functions will be denoted γh in the following) in order to improve the description of hydrogen bonding interaction. Within this correction, the short-range contribution, S, to the function γAl,Bl′ (eq 22) is damped by the factor

(19)

where δσ = 2δσ↑ − 1 and γAl , Bl ′ = (FAl|fhxc [ρ0 ]|FBl ′)

(20)

⎛ δ 2Exc[ρ , m] ⎜ WAl , l ′ = ⎜FAl δm(r)δm(r′) ⎜ ⎝

⎞ ⎟ FAl ′⎟ ⎟ ⎠

ρ0 ,0

⎡ ⎛γ + γBl ′ , Bl ′ ⎞ξ 2 ⎤ Al , Al exp⎢ −⎜ ⎟ RAB⎥ ⎢⎣ ⎝ ⎥⎦ ⎠ 2

(21)

The quantities γAl,Bl′ and WAl,l′ are referred to as the γ functions and spin coupling constants, respectively.26 The latter are treated as strictly on-site parameters, whereas γAl,Bl′ is calculated for every atom pair using an interpolation formula that depends on the distance RAB between the atoms A and B and the atomic Hubbard-like parameters γAl,Al and γBl′,Bl′: 1 − S(RAB , γAl , Al , γBl ′ , Bl ′) RAB

γAl , Bl ′ =

if at least one of atoms A and B is a hydrogen atom. The parameter ξ is fitted to reproduce a desired property, generally binding energies and/or proton affinities of hydrogen-bonded complexes. The combined employment of the third order extension and the damping factor 25 is referred to as the DFTB3 method, and it has been shown to further improve hydrogen bond energies compared to the solely application of the hydrogen correction.13 2.4. Onsite Correction. The Mulliken approach has been shown to be a good approximation to evaluate the differential overlap of AO residing on different atoms and has been successfully applied in many scenarios.28 However, if two different orbitals ϕμ and ϕν have a common center, their differential overlap strictly vanishes under this approximation as the overlap matrix element Sμν is zero for orthogonal basis functions. This implies that many three-, two-, and more importantly, one-center integrals are completely neglected in DFTB even though some of the latter may be fairly large. A visible consequence of this neglect was shown recently for the calculation of the absorption spectra of diatomic molecules.16 In section 4, we show that the Mulliken approximation is also responsible for the inaccurate description of hydrogen-bonded systems. To overcome the shortcomings of the Mulliken approach, we take the expansion of the multicenter integrals beyond the Mulliken term by including every onsite exchange-like term (μν|f στ hxc[ρ0]|μν), with μ,ν ∈ A and μ ≠ ν (see Appendix C in ref 16). By doing this, the multicenter integrals can be expressed as

(22)

S is a short-range function ensuring the correct convergence at RAB = 0 and the correct R−1 behavior at large interatomic distances. The quantities 20 and 21 can be easily evaluated by numerical integration. However, in DFTB, these integrals are, instead, approximated total energy derivatives with respect to the orbital occupation numbers using DFT calculations.2 Due to orbital relaxation, Hubbard parameters obtained in this way are roughly 10−20% smaller than the ones computed from a direct integral evaluation.16 It has been pointed out that use of a low-quality basis set leads to overestimated Hubbard parameters as the missing basis functions would otherwise screen its value.27 Thus, the approximation used for determining γAl,Al in DFTB partly compensates for the error introduced by the employed minimal basis. In part because of this error compensation, DFTB often returns better results than DFT itself when employing the same minimal basis set. 2.3. Third-Order Corrections. DFTB3. For the evaluation of the third-order energy, the same approximations as for the second order are applied, that is, the Mulliken and monopolar approaches. This leads to E(3) =

=

1 6 1 3

∑

ΔqAl ΔqBl ′

ABC , ll ′ l ″

∑

2 ΔqAl ΔqBl ′

∂γAl , Bl ′ ∂qCl ″

ΔqCl ″

στ στ στ (μν|f hxc [ρ0 ]|κλ) ≈ (μν|f hxc [ρ0 ]|κλ)mull + (μν|f hxc [ρ0 ]|κλ)ons

qCl0 ″

∂γAl , Bl ′

AB , ll ′

∂qAl

where 0 qAl

(23)

∂qAl

= 0 qAl

0 qAl

1 4 +

(26)

is given by eq 13 and

α≠μ στ Sαν(μα|f hxc [ρ0 ]|μα)[Sαλδμκ + Sακδμλ] α∈A β≠ν 1 στ Sβμ(βν|f hxc [ρ0 ]|βν)[Sβλδνκ + Sβκδνλ] 4 β∈B

∑

∑

1 στ SκνSμλ(μκ |f hxc [ρ0 ]|μκ )δAC(1 − δμκ ) 4 1 στ + SλνSμκ(μλ|f hxc [ρ0 ]|μλ)δAD(1 − δμλ) 4 1 στ + SκμSνλ(νκ |f hxc [ρ0 ]|νκ )δBC(1 − δνκ ) 4 1 στ + SλμSνκ(νλ|f hxc [ρ0 ]|νλ)δBD(1 − δνλ) 4 +

∂γAl , Bl ′ ∂γAl , Al ∂γAl , Al ∂qAl

mull (μν|f στ hxc[ρ0]|κλ)

στ (μν|f hxc [ρ0 ]|κλ)ons =

where the derivatives of the γ functions are evaluated at the reference atomic charges.13 Please, note that we returned here to a spin-restricted analysis to simplify the appearing equations. A generalization to a spin-unrestricted scheme can be easily formulated. The derivatives in eq 23 can be expressed in terms of the derivatives of the atomic gamma (Hubbard parameters) ∂γAl , Bl ′

(25)

(24) 3538

(27)


Article

The Journal of Physical Chemistry A where it is assumed that μ ∈A, ν ∈B, κ ∈C and λ ∈D. It should ons be noted that the integral (μν|f στ is nonzero if at hxc[ρ0]|κλ) least two of the four orbitals ϕμ, ϕν, ϕκ and ϕλ share a common center, which means that the onsite correction strictly excludes four-center integrals. Moreover, as it should be expected, integrals of the type (μμ|f στ hxc[ρ0]|νν), which are given exactly within the Mulliken approach, are not affected by our correction. After substituting the onsite-corrected integrals 26 into 11 and doing some lengthy but straightforward algebra, it is possible to obtain a refined expression for the DFTB secondorder energy 1 2

E(2) =

In the standard DFTB approach, this condition is met as both terms on the left-hand side of eq 30 are approximated to the parameter Γστ Ap,Ap, while the integral on the right-hand side is neglected. Within the new approach, the right-hand side of eq 30 is no longer zero, and hence, the monopolar approach does not ensure RI. One way to retain this property consists of evaluating every onsite integral exactly, so that eq 30 holds automatically. However, as mentioned above, the use of the screened Γ is convenient as it countervails the modest quality of the basis set employed in DFTB. To preserve RI while still compensating for the employment of a minimal basis set, we proceed as follows: we obtain in practice only the exchange-like integrals directly from the wave functions, whereas the coulombστ like integrals, (pp|f στ hxc[ρ0]|pp) and (pp|f hxc[ρ0]|p′p′), are evaluated using the identity 30 and the known (approximate) Γ:

στ ̃τ ∑ ∑ ΔPμμ̃σ (μμ|f hxc [ρ0 ]|νν)ΔPνν στ

μν μ≠ν

+

∑∑ ∑ A

στ

̃σ (μν|f στ [ρ ]|μν)ΔPμν ̃τ ΔPμν hxc 0

μν ∈ A

4 στ (pp′|f hxc [ρ0 ]|pp′) 3 2 στ στ (pp|f hxc [ρ0 ]|p′p′) = Γ στ (pp′|f hxc [ρ0 ]|pp′) Ap , Ap − 3

στ (pp|f hxc [ρ0 ]|pp) = Γ στ Ap , Ap +

(28)

Equation 28 differs from the conventional formula 14 in the additional second term, which depends on the fluctuations of the off-diagonal elements of the dual density matrix, P̃σμν, connecting orbitals ϕμ and ϕν residing at the same atom. Similarly, by substituting 26 into 12, it is possible to obtain the onsite-corrected Hamiltonian Hnew= Hold+ Hons, where Hold is given by eq 18, and Hons is given by

Apart from the improved accuracy, the use of the traditional Γ keeps the modifications of the original method as small as possible, while RI is still exactly preserved. Additionally, in this way, the simplicity of the scheme is not destroyed, and the number of parameters to add is small. For onsite integrals involving d-orbitals, a similar identity holds: στ στ στ (dd|f hxc [ρ0 ]|dd) − (dd|f hxc [ρ0 ]|d′d′) = 2(dd′|f hxc [ρ0 ]|dd′)

α≠μ ons Hμνσ =

στ ̃τ ∑ ∑ Sαν(μα|f hxc [ρ0 ]|μα)ΔPμα τ

+

(32)

α∈A β≠ν στ ̃τ ∑ ∑ Sβμ(βν|f hxc [ρ0 ]|βν)ΔPβν τ

β∈B

(31)

στ but in this case (dd|f στ hxc[ρ0]|d′d′) and (dd′ |f hxc[ρ0]|dd′) stand for the mean values of the coulomb and exchange-like onsite integrals, respectively, formed by every possible combination (ten in total) of orbitals d and d′, with d ≠ d′. In a similar fashion as for p-orbitals, only the integral (dd′ |f στ hxc[ρ0]|dd′) needs to be computed, whereas the rest of onsite integrals are obtained by using eq 32 and 16:

(29)

Unlike one might deceptively think, the onsite contribution to the Hamiltonian has nonzero off-site elements, that is, our correction affects also matrix elements, Hμνσ, with orbitals ϕμ and ϕν placed at different atoms. The major effect relies, however, on the atomic blocks of the Hamiltonian matrix. Within the new formalism, the Hamiltonian is self-consistent in terms of the whole dual density matrix and not only of its diagonal elements (Mulliken population). Within the onsitecorrected DFTB, convergence is reached when a converged dual density matrix is obtained. As this represents a somewhat heavier convergence criterion, our scheme is in principle expected to moderately increase the number of SCC cycles with respect to the original scheme.

8 στ (dd′|f hxc [ρ0 ]|dd′) 5 2 στ [ρ0 ]|dd′) − (dd′|f hxc 5

στ (dd|f hxc [ρ0 ]|dd) = Γ στ Ad , Ad + στ (dd|f hxc [ρ0 ]|d′d′) = Γ στ Ad , Ad

(33)

As a last requirement for a basis set containing up to d functions, it is necessary to set also every (pp|f στ hxc[ρ0]|dd) integral to its averaged value, which is given by the exact ΓAp,Ad. However, RI fulfillment does not depend on the value of this quantity and so the usual screened parameter is also employed in this case. Off-site integrals are exempt from the RI issue and can be approximated as usual using eq 16. Splitting the first term of eq 28 into onsite and off-site contributions while applying eq 16 for the off-site component and 31 and 33 for the onsite part, one obtains

3. ONSITE CORRECTION IN PRACTICE The usual procedure for the evaluation of the first term of eq 28 within DFTB is to apply the monopole approximation (eq 16). However, with the inclusion of the exchange-like integrals in the second term, the spherical averaging over AO products will in general not lead to an expression that is invariant under a rotation of the molecular frame. A similar problem with the rotational invariance (RI) of the system was pointed out for the INDO model in a paper by Figeys and co-workers.29 In this work, the authors state that, for example, for two orbitals with angular momentum p centered at the same atom A, the following identity must hold to preserve RI:

E(2) =

1 2

∑∑

σ τ ΔqAl Γ στ Al , Bl ′ ΔqBl ′

στ AB , ll ′ μ≠ν

+

∑∑

̃σ Λστ ̃τ ΔPμν Al , Al ′ ΔPμν

∑

στ A , ll ′ μ ∈ A , l ; ν ∈ A , l ′

1 3

∑∑ ∑

1 + 5

∑∑ ∑

+

στ στ στ (pp|f hxc [ρ0 ]|pp) − (pp|f hxc [ρ0 ]|p′p′) = 2(pp′|f hxc [ρ0 ]|pp′), ∀ p , p′ ∈ A

(30)

στ

στ

A

A

̃σ Λστ ̃τ (3δμν − 1)ΔPμμ Ap , Ap ′ ΔPνν

μν ∈ A , p

̃σ Λστ ̃τ (5δμν − 1)ΔPμμ Ad , Ad ′ ΔPνν

μν ∈ A , d

(34) 3539


Article


Table 1. Comparison of Binding Energies As Obtained with the DFTB Method at Different Levels of Approximation and G3B3a system

DFTB2

DFTBh2

OC-DFTB2

DFTB3

OC-DFTB3

DFTBh3

G3B3

(H2O)2 (H2O)3 (H2O)4 (H2O)5 H+(H2O)2 H+(H2O)3 H+(H2O)4 H+(H2O)5 OH−(H2O) OH−(H2O)2 OH−(H2O)3 OH−(H2O)4 NH3 (H2O) NH4+ (H2O) (H2O)6 [book] (H2O)6 [cage] (H2O)6 [prism] (H2O)6 [ring] methylimidazole(-H+)(H2O) methylimidazole (H2O) [donor] methylimidazole (H2O) [acceptor] methylimidazole H+ (H2O) MUD MSD MAX

1.6 5.5 9.7 13.3 4.5 10.4 13.9 18.3 −5.1 −2.6 0.3 6.1 3.2 0.6 16.7 17.2 17.6 16.5 4.1 2.4 3.5 3.3 8.0 7.3 18.3

0.0 −0.6 0.6 1.4 −2.0 −0.1 1.1 1.8 −12.8 −17.0 −17.5 −18.2 2.1 −3.4 1.2 0.3 0.0 1.8 2.0 1.4 2.6 1.2 4.0 −2.5 18.2

0.2 −0.4 −1.4 −1.4 −4.2 −2.4 −1.0 0.0 −10.4 −8.7 −10.5 −9.2 1.8 −4.6 −1.3 −0.8 0.0 −1.6 −0.6 1.2 1.4 −0.4 2.9 −2.5 10.5

1.5 5.4 9.4 12.5 5.9 11.6 15.0 19.7 1.5 5.3 9.0 14.2 3.1 1.4 16.5 17.6 18.0 15.3 3.2 2.6 2.8 3.9 8.9 8.9 19.7

0.3 0.0 −1.1 −1.5 −2.3 0.1 1.9 3.7 −2.3 −0.4 −0.2 0.8 1.7 −3.4 −0.7 0.5 1.1 −2.1 −1.9 1.7 0.5 0.5 1.3 −0.1 3.7

0.0 −0.3 0.8 1.3 0.9 3.7 5.0 6.2 −5.9 −8.4 −7.2 −7.8 2.1 −1.3 1.7 1.5 1.3 1.5 1.2 1.9 1.9 2.3 2.9 0.1 8.4

−4.9 −15.1 −27.4 −36.3 −33.9 −57.3 −77.2 −91.9 −27.4 −48.6 −66.7 −86.3 −6.6 −20.4 −45.8 −46.6 −47.2 −44.7 −15.9 −6.2 −8.2 −16.0

For all DFTB variants, the reported values are the deviations from the G3B3 results. (H2O)n, H+(H2O)n, and OH−(H2O)n denote a neutral, protonated, and deprotonated (hydroxide) water n-mer, respectively. The structures of the four water hexamers are depicted in Figure 1. Methylimidazole(H+)(H2O) and methylimidazole(-H+)(H2O) denote, respectively, a protonated and deprotonated methylimidazole complexed with water, whereas methylimidazole (H2O) [donor] ([acceptor]) stand for the neutral methylimidazole complexed with water as a hydrogen-bond donor (acceptor). All energies are expressed in kcal/mol. a

στ where the shorthand Λστ Al,Al′ = (ll′ |f hxc[ρ0]|ll′) is used. The first term in eq 34 accounts for the noncorrected energy (eq 17). The second term corrects partly for the employed Mulliken approximation applied in the first term. It accounts for the interaction of the fluctuations of off-diagonal elements of the dual density matrix. Finally, third and fourth terms correct for the RI loss due to the concurrent application of the monopolar approximation in the first term and the onsite correction (second term). It is worth stressing that RI is satisfied using eq 34 regardless of the employed Γ and Λ parameters. The latter (which will be referred to as the onsite parameters) are in our approach neither freely adjustable nor fitted to experiments. Instead, they are computed numerically at the PBE level (see ref 16 for details on their evaluation). Due to symmetry, only 10 onsite parameters per element are needed for angular momenta up to l = 2:

4. RESULTS As mentioned above, hydrogen bond energies were substantially improved with the introduction of a short-range correction to the γ functions, and further enhanced with DFTB3, which combines this empirical correction with a thirdorder energy expansion. In this section, we assess the performance of the onsite-corrected DFTB (OC-DFTB) method for the description of hydrogen bonds. To this end, we calculated the binding energies of 22 small systems containing water, including neutral, protonated, and deprotonated water clusters. This benchmark set has been employed elsewhere for the validation of different variants of DFTB3.12,13 The results are shown in Table 1. Different DFTB approaches have been employed. For ease of reference, we distinguish in the following each DFTB variant by a subscript indicating the expansion order, whereas the superscript h will indicate that the empirical correction for hydrogen is applied. Thus, the DFTB3 method, for example, is denoted DFTBh3. The investigated DFTB approaches include the original second-order scheme using the standard and modified γ functions (DFTB2 and DFTBh2, respectively) as well as using third order (DFTB3) and onsite (OC-DFTB2) corrections. It is also interesting to assess the performance of the combination of both the onsite and third order corrections (OC-DFTB3) and so these results are included for comparison. We finally include DFTBh3 results using the Hubbard derivatives obtained from PBE calculations and the parameter ξ fitted to reproduce the binding energy of the neutral water dimer. G3B331 results are employed here as

↑↓ ↑↑ ↑↓ ↑↑ ↑↓ ↑↑ Λ↑↑ As , Ap, ΛAs , Ap, ΛAp , Ap ′, ΛAp , Ap ′, ΛAs , Ad , ΛAs , Ad , ΛAp , Ad , ↑↑ ↑↓ Λ↑↓ Ap , Ad , ΛAd , Ad ′, ΛAd , Ad ′

(35)

These integrals are calculated for every atom type, stored in a file and read during the calculation. In Appendix A, calculated onsite parameters are reported for some elements. Hydrogen has vanishing Λ parameters because the employed basis set contains only s-functions residing on this element. Therefore, no onsite correction applies for hydrogen atoms. The present DFTB scheme has been implemented in a development version of the DFTB+ code.30 3540


Article


Figure 1. Four water hexamers investigated in Table 1.

the description of the water hexamers and methylimidazole− water complexes, with errors under 2 kcal/mol per hydrogen bond. When combining the onsite correction with the thirdorder extension, the results are even better. The MUD is in this case reduced to 1.3 kcal/mol with a maximum individual deviation of 3.7 kcal/mol. The MSD is also small (−0.1 kcal/mol), which indicates that the obtained values are scattered around the reference data. Among all tested DFTB variants in Table 1, the OC-DFTB3 scheme is the only one for which the errors concerning the hydroxide water clusters are really consistent with those related to the neutral and positive charged systems. This indicates that a higher parameter transferability is achieved with such a scheme. More importantly, both OC-DFTB2 and OC-DFTB3 methods are totally free of empirical or semiempirical parameters, which is a desired feature. In our calculations, the computational effort required by every DFTB variant was similar, with the total elapsed (wall clock) time employed for the single-point energy calculation of the 22 systems ranging from 5.8 to 5.9 s, using an Intel Xeon processor E3−1225 at 3.10 GHz. Although the use of the modified γh improves considerably the performance of DFTB for the description of hydrogen bonds with respect to the standard method, the obtained results are yet not satisfactory. This is especially observed for larger systems, and it has been pointed out by Truhlar et al.18 In a recent contribution, they tested the accuracy of 85 DFT and semiempirical methods to reproduce the binding energies of five water hexadecamers (16-mers) (see Figure 2) obtained by Yoo and co-workers,32 using coupled cluster theory with quasiperturbative triplet excitations [CCSD(T)].33,34 The geometries of the water nanoparticles and monomer were optimized at the MP2/aug-cc-pVTZ level. As shown in that work, DFTB2 (Note: In ref 18, standard self-consistent-charge DFTB (DFTB2) is termed SCC-DFTB, whereas DFTBh2 is referred to as SCC-DFTB-γh.) highly underestimates the strength of the hydrogen bonding interaction with a mean binding energy error of 68 kcal/mol. This issue is importantly alleviated with DFTBh2, although absolute values of the binding energy are still substantially underestimated. Furthermore, Truhlar and co-workers assessed the accuracy of the methods in terms of relative energies of the water 16-mers. For this property, the tested DFTB variants perform well, with MUD of 1.2 and 1.7 kcal/mol when employing the standard and modified γ, respectively. In order to assess quantitatively the ability of the investigated approaches to reproduce both relative energies and absolute binding energies, Truhlar et al. defined a characteristic error (CE) as follows:

reference. Except for the onsite corrected schemes, OC-DFTB2 and OC-DFTB3, the binding energies were taken from ref 13. All geometries, optimized at the B3LYP/6-31G(d) level of theory, were also extracted from the mentioned work. Binding energies are defined as the energy difference between the complex and the isolated molecules, so that a negative value designates a bonded system. For all DFTB variants, the deviation with respect to G3B3 results is reported, where the DFTB energies are taken as the minuend, that is, a positive value indicates an underestimation of the absolute value of the binding energy (and hence, an underestimation of the strength of the hydrogen bonding interaction). At the bottom of the table, the mean unsigned and signed deviation or error (MUD and MSD, respectively) as well as the maximum deviation from G3B3 results (MAX) are given for every tested method. As shown in Table 1, except for the deprotonated systems, the standard DFTB2 method underestimates the strength of the hydrogen bonding interaction with an error of about 2−4 kcal/mol per hydrogen bond. This underestimation is especially noticeable for the protonated systems, for which errors are large. In contrast, the underestimation for the hydroxide water tetramer and pentamer is rather small, whereas for the hydroxide dimer and trimer, the strength of the hydrogen bonding interaction is overestimated. Use of the empirical γh functions substantially improves upon binding energies for the neutral and protonated clusters but clearly worsens the description of the deprotonated ones with a remarkable overestimation of the strength of the hydrogen bonding interaction. The overall performance of this correction is however adequate with a MUD reduction of a half. Third-order correction in combination with the standard γ functions does not improve the results overall. However, using both the modified γh and the third-order energy expansion results in a reduced MUD of 2.9 kcal/mol. As can be seen in the table, DFTBh3 retains the performance of DFTBh2 for the neutral and protonated water clusters while systematically improving the description of the deprotonated systems. The overall enhancement of DFTBh3 over the traditional scheme should be, however, mainly attributed to the empirical correction for hydrogen and only in part to the combined application of this modification and the thirdorder extension. A global improvement is also obtained when the onsite correction is applied, resulting in a MUD (2.9 kcal/mol) similar to that obtained with the DFTBh3 method. In contrast to the standard DFTB2 case, our formalism tends to overestimate the strength of the hydrogen bonding interaction. This overestimation is especially important for the hydroxide water clusters. The errors for these systems are smaller than those obtained with the DFTBh2 variant but considerably larger than for the standard approach. Particularly outstanding is, in contrast,

⎛ MUDX ⎞ ⎛ MUDX ⎞ BE RE ⎟ ⎟ + 0.5⎜⎜ CE = 0.5⎜⎜ med ⎟ med ⎟ ⎝ MUDBE ⎠ ⎝ MUDRE ⎠ 3541

(36)


Article


provide the MUDBE, MUDRE and Truhlar’s CE for every DFTB variant. Our results reveal that the inclusion of onsite corrections cures the failure of standard DFTB for the description of these systems. OC-DFTB2 and OC-DFTB3 exhibit MUDBE of only 10.5 and 14.0 kcal/mol, respectively. Third-order corrections alone do not improve the performance of standard DFTB. However, DFTBh3 does return binding energies in better agreement with the CCSD(T) findings. OC-DFTB2 also keeps a low MUDRE, and as a result, it exhibits a Truhlar’s CE of 0.56. With this result, our method would occupy the 16th place in Truhlar’s table, above DFT approaches such as PBE (CE = 0.77) and hybrid schemes such as HSE06 (CE = 0.70) and PBE0 (CE = 0.79). DFTBh3 would also obtain a distinguished position in the table, with a CE of 0.61 when using third-order parameters computed by DFT. Combination of third-order and onsite corrections apparently is not a good formula here; such formalism scores a CE close to 1.0 and would be considered as an “average“ method by Truhlar and co-workers. It is important to note that the aforementioned DFT approaches tested by Truhlar do not have any post-SCF correction. However, the employed basis set (jun-cc-pVTZ) is of sufficient quality, and the BSSE is expected to be small. The smaller CE of onsite-corrected DFTB compared to PBE is due to the slightly better description of the relative energies within DFTB and not because of better absolute values of the binding energies. Indeed, according to the latter indicator, PBE outperforms OC-DFTB2 with a MUD of 5.9 kcal/mol, which is considerably smaller than for our formalism. PBE also returns correct relative energies for the 16-mer boat and antiboat structures (binding energies of boat-a: −166.67 kcal/mol, boat-b: −166.46 kcal/mol, antiboat: −166.10 kcal/mol), but it predicts too low absolute values of the binding energies for the 4444-a and 4444-b clusters (−162,67 and −162.24 kcal/mol, respectively) in comparison to the other structures, although still in better agreement with CCSD(T) values than OC-DFTB2. This results in a MUDRE of 2.5 kcal/mol, compared to the 1.2 kcal/mol value within OC-DFTB2. The reasonably good performance for the relative energies is not a particular feature of the onsite correction, but rather of DFTB in almost any of its variants, as it can be seen in Table 2. Only the concurrent application of the onsite and third-order corrections worsens this result, which approaches the PBE value. Regarding the irregular performance of the OC-DFTB3 variant compared to OC-DFTB2, we would like to mention that the onsite correction constitutes an important step toward

Figure 2. Five water hexadecamers investigated by Yoo et al.32

MUDXBE and MUDXRE are, respectively, the mean unsigned errors in the five binding energies and ten relative energies as computed with method X. MUDmed BE = 19.9 kcal/mol and X MUDmed RE = 2.0 kcal/mol are the median of all of the MUDBE X and MUDRE values, respectively. In ref 18, Truhlar’s group sorted the 85 tested methods by their CE. The authors consider as satisfactory those approaches with a CE well below 1.0, whereas those with CE values close to or greater than 1.0 are disapproved. DFTB2 belongs to the latter group of approaches with a CE = 2.02. If the empirical γh function is employed, the value decreases to 1.31, but the use of the method is still discouraged. We would like to use this indicator to assess our method and compare our results to those obtained by Truhlar. First, we show in Table 2 the binding energies for the five water 16-mers as obtained with OC-DFTB2, DFTB3, OC-DFTB3 and DFTBh3. The same MP2 geometries for the water 16-mer and monomer employed in ref 18 are used throughout. We also include the results for DFTB2 and DFTBh2 as well as the reference CCSD(T) energies from ref 18. At the bottom of the table, we

Table 2. Comparison of Binding Energies of Five Water 16-Mers As Obtained with the DFTB Method at Different Levels of Approximation and CCSD(T)a structure

DFTB2

DFTBh2

OC-DFTB2

DFTB3

OC-DFTB3

DFTBh3

CCSD(T)

4444-a 4444-b antiboat boat-a boat-b MUDBE MUDRE CE

67.03 66.90 68.14 69.04 68.98 68.0 1.2 2.02

33.99 33.45 35.80 36.48 36.46 35.2 1.7 1.31

12.01 10.97 9.50 9.91 10.10 10.5 1.2 0.56

68.26 68.20 67.82 68.73 68.67 68.3 0.5 1.83

16.58 15.67 12.22 12.61 12.77 14.0 2.4 0.95

13.63 12.99 14.22 14.86 14.93 14.1 1.0 0.61

−171.06 −170.52 −170.55 −170.80 −170.64

For all DFTB variants, the reported values are the deviations from the CCSD(T) results. The five water configurations are depicted in Figure 2. At the bottom of the table, we provide the mean unsigned errors in the five binding energies (MUDBE) as well as the mean unsigned errors in the 10 relative energies of every combination pair (MUDRE). Truhlar’s characteristic error (CE) is also given for every DFTB variant. All energies are expressed in kcal/mol. a

3542


Article


This scheme substantially improves upon hydrogen bond energies for a set of 22 water-containing systems. The overall accuracy is comparable to that obtained with the DFTB3 method, which employs empirical parameters fitted to reproduce this property. If, in addition, a third-order extension is combined with our correction, hydrogen bond energies are further improved. Within such a scheme, the transferability to treat different charge states in hydrogen-bonded systems clearly surpasses that of DFTB3, whereas no empirical parameters are necessary. The description of larger neutral water clusters is also improved with respect to standard DFTB results when employing onsite corrections. In this case, however, the combination of thirdorder extensions with our refinement does not outperform the onsite-corrected DFTB. We additionally showed that DFTB3 also describe these systems accurately. The onsite correction has again the advantage of requiring no fitted parameters. We would finally like to mention that all onsite-corrected results reported here were obtained with the existing set of parameters (mio set) generated for the conventional secondorder DFTB. In the future, the DFTB parameters will be adapted for the onsite-corrected scheme and, hence, an improved accuracy is expected. Especially, the generation of adapted repulsive potentials will allow us to assess the performance of the onsite correction for the description of structural properties of molecules and solids and run onsitecorrected molecular dynamic simulations.

the refinement of second-order DFTB, but there is certainly still room for improvements. Therefore, inclusion of higherorder corrections might not necessarily always lead to improvements as they are expected to have a small effect compared to the yet missing second-order corrections. One example of this kind of corrections is the further refinement of the multicenter integrals by including terms of the form (μν|f hxc|κλ), where μ,ν ∈A and κ,λ ∈B with A ≠ B. These terms are neglected within the Mulliken approach (and the onsite correction) if at least one of the conditions μ ≠ ν and κ ≠ λ is met. This future prospect is expected to lead to a further enhanced characterization of hydrogenbonded complexes. For example, apart from the interaction of hydrogen s orbitals with oxygen s or p orbitals considered within the Mulliken approach, (sOsO|f hxc|sHsH) and (pOpO|f hxc|sHsH), the aforementioned two-center-like correction would also be able to model interactions of the type (sOpO|f hxc|sHsH). Another important point to analyze is the fact that the onsite correction improves upon hydrogen bonding despite that hydrogen atoms are not directly affected by this correction. More importantly, both this and the empirical correction γh seem to have a similar effect on the results. In the onsite correction, the self-interaction of the off-diagonal elements of the dual density matrix is modeled. Although no correction applies for the dual density matrix corresponding to H, the dual density matrix elements P̃ pOp′O= 1/2 ∑k(PpOkSkp′O + SpOkPkp′O) of an oxygen atom interacting with the H atom carry information on this interaction through the overlap matrix elements SsHp′O and SpOsH, where sH donotes the s orbital of the H atom. Thus, the onsite correction for O does indirectly strengthen the hydrogen bond interaction. The γh correction, in contrast, directly alters the γO−H function, which mediates the interaction between the oxygen and sH electrons, thereby making it more repulsive in the short-range. Both corrections then lead to the increase of the water molecule polarization compared to the DFTB2 results. The dipole moment of a single water molecule, calculated according to Mulliken population analysis, is 1.688 D within the original approach, which is too small compared to the experimental value 1.855 D35 or the QCISD/aug-cc-PVQZ value 1.925 D.36 The γh and onsite corrections return molecular dipole moments of 1.830 and 1.927 D, respectively, which are more in line with the experimental and ab initio results. The increase of the water polarization in both cases leads to the desired strengthening of the hydrogen bond interaction within DFTB. As a final check, we computed the binding energies for the two water heptadecamers studied in ref 18 using OC-DFTB2, OC-DFTB3, and DFTBh3. These methods also perform fairly well in this case, with MUDBE of 10.1, 13.7, and 15.6 kcal/mol and relative energy errors of 0.26, 1.39, and 0.27 kcal/mol, respectively. Our results thus suggest that both the onsitecorrected DFTB and DFTB3h methods may be reliably employed for the study of some properties of neutral bulk water at a little computational cost.

A. ONSITE PARAMETER FOR SOME ELEMENTS In Appendix Table 3, we report calculated onsite parameters for some elements. Table 3. Onsite Parameters for Some Chemical Elements Calculated at the PBE Level of Theory element H C N O S Ti Au

■

5. CONCLUSIONS A formulation of the onsite-corrected DFTB method for the calculation of ground-state properties has been presented in detail. This approach deals with the refinement of one of the most important approximations within DFTB, namely, the Mulliken approximation. By restoring the disregarded onecenter exchange-like integrals into the formalism, we arrive at a scheme that treats in a self-consistent way the fluctuations of the whole dual density matrix and not only its diagonal elements (charges), as in the original approach.

sp(↑↑)

pp′(↑↑)

sd(↑↑)

pd(↑↑)

dd′(↑↑)

sp(↑↓)

pp′(↑↓)

sd(↑↓)

pd(↑↓)

dd′(↑↓)

0.00000 0.00000 0.04973 0.10512 0.06816 0.12770 0.08672 0.14969 0.07501 0.11653 0.02659 0.06881 0.03752 0.06928

0.00000 0.00000 −0.01203 0.02643 −0.00879 0.03246 −0.00523 0.03834 0.00310 0.03058 −0.01297 0.01640 −0.00505 0.01677

0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00398 0.03915 −0.00587 0.01239 0.00073 0.01339

0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.01100 0.04979 −0.00523 0.01144 −0.00002 0.01228

0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 −0.01792 0.01582 −0.00750 0.02604 0.00531 0.02519

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support from the Deutsche Forschungsgemeinschaft (GRK 1570 and SPP 1243) and the DAAD is gratefully acknowledged. Dr. B. Aradi is acknowledged for fruitful discussions 3543


Article


(18) Leverentz, H. R.; Qi, H. W.; Truhlar, D. G. Assessing the Accuracy of Density Functional and Semiempirical Wave Function Methods for Water Nanoparticles: Comparing Binding and Relative Energies of (H2O)16 and (H2O)17 to CCSD(T) Results. J. Chem. Theory Comput. 2013, 9, 995−1006. (19) Choi, T. H.; Liang, R.; Maupin, C. M.; Voth, G. A. Application of the SCC-DFTB Method to Hydroxide Water Clusters and Aqueous Hydroxide Solutions. J. Phys. Chem. B 2013, 117, 5165−5179. (20) Liang, R.; Swanson, J. M. J.; Voth, G. A. Benchmark Study of the SCC-DFTB Approach for a Biomolecular Proton Channel. J. Chem. Theory Comput. 2014, 10, 451−462. (21) Miro, P.; Cramer, C. J. Water Clusters to Nanodrops: A TightBinding Density Functional Study. Phys. Chem. Chem. Phys. 2013, 15, 1837−1843. (22) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553. (23) Foulkes, W. M. C.; Haydock, R. Tight-Binding Models and Density-Functional Theory. Phys. Rev. B 1989, 39, 12520−12536. (24) Elstner, M. Weiterentwicklung quantenmechanischer Rechenverfahren für organische Moleküle und Polymere. Ph.D. thesis, Universität Gesamthochschule Paderborn, 1998. (25) Han, M. J.; Ozaki, T.; Yu, J. O(N) LDA+U Electronic Structure Calculation Method Based on the Nonorthogonal Pseudoatomic Orbital Basis. Phys. Rev. B 2006, 73, 45110. (26) Frauenheim, T.; Seifert, G.; Elstner, M.; Niehaus, T.; Kohler, C.; Amkreutz, M.; Sternberg, M.; Hajnal, Z.; Di Carlo, A.; Suhai, S. Atomistic Simulations of Complex Materials: Ground-State and Excited-State Properties. J. Phys.: Condens. Matter 2002, 14, 3015− 3047. (27) Cox, B. N.; Coulthard, M. A.; Lloyd, P. A Calculation of the Coulomb Correlation Energy, U, for Transition Metals in Hubbard’s Model. J. Phys. F: Met. Phys. 1974, 4, 807−820. (28) Ayed, O.; Bernard, E.; Silvi, B. On the Mulliken Approximation of Multicenter Integrals. J. Mol. Struct.: THEOCHEM 1986, 135, 159− 168. (29) Figeys, H. P.; Geerlings, P.; Van Alsenoy, C. Rotational Invariance of INDO Theories Including d-Orbitals into the Basis Set. Int. J. Quantum Chem. 1977, 11, 705−713. (30) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method. J. Phys. Chem. A 2007, 111, 5678−5684. (31) Baboul, A. G.; Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-3 Theory Using Density Functional Geometries and ZeroPoint Energies. J. Chem. Phys. 1999, 110, 7650−7657. (32) Yoo, S.; Aprà, E.; Zeng, X. C.; Xantheas, S. S. High-Level Ab Initio Electronic Structure Calculations of Water Clusters (H2O)16 and (H2O)17: A New Global Minimum for (H2O)16. Phys. Chem. Lett. 2010, 1, 3122−3127. (33) Purvis, G. D.; Bartlett, R. J. A Full Coupled Cluster Singles and Doubles Model: The Inclusion of Disconnected Triples. J. Chem. Phys. 1982, 76, 1910−1918. (34) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A Fifth-Order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 1989, 157, 479−483. (35) Clough, S. A.; Beers, Y.; Klein, G. P.; Rothman, L. S. Dipole Moment of Water from Stark Measurements of H2O, HDO, and D2O. J. Chem. Phys. 1973, 59, 2254−2259. (36) Martin, F.; Zipse, H. Charge Distribution in the Water Molecule - A Comparison of Methods. J. Comput. Chem. 2005, 26, 97−105.

and his assistance during the implementation of the OC-DFTB formalism.

■

REFERENCES

(1) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. Construction of Tight-Binding-like Potentials on the Basis of DensityFunctional Theory: Application to Carbon. Phys. Rev. B 1995, 51, 12947−12957. (2) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260−7268. (3) Frauenheim, T.; Seifert, G.; Elstner, M.; Hajnal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R. A Self-Consistent Charge Density-Functional based Tight-Binding Method for Preditive Materials Simulations in Physics, Chemistry and Biology. Phys. Status Solidi B 2000, 217, 41−62. (4) Köhler, C.; Seifert, G.; Gerstmann, U.; Elstner, M.; Overhof, H.; Frauenheim, T. Approximate Density Functional Calculations of Spin Densities in Large Molecular Systems and Complex Solids. Phys. Chem. Chem. Phys. 2001, 3, 5109−5114. (5) Witek, H. A.; Kohler, C.; Frauenheim, T.; Morokuma, K.; Elstner, M. Relativistic Parametrization of the Self-Consistent-Charge DensityFunctional Tight-Binding Method. 1. Atomic Wave Functions and Energies. J. Phys. Chem. A 2007, 111, 5712−5719. (6) Wahiduzzaman, M.; Oliveira, A. F.; Philipsen, P.; Zhechkov, L.; van Lenthe, E.; Witek, H. A.; Heine, T. DFTB Parameters for the Periodic Table: Part 1, Electronic Structure. J. Chem. Theory Comput. 2013, 9, 4006−4017. (7) Seabra, G. M.; Walker, R. C.; Elstner, M.; Case, D. A.; Roitberg, A. E. Implementation of the SCC-DFTB Method for Hybrid QM/ MM Simulations within the Amber Molecular Dynamics Package. J. Phys. Chem. A 2007, 111, 5655−5664. (8) Woodcock, H. L.; Hodoscek, M.; Brooks, B. R. Exploring SCCDFTB Paths for Mapping QM/MM Reaction Mechanisms. J. Phys. Chem. A 2007, 111, 5720−5728. (9) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen Bonding and Stacking Interactions of Nucleic Acid Base Pairs: A Density-Functional-Theory based Treatment. J. Chem. Phys. 2001, 114, 5149. (10) Zhechkov, L.; Heine, T.; Patchkovskii, S.; Seifert, G.; Duarte, H. A. An Efficient a Posteriori Treatment for Dispersion Interaction in Density-Functional-Based Tight Binding. J. Chem. Theory Comput. 2005, 1, 841−847. (11) Elstner, M. The SCC-DFTB Method and its Application to Biological Systems. Theor. Chem. Acc. 2006, 116, 316−325. (12) Yang, Y.; Yu, H.; York, D.; Cui, Q.; Elstner, M. Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method: Third-Order Expansion of the Density Functional Theory Total Energy and Introduction of a Modified Effective Coulomb Interaction. J. Phys. Chem. A 2007, 111, 10861−10873. (13) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the SelfConsistent-Charge Density-Functional Tight-Binding Method (SCCDFTB). J. Chem. Theory Comput. 2011, 7, 931−948. (14) Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. J. Chem. Theory Comput. 2013, 9, 338−354. (15) Bodrog, Z.; Aradi, B. Possible Improvements to the SelfConsistent-Charges Density-Functional Tight-Binding Method within the Second Order. Phys. Status Solidi B 2012, 249, 259−269. (16) Domínguez, A.; Aradi, B.; Frauenheim, T.; Lutsker, V.; Niehaus, T. A. Extensions of the Time-Dependent Density Functional Based Tight-Binding Approach. J. Chem. Theory Comput. 2013, 9, 4901− 4914. (17) Maupin, C. M.; Aradi, B.; Voth, G. A. The Self-Consistent Charge Density Functional Tight Binding Method Applied to Liquid Water and the Hydrated Excess Proton: Benchmark Simulations. J. Phys. Chem. B 2010, 114, 6922−6931. 3544


Accurate Hydrogen Bond Energies within the Density Functional Tight

Recommend Documents