Adiabatic Connection without Coupling Constant Integration - Journal

---

Subscriber access provided by Kaohsiung Medical University

Quantum Electronic Structure

Adiabatic connection without coupling constant integration Jefferson E Bates, Niladri Sengupta, Jonathon Sensenig, and Adrienn Ruzsinszky J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00067 • Publication Date (Web): 07 May 2018 Downloaded from http://pubs.acs.org on May 10, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Adiabatic connection without coupling constant integration Jefferson E. Bates,∗,† Niladri Sengupta,‡ Jonathon Sensenig,‡ and Adrienn Ruzsinszky‡ †Department of Chemistry, Appalachian State University, Boone, North Carolina 28607, United States ‡Department of Physics, Temple University, Philadelphia, Pennsylvania 19122, United States E-mail: [email protected]

Abstract Using a second-order approximation to Random Phase Approximation renormalized (RPAr) many-body perturbation theory for the interacting density-density response function, we have developed a so-called higher-order terms (HOT) approximation for the correlation energy. In combination with the first-order RPAr correction, our new method faithfully captures the infinite-order correlation for a given exchangecorrelation kernel, yielding errors of the total correlation energy on the order of 1% or less for most systems. For exchange-like kernels, our new method has the further benefit that the coupling-strength integration can be completely eliminated resulting in a modest reduction in computational cost compared to the traditional approach. When the correlation energy is accurately reproduced by the HOT approximation, structural properties and energy differences are also accurately reproduced, as we demonstrate for several periodic solids and some molecular systems. Energy differences involving

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

fragmentation are challenging for the HOT method, however, due to errors that may not cancel between a composite system and its constituent pieces.

1

Introduction

The Random Phase Approximation (RPA) has established iteself as an increasingly common Density Functional Theory 1,2 (DFT) based correlation method. 3–6 RPA is itself the simplest approximation within the adiabatic connection fluctutation-dissipation (ACFD) theorem formulation of DFT, 7–9 which is an exact formulation of the correlation energy. Since the ACFD approach accounts for density fluctuations, dispersion interactions are naturally included. ACFD methods are also inherently self-interaction error free in their exchange energies, thus the ACFD is a straightforward route to eliminating two major challenges of approximate semilocal density functionals. 5,6,10 The catch for these improvements is the increased computational cost associated with the algorithms for evaluating the correlation energy, however the difference in formal scaling between ACFD and semilocal DFT calculations can be made very close. 11–16 RPA has been demonstrated to deliver improved predictions compared to the semilocal functional used to evaluate it for thermochemistry, structural properties, and weak-interactions. 17–32 In spite of the accuracy of RPA for a wide range of problems, there are situations where it is unable to deliver high accuracy. Non-isogyric processes that break the number of electron pairs are one such class of energy differences where RPA struggles due to its improper treatment of short-ranged correlation. 33–36 Beyond-RPA (bRPA) methods are needed in such cases, and within the ACFD formalism can be easily obtained through the introduction of an exchange-correlation (xc) kernel from time dependent DFT into the solution of the interacting density response function. 9,37–40 Approaches to include exchange effects beyond RPA have also been developed from Coupled Cluster (CC) and other wavefunction theories, 41 due to the connection between CC doubles and the adiabatic connection. 42,43 The

2

ACS Paragon Plus Environment

Page 2 of 40

Page 3 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

most widely used of these methods include SOSEX 44,45 and AC-SOSEX, 40,43,46,47 which are closely related, and can be cast as an approximation to the first-order RPA renormalized perturbation theory developed in Ref. 40. In this work we will focus on corrections from xc kernels within the ACFD. Many exchange-correlation kernels have been suggested in the literature, 37,40,46,48–55 and the major distinction between most of them is whether they are derived from a model system such as the electron gas, or from many-body theory. 44 A further distinction is based on the coupling-strength dependence: for exchange-like kernels, uniform coordinate scaling shows these are linear in the coupling-strength, while for exchange-correlation kernels the dependence is non-linear. 37 Herein we will make use of several exchange-like kernels originating from models of the homogeneous electron gas (HEG), since they were demonstrated to yield accurate structural properties and energy differences at a fraction of the cost of an exchangecorrelation kernel calculation for bulk solids. 56–58 The non-local, energy-optimized (NEO) kernel 46 and the spatially renormalized adiabatic DFT (rADFT) kernels 52,53 are among the simplest exchange-like kernels available and will be combined with the ACFD herein to study both model and real systems. Instead of utilizing the traditional ACFD approach for computing the bRPA correlation energy, RPA renormalized (RPAr) perturbation theory 40,46,59 was shown to deliver a robust treatment of bRPA correlation with several advantages compared to the infinite-order method. For exchange-like kernels, an analytic integral exists for the first-order (RPAr1) correction. 40,46 RPAr1 also guarantees the elimination of divergences due to electronic instabilities. 27,40,60 Compared to the infinite-order method, we showed previously that RPAr1 systematically underestimates the bRPA correction by approximately 10%, 46 which can be straightforwardly accounted for through the inclusion of higher-orders in the RPAr series. 59 The drawback of higher-orders in the series is that an analytic coupling-strength integral is not currently known for those terms and thus the coupling-strength integration must be reintroduced just to compute a small fraction of the energy.

3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

To circumvent the coupling-strength integration entirely, we have developed a higherorder terms (HOT) approximation that seeks to account for the missing correlation beyond RPAr1. By studying the coupling-strength dependence of second and higher orders of RPAr in the electron gas, as well as the kinetic and potential contributions to correlation for the interacting system at full-coupling, we find that a one-parameter expression for the beyondRPAr1 correlation reproduces the infinite-order result to within about 1%. When the HOT approximation accurately reproduces the infinite-order bRPA correlation for a given system and kernel, structural properties and energy differences are also accurate. For open-shell systems our approximation is less robust, but we are able to understand the origin of this deficiency in terms of the behavior of the RPAr series in these cases. The remainder of the paper is structured as follows. In Sec. 2 we give a brief review of the ACFD and RPA renormalization formalisms, as well as introduce the HOT approximation. In Sec. 3 we analyze the behavior of the HOT approximation within two contexts for the electron gas and follow with results for real solids. We also present our results for a few select molecular systems, with a focus on fragmentation processes, to illustrate the challenges for future developments. Finally, in Sec. 4 we give a brief overview of the results and discuss their implications for beyond-RPA correlation methods.

2 2.1

Theory Adiabatic Connection Fluctuation-Dissipation DFT

Conceptually, the adiabatic connection approach in DFT links the non-interacting, KohnSham (KS) system to the interacting system through the introduction of a coupling-strength, λ, associated with the electron-electron interaction. 5,61,62 The total energy can then be obtained from the Hellman-Feynman theorem by integrating over all intermediate systems between the KS system at λ = 0 and the interacting system at λ = 1. 63 The exchangecorrelation (xc) energy can be isolated from the kinetic and external potential contributions, 4

ACS Paragon Plus Environment

Page 4 of 40

Page 5 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

and expressed as an integral over the potential contribution to exchange and correlation. 1

Z Exc [ρ] =

dλ Uxc [ρ](λ) .

(1)

0

Furthermore, the exchange piece is straightforward to evaluate and separate. The resulting expression for the total energy, neglecting correlation, is equivalent to the exact exchange (EXX), or Hartree-Fock, total energy, evaluated using orbitals from a semilocal DFT calculation. An exact expression for the correlation energy can then be derived based on the fluctuationdissipation theorem at zero temperature. 3,7 The correlation potential provides the connection between density fluctuations described by time dependent DFT (TDDFT) to the ground state pair-density which describes electron correlation 5 Z Ec [ρ] = −

1

Z

∞

dλ Im 0

0

dω hV (χλ (ω) − χ0 (ω))i 2π

(2)

where V is the direct Coulomb (Hartree) interaction, Im indicates the imaginary part, hAi is the trace of matrix A, and atomic units are used unless otherwise specified. The densitydensity response functions χλ and χ0 satisfy a Dyson-type equation 0

0

Z

χλ (x, x ; ω) =χ0 (x, x ; ω) +

dx1 dx2 χ0 (x, x1 ; ω)

  λ × Vλ (x1 , x2 ) + fxc (x1 , x2 ; ω) χλ (x2 , x0 ; ω) ,

(3)

where x = {σ, r} is short for a set of spin and spatial coordinates, χ0 (ω) is the non-interacting KS response function, and fxc (ω) is the (exact) frequency-dependent exchange-correlation (xc) kernel. 64 The Coulomb interaction is linear in the coupling strength, Vλ = λV , and the behavior of the xc-kernel can be determined from uniform coordinate scaling, 37 though we focus here on kernels which scale linearly with λ. Once the kernel and the KS response function are available, the interacting response function can be extracted by solving Eq. (3) 5

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 40

and the correlation energy computed from Eq. (2). We will call a direct solution of Eq. (3) and the corresponding correlation energy including the xc kernel as the “infinite-order” method in this work, since by inversion
−1 λ χλ (ω) = χ−1 . 0 (ω) − Vλ − fxc (ω)

