Second derivatives for approximate spin projection methods Lee M. Thompson and Hrant P. Hratchian Citation: The Journal of Chemical Physics 142, 054106 (2015); doi: 10.1063/1.4907269 View online: http://dx.doi.org/10.1063/1.4907269 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/5?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Semi-quartic force fields retrieved from multi-mode expansions: Accuracy, scaling behavior, and approximations J. Chem. Phys. 142, 154118 (2015); 10.1063/1.4918587 Dynamical second-order Bethe-Salpeter equation kernel: A method for electronic excitation beyond the adiabatic approximation J. Chem. Phys. 139, 154109 (2013); 10.1063/1.4824907 Communication: An efficient analytic gradient theory for approximate spin projection methods J. Chem. Phys. 138, 101101 (2013); 10.1063/1.4795429 Vibrational quenching of excitonic splittings in H-bonded molecular dimers: Adiabatic description and effective mode approximation J. Chem. Phys. 137, 184312 (2012); 10.1063/1.4763979 The ab initio model potential method: Third-series transition metal elements J. Chem. Phys. 110, 784 (1999); 10.1063/1.478046
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
THE JOURNAL OF CHEMICAL PHYSICS 142, 054106 (2015)
Second derivatives for approximate spin projection methods Lee M. Thompson and Hrant P. Hratchiana) Chemistry and Chemical Biology, University of California, Merced, California 95343, USA
(Received 27 November 2014; accepted 16 January 2015; published online 5 February 2015) The use of broken-symmetry electronic structure methods is required in order to obtain correct behavior of electronically strained open-shell systems, such as transition states, biradicals, and transition metals. This approach often has issues with spin contamination, which can lead to significant errors in predicted energies, geometries, and properties. Approximate projection schemes are able to correct for spin contamination and can often yield improved results. To fully make use of these methods and to carry out exploration of the potential energy surface, it is desirable to develop an efficient second energy derivative theory. In this paper, we formulate the analytical second derivatives for the Yamaguchi approximate projection scheme, building on recent work that has yielded an efficient implementation of the analytical first derivatives. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4907269]
I. INTRODUCTION
Density Functional Theory (DFT) is the most widespread class of methods used in quantum chemical studies of medium and large sized systems.1,2 Due to a mean-field cost and broad applicability, DFT strikes an excellent balance of computational cost with accuracy, allowing researchers to compute reaction energies, activation energies, and molecular properties with reasonable confidence across a broad set of chemical system classes. Despite the success of current-day models, DFT does suffer from a few well known systematic failings. Of specific interest here is spin contamination in open-shell species, which results from the single-determinant framework underlying Kohn–Sham DFT. Spin contamination can be particularly severe in systems exhibiting electronic strain, spin polarization, or broken symmetry electronic structures.3–10 It should be noted that, in principle, this limitation of DFT can be overcome using sophisticated multireference electronic structure methods. However, the steep computational scaling of such models rapidly renders them unfeasible for many applications.11 To solve spin-contamination issues within a singledeterminant model, it is possible to use spin projection operators to annihilate higher-order spin states.12–14 This forms a perturbation series, where the first order annihilates the (s + 1)th state, the second order annihilates the (s + 2)th state, etc.15,16 Thus, even at first order, analytical spin projection requires evaluation of double excitation integrals and so increases the cost beyond that of a mean-field approach. An alternative which avoids this issue is Approximate Projection (AP), such as those of Noodleman17 and Yamaguchi.18–20 Whereas AP energies have been used in studies for the last couple of decades, approximate AP gradient formulations were only reported in the mid-2000s.21–24 Most recently, our group reported21 the first complete and efficient gradient theory based on the z-vector scheme,25 allowing efficient a)Electronic mail: [email protected]
0021-9606/2015/142(5)/054106/8/$30.00
exploration of the potential energy surface (PES) in its spin pure state. While the availability of analytic forces enhances the study of contaminated and broken-symmetry states, second derivatives will be necessary for efficiently characterizing stationary points and carrying out geometry optimization on flat PESs. In this paper, we report a complete analytic second derivative theory and implementation for the AP model. This implementation is derived using the z-vector derivative scheme and is thus quite efficient, as it requires calculation of the molecular orbital (MO) density perturbation response to first order only. This exhibits the same overall cost scaling as self consistent field (SCF) second derivatives and requiring only minimal additional costs.
II. METHODS
AP approaches approximate the broken symmetry solutions of the Kohn-Sham or Hartree-Fock equations to be eigenfunctions of the Heisenberg spin Hamiltonian, which resolves the energy of the pure spin states contaminating the broken symmetry solution. In the limiting case, where there is a single contaminating spin state, which represents the current level of development for these methods, this requires two SCF calculations—one of the broken symmetry state and one of the higher energy contaminating state. Extending this beyond two states would require an extra SCF calculation for each additional contaminating state. These methods have the advantage over formal analytical projection methods, which involve factorial cost increase with the number of projected states. The AP method removes the (s + 1)th contaminant from the spin state of interest, referred to as the low spin (LS) state, using the energy expression EAP = αELS + (1 − α)EHS,
(1)
where ELS is the LS energy, EHS is the energy of the contaminating state, referred to as the high spin (HS) state, and α is 142, 054106-1
© 2015 AIP Publishing LLC
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-2
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
obtained from
2 x ⟨SLS ⟩ =
α= 2 ⟨SHS ⟩
2 ⟨SHS ⟩ − Sz,LS(Sz,LS + 1) 2 2 ⟨SHS ⟩ − ⟨SLS ⟩
µν
,
(2)
where and are the total electron spin squared expectation values of the HS and LS states, respectively, and Sz,LS is the z component of the total electron spin of the LS state.
+ −
−
−
as x x x , + (1 − α)EHS = α E˜LS EAP (E − E ) LS HS x x ⟨S 2 ⟩ x , E˜LS = ELS + 2 2 ⟨SHS⟩ − ⟨SLS ⟩ LS
( )x 2 x ⟨SLS ⟩ = Sz,LS(Sz,LS + 1) + n β − |Si j¯|2 ,
(5)
2 x ⟨SLS ⟩
µν
x Ω µν Sµν ,
(7)
where the intermediate matrices Υ, Λ, and Ω are given by Λ µν = Sµλ P¯λσ Sσν , (8) Υµν =
λσ
) x x Cν j¯ , Cµi Sµν Cν j¯ + Cµi Sµν Cνxj¯ + Cµi Sµν
(9)
Pµλ Sλσ P¯σν ,
(10)
λσ
where the over-bar indicates the β density matrix and its absence indicates the α density matrix. Expression (4) can now be re-written in a standard postSCF manner, which is amenable to z-vector substitution, x AP AP x E˜LS = Γµνλσ (µν|λσ) x + Pµν H µν µν
µνλσ
+
AP x W µν Sµν + VNx N ,
(11)
µν
where P AP is the sum of the SCF density matrix P and the AP correction density matrix P∆. Γ AP is given by AP ∆ ∆ Γµνλσ = Pµν Pλσ + Pµσ Pλν ,
(12)
AP ∆ ¯ ∆ ¯ Γ¯ µνλσ = P¯µν Pλσ + P¯µσ Pλν .
(13) AP
The energy weighted density matrix W is obtained AP from the sum of SCF and AP contributions, W µν = W µν ∆ AP − Ω µν + W µν . The term Ω is summed into W µν . This term is not common to post-SCF gradient theorems and is effectively a non-separable term arising from derivation in the AO basis. The AP energy weighted density matrix correction is most conventionally defined in the MO basis ∆ (p j||qi), (14) Wi∆j = − Λ µν Cµi Cν j − Ppq µν
pq
(6) where Sµν denotes the AO overlap matrix elements and C is the matrix of MO coefficients. We use the conventional notation— greek subscripts refer to AOs; roman subscripts refer to MOs; p, q, r, s . . . refer to all MOs; i, j, k, l . . . indicate occupied MOs; and a, b, c, d . . . indicate virtual MOs.26 Solving Eq. (6) requires the first order response of the density to the external perturbation(s), however, solution of the first-order coupled perturbed Hartree-Fock (CPHF) equations can be avoided by reformulating it as AO density matrix derivatives P x , back-transforming to the MO basis, and separating summations into occupied-occupied (oo), virtualvirtual (vv), occupied-virtual (ov), and virtual-occupied (vo) blocks. Then, using identities from CPHF theory (Pixj = −Sixj x x 27 and Pai = Pia ), we obtain
Sµσ Pσλ Sλν ,
λσ
Ω µν = 2
where n β is the number of β electrons in the LS system and Si j¯ is the α– β overlap matrix of the LS system, where the over-bar indicates a β spin orbital. Given that the first two terms in Eq. (5) are constants, only the last term requires special care. To obtain the α– β overlap derivative, we formulate it in the atomic orbital (AO) basis, giving
i j¯
Υµν (Cµ i¯Cν a¯ + Cµ a¯ Cν i¯)Pax¯ i¯
i¯a¯
(4)
i j¯
(
ia
µν
(3)
where the superscript denotes differentiation with respect to an external perturbation x. Note that it is assumed that the 2 HS state is spin pure, and ⟨SHS ⟩ is unchanged with geometric distortion. 2 x Thus, ⟨SLS ⟩ needs to be solved and depends upon the first order density response to the perturbation. While the value of ⟨S 2⟩ for DFT requires the two-particle density matrix, however, it has been shown that it is possible to approximate it using the same approach as Hartree-Fock (HF),
Si j¯
x Λ µν (Cµi Cνa + CµaCνi )Pai
µν
Following our previous work, the gradient can be written
Υµν Cµ i¯Cν j¯Si¯xj¯
i¯ j¯
µν
A. AP first derivatives
= −2
ij
µν
2 ⟨SLS ⟩
Λ µν Cµi Cν j Sixj
Wi¯∆j¯
=−
Υµν Cµ i¯Cν j¯ −
µν
¯ q¯i), ¯ Pp∆¯ q¯ ( p¯ j||
(15)
p¯ q¯
∆ ∆ Wai = −Pai ϵ i,
(16)
−Pa∆¯ i¯ϵ i¯,
(17)
Wa∆¯ i¯ ∆ Wab Wa∆¯ b¯
=
= 0,
(18)
= 0,
(19)
and once all blocks of W∆ are found, it is rotated to the AO basis ∆ ∆ W µν = Cµ pCνqW pq + Cµ p¯ Cν q¯ W p∆¯ q¯ . (20) pq
Pi∆j ,
p¯ q¯
Pi∆¯ j¯,
∆ Pab ,
Finally, the and Pa∆¯ b¯ terms are zero and the remaining ov and vo terms are obtained by invoking the
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-3
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
z-vector approach25 ∆ = L ai , (21) [(i j||ab) − (ib|| ja)]Pb∆j + (ϵ a − ϵ i )Pai
Cqx µ according to Cqx µ =
bj
u xpqCp µ .
(30)
p
¯ − (i¯b|| ¯ j¯a)]P ¯ a¯ b) [(i¯ j|| ¯ b∆¯ j¯ + (ϵ a¯ − ϵ i¯)Pa∆¯ i¯ = L a¯ i¯, (22) b¯ j¯
for Eqs. (21) and (22), ϵ p is the energy of MO p and the AP Lagrangian matrices L are L ai = − Λ µν (Cµi Cνa + CµaCνi ), (23)
The u x vectors are obtained as solutions to the first-order CPHF equations and used to derive the MO density derivative matrix with the relations Pixj = uix†j + uixj = −Sixj , x Pia x Pab
µν
L a¯ i¯ = −
Υµν (Cµ i¯Cν a¯ + Cµ a¯ Cν i¯).
(24)
µν
=
x uai ,
(32)
= 0.
(33)
The AP derivative density matrix can then be obtained from AP x x ∆x Pµν = Pµν + Pµν .
B. AP second derivatives
The second derivative term can be formulated in a similar manner to Eqs. (1) and (3) as xy EAP
xy xy = α E˜LS + (1 − α)EHS ,
(25)
where E˜ x y is given by xy xy E˜LS = E˜LS +
x x ( E˜LS − EHS ) 2 y ⟨S ⟩ , 2 2 ⟨SHS⟩ − ⟨SLS⟩ LS
Similarly, the AP effective two-particle density matrix derivative is given by AP x ∆x ∆ x ∆x ∆ x Γµνλσ = Pµν Pλσ + Pµν Pλσ − Pµσ Pλµ − Pµσ Pλν ,
(35)
∆ ¯x ∆x ¯ ∆ ¯x AP x ∆x ¯ Pλν , Pλµ − P¯µσ Pλσ − P¯µσ Pλσ + P¯µν Γ¯ µνλσ = P¯µν
(36)
xy xy E˜LS = ELS +
AP x x ∆x W µν = W µν − Ω xµν + W µν .
(26)
In Eq. (37), Ω is given by x x x ¯ Ω xµν = 2 Pµλ Sλσ P¯σν + Pµλ Sλσ P¯σν + Pµλ Sλσ Pσν ,
y y (ELS − EHS )
(ELS − EHS) 2 x y + 2 ⟨S ⟩ . 2 ⟨SHS⟩ − ⟨SLS ⟩ LS
and the components of W∆x are Wi∆x Λ xµν Cµi Cν j j =− µν
(ELS − EHS) 2 x y + 2 ⟨S ⟩ . 2 ⟨SHS⟩ − ⟨SLS ⟩ LS
+ +
µνλσ
APy x Pµν H µν +
µν
µν
AP x y H µν + Pµν
µν AP x y W µν Sµν
+ VNx yN .
−
µν
APy x W µν Sµν
µν
(29)
In addition to terms defined in the AP first derivatives, we now require the terms ΓAPy, PAPy, WAPy, as well as the one-electron, two-electron, and overlap integral second derivatives. These terms involve the vectors u xpq that describe the rotation of unperturbed orbitals Cp µ to perturbed orbitals
Λ µν (u xpi Cµ pCν j + Cµi u xp j Cν p )
p
∆x ∆ Ppq (ip|| jq) + Ppq (ip|| jq) x ,
(39)
pq
Wi¯∆x =− j¯
x Υµν Cµ i¯Cν j¯
µν
−
−
µν
(28)
In order to avoid the need to compute the second-order response of the MO density to a perturbation via the secondorder CPHF equations, Eq. (28) can be cast in a post-SCF second derivative form and z-vector substitution carried out, APy x y′ AP (µν|λσ) x + Γµνλσ (µν|λσ) x y E˜LS = Γµνλσ
− (27)
The first term in Eq. (27) is the second derivative of the spin-contaminated LS state, and the final term is dependent on the second derivative of the LS spin-squared expectation value
µνλσ
(38)
λσ
⟨S 2 ⟩ x 2 2 ⟨SHS ⟩ − ⟨SLS ⟩ LS (ELS − EHS) 2 y ⟩ + ⟨S 2 ⟩ x ⟨SLS 2 2 2 LS (⟨SHS ⟩ − ⟨SLS ⟩)
=
(37)
x
xy x where E˜LS is determined by Eq. (4) and E˜LS is found from direct differentiation of Eq. (4) to be
xy ELS
(34)
and the energy weighted density matrix derivative is
LS
x y′ E˜LS
(31)
Υµν (u xp¯ i¯Cµ p¯ Cν j¯ + Cµ i¯u xp¯ j¯Cν p¯ )
p¯
¯ ¯ j¯q) Pp∆x ¯ + Pp∆¯ q¯ (i¯p|| ¯ j¯q) ¯ x, ¯ q¯ (i p||
(40)
p¯ q¯ ∆x ∆x ∆ Wai = −Pai ϵ i − Pai (Fiix − Siix ϵ i ),
Wa∆x ¯ i¯
=
−Pa∆x ϵ − ¯ i¯ i¯
Pa∆¯ i¯(Fi¯xi¯ −
Si¯xi¯ϵ i¯),
(41) (42)
where we have used the relation ϵ ix = Fiix − ϵ i Siix . The derivative matrices Λ x and Υ x are given by x ¯ x x ¯x Λ xµν = Sµλ Pλσ Sσν + Sµλ P¯λσ Sσν + Sµλ Pλσ Sσν , (43) x Υµν
=
λσ
x x x x Sµλ Pλσ Sσν + Sµλ Pλσ Sσν + Sµλ Pλσ Sσν .
(44)
λσ ∆x Finally, the Pµν terms are determined by solving the set of response-like equations, ∆x [(i j||ab) − (ib|| ja)]Pb∆xj + (ϵ a − ϵ i )Pai = L ∆x (45) ai , bj
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-4
L. M. Thompson and H. P. Hratchian
x L ∆x ai = L ai −
J. Chem. Phys. 142, 054106 (2015)
x ∆ − Fixj )Pai , [(i j||ab) x − (ib|| ja) x ]Pb∆j − (Fab bj
(46)
¯ − ¯ a¯ b) [(i¯ j||
¯ j¯a)]P (i¯b|| ¯ b∆x ¯ j¯
+ (ϵ a¯ −
ϵ i¯)Pa∆x ¯ i¯
=
L ∆x , a¯ i¯
(47)
b¯ j¯
L ∆x = L ax¯ i¯ − a¯ i¯
¯ x − (i¯b|| ¯ j¯a) [(i j|| a¯ b) ¯ x ]Pb∆¯ j¯ − (Fa¯xb¯ − Fi¯xj¯)Pa∆¯ i¯, b¯ j¯
(48) where the Lagrangian derivative matrices are x L ai =− Λ xµν (Cµi Cνa + CµaCνi ) −
µν
Λ µν (u xpi Cµ pCνa
L ax¯ i¯
=−
+ Cµi u xpaCν p
+ Cµau xpi Cν p ),
x Υµν (Cµ i¯Cν a¯
(49)
+ Cµ a¯ Cν i¯)
µν
−
µν
Υµν (u xp¯ i¯Cµ p¯ Cν a¯ + Cµ i¯u xp¯ a¯ Cν p¯
p¯
+ u xp¯ a¯ Cµ p¯ Cν i¯ + Cµ a¯ u xp¯ i¯Cν p¯ ).
x L ai ) ( x† t t = − u C Λ C + C Λ C µ p µν νq µq µν ν p ∆x pi Wab p µν µν ( ) t t − u xt pa C µq Λ µν Cν p + C µ p Λ µν Cνq ,
(51)
p
+ u xpaCµ pCνi
Wi∆x j L x ia
p
µν
with u xpi at O(MON 2) cost. Focusing only on the second term of each equation, it is seen that the supermatrix (for the α spin case) is
(50)
Thus, by utilization of z-vector substitution, we only require the first-order response of the MO derivative coefficients to the perturbation through a set of CPHF equations in order to obtain AP second derivatives. This results in significant cost savings compared to the explicit calculation of the secondorder MO density response. C. Efficiency considerations for AP second derivatives
We note that implementation of Eqs. (25)–(50) as written above yields significant computation bottlenecks. As a result, the cost of an AP force constant calculation can require substantially more time than the combined time necessary for the constituent HS and LS frequency calculations. However, employing straightforward factorization techniques provides important cost improvements (vide infra). In what follows, we define M, N, O, and V to be the number of perturbations, AO basis functions, occupied MOs, and virtual MOs, respectively. As written, the second term of Eqs. (39), (40), (49), and (50) involves O(MON 4) work. This is because these equations represent the oo block (Eqs. (39) and (40)) and the ov block (Eqs. (49) and (50)) of two supermatrices, one for α and one for β spins. Due to the symmetry of these matrices and that the vv block is zero, we only need to solve for the oo and one of the ov blocks, leading to the O(ON) scaling on top of the O(M N 3) for matrix multiplications. This step can be reduced to a single O(N 3) step followed by individual O(ON 2) steps for each perturbation M. Furthermore, the number of matrix multiplications with u x terms can be reduced by recognizing the symmetry of Λ and Υ. The remaining terms are independent of perturbation and correspond to transformation of Λ and Υ into the MO basis. Thus, these matrices only require a single evaluation at O(N 3) cost, before looping over each perturbation and contracting
where, as described above, the symmetry of the matrices involved requires only explicit formation of the first summa∆x 1 x tion over p in Eq. (51) to obtain Wi∆x j , Wab , and /2 L ia . Second, using techniques derived in the solution of other post-SCF derivations,28,29 we avoid the need for the second interchange suggested by Eqs. (45) and (47). Although it is ∆x possible to directly form Pai , inspection of the terms with ∆x which Pai is contracted shows that they form the right-side of y the CPHF system of equations (commonly denoted as Q ai in CPHF theory), ( y y ∆x Pai hai − ϵ i Sai − Skyl (al||ik) kl
+
) ∆x y Q ai . Cµi Cνa Pλσ (µλ||νσ) y = Pai
(52)
µνλσ
The terms of Eqs. (45) and (47) can be regrouped to form the left-hand side of the CPHF equations (commonly denoted as A in CPHF theory) ) ( L ∆y ( ja||ab) − (ib|| ja) ai ∆x = , (53) Pai 1+ ϵa − ϵi ϵa − ϵi ∆x Pai
L ∆y ai = [1 + A]−1, ϵa − ϵi
(54)
y such that by multiplying through by Q ai , we are able to x generate an expression in terms of uai , ∆x y Pai Q ai
=
L ∆y ai [1
+ A]
−1
Q by j ϵb − ϵ j
∆x y x Pai Q ai = L ∆y ai uai ,
,
(55) (56)
x where L ∆y ai and uai have already been calculated prior to ∆x solving for Pai and can be contracted directly. This technique ∆ is the reverse of the interchange invoked to obtain Pai in the gradient derivation and avoids solution of an additional set of response equations in the calculation of force constants at the ∆x expense of never explicitly forming Pai .
III. NUMERICAL TESTS A. Computational scaling and timing
To illustrate the efficiency of our AP second derivative theory, this section describes numerical tests of computational scaling of the method with respect to system size. These calculations have been carried out on the diradical intermediate in the [2 + 2] cycloaddition of singlet oxygen with a terminal alkene chain (Fig. 1). The derivation
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-5
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
FIG. 1. Diradical intermediate in the [2 + 2] cycleaddition of singlet oxygen to a terminal alkene chain, with the α and β spin localization indicated by red and blue halos.
described above was implemented in a local development version of Gaussian30 and calculations were carried out using a production cluster equipped with dual 6-core 2.60 GHz Intel Xeon processors. Except where noted, each calculation was carried out in serial mode (using a single core) using up to 32 GB RAM. In order for our derivation to be useful, it is important that scaling of the method with respect to system size is kept to mean-field order and to minimize overhead relative to the sum of the two constituent SCF calculations. The minimum cost of an AP second derivative calculation is found in the limiting case that α is invariant to the perturbation(s). This xy xy corresponds to Eq. (25) in which E˜LS = ELS and so only SCF frequency calculations at both LS and HS multiplicities are required. In these calculations, the main expense is the solution of the first-order CPHF equations which grows as O(M 2OV ) in a direct formalism.31,32 Calculating the change in α with respect to the perturbation(s) involves computing the additional post-SCF terms, which increases the cost beyond that of the sum of the LS and HS calculations. These extra terms include the assembly of the Lagrangian derivative and calculation of the response density. Additionally, hybridDFT has additional cost beyond HF due to the evaluation of two electron integral derivatives contracted with the density over grid points.
FIG. 2. Comparison of elapsed time for AP second derivative calculation (red triangles) against LS + HS frequency calculations (blue squares) for varying number of basis functions. The number of perturbations is constant over all points. Calculations were carried out using B3LYP in serial mode on a production computer equipped with dual 6 core 2.60 GHz Intel Xeon processors.
To demonstrate the scaling of our program in terms of the number of perturbations, M, and the number of basis functions, N, we report two different sets of calculations. First, six different basis sets have been employed to calculate AP force constants for the compound shown in Fig. 1. This group of basis sets used yields a range of basis set size from 66 to 590, corresponding to a size range of the OV space in the LS sub-calculation from 14 400 to 59 904. Second, the 6-311G(d,p) basis set is used with a series of compounds based on that shown in Fig. 1. Each member of this set of compounds extends the alkane chain of the parent molecule by two carbon centers. The resulting group of molecules ranges from 26 to 68 atomic centers, corresponding to a range of 81–214 nuclear displacement perturbations included in force constant calculations. Calculation times for both sets of jobs using HF and B3LYP are reported in Table I. Figure 2 shows the scaling of AP second derivatives with number of basis functions for the molecule shown in Fig. 1 (red triangles) versus the scaling of the component SCF second derivative calculations (blue squares). This data are compiled from the top part of Table I, which keeps the
TABLE I. Parameters and elapsed times in seconds for frequency calculations using AP and the sum of the two constituent SCF second derivative calculations. Values are taken from calculations carried out using a single core. HS + LS SCF freq. Natoms
AP freq.
Basis set
Nbasis
M
OVLS
HF
B3LYP
HF
B3LYP
26
STO-3G 6-31G 6-311G 6-311G(d,p) 6-311++G (d,p) cc-pVTZ
66 122 178 286 342 590
81 81 81 81 81 81
1040 3280 5520 9840 12 080 22 000
11.6 61.2 218.9 1468.9 2856.6 35 972.1
306.3 756.6 1405.3 3787.3 7553.7 39 599.4
12.8 68.0 243.4 1708.7 3276.4 39 453.9
423.5 1014.5 1880.4 5101.3 9256.9 45 518.4
32 38 44 50 56 62 68
6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p)
348 410 472 534 596 658 720
96 114 132 150 168 196 214
14 400 19 824 26 112 33 264 41 280 50 160 59 904
3514.5 7846.6 13 101.8 19 710.2 28 203.7 36 050.7 46 398.8
6799.5 11 811.2 17 429.2 24 721.3 33 488.9 42 496.4 54 897.4
4126.7 8995.6 14 909.5 22 397.2 32 191.0 41 570.5 53 290.9
8870.6 14 682.7 21 834.8 31 095.7 42 346.1 50 903.9 63 721.2
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-6
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
TABLE II. Parameters and elapsed times in seconds for AP and SCF of sample calculations in parallel mode using a production computer equipped with dual 6-core 2.60 GHz Intel Xeon processors. HS + LS SCF freq. Natoms
26
AP freq.
Basis set
Nbasis
M
OVLS
HF
B3LYP
HF
B3LYP
STO-3G 6-31G 6-311G 6-311 + +G (d,p) cc-pVTZ
66 122 178 342 590
81 81 81 81 81
1040 3280 5520 12 080 22 000
10.1 16.4 39.5 390.3 4720.4
37.0 85.4 151.7 793.9 4859.9
10.6 21.8 47.8 704.7 5252.9
48.6 108.1 198.6 1078.0 5614.1
number of perturbations constant. From this figure, it can be seen that the overhead of the AP second derivatives is small compared to the underlying SCF calculations. In fact, the overhead as a fraction of the SCF second derivative calculations, which ranges from 10% to 40%, decreases as the OV space increases. This is due to the scaling of the CPHF solutions with the number of virtual orbitals, which grows with the size of basis. As this is common to both sets of calculations, its expense dominates the calculation and the relative cost of the post-SCF like terms rapidly decreases. For example, with 178 basis functions, solution of the CPHF equations in the B3LYP calculation accounts for 27% of job time. In contrast, at 590 basis functions solution of the CPHF work accounts for 65% of total CPU time. In the lower part of Table I, we show the scaling for calculations on the molecule in Fig. 1 using the 6-311G(d,p) basis set as the number of carbon atoms in the alkyl chain is increased by two carbons from 6 to 20. Therefore, this indicates how our method scales with regard to an increase in the number of perturbations, as well as the number of basis functions and OV size. This set of results indicates that the overhead of the AP method does not depend on the number of perturbations, giving the same trend as the results in which only the number of basis functions were increased. In all sets of results, the HF calculations were observed to have smaller overhead than the corresponding B3LYP calculations. This is because the formation of the Lagrangian derivative requires the contraction of the integral derivatives with the response density evaluated over grid points. In order to demonstrate that our AP second-derivative implementation performs well in parallel computing environments, Table II shows job times for sample calculations carried out using 12 processors. The total AP wall time speedup going from 1 to 12 processors ranges from 68% to 79% efficiency for B3LYP, which is comparable to standard DFT frequency calculations. The efficiency is somewhat lower for HF as this does not benefit from the highly parallelized exchange-correlation integration over grid points. Despite this, HF also gives valuable cost reduction when run in parallel mode. This indicates that the post-SCF components of the AP second derivative calculation parallelize as well as the components from SCF calculations. Figure 3 shows the same trend as the serial calculations for the overhead with regard to system size, indicating that in this case, it is also the solution of the CPHF equations which dominate at large N and OV .
B. Example calculations
In order to demonstrate how the application of AP second derivatives modifies the calculated force constants, we have run several calculations on potentially challenging radical triatomic systems (Table III).33 These results show that the broken symmetry solutions of LiO2 and NO2 have small amounts of spin contamination. For this reason, application of the AP method results in very small changes in the calculated vibrational frequencies. This is satisfying as both the broken symmetry and AP solutions give respectable agreement with Coupled Cluster Singles & Doubles (CCSD) calculations using the same basis set. However, the broken 2 ⟩ symmetry C+3 is more significantly contaminated, with ⟨SLS equal to 1.022. This causes UB3LYP to give shorter bond lengths than the equivalent CCSD calculation, as well as a significant deviation in the ω1(a1) vibrational mode. This mode corresponds to the bending mode, and application of AP causes a shift to higher frequency of 116 cm−1. This suggests that the gap between the contaminating state and the ground state varies along this coordinate and is close enough that this variance has an effect on the curvature of the broken symmetry solution. In the case of the other modes, either the change in energy gap is small, or the gap is too large to have any effect. Further investigation into the effects of second derivatives on vibrational modes and spectra is ongoing in our lab and will be reported in due course.
FIG. 3. Comparison of elapsed time for AP second derivative calculation (red triangles) against LS + HS frequency calculations (blue squares) for varying number of basis functions. The number of perturbations is constant over all points. Calculations were carried out in parallel mode using 12 cores on production computer equipped with dual 6-core 2.60 GHz Intel Xeon processors.
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-7
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
TABLE III. Geometric parameters and vibrational modes of example triatomic molecules calculated using UB3LYP and AP-B3LYP with cc-pVTZ compared to CCSD. 2A 2
LiO2
Method
RLiO
θ OLiO
2 ⟩ ⟨SLS
UB3LYP AP-B3LYP CCSD
1.747 1.746 1.757
45.2 45.3 44.6
0.757 ... ... 2B 2
ω 3 (a 1)
535 537 536
761 762 763
1170 1167 1204
C+3
RCC
θ CCC
2 ⟩ ⟨SLS
ω 1 (a 1)
ω 2 (b2)
ω 3 (a 1)
UB3LYP AP-B3LYP CCSD
1.309 1.315 1.317
66.6 65.2 65.4
1.022 ... ...
742 958 877
1238 1302 1297
1660 1692 1666
1
NO2
Method
RNO
θ ONO
2 ⟩ ⟨SLS
ω 1 (a 1)
ω 2 (a 1)
ω 3 (b 2)
UB3LYP AP-B3LYP CCSD
1.190 1.190 1.183
134.4 134.4 134.8
0.754 ... ...
769 769 786
1394 1394 1434
1696 1694 1763
IV. CONCLUSION
By casting the equations into the form of a postSCF model, we have been able to derive analytical second derivatives of the AP method. Application of interchange and reverse-interchange techniques has enabled implementation of an efficient algorithm that we have demonstrated has a small amount of overhead relative to two SCF frequency calculations. These second derivatives include the orbital relaxation terms that account for the change in the alphabeta overlap matrix under the perturbations considered and should be significantly more accurate than derivations which ignore this term. The availability of AP second derivatives will allow the study of spin pure states in a number of systems for which spin contamination adversely affects results and for which analytical second derivatives are desired. ACKNOWLEDGMENTS
We gratefully acknowledge Dr. Michael Frisch (Gaussian, Inc.), Dr. Giovanni Scalmani (Gaussian, Inc.), and Professor Marco Caricato (University of Kansas) for helpful discussions and the University of California, Merced for financial support. Burke, J. Chem. Phys. 136, 150901 (2012). D. Becke, J. Chem. Phys. 140, 18A301 (2014). 3J. L. Sonnenberg, H. B. Schlegel, and H. P. Hratchian, in Computational Inorganic and Bioinorganic Chemistry, edited by E. I. Solomon, R. A. Scott, and B. A. King (John Wiley & Sons, Ltd., 2009). 4F. Neese, Coord. Chem. Rev. 253, 526 (2009). 5S. E. Waller, J. E. Mann, D. W. Rothgeb, and C. C. Jarrold, J. Phys. Chem. A 116, 9639 (2012). 6D. Aradilla and C. Aleman, J. Phys. Org. Chem. 27, 867 (2014). 7P. Jost and C. van Wüllen, Phys. Chem. Chem. Phys. 15, 16426 (2013). 8J. A. Pople, P. M. W. Gill, and N. C. Handy, Int. J. Quantum Chem. 56, 303 (1995). 9U. Bozkaya, J. Chem. Theory Comput. 10, 4389 (2014). 10U. Bozkaya, J. Chem. Phys. 135, 224103 (2011). 2A.
ω 2 (a 1)
Method
2A
1K.
ω 1 (b2)
11H.
B. Schlegel and M. Frisch, in Theoretical and Computational Models for Organic Chemistry, edited by S. J. Formosinho, I. Csizmadia, and L. Arnaut (Kluwer Academic Publishers, 1991). 12P. O. Löwdin, Phys. Rev. 97, 1509 (1955). 13C. Sosa and H. B. Schlegel, Int. J. Quantum Chem. 32, 267 (1987). 14N. C. Handy and P. J. Knowles, J. Chem. Phys. 88, 6991 (1988). 15N. C. Handy and P. J. Knowles, J. Phys. Chem. 92, 3097 (1988). 16L. M. Thompson and H. P. Hratchian, J. Chem. Phys. 141, 034108 (2014). 17L. Noodleman, C. Y. Peng, D. A. Case, and J. M. Mouesca, Coord. Chem. Rev. 144, 199 (1995). 18K. Yamaguchi, F. Jensen, A. Dorigo, and K. N. Houk, Chem. Phys. Lett. 149, 537 (1988). 19K. Yamaguchi, M. Okumura, W. Mori, and J. Maki, Chem. Phys. Lett. 210, 201 (1993). 20S. Yamanaka, M. Okumura, M. Nakano, and K. Yamaguchi, J. Mol. Struct.: THEOCHEM 310, 205 (1994). 21H. P. Hratchian, J. Chem. Phys. 138, 101101 (2013). 22Y. Kitagawa, T. Saito, M. Ito, M. Shoji, K. Koizumi, S. Yamanaka, T. Kawakami, M. Okumura, and K. Yamaguchi, Chem. Phys. Lett. 442, 445 (2007). 23T. Saito and W. Thiel, J. Phys. Chem. A 116, 10864 (2012). 24T. Saito, S. Nishihara, Y. Kataoka, Y. Nakanishi, T. Matsui, Y. Kitagawa, T. Kawakami, M. Okumura, and K. Yamaguchi, Chem. Phys. Lett. 483, 168 (2009). 25N. C. Handy and H. F. Schaefer, J. Chem. Phys. 81, 5031 (1984). 26J. A. Pople and D. L. Beveridge, Approximate Molecular Orbital Theory (McGraw-Hill, Inc., 1970). 27K. Raghavachari, J. A. Pople, and H. B. Schlegel, Int. J. Quantum Chem. 16, 225 (1979). 28R. Cammi, B. Mennucci, C. Pomelli, C. Cappelli, S. Corni, L. Frediani, G. W. Trucks, and M. Frisch, Theor. Chem. Acc. 111, 66 (2004). 29G. W. Trucks, M. Frisch, M. Head-Gordon, H. B. Schlegel, E. A. Salter, and J. L. Andres, “An efficient theory and implementation of MP2 second derivatives” (unpublished). 30M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, B. Janesko, F. Lipparini, G. Zheng, J. L. Sonnenberg, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin,
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-8
L. M. Thompson and H. P. Hratchian
K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, P. V. Parandekar, N. J. Mayhall, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, “Gaussian development version, Revision H.38+,” Gaussian, Inc., 2010.
J. Chem. Phys. 142, 054106 (2015) 31M. Head-Gordon, J. A. Pople, and M. Frisch, Chem. Phys. 141, 189 (1990). 32L.
M. Thompson, A. Lasoroski, P. M. Champion, J. T. Sage, M. Frisch, J. J. van Thor, and M. J. Bearpark, J. Chem. Theory Comput. 10, 751 (2014). 33U. Bozkaya and C. D. Sherrill, J. Chem. Phys. 138, 184103 (2013).
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
THE JOURNAL OF CHEMICAL PHYSICS 142, 054106 (2015)
Second derivatives for approximate spin projection methods Lee M. Thompson and Hrant P. Hratchiana) Chemistry and Chemical Biology, University of California, Merced, California 95343, USA
(Received 27 November 2014; accepted 16 January 2015; published online 5 February 2015) The use of broken-symmetry electronic structure methods is required in order to obtain correct behavior of electronically strained open-shell systems, such as transition states, biradicals, and transition metals. This approach often has issues with spin contamination, which can lead to significant errors in predicted energies, geometries, and properties. Approximate projection schemes are able to correct for spin contamination and can often yield improved results. To fully make use of these methods and to carry out exploration of the potential energy surface, it is desirable to develop an efficient second energy derivative theory. In this paper, we formulate the analytical second derivatives for the Yamaguchi approximate projection scheme, building on recent work that has yielded an efficient implementation of the analytical first derivatives. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4907269]
I. INTRODUCTION
Density Functional Theory (DFT) is the most widespread class of methods used in quantum chemical studies of medium and large sized systems.1,2 Due to a mean-field cost and broad applicability, DFT strikes an excellent balance of computational cost with accuracy, allowing researchers to compute reaction energies, activation energies, and molecular properties with reasonable confidence across a broad set of chemical system classes. Despite the success of current-day models, DFT does suffer from a few well known systematic failings. Of specific interest here is spin contamination in open-shell species, which results from the single-determinant framework underlying Kohn–Sham DFT. Spin contamination can be particularly severe in systems exhibiting electronic strain, spin polarization, or broken symmetry electronic structures.3–10 It should be noted that, in principle, this limitation of DFT can be overcome using sophisticated multireference electronic structure methods. However, the steep computational scaling of such models rapidly renders them unfeasible for many applications.11 To solve spin-contamination issues within a singledeterminant model, it is possible to use spin projection operators to annihilate higher-order spin states.12–14 This forms a perturbation series, where the first order annihilates the (s + 1)th state, the second order annihilates the (s + 2)th state, etc.15,16 Thus, even at first order, analytical spin projection requires evaluation of double excitation integrals and so increases the cost beyond that of a mean-field approach. An alternative which avoids this issue is Approximate Projection (AP), such as those of Noodleman17 and Yamaguchi.18–20 Whereas AP energies have been used in studies for the last couple of decades, approximate AP gradient formulations were only reported in the mid-2000s.21–24 Most recently, our group reported21 the first complete and efficient gradient theory based on the z-vector scheme,25 allowing efficient a)Electronic mail: [email protected]
0021-9606/2015/142(5)/054106/8/$30.00
exploration of the potential energy surface (PES) in its spin pure state. While the availability of analytic forces enhances the study of contaminated and broken-symmetry states, second derivatives will be necessary for efficiently characterizing stationary points and carrying out geometry optimization on flat PESs. In this paper, we report a complete analytic second derivative theory and implementation for the AP model. This implementation is derived using the z-vector derivative scheme and is thus quite efficient, as it requires calculation of the molecular orbital (MO) density perturbation response to first order only. This exhibits the same overall cost scaling as self consistent field (SCF) second derivatives and requiring only minimal additional costs.
II. METHODS
AP approaches approximate the broken symmetry solutions of the Kohn-Sham or Hartree-Fock equations to be eigenfunctions of the Heisenberg spin Hamiltonian, which resolves the energy of the pure spin states contaminating the broken symmetry solution. In the limiting case, where there is a single contaminating spin state, which represents the current level of development for these methods, this requires two SCF calculations—one of the broken symmetry state and one of the higher energy contaminating state. Extending this beyond two states would require an extra SCF calculation for each additional contaminating state. These methods have the advantage over formal analytical projection methods, which involve factorial cost increase with the number of projected states. The AP method removes the (s + 1)th contaminant from the spin state of interest, referred to as the low spin (LS) state, using the energy expression EAP = αELS + (1 − α)EHS,
(1)
where ELS is the LS energy, EHS is the energy of the contaminating state, referred to as the high spin (HS) state, and α is 142, 054106-1
© 2015 AIP Publishing LLC
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-2
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
obtained from
2 x ⟨SLS ⟩ =
α= 2 ⟨SHS ⟩
2 ⟨SHS ⟩ − Sz,LS(Sz,LS + 1) 2 2 ⟨SHS ⟩ − ⟨SLS ⟩
µν
,
(2)
where and are the total electron spin squared expectation values of the HS and LS states, respectively, and Sz,LS is the z component of the total electron spin of the LS state.
+ −
−
−
as x x x , + (1 − α)EHS = α E˜LS EAP (E − E ) LS HS x x ⟨S 2 ⟩ x , E˜LS = ELS + 2 2 ⟨SHS⟩ − ⟨SLS ⟩ LS
( )x 2 x ⟨SLS ⟩ = Sz,LS(Sz,LS + 1) + n β − |Si j¯|2 ,
(5)
2 x ⟨SLS ⟩
µν
x Ω µν Sµν ,
(7)
where the intermediate matrices Υ, Λ, and Ω are given by Λ µν = Sµλ P¯λσ Sσν , (8) Υµν =
λσ
) x x Cν j¯ , Cµi Sµν Cν j¯ + Cµi Sµν Cνxj¯ + Cµi Sµν
(9)
Pµλ Sλσ P¯σν ,
(10)
λσ
where the over-bar indicates the β density matrix and its absence indicates the α density matrix. Expression (4) can now be re-written in a standard postSCF manner, which is amenable to z-vector substitution, x AP AP x E˜LS = Γµνλσ (µν|λσ) x + Pµν H µν µν
µνλσ
+
AP x W µν Sµν + VNx N ,
(11)
µν
where P AP is the sum of the SCF density matrix P and the AP correction density matrix P∆. Γ AP is given by AP ∆ ∆ Γµνλσ = Pµν Pλσ + Pµσ Pλν ,
(12)
AP ∆ ¯ ∆ ¯ Γ¯ µνλσ = P¯µν Pλσ + P¯µσ Pλν .
(13) AP
The energy weighted density matrix W is obtained AP from the sum of SCF and AP contributions, W µν = W µν ∆ AP − Ω µν + W µν . The term Ω is summed into W µν . This term is not common to post-SCF gradient theorems and is effectively a non-separable term arising from derivation in the AO basis. The AP energy weighted density matrix correction is most conventionally defined in the MO basis ∆ (p j||qi), (14) Wi∆j = − Λ µν Cµi Cν j − Ppq µν
pq
(6) where Sµν denotes the AO overlap matrix elements and C is the matrix of MO coefficients. We use the conventional notation— greek subscripts refer to AOs; roman subscripts refer to MOs; p, q, r, s . . . refer to all MOs; i, j, k, l . . . indicate occupied MOs; and a, b, c, d . . . indicate virtual MOs.26 Solving Eq. (6) requires the first order response of the density to the external perturbation(s), however, solution of the first-order coupled perturbed Hartree-Fock (CPHF) equations can be avoided by reformulating it as AO density matrix derivatives P x , back-transforming to the MO basis, and separating summations into occupied-occupied (oo), virtualvirtual (vv), occupied-virtual (ov), and virtual-occupied (vo) blocks. Then, using identities from CPHF theory (Pixj = −Sixj x x 27 and Pai = Pia ), we obtain
Sµσ Pσλ Sλν ,
λσ
Ω µν = 2
where n β is the number of β electrons in the LS system and Si j¯ is the α– β overlap matrix of the LS system, where the over-bar indicates a β spin orbital. Given that the first two terms in Eq. (5) are constants, only the last term requires special care. To obtain the α– β overlap derivative, we formulate it in the atomic orbital (AO) basis, giving
i j¯
Υµν (Cµ i¯Cν a¯ + Cµ a¯ Cν i¯)Pax¯ i¯
i¯a¯
(4)
i j¯
(
ia
µν
(3)
where the superscript denotes differentiation with respect to an external perturbation x. Note that it is assumed that the 2 HS state is spin pure, and ⟨SHS ⟩ is unchanged with geometric distortion. 2 x Thus, ⟨SLS ⟩ needs to be solved and depends upon the first order density response to the perturbation. While the value of ⟨S 2⟩ for DFT requires the two-particle density matrix, however, it has been shown that it is possible to approximate it using the same approach as Hartree-Fock (HF),
Si j¯
x Λ µν (Cµi Cνa + CµaCνi )Pai
µν
Following our previous work, the gradient can be written
Υµν Cµ i¯Cν j¯Si¯xj¯
i¯ j¯
µν
A. AP first derivatives
= −2
ij
µν
2 ⟨SLS ⟩
Λ µν Cµi Cν j Sixj
Wi¯∆j¯
=−
Υµν Cµ i¯Cν j¯ −
µν
¯ q¯i), ¯ Pp∆¯ q¯ ( p¯ j||
(15)
p¯ q¯
∆ ∆ Wai = −Pai ϵ i,
(16)
−Pa∆¯ i¯ϵ i¯,
(17)
Wa∆¯ i¯ ∆ Wab Wa∆¯ b¯
=
= 0,
(18)
= 0,
(19)
and once all blocks of W∆ are found, it is rotated to the AO basis ∆ ∆ W µν = Cµ pCνqW pq + Cµ p¯ Cν q¯ W p∆¯ q¯ . (20) pq
Pi∆j ,
p¯ q¯
Pi∆¯ j¯,
∆ Pab ,
Finally, the and Pa∆¯ b¯ terms are zero and the remaining ov and vo terms are obtained by invoking the
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-3
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
z-vector approach25 ∆ = L ai , (21) [(i j||ab) − (ib|| ja)]Pb∆j + (ϵ a − ϵ i )Pai
Cqx µ according to Cqx µ =
bj
u xpqCp µ .
(30)
p
¯ − (i¯b|| ¯ j¯a)]P ¯ a¯ b) [(i¯ j|| ¯ b∆¯ j¯ + (ϵ a¯ − ϵ i¯)Pa∆¯ i¯ = L a¯ i¯, (22) b¯ j¯
for Eqs. (21) and (22), ϵ p is the energy of MO p and the AP Lagrangian matrices L are L ai = − Λ µν (Cµi Cνa + CµaCνi ), (23)
The u x vectors are obtained as solutions to the first-order CPHF equations and used to derive the MO density derivative matrix with the relations Pixj = uix†j + uixj = −Sixj , x Pia x Pab
µν
L a¯ i¯ = −
Υµν (Cµ i¯Cν a¯ + Cµ a¯ Cν i¯).
(24)
µν
=
x uai ,
(32)
= 0.
(33)
The AP derivative density matrix can then be obtained from AP x x ∆x Pµν = Pµν + Pµν .
B. AP second derivatives
The second derivative term can be formulated in a similar manner to Eqs. (1) and (3) as xy EAP
xy xy = α E˜LS + (1 − α)EHS ,
(25)
where E˜ x y is given by xy xy E˜LS = E˜LS +
x x ( E˜LS − EHS ) 2 y ⟨S ⟩ , 2 2 ⟨SHS⟩ − ⟨SLS⟩ LS
Similarly, the AP effective two-particle density matrix derivative is given by AP x ∆x ∆ x ∆x ∆ x Γµνλσ = Pµν Pλσ + Pµν Pλσ − Pµσ Pλµ − Pµσ Pλν ,
(35)
∆ ¯x ∆x ¯ ∆ ¯x AP x ∆x ¯ Pλν , Pλµ − P¯µσ Pλσ − P¯µσ Pλσ + P¯µν Γ¯ µνλσ = P¯µν
(36)
xy xy E˜LS = ELS +
AP x x ∆x W µν = W µν − Ω xµν + W µν .
(26)
In Eq. (37), Ω is given by x x x ¯ Ω xµν = 2 Pµλ Sλσ P¯σν + Pµλ Sλσ P¯σν + Pµλ Sλσ Pσν ,
y y (ELS − EHS )
(ELS − EHS) 2 x y + 2 ⟨S ⟩ . 2 ⟨SHS⟩ − ⟨SLS ⟩ LS
and the components of W∆x are Wi∆x Λ xµν Cµi Cν j j =− µν
(ELS − EHS) 2 x y + 2 ⟨S ⟩ . 2 ⟨SHS⟩ − ⟨SLS ⟩ LS
+ +
µνλσ
APy x Pµν H µν +
µν
µν
AP x y H µν + Pµν
µν AP x y W µν Sµν
+ VNx yN .
−
µν
APy x W µν Sµν
µν
(29)
In addition to terms defined in the AP first derivatives, we now require the terms ΓAPy, PAPy, WAPy, as well as the one-electron, two-electron, and overlap integral second derivatives. These terms involve the vectors u xpq that describe the rotation of unperturbed orbitals Cp µ to perturbed orbitals
Λ µν (u xpi Cµ pCν j + Cµi u xp j Cν p )
p
∆x ∆ Ppq (ip|| jq) + Ppq (ip|| jq) x ,
(39)
pq
Wi¯∆x =− j¯
x Υµν Cµ i¯Cν j¯
µν
−
−
µν
(28)
In order to avoid the need to compute the second-order response of the MO density to a perturbation via the secondorder CPHF equations, Eq. (28) can be cast in a post-SCF second derivative form and z-vector substitution carried out, APy x y′ AP (µν|λσ) x + Γµνλσ (µν|λσ) x y E˜LS = Γµνλσ
− (27)
The first term in Eq. (27) is the second derivative of the spin-contaminated LS state, and the final term is dependent on the second derivative of the LS spin-squared expectation value
µνλσ
(38)
λσ
⟨S 2 ⟩ x 2 2 ⟨SHS ⟩ − ⟨SLS ⟩ LS (ELS − EHS) 2 y ⟩ + ⟨S 2 ⟩ x ⟨SLS 2 2 2 LS (⟨SHS ⟩ − ⟨SLS ⟩)
=
(37)
x
xy x where E˜LS is determined by Eq. (4) and E˜LS is found from direct differentiation of Eq. (4) to be
xy ELS
(34)
and the energy weighted density matrix derivative is
LS
x y′ E˜LS
(31)
Υµν (u xp¯ i¯Cµ p¯ Cν j¯ + Cµ i¯u xp¯ j¯Cν p¯ )
p¯
¯ ¯ j¯q) Pp∆x ¯ + Pp∆¯ q¯ (i¯p|| ¯ j¯q) ¯ x, ¯ q¯ (i p||
(40)
p¯ q¯ ∆x ∆x ∆ Wai = −Pai ϵ i − Pai (Fiix − Siix ϵ i ),
Wa∆x ¯ i¯
=
−Pa∆x ϵ − ¯ i¯ i¯
Pa∆¯ i¯(Fi¯xi¯ −
Si¯xi¯ϵ i¯),
(41) (42)
where we have used the relation ϵ ix = Fiix − ϵ i Siix . The derivative matrices Λ x and Υ x are given by x ¯ x x ¯x Λ xµν = Sµλ Pλσ Sσν + Sµλ P¯λσ Sσν + Sµλ Pλσ Sσν , (43) x Υµν
=
λσ
x x x x Sµλ Pλσ Sσν + Sµλ Pλσ Sσν + Sµλ Pλσ Sσν .
(44)
λσ ∆x Finally, the Pµν terms are determined by solving the set of response-like equations, ∆x [(i j||ab) − (ib|| ja)]Pb∆xj + (ϵ a − ϵ i )Pai = L ∆x (45) ai , bj
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-4
L. M. Thompson and H. P. Hratchian
x L ∆x ai = L ai −
J. Chem. Phys. 142, 054106 (2015)
x ∆ − Fixj )Pai , [(i j||ab) x − (ib|| ja) x ]Pb∆j − (Fab bj
(46)
¯ − ¯ a¯ b) [(i¯ j||
¯ j¯a)]P (i¯b|| ¯ b∆x ¯ j¯
+ (ϵ a¯ −
ϵ i¯)Pa∆x ¯ i¯
=
L ∆x , a¯ i¯
(47)
b¯ j¯
L ∆x = L ax¯ i¯ − a¯ i¯
¯ x − (i¯b|| ¯ j¯a) [(i j|| a¯ b) ¯ x ]Pb∆¯ j¯ − (Fa¯xb¯ − Fi¯xj¯)Pa∆¯ i¯, b¯ j¯
(48) where the Lagrangian derivative matrices are x L ai =− Λ xµν (Cµi Cνa + CµaCνi ) −
µν
Λ µν (u xpi Cµ pCνa
L ax¯ i¯
=−
+ Cµi u xpaCν p
+ Cµau xpi Cν p ),
x Υµν (Cµ i¯Cν a¯
(49)
+ Cµ a¯ Cν i¯)
µν
−
µν
Υµν (u xp¯ i¯Cµ p¯ Cν a¯ + Cµ i¯u xp¯ a¯ Cν p¯
p¯
+ u xp¯ a¯ Cµ p¯ Cν i¯ + Cµ a¯ u xp¯ i¯Cν p¯ ).
x L ai ) ( x† t t = − u C Λ C + C Λ C µ p µν νq µq µν ν p ∆x pi Wab p µν µν ( ) t t − u xt pa C µq Λ µν Cν p + C µ p Λ µν Cνq ,
(51)
p
+ u xpaCµ pCνi
Wi∆x j L x ia
p
µν
with u xpi at O(MON 2) cost. Focusing only on the second term of each equation, it is seen that the supermatrix (for the α spin case) is
(50)
Thus, by utilization of z-vector substitution, we only require the first-order response of the MO derivative coefficients to the perturbation through a set of CPHF equations in order to obtain AP second derivatives. This results in significant cost savings compared to the explicit calculation of the secondorder MO density response. C. Efficiency considerations for AP second derivatives
We note that implementation of Eqs. (25)–(50) as written above yields significant computation bottlenecks. As a result, the cost of an AP force constant calculation can require substantially more time than the combined time necessary for the constituent HS and LS frequency calculations. However, employing straightforward factorization techniques provides important cost improvements (vide infra). In what follows, we define M, N, O, and V to be the number of perturbations, AO basis functions, occupied MOs, and virtual MOs, respectively. As written, the second term of Eqs. (39), (40), (49), and (50) involves O(MON 4) work. This is because these equations represent the oo block (Eqs. (39) and (40)) and the ov block (Eqs. (49) and (50)) of two supermatrices, one for α and one for β spins. Due to the symmetry of these matrices and that the vv block is zero, we only need to solve for the oo and one of the ov blocks, leading to the O(ON) scaling on top of the O(M N 3) for matrix multiplications. This step can be reduced to a single O(N 3) step followed by individual O(ON 2) steps for each perturbation M. Furthermore, the number of matrix multiplications with u x terms can be reduced by recognizing the symmetry of Λ and Υ. The remaining terms are independent of perturbation and correspond to transformation of Λ and Υ into the MO basis. Thus, these matrices only require a single evaluation at O(N 3) cost, before looping over each perturbation and contracting
where, as described above, the symmetry of the matrices involved requires only explicit formation of the first summa∆x 1 x tion over p in Eq. (51) to obtain Wi∆x j , Wab , and /2 L ia . Second, using techniques derived in the solution of other post-SCF derivations,28,29 we avoid the need for the second interchange suggested by Eqs. (45) and (47). Although it is ∆x possible to directly form Pai , inspection of the terms with ∆x which Pai is contracted shows that they form the right-side of y the CPHF system of equations (commonly denoted as Q ai in CPHF theory), ( y y ∆x Pai hai − ϵ i Sai − Skyl (al||ik) kl
+
) ∆x y Q ai . Cµi Cνa Pλσ (µλ||νσ) y = Pai
(52)
µνλσ
The terms of Eqs. (45) and (47) can be regrouped to form the left-hand side of the CPHF equations (commonly denoted as A in CPHF theory) ) ( L ∆y ( ja||ab) − (ib|| ja) ai ∆x = , (53) Pai 1+ ϵa − ϵi ϵa − ϵi ∆x Pai
L ∆y ai = [1 + A]−1, ϵa − ϵi
(54)
y such that by multiplying through by Q ai , we are able to x generate an expression in terms of uai , ∆x y Pai Q ai
=
L ∆y ai [1
+ A]
−1
Q by j ϵb − ϵ j
∆x y x Pai Q ai = L ∆y ai uai ,
,
(55) (56)
x where L ∆y ai and uai have already been calculated prior to ∆x solving for Pai and can be contracted directly. This technique ∆ is the reverse of the interchange invoked to obtain Pai in the gradient derivation and avoids solution of an additional set of response equations in the calculation of force constants at the ∆x expense of never explicitly forming Pai .
III. NUMERICAL TESTS A. Computational scaling and timing
To illustrate the efficiency of our AP second derivative theory, this section describes numerical tests of computational scaling of the method with respect to system size. These calculations have been carried out on the diradical intermediate in the [2 + 2] cycloaddition of singlet oxygen with a terminal alkene chain (Fig. 1). The derivation
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-5
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
FIG. 1. Diradical intermediate in the [2 + 2] cycleaddition of singlet oxygen to a terminal alkene chain, with the α and β spin localization indicated by red and blue halos.
described above was implemented in a local development version of Gaussian30 and calculations were carried out using a production cluster equipped with dual 6-core 2.60 GHz Intel Xeon processors. Except where noted, each calculation was carried out in serial mode (using a single core) using up to 32 GB RAM. In order for our derivation to be useful, it is important that scaling of the method with respect to system size is kept to mean-field order and to minimize overhead relative to the sum of the two constituent SCF calculations. The minimum cost of an AP second derivative calculation is found in the limiting case that α is invariant to the perturbation(s). This xy xy corresponds to Eq. (25) in which E˜LS = ELS and so only SCF frequency calculations at both LS and HS multiplicities are required. In these calculations, the main expense is the solution of the first-order CPHF equations which grows as O(M 2OV ) in a direct formalism.31,32 Calculating the change in α with respect to the perturbation(s) involves computing the additional post-SCF terms, which increases the cost beyond that of the sum of the LS and HS calculations. These extra terms include the assembly of the Lagrangian derivative and calculation of the response density. Additionally, hybridDFT has additional cost beyond HF due to the evaluation of two electron integral derivatives contracted with the density over grid points.
FIG. 2. Comparison of elapsed time for AP second derivative calculation (red triangles) against LS + HS frequency calculations (blue squares) for varying number of basis functions. The number of perturbations is constant over all points. Calculations were carried out using B3LYP in serial mode on a production computer equipped with dual 6 core 2.60 GHz Intel Xeon processors.
To demonstrate the scaling of our program in terms of the number of perturbations, M, and the number of basis functions, N, we report two different sets of calculations. First, six different basis sets have been employed to calculate AP force constants for the compound shown in Fig. 1. This group of basis sets used yields a range of basis set size from 66 to 590, corresponding to a size range of the OV space in the LS sub-calculation from 14 400 to 59 904. Second, the 6-311G(d,p) basis set is used with a series of compounds based on that shown in Fig. 1. Each member of this set of compounds extends the alkane chain of the parent molecule by two carbon centers. The resulting group of molecules ranges from 26 to 68 atomic centers, corresponding to a range of 81–214 nuclear displacement perturbations included in force constant calculations. Calculation times for both sets of jobs using HF and B3LYP are reported in Table I. Figure 2 shows the scaling of AP second derivatives with number of basis functions for the molecule shown in Fig. 1 (red triangles) versus the scaling of the component SCF second derivative calculations (blue squares). This data are compiled from the top part of Table I, which keeps the
TABLE I. Parameters and elapsed times in seconds for frequency calculations using AP and the sum of the two constituent SCF second derivative calculations. Values are taken from calculations carried out using a single core. HS + LS SCF freq. Natoms
AP freq.
Basis set
Nbasis
M
OVLS
HF
B3LYP
HF
B3LYP
26
STO-3G 6-31G 6-311G 6-311G(d,p) 6-311++G (d,p) cc-pVTZ
66 122 178 286 342 590
81 81 81 81 81 81
1040 3280 5520 9840 12 080 22 000
11.6 61.2 218.9 1468.9 2856.6 35 972.1
306.3 756.6 1405.3 3787.3 7553.7 39 599.4
12.8 68.0 243.4 1708.7 3276.4 39 453.9
423.5 1014.5 1880.4 5101.3 9256.9 45 518.4
32 38 44 50 56 62 68
6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p)
348 410 472 534 596 658 720
96 114 132 150 168 196 214
14 400 19 824 26 112 33 264 41 280 50 160 59 904
3514.5 7846.6 13 101.8 19 710.2 28 203.7 36 050.7 46 398.8
6799.5 11 811.2 17 429.2 24 721.3 33 488.9 42 496.4 54 897.4
4126.7 8995.6 14 909.5 22 397.2 32 191.0 41 570.5 53 290.9
8870.6 14 682.7 21 834.8 31 095.7 42 346.1 50 903.9 63 721.2
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-6
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
TABLE II. Parameters and elapsed times in seconds for AP and SCF of sample calculations in parallel mode using a production computer equipped with dual 6-core 2.60 GHz Intel Xeon processors. HS + LS SCF freq. Natoms
26
AP freq.
Basis set
Nbasis
M
OVLS
HF
B3LYP
HF
B3LYP
STO-3G 6-31G 6-311G 6-311 + +G (d,p) cc-pVTZ
66 122 178 342 590
81 81 81 81 81
1040 3280 5520 12 080 22 000
10.1 16.4 39.5 390.3 4720.4
37.0 85.4 151.7 793.9 4859.9
10.6 21.8 47.8 704.7 5252.9
48.6 108.1 198.6 1078.0 5614.1
number of perturbations constant. From this figure, it can be seen that the overhead of the AP second derivatives is small compared to the underlying SCF calculations. In fact, the overhead as a fraction of the SCF second derivative calculations, which ranges from 10% to 40%, decreases as the OV space increases. This is due to the scaling of the CPHF solutions with the number of virtual orbitals, which grows with the size of basis. As this is common to both sets of calculations, its expense dominates the calculation and the relative cost of the post-SCF like terms rapidly decreases. For example, with 178 basis functions, solution of the CPHF equations in the B3LYP calculation accounts for 27% of job time. In contrast, at 590 basis functions solution of the CPHF work accounts for 65% of total CPU time. In the lower part of Table I, we show the scaling for calculations on the molecule in Fig. 1 using the 6-311G(d,p) basis set as the number of carbon atoms in the alkyl chain is increased by two carbons from 6 to 20. Therefore, this indicates how our method scales with regard to an increase in the number of perturbations, as well as the number of basis functions and OV size. This set of results indicates that the overhead of the AP method does not depend on the number of perturbations, giving the same trend as the results in which only the number of basis functions were increased. In all sets of results, the HF calculations were observed to have smaller overhead than the corresponding B3LYP calculations. This is because the formation of the Lagrangian derivative requires the contraction of the integral derivatives with the response density evaluated over grid points. In order to demonstrate that our AP second-derivative implementation performs well in parallel computing environments, Table II shows job times for sample calculations carried out using 12 processors. The total AP wall time speedup going from 1 to 12 processors ranges from 68% to 79% efficiency for B3LYP, which is comparable to standard DFT frequency calculations. The efficiency is somewhat lower for HF as this does not benefit from the highly parallelized exchange-correlation integration over grid points. Despite this, HF also gives valuable cost reduction when run in parallel mode. This indicates that the post-SCF components of the AP second derivative calculation parallelize as well as the components from SCF calculations. Figure 3 shows the same trend as the serial calculations for the overhead with regard to system size, indicating that in this case, it is also the solution of the CPHF equations which dominate at large N and OV .
B. Example calculations
In order to demonstrate how the application of AP second derivatives modifies the calculated force constants, we have run several calculations on potentially challenging radical triatomic systems (Table III).33 These results show that the broken symmetry solutions of LiO2 and NO2 have small amounts of spin contamination. For this reason, application of the AP method results in very small changes in the calculated vibrational frequencies. This is satisfying as both the broken symmetry and AP solutions give respectable agreement with Coupled Cluster Singles & Doubles (CCSD) calculations using the same basis set. However, the broken 2 ⟩ symmetry C+3 is more significantly contaminated, with ⟨SLS equal to 1.022. This causes UB3LYP to give shorter bond lengths than the equivalent CCSD calculation, as well as a significant deviation in the ω1(a1) vibrational mode. This mode corresponds to the bending mode, and application of AP causes a shift to higher frequency of 116 cm−1. This suggests that the gap between the contaminating state and the ground state varies along this coordinate and is close enough that this variance has an effect on the curvature of the broken symmetry solution. In the case of the other modes, either the change in energy gap is small, or the gap is too large to have any effect. Further investigation into the effects of second derivatives on vibrational modes and spectra is ongoing in our lab and will be reported in due course.
FIG. 3. Comparison of elapsed time for AP second derivative calculation (red triangles) against LS + HS frequency calculations (blue squares) for varying number of basis functions. The number of perturbations is constant over all points. Calculations were carried out in parallel mode using 12 cores on production computer equipped with dual 6-core 2.60 GHz Intel Xeon processors.
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-7
L. M. Thompson and H. P. Hratchian
J. Chem. Phys. 142, 054106 (2015)
TABLE III. Geometric parameters and vibrational modes of example triatomic molecules calculated using UB3LYP and AP-B3LYP with cc-pVTZ compared to CCSD. 2A 2
LiO2
Method
RLiO
θ OLiO
2 ⟩ ⟨SLS
UB3LYP AP-B3LYP CCSD
1.747 1.746 1.757
45.2 45.3 44.6
0.757 ... ... 2B 2
ω 3 (a 1)
535 537 536
761 762 763
1170 1167 1204
C+3
RCC
θ CCC
2 ⟩ ⟨SLS
ω 1 (a 1)
ω 2 (b2)
ω 3 (a 1)
UB3LYP AP-B3LYP CCSD
1.309 1.315 1.317
66.6 65.2 65.4
1.022 ... ...
742 958 877
1238 1302 1297
1660 1692 1666
1
NO2
Method
RNO
θ ONO
2 ⟩ ⟨SLS
ω 1 (a 1)
ω 2 (a 1)
ω 3 (b 2)
UB3LYP AP-B3LYP CCSD
1.190 1.190 1.183
134.4 134.4 134.8
0.754 ... ...
769 769 786
1394 1394 1434
1696 1694 1763
IV. CONCLUSION
By casting the equations into the form of a postSCF model, we have been able to derive analytical second derivatives of the AP method. Application of interchange and reverse-interchange techniques has enabled implementation of an efficient algorithm that we have demonstrated has a small amount of overhead relative to two SCF frequency calculations. These second derivatives include the orbital relaxation terms that account for the change in the alphabeta overlap matrix under the perturbations considered and should be significantly more accurate than derivations which ignore this term. The availability of AP second derivatives will allow the study of spin pure states in a number of systems for which spin contamination adversely affects results and for which analytical second derivatives are desired. ACKNOWLEDGMENTS
We gratefully acknowledge Dr. Michael Frisch (Gaussian, Inc.), Dr. Giovanni Scalmani (Gaussian, Inc.), and Professor Marco Caricato (University of Kansas) for helpful discussions and the University of California, Merced for financial support. Burke, J. Chem. Phys. 136, 150901 (2012). D. Becke, J. Chem. Phys. 140, 18A301 (2014). 3J. L. Sonnenberg, H. B. Schlegel, and H. P. Hratchian, in Computational Inorganic and Bioinorganic Chemistry, edited by E. I. Solomon, R. A. Scott, and B. A. King (John Wiley & Sons, Ltd., 2009). 4F. Neese, Coord. Chem. Rev. 253, 526 (2009). 5S. E. Waller, J. E. Mann, D. W. Rothgeb, and C. C. Jarrold, J. Phys. Chem. A 116, 9639 (2012). 6D. Aradilla and C. Aleman, J. Phys. Org. Chem. 27, 867 (2014). 7P. Jost and C. van Wüllen, Phys. Chem. Chem. Phys. 15, 16426 (2013). 8J. A. Pople, P. M. W. Gill, and N. C. Handy, Int. J. Quantum Chem. 56, 303 (1995). 9U. Bozkaya, J. Chem. Theory Comput. 10, 4389 (2014). 10U. Bozkaya, J. Chem. Phys. 135, 224103 (2011). 2A.
ω 2 (a 1)
Method
2A
1K.
ω 1 (b2)
11H.
B. Schlegel and M. Frisch, in Theoretical and Computational Models for Organic Chemistry, edited by S. J. Formosinho, I. Csizmadia, and L. Arnaut (Kluwer Academic Publishers, 1991). 12P. O. Löwdin, Phys. Rev. 97, 1509 (1955). 13C. Sosa and H. B. Schlegel, Int. J. Quantum Chem. 32, 267 (1987). 14N. C. Handy and P. J. Knowles, J. Chem. Phys. 88, 6991 (1988). 15N. C. Handy and P. J. Knowles, J. Phys. Chem. 92, 3097 (1988). 16L. M. Thompson and H. P. Hratchian, J. Chem. Phys. 141, 034108 (2014). 17L. Noodleman, C. Y. Peng, D. A. Case, and J. M. Mouesca, Coord. Chem. Rev. 144, 199 (1995). 18K. Yamaguchi, F. Jensen, A. Dorigo, and K. N. Houk, Chem. Phys. Lett. 149, 537 (1988). 19K. Yamaguchi, M. Okumura, W. Mori, and J. Maki, Chem. Phys. Lett. 210, 201 (1993). 20S. Yamanaka, M. Okumura, M. Nakano, and K. Yamaguchi, J. Mol. Struct.: THEOCHEM 310, 205 (1994). 21H. P. Hratchian, J. Chem. Phys. 138, 101101 (2013). 22Y. Kitagawa, T. Saito, M. Ito, M. Shoji, K. Koizumi, S. Yamanaka, T. Kawakami, M. Okumura, and K. Yamaguchi, Chem. Phys. Lett. 442, 445 (2007). 23T. Saito and W. Thiel, J. Phys. Chem. A 116, 10864 (2012). 24T. Saito, S. Nishihara, Y. Kataoka, Y. Nakanishi, T. Matsui, Y. Kitagawa, T. Kawakami, M. Okumura, and K. Yamaguchi, Chem. Phys. Lett. 483, 168 (2009). 25N. C. Handy and H. F. Schaefer, J. Chem. Phys. 81, 5031 (1984). 26J. A. Pople and D. L. Beveridge, Approximate Molecular Orbital Theory (McGraw-Hill, Inc., 1970). 27K. Raghavachari, J. A. Pople, and H. B. Schlegel, Int. J. Quantum Chem. 16, 225 (1979). 28R. Cammi, B. Mennucci, C. Pomelli, C. Cappelli, S. Corni, L. Frediani, G. W. Trucks, and M. Frisch, Theor. Chem. Acc. 111, 66 (2004). 29G. W. Trucks, M. Frisch, M. Head-Gordon, H. B. Schlegel, E. A. Salter, and J. L. Andres, “An efficient theory and implementation of MP2 second derivatives” (unpublished). 30M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, B. Janesko, F. Lipparini, G. Zheng, J. L. Sonnenberg, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin,
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
054106-8
L. M. Thompson and H. P. Hratchian
K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, P. V. Parandekar, N. J. Mayhall, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, “Gaussian development version, Revision H.38+,” Gaussian, Inc., 2010.
J. Chem. Phys. 142, 054106 (2015) 31M. Head-Gordon, J. A. Pople, and M. Frisch, Chem. Phys. 141, 189 (1990). 32L.
M. Thompson, A. Lasoroski, P. M. Champion, J. T. Sage, M. Frisch, J. J. van Thor, and M. J. Bearpark, J. Chem. Theory Comput. 10, 751 (2014). 33U. Bozkaya and C. D. Sherrill, J. Chem. Phys. 138, 184103 (2013).
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: 128.6.218.72 On: Wed, 12 Aug 2015 06:56:27
