The relationship between the random-phase-approximation (RPA) correlation energy and the continuous algebraic Riccati equation is examined and the importance of a stabilizing solution is emphasized. The criterion to distinguish this from non-stabilizing solutions can be used to ensure that physical, smooth potential energy surfaces are obtained. An implementation of analytic RPA molecular gradients is presented using the Lagrangian technique. Illustrative calculations indicate that RPA with Hartree-Fock reference orbitals delivers an accuracy similar to that of second-order Møller–Plesset perturbation theory.

Theories based on the random-phase approximation (RPA) have received a great deal of recent attention, see Refs. 1–3 for recent reviews. In the context of hybrid density-functional theories, RPA methods offer the tantalizing prospect of a relatively low-cost non-local orbital-dependent correlation energy applicable to metallic and insulating systems alike. The connection of RPA methods to ring-coupled-cluster-doubles (rCCD) methods was recently clarified and inspired a number of related theories that attempt to include additional diagrams in the rCCD approach to improve the accuracy of the method further, while maintaining its attractive low cost.

The triplet instability problems4 that afflict full RPA are ameliorated by using direct (non-anti-symmetrized) integrals in the computation of the correlation energy. Two stable RPA variants have become particularly popular: direct-RPA (dRPA)5 and dRPA with second-order screened exchange (SOSEX).6 Here, we examine the relationship between the rCCD-based implementations of these approaches and their associated continuous algebraic Riccati equations (CAREs).7 We highlight the need for care in selecting the appropriate solution to these equations and demonstrate the use of a criterion that identifies the desired solution. We show how the Lagrangian approach to molecular properties8 may be used to derive analytic molecular gradients of the dRPA and SOSEX energies. The quality of dRPA and SOSEX geometries is assessed by comparison with accurate geometries for small covalently bound molecules and rare-gas (RG) dimers. Hartree-Fock (HF) reference orbitals are used throughout this work.

The RPA ground-state correlation energy is given by the plasmon formula9 

\begin{equation}E_{\rm{{c}}} = \frac{1}{2} \sum _{K=1}^N \big( \omega _K - \omega _K^{\rm{{TDA}}} \big),\end{equation}
Ec=12K=1NωKωK TDA ,
(1)

where N is the number of singly-excited determinants with respect to the HF reference state and ωK and

$\omega _K^{\rm{{TDA}}}$
ωK TDA are excitation energies within the RPA and the Tamm-Dancoff approximation (TDA), respectively. The neglect of correlated exchange in dRPA implies that only singlet states contribute to the correlation energy. We here consider a spin-adapted formulation of dRPA with a closed-shell HF reference state. The plasmon formula in Eq. (1) is inconvenient for computation of the energy and its derivatives, and instead we set up a variational Lagrangian for the total energy.8 

To establish notation, we recapitulate canonical HF theory for closed-shell ground states.4 Molecular orbitals (MOs) are determined from the canonical HF conditions

\begin{equation}F_{pq} = \delta _{pq} \varepsilon _p.\end{equation}
Fpq=δpqɛp.
(2)

Here, ɛp is the orbital energy of the pth MO (indices p, q, r, s, t are used for general MOs, i, j for occupied MOs, and a, b for virtual MOs) and F is the Fock matrix with elements

\begin{equation}F_{pq} = h_{pq} + \sum _i [ 2(pq \vert ii) - (pi \vert iq) ] ,\end{equation}
Fpq=hpq+i[2(pq|ii)(pi|iq)],
(3)

where hpq and (pq|rs) are one- and two-electron (Mulliken notation) Hamiltonian integrals, respectively. The HF ground-state energy may then be written as

\begin{equation}E_{\rm{{HF}}} = 2\sum _i h_{ii} + \sum _{ij} [ 2(ii \vert jj) - (ij \vert ji) ].\end{equation}
E HF =2ihii+ij[2(ii|jj)(ij|ji)].
(4)

For a closed-shell HF reference state, the dRPA correlation energy is given by7 