(4)

Since fxc appears in the inverse it has been included in the response-function, and consequently the correlation energy, through infinite-order in the traditional many-body perturbation theory sense.

2.2

RPA renormalization

First introduced in Ref. 40, RPA renormalization is an alternative formulation of many-body perturbation theory that naturally includes screening effects in bRPA correlation through its dependence on the RPA response function. 46,59 By a slight rearrangement in Eq. (4), the Dyson-type equation can be expressed as

λ χλ (ω) = χˆλ (ω) + χˆλ (ω)fxc (ω)χλ (ω) ,

with the RPA response function χˆλ = χ−1 0 − λV


−1

(5)

= (1 − χ0 λV )−1 χ0 . If the xc kernel is

kept to infinite-order in Eq. (5), the resulting correlation energy is equivalent to that from the traditional approach. The utility of RPAr therefore lies in approximations to Eq. (5), which we can obtain order-by-order just by expanding the inverse in powers of the product χf ˆ xc λ −1 ) χˆλ χλ = (1 − χˆλ fxc λ λ λ χˆλ + . . . ≈ χˆλ + χˆλ fxc χˆλ + χˆλ fxc χˆλ fxc

6

ACS Paragon Plus Environment

(6)

Page 7 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

The correlation energy is then obtained as a sum of the RPA correlation energy, which is zeroth-order in RPAr, and all bRPA contributions

Ec [fxc ] = EcRPA +

∞ X

∆EcRPAr-n [fxc ] ,

(7)

n=1

where the n-th order RPAr correction for a given kernel is

∆EcRPAr-n [fxc ]

Z =−

1

Z dλ

0

0

∞

 du

λ V (χˆλ (iu)fxc (iu))n χˆλ (iu) 2π

(8)

and the total (infinite-order) bRPA correction can be computed as 46

∆EcbRPA [fxc ] = Ec [fxc ] − EcRPA Z 1 Z ∞  du

λ V χˆλ (iu)fxc (iu) χλ (iu) . =− dλ 2π 0 0

(9)

Keeping with our previous naming conventions, 46,59 when a given kernel is referenced directly as a correlation method we mean the infinite-order method, Eq. (9). The acronym XACFD is also used to mean the exact ACFD correlation for a given kernel. We showed previously that RPAr1 systematically recovers more than 90% of the bRPA correlation energy for the NEO kernel, 46 implying that second-order and higher corrections make up the remaining 10%, provided the series expansion converges. For example, Figure 1 illustrates the breakdown of correlation contributions within RPAr. The kinetic, potential, and total correlation energy contributions can all be extracted from such a plot. While we did not analytically prove the convergence of the RPAr series, we demonstrated numerically that it does tend to converge to the exact result in both model and real solids, as well as for several finite systems using exchange-like kernels. 59 Additionally the convergence is found to be monotonic from below for the systems we studied because each correction was positive and smaller than the previous order. This trend was rooted in an analytical analysis of the correlation potential which shows RPA renormalization systematically sums subsets of terms 7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

∆U

c bRP A

0

∆EcbRPAr1

EcXACFD ∆EcbRPA −TcRPA

0

∆EcRPAr1

RPA Exact ACFD

Exact ACFD RPAr1 λ

1

0

UcRPA

−∆TcbRPA

∆UcbRPA (λ)

Uc (λ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 40

0

λ

1

Figure 1: Cartoon of the coupling strength integrands for the adiabatic-connection showing the RPA and bRPA components. RPA overestimates the correlation energy, Eq. (2), (left) and so the adiabatic connection is too negative in comparison to the exact ACFD (XACFD). The shaded region between the RPA and the XACFD curves corresponds to the infinite-order bRPA correlation energy, Eq. (9), whose adiabatic-connection is plotted on the right. RPA renormalization to first order (RPAr1, Eq. (8)) systematically recovers the dominant bRPA correlation (right), but underestimates the infinite-order method (Eq. (9)) by the shaded region between the solid and dashed lines, which can be exactly computed from Eq. (10). The HOT correction is intended to make up this remaining piece of the correlation energy. The hatched region corresponds to −Tc and −∆TcbRPA for the right and left plots, respectively. in many-body perturbation theory (MBPT) that all have the same signed contribution to the correlation energy. 59 This is quite impressive for a perturbation expansion since traditional MBPT expansions based on non-interacting references are not gauranteed to converge at higher-orders and the convergence is not necessarily monotonic. 65–68

2.3

Higher-Order Terms Approximation

For RPA and RPAr1 with exchange-like kernels 46 the coupling-strength integration can be performed analytically so it is counter-productive to evaluate higher-orders exactly for real systems because numerical integrations of λ are required. To avoid this we have studied the coupling-strength dependence of the integrands in Eq. (8) and find that an accurate approximation can be obtained from a second-order correction to the response function.

8

ACS Paragon Plus Environment

Page 9 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Starting with the exact definition of correlation beyond RPAr1,

∆EcbRPAr1 [fxc ] =EcbRPA [fxc ] − EcRPAr1 [fxc ] Z 1 Z ∞ du λ λ dλ =− hV χˆλ fxc χˆλ fxc χλ i , 2π 0 0

(10)

we decompose this correlation energy into kinetic and potential contributions 39 and write the total correction in terms of the potential contribution

∆EcbRPAr1 [fxc ] = (1 − ˆb)∆UcbRPAr1 [fxc ] , bRPA

(11)

RPAr1

− ∆Tc ˆb = − ∆Tc , bRPA ∆Uc − ∆UcRPAr1

(12)

in analogy to the adiabatic connection for the total correlation energy. 69 Such a construction corresponds to a two-node trapezoidal rule quadrature for integrating the λ-depdendence, which itself implies ˆb = 1/2. 70 These terms can be identified from Fig. 1, and are defined in the supporting information. 71 To obtain the HOT correction we approximate two quantities, ˆb ≈ 1/2 and ∆U bRPAr1 ≈ ∆U RPAr2 , arriving at the following expression for correlation beyond c c RPAr1 1 ∆EcbRPAr1 [fxc ] ≈ ∆EcHOT [fxc ] = ∆UcRPAr2 [fxc ] . 2

(13)

While ˆb can in principle be computed correctly for any system, we find that approximating it as ˆb = 1/2 is surprisingly robust and yields acceptable accuracy. In fact for several systems we find that the exact ˆb defined in this way for approximate kernels is around 0.5. 71 Furthermore, by comparing to the properly integrated RPAr-n corrections of Eq. (8) this approximation for ˆb seems to compensate for neglecting third and higher order terms in the expansion of the RPAr correlation energy. The total correlation energy within the HOT

9

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 40

approximation is computed as 1 EcHOT [fxc ] = EcRPA + ∆EcRPAr1 [fxc ] + ∆UcRPAr2 [fxc ] , 2

(14)

with

∆UcRPAr2 [fxc ]

Z

∞

=− 0

du hV χf ˆ xc χf ˆ xc χi ˆ . 2π

(15)

We note a fixed choice of ˆb could introduce size-consistency errors since the exact ˆb is a function of the molecular configuration and changes as the configuration of atoms changes for a given kernel. However, using ∆UcbRPAr1 ≈ ∆UcRPAr2 with a fixed value of ˆb will result in a size-consistent energy since ∆UcRPAr2 is itself size-consistent. This in addition to the fact that RPA and RPAr1 are themselves size-consistent. The fixed value of ˆb = 1/2 biases the HOT approximation towards lower densities of the electron gas for readily available kernels, see Fig. S1. Additionally, by rearranging Eq. (11), one can show that in the case of very-low densities or strong static correlation, ˆb → 0 as

∆EcbRPAr1 ∆UcbRPAr1

→ 1, which stems from the kinetic

contributions to correlation going to zero. 69,72 For spin-unpolarized systems, or spin-independent kernels, the coupling strength integral for ∆EcRPAr1 can be taken analytically for exchange-like kernels, 40,46 ∆EcRPAr1 [fx ]

Z =− 0

∞

 du −1  V fx (1 − χ0 V )−1 χ0 V + ln (1 − χ0 V ) , 2π

(16)

meaning only one diagonalization of χˆλ at λ = 1 is needed to compute the correlation energy, thus avoiding instabilities that come from having to diagonalize a response function that directly includes fxc . 40,59,60 The reduction in computational cost for periodic boundary conditions is related to the reduced amount of time spent needed for diagonalization of χˆλ in the reciprocal lattice basis, which is not rate-determining as the cost to build χ0 is quartic in a plane-wave basis. However, for systems such as metals where dense k-mesh samplings are

10

ACS Paragon Plus Environment

Page 11 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

needed for converged results, the reduction can be appreciable since the savings is per q-point. As stressed in Ref. 59, the RPA and beyond RPA corrections are obtained simultaneously and one does not need to perform two separate calculations in order to determine the behavior of each method. Instead, the total energy is a sum of separate pieces and the impact of the kernel compared to RPA can be directly determined. Exchange-correlation kernels or spin-polarized systems with any kernel still require numerical integration over λ within RPAr as implemented using the density-density response function in a plane-wave basis, so there is no reduction in computational cost compared to the infinite-order method, however instabilities are still avoided since only χˆλ needs to be diagonalized. For exchange-correlation kernels, there is no analytic formula for RPAr1 in terms of integrating out the λ dependence. Thus when using an exchange-correlation kernel, numerical integration over the coupling strength is always required because a way to eliminate it analytically has not been reported. For spin-polarized systems, the Coulomb interaction is spin-independent so RPA calculations can take advantage of the spin-summed response function to reduce the computational cost to that of a spin-independent system. However, once an exchange-correlation kernel is introduced, this is no longer possible because the kernel effects α and β spins differently. 27,56 Thus, for spin-dependent systems with any type of xc kernel, a numerical integration over λ is needed because one must work with the λ- and spin-dependent Dyson equation for the response function. Though we have utilized a plane-wave implementation, for localized basis sets the savings introduced by avoiding the repeated diagonalization of the RPA response function (or supermatrix 8 ) is actually significant. The diagonalization of matrices in the particle-hole basis scales as N6 , where N is a measure of the basis set size, and the diagonalization must be carried out for every coupling-strength point used to integrate the correlation energy (typically between 6 and 8 points). Thus one is able to avoid repeatedly performing the rate-determining step by eliminating the numerical integration over the coupling strength. If the electron repulsion integrals are constructed using resolution of the identity methods in

11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

an integral-direct fashion, the cost to perform a HOT calculation via numerical integration of the frequency along the imaginary axis will be the same as that for RPA 11,28,73 or RPAr1, 40 approximately O(N 4 ) with a prefactor that depends on the details of the implementation.

2.4

Computational Details

Using a modified version of the gpaw code, 74–76 we have implemented the HOT approximation in conjunction with the spin-dependent (rADFT) and spin-independent (rADFTns) forms of the exchange-like kernels rALDA 52 and rAPBE. 53 When necessary we use the notation RPAr@kernel to denote which kernel has been used with a particular RPA renormalization approximation. The rALDA and rAPBE kernels have been demonstrated to improve upon RPA for energy differences and structural properties, and preserve the accurate description of dispersion interactions for atoms, molecules, and solids. 52,53,56–58 Since the combination of these kernels and RPAr resulted in analogous trends, 59 here we will focus almost solely on the rAPBE kernel. In addition to specifying the kernel within gpaw, one must specify a reciprocal lattice vector basis averaging scheme to compute the kernel since there are multiple ways to extend homogeneous electron gas kernels to inhomogeneous systems. 52,54,57,77 We have used the socalled wavevector symmetrization throughout to ensure that the kernel is symmetric in the reciprocal lattice basis because of its computational advantages, as discussed in Ref. 57. There is the added bonus that the wavevector symmetrization also preserves the proper divergence behaviors of the kernel in the q → 0 limit. 57 See Ref. 57 for an overview of the averaging schemes. Calculations were performed within the projector-augmented wave formalism 78 using the 0.9.20000 gpaw datasets, which treat the 4s and 3d shells as a part of the valence for transition metal atoms, as well as treating the 3p semicore states for Ni. Results for real systems were obtained using PBE orbitals to construct the Kohn-Sham reference determinant, and gamma-centered Monkhorst-Pack k-point meshes 79 were used throughout. The maximum cutoff for the response function was chosen between 300 and 12

ACS Paragon Plus Environment

Page 12 of 40

Page 13 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

400 eV, the number of bands was chosen to be equal to the number of plane-waves, and the perturbative approach from Ref. 80 was used to treat the divergence of the Coulomb interaction at small wavevectors. The frequency integral was performed as in Ref. 17 using a 16 point Gauss-Legendre quadrature, and with a frequency scale of 2 for non-metals, and 2.5 for metals. 27 Extrapolation of the correlation energies to the basis set limit were performed using the Harl-Kresse 17 method with at least four cutoffs below the maximum. The cutoff used to generate the wavefunctions for the response function was 600 eV. FermiDirac occupations corresponding to an electronic temperature of 0.01 eV were used for all periodic systems. A Wigner-Seitz truncation scheme 81 was used to treat the small wavevector divergence of the Coulomb interaction in the exchange energy. For atomic and molecular systems we used rectangular boxes with unequal side lengths to break spatial symmetry, and extrapolated the correlation energy from calculations at planewave cutoffs of 250 and 300 eV for the response function. These cutoffs were previously shown to deliver small errors compared to extrapolations with higher cutoffs for RPA. 27 Larger plane-wave cutoffs, box sizes, or volume based extrapolations of the correlation energy 17 would be needed to obtain fully converged results for each kernel, however the relative performance of RPAr to the infinite-order method is usually independent of the basis for cutoffs larger than 250 eV. For atoms, a larger variation in finite cutoff and extrapolated values can result in slower convergence of RPAr for the extrapolated results. For molecules and bulk systems the differences between finite cutoff and extrapolated results tends to be negligible for any k-mesh and cutoffs above 200 eV. The computational settings needed to reproduce our results have been included in the supporting information. 71 Results for the electron gas were obtained using an in-house python code. GaussLegendre quadratures of 12-20 points were used for the frequency and coupling-strength integration, while the q-integration is implemented using the rectangle method over 0 < q ≤ 15 a.u. with 3000 points. The exchange-like nonlocal, energy-optimized (NEO) kernel 46 was also implemented in conjunction with RPAr for the electron gas, and results are presented

13

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

below for comparison to rALDA. Note that for the electron gas rALDA and rAPBE are equivalent. 58

3

Results

First we will discuss the performance of the HOT approximation for predicting total bRPA correlation energies, as well as present an alternative interpretation for its accuracy based on the behavior of the adiabatic connection integrand. After establishing a trend in behavior, we discuss the application of the HOT approximation to predicting the structural properties of several simple solids as well as for the pressure induced phase transition of Silicon Carbide. We wrap up the results section by discussing some of the challenges facing the HOT approximation for molecular and atomic systems which are often needed to compute fragmentation reactions such as the atomization energy or binding energy of dispersion-bound complexes.

3.1

HOT Approximation Compared to Finite-Order RPAr

The HOT approximation can be recast as a linear function of λ whose integral from 0 to 1 yields the prefactor 1/2 in Eq. (13). In Figure 2 we plot the coupling strength dependence of the HOT approximation and the exact bRPAr1 coupling strength integrand in the spinunpolarized HEG and diamond carbon for perspective. At small values of λ, the HOT approximation overestimates the the exact bRPAr1 correlation contribution, while for larger values of λ closer to 1 it is an underestimate. In spite of this, the total area of the HOT correction is quite accurate for, e.g., rs = 2.07 due to the error compensation between small and large values of λ. The total difference in the areas underneath these curves is less than 3 meV per electron. The error is larger for rs = 4 with NEO or for rAPBE and diamond where the HOT approximation is an underestimate and overestimate, respectively. One immediate consequence of this behavior is that the systematic underestimation of beyondRPA correlation delivered by properly λ-integrated RPA renormalization is not retained by

14

ACS Paragon Plus Environment

Page 14 of 40

Page 15 of 40

NEO - HEG rs = 2.07 bRPAr1

NEO - HEG rs = 4 bRPAr1

40

40

30

30

20

20

10

10

UCbRPAr1( ) (meV/e )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

00

bRPAr1 HOT

1 00

rAPBE - C(A4) bRPAr1

70 60 50 40 30 20 10 1 00

Figure 2: Adiabatic connection integrand for the beyond RPAr1 correlation energy per particle, Eq. (10), with the NEO kernel for the HEG at rs = 2.07 and rs = 4, and the rAPBE kernel for carbon in the diamond (A4) phase with a = 3.566 ˚ A. The HOT approximation is also plotted as a linear function of λ for comparison. The shaded region represents the difference between the HOT and exact bRPAr1 corrections. The accuracy of the HOT approximation can be seen as a compensation of errors between the small and large λ behavior since for the former HOT is an overestimate, but for the latter an underestimate. the HOT approximation, and instead the HOT correction can lead to either an over or underestimation of beyond RPA correlation. This behavior is more clearly understood by comparing the integrated beyond RPAr1 corrections. Figure 3 shows the comparison between the integrated first four orders of RPAr and the HOT approximation. The contribution of each higher order of RPAr is smaller by nearly an order of magnitude than the previous order. The fourth-order correction is not even noticable on this scale, except for rALDA where it is clearly a tiny fraction of the total. We have also included results for the Constantin-Pitarke exchange-correlation kernel 49 (CP07, non-linear in λ) for comparison. The HOT approximation understimates the total beyond RPAr1 correlation energy at this density by 2.6 meV per electron for the NEO kernel and 0.8 meV per electron for rALDA, while for the CP07 kernel the HOT method overestimates by 3 meV per electron. As a function of rs , the HOT method is typically an overestimate at high densities, but becomes an underestimate in the low-density limit where contributions from RPAr beyond second-order are more important. A plot illustrating this behavior is

15

ACS Paragon Plus Environment

UcRPAr2

1

Journal of Chemical Theory and Computation

Homogeneous Electron Gas

Real Systems - rAPBE kernel RPAr1 RPAr2 RPAr3

3.2

RPAr1 RPAr2 RPAr3 RPAr4 HOT

3.0

RPAr4 HOT

EcbRPA (eV)

2.8 2.6 2.4 2.2 O* -B1 Mg

A4 Si-

C -FC Fe

( 3P O*

(1 CO

)

)

7 CP 0

DA

NE

O

2.0 rAL

420 400 380 360 340 320 300 280 260

EcbRPA (meV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 40

Figure 3: RPA renormalization contributions to the beyond RPA correlation energy for first through fourth order in the spin-unpolarized HEG for rs = 4 and for a molecule, atom, metal, semiconductor, and insulator. The correlation energy per particle is plotted for the HEG, while correlation energy per unit cell is plotted for the real systems. The HOT correction is also plotted for comparison with RPAr2-4. The bulk of the bRPA correlation is accounted for through RPAr1 and adding the HOT correction yields an error typically less than 1%. The MgO correlation energies have been reduced by 2 eV and the O energies increased by 1 eV to improve the comparison. given in the supporting information, Fig. S1. These errors are exceedingly small, and in terms of percentage, the HOT method yields a relative error within 1% compared to the infinite-order bRPA correlation energy for these kernels over a wide range of densities. For the real systems with rAPBE and rALDA kernels, the HOT approximation yields equally accurate beyond RPAr1 correlation energies. For the molecular systems the HOT approximation is an underestimate, though more so for oxygen atom than CO. For the bulk solids, the HOT approximation is a slight overestimate of bRPAr1 correlation, but the error is still less than 1% of the total beyond RPA correlation energy. Given that the HOT approximation accurately reproduces the total correlation energy, any computable property that depends on the correlation energy should also be accurately reproduced. We explored this point through structural parameters of several simple solids, energy differences of weakly bound molecular complexes, and some cohesive and atomization energies.

16

ACS Paragon Plus Environment

Page 17 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

3.2

Structural and Energetic Properties

The lattice constants and bulk moduli of 13 solids were computed using the HOT approximation and the infinite-order method for correlation beyond RPA. We chose several metals (Cu, Ag, Pd, Na, Al), as well several gapped materials (C, Si, SiC, LiF, GaN, AlN, NaF, MgO). We used between 7 and 10 volume points around the minimum and fit the EV data to the stabilized jellium equation of state (SJEOS). 82 The SJEOS equation of state requires four fit parameters to determine the energy as a function of volume: the energy at the minimum E0 , the volume at the minimum V0 , the bulk modulus B = −V (∂P/∂V ), and its pressure derivative B 0 = (∂B/∂P )P0 . Computational settings were similar to those from Refs. 83 and 84, and are tabulated in the supporting information along with the numerical results for RPA and the kernel-corrected methods. ACFD methods are known to yield accurate structural properties in comparison to experiment, 26,27,83,84 though we stress here the degree to which the HOT approximation is capable of reproducing the infinite-order results. The agreement is excellent, lattice constants are reproduced to within 0.001 ˚ A on average and bulk moduli to within 1 GPa, Figure 4. The difference between RPAr1 and XACFD for the NEO kernel 46 was previously demonstrated to be small, and the same behavior persists for the rAPBE kernel. Though there is a small error in the energy on the order of 1% in the correlation energy, the HOT approximation reproduces both the minimum and curvature of the infinite-order EV fits. Energy differences can prove to be more challenging for the HOT approximation, especially when open-shell atoms are involved. To test the performance of the HOT approximation for energy differences with the rAPBE kernel, we computed 6 cohesive energies of bulk solids (C, Si, AlN, SiC, Al, Pd), 9 atomization energies of small molecules, the binding energies of H2 O and NH3 dimers from the S22 test set, 85,86 and the interlayer binding energy of 2H MoS2 . Previous works have made it clear that RPA struggles for energy differences that break the number of electron pairs, 5,10,27 and that inclusion of an xc kernel remarkably improves the agreement with reference values. 40,53,56,87 For weak interactions, RPA can be acceptably accurate and addition of the kernel is only 17

ACS Paragon Plus Environment

Signed Error (GPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Signed Error (pm)

Journal of Chemical Theory and Computation

2 0 2 4 6

12.5 10.0 7.5 5.0 2.5 0.0 2.5 5.0

Lattice Constant Deviations

RPA RPAr1 HOT

Bulk Modulus Deviations

C Si AlN GaN SiC MgO LiF NaF Na Al Cu Pd Ag

Figure 4: Deviations from the infinite-order (XACFD) rAPBE results of the lattice constants and bulk moduli for 13 solids. Numerical results for a0 and B are given in the supporting information. RPA generally agrees with the kernel-corrected results, except for ionic materials, such as LiF and NaF, where short-ranged forces are needed to reduce the lattice constant. Overall, the HOT method quantitatively reproduces the infinite-order structural parameters.

18

ACS Paragon Plus Environment

Page 18 of 40

Page 19 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

required for cases where there may still be some residual short-ranged interactions. For the cohesive (Table 1) and atomization energies (Table 2), RPAr1 systematically underestimates the infinite-order energy differences on average by 0.25 eV and 0.44 eV, respectively, as expected from its underestimation of the total beyond-RPA correlation. 59 Our PBE, EXX, RPA, and XACFD@rAPBE results agree reasonably well with those in the literature, 26,27,53,56 with differences likely stemming from different computational parameters (number of kpoints, cutoffs, etc). The HOT approximation typically increases the RPAr1 energy difference, reducing the average difference of the cohesive and atomization energies compared to the infinite-order method to an underestimation of 0.18 eV and 0.28 eV, respectively. For molecular and bulk systems, even though the HOT approximation makes a small error in the total energy on the order of 1-2%, the larger (in magnitude) error made for the spin-polarized atoms dominates the error in the cohesive or atomization energy. 59 For instance, if we use the infinite-order correlation energies of the atoms with the HOT approximation for the molecule or bulk solid, the average difference in cohesive and atomization energies drops to an underestimation of 0.01 eV and 0.03 eV, respectively, which is the same order of magnitude as the HOT method’s error in the total energy for the molecular or bulk system, see Tables S2 and S3 in the supporting information. Thus we recommend to replace the HOT approximation with XACFD when computing open-shell atomic correlation energies as a remedy for the poor performance of the bare HOT approximation for fragmentation reactions. A similar shortcoming of the HOT approximation is also encountered for the two weakly bound dimers, Table 3. In these systems, the error in the total correlation energy for both the monomer and dimer is roughly equal, but of opposite sign. As a result, the errors do not cancel and the HOT approximation yields interaction energies that are too large in comparison to the infinite-order method. This behavior can be traced back to the approximation ˆb ≈ 1/2 for both systems, which results in a beyond RPAr1 correction that is too large for the dimers and too small for the monomers. Fortunately for these systems, however, RPAr1

19

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 1: Cohesive energies (∆Ecoh , eV) for 6 simple solids. The HOT approximation recovers most of the remaining binding missed by RPAr1 in comparison to the infinite-order method (the exact ACFD, XACFD). A4 corresponds to the diamond structure while B3 is zincblende. Our XACFD@rAPBE results agree to within 0.15 eV of those in Ref. 53. Mean signed errors and mean absolute errors with respect to experiment are also included. ∆Ecoh (eV) Pd (FCC) C (A4) Si (A4) Al (FCC) AlN (B3) SiC (B3) MSE MAE

PBE EXX RPA RPAr1 HOT XACFD Expt. 3.68 –1.38 3.81 3.46 3.45 3.49 3.94 7.71 5.15 7.02 7.09 7.20 7.49 7.55 4.55 2.83 4.33 4.39 4.47 4.75 4.68 3.43 1.34 3.17 3.26 3.35 3.49 3.43 5.68 3.66 5.54 5.42 5.47 5.54 5.85 6.37 4.32 5.99 6.02 6.12 6.40 6.48 –0.08 –2.67 –0.35 –0.38 –0.31 –0.13 0.14 2.67 0.35 0.38 0.31 0.17

Table 2: Atomization energies (∆EAE , eV) for 9 small molecules. The HOT approximation typically increases the RPAr1@rAPBE result by 0.16 eV on average, improving the agreement with both experiment and the infinite-order result (XACFD). Our XACFD@rAPBE results agree to within 5% of those in Ref. 53. Mean signed errors and mean absolute errors with respect to experiment are also included. ∆Ecoh (eV) PBE EXX RPA RPAr1 C2 H2 17.97 12.56 16.64 16.95 C2 H4 24.78 18.43 23.36 23.63 CO2 18.04 9.89 15.74 16.24 CO 11.65 7.19 10.61 10.89 H2 4.53 3.64 4.72 4.69 H2 O 10.16 6.64 9.59 9.73 N2 10.56 4.57 9.79 9.79 NH3 13.10 8.66 12.61 12.58 O2 6.25 1.04 4.96 5.08 MSE 0.37 -3.72 -0.29 -0.21 MAE 0.44 3.72 0.29 0.21

20

HOT XACFD Expt. 17.20 17.77 17.56 23.93 24.51 24.46 16.52 17.08 16.87 11.06 11.47 11.23 4.72 4.72 4.73 9.86 10.00 10.10 9.83 9.82 9.89 12.68 12.69 12.88 5.20 5.47 5.20 -0.11 0.02 0.11 0.15

ACS Paragon Plus Environment

Page 20 of 40

Page 21 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

recovers the dominant beyond-RPA contribution to the binding energy and yields a small error compared to the infinite-order method. Table 3: Binding energies (kcal/mol) for the ammonia and water dimers from the S22 test set, and inter-layer binding energy for MoS2 (meV/˚ A2 ). Reference values were taken from Takatani et al 86 for the S22 dimers, and from Ref. 21 for MoS2 . The reported results for RPA are within 1 kcal/mol to those of Ref. 56, from which the rALDA results were taken for comparison. RPAr1@rAPBE agrees with the infinite-order results (XACFD@rAPBE) to within 0.5 kcal/mol, but addition of the HOT approximation grossly overestimates the total binding energy due to a poor cancellation of errors. The binding energy of MoS2 , on the other hand, is accurately reproduced by the HOT method since it involves “bulk-like” energy differences. (NH3 )2 (H2 O)2 MoS2

PBE EXX 2.60 1.14 4.61 3.81 0.63 –14.9

RPA RPAr1 HOT XACFD rALDA 56 2.02 2.99 4.52 3.12 2.31 3.28 4.24 6.21 4.46 3.92 22.4 20.8 20.8 21.9 –

Ref. 3.15 5.07 20.5

To understand the challenge facing the HOT method for these energy differences, Figure 5 shows the error of the HOT approximation in comparison to the infinite-order beyond RPAr1 correction as a function of ˆb for the NH3 dimer, NH3 , and the binding energy. If ˆb is chosen to be the same value for both NH3 and NH3 dimer, then the value that yields the smallest error is close to 1, i.e. RPAr1 without correction since ∆EcHOT ∝ (1 − ˆb). If ˆb can be chosen independently for each system, the error in the binding energy is minimized when the error in the total correlation energies for each system is zero; a value of ˆb = 0.55 for NH3 is needed to accomplish this, while a value of 0.5 was already quite accurate for the NH3 dimer. The same challenge also applies to the cohesive and atomization energies in that choosing one value of ˆb for both atoms and composite systems is limited in accuracy due to poor cancellation of errors between the atomic errors and the bulk errors.

3.3

Transition Pressure of SiC and Binding Energy of MoS2

Though the HOT approximation struggles for energy differences describing fragmentation into monomer or atomic species, energy differences that involve only bulk or surface species can still be accurately reproduced. This property of the method is important for, e.g., 21

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

NH3 Correlation (NH3)2 Correlation Binding Energy

EcbRPAr1 (kcal/mol)

6

EcHOT

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 40

4 2 0 2 4 0.0

0.2

0.4

b

0.6

0.8

1.0

Figure 5: Error of the HOT approximation (∆EcHOT ) compared to the exact beyond RPAr1 (∆EcbRPAr1 ) correlation energy for NH3 , NH3 dimer, and the binding energy correction as a function of ˆb. The optimal fixed value of ˆb for both systems is close to 1, meaning RPAr1 without correction. If ˆb can vary from system to system, the binding energy error is near zero for ˆb = 0.49 for NH3 dimer and 0.55 for NH3 . The dashed vertical line indicates ˆb = 0.5. predicting the relative phase ordering of different materials or organic molecules. 25,88–93 Two simple examples are the pressure induced phase transition between two phases of a material and the interlayer binding energy of 2H MoS2 . The transition pressure (Pt ) at which the phase transition begins can be computed from the common tangent between the EV curves of each material or through the equal enthalpy condition, H1 (V1 , Pt ) = H2 (V2 , Pt ), where H(V, P ) = E(V ) + P V (P ) is the enthalpy. Using a parametrized EOS to determine the energy as a function of volume and the volume as a function of pressure, Pt can be determined through numerical optimization. The biggest factors that determine the transition pressure are the energy and volume differences between phases since one can think of the common tangent as being Pt ≈ (E1 − E2 )/(V1 − V2 ) for the two phases. An accurate prediction of the energy and volume difference between phases is therefore a central part of predicting the transition pressure. As an example of such a process, we have computed the transition pressure of Silicon Carbide from the zincblende (ZB) phase to the rocksalt (RS) phase. Both phases are semiconductors, and experimentally the transition pressure is observed to be on the order of 100 22

ACS Paragon Plus Environment

Page 23 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

GPa. 94,95 Semilocal functionals underestimate the experimental transition pressure due to a too small energy difference between the phases. 95–98 A hybrid functional, such as B3LYP, increases the energy difference and subsequently the transition pressure between the phases due to the addition of exact exchange. 99 RPA has been applied previously for other phase transitions, such as in Si 14 and SiO2 , 22 however it does not always yield smaller errors in comparison to experiment. In such cases a kernel correction could improve the result, and we found for many simple solids this is indeed the case. 100 Here we use PBE input orbitals and the rAPBE kernel to evaluate RPAr1, HOT, and infinite-order method. The computational parameters can be found in the supporting information. For the ZB→RS transition in SiC, our PBE results agree well with previously reported LDA and GGA values in the literature, 84,95–97,101 Table 4, for both the structural parameters and the transition pressure. Though not discussed in previous reports of RPA transition pressures, the EXX and RPA correlation contributions to the transition pressure are often large and of opposite sign. In this case, EXX yields a transition pressure of 114 GPa and an energy difference of 2.66 eV/SiC between the phases, while the addition of RPA correlation reduces Pt by 40 GPa and the energy difference by 1 eV/SiC. Adding a kernel correction from RPAr yields a relatively small shift in the RPA results, reducing the transition pressure by ∼ 4 GPa and the energy difference by 0.07 eV/SiC. All of the kernel corrected methods yield essentially equivalent transition pressures and energy differences, though there is a 30 meV reduction in the energy difference from RPAr1 to HOT or the infinite-order method. Thus for this system it appears that the RPA results are quite acceptable for comparisons with experiment and the kernel introduces a relatively small shift. Though the predicted pressures appear too small in comparison to experiment, Ref. 97 argues that the stiffness of this material leads to a larger experimental transition pressure in order to overcome the transition barrier. By performing phonon calculations as a function of pressure, Ref. 97 reported that the ZB phase becomes unstable around 120 GPa, in better agreement with the experimental value.

23

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Looking at the performance of RPA renormalization, the first-order approximation already yields accurate transition pressure parameters compared to the infinite-order method. In contrast to the weak interactions in the ammonia and water dimers, the accuracy of the HOT approximation quantitatively reproduces the infinite-order volume of each phase at the phase transition, the transition pressure, and the energy difference between phases. This is somewhat expected given the degree with which the HOT approximation is able to reproduce the EOS fits for the solids discussed above. The excellent performance of the HOT approximation for this type of energy difference stems from the consistent overestimation of the bRPAr1 correlation in each phase by approximately 30 meV. Since the error has the same sign in both phases it drops out when taking the energy difference, thus the HOT approximation is capable of reproducing energy differences between phases of SiC with high fidelity. The same cancellation of errors can also be found for the interlayer binding energy of 2H MoS2 . Dominated by the weak interaction between layers, the interlayer binding energy is difficult to predict with bare semilocal functionals and typically a dispersion correction is needed for accurate results. 21,103 RPA has been used as a benchmark in these systems due to absence of high quality experimental data or other high level calculations, 21 yielding roughly 20 meV/˚ A2 for the binding energy of MoS2 . Our kernel-corrected calculations further confirm the accuracy of RPA for these systems since there is little difference between any of our RPAr results with rAPBE and the original RPA value. RPAr1 itself seems to capture all but 1 meV of the binding energy compared to the infinite-order method, and the addition of the HOT correction changes the RPAr1 result by very little. Just as for the SiC transition pressure, the error of the HOT approximation for mono- and bilayer exhibits the same sign and so the error for the difference drops out yielding high fidelity to the XACFD result. Thus the HOT approximation is likely capable of predicting accurate exfoliation and surface properties since these involve differences of bulk-like structures. Furthermore, the adsorption of molecules on surfaces should also be

24

ACS Paragon Plus Environment

Page 24 of 40

Page 25 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 4: SJEOS equation of state parameters and phase transition parameters for the ZB→RS transition in SiC. ∆E0 = E0 (RS) − E0 (ZB) at the equilibrium volume of each phase, Pt is the transition pressure, ∆V = (VtZB − VtRS )/VtZB is the percent reduction in volume at the transition, and Vt is the volume at the transition pressure computed for a given method. EXX yields a large estimate for the transition pressure, which is greatly reduced by the addition of RPA correlation. The addition of a kernel correction produces a relatively small shift from the RPA results. RPA RPAr1d Zincblende ˚ a (A) 4.390 4.363 4.380 4.379 B0 (GPa) 211 252 222 221 0 B0 3.91 3.61 3.86 3.93 Rocksalt a (˚ A) 4.071 3.980 4.049 4.047 B0 (GPa) 250 372 287 289 0 4.16 3.73 3.87 3.88 B0 ∆E0 (eV/SiC) 1.47 2.66 1.65 1.61 Pt (GPa) 66 115 74 71 ∆V (%) 18 19 18 18 ZB ˚3 Vt (A /SiC) 17.180 15.748 16.866 16.966 VtRS (˚ A3 /SiC) 14.076 12.796 13.840 13.897 PBE

a b c d

EXX

HOTd

XACFDd

Ref.

4.376 222 3.94

4.375 223 3.93

4.36a 224a 3.57a

4.043 292 3.87 1.58 70 18 16.987 13.906

4.042 293 3.87 1.58 70 18 16.981 13.903

4.04b 252b 4.26b – 100a 20a 15.685c 12.494c

Experimental data from Ref. 94 for ZB and the phase transition. Rocksalt LDA values from Ref. 102. Back calculated from Refs. 94 and 97. Kernel corrected methods used the rAPBE kernel.

25

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

accurately reproduced by the HOT approximation since both the surface and the molecule tend to be spin-unpolarized species and the HOT approximation is accurate in those cases.

4

Discussion and Conclusions

In order to accurately treat short-ranged correlation within the ACFD, an exchange-correlation kernel should be added to RPA. RPA renormalization provides a systematic route to compute such corrections, with the added benefit that both RPA and beyond RPA contributions to the correlation energy can be computed simultaneously. To recover the missing correlation beyond RPAr1, we developed the HOT approximation to account for second-order and higher contributions from RPA renormalization. Though the HOT approximation is no longer a systematic underestimate of the beyond RPA correlation energy as is RPAr1, it delivers total correlation energies typically within 1% of the infinite-order method. Such kernel-corrected calculations can be used for highly accurate predictions, or as an alternative benchmark for semilocal functionals in systems without other reference values. We also note that, in general, the RPAr1 and HOT corrections converge more rapidly than RPA with respect to k-mesh and cutoff choices, 57 and so often a kernel corrected calculation may be more well converged with reduced computational settings compared to a full RPA calculation in a plane-wave basis. Though the total correlation energy exhibits a small error, energy differences can be adversely affected by the non-systematic behavior of the HOT approximation. Atomization and cohesive energies, as well as the weak interactions in two molecular dimers, exhibit a weak cancellation of errors between a composite system and its fragments. The majority of this error stems from the error of the HOT approximation for the fragments such as in the open-shell oxygen and carbon atoms. We note that our current approximation ˆb = 1/2 is also a poor choice for strongly correlated systems, in addition to spin-polarized atoms. In spite of these challenges, the HOT approximation quantitatively reproduces the structural

26

ACS Paragon Plus Environment

Page 26 of 40

Page 27 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

parameters of several bulk solids, and the pressure induced phase transition of SiC from zincblende to rocksalt structures. We thus view the HOT approximation as an accurate beyond-RPA method for exploring potential energy surfaces and isogyric processes (processes that conserve the number of electron pairs), since for these situations the HOT approximation can be expected to faithfully reproduce the XACFD result, while simultaneously reducing the computational cost. For fragmentation processes such as atomization and cohesive energies, we can recommend the HOT approximation for the bulk-like or molecular species, but that the infinite-order (XACFD) approach be utilized for any open-shell atomic species to eliminate the deficiencies of the HOT approximation for spin-polarized atoms. For weak interactions, on the other hand, RPAr1 is likely sufficient for all species and the HOT correction is not needed. To remedy the challenges facing the HOT approximation, one could fit ˆb to minimize an error function, or by determining ˆb as a global function of spin-polarization to make up for the poor choice of approximating the bRPAr1 correlation potential from RPAr2. Further work is underway in these directions, as well as exploring the HOT approximation in combination with non-linear exchange-correlation kernels.

Acknowledgements We thank K. Burke and J. P. Perdew for useful discussions. J.E.B. and A.R. acknowledge the National Science Foundation (NSF) for financial support. This project was supported by the NSF Division of Materials Research under Grant No. DMR–1553022, and in part by the NSF through major research instrumentation Grant No. CNS–09–58854. Figures were created using matplotlib. 104

References (1) Hohenberg, P.; Kohn, W. Inhomogeneous electron gas. Phys. Rev. 1964, 136, B864. 27

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(2) Kohn, W.; Sham, L. J. Self-consistent equations including exchange and correlation effects. Phys. Rev. 1965, 140, A1133. (3) Langreth, D. C.; Perdew, J. P. Exchange-correlation energy of a metallic surface: Wave-vector analysis. Phys. Rev. B 1977, 15, 2884. (4) Heßelmann, A.; G¨orling, A. Random-phase approximation correlation methods for molecules and solids. Mol. Phys. 2011, 109, 24732500. (5) Eshuis, H.; Bates, J. E.; Furche, F. Electron correlation methods based on the random phase approximation. Theor. Chem. Acc. 2012, 131, 1084. (6) Ren, X.; Rinke, P.; Joas, C.; Scheffler, M. Random-phase approximation and its applications in computational chemistry and materials science. J. Mater. Sci. 2012, 47, 7447–7471. (7) Langreth, D. C.; Perdew, J. P. The exchange-correlation energy of a metallic surface. Solid State Commun. 1975, 17, 1425. (8) Furche, F. Molecular tests of the random phase approximation to the exchangecorrelation energy functional. Phys. Rev. B 2001, 64, 195120. (9) Furche, F.; Van Voorhis, T. Fluctuation-dissipation theorem density-functional theory. J. Chem. Phys. 2005, 122, 164106. (10) Paier, J.; Ren, X.; Rinke, P.; Scuseria, G.; Gr¨ uneis, A.; Kresse, G.; Scheffler, M. Assessment of correlation energies based on the random-phase approximation. New J. of Phys. 2012, 14, 043002. (11) Eshuis, H.; Yarkony, J.; Furche, F. Fast computation of molecular random phase approximation correlation energies using resolution of the identity and imaginary frequency integration. J. Chem. Phys. 2010, 132, 234114.

28

ACS Paragon Plus Environment

Page 28 of 40

Page 29 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(12) Del Ben, M.; Hutter, J.; VandeVondele, J. Electron Correlation in the Condensed Phase from a Resolution of Identity Approach Based on the Gaussian and Plane Waves Scheme. J. Chem. Theory Comput. 2013, 9, 2654–2671. (13) Moussa, J. E. Cubic-scaling algorithm and self-consistent field for the random-phase approximation with second-order screened exchange. J. Chem. Phys. 2014, 140, 014107. (14) Kaltak, M.; Klimeˇs, J.; Kresse, G. Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in Si. Phys. Rev. B 2014, 90, 054115. (15) K´allay, M. Linear-scaling implementation of the direct random-phase approximation. J. Chem. Phys. 2015, 142, 204105. (16) Wilhelm, J.; Seewald, P.; Del Ben, M.; Hutter, J. Large-Scale Cubic-Scaling Random Phase Approximation Correlation Energy Calculations Using a Gaussian Basis. J. Chem. Theory Comput. 2016, 12, 5851–5859. (17) Harl, J.; Kresse, G. Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory. Phys. Rev. B 2008, 77, 045136. (18) Harl, J.; Kresse, G. Accurate Bulk Properties from Approximate Many-Body Techniques. Phys. Rev. Lett. 2009, 103, 056401. ´ (19) Leb`egue, S.; Harl, J.; Gould, T.; Angy´ an, J. G.; Kresse, G.; Dobson, J. F. Cohesive properties and asymptotics of the dispersion interaction in graphite by the random phase approximation. Phys. Rev. Lett. 2010, 105, 196401. (20) Schimka, L.; Harl, J.; Stroppa, A.; Gr¨ uneis, A.; Marsman, M.; Mittendorfer, F.; Kresse, G. Accurate surface and adsorption energies from many-body perturbation theory. Nat. Mater. 2010, 9, 741–744.

29

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(21) Bj¨orkman, T.; Gulans, A.; Krasheninnikov, A. V.; Nieminen, R. M. van der Waals Bonding in Layered Compounds from Advanced Density-Functional First-Principles Calculations. Phys. Rev. Lett. 2012, 108, 235502. (22) Xiao, B.; Sun, J.; Ruzsinszky, A.; Feng, J.; Perdew, J. P. Structural phase transitions in Si and SiO 2 crystals via the random phase approximation. Phys. Rev. B 2012, 86, 094109. (23) Yan, J.; Nørskov, J. K. Calculated formation and reaction energies of 3d transition metal oxides using a hierarchy of exchange-correlation functionals. Phys. Rev. B 2013, 88, 245204. (24) Karlick´ y, F.; Lazar, P.; Dubeck´ y, M.; Otyepka, M. Random Phase Approximation in Surface Chemistry: Water Splitting on Iron. J. Chem. Theory Comput. 2013, 9, 3670–3676. (25) Peng, H.; Lany, S. Polymorphic energy ordering of MgO, ZnO, GaN, and MnO within the random phase approximation. Phys. Rev. B 2013, 87, 174113. (26) Schimka, L.; Gaudoin, R.; Klimeˇs, J.; Marsman, M.; Kresse, G. Lattice constants and cohesive energies of alkali, alkaline-earth, and transition metals: Random phase approximation and density functional theory results. Phys. Rev. B 2013, 87, 214102. (27) Olsen, T.; Thygesen, K. S. Random phase approximation applied to solids, molecules, and graphene-metal interfaces: From van der Waals to covalent bonding. Phys. Rev. B 2013, 87, 075111. (28) Burow, A. M.; Bates, J. E.; Furche, F.; Eshuis, H. Analytical First-Order Molecular Properties and Forces within the Adiabatic Connection Random Phase Approximation. J. Chem. Theory Comput. 2014, 10, 180–194.

30

ACS Paragon Plus Environment

Page 30 of 40

Page 31 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(29) Bao, J. L.; Yu, H. S.; Duanmu, K.; Makeev, M. A.; Xu, X.; Truhlar, D. G. Density Functional Theory of the Water Splitting Reaction on Fe(0): Comparison of Local and Nonlocal Correlation Functionals. ACS Catal. 2015, 5, 2070–2080. (30) Mezei, P. D.; Csonka, G. I.; K´allay, M. Accurate DielsAlder Reaction Energies from Efficient Density Functional Calculations. J. Chem. Theory Comput. 2015, 11, 2879– 2888. (31) Waitt, C.; Ferrara, N. M.; Eshuis, H. Thermochemistry and Geometries for TransitionMetal Chemistry from the Random Phase Approximation. J. Chem. Theory. Comput. 2016, 12, 5350–5360. (32) Bates, J. E.; Mezei, P. D.; Csonka, G. I.; Sun, J.; Ruzsinszky, A. Reference Determinant Dependence of the Random Phase Approximation in 3d Transition Metal Chemistry. J. Chem. Theory Comput. 2017, 13, 100–109. (33) Bohm, D.; Pines, D. A collective description of electron interactions: II. Collective vs individual particle aspects of the interactions. Phys. Rev. 1952, 85, 338. (34) Fuchs, M.; Niquet, Y. M.; Gonze, X.; Burke, K. Describing static correlation in bond dissociation by Kohn–Sham density functional theory. J. Chem. Phys. 2005, 122, 094116. (35) Jiang, H.; Engel, E. Random-phase-approximation-based correlation energy functionals: Benchmark results for atoms. J. Chem. Phys. 2007, 127, 184108. ´ (36) Mussard, B.; Reinhardt, P.; Angy´ an, J. G.; Toulouse, J. Spin-unrestricted randomphase approximation with range separation: Benchmark on atomization energies and reaction barrier heights. J. Chem. Phys. 2015, 142, 154123. (37) Lein, M.; Gross, E. K. U.; Perdew, J. P. Electron correlation energies from scaled

31

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

exchange-correlation kernels: Importance of spatial versus temporal nonlocality. Phys. Rev. B 2000, 61, 13431. (38) Heßelmann, A.; G¨orling, A. Random phase approximation correlation energies with exact KohnSham exchange. Mol. Phys. 2010, 108, 359–372. ´ (39) Angy´ an, J. G.; Liu, R.; Toulouse, J.; Jansen, G. Correlation Energy Expressions from the Adiabatic-Connection Fluctuation-Dissipation Theorem Approach. J. Chem. Theory Comput. 2011, 7, 3116. (40) Bates, J. E.; Furche, F. Communication: Random phase approximation renormalized many-body perturbation theory. J. Chem. Phys. 2013, 139, 171103. (41) Lotrich, V.; Bartlett, R. J. External coupled-cluster perturbation theory: Description and application to weakly interaction dimers. Corrections to the random phase approximation. J. Chem. Phys. 2011, 134, 184108. (42) Scuseria, G. E.; Henderson, T. M.; Sorensen, D. C. The ground state correlation energy of the random phase approximation from a ring coupled cluster doubles approach. J. Chem. Phys. 2008, 129, 231101. ´ (43) Jansen, G.; Liu, R.-F.; Angy´ an, J. G. On the equivalence of ring-coupled cluster and adiabatic connection fluctuation-dissipation theorem random phase approximation correlation energy expressions. J. Chem. Phys. 2010, 133, 154106. (44) Gr¨ uneis, A.; Marsman, M.; Harl, J.; Schimka, L.; Kresse, G. Making the random phase approximation to electronic correlation accurate. J. Chem. Phys. 2009, 131, 154115. (45) Henderson, T. M.; Scuseria, G. E. The connection between self-interaction and static correlation: a random phase approximation perspective. Mol. Phys. 2010, 108, 2511. (46) Bates, J. E.; Laricchia, S.; Ruzsinszky, A. A Non-Local, Energy-Optimized Kernel:

32

ACS Paragon Plus Environment

Page 32 of 40

Page 33 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Recovering Second-Order Exchange in the Homogeneous Electron Gas. Phys. Rev. B 2016, 93, 045119. (47) Dixit, A.; Claudot, J.; Leb´egue, S.; Rocca, D. Improving the Efficiency of Beyond-RPA Methods within the Dielectric Matrix Formulation: Algorithms and Applications to the A24 and S22 Test Sets. J. Chem. Theory Comput. 2017, 13, 5432–5442. (48) Corradini, M.; Del Sole, R.; Onida, G.; Palummo, M. Analytical expressions for the local-field factor G(q) and the exchange-correlation kernel Kxc (r) of the homogeneous electron gas. Phys. Rev. B 1998, 57, 14569. (49) Constantin, L. A.; Pitarke, J. M. Simple dynamic exchange-correlation kernel of a uniform electron gas. Phys. Rev. B 2007, 75, 245127. (50) Hesselmann, A.; G¨orling, A. Correct description of the bond dissociation limit without breaking spin symmetry by a random-phase-approximation correlation functional. Phys. Rev. Lett. 2011, 106, 093001. (51) Gould, T. Communication: Beyond the random phase approximation on the cheap: Improved correlation energies with the efficient radial exchange hole kernel. J. Chem. Phys. 2012, 137, 111101. (52) Olsen, T.; Thygesen, K. S. Extending the random-phase approximation for electronic correlation energies: The renormalized adiabatic local density approximation. Phys. Rev. B 2012, 86, 081103(R). (53) Olsen, T.; Thygesen, K. S. Accurate Ground-State Energies of Solids and Molecules from Time-Dependent Density-Functional Theory. Phys. Rev. Lett. 2014, 112, 203001. (54) Lu, D. Evaluation of model exchange-correlation kernels in the adiabatic connection fluctuation-dissipation theorem for inhomogeneous systems. J. Chem. Phys. 2014, 140, 18A520. 33

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(55) Erhard, J.; Bleiziffer, P.; G¨orling, A. Power Series Approximation for the Correlation Kernel Leading to Kohn-Sham Methods Combining Accuracy, Computational Efficiency, and General Applicability. Phys. Rev. Lett. 2016, 117, 143002. (56) Olsen, T.; Thygesen, K. S. Beyond the random phase approximation: Improved description of short-range correlation by a renormalized adiabatic local density approximation. Phys. Rev. B 2013, 88, 115131. (57) Patrick, C. E.; Thygesen, K. S. Adiabatic-connection fluctuation-dissipation DFT for the structural properties of solidsThe renormalized ALDA and electron gas kernels. J. Chem. Phys. 2015, 143, 102802. (58) Jauho, T. S.; Olsen, T.; Bligaard, T.; Thygesen, K. S. Improved description of metal oxide stability: Beyond the random phase approximation with renormalized kernels. Phys. Rev. B 2015, 92, 115140. (59) Bates, J. E.; Sensenig, J.; Ruzsinszky, A. Convergence behavior of the random phase approximation renormalized correlation energy. Phys. Rev. B 2017, 95, 195158. (60) Colonna, N.; Hellgren, M.; de Gironcoli, S. Correlation energy within exact-exchange adiabatic connection fluctuation-dissipation theory: Systematic development and simple approximations. Phys. Rev. B 2014, 90, 125150. (61) G¨orling, A.; Levy, M. Correlation-energy functional and its high-density limit obtained from a coupling-constant perturbation expansion. Phys. Rev. B 1993, 47, 13105. (62) G¨orling, A.; Levy, M. Exact Kohn-Sham scheme based perturbation theory. Phys. Rev. A 1994, 50, 196. (63) K. Burke et al., The ABC of DFT ; 2007. (64) Petersilka, M.; Gossmann, U. J.; Gross, E. K. U. Excitation energies from timedependent density-functional theory. Phys. Rev. Lett. 1996, 76, 1212. 34

ACS Paragon Plus Environment

Page 34 of 40

Page 35 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(65) Laidig, W. D.; Fitzgerald, G.; Bartlett, R. J. Is fifth-order MBPT enough? Chem. Phys. Lett. 1985, 113, 151 – 158. (66) Cremer, D.; He, Z. Sixth-Order Møller- Plesset Perturbation Theory On the Convergence of the MP n Series. J. Phys. Chem. 1996, 100, 6173–6188. (67) Seidl, M.; Perdew, J. P.; Kurth, S. Simulation of All-Order Density-Functional Perturbation Theory, Using the Second Order and the Strong-Correlation Limit. Phys. Rev. Lett. 2000, 84, 5070–5073. (68) Manby, F.; Knowles, P. J. An exchange functional for accurate virtual orbital energies. J. Chem. Phys. 2000, 112, 7002. (69) Burke, K.; Ernzerhof, M.; Perdew, J. P. The adiabatic connection method: a nonempirical hybrid. Chem. Phys. Lett. 1997, 265, 115 – 120. (70) This was graciously pointed out by one of the Referees. (71) See the supplemental material available free of charge via http://pubs.acs.org for a compressed file containing some raw data, as well as a pdf describing the contents of the compressed file and computational details. (72) Ernzerhof, M.; Perdew, J. P.; Burke, K. Coupling-constant dependence of atomization energies. Int. J. Quant. Chem. 1997, 64, 285. (73) Grundei, M. M. J.; Burow, A. M. Random Phase Approximation for Periodic Systems Employing Direct Coulomb Lattice Summation. J. Chem. Theory Comput. 2017, 13, 1159–1175. (74) Walter, M.; H¨akkinen, H.; Lehtovaara, L.; Puska, M.; Enkovaara, J.; Rostgaard, C.; Mortensen, J. J. Time-dependent density-functional theory in the projector augmented-wave method. J. Chem. Phys. 2008, 128, 244101.

35

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(75) Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Real-space grid implementation of the projector augmented wave method. Phys. Rev. B 2005, 71, 035109. (76) Bahn, S. R.; Jacobsen, K. W. An object-oriented scripting interface to a legacy electronic structure code. Comput. Sci. Eng. 2002, 4, 56–66. (77) Trevisanutto, P. E.; Terentjevs, A.; Constantin, L. A.; Olevano, V.; Sala, F. D. Optical spectra of solids obtained by time-dependent density functional theory with the jellium-with-gap-model exchange-correlation kernel. Phys. Rev. B 2013, 87, 205143. (78) Bl¨ochl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953–17979. (79) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188–5192. (80) Yan, J.; Mortensen, J. J.; Jacobsen, K. W.; Thygesen, K. S. Linear density response function in the projector augmented wave method: Applications to solids, surfaces, and interfaces. Phys. Rev. B 2011, 83, 245122. (81) Sundararaman, R.; Arias, T. A. Regularization of the Coulomb singularity in exact exchange by Wigner-Seitz truncated interactions: Towards chemical accuracy in nontrivial systems. Phys. Rev. B 2013, 87, 165122. (82) Alchagirov, A. B.; Perdew, J. P.; Boettger, J. C.; Albers, R. C.; Fiolhais, C. Energy and pressure versus volume: Equations of state motivated by the stabilized jellium model. Phys. Rev. B 2001, 63, 224115. (83) Harl, J.; Schimka, L.; Kresse, G. Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids. Phys. Rev. B 2010, 81, 115126.

36

ACS Paragon Plus Environment

Page 36 of 40

Page 37 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(84) Patrick, C. E.; Thygesen, K. S. Hubbard-U corrected Hamiltonians for non-selfconsistent random-phase approximation total-energy calculations: A study of ZnS, TiO2 , and NiO. Phys. Rev. B 2016, 93, 035133. (85) Jurecka, P.; Sponer, J.; Cerny, J.; Hobza, P. Benchmark database of accurate (MP2 and CCSD(T) complete basis set limit) interaction energies of small model complexes, DNA base pairs, and amino acid pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985. (86) Takatani, T.; Hohenstein, E. G.; Malagoli, M.; Marshall, M. S.; Sherrill, C. D. Basis set consistent revision of the S22 test set of noncovalent interaction energies. J. Chem. Phys. 2010, 132, 144104. (87) Ren, X.; Rinke, P.; Scuseria, G. E.; Scheffler, M. Renormalized second-order perturbation theory for the electron correlation energy: Concept, implementation, and benchmarks. Phys. Rev. B 2013, 88, 035120. (88) Pantelides, C. C.; Adjiman, C. S.; Kazantsev, A. V. General Computational Algorithms for Ab Initio Crystal Structure Prediction for Organic Molecules. Top. Curr. Chem. 2014, 345, 25. (89) Price, S. L. Predicting crystal structures of organic compounds. Chem. Soc. Rev. 2014, 43, 2098–2111. (90) Szalewicz, K. Determination of Structure and Properties of Molecular Crystals from First Principles. Acc. Chem. Res. 2014, 47, 3266–3274. (91) Cui, Z.-H.; Wu, F.; Jiang, H. First-principles study of relative stability of rutile and anatase TiO 2 using the random phase approximation. Phys. Chem. Chem. Phys. 2016, 18, 29914–29922. (92) Brandenburg, J. G.; Bates, J. E.; Sun, J.; Perdew, J. P. Benchmark tests of a strongly

37

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

constrained semilocal functional with a long-range dispersion correction. Phys. Rev. B 2016, 94, 115144. (93) Whittleton, S. R.; Otero-de-la Roza, A.; Johnson, E. R. The exchange-hole dipole dispersion model for accurate energy ranking in molecular crystal structure prediction. J. Chem. Theory Comput. 2017, 13, 441. (94) Yoshida, M.; Onodera, A.; Ueno, M.; Takemura, K.; Shimomura, O. Pressure-induced phase transition in SiC. Phys. Rev. B 1993, 48, 10587–10590. (95) Mujica, A.; Rubio, A.; Mu˜ noz, A.; Needs, R. J. High-pressure phases of group-IV, III–V, and II–VI compounds. Rev. Mod. Phys. 2003, 75, 863–912. (96) Karch, K.; Pavone, P.; Windl, W.; Sch¨ utt, O.; Strauch, D. Ab initio calculation of structural and lattice-dynamical properties of silicon carbide. Phys. Rev. B 1994, 50, 17054–17063. (97) Lu, Y.-P.; He, D.-W.; Zhu, J.; Yang, X.-D. First-principles study of pressure-induced phase transition in silicon carbide. Physica B: Condensed Matter 2008, 403, 3543 – 3546. (98) Sarasamak, K.; Kulkarni, A. J.; Zhou, M.; Limpijumnong, S. Stability of wurtzite, unbuckled wurtzite, and rocksalt phases of SiC, GaN, InN, ZnO, and CdSe under loading of different triaxialities. Phys. Rev. B 2008, 77, 024104. (99) Catti, M. Orthorhombic Intermediate State in the Zinc Blende to Rocksalt Transformation Path of SiC at High Pressure. Phys. Rev. Lett. 2001, 87, 035504. (100) Sengupta, N.; Bates, J. E.; Ruzsinszky, A. From Semilocal Density Functionals to Random Phase Approximation Renormalized Perturbation Theory: A Methodological Assessment of Structural Phase Transitions. submitted 2018,

38

ACS Paragon Plus Environment

Page 38 of 40

Page 39 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(101) Karch, K.; Bechstedt, F.; Pavone, P.; Strauch, D. Pressure-dependent properties of SiC polytypes. Phys. Rev. B 1996, 53, 13400–13413. (102) Kohn, W.; Becke, A. D.; Parr, R. G. Density functional theory of electronic structure. J. Phys. Chem. 1996, 100, 12974. (103) Peng, H.; Yang, Z.-H.; Perdew, J. P.; Sun, J. Versatile van der Waals Density Functional Based on a Meta-Generalized Gradient Approximation. Phys. Rev. X 2016, 6, 041005. (104) Hunter, J. D. Matplotlib: A 2D graphics environment. Comput. Sci. Eng. 2007, 9, 90–95.

39

ACS Paragon Plus Environment

RPA Renormalized (RPAr) Correlation Energies

Journal of Chemical 3.2 Theory and Computation Page 40 RPAr1 40 of RPAr3 RPAr2

HOT

1 2 3 4 5 6 7

EcbRPA (eV)

3.0

2.8

2.6

fxc

ACS Paragon2.4Plus Environment fxc

2.2

CO (1 )

Si2 - diamond