Skip to Main Content

Skip Nav Destination

Article navigation

Rapid Communication| August 27 2013

Communication: Analytic gradients in the random-phase approximation

Johannes Rekkedal;

Johannes Rekkedal

1Centre for Theoretical and Computational Chemistry, Department of Chemistry,

University of Oslo

, P.O. Box 1033 Blindern, N-0315 Oslo,

Norway

Search for other works by this author on:

PubMed

Google Scholar

Sonia Coriani;

Sonia Coriani

2Dipartimento di Scienze Chimiche e Farmaceutiche,

Università degli Studi di Trieste

, Via Licio Giorgieri 1, IT-34127 Trieste,

Italy

Search for other works by this author on:

PubMed

Google Scholar

Maria Francesca Iozzi;

Maria Francesca Iozzi

1Centre for Theoretical and Computational Chemistry, Department of Chemistry,

University of Oslo

, P.O. Box 1033 Blindern, N-0315 Oslo,

Norway

Search for other works by this author on:

PubMed

Google Scholar

Andrew M. Teale;

Andrew M. Teale

3School of Chemistry, University of Nottingham,

University Park

, Nottingham NG7 2RD,

United Kingdom

Search for other works by this author on:

PubMed

Google Scholar

Trygve Helgaker;

Trygve Helgaker

1Centre for Theoretical and Computational Chemistry, Department of Chemistry,

University of Oslo

, P.O. Box 1033 Blindern, N-0315 Oslo,

Norway

Search for other works by this author on:

PubMed

Google Scholar

Thomas Bondo Pedersen

Thomas Bondo Pedersen

1Centre for Theoretical and Computational Chemistry, Department of Chemistry,

University of Oslo

, P.O. Box 1033 Blindern, N-0315 Oslo,

Norway

Search for other works by this author on:

PubMed

Google Scholar

J. Chem. Phys. 139, 081101 (2013)

https://doi.org/10.1063/1.4819399

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.

I. INTRODUCTION AND 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 problems⁴ 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)⁵ and dRPA with second-order screened exchange (SOSEX).⁶ Here, we examine the relationship between the rCCD-based implementations of these approaches and their associated continuous algebraic Riccati equations (CAREs).⁷ 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 properties⁸ 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 formula⁹

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

E_{c} = \frac{1}{2} \sum_{K = 1}^{N} (ω_{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. is inconvenient for computation of the energy and its derivatives, and instead we set up a variational Lagrangian for the total energy.⁸

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

\begin{equation}F_{pq} = \delta _{pq} \varepsilon _p.\end{equation}

F_{p q} = δ_{p q} ɛ_{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}

F_{p q} = h_{p q} + \sum_{i} [2 (p q | i i) - (p i | i q)],

(3)

where h_pq 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} = 2 \sum_{i} h_{i i} + \sum_{i j} [2 (i i | j j) - (i j | j i)] .

(4)

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

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

E_{c} = \frac{1}{2} Tr (ZB) = \frac{1}{2} \sum_{a i b j} Z_{a i b j} B_{b j a i},

(5)

where the amplitude matrix

${\bm Z}=\bm{YX}^{-1}$

Z = {YX}^{- 1}

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}

(\begin{matrix} A & B \\ B & A \end{matrix}) (\begin{matrix} X \\ Y \end{matrix}) = (\begin{matrix} X ω \\ - Y ω \end{matrix}) .

(6)

Here, ω = diag(ω₁, ω₂, ⋅⋅⋅, ω_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}

A_{a i b j} = F_{b a} δ_{i j} - F_{i j} δ_{a b} + B_{a i b j}, B_{a i b j} = 2 (a i | b j) .

(7)

The positive definiteness of B guarantees the existence of

${\bm X}^{-1}$

X^{- 1}

and the negative definiteness of Z.⁷ 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,⁶ 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}

E_{c} = \frac{1}{2} \sum_{a i b j} Z_{a i b j} (B_{b j a i} - \frac{1}{2} B_{b i a j}) .

(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.¹⁰

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$

{\tilde{ϕ}}

⁠, Eq. determines the orbital rotation parameters κ_pq = −κ_qp such that

$\phi _p = \sum _q \tilde{\phi }_q [ \exp ( -\bm{\kappa } ) ]_{qp}$

ϕ_{p} = \sum_{q} {\tilde{ϕ}}_{q} {[\exp (- κ)]}_{q p}

⁠. Using Eqs. as constraints and introducing the Lagrange multipliers

$\bar{{\bm Z}}$

\bar{Z}

and

$\bar{\bm{\kappa }}$

\bar{κ}

⁠, 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}

\begin{matrix} L = & E_{HF} (κ) + E_{c} (κ, Z) \\ + \sum_{a i \geq b j} {\bar{Z}}_{a i b j} R_{a i b j} (κ, Z) + \sum_{p > q} {\bar{κ}}_{p q} F_{p q} (κ) . \end{matrix}

(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}}$

\bar{Z}