\begin{equation}E_{\rm{{c}}} = \frac{1}{2}{\rm Tr}( \bm{ZB} ) = \frac{1}{2} \sum _{aibj} Z_{aibj}B_{bjai},\end{equation}
Ec=12 Tr (ZB)=12aibjZaibjBbjai,
(5)

where the amplitude matrix

${\bm Z}=\bm{YX}^{-1}$
Z=YX1 is defined in terms of the excitation (X) and de-excitation (Y) parts of the eigenvectors of the symplectic eigenvalue problem

\begin{equation}\left( \begin{array}{r@\quad r}{\bm A} & {\bm B} \\[6pt]{\bm B} & {\bm A} \end{array} \right) \left( \begin{array}{r}{\bm X} \\[6pt]{\bm Y} \end{array} \right) = \left( \begin{array}{rr}\bm{X\omega } \\[6pt]-\bm{Y\omega } \end{array} \right).\end{equation}
ABBAXY=XωYω.
(6)

Here, ω = diag(ω1, ω2, ⋅⋅⋅, ωN) are the dRPA singlet excitation energies and

\begin{equation}A_{aibj} = F_{ba}\delta _{ij} - F_{ij}\delta _{ab} + B_{aibj}, \quad B_{aibj} = 2(ai \vert bj).\end{equation}
Aaibj=FbaδijFijδab+Baibj,Baibj=2(ai|bj).
(7)

The positive definiteness of B guarantees the existence of

${\bm X}^{-1}$
X1 and the negative definiteness of Z.7 As shown in Ref. 7, Z may be obtained directly from the CARE

\begin{equation}{\bm R}({\bm Z}) = {\bm B}+\bm{AZ}+\bm{ZA}+\bm{ZBZ} = {\bm 0}.\end{equation}
R(Z)=B+AZ+ZA+ZBZ=0.
(8)

In SOSEX,6 exchange screening is reintroduced by calculating the energy from the dRPA amplitudes via

\begin{equation}E_{\rm{{c}}} = \frac{1}{2} \sum _{aibj} Z_{aibj} \bigl ( B_{bjai}-\textstyle{\frac{1}{2}}B_{biaj} \bigr ).\end{equation}
Ec=12aibjZaibjBbjai12Bbiaj.
(9)

The total dRPA [SOSEX] energy is defined as the sum of Eqs. (4) and (5) [(9)]. While SOSEX eliminates self-correlation for one-electron systems, the static correlation captured by dRPA is ruined.10 

To compute derivatives of dRPA/SOSEX energies, we must consider the variations in Eqs. (2) and (8) due to variations in the orbitals. For a set of orthonormal orbitals

$\lbrace \tilde{\phi } \rbrace$
{ϕ̃}⁠, Eq. (2) determines the orbital rotation parameters κpq = −κqp such that
$\phi _p = \sum _q \tilde{\phi }_q [ \exp ( -\bm{\kappa } ) ]_{qp}$
ϕp=qϕ̃q[exp(κ)]qp
. Using Eqs. (2) and (8) as constraints and introducing the Lagrange multipliers
$\bar{{\bm Z}}$
Z¯
and
$\bar{\bm{\kappa }}$
κ¯
, the dRPA/SOSEX Lagrangian becomes

\begin{align}L = &E_{\rm{{HF}}}(\bm{\kappa }) + E_{\rm{{c}}}(\bm{\kappa },{\bm Z}) \nonumber \\&+ \sum _{ai \ge bj} \bar{Z}_{aibj} R_{aibj}(\bm{\kappa },{\bm Z}) + \sum _{p>q} \bar{\kappa }_{pq} F_{pq}(\bm{\kappa }).\end{align}
L=E HF (κ)+Ec(κ,Z)+aibjZ¯aibjRaibj(κ,Z)+p>qκ¯pqFpq(κ).
(10)

To obtain the dRPA/SOSEX Lagrangian for a Kohn–Sham (KS) reference, simply replace the Fock matrix of the last term with the KS matrix. While Z and

