Random phase approximation with exchange for an accurate description of crystalline polymorphism

We determine the correlation energy of BN, SiO$_2$ and ice polymorphs employing a recently developed RPAx (random phase approximation with exchange) approach. The RPAx provides larger and more accurate polarizabilities as compared to the RPA, and captures effects of anisotropy. In turn, the correlation energy, defined as an integral over the density-density response function, gives improved binding energies without the need for error cancellation. Here, we demonstrate that these features are crucial for predicting the relative energies between low- and high-pressure polymorphs of different coordination number as, e.g., between $\alpha$-quartz and stishovite in SiO$_2$, or layered and cubic BN. Furthermore, a reliable (H$_2$O)$_2$ potential energy surface is obtained, necessary for describing the various phases of ice. The RPAx gives results comparable to other high-level methods such as coupled cluster and quantum Monte Carlo, also in cases where the RPA breaks down. Although higher computational cost than RPA we observe a faster convergence with respect to the number of eigenvalues in the response function.


I. INTRODUCTION
The development of an accurate, yet computationally efficient, electronic structure approach able to treat electron correlation in solids remains an important task in condensed matter physics. Density functional theory (DFT) provides a rigorous and computationally appealing framework [1,2] that has been widely successful, but available approximations still suffer a number of drawbacks. The van der Waals (vdW) interactions, ubiquitous in systems ranging from layered materials to superconducting hydrides, are challenging to include [3], and most exchange-correlation (xc) functionals do not achieve the accuracy needed to satisfactorily describe the energetics of phase transformations and polymorphism, even in common systems like water [4] and silica [5].
The random phase approximation (RPA) represents the highest level of sophistication currently applied in materials science [6][7][8]. Being based on the exact adiabatic-connection fluctuation-dissipation (ACFD) formula [9,10], first-order Hartree-Fock exchange is exactly incorporated and vdW forces are seamlessly built in, providing an overall improved accuracy. However, the performance of the RPA strongly relies on error cancellation, leading to unpredictable errors when dealing with crystals of low symmetry [11,12]. Futhermore, vdW forces are underestimated by an average of 20% [13][14][15]. The origin of these errors can, however, easily be traced to the approximate Hartree response function used for constructing the correlation energy.
The combination of the ACFD formula and timedependent (TD) DFT has opened a path for systematic improvements of the RPA [16,17]. Within TDDFT, xc effects are incorporated via the xc kernel that can be defined as the second density variation of the xc action functional [18][19][20]. Starting from the RPA, i.e., the Hartree kernel, further refinements then naturally proceed via the local density approximation (LDA) and generalized gradient approximations (GGAs). However, first results found that these approximations can worsen the performance with respect to the RPA [21]. A renormalized LDA kernel that introduces spacial nonlocality was, therefore, formulated in Ref. [22]. This kernel, as well as a renormalized PBE (Perdew-Burke-Ernzerhof) kernel, have shown to improve the RPA correlation energy, leading to a better performance in several cases [12,23,24].
Within an approach that combines many-body perturbation theory (MBPT) and TDDFT the first step beyond RPA is to include the full Fock exchange term in the density response function via the nonlocal and frequency dependent exact-exchange (EXX) kernel. This generates the random phase approximation with exchange (RPAx), that exactly includes the second order exchange diagram as well as higher order exchange effects, in addition to the RPA ring series of diagrams [25][26][27][28][29]. Various partial resummations of the RPAx correlation terms are thus equal to other advanced expressions for the correlation energy defined within MBPT such as, e.g., SOSEX (second order screened exchange) [15,[30][31][32][33][34]. First tests on atoms and molecules have shown that not only accurate correlation energies are obtained with RPAx [15,35], but also xc potentials, i.e., electronic densities, and polarisabilities [28,36]. Furthermore, studies on the homogeneous electron gas, the simplest model of a metal, have shown promising results [34].
In this work, we extend the scope of applications to solids and show that it is possible to reach results of similar quality as for molecules. We calculate the relative energies between various crystalline phases within the BN, SiO 2 and ice polymorphs. In particular, we investigate cases where the RPA has shown insufficient, as, e.g., for the α-quartz-stishovite energy difference in SiO 2 , for which the change in coordination number necessitate accurate correlation energies [5,11,12]. A similar problem occurs in cubic and layered BN, inverting their stability order [37,38]. In both SiO 2 and BN the vdW forces play an important role [3,5,39]. Here we investigate the effect of the EXX kernel on the binding energy of arXiv:2101.07052v2 [cond-mat.mtrl-sci] 31 Aug 2021 layered BN as well as the energy difference between αquartz and cristobalite. A study of purely vdW bonded solids, such as Ar and Kr, confirms that the vdW bond is well described with RPAx, yielding an accuracy similar to that found for molecular dimers. Finally, we perform a detailed analysis of the water dimer potential energy surface and ice polymorphism. The delicate balance between Pauli repulsion, static and dynamic correlation has made water a long-standing problem for electronic structure methods.
Our results on this broad and challenging set of systems show promise for future applications. The RPAx stays within the simple computational framework of the RPA, but has an accuracy comparable with more sophisticated methods such as quantum Monte Carlo (QMC) and coupled cluster (CC).

II. RPA WITH EXCHANGE
In this section we summarize the equations of the RPAx approximation for periodic solids. These have been introduced previously [27,28] but are here adapted to the specific implementation employed in this work [15,34,35].
The RPAx approximation for the correlation energy is based on the ACFD (adiabatic connection fluctuationdissipation) formula, in which the exact correlation energy is expressed in terms of the linear density response function, χ λ , multiplied by the Coulomb interaction, v, and integrated over the frequency and the interaction strength λ The underlying Hamiltonian is given by where T is the kinetic energy operator, W is the electronelectron interaction operator, and V λ is a one-body potential operator composed of the external potential and a fictitious potential that fix the density at full interaction strength for every λ. With this definition, the system at λ = 0 already produces the exact density, n, via the exact Kohn-Sham (KS) xc potential, v xc . The KS independent-particle response function, subtracted in Eq.
(1), is denoted χ s . Within TDDFT, χ λ becomes a functional of the ground-state density, n 0 , via the xc kernel, f λ xc = δv λ xc /δn| n0 . This can be expressed in a Dyson-like equation similar to the RPA (or Hartree) equation (f xc = 0). A systematic approach for generating xc kernels based on many-body perturbation theory was introduced in Ref. 20. Within this approach the simplest approximation corresponds to including the full Fock exchange term via a time-dependent optimized effective potential scheme. This results in the EXX kernel, f x , defined by the following equation The right hand side can be interpreted diagrammatically in terms of KS Green functions, being the sum of the first order vertex diagram, R V (ω), and the first order self-energy correction diagrams R Σ (ω) (see Fig. 1). Although f x is first order in the Coulomb interaction, the correlation energy, via Eq. (1), contains all orders in v. This approximation to the correlation energy is what we call RPAx and is the same as the one used in Ref. [29]. Other variants based on a dielectric formulation [40], or Hartree-Fock orbitals and range-separation have also been proposed [41]. We emphasize that the RPAx, as defined here, is based on TDDFT and KS orbitals, thus falling in the category of orbital functionals in DFT (i.e. functionals depending on the (un)occupied KS orbitals and eigenvalues). The total energy functional can then be minimized via a local KS potential in the usual way [28]. The equation for the EXX kernel, Eq. (4), is usually expressed in terms of KS orbitals (occupied and unoccupied), ϕ n , and KS eigenvalues, ε n . In Ref. 34 it was shown that the same expression can be reformulated in terms of occupied KS orbitals and their responses. For periodic systems this expression reads where with Σ k (v x ) the nonlocal(local) Fock exchange potential and V α q a local potential. In this form one can apply the linear response techniques of density functional perturbation theory to determine Eq. (4) and, hence, avoid the generation of unoccupied KS states.
By solving the generalized eigenvalue problem where R = R V + R Σ + χ s vχ s , the correlation energy can be computed from where N ν is the total number of eigenvalues and N q is the total number of q-points in the Brillouin zone. The latter can be reduced exploiting crystal symmetries. The λ-integration is done analytically, while the frequency integral has to be carried out numerically. We immediately see that this expression is reduced to RPA by setting R V = 0 and R Σ = 0 in Eq. (9). Equations (5)- (10) have been implemented within the plane-wave and pseudopotential framework, but have so far only been applied to molecular systems [15,34,35]. Binding energies have been shown to converge quickly with respect to the number of eigenvalues (≈ 10 times the number of electrons (N e )). Due to the additional kpoint sum in R V , RPAx scales as N 3 k with respect to the number of k-points, which can be compared to N 2 k with RPA. On the other hand, as we will see later, the RPAx converges, in many cases, faster than RPA with respect to N ν .
The underlying many-body formulation of the f xkernel allows also a variant of SOSEX [30][31][32][33] to be generated within the present formalism, making use of an alternative resummation of Dyson's equation (Eq. (3)), that only includes a subset of terms [15]. In the following we will present both RPAx and SOSEX results.
Since a fully self-consistent implementation is still not developed, all calculations are done on top of a PBE ground-state [43], unless otherwise stated. The electronion interaction is treated with standard norm-conserving pseudopotentials [44][45][46].

III. NOBLE GAS SOLIDS
The ACFD formula for the correlation energy provides a natural starting-point for including the long-range vdW force, i.e., the attractive force between two charge neutral, closed-shell atoms A and B due to fluctuating dipoles present in a correlated wave function. At large separation R, the interaction energy can, to leading order, be written as [47] where is the isotropic dynamical polarizability of atom A, and C 6 is the vdW coefficient. In Ref. 48 Cohesive energy (eV) as a function of nearest neighbour distance (Å) in solid Ar (a) and Kr (b). The correlation energy (Ec-RPAx) is also shown, together with a fitted y = k/x 6 curve. Green and blue squares correspond to CCSD(T) results [42].
that the RPA correlation energy exactly reproduces Eq. (11), but with the exact polarizability replaced by the RPA polarizability. However, as shown in several previous works, the RPA C 6 -coefficients are underestimated (up to 40% in highly polarizable atoms) [28,33,[49][50][51][52]. For noble gas atoms the error lies at about 17%, which induces an error of approximately 30% in the binding energy of dimers [35,49]. Thus, while the RPA captures the vdW forces in a qualitatively correct way, the quantitative errors are large.
Including the xc kernel, it is not straightforward to derive an expression similar to Eq. 11. In fact, in Ref. [53] it was shown that semilocal kernels produce an additional 1/R 6 -term. It remains to be proven analytically whether the nonlocal RPAx reduces to Eq. 11.
In Fig. 2 we have plotted the cohesive energy of solid Ar and Kr as a function of nearest neighbour (NN) distance. In RPA and RPAx a reduced 3×3×3 shifted k-point grid is used for the correlation energy giving results within 4 meV of the fully converged values (see Appendix C).
In the RPA, the cohesive energy is underestimated by 25%. Our results thus differ somewhat from those obtained in Ref. [54], but are more consistent with the results found for dimers (see Refs. [35,49]). As expected, based on the accurate C 6 -coefficients and dimer binding energies, we find RPAx to be in good agreement with both CCSD(T) [42] and experimental results [55]. SO-SEX gives only slightly larger binding energies.
In Fig. 2 we also display the correlation energy (E c -RPAx) as a function of NN distance, exhibiting a perfect 1/R 6 behaviour. Given that Ar and Kr have 12 NN the fitted coefficients 67*12 and 135*12 (atomic units) compare well with the dimer C 6 -coefficients of 64 and 130, respectively. Similar accuracy for the C 6 -coefficients are found with other variants of RPAx [33,52,56,57] and the renormalized LDA kernel [51]. We note that the SO-SEX defined in Ref. [33] is different from the one defined here since it does not contain self-energy corrections (here included via R Σ in Eq. (4) [15]). The self-energy corrections are similar to the 'single excitations' in Ref. [33].

IV. BORON NITRIDE AND SILICA
Having established that vdW bonded solids are accurately reproduced with RPAx we now turn to systems were both vdW forces and covalent/ionic bonding are important. The polymorphs of BN consist of the highpressure cubic (c) and wurtzite (w) structures with fourfold coordinated B/N atoms, and several low-pressure layered structures, such as hexagonal (h) and rhombohedral (r) BN, which are threefold coordinated (see Fig.  3). The density, symmetry and ionic character of the B-N bond change from one phase to the other, with the largest difference between fourfold and threefold coordinated crystals. In these cases, functionals that largely rely on error cancellation become unreliable, producing substantial errors in energy differences that can invert the stability order. Although it is well known that the layered polymorphs need a good description of the vdW forces [39], this is not the only missing component. Inaccuracies in describing the covalent bond will show up if errors fail to cancel.
In Fig. 4 we plot the interlayer binding energy of r-BN as a function of interlayer distance, keeping the intralayer B-N distance fixed to the experimental value. In general, the interlayer bonding is difficult to model with vdW corrected DFT. Binding energies are typically largely overestimated, often twice as large as RPA values [3,58]. On the other hand, RPA underestimates the vdW forces, but it is not clear to what extent in layered materials. We have verified that we obtain the same RPA result as in Ref. [3], taking into account that h-BN is 5 meV lower in energy than r-BN [37]. The RPAx increases the binding energy by 7 meV (or 9%). This correction is smaller than the correction in purely vdW bonded solids such as Ar and Kr analyzed above. The RPA thus performs better than expected, at least for the BN layered system. We note that, except for some differences in the dissociation tail, the SOSEX and RPAx give essentially the same result.
In the lower panel of Fig. 5 we present results for the energy difference between c-BN and w-BN, and between c-BN and r-BN. c-BN is the reference level at zero in the figure. Calculations are done on experimental structures [59]. We compare our RPA and RPAx results to recent CCSD and CCSD(T) results [38], to LDA and PBE0 (i.e. a hybrid functional with 25% of Fock exchange). We see a striking difference in the performance on the different polymorphs. Both the wurtzite and cubic structures are high-pressure phases with the same B-N coordination and similar volume. This is thus the ideal situation for approximations to benefit from error cancellation. Indeed, almost all approximations give the same value within a couple of meV. The situation is very different when looking at the r-BN and c-BN energy difference. Here, the volume differ by 35% and the B-N bond change character. LDA is seen to hugely overestimate the energy difference (105 meV) and to predict c-BN to be lowest in energy. Experimental results are mixed but CC results revert the energy order (see Ref. [38] and experimental references therein [60,61]). Additionally, PBE0 favours this energy order, predicting r-BN to be lowest in energy, although with a too large energy difference (-67 meV). The RPA, that is known to rely on error cancellation, predicts the same energy order as LDA, but with a smaller energy difference (25 meV). We note that our result is in good agreement with the RPA calculation in Ref. [37]. The RPAx, known to give accurate absolute energies, and not only energy differences, gives a result very close to the CCSD result with an energy difference of -15 meV, and thus reverts the stability order. SOSEX gives a slightly larger energy difference (-23 meV). In Ref. [12] it was shown that the renormalized PBE kernel also provides an improvement upon RPA (-48-(-38) meV). We note that the SOSEX defined here is different from the ACSOSEX defined in Refs. [12,62]. The latter is based on the renormalized PBE kernel whereas the one here is based on many-body perturbation theory [15].
The difference between RPA and RPAx cannot solely be explained by the difference in how they describe the vdW forces. This suggests that the RPAx captures the covalent B-N bond better than RPA. Our result provides further confirmation that the correct order places the layered low-density structures at a slightly lower energy than the cubic form, even at zero temperature. The difference in zero-point (vibrational) energy (ZPE), has been estimated to -10 meV [37]. We also mention that we tried various DFT functionals employing a vdW correction, again with very mixed results.
In the upper panel of Fig. 5 we show the convergence of the correlation energy difference between c-BN and r-BN (∆E r−c c ) with respect to N ν /N e . We note that when comparing the correlation energy of two systems the same number of eigenvalues per electron has to be used [15]. Since r-BN and c-BN have the same number of electrons in the unit cell this implies the same number of eigenvalues. w-BN has twice the number of electrons and, therefore, needs twice the number of eigenvalues. However, while the correlation energy difference between w-BN and c-BN converges quickly with respect to the number of eigenvalues (< 10 × N e ) the RPA exhibits a slow convergence for the r-BN-c-BN energy difference. We used up to 512 eigenvalues in order to produce results within 2 meV (marked by the dashed horizontal lines in Fig. 5), which is more than 60 times the number of electrons and hence 6 times more than usual. The problem of not cancelling errors is thus visible also in this technical aspect of our approach. The absence of this problem is remarkably evident in the convergence   of the RPAx correlation energy difference. In RPAx we needed < 64 eigenvalues, which is similar to the number used for other systems. Further convergence studies can be found in Appendix B.
A similar effect, as studied above in BN, has previously been observed in silica (SiO 2 ) polymorphs [5,11,12]. Silica, similar to BN, has polar covalent bonds and polymorphs of different coordination number. In α-quartz (space group P3 1 21), the most stable crystal at ambient pressure, silicon is fourfold coordinated by oxygen, forming a tetrahedral unit. These units are connected via a flexible Si-O-Si bond, sensitive to vdW forces (see Fig. 3). Cristobalite (space group P4 1 2 1 2) has a similar local structure and only 5% larger volume. The energy difference between α-quartz and cristobalite is found to be approximately 20 meV in experiment [64,65]. When the ZPE contribution is subtracted the energy difference doubles [5]. While LDA predicts the correct energy order PBE erroneously predicts cristobalite to be lower in energy due to missing vdW forces [5]. A similar effect was recently observed in the borate (B 2 O 3 ) polymorphs [66]. Stishovite (space group P4 2 /mnm) is a high-pressure phase with 37% higher density, having sixfold coordinated Si atoms. The energy difference between stishovite and α-quartz is experimentally determined to 0.51-0.54 eV (see Ref. [67] and experimental references therein [68,69]). LDA underestimates this value by more than 0.4 eV, whereas PBE here predicts the correct en- ergy difference [67].
In the lower panels of Fig. 6 we illustrate the results for the α-quartz-cristobalite (red) and α-quartzstishovite (yellow) energy differences. We optimized all structures with PBE plus a vdW correction at the Tkatchenko-Scheffler (TS) level [70]. In the case of αquartz-cristobalite, LDA gives a value of 33 meV and PBE0 a value of 12 meV. As discussed in Ref. [5] a vdW correction improves this result. RPA gives a value of 31 meV, which is underestimated by 15 meV according to QMC calculations. As the Si-O bond remains similar in α-quartz and cristobalite, we interpret this relatively small error coming from the underestimated vdW forces with RPA. RPAx and SOSEX bring the result in good agreement with the QMC value of 43 meV.
The situation is different in the case of stishovite which has sixfold coordinated Si atoms. Here the variations between the different methods are an order of magnitude larger. Subtracting the ZPE, calculated in Ref. [5], QMC predicts an energy difference of 0.51 eV [63], similar to the PBE0 result (0.52 eV) and experiment. With RPAx we find a value of 0.49 eV which is 80 meV larger than the RPA result (0.41 eV). Although still somewhat underestimated we see that the exchange kernel is responsible for the improved accuracy. Self-consistency could further improve this result. That an xc kernel can improve the α-quartz-stishovite energy difference was also found in Ref. [12] using the renormalized PBE kernel (0.51-0.54

eV).
Our RPA result differs from the result presented in Ref. [11] by 20 meV. The origin of this difference, occurring already at the PBE level, can be traced to the use of PAW (projected augmented wave) potentials [71] in Ref. [11]. A larger difference can be found with respect to the RPA result presented in Ref. [12], a difference that can be traced mainly to the EXX part of the energy (see Appendix C).
From the technical point of view the situation is similar to the BN polymorphs. In the upper panel of Fig. 6 we show the convergence of the α-quartz-stishovite correlation energy difference (∆E q−s c ) with respect to N ν . Again, we find a slower convergence in RPA as compared to RPAx. The RPA and RPAx values were converged within 10 meV using 30 × N e and 10 × N e eigenvalues, respectively.

V. ICE POLYMORPHISM
The problem in predicting the relative energies between high-and low-pressure phases is relevant also for molecular solids such as the ice polymorphs. At ambient pressure the oxygen atoms in water form an hexagonal lattice held together by hydrogen bonds. By applying pressure the average O-O distance reduces such that the oxygen atoms are approached in non-hydrogen bonded configurations, which increases the coordination number from four up to eight in some dense high-pressure structures [4]. Calculating the energy difference between these phases with standard functionals in DFT, again produces large errors. One of the failures has been attributed to the lack of vdW forces, which are expected to play an important role in high-pressure phases as well as liquid water [75][76][77][78]. In addition, a fraction of HF exchange has shown to be important for the description of the hydrogen bond [77,78]. However, combining a hybrid functional with a vdW correction still does not give a complete description. The strongly constrained and appropriately normed (SCAN) meta-GGA functional [79] has shown that it is possible to obtain good energy differences between low-and high-density phases using a semi-local functional. On the other hand, lattice energies, i.e., the energy gain per monomer in the solid, remain overestimated with SCAN [80]. The more advanced RPA also provides improved energy differences, but in this case with underestimated lattice energies [81], due to the poor description of the hydrogen bond within RPA [13][14][15]. Calculations using second order Møller-Plesset perturbation theory have, indeed, shown that higher order exchange terms are important [82].
We will now investigate the performance of the RPAx, and start by exposing a fundamental problem of the RPA when calculating the polarizability of the water molecule. Since the static polarizability α(0) (Eq. (12)) is not a variational quantity like the total energy, it has a stronger dependency on the input orbitals (and eigenvalues) [83]. It can be argued that the correct orbitals should be those coming from a self-consistent RPA/RPAx calculation [28]. However, instead of these orbitals, we will use both PBE orbitals and orbitals from an optimal hybrid functional based on a local KS potential [84]. The optimization of the hybrid functional is carried out by minimizing the total RPA/RPAx energy with respect to the fraction of exact-exchange in hybrid functional used to generate the input orbitals (see further details in Appendix A and Ref. [85]). In order to study effects of anisotropy we calculate all diagonal elements of the polarizability. The results, obtained using the geometry of Refs. 73 and 74, are presented in Table I. We see that RPA clearly underestimates the polarizability and the magnitude of anisotropy when evaluated with optimal hybrid orbitals (20% of exchange in RPA). The isotropic RPA polarizability improves with PBE orbitals, but the anisotropy is then completely lost. The origin of any anisotropy within RPA is, thus, solely due to an effect of exchange in the description of the orbitals. In contrast, the RPAx exhibits anisotropy already with PBE orbitals stemming from the exchange kernel. This effect is enhanced using optimal hybrid orbitals (35% of exchange  in RPAx). The RPAx agrees very well with CCSD results [73], and are not so far from CCSD(T) results [74]. These results highlight the important role of the exchange kernel for describing not only the water molecule, but also the long-range interaction between water molecules.
In order to study the water intermolecular interaction in more detail we determine the binding energy of the ten Smith stationary points (SPs) [86]. In addition to the global minimum, the SPs correspond to transition states or higher order saddle points on the water dimer potential energy surface [72]. The geometries of the SPs, depicted Appendix C are relevant for liquid water [87,88], and some are found in high-pressure ice polymorphs. The lowest energy configuration (SP1) corresponds to the optimal hydrogen bond geometry, included in many molecular test-sets. Previous work has shown that RPAx and related methods capture the binding energy of this configuration rather well, while RPA fails by around 15% [13][14][15]. In Fig. 7 we have summarized the results for all ten SPs using structures optimized at the CCSD(T) level [72]. The HYB opt corresponds to the hybrid functional optimised with RPAx (35% of exchange). The HYB opt orbitals were used for evaluating the energy in both RPA and RPAx. For the total energy the choice of orbitals has a small effect. The largest difference is found for the SP1 configuration where the RPAx binding energy change from 212 meV to 216 meV when going from PBE to HYB opt orbitals. We have also included a vdW correction at the TS level in combination with the same hybrid (HYB opt +vdW). Errors are determined with respect CCSD(T) reference values [72]. It is clear that RPA fails not only for SP1 but for all configurations. The error is, however, rather constant suggesting that the energy differences are still quite well described. The hybrid functional performs very well for the first three hydrogen bonded configurations but fails, similarly to RPA, for the others. Energy differences will, in this case, not benefit from error cancellation. Adding the TS-vdW correction shifts all values with a similar constant (except for SP8). The vdW correction thus improves results for some configurations and worsens the results for others, but the quality of the energy difference will be similar to that of the hybrid functional without a vdW correction. These results show how complex and challenging it is to capture the full (H 2 O) 2 potential energy surface. RPAx is the only approximation that performs well for all SPs and produces results comparable to reference values at the CCSD(T) level.
In Fig. 8 we plot the dissociation energy curves of the SP1 and SP6 configurations. The dimers are stretched with fixed intramolecular geometries. Given the accuracy of the RPAx we can now use this approximation as a reference to compare the other methods too. It is interesting to note that the hybrid functional which gives a good value at the SP1 bond midpoint also captures the dissociation tail. RPAx and HYB opt perfectly coincides up to 5Å. Adding a vdW correction improves the binding energy of the SP6 configuration. However, from the full dissociation curve, we see that this is at the expense of producing the wrong tail. Even in this non-hydrogen bonded configuration, RPAx and HYB opt are closer at stretched geometries.
We are now ready to study solid forms of water. In Fig. 9 we illustrate four different phases of ice in the order of decreasing volume per molecular unit: XI, II, IX and VIII. All structures are relaxed with the PBE+TS functional. Due to the smaller unit cell, we have chosen the ferroelectrically proton-ordered ice XI (Cmc2 1 ) instead of the ordinary disordered I h [89,90]. The energy of XI is three meV lower than I h in RPA [81]. Both are purely hydrogen bonded solids with water molecules arranged in SP(1-3)-like configurations. In clusters and solids the strength of the hydrogen bond is enhanced as compared to the isolated dimer bond strength. This effect is captured by most functionals. In Table II we present lattice energies using HYB opt , HYB opt +vdW, RPA, RPAx and SOSEX. The results are compared to CCSD(T) [82] and recent QMC results [91]. We see that both HYB opt and II XI IX VIII RPAx produce a result for XI of similar quality as for the SP1 dimer. In the same way, RPA underestimates the lattice energy. We note that for the solids we have evaluated both RPA and RPAx with PBE orbitals in order to compare with previous RPA results [81]. The use of HYB opt orbitals (always via a local potential) increases the XI result by 10 meV (< 2%).
Ice IX and II have similar volumes, approximately 20% smaller than XI. The smaller volume makes nonhydrogen bonded configurations more relevant. In Fig. 9 we have identified pairs of H 2 O monomers for which the magnitude of the O-O distance is approaching that of hydrogen-bonded configurations. The ability of a functional to capture the energies of different dimer orientations thus increases in importance for the IX and II phases. In Table II we see that the good performance of HYB opt is now lost, producing underestimated lattice energies. Adding the vdW correction overcorrects the lattice energies but recovers the near degeneracies between XI, II and IX. This result was highlighted in Ref. [77] suggesting that the vdW forces play an important role in the II and IX phases. We note that the most recent QMC calculation predicts ice II and I h to be nearly degenerate at 613 and 615 meV, respectively (with a 5 meV error bar) [91], while an earlier QMC calculation from 2011 predicts ice II to be lower in energy than I h by 4 meV (609 and 605 meV respectively) [77]. The corresponding experimental results have been extracted in Ref. [92], yielding 609 and 610 meV respectively. Our RPA calculation gives an II-XI energy difference of 9 meV, in good agreement with previous calculations [81]. RPAx clearly reduces this energy difference making them degenerate. SOSEX lowers the energy of ice II with respect XI even further, reverting the energy order. Given that ice XI is expected to be a couple of meV lower in energy than I h , even by experiment [93], our results indicate that the II-I h energy difference should be ≤ 0 meV. We also note that HYB opt and HYB opt +vdW predict ice IX to be lower in energy than ice II, in contrast to the present (beyond-)RPA methods and the previously mentioned SCAN functional [80].
Ice VIII has a volume 40% smaller than ice IX and is stable above 2 GPa. It corresponds to a proton-ordered phase of VII [94,95]. In this structure the O atoms are eightfold coordinated and it is easy to identify both SP6-like and SP8-like dimers, in addition to the hydrogen bonded configurations. As pointed out in Ref. [4] the smallest O-O distance is actually between monomers in non-hydrogen bonded configurations. The experimental result for the lattice energy was initially determined to 577 meV [92]. Later, this has been revised to 595 meV [91]. Again HYB opt fails, while HYB opt +vdW only slightly overcorrects in this case. The energy difference with respect to ice XI remains, however, largely overestimated. The error of the RPA is slightly larger than for the other polymorphs. RPAx is again very accurate lying between CCSD(T) and QMC.
Overall the RPA results for the ice polymorphs are consistent with the analysis of the SPs above. Lattice energies are underestimated by approximately 13%, which is very similar to the error found for the SPs. Relative energies are much better described, benefiting from error cancellation. RPAx, which captures the essential correlation effects, gets both lattice energies and relative energies in good agreement with experiment and more sophisticated methods. SOSEX slightly overestimates the hydrogen bond and thereby overestimates the lattice energies for all polymorphs.

VI. CONCLUSIONS
The RPA is becoming an important tool in materials science for calculating ground-state properties of solids. There are, however, some inherent limitations of the RPA which can cause significant errors in certain systems or situations. In this work we have analyzed some of these cases and demonstrated how exact-exchange in the response function provides a theoretically well-defined, accurate and reliable improvement.
In general, the exchange term increases atomic polarizabilities leading to, e.g., improved vdW coefficients. Here Table II. Lattice energies (meV) of ice XI, II, IX and VIII with the optimized hybrid functional HYBopt (35% exchange), the same hybrid functional with a TS-vdW correction, RPA, RPAx, SOSEX, CCSD(T) [82], QMC [91] and experiment [92]. The relative energies with respect to XI are also presented. For CCSD(T) and QMC XI should be replaced by I h .

Ice
HYBopt HYBopt+vdW RPA RPAx SOSEX CCSD(T) QMC Exp. XI (Cmc21)  614  676  551  609  634  601  615  610  II  559  671  542  609  639  601  613  609  IX  570  673  538  598  625  --606  VIII  471  606  516  584  613  574  594  595  II-XI  55  5  9  0  -5  0  2  1  IX-XI  44  3  13  11  9  --4  VIII-XI  143  70  35  25  21 27 21 15 we have shown that for the cohesive energy of purely vdW bonded solids, such as Ar and Kr, RPAx produces very accurate results, improving the RPA by more than 20%. vdW forces also play an important role in the SiO 2 and BN polymorphs. We have shown that the interlayer binding energy of r-BN is enhanced with RPAx. However, with a smaller correction than for purely vdW bonded systems. The α-quartz-cristobalite energy difference is also enhanced with RPAx, and now agrees with reference QMC calculations. Earlier work have shown that the RPAx corrects the overestimated RPA correlation energy. As often pointed out, this error cancels to a large extent when looking at energy differences of similar systems. Which systems are similar enough is, however, not well-defined. In this work we have highlighted that when the coordination number changes, as between high-and low-pressure phases, errors do not cancel to a satisfactory degree. This has been demonstrated for the r-BN-c-BN and the α-quartzstishovite energy differences. In the first case, RPA and RPAx predict different energy ordering, and in the second case, RPAx strongly enhances the energy difference. In both cases RPAx produces results in line with highly accurate methods such as CC or QMC.
We have also studied the (H 2 O) 2 potential energy surface and ice polymorphism. The errors of RPA (and hybrid functionals) in describing the lattice energies of ice IX, II, XI and VIII were analysed in terms of the polarizability of the water molecule and the 10 SPs. The RPA clearly underestimates both magnitude and anisotropy of the polarizability, but captures rather well the energy differences between different H 2 O dimer configurations thanks to a good error cancellation. This implies under-estimated lattice energies of ice but good relative energies. The RPAx, which includes higher order exchange effects, captures the correct strength of the hydrogen bond and is able to describe the full anisotropic interaction between H 2 O monomers. Not only relative energies but also lattice energies are, thereby, in very good agreement with experiment and more sophisticated methods.
Regarding the computational cost and feasibility of the RPAx the main constraint is the N 3 k scaling, which can make k-point demanding systems several times more expensive than RPA. On the other hand, the convergence with respect to eigenvalues in the response function is in many cases faster with RPAx. Other parameters such as the total number of k-points or frequency points have a similar behaviour in RPA and RPAx. For the systems studied here, the cost of the RPAx calculation is 3-10 times that of the RPA calculation.
All calculations were done non-self-consistently, on top of a PBE ground-state. For the water dimers and ice we showed that the use of hybrid orbitals, which can be considered as a step towards self-consistency, increases the energy differences by a few meV only. Further developments in this direction are, however, interesting for eventually calculating forces and lattice vibrations within the same framework.
To conclude, we have shown that the RPAx opens up for highly accurate density functional calculations on solids that can address structural phase transitions involving distinct but nearly degenerate phases, and provide more reliable reference results than RPA when experimental data are missing or difficult to access.

ACKNOWLEDGMENTS
The authors would like to thank Lorenzo Paulatto, Guillaume Ferlat and Michele Casula for discussions.
The work was performed using HPC resources from GENCI-TGCC/CINES/IDRIS ( Figure 10. The α-parameter of the hybrid functional was chosen so as to minimize the total energy of the water molecule. The minimum occurs at 20% in RPA and at 35% in RPAx. (N k ). We used an energy cut-off of 80 Ry in all calculations except for the SiO 2 polymorphs where we converged the correlation energy with 60 Ry.
The convergence with respect to N ν for the α-quartzstishovite and c-BN-r-BN energy difference is presented in the main text. Here, in Fig. 11, we present the N νconvergence for the cohesive energy of Ar, the interlayer binding energy of r-BN, the α-quartz-cristobalite energy difference and the binding energy of the water dimer. For these systems a value N ν < 10 * N e is sufficient to converge within 2 meV for both RPA and RPAx.
In Fig. 12 we present the convergence of the correlation energy with respect k-points for Ar, c-BN, stishovite and Ice-XI. We used shifted grids for faster convergence. Ar and Kr were converged within 2 meV using a 5×5×5 grid. For c-BN the last point corresponds to a 6×6×6 grid. For stishovite the last point corresponds to a 4 × 4 × 4 grid. In order to save computing hours, the last RPAx point for stishovite is obtained by extrapolation, using the fact that the ratio between the RPA and RPAx correlation energy differences is essentially a constant. The last point for Ice-XI corresponds to a 3 × 3 × 2 grid.