are symmetric, κ and

$\bar{\bm{\kappa }}$

\bar{κ}

are antisymmetric and symmetric, respectively, and

$\bar{\kappa }_{pp}=0$

{\bar{κ}}_{p p} = 0

⁠. The variational conditions on the multipliers give Eqs. . The variational condition on Z gives the Lyapunov equation for

$\bar{{\bm Z}}$

\bar{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}

G \bar{Z} + \bar{Z} G^{T} = - 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 }}$

\bar{κ}

⁠,

\begin{equation}\bar{\bm{\kappa }} \bm{\cal {A}} = -\bm{\eta },\end{equation}

\bar{κ} A = - η,

(12)

where

$\bm{\cal {A}}$

A

depends on orbital energies and two-electron integrals, while η depends also on Z and

$\bar{{\bm Z}}$

\bar{Z}

⁠. To construct the Lagrangian, we compute κ from Eq. and then Z from Eq. ; next, we obtain

$\bar{{\bm Z}}$

\bar{Z}

from Eq. and finally

$\bar{\bm{\kappa }}$

\bar{κ}

from Eq. . The same procedure applies to SOSEX, with modified right-hand sides of Eqs. .

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 = \sum_{p q} D_{p q} h_{p q} + \frac{1}{2} \sum_{p q r s} d_{p q r s} (p q | r s),

(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,⁸

\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}

\begin{matrix} L^{(1)} & = & \sum_{μ ν} D_{μ ν} h_{μ ν}^{(1)} + \frac{1}{2} \sum_{μ ν κ λ} d_{μ ν κ λ} {(μ ν | κ λ)}^{(1)} \\ - \sum_{μ ν} F_{μ ν}^{eff} S_{μ ν}^{(1)}, \end{matrix}

(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_μpD_pqC_νq and d_μνκλ = ∑_pqrsC_μpC_νqd_pqrsC_κrC_λs where C contains the canonical MO coefficients at the reference geometry. As in CC gradient theory,¹³ 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}

F_{p q}^{eff} = \sum_{t} D_{p t} h_{t q} + \frac{1}{2} \sum_{t r s} (d_{p t r s} + d_{t p r s}) (t q | r s) .

(15)

The solution of the CARE, Eq. (8), can be challenging.¹⁰ 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. 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}

(\begin{matrix} A & B \\ B & A \end{matrix}) (\begin{matrix} Y \\ X \end{matrix}) = - (\begin{matrix} Y ω \\ - X ω \end{matrix})

(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}$

\tilde{Z} = {XY}^{- 1} = Z^{- 1}

⁠. Both Z and

${\bm Z}^{-1}$

Z^{- 1}

are guaranteed to exist in dRPA since B is positive definite.⁷ If (as assumed) Z is stabilizing, then

$\bm{\tilde{Z}}$

\tilde{Z}

is not stabilizing since the eigenvalues of

${\bm G}(\bm{\tilde{Z}}) = -\bm{Y\omega Y}^{-1}$

G (\tilde{Z}) = - {Y ω Y}^{- 1}

are all negative. Use of the inverse solution to compute the correlation energy from Eq. will not reproduce the plasmon formula Eq. . 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}$

E_{c} = - Tr (\tilde{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)¹⁴ to accelerate convergence. The iterative procedure is defined by¹⁵

\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)})$

G^{T} (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}

N_{a i b j}^{(k)} = - \frac{R_{a i b j} (Z^{(k)})}{G_{a i a i} (Z^{(k)}) + G_{b j b j} (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}

Z_{a i b j}^{(0)} = - \frac{B_{a i b j}}{ɛ_{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 guess¹⁰ 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}

Z_{a i b j}^{(0)} = - \frac{B_{a i b j}}{A_{a i a i} + A_{b j b j}} .

(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. .

II. RESULTS

All calculations have been performed with a development version of DALTON¹⁶ using HF reference orbitals.

Potential-energy curves for H₂ 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 < 8a₀, when the MP2 start guess, Eq. (19), is used, albeit more iterations (up to 33) are required for r > 4 a₀. For r > 8 a₀, 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.

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).

View large Download slide

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).

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 (NH₃, CH₂, F₂, C₂H₄, H₂, HCN, HNC, HNO, CO, CO₂, CH₄, HOF, N₂H₂, N₂, CH₂O, H₂O) with accurate reference geometries. These were obtained by Pawłowski et al.¹⁷ 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.

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

View large Download slide

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

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 F₂ 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 Hutter¹⁸ investigated the performance of some popular nonlocal dispersion corrections—namely, the van der Waals density functionals (vdW-DF¹⁹ and vdW-DF2²⁰) of Langreth, Lundqvist and co-workers, the D3 dispersion correction of Grimme,²¹ and Björkman's revised Vydrov-Van Voorhis functional (rVV10).^22,23 Equilibrium distances and binding energies for He₂, Ne₂, and Ar₂ 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.²⁵

Table I.

Equilibrium distances (r_e, 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.

	He₂		Ne₂		Ar₂
	r_e	ΔE	r_e	ΔE	r_e	Δ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

	He₂		Ne₂		Ar₂
	r_e	ΔE	r_e	ΔE	r_e	Δ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

III. CONCLUDING REMARKS

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(M⁶), 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.

ACKNOWLEDGMENTS

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.

REFERENCES

1.

A.

Heßelmann

and

A.

Görling

,

Mol. Phys.

109

,

2473

(

2011

).

https://doi.org/10.1080/00268976.2011.614282

2.

H.

Eshuis

,

J. E.

Bates

, and

F.

Furche

,

Theor. Chem. Acc.

131

,

1084

(

2012

).

https://doi.org/10.1007/s00214-011-1084-8

3.

X.

Ren

,

P.

Rinke

,

C.

Joas

, and

M.

Scheffler

,

J. Mater. Sci.

47

,

7447

(

2012

).

https://doi.org/10.1007/s10853-012-6570-4

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

).

https://doi.org/10.1103/PhysRev.111.442

6.

A.

Grüneis

,

M.

Marsman

,

J.

Harl

,

L.

Schimka

, and

G.

Kresse

,

J. Chem. Phys.

131

,

154115

(

2009

).

https://doi.org/10.1063/1.3250347

7.

G. E.

Scuseria

,

T. M.

Henderson

, and

D. C.

Sorensen

,

J. Chem. Phys.

129

,

231101

(

2008

).

https://doi.org/10.1063/1.3043729

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

).