$\bar{{\bm Z}}$
Z¯ are symmetric, κ and
$\bar{\bm{\kappa }}$
κ¯
are antisymmetric and symmetric, respectively, and
$\bar{\kappa }_{pp}=0$
κ¯pp=0
. The variational conditions on the multipliers give Eqs. (2) and (8). The variational condition on Z gives the Lyapunov equation for
$\bar{{\bm Z}}$
Z¯
,

\begin{equation}{\bm G}\bar{{\bm Z}} + \bar{{\bm Z}}{\bm G}^T = -{\bm B}, \qquad {\bm G}({\bm Z}) = {\bm A}+\bm{BZ},\end{equation}
GZ¯+Z¯GT=B,G(Z)=A+BZ,
(11)

which is solved by the conjugate-gradient method.11,12 The variational condition on κ gives an equation for

$\bar{\bm{\kappa }}$
κ¯⁠,

\begin{equation}\bar{\bm{\kappa }} \bm{\cal {A}} = -\bm{\eta },\end{equation}
κ¯A=η,
(12)

where

$\bm{\cal {A}}$
A depends on orbital energies and two-electron integrals, while η depends also on Z and
$\bar{{\bm Z}}$
Z¯
. To construct the Lagrangian, we compute κ from Eq. (2) and then Z from Eq. (8); next, we obtain
$\bar{{\bm Z}}$
Z¯
from Eq. (11) and finally
$\bar{\bm{\kappa }}$
κ¯
from Eq. (12). The same procedure applies to SOSEX, with modified right-hand sides of Eqs. (11) and (12).

The variational Lagrangian Eq. (10) may be recast as

\begin{equation}L = \sum _{pq} D_{pq}h_{pq} + \frac{1}{2} \sum _{pqrs} d_{pqrs}(pq \vert rs),\end{equation}
L=pqDpqhpq+12pqrsdpqrs(pq|rs),
(13)

where D and d are the variational one- and two-electron density matrices, chosen to be permutationally symmetric. We may then compute forces in the usual manner,8 

\begin{eqnarray}L^{(1)} &=& \sum _{\mu \nu } D_{\mu \nu } h^{(1)}_{\mu \nu } + \frac{1}{2} \sum _{\mu \nu \kappa \lambda } d_{\mu \nu \kappa \lambda } (\mu \nu \vert \kappa \lambda )^{(1)} \nonumber \\&&- \sum _{\mu \nu } F^{\text{{eff}}}_{\mu \nu } S^{(1)}_{\mu \nu },\end{eqnarray}
L(1)=μνDμνhμν(1)+12μνκλdμνκλ(μν|κλ)(1)μνFμνeffSμν(1),
(14)

where superscript (1) denotes the first derivative with respect to a nuclear coordinate at the reference geometry, S is the atomic-orbital (AO) overlap matrix, and Greek indices denote AOs. The AO density matrices are given by Dμν = ∑pqCμpDpqCνq and dμνκλ = ∑pqrsCμpCνqdpqrsCκrCλs where C contains the canonical MO coefficients at the reference geometry. As in CC gradient theory,13 the reorthonormalization term is computed through a generalized effective Fock matrix backtransformed to AO basis,

\begin{equation}F^{\text{{eff}}}_{pq} = \sum _t D_{pt}h_{tq} + \frac{1}{2} \sum _{trs} ( d_{ptrs}+d_{tprs} ) (tq \vert rs).\end{equation}
Fpqeff=tDpthtq+12trs(dptrs+dtprs)(tq|rs).
(15)

The solution of the CARE, Eq. (8), can be challenging.10 Only one of the multiple solutions corresponds to the ground state of the symplectic eigenvalue problem in Eq. (6). The (physical) desired solution is stabilizing, meaning that the nonsymmetric matrix

${\bm G}({\bm Z})$
G(Z) has only positive eigenvalues. A stabilizing solution to Eq. (8) is not necessarily unique, however.

Among the non-stabilizing solutions, one is of particular interest. Starting from the paired eigenvalue problem,

