Analytic cubic and quartic force fields using density-functional theory

Analytic cubic and quartic force fields using density-functional theory

Magnus Ringholm,¹Dan Jonsson,^1,2Radovan Bast,³Bin Gao,¹Andreas J. Thorvaldsen,¹ Ulf Ekström,⁴Trygve Helgaker,⁴and Kenneth Ruud¹

1Centre for Theoretical and Computational Chemistry (CTCC), Department of Chemistry, University of Tromsø—The Arctic University of Norway, 9037 Tromsø, Norway

2High Performance Computing Group, University of Tromsø—The Arctic University of Norway, 9037 Tromsø, Norway

3Theoretical Chemistry and Biology, School of Biotechnology, Royal Institute of Technology,

AlbaNova University Center, S-10691 Stockholm, Sweden and PDC Center for High Performance Computing, Royal Institute of Technology, S-10044 Stockholm, Sweden

4Center for Theoretical and Computational Chemistry (CTCC), Department of Chemistry, University of Oslo, P.O. Box 1033, Blindern, 0315 Oslo, Norway

(Received 16 October 2013; accepted 19 December 2013; published online 15 January 2014) We present the first analytic implementation of cubic and quartic force constants at the level of Kohn–Sham density-functional theory. The implementation is based on an open-ended formalism for the evaluation of energy derivatives in an atomic-orbital basis. The implementation relies on the availability of open-ended codes for evaluation of one- and two-electron integrals differentiated with respect to nuclear displacements as well as automatic differentiation of the exchange–correlation ker- nels. We use generalized second-order vibrational perturbation theory to calculate the fundamental frequencies of methane, ethane, benzene, and aniline, comparing B3LYP, BLYP, and Hartree–Fock results. The Hartree–Fock anharmonic corrections agree well with the B3LYP corrections when cal- culated at the B3LYP geometry and from B3LYP normal coordinates, suggesting that the inclusion of electron correlation is not essential for the reliable calculation of cubic and quartic force constants.

© 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.[http://dx.doi.org/10.1063/1.4861003]

I. INTRODUCTION

Vibrational spectroscopy provides a rich and diverse source of information about molecular structure and func- tionality. For this reason, methods for calculating molecular vibrational spectra were developed already in the early days of quantum-chemical computations, to help elucidate molec- ular structure and to provide insight into experimental obser- vations. Already in 1958, Bratoz recognized the benefits of calculating the forces acting on the nuclei of a molecule in an analytic manner.¹ However, the breakthrough in terms of computing molecular gradients and force fields came with the efficient implementation of molecular gradients by Pulay in 1969.²

Molecular gradients are the first-order derivatives of the molecular energy with respect to nuclear displacements and can be determined from the unperturbed electron density and from differentiated one- and two-electron integrals. The com- plexity increases significantly when going to the molecu- lar force fields or molecular Hessians (corresponding to the second-order derivatives of the molecular energy with respect to nuclear displacements) as also the perturbed electron den- sity matrix is needed in this case. Thomson and Swanstrøm presented the first implementation of molecular Hessians in 1973.³In the early 1980s, implementations of molecular Hes- sians for unrestricted and restricted open-shell Hartree–Fock (HF) wave functions were published by Yamaguchi, Schaefer, and co-workers.^4,5 In the following years, second-order geo- metrical derivatives were derived and implemented for a wide

range of correlated wave functions.^6–15More recently, imple- mentations of molecular Hessians have been presented at the level of density-functional theory (DFT).^16–18 For a detailed historical account, see recent reviews of analytic derivative techniques for molecular properties in general and molecular force fields in particular.^19–21

In parallel with the development of force-field calcu- lations for correlated wave functions, Schaefer, Handy, and co-workers extended the evaluation of geometrical deriva- tives of the HF energy to third^22,23 and fourth²⁴ orders.

Despite the importance of cubic and quartic force fields for determining, for instance, anharmonic corrections to vi- brational frequencies, Fermi resonances,²⁵ rotation–vibration constants,²⁶ vibrationally averaged geometries,^27–29 and l- doubling constants,^25,30,31 these implementations have seen little use in the literature. Part of the reason for this may be that anharmonic corrections are important only in high- accuracy studies,^32–34 where the HF approximation may not be sufficient. Also, the implementation of efficient schemes for obtaining anharmonic force fields by numerical differenti- ation of forces and force constants at highly correlated lev- els of theory has proven feasible, even for relatively large molecules.^35,36

To some extent, the low cost of Kohn–Sham DFT has changed this picture. In a series of studies, Barone and co- workers have shown that high-quality harmonic force fields in combination with DFT anharmonic corrections provide reliable estimates of anharmonic force fields^37,38 and of

0021-9606/2014/140(3)/034103/11 140, 034103-1 © Author(s) 2014

anharmonic corrections to intensities,³⁹ thus demonstrating that anharmonic effects can now be studied straightforwardly for large and complex molecular systems. In their studies, ge- ometrical derivatives beyond second order were determined by finite differences.^32,40,41 The high cost of such finite- difference schemes and their numerical instability provide a strong motivation for developing analytic methods for third- and fourth-order geometrical derivatives at the DFT level.

Also, fifth- and sixth-order geometrical derivatives are needed in fourth-order vibrational perturbation theory (VPT4); re- cently, we have calculated analytically quintic and sextic force fields at the HF level of theory.⁴²

In this work, we describe the analytic calculation of cu- bic and quartic force constants at the level of Kohn–Sham DFT, using generalized-gradient-approximation (GGA) and hybrid functionals. This work builds on several developments in our groups in recent years. In particular, we use the frame- work of an atomic-orbital (AO) based, open-ended quasi- energy response-theory formalism described by Thorvaldsen et al.,⁴³ which for force fields reduces to regular energy- derivative theory in the AO basis,⁴⁴extended to fourth-order derivatives. The geometrical derivatives of the one-electron integrals arising from the geometry dependence of the AOs are evaluated using the one-electron integral framework of Gao, Thorvaldsen, and Ruud.⁴⁵ The evaluation of geometri- cal derivatives of the two-electron repulsion integrals follows the approach of Reine, Tellgren, and Helgaker,⁴⁶ expanding solid-harmonic Gaussians directly in Hermite Gaussians. We furthermore extend automatic differentiation of exchange–

correlation kernels⁴⁷ to include corrections arising from the dependence of the AOs on the nuclear positions. Finally, we demonstrate the usefulness of the code by evaluating the cubic and quartic force constants and anharmonic force-field correc- tions for selected molecules.

The bulk of this paper is organized as follows. In Sec. II, we give a brief account of the AO-based energy- derivative framework used by us and give the expressions for the cubic and quartic force constants. A brief description of the evaluation of the exchange–correlation contribution is also given. Section III contains computational details of the cal- culations, while Sec. IV presents and discusses the results.

Finally, in Sec.V, we give some concluding remarks and per- spectives for the analytic calculation of higher-order proper- ties that involve geometrical distortions.

II. THEORY

We here present the theory behind our AO-based imple- mentation of DFT cubic and quartic force fields, the AO- formulation ensuring that the approach is suitable for linear- scaling methodology.⁴⁸ The approach builds on the general AO-based framework for time- and perturbation-dependent basis sets by Thorvaldsen et al.,⁴³ here applied to time- independent perturbations. We note that, even though ex- plicit equations are given for the evaluation of the cubic and quartic force fields, our implementation uses a recursive scheme, for which explicit expressions for energy derivatives are not needed.⁴⁹Compared with the earlier implementations by Schaefer, Handy, and co-workers,^22–24our implementation

is extended to DFT and formulated entirely in the AO basis, although at present a molecular-orbital-based response solver is used to determine the perturbed density matrices.

In Sec.II A, we review the basic theory of the analytic calculation of geometrical derivatives at the DFT level, pro- viding explicit expressions for the cubic and quartic force constants; next, in Sec. II B, we describe the evaluation of exchange–correlation contributions to the cubic and quartic force constants, combining the perturbation dependence of the overlap distributions with the use of automatic differentia- tion to evaluate the higher-order exchange–correlation kernel derivatives.⁴⁷

A. AO-based energy derivative theory

We follow the notation of the AO-based response the- ory for self-consistent-field (SCF) methods with time- and perturbation-dependent basis sets by Thorvaldsen et al.,⁴³ specialized to static perturbations. In Kohn–Sham DFT, the energy can in the AO basis be written as

E(D)=^TrhD+1

2G^γ (D)D+E_xc[ρ(D)]+h_nuc. (1) In this expression, “Tr” indicates that the trace of the ma- trix products on the right-hand side of the equation is taken, E_xc[ρ(D)] is the exchange–correlation energy, which is a functional of the generalized density vectorρ(see Sec.II B), hnucis the nuclear repulsion energy, andDis the AO density matrix. The density matrix fulfills the idempotency relation

0=DSD−D. (2) We have in Eqs.(1)and(2)also introduced the one-electron integral matrix, h, the two-electron integral matrix con- structed fromDwithγ fractional exchange,G^γ(D), and the overlap matrix,S, whose elements are given by

hμν= χμ| − 1

2∇²−

K

ZK

|RK−r||χν, (3)

G^γ_μν(M)=

αβ

M_βα(gμναβ−γ g_μβαν), (4)

Sμν = χμ|χν. (5) The χμ are spherical-harmonic Gaussian AOs and the sum- mation in Eq.(3)is over atomic nuclei atRK and with charge ZK. The two-electron integrals are defined in the conventional manner as

gμνρσ =

dx₁dx₂χ_μ^∗(x₁)χν(x₁)r₁₂⁻¹χ_ρ^∗(x₂)χσ(x₂), (6) with integration over all spin and spatial coordinates.

The optimized Kohn–Sham density fulfills the SCF condition

FDS−SDF=0, (7) where the Kohn–Sham matrixFis defined as

F= ∂E

∂D^T =h+G^γ(D)+F_xc, (8)

with

Fxc=

dr ∂E_xc(r)

∂ρ(r)

∂ρ(r)

∂D^T . (9) Following Refs. 43, 49, and 50, we take as our start- ing point for the generation of higher-order derivatives the energy-gradient LagrangianL^adefined as

L^a=L^a(D,ζa,λa)=^Tr ∂E(D)

∂εa

−S^aW−λaY−ζaZ, (10) where the superscriptaindicates differentiation with respect to an applied perturbation of strengthε_a. In Eq.(10)we have introduced the energy-weighted density matrix

W=DFD, (11)

as well as matrices that represent the constraints on the un- perturbed reference state—in particular, the idempotency and SCF-state matrices, respectively,

Z=DSD−D, (12)

Y=FDS−SDF. (13) For an optimized SCF state, Eqs.(2) and(7) can be written compactly asZ=0andY=0. Additionally, we have intro- duced the Lagrange multipliers λa andζa, respectively, for these constraints. In Ref.43, it is shown that Eq.(10)is vari- ational inDif the zeroth-order multipliers are defined as

λa =D^aSD−DSD^a, (14)

ζa=F^aDS+SDF^a−FDS^a−S^aDF−F^a. (15) The subscriptaon the multipliers does not indicate differen- tiation, merely the relation toL^a.

Since we require Eq. (10)to be variational in the den- sityDand the multipliersλa andζa, we can take advantage of the 2n +1 rule for the density and the 2n +2 rule for the multipliers^50,51 when differentiating the energy gradient Lagrangian.

1. Molecular gradient

Let us first consider the first-order geometrical deriva- tives of the molecular energy. In this case, Eq.(10)simplifies to Pulay’s expression from 1969²

dE

dεa =L^a=^Tr ∂E

∂εa −S^aW, (16) where no derivatives of the Lagrange multipliers are required because of the 2n+2 rule.⁵⁰We note that the molecular gra- dient can be determined from a knowledge of the unperturbed density alone, in accordance with the 2n+1 rule.

2. Molecular Hessian

Differentiating Eq.(10)with respect tob, keeping only terms that fulfill the 2n+1 and 2n+2 rules, we obtain for

the molecular Hessian the expression d²E

dεadεb = dL^a

dεb =L^ab=^Tr ∂²E

∂εa∂εb + ∂²E

∂εa∂D^TD^b

−S^abW−S^aW^b. (17) As for the gradient, no zeroth-order multipliers are required.

However, the first-order perturbed density matrix is needed to calculate the molecular Hessian; we return to the evaluation of the perturbed densities in Sec.II A 5. Henceforth, we adopt the subscript notation of Ref.43, writing the Hessian as

d²E

dεadεb =L^ab_0,1 =^TrE_0,1^ab −(SW)^ab_1W, (18) where superscripts denote total derivatives. The subscripts (k, n) specify the maximum order of the perturbed density matrix D: to orderkfor collections of perturbations involving pertur- bation a, and to ordern for collections of perturbations not involving perturbation a. The notationnW for the SWterm specifies in a similar manner the maximum order of differen- tiation ofWin this term, as dictated by the value ofn.

The Hessian expression in Eq.(17)is not explicitly sym- metric in aandb (the numerical values, of course, are). As shown by Sellers, an explicitly symmetric formula can be ad- vantageous from a numerical point of view.⁵²

3. Cubic force constants

For the calculation of cubic force constants, the third- order energy derivative is needed. Proceeding as above, we differentiate the gradient Lagrangian twice, keeping only terms that fulfill the 2n+1 and 2n+2 rules

L^abc_1,1 =^Tr(E^a)^bc₁ −S^abcW−S^acW^b−S^abW^c−S^aW^bc₁

−λaY^bc₁ −ζaZ^bc₁ , (19) where the subscript 1 on the right-hand side indicates that only terms with density matrices up to first order are to be included. For example, in the expressions

W^bc₁ =D^bF^cD+D^bFD^c+DF^bD^c+DF^bc₁ D, (20) Z^bc₁ =D^bS^cD+D^bSD^c+DS^bD^c+DS^bcD (21) there are no terms containingD^bc. In Ref.43it is shown how to rewrite the above expression in a more symmetric form

L^abc_1,1 =^TrE_1,1^abc−(SW)^abc_1W −S^aW^bc₁ −λaY^bc₁ −ζaZ^bc₁. (22) Here, a prime on the subscript means that the subscript refers to the maximum order of differentiation ofS,D, andF(rather than the order ofD), for example,

W^bc₁ =D^bF^cD+D^bFD^c+DF^bD^c, (23) Z^bc₁ =D^bS^cD+D^bSD^c+DS^bD^c. (24) Even though the two first termsE_1,1^abcand (SW)^abc_1W in Eq.(22) are explicitly symmetric inabc, the remaining terms are not.

It is possible to symmetrize the expression but this does not give any computational benefits.

4. Quartic force constants

Following the same procedure as for the cubic force con- stants, it is possible to derive the following expression for the fourth-order energy derivative:⁴³

L^abcd_2,1 =^TrE_2,1^abcd−(SW)^abcd_1W

−S^adW^bc₁ −S^acW^bd₁ −S^abW^cd₁ −S^aW^bcd₁

−λ^daY^bc₁ −λ^caY^bd₁ −λ^baY^cd₁ −λaY^bcd₁

−ζ^daZ^bc₁ −ζ^caZ^bd₁ −ζ^baZ^cd₁ −ζaZ^bcd₁ . (25) As expected, we need up to second-order perturbed density matrices and first-order multipliers. By retaining the second- order density matrices involving the perturbations bcd, it is possible to eliminate the first-order multipliers⁴³

L^abcd_1,2 =^TrE_1,2^abcd−(SW)^abcd_2W −S^aW^bcd₂ −λaY^bcd₂ −ζaZ^bcd₂ . (26) In this notation, the latter expression appears more com- pact, but if one expands the different terms, one finds that Eqs.(25)and(26)are of similar complexity. Finally, we note that neither expression is explicitly symmetric in the pertur- bation labels.

5. Perturbed density matrices

To evaluate the expressions for the energy derivatives, we need first- and second-order perturbed density matrices. Since the unperturbed density matrix satisfies the idempotency con- dition and the SCF equations, the perturbed density can be obtained by differentiating Eqs.(2) and(7). In other words, the matrix D^b is the solution to the simultaneous equations Z^b =0andY^b=0:

DSD^b+DS^bD+D^bSD−D^b=0, (Z^b=0), (27)

FDS^b+FD^bS+F^bDS−SDF^b−SD^bF−S^bDF=0,

(Y^b =0). (28)

We rewrite Eq. (27)by collecting only terms containing D^b on one side, yielding

DSD^b+D^bSD−D^b =N, N=Z^b|D^b=0, (29) which has a solution of the general form

D^b =Dp+Dh, (30) Dp=NSD+DSN−N, (31) D_h =XSD−DSX, (32) whereD^bhas been partitioned into a particular partDpand a homogeneous partD_h. In the homogeneous equation,N=0 and the equation is automatically satisfied by the ansatz D_h=XSD−DSX.

To determine X, we use the differentiated SCF equa- tion. Inserting Eq.(30)into Eq.(28)and collecting all terms containing Xon the left, we arrive at the coupled-perturbed

Kohn–Sham equations (or more generally the linear response equations)

E^[2]X=M, (33) where we can identify the electronic HessianE^[2] in the AO basis

E^[2]X= ∂F

∂D^TDh

DS−SD ∂F

∂D^TDh

+FDhS−SDhF, (34) and the right-hand side

M=Y^b|D^b=D^b_p. (35) In the same way, the second-order perturbed density ma- trixD^bc can be determined from the equationsZ^bc =0 and Y^bc=0. The resulting equations have the same structure as in Eqs.(30)and(33), the matricesNandMnow being

N=Z^bc|D^bc=0, (36)

M=Y^bc|D^bc=D^bc_p. (37) To summarize, we must solve one set of linear response equa- tions for each perturbed density matrix, where the right-hand side depends on (perturbed) density matrices of lower orders.

We refer to Ref.43for further details.

B. Evaluation of exchange–correlation contributions We employ an exchange–correlation energyE_xcdefined as the integral over a local function_xc(r) that depends on the densityn(r) and its Cartesian gradient∇n(r):

E_xc=

dr_xc(n(r),∇n(r))=

dr_xc(ρ(r)), (38) where

n(r)=^Tr(r)D, (39)

∇n(r)=^Tr(∇(r))D. (40) To simplify notation, we henceforth collect the density vari- ablesn(r) and∇n(r) in a generalized density vectorρ(r). This notation also simplifies a generalization of our implementa- tion to other density variables such as the kinetic-energy den- sity in meta-GGA functionals and density variables in spin DFT and current DFT.

The exchange–correlation energy and the exchange–

correlation potential matrix are integrated on a numerical grid defined by a set of suitably chosen grid points ri and grid weightswi, according to

E_xc≈

i

wi_xc(ρ(ri)), (41)

(Fxc)μν ≈

i

wi

∂_xc(ρ(ri))

∂ρ(ri)

∂ρ(ri)

∂D_νμ

=

i

w_iv_xc(ri)(ρ)μν(ri). (42)

When differentiating the exchange–correlation energy and potential matrix we ignore the contribution from the grid- weight derivatives. The importance of grid-weight deriva- tives in the evaluation of geometrical derivatives at the DFT level has been discussed by Baker et al.⁵³ and Johnson and Frisch.¹⁷The extension to higher order is not straightforward, and for this reason we use very large grids in order to mini- mize the errors arising from the lack of grid-weight deriva- tive contributions, and the quality of the results has been verified by test calculations against numerically calculated derivatives.

For the implementation of DFT analytic cubic and quar- tic force constants, we need up to fourth-order geometri- cal derivatives of the exchange–correlation energy density (_xc^a, ^ab_xc, _xc^abc, and _xc^abcd) and up to second-order geomet- ric derivatives of the exchange–correlation potential matrix contributions (v^a_xc and v^ab_xc). The exchange–correlation en- ergy density derivatives are evaluated using the following expressions:

_xc^a =∂_xc

∂ρ ρ^a, (43)

_xc^ab= ∂_xc

∂ρ ρ^ab+∂²_xc

∂ρ² ρ^aρ^b, (44)

^abc_xc = ∂_xc

∂ρ ρ^abc+∂²_xc

∂ρ² [ρ^aρ^bc+ρ^bρ^ac+ρ^cρ^ab] + ∂³_xc

∂ρ³ ρ^aρ^bρ^c, (45)

_xc^abcd =∂_xc

∂ρ ρ^abcd+∂²_xc

∂ρ² [ρ^aρ^bcd+ρ^bρ^acd + ρ^cρ^abd+ρ^dρ^abc]

+ ∂³_xc

∂ρ³ [ρ^aρ^bρ^cd+ρ^aρ^cρ^bd+ρ^aρ^dρ^bc + ρ^bρ^cρ^ad+ρ^bρ^dρ^ac+ρ^cρ^dρ^ab] + ∂⁴_xc

∂ρ⁴ ρ^aρ^bρ^cρ^d, (46) where the arguments of densities, functional derivatives, and overlap distributions have been omitted for notational clarity.

In our code, the contractions of the functional deriva- tive vectors with the perturbed generalized density vectors are not explicitly programmed. Instead, we obtain the perturbed exchange–correlation energy densities _xc directly from the XCFUNprogram^47,54 by forming a generalized density Tay- lor series expansion (ρ, ρ^a, ρ^b, ρ^ab, . . . ), which is inter- nally contracted with the density functional Taylor expan- sion. This approach significantly reduces the complexity of the exchange–correlation integrator.

If the total energy had been variational with respect to the density ρ, then, according to the 2n +1 rule, we would only need the first-order (second-order) perturbed densities for the cubic (quartic) force field. In our case, the energy is

variational with respect to the AO density matrix D, which means that we still need the second- and third-order (third- and fourth-order) perturbed densities.

This dependence is given as

ρ^a=^TraρD+ρD^a, (47) ρ^ab=^Trabρ D+^aρD^b+^bρD^a+ρD^ab, (48)

ρ^abc=^Trabcρ D

+^abρ D^c+^acρD^b+^bcρD^a +^aρD^bc+^bρD^ac+^cρD^ab

+ρD^abc, (49) ρ^abcd =^Trabcdρ D

+^abcρ D^d+^abdρ D^c+^acdρ D^b+^bcdρ D^a +^abρ D^cd+^acρD^bd+^adρ D^bc+^bcρD^ad +^bdρ D^ac+^cdρ D^ab

+^aρD^bcd+^bρD^acd+^cρD^abd +^dρD^abc +ρD^abcd. (50) Depending on the perturbation order, many of the above terms are omitted when applying the 2n + 1 rule to the density matrix. Note that the omitted terms are not zero by themselves, but only in combination with non-exchange–

correlation terms containing the same density matrices. Here we have used

^aρ = −2^a,0ρ , (51) ^abρ =2

^ab,0ρ +^a,bρ

, (52)

^abcρ = −2

^abc,0ρ +^ab,cρ +^ac,bρ +^a,bcρ

, (53)

^abcdρ =2

^abcd,0ρ +^abc,dρ +^abd,cρ +^ab,cdρ

+ ^acd,bρ +^ac,bdρ +^ad,bcρ +^a,bcdρ

, (54) collecting

(^p,q)μν =χ_μ^∗pχ_ν^q, (55) (∇^p,q)μν =

∇χ_μ^∗p χ_ν^q+χ_μ^∗p

∇χ_ν^q , (56) into the generalized overlap distribution vector (^p,qρ )μν.

Having discussed exchange–correlation energy density contributions, we now turn to the exchange–correlation po- tential matrix contributions. The perturbations can either act on the generalized overlap distribution or on the functional derivative term, giving

[vxc(ρ)μν]^a=v_xc^a(ρ)μν+v_xc

^a_{ρ μν}, (57) [vxc(ρ)μν]^ab=v^ab_xc(ρ)μν+v_xc^a

^b_{ρ μν} +v^b_xc

^a_{ρ μν}+v_xc

^ab_{ρ μν}, (58)

where we have used

v_xc^a = ∂²_xc

∂ρ² ρ^a, (59)

v_xc^ab=∂³_xc

∂ρ³ ρ^aρ^b+∂²_xc

∂ρ² ρ^ab. (60) Finally, we note that an efficient implementation of the density evaluation and matrix distribution routines is essen- tial, bearing in mind the large number of terms that need to be evaluated. We evaluate both the densities and the matrix el- ements in a blocked manner, allowing mathematical matrix–

matrix multiplication libraries to be used in conjunction with efficient pre-screening techniques.

III. COMPUTATIONAL DETAILS

To calculate the cubic and quartic force constants, the re- cursive implementation⁴⁹ of the open-ended response-theory framework by Thorvaldsen et al.⁴³ has been used, as pro- vided by the OPENRSP program package. We use the DAL-

TONprogram package⁵⁵ as a backend for the calculation of undifferentiated integrals and the unperturbed and perturbed density matrices, which are obtained with the linear response solver of Jørgensenet al.⁵⁶The calculation of properties as- sociated with one-electron integrals was carried out using the GEN1INTlibrary,⁵⁷ building on the flexible integral evalua- tion scheme of Gao and co-workers.⁴⁵The differentiated two- electron integrals were mainly calculated using Thorvaldsen’s CGTO-DIFF-ERI code,⁵⁸ which uses the scheme of Reine et al.for the evaluation of differentiated two-electron integrals using solid-harmonic Gaussians,⁴⁶ but some of the lower- order contributions were calculated using DALTON. The dif- ferentiated exchange–correlation energy and potential contri- butions needed for the cubic and quartic force constants were calculated using the XCFUNlibrary,⁵⁴ which uses automatic differentiation for evaluating the derivatives of the exchange–

correlation energy.⁴⁷ We have used an in-house integra- tor to perform the integration of the exchange–correlation contributions.

Cubic and quartic force constants in the Cartesian ba- sis were calculated at the HF and DFT levels of theory for methane, ethane, benzene, and aniline. For methane, we per- formed a basis-set convergence study using the 6-31G⁵⁹ and the correlation-consistent basis sets⁶⁰ of double-, triple-, and quadruple-zeta quality (cc-pVDZ, cc-pVTZ, and cc-pVQZ).

For the other molecules, we have used the cc-pVTZ basis set for ethane and the cc-pVDZ basis set for benzene and aniline.

In order to explore the sensitivity of the results to the choice of exchange–correlation functional, both the BLYP^61–64 and the B3LYP⁶⁵functionals have been used.

For the HF and B3LYP calculations, the geometry was optimized and the molecular Hessian was calculated at the DFT (B3LYP) level of theory with the DALTON program package,⁵⁵ using the same basis as in the anharmonic force field calculations. The B3LYP Hessian was used in the vibra- tional analysis—both for evaluating the harmonic vibrational frequencies and for transforming the anharmonic force con- stants to a reduced normal coordinate basis²⁶ before evalu-

ating the fundamental frequencies (vide infra). Although not consistent, this approach circumvents the well-known defi- ciencies of the HF method for harmonic frequencies and al- lows us to get a better impression of the quality of the HF cubic and quartic force constants. For the calculations involv- ing the BLYP functional, the geometry optimization, the vi- brational analysis, and the cubic and quartic force constants were calculated using this functional, allowing us to com- pare directly the results obtained using the BLYP and B3LYP functionals.

In the calculations, we have converged the coupled- perturbed Kohn–Sham equations to a relative norm of 10⁻⁶, observing no problems with convergence of the response equations. To reduce the errors arising from the lack of grid- weight derivative contributions, we have used an ultrafine grid with a radial quadrature accuracy of 2 ×10⁻¹⁵ and with an angular expansion order of 64.

From the cubic and quartic force constants, anharmonic frequency corrections were calculated using the generalized vibrational second-order perturbation (GVPT2) model,^66,67in which terms that are too large because of Fermi resonances are excluded from the perturbational treatment⁶⁸ and treated variationally.⁴¹ The threshold criteria for the identification of Fermi resonances are those used by Bloino and Barone⁶⁹ ex- cept for ethane, where the threshold for the Martin parameters was increased to 1.5 cm⁻¹from the default value of 1 cm⁻¹, to avoid a splitting of degenerate modes due to unevenly dis- tributed interactions between these modes and a different set of two degenerate modes, which would otherwise lead to an unsymmetric identification of Fermi resonances. We refer to Refs.41and69for more details about the GVPT2 model and the treatment of Fermi resonances. All rotational effects, as described by the rotational constants and the Coriolis cou- pling constants in the GVPT2 scheme, are disregarded in the present work.

IV. RESULTS AND DISCUSSION

In TablesI–IV, we have listed the calculated fundamen- tal frequencies of methane, ethane, benzene, and aniline, re- spectively. Regarding the basis-set dependence of the anhar- monic corrections, the methane results in Table Iindicate it is rather weak, with small differences between 6-31G and cc- pVQZ anharmonic corrections, the largest difference between the HF/6-31G and HF/cc-pVQZ results being 6 cm⁻¹(4%).

Regarding the differences between the various levels of electronic-structure theory, we see from TableIthat the dif- ference between the HF and the B3LYP corrections (both calculated using B3LYP geometries, harmonic frequencies, and normal coordinates) are not large for methane, the ab- solute value of the B3LYP corrections being on average smaller than 10%. From the results for the larger molecules in Tables II–IV, this behaviour appears to be a general trend, but with some discrepancies being slightly larger than 10%.

Also, in a few cases, the B3LYP anharmonic corrections are larger than the corresponding HF corrections—for modes 8 and 12 in benzene and for mode 13 in aniline, for instance, the B3LYP correction is substantially more negative than the HF correction. These differences are mostly the result of

TABLE I. Harmonic fundamental vibrational frequenciesω, corrected fundamental frequenciesν, and anhar- monic vibrational correctionsδfor methane. All values are in cm⁻¹.

Mode ωB3LYP νHF δHF νB3LYP δB3LYP νBLYP δBLYP ωBLYP ωexp νexpa

6-31G

1 3165 2999 −166 3011 −154 2933 −157 3090

2 3043 2910 −132 2920 −122 2847 −126 2973

3 1601 1552 −48 1557 −44 1521 −45 1566

4 1403 1357 −45 1362 −41 1325 −42 1367

cc-pVDZ

1 3146 2977 −169 2988 −158 2906 −162 3068

2 3025 2887 −138 2892 −133 2817 −137 2954

3 1530 1484 −46 1488 −42 1451 −43 1494

4 1309 1264 −46 1268 −41 1233 −42 1275

cc-pVTZ

1 3129 2971 −158 2981 −148 2906 −152 3058

2 3027 2900 −127 2904 −122 2836 −127 2963

3 1559 1511 −48 1514 −44 1481 −45 1526

4 1341 1294 −47 1298 −43 1267 −44 1311

cc-pVQZ

1 3127 2967 −160 2979 −148 2903 −152 3055 3156.8 3022.5

2 3025 2896 −129 2902 −122 2833 −127 2960 3025.5 2920.9

3 1558 1510 −48 1514 −44 1481 −44 1524 1582.7 1532.4

4 1340 1293 −47 1298 −43 1267 −43 1310 1367.4 1308.4

aExperimental data taken from Ref.74and ordered by decreasing frequency.

differences between the HF and B3LYP values for the asso- ciated diagonal (iiii) quartic force constants; for mode 12 in benzene, differences in the semidiagonal (iijj) quartic force constants are also important. Overall, the BLYP anharmonic corrections are in good agreement with the B3LYP correc- tions, but with some exceptions. For mode 2 in benzene, for example, the anharmonic correction is substantially less negative with the BLYP exchange–correlation functional, be- cause of different identifications of Fermi resonances at dif- ferent levels of theory. In general, however, the BLYP anhar- monic corrections are slightly more negative than the B3LYP corrections.

Among the different levels of theory applied here, the B3LYP results are in best agreement with the experimental fundamental frequencies. The listed HF results are of com- parable quality but have been obtained at the B3LYP geome- try and are based on the B3LYP harmonic vibrational analy- sis; they would have been considerably worse had they been based on HF quantities alone. In any case, the calculation of HF anharmonic corrections based on DFT geometries and harmonic frequencies should be a viable approach in many cases. Likewise, we expect harmonic frequencies calculated at higher levels of theory—for instance, at the coupled-cluster level of theory—to perform well in combination with SCF an- harmonic corrections; indeed, such an approach has been used in earlier works.^38,70,71

The BLYP results consistently show the poorest agree- ment with experiment. However, much of the discrepancy arises from inaccurate harmonic frequencies. For all sys- tems, the BLYP corrections are mostly close to the HF and B3LYP corrections—however, for methane, the BLYP cor- respondence to experiment for the harmonic frequencies is clearly inferior to the B3LYP correspondence. This is further

accentuated when we note that, for the high-frequency modes in methane, the derived experimental anharmonic corrections are in general smaller than the calculated ones, suggesting that the experimental harmonic frequencies are underestimated. In general, the agreement between the B3LYP and BLYP anhar- monic corrections is good, supporting the notion that it is the harmonic frequencies that are poorly described by the BLYP functional.

In order to get a better global understanding of the per- formance of the different computational levels, we have in TablesII–IValso collected the mean absolute errors (MAEs) for the different computational levels compared to experi- ment. We note that the harmonic B3LYP frequencies in gen- eral are about 2.5% off the experimental frequencies, in line with the recommended scaling factors often used for B3LYP calculations of 0.9679.⁷² We note that for aniline, the MAE is somewhat larger than for the other molecules, but this is largely due to mode 30 which displays a MAE of almost 25%. Excluding this mode, the MAE for the B3LYP harmonic frequencies is 2.4%. Including anharmonic corrections to the B3LYP harmonic frequencies, either at the HF or B3LYP lev- els of theory, brings the MAE down to about 1% with the HF error being slightly larger and the B3LYP error slightly smaller than 1%, the differences in general being small. How- ever, once again mode 30 in aniline is an interesting case, clearly showing the superiority of the B3LYP method over the HF method for difficult cases, as the HF fundamental fre- quency for this mode is off by 23% from the experimental value, B3LYP instead being only 0.8% off the experimental data.

The BLYP model in contrast provides a much bet- ter MAE than the B3LYP model for harmonic frequencies.

However, this agreement is fortuitous and when adding the

TABLE II. Harmonic fundamental vibrational frequenciesω, corrected fundamental frequenciesν, and anhar- monic vibrational correctionsδfor ethane using the cc-pVTZ basis set. All values are in cm⁻¹.

Mode ωB3LYP νHF δHF νB3LYP δB3LYP νBLYP δBLYP ωBLYP νexpa

1 3093 2945 −148 2953 −140 2875 −145 3020 2977.7

2 3068 2923 −145 2932 −136 2854 −141 2994 2955.0

3 3025 2867 −159 2870 −155 2800 −158 2958 2920

4 3024 2867 −158 2868 −156 2797 −159 2956 2915

5 1507 1458 −48 1462 −45 1427 −46 1473 1471.6

6 1503 1452 −51 1456 −47 1422 −48 1469 1468.1

7 1423 1387 −36 1391 −32 1352 −33 1385 1388.4

8 1413 1376 −37 1379 −34 1346 −35 1381 1379.2

9 1223 1188 −34 1191 −31 1159 −32 1191 1190

10 995 969 −27 972 −23 934 −25 958 994.8

11 827 823 −4 821 −5 802 −6 809 821.6

12 305 267 −38 273 −32 265 −32 297 289

Mean absolute error relative to experiment in percent

2.81 1.55 1.22 3.80 1.17

aExperimental data taken from Ref.75and ordered by decreasing frequency.

anharmonic corrections to the BLYP data, the MAE actually increases, from about 1% to 3%–4%. If we instead add the BLYP corrections to the B3LYP harmonic frequencies, the MAE becomes comparable to that obtained with HF (about 1%).

The differences between the anharmonic corrections ob- tained at various levels of theory are relatively small. How- ever, for certain spectroscopic processes that have recently received increased attention—for example, the doubly vi-

brationally enhanced four-wave-mixing using two incident infrared lasers discussed in Ref. 73—the principal two- dimensional spectroscopic features may consist of closely spaced peaks separated by a distance related to the an- harmonic coupling between the modes involved, in addi- tion to lower-order contributions. The shape of such fea- tures can be very sensitive to the values of the anharmonic corrections, putting higher demands on the accuracy in the calculated anharmonic corrections for two-dimensional than

TABLE III. Harmonic fundamental vibrational frequenciesω, corrected fundamental frequenciesν, and anhar- monic vibrational correctionsδfor benzene using the cc-pVDZ basis set. All values are in cm⁻¹.

Mode ωB3LYP νHF δHF νB3LYP δB3LYP νBLYPa δBLYP ωBLYP νexpb

1 3198 3047 −151 3054 −144 2960 −155 3115 3073.942

2 3187 3029 −159 3034 −153 2975 −130 3104 (3057)

3 3171 3027 −144 3035 −136 2944 −144 3088 3056.7

4 3161 2991 −169 2996 −165 2901 −177 3078 3064.3674

5 1645 1596 −49 1599 −46 1538 −45 1583 1600.9764

6 1506 1475 −32 1476 −30 1430 −31 1462 1483.9854

7 1364 1331 −33 1333 −31 1294 −33 1327 (1350)

8 1356 1346 −10 1330 −26 1298 −30 1328 1309.4

9 1186 1171 −15 1171 −15 1140 −16 1156 1177.776

10 1162 1152 −10 1150 −12 1123 −12 1136 1149.7

11 1059 1038 −21 1039 −20 1007 −21 1028 1038.2670

12 1022 993 −29 982 −39 943 −46 989 993.071

13 1019 1002 −17 1004 −15 971 −16 986 (1010)

14 1013 997 −16 999 −14 973 −15 987 (990)

15 986 961 −26 959 −28 920 −31 952 (967)

16 866 844 −22 843 −23 814 −25 839 847.1

17 723 709 −14 705 −18 681 −21 702 (707)

18 691 678 −13 677 −14 656 −16 672 673.97465

19 618 612 −6 611 −7 596 −7 603 608.13

20 414 404 −10 404 −10 391 −11 402 (398)

Mean absolute error relative to experiment in percent

2.43 0.83 0.76 3.30 1.13

aModes 7 and 8 and modes 13 and 14 were switched after an analysis of the normal coordinates to agree with the B3LYP ordering.

bExperimental data taken from Ref.76and ordered by decreasing frequency.