Double-Hybrid DFT Functionals for the Condensed Phase: Gaussian and Plane Waves Implementation and Evaluation

Intermolecular interactions play an important role for the understanding of catalysis, biochemistry and pharmacy. Double-hybrid density functionals (DHDFs) combine the proper treatment of short-range interactions of common density functionals with the correct description of long-range interactions of wave-function correlation methods. Up to now, there are only a few benchmark studies available examining the performance of DHDFs in condensed phase. We studied the performance of a small but diverse selection of DHDFs implemented within Gaussian and plane waves formalism on cohesive energies of four representative dispersion interaction dominated crystal structures. We found that the PWRB95 and ωB97X-2 functionals provide an excellent description of long-ranged interactions in solids. In addition, we identified numerical issues due to the extreme grid dependence of the underlying density functional for PWRB95. The basis set superposition error (BSSE) and convergence with respect to the super cell size are discussed for two different large basis sets.


Introduction
Electronic structure calculations for realistic condensed-phase systems are generally more involved than those for molecules. The former include more atoms and are performed under periodic boundary conditions (PBC), implying interactions between periodic images. Therefore, condensed-phase electronic structure modelling often relies on simple approximations. Tight-binding approaches-semiempirical methods, density functional based tight-binding (DFTB)-used to be the work horse in the field. With increased computational power Kohn-Sham density functional theory (KS DFT) [1] became a standard approach. Recently, implementations of wave function theories (WFT) became available, although their application is far from routine.