\begin{equation}\left( \begin{array}{r@\quad r}{\bm A} & {\bm B} \\[6pt]{\bm B} & {\bm A} \end{array} \right) \left( \begin{array}{r}{\bm Y} \\[6pt]{\bm X} \end{array} \right) = - \left( \begin{array}{rr}\bm{Y\omega } \\[6pt]-\bm{X\omega } \end{array} \right)\end{equation}
ABBAYX=YωXω
(16)

and repeating the derivation of Ref. 7, we obtain the CARE for the inverse amplitudes

$\bm{\tilde{Z}}=\bm{XY}^{-1}={\bm Z}^{-1}$
Z̃=XY1=Z1⁠. Both Z and
${\bm Z}^{-1}$
Z1
are guaranteed to exist in dRPA since B is positive definite.7 If (as assumed) Z is stabilizing, then
$\bm{\tilde{Z}}$
Z̃
is not stabilizing since the eigenvalues of
${\bm G}(\bm{\tilde{Z}}) = -\bm{Y\omega Y}^{-1}$
G(Z̃)=YωY1
are all negative. Use of the inverse solution to compute the correlation energy from Eq. (5) will not reproduce the plasmon formula Eq. (1). Instead, the correct energy may be obtained directly from the inverse solution according to
$E_{\rm{{c}}} = -\rm{Tr}(\bm{\tilde{Z}}{\bm B})/2 -\rm{Tr}{\bm A}$
Ec= Tr (Z̃B)/2 Tr A
, which may be used to set up an alternative Lagrangian for this particular non-stabilizing solution. In general, a range of non-stabilizing solutions, whose amplitudes are not the inverse of the physical solution, may be obtained and the positive definiteness of
${\bm G}({\bm Z})$
G(Z)
is the only useful criterion for identifying a stabilizing CARE solution.

To solve Eq. (8), we use an iterative Newton-like algorithm combined with direct inversion in the iterative subspace (DIIS)14 to accelerate convergence. The iterative procedure is defined by15 

\begin{equation}{\bm Z}^{(k+1)} = {\bm Z}^{(k)} + {\bm N}^{(k)},\end{equation}
Z(k+1)=Z(k)+N(k),
(17)

where

${\bm N}^{(k)}$
N(k) is an approximate solution of the Lyapunov equation
${\bm G}^T({\bm Z}^{(k)}){\bm N}^{(k)} + {\bm N}^{(k)}{\bm G}({\bm Z}^{(k)}) = -{\bm R}({\bm Z}^{(k)})$
GT(Z(k))N(k)+N(k)G(Z(k))=R(Z(k))
, obtained by assuming that
${\bm G}({\bm Z}^{(k)})$
G(Z(k))
is diagonal dominant,

\begin{equation}N^{(k)}_{aibj} = -\frac{R_{aibj}({\bm Z}^{(k)})}{G_{aiai}({\bm Z}^{(k)})+G_{bjbj}({\bm Z}^{(k)})}.\end{equation}
Naibj(k)=Raibj(Z(k))Gaiai(Z(k))+Gbjbj(Z(k)).
(18)

The amplitudes are then refined by means of DIIS, and the iterative procedure converges when the norm of the dRPA vector function is below a given tolerance. A unique stabilizing solution is not guaranteed but we may check stability by diagonalizing

${\bm G}({\bm Z})$
G(Z)⁠.

The initial trial vector

${\bm Z}^{(0)}$
Z(0) may be chosen in different ways. The standard choice in CC theory for insulators is the direct second-order Møller–Plesset (MP2) start guess

\begin{equation}Z^{(0)}_{aibj} = -\frac{B_{aibj}}{\varepsilon _a-\varepsilon _i+\varepsilon _b-\varepsilon _j},\end{equation}
Zaibj(0)=Baibjɛaɛi+ɛbɛj,
(19)

which, however, becomes numerically unstable with a small or vanishing gap between the occupied and virtual orbital energies. A better start guess10 is an approximate solution of the Lyapunov equation

