Analytic cubic and quartic force fields using density-functional theory.

Analytic cubic and quartic force fields using density-functional theory Magnus Ringholm, Dan Jonsson, Radovan Bast, Bin Gao, Andreas J. Thorvaldsen, Ulf Ekström, Trygve Helgaker, and Kenneth Ruud Citation: The Journal of Chemical Physics 140, 034103 (2014); doi: 10.1063/1.4861003 View online: http://dx.doi.org/10.1063/1.4861003 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/3?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Analytic calculations of hyper-Raman spectra from density functional theory hyperpolarizability gradients J. Chem. Phys. 141, 134107 (2014); 10.1063/1.4896606 Communication: Rationale for a new class of double-hybrid approximations in density-functional theory J. Chem. Phys. 135, 101102 (2011); 10.1063/1.3640019 Double-hybrid density-functional theory made rigorous J. Chem. Phys. 134, 064113 (2011); 10.1063/1.3544215 Resonance Raman spectra of uracil based on Kramers–Kronig relations using time-dependent density functional calculations and multireference perturbation theory J. Chem. Phys. 120, 11564 (2004); 10.1063/1.1697371 Fundamental vibrational frequencies of small polyatomic molecules from density-functional calculations and vibrational perturbation theory J. Chem. Phys. 118, 7215 (2003); 10.1063/1.1561045

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.5.159.110 On: Tue, 07 Oct 2014 00:49:17

THE JOURNAL OF CHEMICAL PHYSICS 140, 034103 (2014)

Analytic cubic and quartic force fields using density-functional theory Magnus Ringholm,1 Dan Jonsson,1,2 Radovan Bast,3 Bin Gao,1 Andreas J. Thorvaldsen,1 Ulf Ekström,4 Trygve Helgaker,4 and Kenneth Ruud1 1

Centre for Theoretical and Computational Chemistry (CTCC), Department of Chemistry, University of Tromsø—The Arctic University of Norway, 9037 Tromsø, Norway 2 High Performance Computing Group, University of Tromsø—The Arctic University of Norway, 9037 Tromsø, Norway 3 Theoretical 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 4 Center 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 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. © 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 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 0021-9606/2014/140(3)/034103/11

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 ldoubling 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 highaccuracy 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 coworkers have shown that high-quality harmonic force fields in combination with DFT anharmonic corrections provide reliable estimates of anharmonic force fields37, 38 and of

140, 034103-1

© Author(s) 2014


034103-2

Ringholm et al.

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 finitedifference schemes and their numerical instability provide a strong motivation for developing analytic methods for thirdand 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 quasienergy response-theory formalism described by Thorvaldsen et al.,43 which for force fields reduces to regular energyderivative 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 energyderivative 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 AOformulation ensuring that the approach is suitable for linearscaling 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 timeindependent 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

J. Chem. Phys. 140, 034103 (2014)

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 1 Tr E (D) = hD + Gγ (D) D + Exc [ρ (D)] + hnuc . (1) 2 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 0 = DSD − D.

(2)

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 ZK 1 |χν , hμν = χμ | − ∇ 2 − (3) 2 |RK − r| K Gγμν (M) =

Mβα (gμναβ − γ gμβαν ),

(4)

αβ

Sμν = χμ |χν .

(5)

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 −1 ∗ dx1 dx2 χμ∗ (x1 ) χν (x1 ) r12 χρ (x2 ) χσ (x2 ) , (6) gμνρσ = with integration over all spin and spatial coordinates. The optimized Kohn–Sham density fulfills the SCF condition FDS − SDF = 0, where the Kohn–Sham matrix F is defined as ∂E = h + Gγ (D) + Fxc , F= ∂DT

(7)

(8)


034103-3

Ringholm et al.

with

J. Chem. Phys. 140, 034103 (2014)

the molecular Hessian the expression

Fxc =

dr

∂Exc (r) ∂ρ(r) . ∂ρ(r) ∂DT

(9)

Following Refs. 43, 49, and 50, we take as our starting point for the generation of higher-order derivatives the energy-gradient Lagrangian La defined as Tr

La = La (D, ζ a , λa ) =

∂E(D) − Sa W − λa Y − ζ a Z, (10) ∂εa

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 W = DFD,