https://doi.org/10.1103/RevModPhys.36.844

10.

T. M.

Henderson

and

G. E.

Scuseria

,

Mol. Phys.

108

,

2511

(

2010

).

https://doi.org/10.1080/00268976.2010.507227

11.

J.

Olsen

,

H. J. A.

Jensen

, and

P.

Jørgensen

,

J. Comput. Phys.

74

,

265

(

1988

).

https://doi.org/10.1016/0021-9991(88)90081-2

12.

O.

Christiansen

,

A.

Halkier

,

H.

Koch

,

P.

Jørgensen

, and

T.

Helgaker

,

J. Chem. Phys.

108

,

2801

(

1998

).

https://doi.org/10.1063/1.475671

13.

K.

Hald

,

A.

Halkier

,

P.

Jørgensen

,

S.

Coriani

,

C.

Hättig

, and

T.

Helgaker

,

J. Chem. Phys.

118

,

2985

(

2003

).

https://doi.org/10.1063/1.1531106

14.

P.

Pulay

,

Chem. Phys. Lett.

73

,

393

(

1980

).

https://doi.org/10.1016/0009-2614(80)80396-4

15.

C.-H.

Guo

and

A. J.

Laub

,

SIAM J. Matrix Anal. Appl.

21

,

694

(

2000

).

https://doi.org/10.1137/S0895479898348519

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

).

https://doi.org/10.1063/1.1459782

18.

F.

Tran

and

J.

Hutter

,

J. Chem. Phys.

138

,

204103

(

2013

).

https://doi.org/10.1063/1.4807332

19.

M.

Dion

,

H.

Rydberg

,

E.

Schröder

,

D. C.

Langreth

, and

B. I.

Lundqvist

,

Phys. Rev. Lett.

92

,

246401

(

2004

);

https://doi.org/10.1103/PhysRevLett.92.246401

[PubMed]

M.

Dion

,

H.

Rydberg

,

E.

Schröder

,

D. C.

Langreth

, and

B. I.

Lundqvist

,

Phys. Rev. Lett.

95

,

109902

–

E

(

2005

) (Erratum).

https://doi.org/10.1103/PhysRevLett.95.109902

20.

K.

Lee

,

É. D.

Murray

,

L.

Kong

,

B. I.

Lundqvist

, and

D. C.

Langreth

,

Phys. Rev. B.

82

,

081101

(R) (

2010

).

https://doi.org/10.1103/PhysRevB.82.081101

21.

S.

Grimme

,

J.

Antony

,

S.

Ehrlich

, and

H.

Krieg

,

J. Chem. Phys.

132

,

154104

(

2010

).

https://doi.org/10.1063/1.3382344

22.

O. A.

Vydrov

and

T.

Van Voorhis

,

J. Chem. Phys.

133

,

244103

(

2010

).

https://doi.org/10.1063/1.3521275

23.

T.

Björkman

,

Phys. Rev. B

86

,

165109

(

2012

).

https://doi.org/10.1103/PhysRevB.86.165109

24.

K. T.

Tang

and

J. P.

Toennis

,

J. Chem. Phys.

118

,

4976

(

2003

).

https://doi.org/10.1063/1.1543944

25.

S. M.

Cybulski

and

R. R.

Toczyłowski

,

J. Chem. Phys.

111

,

10520

(

1999

).

https://doi.org/10.1063/1.480430

© 2013 AIP Publishing LLC.

2013

AIP Publishing LLC