$\bm{AZ}+\bm{ZA}=-{\bm B}$
AZ+ZA=B⁠, obtained by assuming that A is diagonal dominant,

\begin{equation}Z^{(0)}_{aibj} = -\frac{B_{aibj}}{A_{aiai}+A_{bjbj}}.\end{equation}
Zaibj(0)=BaibjAaiai+Abjbj.
(20)

Note that this guess is identical to the amplitudes of the first iteration when the zero vector is chosen as the initial trial vector. The zero vector has the advantage that

${\bm G}({\bm 0})={\bm A}$
G(0)=A is symmetric positive definite and therefore a formally valid guess for a stabilizing solution. This is not necessarily the case for the MP2 start guess of Eq. (19).

All calculations have been performed with a development version of DALTON16 using HF reference orbitals.

Potential-energy curves for H2 computed at the (restricted) HF, MP2, CCD, coupled-cluster singles and doubles (CCSD), dRPA, and SOSEX levels of theory in the aug-cc-pVQZ basis set are shown in Fig. 1. The stabilizing CARE solution is obtained in ten or less iterations at all bond distances r with Eq. (20) as initial guess. The same stabilizing solution is obtained for r < 8a0, when the MP2 start guess, Eq. (19), is used, albeit more iterations (up to 33) are required for r > 4 a0. For r > 8 a0, the MP2 start guess typically does not lead to a stabilizing solution. The dRPA and SOSEX energies corresponding to these non-stabilizing solutions are plotted as points in Fig. 1. These points form alternative continuous curves, suggesting that these solutions are present at all distances. An analysis reveals that they are not the inverse of the stabilizing solutions—that is, they do not correspond to the dRPA symplectic eigenvalue problem and the resulting energies are unphysical. The start guess in Eq. (20) is therefore used in this work.

FIG. 1.

Potential energy curves computed in the aug-cc-pVQZ basis set. dRPA(MP2) and SOSEX(MP2) refer to results obtained using the MP2 start guess, Eq. (19).

FIG. 1.

Potential energy curves computed in the aug-cc-pVQZ basis set. dRPA(MP2) and SOSEX(MP2) refer to results obtained using the MP2 start guess, Eq. (19).

Close modal

Figure 2 compares the results of analytic gradient-based geometry optimizations at the MP2, dRPA, and SOSEX levels of theory in the cc-pVQZ basis for 16 small molecules (NH3, CH2, F2, C2H4, H2, HCN, HNC, HNO, CO, CO2, CH4, HOF, N2H2, N2, CH2O, H2O) with accurate reference geometries. These were obtained by Pawłowski et al.17 from experimental rotational constants and highly accurate vibration-rotation interaction constants computed at the CCSD with perturbative triples (CCSD(T)) level in the cc-pVQZ basis. All bond lengths and non-trivial bond angles (i.e., those not determined by symmetry) are included in Fig. 2.

FIG. 2.

(a) Computed covalent bond lengths in pm and (b) bond angles in degrees versus semiempirical reference values.17 

FIG. 2.

(a) Computed covalent bond lengths in pm and (b) bond angles in degrees versus semiempirical reference values.17 

Close modal

The geometries obtained with the Coulomb-only dRPA correlation energy are generally more accurate than those obtained with SOSEX, which does not bring the results closer to the MP2 values. The mean absolute bond-distance error is 0.5 pm for MP2, 1.6 pm for dRPA, and 2.0 pm for SOSEX, while the mean absolute bond-angle error is 0.3° for MP2, 0.4° for dRPA, and 0.8° for SOSEX. The largest absolute bond-distance errors occur for the longest bonds (the F–F bond in F2 at 141.3 pm and the O–F bond in HOF at 143.4 pm), which contribute significantly to the mean relative errors of 0.4%, 1.4%, and 1.8% for MP2, dRPA, and SOSEX, respectively. The mean relative bond-angle error is less than 1% for all three methods. Although the errors are somewhat greater, dRPA and SOSEX geometries are comparable to those of MP2.

Tran and Hutter18 investigated the performance of some popular nonlocal dispersion corrections—namely, the van der Waals density functionals (vdW-DF19 and vdW-DF220) of Langreth, Lundqvist and co-workers, the D3 dispersion correction of Grimme,21 and Björkman's revised Vydrov-Van Voorhis functional (rVV10).22,23 Equilibrium distances and binding energies for He2, Ne2, and Ar2 computed at the MP2, dRPA, and SOSEX levels of theory are listed along with (selected) results of Tran and Hutter in Table I and compared with accurate reference values. The accuracy of dRPA/SOSEX is similar to that of MP2 for equilibrium distances and binding energies of the RG dimers, with SOSEX increasing the difference between dRPA and MP2. Equilibrium distance errors from dRPA and SOSEX are less than 5%, while binding energy errors range from 15% to almost 100%. None of the DFT methods in the table are consistently better than dRPA/SOSEX, and none can be considered highly accurate relative to CCSD(T), which agrees with the reference values within 0.5% for distances and 3% for energies of the RG dimers.25 

Table I.

Equilibrium distances (re, in Å) and binding energies (ΔE, in meV) for RG dimers in the aug-cc-pV5Z basis without counterpoise correction. DFT results from Ref. 18, reference results from Ref. 24.

 He2Ne2Ar2
 reΔEreΔEreΔE
vdW-DF 2.82 6.6 3.07 14.1 3.92 23.1 
vdW-DF2 2.75 2.8 2.95 9.2 3.75 18.3 
PBE-D3 2.66 5.7 3.01 9.9 3.88 15.3 
revPBE-D3 2.90 3.0 3.20 5.6 3.93 12.8 
B97D-D3 3.01 0.8 3.33 4.3 3.99 11.3 
rVV10 2.92 0.9 3.01 5.6 3.73 13.9 
MP2 3.07 0.61 3.15 3.12 3.63 27.43 
dRPA 3.13 0.46 3.14 3.10 3.72 21.96 
SOSEX 3.13 0.45 3.17 2.51 3.78 16.32 
Reference 2.97 0.95 3.09 3.65 3.76 12.35 
 He2Ne2Ar2
 reΔEreΔEreΔE
vdW-DF 2.82 6.6 3.07 14.1 3.92 23.1 
vdW-DF2 2.75 2.8 2.95 9.2 3.75 18.3 
PBE-D3 2.66 5.7 3.01 9.9 3.88 15.3 
revPBE-D3 2.90 3.0 3.20 5.6 3.93 12.8 
B97D-D3 3.01 0.8 3.33 4.3 3.99 11.3 
rVV10 2.92 0.9 3.01 5.6 3.73 13.9 
MP2 3.07 0.61 3.15 3.12 3.63 27.43 
dRPA 3.13 0.46 3.14 3.10 3.72 21.96 
SOSEX 3.13 0.45 3.17 2.51 3.78 16.32 
Reference 2.97 0.95 3.09 3.65 3.76 12.35 

We have presented an implementation of analytic molecular gradients using the Lagrangian approach for the dRPA and SOSEX methods. Our pilot implementation scales with system size M as O(M6), which may be reduced using the techniques of linear-scaling CC theories. The importance of selecting the proper stabilizing solution to the dRPA CARE has been highlighted. Illustrative calculations indicate that dRPA and SOSEX deliver an accuracy on par with MP2. SOSEX does not always improve the agreement with MP2.

This work was supported by the Research Council of Norway (RCN) through CoE Grant No. 179568/V30, by the Norwegian Supercomputing Program (Grant No. NN4654K), and by the European Research Council under the European Union Seventh Framework Program through the Advanced Grant ABACUS, ERC Grant Agreement No. 267683. A.M.T. was supported by the RCN (Grant No. 197446/V30) and a Royal Society university research fellowship. S.C. was supported by the Italian Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR), PRIN2009 Grant No. 2009C28YBF_001.

1.
A.
Heßelmann
and
A.
Görling
,
Mol. Phys.
109
,
2473
(
2011
).
2.
H.
Eshuis
,
J. E.
Bates
, and
F.
Furche
,
Theor. Chem. Acc.
131
,
1084
(
2012
).
3.
X.
Ren
,
P.
Rinke
,
C.
Joas
, and
M.
Scheffler
,
J. Mater. Sci.
47
,
7447
(
2012
).
4.
T.
Helgaker
,
P.
Jørgensen
, and
J.
Olsen
,
Molecular Electronic-Structure Theory
(
Wiley
,
Chichester
,
2000
).
5.
P.
Nozières
and
D.
Pines
,
Phys. Rev.
111
,
442
(
1958
).
6.
A.
Grüneis
,
M.
Marsman
,
J.
Harl
,
L.
Schimka
, and
G.
Kresse
,
J. Chem. Phys.
131
,
154115
(
2009
).
7.
G. E.
Scuseria
,
T. M.
Henderson
, and
D. C.
Sorensen
,
J. Chem. Phys.
129
,
231101
(
2008
).
8.
T.
Helgaker
, and
P.
Jørgensen
, in
Methods in Computational Molecular Physics
, edited by
S.
Wilson
and
G. H. F.
Diercksen
(
Plenum
,
New York
,
1992
), pp.
353
421
.
9.
A. D.
McLachlan
and
M. A.
Ball
,
Rev. Mod. Phys.
36
,
844
(
1964
).
10.
T. M.
Henderson
and
G. E.
Scuseria
,
Mol. Phys.
108
,
2511
(
2010
).
11.
J.
Olsen
,
H. J. A.
Jensen
, and
P.
Jørgensen
,
J. Comput. Phys.
74
,
265
(
1988
).
12.
O.
Christiansen
,
A.
Halkier
,
H.
Koch
,
P.
Jørgensen
, and
T.
Helgaker
,
J. Chem. Phys.
108
,
2801
(
1998
).
13.
K.
Hald
,
A.
Halkier
,
P.
Jørgensen
,
S.
Coriani
,
C.
Hättig
, and
T.
Helgaker
,
J. Chem. Phys.
118
,
2985
(
2003
).
14.
P.
Pulay
,
Chem. Phys. Lett.
73
,
393
(
1980
).
15.
C.-H.
Guo
and
A. J.
Laub
,
SIAM J. Matrix Anal. Appl.
21
,
694
(
2000
).
16.
DALTON, A Molecular Electronic Structure Program, Release Dalton2011, 2011; see http://www.daltonprogram.org.
17.
F.
Pawłowski
,
P.
Jørgensen
,
J.
Olsen
,
F.
Hegelund
,
T.
Helgaker
,
J.
Gauss
,
K. L.
Bak
, and
J. F.
Stanton
,
J. Chem. Phys.
116
,
6482
(
2002
).
18.
F.
Tran
and
J.
Hutter
,
J. Chem. Phys.
138
,
204103
(
2013
).
19.
M.
Dion
,
H.
Rydberg
,
E.
Schröder
,
D. C.
Langreth
, and
B. I.
Lundqvist
,
Phys. Rev. Lett.
92
,
246401
(
2004
);
[PubMed]
M.
Dion
,
H.
Rydberg
,
E.
Schröder
,
D. C.
Langreth
, and
B. I.
Lundqvist
,
Phys. Rev. Lett.
95
,
109902
E
(
2005
) (Erratum).
20.
K.
Lee
,
É. D.
Murray
,
L.
Kong
,
B. I.
Lundqvist
, and
D. C.
Langreth
,
Phys. Rev. B.
82
,
081101
(R) (
2010
).
21.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
22.
O. A.
Vydrov
and
T.
Van Voorhis
,
J. Chem. Phys.
133
,
244103
(
2010
).
23.
T.
Björkman
,
Phys. Rev. B
86
,
165109
(
2012
).
24.
K. T.
Tang
and
J. P.
Toennis
,
J. Chem. Phys.
118
,
4976
(
2003
).
25.
S. M.
Cybulski
and
R. R.
Toczyłowski
,
J. Chem. Phys.
111
,
10520
(
1999
).