Theoretical Background
In the following, a, b, . . . are virtual orbital indices, i, j, . . . occupied orbital indices, p, q, . . . general orbital indices, and P, Q, . . . auxiliary function indices. In DFT, the total energy is given as a functional of the total ground-state density n( r): where E DFT [n] is the total energy functional, T 0 [n] is the kinetic energy of a reference system of non-interacting electrons, E ne [n] is the nuclei-electron interaction energy, E H [n] is the Hartree energy describing the classical electron-electron interaction energy, and E XC [n] is the exchange-correlation energy describing the quantum mechanical contributions of the electron-electron interaction.
The ground-state density is expressed in terms of orbital functions ψ i ( r) n( r) = ∑ i |ψ i ( r)| 2 (2) where i runs over all occupied orbitals. The orbital functions fulfill the orthonormality constraint: with the Kronecker delta δ ij . The orbitals are solutions of the Kohn-Sham (KS) equation: with the potential arising from the nuclei v ne ( r), the Hartree potential v H [n]( r), the exchangecorrelation (XC) potential v XC [n]( r), and the orbital energy i of orbital i. In this article, we will consider Gaussian functions centered at the atoms only. Because the total energy functional is not known explicitly in terms of the ground-state density, we rely on approximations of the XC functional. These approximate energy functionals are given as integrals of a function explicitly depending on the ground-state density, its gradient and its Laplacian. For convenience, the XC functional is split into an exchange functional E X [n] and a correlation functional E C [n].
The more complex hybrid density functionals (HDFs) [16] include explicit information of the occupied orbitals. They modify the exchange functional by including a certain amount α X,HF of Hartree-Fock (HF) exchange E X,HF [n] providing the exchange functional We introduce the amount α X,DFT of DFT exchange E X,DFT [n] to reflect that the DFT exchange functional is an already known GGA or meta-GGA functional (compare [15,16]). Using the Mulliken notation (chemists' notation) for electron repulsion integrals the HF exchange energy can be written as Non-HDFs suffer from self-interaction errors [46]. These are reduced in HDFs but usually not fully cancelled since α X,HF = 1 in general case. This self-interaction error results in erroneous description of charge-separation processes and transition states. But even hybrid methods and HF lack a reasonable description of dispersion interactions decaying like R −6 with R being a measure of charge separation.
For increased flexibility, we can further split the exchange functional in a long-range and a short-range functional and describe both with a given mixture of HF theory and DFT resulting in range-separated HDFs [19].
The highest flexibility is achieved by including virtual orbitals ψ a ( r). Double-hybrid density functionals (DHDFs) are HDFs in which the correlation functional is composed of a mixture of a DFT correlation functional E C,DFT [n] with ratio α C,DFT and correlation energy E C,WFT [n] of a WFC method with ratio α C,WFT providing a functional Because WFC methods are computationally more demanding than HDFs or standard DFT functionals, most DHDFs exploit the MP2 theory, the SOS-MP2 theory, or the RPA method.
The correlation energy within the MP2 theory for closed-shell systems is The computationally most expensive step of the MP2 method is given by the transformation of the electron interaction integrals from atom orbital basis to molecular orbital basis leading to a O(N 5 ) scaling with N being a measure of system size. The prefactor can be reduced by the resolution-of-the-identity (RI) approach introducing an auxiliary basis in which densities are expanded giving the equation This method is called RI-MP2 [47,48]. A simplified version of the RI-MP2 method is the Scaled-Opposite-Spin(SOS)-MP2 method [49] given by with The integration is carried out numerically using a Minimax quadrature. The RI-SOS-MP2 method scales like O(N 4 ).
Another correlation method with increasing popularity is the Random Phase Approximation (RPA) method [50,51] within the RI approximation with RI-RPA scales like O(N 4 ). As with the RI-SOS-MP2 method, the integration is carried out numerically using a Clenshaw-Curtis grid [52] or a Minimax grid [53,54]. All WFC methods and all DHDFs correctly reproduce the R −6 energy behaviour of long-range interactions. Comparable to range-separated HDFs, there are DHDFs with range-separated exchange functionals like the ωB97X-2 functional [55]. Further, there are DHDFs with range-separated correlation functionals [56]. In this article, we will not focus on DHDFs with range-separated correlation functionals and refer to the literature [40,[57][58][59][60] for more details.

Gaussian and Plane Waves Method (GPW) and Integral Evaluation
The Gaussian and plane waves method (GPW) [61] allows for efficient periodic calculations with Gaussian basis sets using a dual representation of the electronic density and molecular orbitals. It assumes the use of a primary Gaussian basis for the expansion of matrix quantities (density matrix, KS matrix) and an auxiliary plane waves (PW) basis for the evaluation of the Hartree potential and the numerical integration of density functionals. To converge GPW calculations, one has to pay attention to both the size and quality of the Gaussian basis and the energy cutoff for the PWs. In the current implementation, GPW is used for the calculation of the Hartree potential, XC functionals, and two and three center integrals necessary for the RI-MP2 and RI-RPA methods. Exchange integrals are computed analytically using a truncated Coulomb potential [62].

Test Systems
Because we are interested in the description of intermolecular interactions, we are testing the functionals on molecular crystals (NH3, HCN) and rare-gase crystals (Ar, Ne) which have been studied by Sansone et al. [63]. Structural information of the unit cells are summarized in Table 1. In case of molecular solids, we were using structures reoptimized at the B3LYP-D* level [64].
Molecular crystals represent systems with a mixture of covalent bonding and dispersion interactions. NH 3 and HCN crystals additionally contain hydrogen bonds which are crucial for the discription of proteins. In contrast to that, there are only dispersion interactions within the rare-gas crystals. This results in low cohesive energies and the need for well-balanced functionals. Table 1. Structural information about the bulk structures used in this study. n f u is the number of formula units per unit cell. References for the geometrical information of the respective system are provided in the last column. Please note that there was a mistake in the cell parameters of CO 2 provided in reference [63].

Parameters of the Calculations
All calculations have been carried using a development version 8.0 of CP2K [42]. To ensure convergence with respect to the density cutoff, we were using high cutoffs of 1500 Ry for all RPA and MP2 calculations, 4000 Ry for remaining calculations of the molecular crystals NH 3 and HCN, and 10, 000 Ry for the rare-gas crystals Ar and Ne (see Section 4.1 for more details) and a relative cutoff of 50 Ry. For the rare-gas crystals, we set the parameters EPS_DEFAULT, EPS_PGF_ORB, EPS_SCF, and EPS_SCHWARZ in the HF section to 10 −30 , 10 −50 , 10 −5 , and 10 −10 , respectively, for the molecular crystals, we were using for the same parameters 10 −20 , 10 −40 , 10 −5 , and 10 −9 , respectively (see the CP2K manual for the meaning of these parameters). HF calculations for the bulk systems were using a truncated Coulomb potential with a cutoff radius of roughly half the super cell size. All densities have been smoothed using the NN10 method.
RI-MP2, RI-SOS-MP2 and RI-RPA calculations have been carried out using the GPW method to determine all integrals with a primary cutoff of 300 Ry and a relative cutoff of 50 Ry. We have exploited an 8-point minimax grid for all RI-RPA and RI-SOS-MP2 calculations.

Choice of Functionals and Implementation
We carried out calculations at the PBE [10], ωB97M-V [72], ωB97X-2 [55], PW6B95 [73], PWRB95 [74], SOS-PBE0-2 [75], RI-MP2 [35,47] and RI-RPA [51] levels of theory. PBE and RI-MP2 are used to compare differences between valence-only calculations of our valence-only calculations and the all-electron calculations of Sansone et al. [63]. PW6B95 is a meta-hybrid functional which performed best for weakly interacting systems with more pronounced dispersion interactions. PWRB95 is its RPA-based DHDF. ωB97M-V is a dispersion-corrected range-separated meta-hybrid functional. ωB97X-2 is an MP2-based DHDF with range-separated exchange functional. SOS-PBE0-2 is a RI-SOS-MP2-based DHDF. With this choice, we cover a large variety of different flavours of meta-hybrid and DHDF theories. Due to very high computational cost, we have restricted ourselves to this small, but representative set of functionals: one for each flavour of DHDFs and a corresponding HDF.
PBE calculations have been carried out using the CP2K implementation of PBE. For the ωB97M-V, and the PW6B95 functionals, we exploited the implementations of the LibXC library [76], version 4.3.4. Since the VV10 dispersion correction is not available in CP2K, we relied on the rVV10 correction and the parametrization suggested by Mardirossian et al. [77]. For the ωB97X-2 and PWRB95 functionals, we implemented the required parameter sets into the LibXC library.

Basis Sets and Pseudopotentials
The MP2 and RPA implementations within CP2K rely on a pseudopotential (PP) approach with Goedecker-Teter-Hutter PPs [78]. For the PBE functional, we used PPs optimized for PBE, for RPA and MP2 calculations, we were using PPs optimized for HF whereas for both HDFs and all DHDFs, we utilized PPs optimized for the PBE0 functional. All PPs have been taken from the Github repository of Jürg Hutter [79].
Correlation-consistent primary basis sets and suitable auxiliary basis sets of double zeta (DZ) and triple zeta (TZ) quality for the elements C, H, N and O have been taken from Del Ben et al. [80] We have optimized appropriate correlation-consistent primary and auxiliary basis sets of the same qualities for Ne and Ar using the polarization functions of the respective Dunning basis sets [48,81,82]. All PPs and primary and auxiliary basis sets are compiled in the Appendxes A-C.

Cohesive Energies and Basis Set Superposition Error
To determine total energies per formula unit, we carried out calculations of 2 × 2 × 2, 3 × 3 × 3 and 4 × 4 × 4 supercells of all given unit cells and used a linear fit of the total energy per formula unit against the inverse of the cell volume.
Because calculations of cohesive energies usually suffer from basis set superposition errors (BSSE), we perform a counterpoise correction [83]. The BSSE-free cohesive energies E coh are calculated according to with the total bulk energy per formula unit extrapolated to infinite cell volume E bulk , the energy of the molecule with ghost atoms E mol+ghost,bulk , the energy of the molecule using the bulk geometry E mol,bulk , and the total energy of the molecule using an optimized gas phase structure E mol,gas . For Ar and Ne, we trivially have E mol,bulk = E mol,gas . The corresponding BSSE is given by For the BSSE calculations, we took the crystalline structures, chose one molecule (or atom for Ar and Ne) surrounded by all ghost atoms within a 3 × 3 × 3 supercell.

General Remarks
We found the convergence of total energies of meta-HDFs PW6B95 and PWRB95 requires very tight energy cutoffs for the auxiliary PW basis of at least 4000 Ry. In contrast to that, calculations with the other meta-HDFs in our benchmark study, ωB97M-V, provided reasonable results with a cutoff of only 1200 Ry. Because the basis functions for the elements argon and neon are more localized than those for hydrogen, carbon and nitrogen, higher cutoffs for the noble gases were needed for an adequate representation of the basis functions of these elements on the grid.
It is well-known that GGA functionals and especially meta-GGA functionals require very tight integration grids for convergence and thus accurate results. Such cutoffs reflect numerical issues and the need for very fine integration grids when using the PW6B95 and PWRB95 functionals. Such grids are not necessary for the ωB97M-V functional which has been optimized with coarser integration grids in mind [72]. Thus, energy differences converged faster with ωB97M-V and PBE. Nevertheless, the total energies were not converged. To remove any possible problems due to incomplete convergence with respect to cutoffs, we utilized unusually high cutoffs for all density functionals.
Furthermore, we have found convergence problems with the PW6B95 and PWRB95 functionals, which can be resolved with density smoothing. Unfortunately, in some cases an increase of the energy cutoff for the density resulted in SCF convergence issues which could not be resolved with tighter filter thresholds. Nevertheless, we were able to achieve convergence by restarting the calculations with a higher cutoff starting from the converged SCF results with a lower cutoff. This was not possible for argon, where we exploited a cutoff of 4000 Ry for the PW6B95 and PWRB95 functional. Thus, some numbers for the PW6B95 and PWRB95 functionals are not fully converged with respect to the density cutoff.
Due to the higher computational costs, we have not carried out calculations of the 4 × 4 × 4 supercells on the TZ level.
All cohesive energies are compiled in Tables 2 and 3.  Table 3. Same as Table 2, but with basis sets of TZ quality. b exploits basis sets of augmented DZ quality.

Convergence with Respect to Super Cell Size
In Figure 1, we compiled the differences in total energies per formula unit relative to the extrapolated total energies. In general, we expect the total energies to decrease with increasing supercell size and the extrapolated value is a lower bound for the total energies of the super cells.
Our results show exactly this behaviour.  An important question is for which supercell size the error becomes negligible. A useful magnitude is given by the chemical accuracy of 4 kJ·mol −1 . For weakly-interacting systems such as rare-gas crystals with cohesive energy of less than chemical accuracy, the order of magnitude is set by the cohesive energy itself. As the error of a method should be not larger than chemical accuracy, the allowed error of the supercell method must be at least one order of magnitude smaller then the methodological error, i.e., not larger than 0.4 kJ·mol −1 . We find that a 3 × 3 × 3 super cell provides sufficient accuracy for all functionals and test systems. This behaviour is in agreement with the literature [80]. Sometimes, the total energy per formula unit of the 4 × 4 × 4 super cell has a higher magnitude than this of the 3 × 3 × 3 supercell, which may be due to numerical issues. For PBE, a cubic fit does not seem to be appropriate, and an exponential fit should be used instead.

Convergence of the BSSE
The BSSEs for the different test systems are compiled in Figure 2. First, we would like to point out that the BSSE is significantly larger for the molecular crystals than for the rare-gas crystals. This might be related to the larger number of atoms per molecule and to the spread of the basis functions. Since the effective core charge of rare-gas atoms is larger than for carbon or nitrogen, the basis functions are more localized which results in weaker overlap with neighbouring atoms. This is supported by the smaller reduction in BSSE for Ar and Ne when we exploit larger basis sets. Thus, augmentation of basis sets must significantly reduce BSSEs of Ar and Ne. Indeed, diffuse basis functions actually improve cohesive energies as shown by Sansone et al. [63]. Molecular crystals are thus more suitable objects to study BSSE than rare-gas crystals. For both molecular crystals in the test set, the non-DHDFs PBE, PW6B95 and ωB97M-V, provide the smallest BSSEs whereas the two WFC methods MP2 and RPA have the largest BSSEs, as expected. The DHDFs have a BSSE between both classes of methods because they employ a mixture of DFT and WFC methods.
Furthermore, we note that the WFC methods in CP2K are implemented within the RI approximation employing an auxiliary basis set. This leads to an additional source of BSSE for RPA, MP2 and all the DHDFs because the addition of the auxiliary functions of the ghost atoms increases the overall accuracy.

Convergence with Respect to Basis Set Size
In numerous studies, it was shown that total energies from DFT calculations converge exponentially with respect to basis set size. In contrast to that, total energies from WFC methods converge cubically with respect to basis set size when employing correlation-consistent basis sets. Thus, most DHDFs are expected to have a cubic convergence with respect to basis set size but with a smaller prefactor. DHDFs employing a long-ranged Coulomb operator only and describing short-ranged interactions with a density functional, converge exponentially [40]. This behaviour is confirmed with our data compiled in  Since larger basis sets systematically reduce total energies, cohesive energies increase. We observe this behaviour for the WFC methods and almost all DHDFs. The slight difference for PWRB95 in case of Ar may be due to not full convergence with respect to super cell size. For the other functionals-PBE, ωB97M-V and PW6B95-the cohesive energy from the TZ basis set is sometimes higher, i.e., the system is weaker bound. One problem might be that the 2 × 2 × 2 super cells are not yet fully converged or the extrapolation scheme using a linear fit of the total energies versus the inverse of the volume is not appropriate and an exponential fit might be more suitable.
Next, we would like to discuss the results obtained for the molecular crystals NH 3 and HCN. They are bound together by covalent bonds, dipole-dipole interactions, and dispersion interactions. For both systems, the results with the RPA and MP2 methods significantly improve the results over GGA DFT functionals, MP2 even achieving chemical accuracy. The ωB97M-V functional also provides very accurate numbers. The PW6B95 functional, as PBE, systematically underestimates the cohesive energies with errors compatible to PBE. The PWRB95 functional significantly improves upon the results of its relative PW6B95, bringing them within 1 kJ·mol −1 from the experiment. The same holds for the ωB97X-2 functional compared with the ωB97M-V functional, although the DHDF is based on the non-meta-GGA HDF ωB97X [86]. One of the worst performing functionals is SOS-PBE0-2.
For the rare-gas crystals, the picture is more complicated because the absolute values of the cohesive energies are of the order of the chemical accuracy. As pointed out by Sansone et al. [63], augmented basis sets are required for these systems. Our cohesive energies from MP2 with a TZ basis set are only slightly lower than their result with a DZ basis set but still much worse than those with an augmented basis set for both Ne and Ar. Consequently, our results for Ne do not allow for an evaluation of the performance of these functionals and further studies employing either quadruple or augmented basis sets (which are to be constructed) are needed. Nevertheless, our results for Ar show that the ωB97X-2 functional provides a good description. The same holds for MP2, PW6B95 and PWRB95, although one needs further investigations with augmented basis sets.
This issue does not apply to the molecular crystals. Indeed, our cohesive energies with a TZ basis set are even lower than those with an augmented basis set. Thus, the use of augmented basis sets is not necessary for the molecular crystals. This result is important for reducing computational costs of HF calculations and low-scaling WFC methods.

Discussion
Because DHDFs can be considered to be a mixture of DFT and WFC methods, the flexibility of DHDF parametrizations can yield approaches more accurate than the parent DFT and WFC functionals. At the same time, they inherit the shortcomings of both classes. Due to the dependence on the grid parameters, the functionals PW6B95 and PWRB95 are more difficult to use than others: care must be taken to check whether the results are converged with respect to the grid parameters, in CP2K, the density cutoff.
As expected, PBE can only provide the order of magnitude for weakly interacting systems, although it converges fast with respect to basis set size and has a low BSSE. MP2 and RPA are more sensitive to the basis set size and exhibit large BSSEs. These methods provide a moderate accuracy for different systems with small basis sets.
Non-DHDFs benefit from lower BSSEs. The PW6B95 functional has high demands on integration grids. Both considered functionals also provide a moderate accuracy and should be favourable over MP2 and RPA with their higher computational costs.
The double-hybrid functionals PWRB95 and ωB97X-2 show excellent performance with moderate BSSEs and lower basis set incompleteness errors. Both have computational costs compatible to full MP2 or RPA calculations and inherit the need of fine integration grids for accurate results, especially for PWRB95.
The non-empirical SOS-MP2 based DHDF, SOS-PBE0-2, does not provide any advantage as compared to the original methods. It was pointed out by different authors [87,88] that non-empirical DHDFs usually perform worse than empirical DHDFs.

Conclusions
In this study, we examined a selection of different HDFs and DHDFs by computing cohesive energies in four different crystal structures. Our results show that DHDFs inherit the shortcomings of the underlying DFT functional (integration grids) and the underlying WFC method (computational costs, BSSE, basis set dependence). We were able to show that the PWRB95 and the ωB97X-2 functionals provide excellent accuracy for molecular and rare-gas crystals. The HDFs ωB97M-V and PW6B95 also provide reasonable accuracy for these systems, whereas the SOS-PBE0-2 functional underperforms and can not be recommended.
The exploited basis sets allow a good description of molecular crystals. For the rare-gas crystals, we showed that non-augmented basis sets are not sufficient to achieve energy convergence with respect to the basis set size. Due to the high computational costs, we leave studies with augmented basis sets (and the construction of those basis sets) as well as benchmarking more range-separated DHDFs for future prospect. Acknowledgments: This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project IDs mr25 and uzh1.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: