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 kernels. 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 calculated 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.
I. INTRODUCTION
Vibrational spectroscopy provides a rich and diverse source of information about molecular structure and functionality. For this reason, methods for calculating molecular vibrational spectra were developed already in the early days of quantum-chemical computations, to help elucidate molecular structure and to provide insight into experimental observations. Already in 1958, Bratoz recognized the benefits of calculating the forces acting on the nuclei of a molecule in an analytic manner.1 However, the breakthrough in terms of computing molecular gradients and force fields came with the efficient implementation of molecular gradients by Pulay in 1969.2
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 complexity increases significantly when going to the molecular 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 density matrix is needed in this case. Thomson and Swanstrøm presented the first implementation of molecular Hessians in 1973.3 In the early 1980s, implementations of molecular Hessians 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 geometrical derivatives were derived and implemented for a wide range of correlated wave functions.6–15 More recently, implementations 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 calculations for correlated wave functions, Schaefer, Handy, and co-workers extended the evaluation of geometrical derivatives of the HF energy to third22,23 and fourth24 orders. Despite the importance of cubic and quartic force fields for determining, for instance, anharmonic corrections to vibrational frequencies, Fermi resonances,25 rotation–vibration constants,26 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 differentiation of forces and force constants at highly correlated levels 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 fields37,38 and of anharmonic corrections to intensities,39 thus demonstrating that anharmonic effects can now be studied straightforwardly for large and complex molecular systems. In their studies, geometrical 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); recently, we have calculated analytically quintic and sextic force fields at the HF level of theory.42
In this work, we describe the analytic calculation of cubic 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 framework of an atomic-orbital (AO) based, open-ended quasi-energy response-theory formalism described by Thorvaldsen et al.,43 which for force fields reduces to regular energy-derivative theory in the AO basis,44 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.45 The evaluation of geometrical derivatives of the two-electron repulsion integrals follows the approach of Reine, Tellgren, and Helgaker,46 expanding solid-harmonic Gaussians directly in Hermite Gaussians. We furthermore extend automatic differentiation of exchange–correlation kernels47 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 corrections 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 calculations, while Sec. IV presents and discusses the results. Finally, in Sec. V, we give some concluding remarks and perspectives for the analytic calculation of higher-order properties that involve geometrical distortions.
II. THEORY
We here present the theory behind our AO-based implementation of DFT cubic and quartic force fields, the AO-formulation ensuring that the approach is suitable for linear-scaling methodology.48 The approach builds on the general AO-based framework for time- and perturbation-dependent basis sets by Thorvaldsen et al.,43 here applied to time-independent perturbations. We note that, even though explicit 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.49 Compared with the earlier implementations by Schaefer, Handy, and co-workers,22–24 our 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, providing 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 differentiation to evaluate the higher-order exchange–correlation kernel derivatives.47
A. AO-based energy derivative theory
We follow the notation of the AO-based response theory for self-consistent-field (SCF) methods with time- and perturbation-dependent basis sets by Thorvaldsen et al.,43 specialized to static perturbations. In Kohn–Sham DFT, the energy can in the AO basis be written as
In this expression, “Tr” indicates that the trace of the matrix products on the right-hand side of the equation is taken, Exc[ρ(D)] is the exchange–correlation energy, which is a functional of the generalized density vector ρ (see Sec. II B), hnuc is the nuclear repulsion energy, and D is the AO density matrix. The density matrix fulfills the idempotency relation
We have in Eqs. (1) and (2) also introduced the one-electron integral matrix, h, the two-electron integral matrix constructed from D with γ fractional exchange, Gγ(D), and the overlap matrix, S, whose elements are given by
The χμ are spherical-harmonic Gaussian AOs and the summation in Eq. (3) is over atomic nuclei at RK and with charge ZK. The two-electron integrals are defined in the conventional manner as
with integration over all spin and spatial coordinates.
The optimized Kohn–Sham density fulfills the SCF condition
where the Kohn–Sham matrix F is defined as
with
Following Refs. 43,49,50, we take as our starting point for the generation of higher-order derivatives the energy-gradient Lagrangian
where the superscript a indicates differentiation with respect to an applied perturbation of strength ɛa. In Eq. (10) we have introduced the energy-weighted density matrix
as well as matrices that represent the constraints on the unperturbed reference state—in particular, the idempotency and SCF-state matrices, respectively,
For an optimized SCF state, Eqs. (2) and (7) can be written compactly as Z = 0 and Y = 0. Additionally, we have introduced the Lagrange multipliers λa and ζa, respectively, for these constraints. In Ref. 43, it is shown that Eq. (10) is variational in D if the zeroth-order multipliers are defined as
The subscript a on the multipliers does not indicate differentiation, merely the relation to
1. Molecular gradient
Let us first consider the first-order geometrical derivatives of the molecular energy. In this case, Eq. (10) simplifies to Pulay's expression from 19692
where no derivatives of the Lagrange multipliers are required because of the 2n + 2 rule.50 We note that the molecular gradient 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 to b, keeping only terms that fulfill the 2n + 1 and 2n + 2 rules, we obtain for the molecular Hessian the expression
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
where superscripts denote total derivatives. The subscripts (k, n) specify the maximum order of the perturbed density matrix D: to order k for collections of perturbations involving perturbation a, and to order n for collections of perturbations not involving perturbation a. The notation nW for the SW term specifies in a similar manner the maximum order of differentiation of W in this term, as dictated by the value of n.
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
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
there are no terms containing Dbc. In Ref. 43 it is shown how to rewrite the above expression in a more symmetric form
Here, a prime on the subscript means that the subscript refers to the maximum order of differentiation of S, D, and F (rather than the order of D), for example,
Even though the two first terms
4. Quartic force constants
Following the same procedure as for the cubic force constants, it is possible to derive the following expression for the fourth-order energy derivative:43
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 multipliers43
In this notation, the latter expression appears more compact, 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 perturbation 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 condition and the SCF equations, the perturbed density can be obtained by differentiating Eqs. (2) and (7). In other words, the matrix Db is the solution to the simultaneous equations Zb = 0 and Yb = 0:
We rewrite Eq. (27) by collecting only terms containing Db on one side, yielding
which has a solution of the general form
where Db has been partitioned into a particular part
To determine X, we use the differentiated SCF equation. Inserting Eq. (30) into Eq. (28) and collecting all terms containing X on the left, we arrive at the coupled-perturbed Kohn–Sham equations (or more generally the linear response equations)
where we can identify the electronic Hessian E[2] in the AO basis
and the right-hand side
In the same way, the second-order perturbed density matrix Dbc can be determined from the equations Zbc = 0 and Ybc = 0. The resulting equations have the same structure as in Eqs. (30) and (33), the matrices N and M now being
To summarize, we must solve one set of linear response equations for each perturbed density matrix, where the right-hand side depends on (perturbed) density matrices of lower orders. We refer to Ref. 43 for further details.
B. Evaluation of exchange–correlation contributions
We employ an exchange–correlation energy
where
To simplify notation, we henceforth collect the density variables n(r) and ∇n(r) in a generalized density vector ρ(r). This notation also simplifies a generalization of our implementation to other density variables such as the kinetic-energy density 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 weights wi, according to
When differentiating the exchange–correlation energy and potential matrix we ignore the contribution from the grid-weight derivatives. The importance of grid-weight derivatives in the evaluation of geometrical derivatives at the DFT level has been discussed by Baker et al.53 and Johnson and Frisch.17 The extension to higher order is not straightforward, and for this reason we use very large grids in order to minimize the errors arising from the lack of grid-weight derivative 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 quartic force constants, we need up to fourth-order geometrical derivatives of the exchange–correlation energy density (
where the arguments of densities, functional derivatives, and overlap distributions have been omitted for notational clarity.
In our code, the contractions of the functional derivative vectors with the perturbed generalized density vectors are not explicitly programmed. Instead, we obtain the perturbed exchange–correlation energy densities
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
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
collecting
into the generalized overlap distribution vector
Having discussed exchange–correlation energy density contributions, we now turn to the exchange–correlation potential matrix contributions. The perturbations can either act on the generalized overlap distribution or on the functional derivative term, giving
where we have used
Finally, we note that an efficient implementation of the density evaluation and matrix distribution routines is essential, bearing in mind the large number of terms that need to be evaluated. We evaluate both the densities and the matrix elements 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 recursive implementation49 of the open-ended response-theory framework by Thorvaldsen et al.43 has been used, as provided by the OpenRSP program package. We use the Dalton program package55 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ørgensen et al.56 The calculation of properties associated with one-electron integrals was carried out using the Gen1Int library,57 building on the flexible integral evaluation scheme of Gao and co-workers.45 The differentiated two-electron integrals were mainly calculated using Thorvaldsen's Cgto-Diff-Eri code,58 which uses the scheme of Reine et al. for the evaluation of differentiated two-electron integrals using solid-harmonic Gaussians,46 but some of the lower-order contributions were calculated using Dalton. The differentiated exchange–correlation energy and potential contributions needed for the cubic and quartic force constants were calculated using the XCFun library,54 which uses automatic differentiation for evaluating the derivatives of the exchange–correlation energy.47 We have used an in-house integrator to perform the integration of the exchange–correlation contributions.
Cubic and quartic force constants in the Cartesian basis were calculated at the HF and DFT levels of theory for methane, ethane, benzene, and aniline. For methane, we performed a basis-set convergence study using the 6-31G59 and the correlation-consistent basis sets60 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 BLYP61–64 and the B3LYP65 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,55 using the same basis as in the anharmonic force field calculations. The B3LYP Hessian was used in the vibrational analysis—both for evaluating the harmonic vibrational frequencies and for transforming the anharmonic force constants to a reduced normal coordinate basis26 before evaluating the fundamental frequencies (vide infra). Although not consistent, this approach circumvents the well-known deficiencies of the HF method for harmonic frequencies and allows us to get a better impression of the quality of the HF cubic and quartic force constants. For the calculations involving the BLYP functional, the geometry optimization, the vibrational analysis, and the cubic and quartic force constants were calculated using this functional, allowing us to compare 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−6, 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−15 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,67 in which terms that are too large because of Fermi resonances are excluded from the perturbational treatment68 and treated variationally.41 The threshold criteria for the identification of Fermi resonances are those used by Bloino and Barone69 except for ethane, where the threshold for the Martin parameters was increased to 1.5 cm−1 from the default value of 1 cm−1, to avoid a splitting of degenerate modes due to unevenly distributed 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. 41,69 for more details about the GVPT2 model and the treatment of Fermi resonances. All rotational effects, as described by the rotational constants and the Coriolis coupling constants in the GVPT2 scheme, are disregarded in the present work.
IV. RESULTS AND DISCUSSION
In Tables I–IV, we have listed the calculated fundamental frequencies of methane, ethane, benzene, and aniline, respectively. Regarding the basis-set dependence of the anharmonic corrections, the methane results in Table I indicate 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−1 (4%).
Harmonic fundamental vibrational frequencies ω, corrected fundamental frequencies ν, and anharmonic vibrational corrections δ for methane. All values are in cm−1.
| 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 |
| 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 |
Experimental data taken from Ref. 74 and ordered by decreasing frequency.
Harmonic fundamental vibrational frequencies ω, corrected fundamental frequencies ν, and anharmonic vibrational corrections δ for ethane using the cc-pVTZ basis set. All values are in cm−1.
| 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 | |||||
| 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 | |||||
Experimental data taken from Ref. 75 and ordered by decreasing frequency.
Harmonic fundamental vibrational frequencies ω, corrected fundamental frequencies ν, and anharmonic vibrational corrections δ for benzene using the cc-pVDZ basis set. All values are in cm−1.
| 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 | |||||
| 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 | |||||
Modes 7 and 8 and modes 13 and 14 were switched after an analysis of the normal coordinates to agree with the B3LYP ordering.
Experimental data taken from Ref. 76 and ordered by decreasing frequency.
Harmonic fundamental vibrational frequencies ω, corrected fundamental frequencies ν, and anharmonic vibrational corrections δ for aniline using the cc-pVDZ basis set. All values are in cm−1.
| Mode . | ωB3LYP . | νHF . | δHF . | νB3LYP . | δB3LYP . | νBLYPa . | δBLYP . | ωBLYP . | νexpb . |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 3626 | 3457 | −169 | 3452 | −174 | 3326 | −183 | 3509 | 3485 |
| 2 | 3528 | 3371 | −156 | 3368 | −160 | 3244 | −169 | 3413 | 3401 |
| 3 | 3198 | 3052 | −146 | 3059 | −139 | 2972 | −143 | 3115 | 3094 |
| 4 | 3179 | 3025 | −154 | 3032 | −148 | 2939 | −157 | 3096 | 3074 |
| 5 | 3173 | 3034 | −140 | 3042 | −132 | 2950 | −140 | 3090 | 3050 |
| 6 | 3157 | 3011 | −146 | 3019 | −138 | 2935 | −139 | 3074 | 3040 |
| 7 | 3156 | 2996 | −160 | 3002 | −154 | 2908 | −165 | 3073 | 3013 |
| 8 | 1665 | 1622 | −43 | 1626 | −40 | 1564 | −40 | 1604 | 1618 |
| 9 | 1635 | 1588 | −48 | 1592 | −43 | 1528 | −43 | 1572 | 1603 |
| 10 | 1635 | 1592 | −43 | 1596 | −39 | 1548 | −39 | 1586 | 1590 |
| 11 | 1530 | 1495 | −34 | 1497 | −32 | 1448 | −33 | 1481 | 1501 |
| 12 | 1496 | 1463 | −33 | 1466 | −31 | 1420 | −31 | 1451 | 1470 |
| 13 | 1369 | 1349 | −19 | 1338 | −31 | 1306 | −35 | 1341 | 1324 |
| 14 | 1349 | 1324 | −25 | 1324 | −24 | 1286 | −26 | 1312 | 1308 |
| 15 | 1311 | 1285 | −26 | 1286 | −25 | 1245 | −25 | 1270 | 1278 |
| 16 | 1188 | 1170 | −18 | 1170 | −18 | 1138 | −19 | 1157 | 1173 |
| 17 | 1166 | 1153 | −12 | 1153 | −12 | 1123 | −15 | 1138 | 1152 |
| 18 | 1134 | 1112 | −22 | 1115 | −19 | 1086 | −19 | 1105 | 1115 |
| 19 | 1068 | 1048 | −21 | 1053 | −15 | 1022 | −16 | 1037 | 1054 |
| 20 | 1048 | 1030 | −18 | 1031 | −17 | 998 | −18 | 1016 | 1028 |
| 21 | 1006 | 986 | −20 | 989 | −17 | 959 | −18 | 977 | 996 |
| 22 | 989 | 963 | −26 | 956 | −33 | 915 | −39 | 954 | 968 |
| 23 | 967 | 942 | −25 | 941 | −26 | 901 | −30 | 930 | 957 |
| 24 | 885 | 858 | −27 | 859 | −25 | 823 | −29 | 851 | 874 |
| 25 | 833 | 809 | −24 | 812 | −21 | 786 | −22 | 808 | 823 |
| 26 | 828 | 806 | −22 | 806 | −22 | 775 | −24 | 799 | 808 |
| 27 | 766 | 740 | −26 | 744 | −21 | 716 | −23 | 739 | 747 |
| 28 | 706 | 692 | −14 | 691 | −15 | 668 | −17 | 685 | 689 |
| 29 | 632 | 623 | −8 | 624 | −8 | 606 | −9 | 615 | 619 |
| 30 | 614 | 379 | −235 | 494 | −120 | 487 | −111 | 598 | (440/490) |
| 31 | 534 | 526 | −8 | 528 | −7 | 511 | −8 | 519 | 526 |
| 32 | 509 | 486 | −23 | 493 | −16 | 478 | −15 | 493 | 501 |
| 33 | 419 | 410 | −9 | 410 | −9 | 397 | −9 | 406 | 415 |
| 34 | 382 | 381 | −1 | 380 | −2 | 370 | −3 | 372 | 390 |
| 35 | 289 | 266 | −23 | 286 | −3 | 288 | −4 | 292 | |
| 36 | 223 | 217 | −6 | 218 | −5 | 210 | −5 | 216 | 233 |
| Mean absolute error relative to experiment in percentc | |||||||||
| 3.06 | 1.68 | 0.94 | 3.75 | 2.12 | |||||
| Mode . | ωB3LYP . | νHF . | δHF . | νB3LYP . | δB3LYP . | νBLYPa . | δBLYP . | ωBLYP . | νexpb . |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 3626 | 3457 | −169 | 3452 | −174 | 3326 | −183 | 3509 | 3485 |
| 2 | 3528 | 3371 | −156 | 3368 | −160 | 3244 | −169 | 3413 | 3401 |
| 3 | 3198 | 3052 | −146 | 3059 | −139 | 2972 | −143 | 3115 | 3094 |
| 4 | 3179 | 3025 | −154 | 3032 | −148 | 2939 | −157 | 3096 | 3074 |
| 5 | 3173 | 3034 | −140 | 3042 | −132 | 2950 | −140 | 3090 | 3050 |
| 6 | 3157 | 3011 | −146 | 3019 | −138 | 2935 | −139 | 3074 | 3040 |
| 7 | 3156 | 2996 | −160 | 3002 | −154 | 2908 | −165 | 3073 | 3013 |
| 8 | 1665 | 1622 | −43 | 1626 | −40 | 1564 | −40 | 1604 | 1618 |
| 9 | 1635 | 1588 | −48 | 1592 | −43 | 1528 | −43 | 1572 | 1603 |
| 10 | 1635 | 1592 | −43 | 1596 | −39 | 1548 | −39 | 1586 | 1590 |
| 11 | 1530 | 1495 | −34 | 1497 | −32 | 1448 | −33 | 1481 | 1501 |
| 12 | 1496 | 1463 | −33 | 1466 | −31 | 1420 | −31 | 1451 | 1470 |
| 13 | 1369 | 1349 | −19 | 1338 | −31 | 1306 | −35 | 1341 | 1324 |
| 14 | 1349 | 1324 | −25 | 1324 | −24 | 1286 | −26 | 1312 | 1308 |
| 15 | 1311 | 1285 | −26 | 1286 | −25 | 1245 | −25 | 1270 | 1278 |
| 16 | 1188 | 1170 | −18 | 1170 | −18 | 1138 | −19 | 1157 | 1173 |
| 17 | 1166 | 1153 | −12 | 1153 | −12 | 1123 | −15 | 1138 | 1152 |
| 18 | 1134 | 1112 | −22 | 1115 | −19 | 1086 | −19 | 1105 | 1115 |
| 19 | 1068 | 1048 | −21 | 1053 | −15 | 1022 | −16 | 1037 | 1054 |
| 20 | 1048 | 1030 | −18 | 1031 | −17 | 998 | −18 | 1016 | 1028 |
| 21 | 1006 | 986 | −20 | 989 | −17 | 959 | −18 | 977 | 996 |
| 22 | 989 | 963 | −26 | 956 | −33 | 915 | −39 | 954 | 968 |
| 23 | 967 | 942 | −25 | 941 | −26 | 901 | −30 | 930 | 957 |
| 24 | 885 | 858 | −27 | 859 | −25 | 823 | −29 | 851 | 874 |
| 25 | 833 | 809 | −24 | 812 | −21 | 786 | −22 | 808 | 823 |
| 26 | 828 | 806 | −22 | 806 | −22 | 775 | −24 | 799 | 808 |
| 27 | 766 | 740 | −26 | 744 | −21 | 716 | −23 | 739 | 747 |
| 28 | 706 | 692 | −14 | 691 | −15 | 668 | −17 | 685 | 689 |
| 29 | 632 | 623 | −8 | 624 | −8 | 606 | −9 | 615 | 619 |
| 30 | 614 | 379 | −235 | 494 | −120 | 487 | −111 | 598 | (440/490) |
| 31 | 534 | 526 | −8 | 528 | −7 | 511 | −8 | 519 | 526 |
| 32 | 509 | 486 | −23 | 493 | −16 | 478 | −15 | 493 | 501 |
| 33 | 419 | 410 | −9 | 410 | −9 | 397 | −9 | 406 | 415 |
| 34 | 382 | 381 | −1 | 380 | −2 | 370 | −3 | 372 | 390 |
| 35 | 289 | 266 | −23 | 286 | −3 | 288 | −4 | 292 | |
| 36 | 223 | 217 | −6 | 218 | −5 | 210 | −5 | 216 | 233 |
| Mean absolute error relative to experiment in percentc | |||||||||
| 3.06 | 1.68 | 0.94 | 3.75 | 2.12 | |||||
Modes 9 and 10 were switched after an analysis of the normal coordinates to correspond to the ordering of the B3LYP normal modes.
Experimental data taken from Ref. 77 and ordered by decreasing frequency.
Excluding mode 35 for which no experimental data exist, and using the experimental value of 490 cm−1 for mode 30.
Regarding the differences between the various levels of electronic-structure theory, we see from Table I that the difference between the HF and the B3LYP corrections (both calculated using B3LYP geometries, harmonic frequencies, and normal coordinates) are not large for methane, the absolute 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 differences between the HF and B3LYP values for the associated 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 corrections, but with some exceptions. For mode 2 in benzene, for example, the anharmonic correction is substantially less negative with the BLYP exchange–correlation functional, because of different identifications of Fermi resonances at different levels of theory. In general, however, the BLYP anharmonic 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 comparable quality but have been obtained at the B3LYP geometry and are based on the B3LYP harmonic vibrational analysis; 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 anharmonic corrections; indeed, such an approach has been used in earlier works.38,70,71
The BLYP results consistently show the poorest agreement with experiment. However, much of the discrepancy arises from inaccurate harmonic frequencies. For all systems, the BLYP corrections are mostly close to the HF and B3LYP corrections—however, for methane, the BLYP correspondence 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 anharmonic 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 performance of the different computational levels, we have in Tables II–IV also collected the mean absolute errors (MAEs) for the different computational levels compared to experiment. We note that the harmonic B3LYP frequencies in general are about 2.5% off the experimental frequencies, in line with the recommended scaling factors often used for B3LYP calculations of 0.9679.72 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 levels 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. However, 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 frequency 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 better MAE than the B3LYP model for harmonic frequencies. However, this agreement is fortuitous and when adding the 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 obtained at various levels of theory are relatively small. However, for certain spectroscopic processes that have recently received increased attention—for example, the doubly vibrationally 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 anharmonic coupling between the modes involved, in addition to lower-order contributions. The shape of such features 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 for one-dimensional spectroscopies. Our results suggest that the inclusion of correlation effects in the calculated anharmonic vibrational corrections by means of DFT may be worthwhile, in spite of the added computational cost. Finally, we note that the rotational effects, which have been neglected in our work, should not be sufficiently large to affect our conclusions.
V. SUMMARY AND OUTLOOK
We have presented the first analytic DFT implementation of cubic and quartic force constants. Our implementation is based on a recursive, open-ended, AO- and density-matrix-based energy derivative approach for SCF methods (HF and DFT). In combination with open-ended schemes for one- and two-electron integrals and for exchange–correlation contributions, this approach allows for a compact and efficient code for the analytic evaluation of anharmonic force constants.
We have demonstrated that the hybrid B3LYP exchange–correlation functional is superior to the generalized-gradient BLYP functional for calculating the fundamental frequencies of several small- and medium-sized molecules. However, this observed superiority of the B3LYP model is mostly the result of improvements in the harmonic force field—the differences between the B3LYP and BLYP fundamental frequencies are much larger than the differences in the anharmonic corrections. This effect is also reflected in the results obtained using the HF model: the HF anharmonic corrections are in good agreement with the B3LYP corrections when the HF calculations are based on the B3LYP structure and harmonic force field.
In future work, we intend to apply the calculated cubic and quartic force constants to obtain fundamental vibrational frequencies for use in various spectroscopic designs and to obtain anharmonic corrections to spectroscopic intensities—in particular with an eye towards multidimensional vibrational spectroscopies, where anharmonic effects are important.
ACKNOWLEDGMENTS
This work has received support from the Research Council of Norway through a Centre of Excellence Grant (Grant No. 179568/V30) and a research grant (Grant No. 191251), the European Research Council through a ERC Starting Grant (Grant No. 279619), and from the Norwegian Supercomputer Program through a grant of computer time. T.H. and U.E. acknowledge support from the European Research Council under the European Union Seventh Framework Program through the Advanced Grant ABACUS, ERC Grant Agreement No. 267683. The work was also supported by COST (Action CoDECS: COnvergent Distributed Environment for Computational Spectroscopy). We are grateful to Julien Bloino and Vincenzo Barone for helpful discussions.