(11)

as well as matrices that represent the constraints on the unperturbed 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 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 λa = D SD − DSD ,

(14)

ζ a = Fa DS + SDFa − FDSa − Sa DF − Fa .

(15)

a

a

The subscript a on the multipliers does not indicate differentiation, merely the relation to La . Since we require Eq. (10) to be variational in the density D and the multipliers λa and ζ a , we can take advantage of the 2n + 1 rule for the density and the 2n + 2 rule for the multipliers50, 51 when differentiating the energy gradient Lagrangian.

− Sab W − Sa Wb .

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 (16)

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

(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 d2 E Tr ab ab = Lab 0,1 = E0,1 − (SW)1W , dεa dεb

(18)

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. The Hessian expression in Eq. (17) is not explicitly symmetric in a and b (the numerical values, of course, are). As shown by Sellers, an explicitly symmetric formula can be advantageous from a numerical point of view.52 3. Cubic force constants

For the calculation of cubic force constants, the thirdorder 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 Tr

a bc abc Labc W − Sac Wb − Sab Wc − Sa Wbc 1,1 = (E )1 − S 1 bc − λa Ybc 1 − ζ a Z1 ,

(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 b c b c b c bc Wbc 1 = D F D + D FD + DF D + DF1 D, b c b c b c bc Zbc 1 = D S D + D SD + DS D + DS D

1. Molecular gradient

dE Tr ∂E = La = − Sa W, dεa ∂εa

2 dLa ∂ 2E d2 E Tr ∂ E = = Lab = + Db dεa dεb dεb ∂εa ∂εb ∂εa ∂DT

(20) (21)

there are no terms containing Dbc . In Ref. 43 it is shown how to rewrite the above expression in a more symmetric form Tr

abc abc a bc bc bc Labc 1,1 = E1,1 − (SW)1W − S W1 − λa Y1 − ζ a Z1 . (22)

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, b c b c b c Wbc 1 = D F D + D FD + DF D ,

(23)

b c b c b c Zbc 1 = D S D + D SD + DS D .

(24)

abc Even though the two first terms E1,1 and (SW)abc 1W in Eq. (22) are explicitly symmetric in abc, the remaining terms are not. It is possible to symmetrize the expression but this does not give any computational benefits.


Ringholm et al.

034103-4

J. Chem. Phys. 140, 034103 (2014)

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 Tr Labcd 2,1 =

abcd E2,1

−

(SW)abcd 1W

ac bd ab cd a bcd − S Wbc 1 − S W1 − S W1 − S W1 ad

c bd b cd bcd − λda Ybc 1 − λa Y1 − λa Y1 − λa Y1 c bd b cd bcd − ζ da Zbc 1 − ζ a Z1 − ζ a Z1 − ζ a Z1 .

(25)

As expected, we need up to second-order perturbed density matrices and first-order multipliers. By retaining the secondorder density matrices involving the perturbations bcd, it is possible to eliminate the first-order multipliers43 Tr

abcd abcd a bcd bcd bcd Labcd 1,2 = E1,2 − (SW)2W − S W2 − λa Y2 − ζ a Z2 . (26) 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: DSDb + DSb D + Db SD − Db = 0,

(Zb = 0),

Kohn–Sham equations (or more generally the linear response equations) E[2] X = M,

(33) [2]

where we can identify the electronic Hessian E in the AO basis ∂F ∂F DS − SD + FDh S − SDh F, D D E[2] X = h h ∂DT ∂DT (34) and the right-hand side M = Yb |Db =Dbp .

(35)

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 N = Zbc |Dbc =0 ,

(36)

M = Ybc |Dbc =Dbcp .

(37)

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 Exc defined as the integral over a local function xc (r) that depends on the density n(r) and its Cartesian gradient ∇n(r): Exc = dr xc (n(r), ∇n(r)) = dr xc (ρ(r)), (38)

(27) where Tr

n(r) = (r)D,

FDSb + FDb S + Fb DS − SDFb − SDb F − Sb DF = 0, (Yb = 0).

(28) b

We rewrite Eq. (27) by collecting only terms containing D on one side, yielding DSDb + Db SD − Db = N,

N = Zb |Db =0 ,

(29)

which has a solution of the general form

b

Db = Dp + Dh ,

(30)

Dp = NSD + DSN − N,

(31)

Dh = XSD − DSX,

(32)

where D has been partitioned into a particular part Dp and a homogeneous part Dh . In the homogeneous equation, N = 0 and the equation is automatically satisfied by the ansatz Dh = XSD − DSX. 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

Tr

∇n(r) = (∇(r))D.

(39) (40)

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 wi xc (ρ(ri )), (41) Exc ≈ i

(Fxc )μν ≈

i

=

wi

∂xc (ρ(ri )) ∂ρ(ri ) ∂ρ(ri ) ∂Dνμ

wi vxc (ri )( ρ )μν (ri ).

(42)

i


034103-5

Ringholm et al.

J. Chem. Phys. 140, 034103 (2014)

When differentiating the exchange–correlation energy and potential matrix we ignore the contribution from the gridweight 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 a ab abc abcd , xc , xc , and xc ) and up to second-order geomet(xc ric derivatives of the exchange–correlation potential matrix a ab and vxc ). The exchange–correlation encontributions (vxc ergy density derivatives are evaluated using the following expressions: ∂xc a ρ , ∂ρ

(43)

∂xc ab ∂ 2 xc a b ρ + ρ ρ , ∂ρ ∂ρ 2

(44)

a = xc

ab xc =

abc xc =

abcd xc =

∂ 3 xc a b c ρ ρ ρ , ∂ρ 3

Tr

(48)

Tr

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

+ aρ Dbc + bρ Dac + cρ Dab + ρ Dabc ,

(49)

Tr

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

+ aρ Dbcd + bρ Dacd + cρ Dabd + dρ Dabc + ρ Dabcd .

(45)

∂ 3 xc a b cd [ρ ρ ρ + ρ a ρ c ρ bd + ρ a ρ d ρ bc ∂ρ 3

+ ρ b ρ c ρ ad + ρ b ρ d ρ ac + ρ c ρ d ρ ab ] ∂ 4 xc a b c d ρ ρ ρ ρ , ∂ρ 4

a b b a ab ρ ab = ab ρ D + ρ D + ρ D + ρ D ,

ρ abc = abc ρ 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 derivative vectors with the perturbed generalized density vectors are not explicitly programmed. Instead, we obtain the perturbed exchange–correlation energy densities xc directly from the XCF UN program47, 54 by forming a generalized density Taylor series expansion (ρ, ρ a , ρ b , ρ ab , . . . ), which is internally contracted with the density functional Taylor expansion. 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

(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

+ ρ c ρ abd + ρ d ρ abc ]

+

(47)

Tr

∂xc abcd ∂ 2 xc a bcd ρ + [ρ ρ + ρ b ρ acd ∂ρ ∂ρ 2

+

ρ a = aρ D + ρ Da ,

ac cd ab + bd ρ D + ρ D

∂xc abc ∂ 2 xc a bc ρ + [ρ ρ + ρ b ρ ac + ρ c ρ ab ] ∂ρ ∂ρ 2 +

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

aρ = −2a,0 ρ ,

(51)

ab,0 + a,b ab ρ = 2 ρ ρ ,

(52)

abc,0 abc , + ab,c + ac,b + a,bc ρ = −2 ρ ρ ρ ρ

(53)

= 2 abcd,0 + abc,d + abd,c + ab,cd abcd ρ ρ ρ ρ ρ acd,b ac,bd ad,bc a,bcd , + ρ + ρ + ρ + ρ

(54)

collecting ( p,q )μν = χμ∗p χνq ,

(55)

(∇ p,q )μν = ∇χμ∗p χνq + χμ∗p ∇χνq ,

(56)

p,q

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 a ( ρ )μν + vxc aρ μν , (57) [vxc ( ρ )μν ]a = vxc b ab a [vxc ( ρ )μν ]ab = vxc ( ρ )μν + vxc ρ μν b aρ μν + vxc ab + vxc ρ μν ,

(58)


034103-6

Ringholm et al.

J. Chem. Phys. 140, 034103 (2014)

where we have used ∂ xc a ρ , ∂ρ 2

(59)

∂ 3 xc a b ∂ 2 xc ab ρ ρ + ρ . ∂ρ 3 ∂ρ 2

(60)

a = vxc

ab vxc =

2

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 O PEN RSP program package. We use the DAL TON 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 G EN 1I NT library,57 building on the flexible integral evaluation scheme of Gao and co-workers.45 The differentiated twoelectron integrals were mainly calculated using Thorvaldsen’s C GTO -D IFF -E RI 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 lowerorder 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 XCF UN 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 evalu-

ating 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 coupledperturbed 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 gridweight 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 and 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 ccpVQZ anharmonic corrections, the largest difference between the HF/6-31G and HF/cc-pVQZ results being 6 cm−1 (4%). 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


034103-7

Ringholm et al.

J. Chem. Phys. 140, 034103 (2014)

TABLE I. 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

1 2 3 4

3165 3043 1601 1403

2999 2910 1552 1357

− 166 − 132 − 48 − 45

3011 2920 1557 1362

6-31G − 154 − 122 − 44 − 41

2933 2847 1521 1325

− 157 − 126 − 45 − 42

3090 2973 1566 1367

1 2 3 4

3146 3025 1530 1309

2977 2887 1484 1264

− 169 − 138 − 46 − 46

2988 2892 1488 1268

cc-pVDZ − 158 − 133 − 42 − 41

2906 2817 1451 1233

− 162 − 137 − 43 − 42

3068 2954 1494 1275

1 2 3 4

3129 3027 1559 1341

2971 2900 1511 1294

− 158 − 127 − 48 − 47

2981 2904 1514 1298

cc-pVTZ − 148 − 122 − 44 − 43

2906 2836 1481 1267

− 152 − 127 − 45 − 44

3058 2963 1526 1311

1 2 3 4

3127 3025 1558 1340

2967 2896 1510 1293

− 160 − 129 − 48 − 47

2979 2902 1514 1298

cc-pVQZ − 148 − 122 − 44 − 43

2903 2833 1481 1267

− 152 − 127 − 44 − 43

3055 2960 1524 1310

a

ωexp

νexp a

3156.8 3025.5 1582.7 1367.4

3022.5 2920.9 1532.4 1308.4

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

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


034103-8

Ringholm et al.

J. Chem. Phys. 140, 034103 (2014)

TABLE II. 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 1 2 3 4 5 6 7 8 9 10 11 12

a

ωB3LYP

νHF

δHF

νB3LYP

δB3LYP

νBLYP

δBLYP

ωBLYP

νexp a

3093 3068 3025 3024 1507 1503 1423 1413 1223 995 827 305

2945 2923 2867 2867 1458 1452 1387 1376 1188 969 823 267

− 148 − 145 − 159 − 158 − 48 − 51 − 36 − 37 − 34 − 27 −4 − 38

2953 2932 2870 2868 1462 1456 1391 1379 1191 972 821 273

− 140 − 136 − 155 − 156 − 45 − 47 − 32 − 34 − 31 − 23 −5 − 32

2875 2854 2800 2797 1427 1422 1352 1346 1159 934 802 265

− 145 − 141 − 158 − 159 − 46 − 48 − 33 − 35 − 32 − 25 −6 − 32

3020 2994 2958 2956 1473 1469 1385 1381 1191 958 809 297

2977.7 2955.0 2920 2915 1471.6 1468.1 1388.4 1379.2 1190 994.8 821.6 289

2.81

1.55

Mean absolute error relative to experiment in percent 1.22 3.80

1.17

Experimental data taken from Ref. 75 and 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 obtained at various levels of theory are relatively small. However, 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 twodimensional 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

TABLE III. 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 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

a b

ωB3LYP

νHF

δHF

νB3LYP

δB3LYP

νBLYP a

δBLYP

ωBLYP

νexp b

3198 3187 3171 3161 1645 1506 1364 1356 1186 1162 1059 1022 1019 1013 986 866 723 691 618 414

3047 3029 3027 2991 1596 1475 1331 1346 1171 1152 1038 993 1002 997 961 844 709 678 612 404

− 151 − 159 − 144 − 169 − 49 − 32 − 33 − 10 − 15 − 10 − 21 − 29 − 17 − 16 − 26 − 22 − 14 − 13 −6 − 10

3054 3034 3035 2996 1599 1476 1333 1330 1171 1150 1039 982 1004 999 959 843 705 677 611 404

− 144 − 153 − 136 − 165 − 46 − 30 − 31 − 26 − 15 − 12 − 20 − 39 − 15 − 14 − 28 − 23 − 18 − 14 −7 − 10

2960 2975 2944 2901 1538 1430 1294 1298 1140 1123 1007 943 971 973 920 814 681 656 596 391

− 155 − 130 − 144 − 177 − 45 − 31 − 33 − 30 − 16 − 12 − 21 − 46 − 16 − 15 − 31 − 25 − 21 − 16 −7 − 11

3115 3104 3088 3078 1583 1462 1327 1328 1156 1136 1028 989 986 987 952 839 702 672 603 402

3073.942 (3057) 3056.7 3064.3674 1600.9764 1483.9854 (1350) 1309.4 1177.776 1149.7 1038.2670 993.071 (1010) (990) (967) 847.1 (707) 673.97465 608.13 (398)

2.43

0.83

Mean absolute error relative to experiment in percent 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.


034103-9

Ringholm et al.

J. Chem. Phys. 140, 034103 (2014)

TABLE IV. 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 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

ωB3LYP

νHF

δHF

νB3LYP

δB3LYP

νBLYP a

δBLYP

ωBLYP

νexp b

3626 3528 3198 3179 3173 3157 3156 1665 1635 1635 1530 1496 1369 1349 1311 1188 1166 1134 1068 1048 1006 989 967 885 833 828 766 706 632 614 534 509 419 382 289 223

3457 3371 3052 3025 3034 3011 2996 1622 1588 1592 1495 1463 1349 1324 1285 1170 1153 1112 1048 1030 986 963 942 858 809 806 740 692 623 379 526 486 410 381 266 217

− 169 − 156 − 146 − 154 − 140 − 146 − 160 − 43 − 48 − 43 − 34 − 33 − 19 − 25 − 26 − 18 − 12 − 22 − 21 − 18 − 20 − 26 − 25 − 27 − 24 − 22 − 26 − 14 −8 − 235 −8 − 23 −9 −1 − 23 −6

3452 3368 3059 3032 3042 3019 3002 1626 1592 1596 1497 1466 1338 1324 1286 1170 1153 1115 1053 1031 989 956 941 859 812 806 744 691 624 494 528 493 410 380 286 218

− 174 − 160 − 139 − 148 − 132 − 138 − 154 − 40 − 43 − 39 − 32 − 31 − 31 − 24 − 25 − 18 − 12 − 19 − 15 − 17 − 17 − 33 − 26 − 25 − 21 − 22 − 21 − 15 −8 − 120 −7 − 16 −9 −2 −3 −5

3326 3244 2972 2939 2950 2935 2908 1564 1528 1548 1448 1420 1306 1286 1245 1138 1123 1086 1022 998 959 915 901 823 786 775 716 668 606 487 511 478 397 370 288 210

− 183 − 169 − 143 − 157 − 140 − 139 − 165 − 40 − 43 − 39 − 33 − 31 − 35 − 26 − 25 − 19 − 15 − 19 − 16 − 18 − 18 − 39 − 30 − 29 − 22 − 24 − 23 − 17 −9 − 111 −8 − 15 −9 −3 −4 −5

3509 3413 3115 3096 3090 3074 3073 1604 1572 1586 1481 1451 1341 1312 1270 1157 1138 1105 1037 1016 977 954 930 851 808 799 739 685 615 598 519 493 406 372 292 216

3485 3401 3094 3074 3050 3040 3013 1618 1603 1590 1501 1470 1324 1308 1278 1173 1152 1115 1054 1028 996 968 957 874 823 808 747 689 619 (440/490) 526 501 415 390

3.06

Mean absolute error relative to experiment in percentc 1.68 0.94 3.75

233

2.12

a

Modes 9 and 10 were switched after an analysis of the normal coordinates to correspond to the ordering of the B3LYP normal modes. b Experimental data taken from Ref. 77 and ordered by decreasing frequency. c Excluding mode 35 for which no experimental data exist, and using the experimental value of 490 cm−1 for mode 30.

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 densitymatrix-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


034103-10

Ringholm et al.

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. 1 S.

Bratoz, Colloq. Int. C.N.R.S. 82, 287 (1958). Pulay, Mol. Phys. 17, 197 (1969). 3 K. Thomson and P. Swanstrøm, Mol. Phys. 26, 735 (1973). 4 Y. Osamura, Y. Yamaguchi, P. Saxe, M. A. Vincent, J. F. Gaw, and H. F. Schaefer, Chem. Phys. 72, 131 (1982). 5 Y. Osamura, Y. Yamaguchi, P. Saxe, D. J. Fox, M. A. Vincent, and H. F. Schaefer, J. Mol. Struct.: THEOCHEM 103, 183 (1983). 6 J. A. Pople, R. Krishnan, H. B. Schlegel, and J. S. Binkley, Int. J. Quantum Chem. 16(S13), 225 (1979). 7 R. Krishnan, H. B. Schlegel, and J. A. Pople, J. Chem. Phys. 72, 4654 (1980). 8 M. R. Hoffmann, D. J. Fox, J. F. Gaw, Y. Osamura, Y. Yamaguchi, R. S. Grev, G. B. Fitzgerald, H. F. Schaefer, P. J. Knowles, and N. C. Handy, J. Chem. Phys. 80, 2660 (1984). 9 P. Jørgensen and J. Simons, J. Chem. Phys. 79, 334 (1983). 10 T. U. Helgaker, J. Almlöf, H. J. A. Jensen, and P. Jørgensen, J. Chem. Phys. 84, 6266 (1986). 11 N. C. Handy, R. D. Amos, J. F. Gaw, J. E. Rice, and E. D. Simandiras, Chem. Phys. Lett. 120, 151 (1985). 12 R. J. Harrison, G. B. Fitzgerald, W. D. Laidig, and R. J. Bartlett, Chem. Phys. Lett. 124, 291 (1986). 13 H. Koch, H. J. A. Jensen, P. Jørgensen, T. Helgaker, G. E. Scuseria, and H. F. Schaefer, J. Chem. Phys. 92, 4924 (1990). 14 J. Gauss and J. F. Stanton, Chem. Phys. Lett. 276, 70 (1997). 15 P. G. Szalay, J. Gauss, and J. F. Stanton, Theor. Chem. Acc. 100, 5 (1998). 16 A. Komornicki and G. B. Fitzgerald, J. Chem. Phys. 98, 1398 (1993). 17 B. G. Johnson and M. J. Frisch, Chem. Phys. Lett. 216, 133 (1993). 18 P. Deglmann, F. Furche, and R. Ahlrichs, Chem. Phys. Lett. 362, 511 (2002). 19 J. F. Stanton and J. Gauss, Int. Rev. Phys. Chem. 19, 61 (2000). 20 R. Bast, U. Ekström, B. Gao, T. Helgaker, K. Ruud, and A. J. Thorvaldsen, Phys. Chem. Chem. Phys. 13, 2627 (2011). 21 T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen, and K. Ruud, Chem. Rev. 112, 543 (2012). 22 J. F. Gaw, Y. Yamaguchi, H. F. Schaefer, and N. C. Handy, J. Chem. Phys. 85, 5132 (1986). 23 S. M. Colwell, D. Jayatilaka, P. E. Maslen, R. D. Amos, and N. C. Handy, Int. J. Quantum Chem. 40, 179 (1991). 24 P. E. Maslen, D. Jayatilaka, S. M. Colwell, R. D. Amos, and N. C. Handy, J. Chem. Phys. 95, 7409 (1991). 2 P.

J. Chem. Phys. 140, 034103 (2014) 25 G.

H. Herzberg, Molecular Spectra and Molecular Structure: III. Electronic Spectra and Electronic Structure of Polyatomic Molecules (Van Nostrand Reinhold, Princeton, 1966). 26 I. M. Mills, in Molecular Spectroscopy: Modern Research, edited by K. N. Rao and C. W. Mathews (Academic, New York, 1972), pp. 115–140. 27 I. M. Mills, J. Phys. Chem. 80, 1187 (1976). 28 P.-O. Åstrand, K. Ruud, and D. Sundholm, Theor. Chem. Acc. 103, 365 (2000). 29 P.-O. Åstrand, K. Ruud, and P. R. Taylor, J. Chem. Phys. 112, 2655 (2000). 30 M. L. Grenier-Besson, J. Phys. Radium 21, 555 (1960). 31 M. L. Grenier-Besson, J. Phys. France 25, 757 (1964). 32 J. Breidung, W. Thiel, J. Gauss, and J. F. Stanton, J. Chem. Phys. 110, 3687 (1999). 33 C. Puzzarini, M. Biczysko, and V. Barone, J. Chem. Theory Comput. 7, 3702 (2011). 34 K. Aarset, A. G. Császár, E. L. Sibert, W. D. Allen, H. F. Schaefer, W. Klopper, and J. Noga, J. Chem. Phys. 112, 4053 (2000). 35 D. A. Matthews and J. F. Stanton, Mol. Phys. 107, 213 (2009). 36 M. E. Harding, J. Vázquez, J. Gauss, J. F. Stanton, and M. Kállay, J. Chem. Phys. 135, 044513 (2011). 37 D. Begue, P. Carbonniere, V. Barone, and C. Pouchan, Chem. Phys. Lett. 415, 25 (2005). 38 C. Puzzarini, M. Biczysko, and V. Barone, J. Chem. Theory Comput. 6, 828 (2010). 39 V. Barone, J. Bloino, C. A. Guido, and F. Lipparini, Chem. Phys. Lett. 496, 157 (2010). 40 V. Barone, J. Chem. Phys. 120, 3059 (2004). 41 V. Barone, J. Chem. Phys. 122, 014108 (2005). 42 J. Gong, J. F. Stanton, M. Ringholm, and K. Ruud (unpublished). 43 A. J. Thorvaldsen, K. Ruud, K. Kristensen, P. Jørgensen, and S. Coriani, J. Chem. Phys. 129, 214108 (2008). 44 H. Larsen, T. Helgaker, J. Olsen, and P. Jørgensen, J. Chem. Phys. 115, 10344 (2001). 45 B. Gao, A. J. Thorvaldsen, and K. Ruud, Int. J. Quantum Chem. 111, 858 (2011). 46 S. Reine, E. I. Tellgren, and T. Helgaker, Phys. Chem. Chem. Phys. 9, 4771 (2007). 47 U. Ekström, L. Visscher, R. Bast, A. J. Thorvaldsen, and K. Ruud, J. Chem. Theory Comput. 6, 1971 (2010). 48 S. Coriani, S. Høst, B. Jansik, L. Thøgersen, J. Olsen, P. Jørgensen, S. Reine, F. Pawlowski, T. Helgaker, and P. Sałek, J. Chem. Phys. 126, 154108 (2007). 49 M. Ringholm, D. Jonsson, and K. Ruud, “A general, recursive and openended response code,” J. Comput. Chem. (submitted). 50 K. Kristensen, P. Jørgensen, A. J. Thorvaldsen, and T. Helgaker, J. Chem. Phys. 129, 214103 (2008). 51 T. Helgaker and P. Jørgensen, Theor. Chim. Acta 75, 111 (1989). 52 H. Sellers, Int. J. Quantum Chem. 30, 433 (1986); 30, 713 (1986) (erratum). 53 J. Baker, J. Andzelm, A. Scheiner, and B. Delley, J. Chem. Phys. 101, 8894 (1994). 54 U. Ekström, “Xcfun library,” 2010, http://www.admol.org/xcfun. 55 “DALTON,” A molecular electronic structure program, Release 2.0 (2005), see http://www.kjemi.uio.no/software/dalton/dalton.html. 56 P. Jørgensen, H. J. A. Jensen, and J. Olsen, J. Chem. Phys. 89, 3654 (1988). 57 B. Gao and A. J. Thorvaldsen, “GEN1INT Version 0.2.1,” 2012, Gen1Int is a library to evaluate the derivatives of one-electron integrals with respect to geometry perturbations, external electric and magnetic fields, and total rotational angular momentum at zero fields with contracted rotational London atomic orbitals, and is released under the GNU Lesser General Public License. 58 A. J. Thorvaldsen, “Cgto-diff-eri,” 2012, A library for the evaluation of geometry-differentiated 2e-integrals (dn /dRnk ), released under the GNU Lesser General Public License. 59 W. J. Hehre, R. Ditchfield, and J. A. Pople, J. Chem. Phys. 56, 2257 (1972). 60 T. H. Dunning, J. Chem. Phys. 90, 1007 (1989). 61 A. D. Becke, J. Chem. Phys. 98, 5648 (1993). 62 A. D. Becke, Phys. Rev. A 38, 3098 (1988). 63 C. T. Lee, W. T. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988). 64 S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). 65 P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, J. Phys. Chem. 98, 11623 (1994). 66 H. H. Nielsen, Rev. Mod. Phys. 23, 90 (1951).


034103-11 67 D.

Ringholm et al.

Papoušek and M. R. Aliev, Molecular Vibrational-Rotational Spectra (Elsevier, Amsterdam, 1981). 68 J. M. L. Martin, T. J. Lee, P. R. Taylor, and J.-P. François, J. Chem. Phys. 103, 2589 (1995). 69 J. Bloino and V. Barone, J. Chem. Phys. 136, 124108 (2012). 70 P. Carbonniere, T. Lucca, C. Pouchan, N. Rega, and V. Barone, J. Comput. Chem. 26, 384 (2005). 71 D. Begue, P. Carbonniere, and C. Pouchan, J. Phys. Chem. A 109, 4611 (2005).

J. Chem. Phys. 140, 034103 (2014) 72 M.

P. Andersson and P. Uvdal, J. Phys. Chem. A 109, 2937 (2005). Kwak, S. Cha, M. Cho, and J. C. Wright, J. Chem. Phys. 117, 5675 (2002). 74 D. L. Gray and A. G. Robiette, Mol. Phys. 37, 1901 (1979). 75 J. L. Duncan, R. A. Kelly, G. D. Nivellini, and F. Tullini, J. Mol. Spectrosc. 98, 87 (1983). 76 L. Goodman, A. G. Ozkabak, and S. N. Thakur, J. Phys. Chem. 95, 9044 (1991). 77 J. C. Evans, Spectrochim. Acta 16, 428 (1960). 73 K.


Analytic solution for a quartic electron mirror.

Trihydrogen cation with neon and argon: structural, energetic, and spectroscopic data from quartic force fields.

Fundamental vibrational frequencies and spectroscopic constants of cis- and trans-HOCS, HSCO, and isotopologues via quartic force fields.

Force-field functor theory: classical force-fields which reproduce equilibrium quantum distributions.

Fitting of dihedral terms in classical force fields as an analytic linear least-squares problem.

Determining polarizable force fields with electrostatic potentials from quantum mechanical linear response theory.

Quartic force fields for excited electronic states: rovibronic reference data for the 1 (2)A' and 1 (2)A″ states of the isoformyl radical, HOC.

Dynamic properties of force fields.

Drag force in bimodal cubic-quintic nonlinear Schrödinger equation.

Frequency-dependent force fields for QMMM calculations.

Coupled modes, frequencies and fields of a dielectric resonator and a cavity using coupled mode theory.

A computational theory of visual receptive fields.

Secondary Structure of Rat and Human Amylin across Force Fields.

Evaluation of Representations and Response Models for Polarizable Force Fields.

Polarizable Force Fields: History, Test Cases, and Prospects.

Simulation of lipid bilayer self-assembly using all-atom lipid force fields.

Variational scheme to compute protein reaction pathways using atomistic force fields with explicit solvent.

Application of the queuing analytic theory in emergency department.

A method to study precision grip control in viscoelastic force fields using a robotic gripper.

Predicting water-to-cyclohexane partitioning of the SAMPL5 molecules using dielectric balancing of force fields.

Analytic theory of photoacoustic wave generation from a spheroidal droplet.

Force fields and scoring functions for carbohydrate simulation.

Improved model of hydrated calcium ion for molecular dynamics simulations using classical biomolecular force fields.

Pure-quartic solitons.