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^{4} 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 properties^{8} 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^{9}

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

^{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

Here, ɛ_{p} is the orbital energy of the *p*th 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

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

For a closed-shell HF reference state, the dRPA correlation energy is given by^{7}

where the amplitude matrix

**) and de-excitation (**

*X***) parts of the eigenvectors of the symplectic eigenvalue problem**

*Y*Here, ** ω** = diag(ω

_{1}, ω

_{2}, ⋅⋅⋅, ω

_{N}) are the dRPA singlet excitation energies and

The positive definiteness of ** B** guarantees the existence of

**.**

*Z*^{7}As shown in Ref. 7,

**may be obtained directly from the CARE**

*Z*In SOSEX,^{6} exchange screening is reintroduced by calculating the energy from the dRPA amplitudes via

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

_{pq}= −κ

_{qp}such that

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

**and**

*κ***gives the Lyapunov equation for**

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

where

**depends also on**

*η***and**

*Z***from Eq. (2) and then**

*κ***from Eq. (8); next, we obtain**

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

where ** D** and

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

*d*^{8}

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*

_{μν}= ∑

_{pq}

*C*

_{μp}

*D*

_{pq}

*C*

_{νq}and

*d*

_{μνκλ}= ∑

_{pqrs}

*C*

_{μp}

*C*

_{νq}

*d*

_{pqrs}

*C*

_{κr}

*C*

_{λs}where

**contains the canonical MO coefficients at the reference geometry. As in CC gradient theory,**

*C*^{13}the reorthonormalization term is computed through a generalized effective Fock matrix backtransformed to AO basis,

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

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

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

**and**

*Z***is positive definite.**

*B*^{7}If (as assumed)

**is stabilizing, then**

*Z*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 by^{15}

where

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

The initial trial vector

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

**is diagonal dominant,**

*A*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

## II. RESULTS

All calculations have been performed with a development version of DALTON^{16} using HF reference orbitals.

Potential-energy curves for H_{2} 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* < 8*a*_{0}, when the MP2 start guess, Eq. (19), is used, albeit more iterations (up to 33) are required for *r* > 4 *a*_{0}. For *r* > 8 *a*_{0}, 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.

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_{3}, CH_{2}, F_{2}, C_{2}H_{4}, H_{2}, HCN, HNC, HNO, CO, CO_{2}, CH_{4}, HOF, N_{2}H_{2}, N_{2}, CH_{2}O, H_{2}O) 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.

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_{2} 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^{18} investigated the performance of some popular nonlocal dispersion corrections—namely, the van der Waals density functionals (vdW-DF^{19} and vdW-DF2^{20}) 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 He_{2}, Ne_{2}, and Ar_{2} 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}

. | He_{2}
. | Ne_{2}
. | Ar_{2}
. | |||
---|---|---|---|---|---|---|

. | 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_{2}
. | Ne_{2}
. | Ar_{2}
. | |||
---|---|---|---|---|---|---|

. | 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*^{6}), 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.