Efficient and accurate description of adsorption in zeolites

TL;DR

This paper compares RPA with singles corrections and MP2 for calculating adsorption energies in zeolites, finding RPA superior in accuracy and computational efficiency for porous materials.

Contribution

It demonstrates that RPA with singles corrections is more accurate and faster than MP2 for adsorption energy calculations in zeolites.

Findings

RPA with singles corrections outperforms MP2 in accuracy.

RPA with singles corrections is computationally more efficient.

RPA is suitable for accurate adsorption energy calculations in porous materials.

Abstract

Accurate theoretical methods are needed to correctly describe adsorption on solid surfaces or in porous materials. The random phase approximation with singles corrections scheme (RPA) and the second order M{\o}ller-Plesset perturbation theory (MP2) are two schemes which offer high accuracy at affordable computational cost. However, there is little knowledge about their applicability and reliability for different adsorbates and surfaces. Here we calculate adsorption energies of seven different molecules in zeolite chabazite to show that RPA with singles corrections is superior to MP2 not only in terms of accuracy but also in terms of computer time. Therefore RPA with singles is suitable scheme for obtaining highly accurate adsorption energies in porous materials and similar systems.

Tables

Table 1. Table 1: Adsorption energies of different molecules with bulk chabazite obtained for different methods. Data are in kJ/mol

System	HF	xMP2	MP2	EXX	RPA	+RSE	+GWSE
Methane	$4.6$	$0.1$	$-25.6$	10.1	$-23.0$	$-26.5$	$-26.4$
Ethane	$10.2$	$0.3$	$-37.0$	18.5	$-33.1$	$-38.5$	$-38.2$
Ethylene	$-3.2$	$-0.9$	$-54.2$	3.8	$-46.9$	$-51.6$	$-52.1$
Acetylene	$-6.4$	$-0.3$	$-49.1$	$-$0.2	$-41.7$	$-45.6$	$-46.3$
Propane	$10.4$	$0.4$	$-47.6$	19.5	$-43.2$	$-49.1$	$-48.8$
CO2	$-10.1$	$-2.9$	$-46.3$	$-$2.6	$-41.0$	$-45.9$	$-45.6$
H2O	$-54.4$	$-3.8$	$-82.7$	$-$45.0	$-76.4$	$-82.8$	$-82.1$

Table 2. Table 2: Adsorption energies of different molecules on the 2T cluster as obtained using RPA, RPA+RSE, MP2, and CCSD(T) which serves as a reference. All the data are in kJ/mol.

System	RPA	+RSE	MP2	CCSD(T)
Methane	$-$8.7	$-$11.1	$-$10.6	$-$11.0
Ethane	$-$9.4	$-$11.8	$-$10.9	$-$11.6
Ethylene	$-$18.7	$-$20.9	$-$23.0	$-$21.3
Acetylene	$-$14.0	$-$16.3	$-$18.4	$-$16.7
Propane	$-$10.0	$-$12.1	$-$11.1	$-$11.8
CO2	$-$12.0	$-$14.2	$-$14.1	$-$14.6
H2O	$-$39.8	$-$43.7	$-$44.7	$-$44.9

Table 3. Table 3: Adsorption energies of different molecules on the 4T cluster as calculated with RPA, RPA+RSE, MP2, and the reference PNO-CCSD(T)(F12*). All the data are in kJ/mol.

System	RPA	+RSE	MP2	PNO-CCSD(T)
Methane	$-12.8$	$-16.2$	$-15.3$	$-14.7$
Ethane	$-15.1$	$-18.8$	$-17.9$	$-17.4$
Ethylene	$-30.7$	$-35.3$	$-38.1$	$-32.7$
Acetylene	$-19.1$	$-22.7$	$-25.5$	$-21.8$
Propane	$-16.4$	$-19.9$	$-19.0$	$-18.2$
CO2	$-23.9$	$-27.7$	$-27.7$	$-28.9$
H2O	$-64.7$	$-71.0$	$-72.0$	$-72.0$

Table 4. Table 4: Adsorption energies of different molecules in bulk chabazite obtained for MP2, RPA+GWSE and two dispersion corrected DFT functionals. Data are in kJ/mol

System	MP2	RPA+GWSE	PBE-TS/IH	PBE-D3
Methane	$-25.6$	$-26.4$	$-$31.4	$-$34.9
Ethane	$-37.0$	$-38.2$	$-$47.1	$-$50.1
Ethylene	$-54.2$	$-52.0$	$-$56.1	$-$64.7
Acetylene	$-49.1$	$-46.3$	$-$50.1	$-$57.8
Propane	$-47.6$	$-48.8$	$-$65.0	$-$65.4
CO2	$-46.3$	$-45.6$	$-$39.8	$-$48.9
H2O	$-82.7$	$-82.1$	$-$83.7	$-$88.7

Table 5. Table 5: Adsorption energies of the tested molecules on the 2T model using the reference CCSD(T) scheme, RPA+RSE and dispersion corrected DFT. Data are in kJ/mol

System	CCSD(T)	RPA+RSE	PBE-TS/IH	PBE-D3	PBE-D4	PBE0-D4
Methane	$-11.0$	$-$11.1	$-$10.9	$-$13.0	$-$13.1	$-$12.6
Ethane	$-11.6$	$-$11.8	$-$11.3	$-$13.2	$-$13.2	$-$12.7
Ethylene	$-21.3$	$-$20.9	$-$22.6	$-$25.3	$-$26.0	$-$25.8
Acetylene	$-16.7$	$-$16.3	$-$17.0	$-$20.1	$-$20.9	$-$20.6
Propane	$-11.8$	$-$12.1	$-$12.2	$-$13.5	$-$13.3	$-$12.8
CO2	$-14.6$	$-$14.2	$-$12.5	$-$14.3	$-$14.6	$-$14.7
H2O	$-44.9$	$-$43.7	$-$43.7	$-$46.0	$-$46.5	$-$47.2

Full text

I Introduction

Adsorption of molecules on solid surfaces or in porous materials is a phenomenon important for heterogeneous catalysis in industry, for (de)mineralization in nature, and for many other processes. Computational modelling at atomic and molecular level greatly helps to identify structures of adsorbates or to elucidate reactions catalyzed by solid surfaces.Carrasco et al. (2009) In particular, Kohn-Sham density functional theory (DFT)Hohenberg and Kohn (1964); Kohn and Sham (1965) has been an extremely useful tool for this task. However, in many systems, substantial contributions to adsorption energies come from long-range electron correlations (dispersion forces).Tuma and Sauer (2006); Hammerschmidt, Müler, and Paulus (2012) These are not accounted for by standard semi-local or hybrid DFT approximations and they can thus lead to large errors for adsorption energies.Chakarova-Käck et al. (2006); Sony et al. (2007); Berland, Einstein, and Hylgaard (2009) Adding dispersion corrections usually reduces the errors but they can remain largeGautier et al. (2015) even for seemingly simple systems, such as semiconductors or oxides.Al-Hamdani et al. (2015); Zen et al. (2016a) This means that the DFT based approaches are still far from being reliable enough to be used to obtain reference quality adsorption energies.

To obtain highly accurate adsorption energies one needs a scheme that includes a high-level description of electron correlation effects. Methods that are capable of this task are quantum Monte Carlo (QMC) techniques, such as diffusion Monte Carlo or the coupled cluster scheme, using at least singles, doubles, and perturbative triples excitations, CCSD(T). For finite cluster calculations these methods have been shown to agree to within few per cent, see, e.g., the recent work in Refs. Dubecký et al., 2014; Al-Hamdani et al., 2017; Brandenburg et al., 2019; Dubecký et al., 2019. However, so far there has been less agreement between adsorption energies calculated within periodic boundary conditions,Brandenburg et al. (2019) even though algorithmic improvements are likely to reduce this issue in the future as well.Zen et al. (2016b); Hummel, Tsatsoulis, and Grüneis (2017); Zhang and Grüneis (2019) As an alternative, embedded cluster techniques can be used to calculate adsorption energies,Lopez et al. (1999); Di Valentin et al. (2002) and the recent combination with low scaling coupled cluster implementations seems particularly promising.Kubas et al. (2016) Another alternative is to use less elaborate methods that could offer high accuracy nevertheless. One of them is the second-order Møller-Plesset perturbation theory approximation (MP2) which has been widely used in quantum chemistry. In fact, it has been also used to study adsorption using embedded finite cluster approach, see, e.g. Ref. Boese and Sauer, 2016 for a recent work. Moreover, periodic implementation of MP2 was used to study adsorption in zeolites.Göltl et al. (2012) Another promising method available within periodic boundary conditions is the random phase approximation (RPA).Harl and Kresse (2008); Lu et al. (2009); Harl and Kresse (2010); Nguyen and de Gironcoli (2009); Olsen et al. (2011) RPA surpasses MP2 in description of electron screening but it lacks second and higher order exchange effects,Grüneis et al. (2009) the relative importance of either of these is still rather unclear in general.Ren et al. (2013) The accuracy of RPA, even though it was shown to be consistent for several systems,Harl, Schimka, and Kresse (2010); Bleiziffer, Heßelmann, and Görling (2012); Garrido Torres et al. (2017) was found to be too low to make it the method of choice. Several ways how to improve the accuracy of RPA have been explored, for example obtaining self-consistent RPA energiesBleiziffer, Heßelmann, and Görling (2012, 2013) or using exchange-correlation kernels.Olsen and Thygesen (2013a, 2014); Bleiziffer, Krug, and Görling (2015) Another approach that has lead to improved accuracy of RPA was to add the so-called singles corrections, either the renormalized singles excitations (RSE)Ren et al. (2011) or the $GW$ singles excitations (GWSE),Klimeš et al. (2015) that do not or do account for electron screening, respectively. The high accuracy of RPA with singles has been demonstrated for systems such as molecular solidsKlimeš (2016); Zen et al. (2018) and water adsorption on sodium chloride, hexagonal boron nitride, or graphene.Klimeš et al. (2015); Al-Hamdani et al. (2017); Brandenburg et al. (2019) However, the studies on adsorption considered only a single adsorbate and substrate type at a time and more general understanding of the accuracy and limitations of RPA with singles is missing.

In this work we want to understand the predictive power of MP2 and RPA with singles for the calculation of adsorption energies in simple systems. To achieve this goal we focus on adsorption of small molecules in zeolites, which are porous aluminosilicates. Previous studies have shown that, rather surprisingly, many DFT based schemes give adsorption energies which are too large.Göltl and Hafner (2011); Bučko et al. (2013); Shang et al. (2014); Bučko et al. (2014) This has even partly stimulated a development of improved dispersion correction methods.Bučko et al. (2014) In contrast, MP2 has been often shown to be close to CCSD(T) for finite cluster calculations of zeolitesTuma and Sauer (2015) and RPA with singles was shown to improve the description of CO adsorption on models of zeolites.Rubeš et al. (2018) Due to the problems in describing adsorption in zeolites and the current availability of both RPA with singles and MP2 it is interesting to ask: Which one of these is more reliable for the description of adsorption? Is it the widely used MP2 or RPA which, until 2001,Furche (2001) has received almost no attention in quantum chemistry? To answer these questions we calculated adsorption energies of seven molecules using MP2 and RPA with singles. While they give results which are surprisingly close to each other for adsorption in bulk, RPA with singles outperforms MP2 for adsorption on finite clusters. Moreover, its computational cost is one order of magnitude smaller than that of MP2. Our data, together with the previous resultsKlimeš et al. (2015); Al-Hamdani et al. (2017); Rubeš et al. (2018); Brandenburg et al. (2019) suggest that RPA with singles is the current method of choice for obtaining nearly reference quality adsorption energies in zeolites as well as other systems at affordable computational cost.

II Systems

The basic tetrahedral binding motif of silica groups gives rise to a large number of different zeolite structures which can have unit cell volumes of several thousands of cubic Ångströms. MP2 and RPA are methods based on perturbation theory and, compared to simpler DFT approaches, have larger memory and computational time requirements. Therefore, to obtain precisely converged adsorption energies, one needs a system with relatively small simulation cell, below approximately 2000 Å3. As was done in previous studies,Göltl and Hafner (2011); Göltl et al. (2012) we picked the chabazite structure with unit cell composition Si11AlO24H.

We calculated adsorption energies of seven molecules: methane, ethane, ethylene, acetylene, propane, CO2, and H2O. To obtain the required geometries we first optimized the chabazite with adsorbed methane using the optB88-vdW functional.Dion et al. (2004); Román-Pérez and Soler (2009); Klimeš, Bowler, and Michaelides (2010) The zeolite framework was held fixed in all the subsequent calculations. This was simply to reduce the number of calculations needed to obtain the results. Moreover, this is not an issue as our primary interest here is to understand the differences between MP2 and RPA adsorption energies. Structural optimization, in some cases combined with molecular dynamics, was used to obtain the adsorption structures of the other molecules. The unit cell of chabazite with adsorbed ethylene is shown in Fig. 1(a).

Furthermore, we have created two finite clusters, one with 2 tetrahedral sites (2T, AlSiO7H7) and one with four tetrahedral sites (4T, AlSi3O13H11) to assess the quality of the adsorption energies using reference quality method. For the clusters, the broken bonds were capped with hydrogens which were subsequently relaxed keeping the positions of all the other atoms fixed. The clusters that we used together with adsorbed ethylene are shown in Fig. 1(b) for the 2T cluster and in panel (c) of the figure for the 4T clusters. All the structures are given in the SI.sup

III Computational setup

The adsorption energies $E^{\rm ads}$ were obtained as $E^{\rm ads}=E_{\rm zeo+mol}-E_{\rm zeo}-E_{\rm mol}$, where $E_{\rm zeo+mol}$, $E_{\rm zeo}$, and $E_{\rm mol}$ are the energies of zeolite (either bulk or finite cluster) with adsorbed molecule, the bare zeolite, and the isolated molecule, respectively. In all three cases, the geometries of the molecule or zeolite are identical, that is we are calculating interaction energy. Moreover, all three calculations used the same simulation cell. This reduces errors due to k-point sampling and basis-set incompleteness. For the bulk calculations, there is a small residual error due to interaction with periodic images for the molecule, this amounts to 0.5 kJ/mol for water and is negligible in other cases. In any case this does not affect our comparison as the error is similar for different methods.

The calculations employing periodic boundary conditions were performed using the VASP program.Kresse and Hafner (1993); Kresse and Furthmüller (1996) VASP uses a plane-wave basis set and the projector-augmented wave (PAW) ansatz for the wavefunction.Blöchl (1994); Kresse and Joubert (1999) The MP2 calculations were performed using the standard implementation in VASP.Marsman et al. (2009); Grüneis, Marsman, and Kresse (2010) While the recent $O(N^{4})$ implementation of MP2 was shown to be more efficient for the systems considered here,Schäfer, Ramberger, and Kresse (2017) we did not use it as the corrections for the Gamma-point singularity of the Coulomb interaction were not yet available.Gajdoš et al. (2006); Liao and Grüneis (2016) We used the cubic scaling code for the RPA and singles calculations.Kaltak, Klimeš, and Kresse (2014a, b); Klimeš et al. (2015) Two sets of PAW potentials were used, first “standard” PAWs for testing and volume convergence for the finite clusters, and, second, “hard” PAWs which treat also semi-core electrons as valence and are almost norm-conserving.Klimeš, Kaltak, and Kresse (2014) The specific PAW potentials that we used are listed in the SI.sup

The total RPA energy is composed of the exact exchange (EXX) and random phase correlation (RPAc) energy components. The singles corrections can be added either in the RSE, or in the GWSE flavor. The Perdew-Burke-Ernzerhof (PBE)Perdew, Burke, and Ernzerhof (1996) functional was used as a starting point for the RPA calculations. The EXX and RSE adsorption energies converge quickly with the basis set size and we have used a 1050 eV basis-set cut-off to obtain their values. This cut-off guarantees a convergence of the adsorption energies to within 0.05 kJ/mol, the individual energies are provided in the SI.sup The RPAc and GWSE energies have a stronger dependence on the basis-set cut-off due to the two-electron cusp. We obtained the energies for several values of cut-offs and used data for cut-offs 850, 950, and 1050 eV to extrapolate to infinite basis-set size limit. In the case of GWSE, we extrapolated the difference between GWSE and RSE corrections. In the extrapolation, we assumed that the basis-set incompleteness error is proportional to $E_{\rm cut}^{-3/2}$. The auxiliary basis-set energy cut-off for response related properties (ENCUTGW tag for VASP) was set to one half of the basis-set cut-off (ENCUT).Harl and Kresse (2008) We note that the adsorption energy converges rather quickly with the basis-set size. Interestingly, this holds both for the adsorption energies obtained using a given combination of ENCUT and ENCUTGW and adsorption energies obtained for a given value of ENCUT and extrapolated with the auxiliary basis set ENCUTGW in a single calculation.Harl and Kresse (2008) This is likely because the change of density upon adsorption is small and the leading contributions to basis-set incompleteness errors for the adsorbed system and the individual parts cancel each other.Gulans (2014); Klimeš, Kaltak, and Kresse (2014) Using a basis-set cut-off of 750 eV and the ENCUTGW extrapolation would lead to errors of less than 1.3% or 0.4 kJ/mol in the adsorption energies, which is acceptable in most of applications.

The total MP2 energy is composed of the Hartree-Fock energy (HF) and the MP2 correlation energy, which itself can be divided into direct MP2 (dMP2) and exchange MP2 (xMP2) parts. The HF part was obtained using a basis-set cut-off of 1050 eV. The MP2 implementation has higher computational requirements compared to the RPA implementation and we were able to obtain data only for basis-set cut-offs of 650 eV and, for bulk only, 750 eV. The infinite basis-set limit of the xMP2 adsorption energies was obtained by extrapolating data calculated with 650 eV basis-set cut-off and cut-offs of 325 and 425 eV for the auxiliary basis set (ENCUTGW). For the bulk, the results are within 0.05 kJ/mol of data for which the second point for extrapolation used 750 eV cut-off for the basis-set and 375 eV cut-off for the auxiliary basis. One exception is CO2 where the difference is 0.14 kJ/mol, which is still small not to affect our conclusions. We note that the xMP2 contribution to energy is below 1 kJ/mol for the hydrocarbons while for CO2 and H2O the contribution is larger, being approximately $-3$ and $-4$ kJ/mol, respectively. For the dMP2 part, we obtained additional data using the RPA code and extrapolated to the complete basis set limit using cut-off values of 850, 950, and 1050 eV. As for the RPA calculations, the auxiliary basis-set cut-off energy was set to one half of the cut-off for the orbitals. In the bulk calculations of EXX and HF energies, a 2$\times$2$\times$2 k-point set was used, for all the other calculations (singles corrections and RPA and MP2 correlation energies) only the gamma point was used for k-point sampling.

Interaction energies on a finite cluster were obtained using VASP for RPA, RSE, and MP2 and using Turbomole, which provided basis set limit CCSD(T) reference and MP2 energies for comparison. In the VASP calculations, the HF and EXX energies were obtained using a 1050 eV basis-set cut-off and a cell with at least 18 Å side. The RPA and MP2 binding energies were obtained by a composite procedure. First, a basis-set converged data were obtained in a cubic cell with 9 Å side, using the same settings for basis-set cut-offs as in bulk. Second, correction to an infinite cell volume was added to obtain data converged with both the basis and cell volume. The finite size correction was calculated as the difference between the binding energy in a 9 Å cell and binding energy extrapolated to infinite volume using cells with up to 12 Å side. Similar correction was used to obtain the infinite cell limit of the RSE correction.

The Turbomole calculations used the explicitly correlated coupled cluster approach CCSD(T)(F12*)Hättig, Tew, and Köhn (2010) and the cc-pVTZ-F12 basis set.Peterson, Adler, and Werner (2008) The Slater-type correlation factor was used with an exponent of 1.0 a${}_{0}^{-1}$, together with the specially optimised RI basis sets.Yousaf and Peterson (2008) We also obtained MP2 energies using the MP2-F12 approachBachorz et al. (2011) and cc-pVTZ-F12 basis set to provide reference data to compare to VASP. There is a satisfactory agreement between the MP2 adsorption energies obtained with VASP and Turbomole, the average absolute difference is only 0.1 kJ/mol. Canonical coupled cluster calculations were possible for all cases except for the 4T clusters, which were performed using the PNO-CCSD(T0)(F12*) approach.Schmitz, Hättig, and Tew (2014); Schmitz and Hättig (2016) We found that a very tight value of $T_{\text{pno}}$ parameter is needed to obtain converged adsorption energies and we used $T_{\text{pno}}=10^{-8}$. To obtain the full contribution of triples, we scaled the T0 value by the T/T0 ratio obtained using the cc-VDZ-F12 basis set, which was possible to compute for all 4T clusters.

IV Results

IV.1 Adsorption in bulk material

The adsorption energies of the different molecules in bulk chabazite are summarized in Table 1 and plotted in Fig. 2. Looking at the graph, the first striking thing is the close agreement between MP2 and RPA with singles corrections, either RSE or GWSE. As the table shows, MP2 and RPA+GWSE adsorption energies differ only by 0.6 kJ/mol and 0.7 kJ/mol for water and carbon dioxide, respectively. The differences are larger for other systems, but they are still close to one or two kJ/mol or few per cent. Given the differences between RPA with singles and MP2, the close agreement points to a fact that they provide highly accurate adsorption energies. Moreover, the agreement between RPA with singles and MP2 is unexpected as the results previously published by Göltl et al. imply that there should be a much larger difference, e.g. at least 8 kJ/mol for propane.Göltl et al. (2012) In passing, we note that for the systems considered here the computational cost of RPA with GWSE corrections is about an order of magnitude smaller than the cost of MP2 and the memory requirements are few times smaller as well.

The differences between the RSE and GWSE corrections are small, typically below one per cent so that the data almost overlap in Fig. 2. Looking at the results in a more detail, two cases emerge, for the first, adding GWSE over RSE leads to a small reduction of binding. This is the case of ethane, propane, CO2, and water, where the binding is reduced by 0.3 kJ/mol for the first three and by 0.8 kJ/mol for water. For ethylene and acetylene, using GWSE increases the RPA+RSE binding by 0.5 kJ/mol and by 0.7 kJ/mol, respectively. We note that increase of binding upon adding the GWSE corrections over RSE was also observed for molecular solids where the constituent molecules contained delocalized $\pi$ electrons.Klimeš (2016) For such systems screening is important and the GWSE corrections are expected to be sizeable.Klimeš et al. (2015)

The physical reason why the singles corrections increase the binding is that they change the density matrix towards the HF state,Klimeš et al. (2015) which is known to bind stronger compared to EXX based on PBE orbitals.Ren et al. (2011) Other properties, such as the electron density, are affected as well upon adding the singles. For example, we have recently shown that the charge density difference ($\Delta\rho$) for noble gas dimers obtained with singles corrections is close to HF and coupled cluster density differences in the bonding region.Modrzejewski, Yourdkhani, and Klimes (2019) We observe similar changes in the density differences for the systems studied here. Specifically, in Fig. 3 we show the charge density difference for adsorption of methane, the line profile runs along the direction of the O-H bond. Clearly, the charge accumulation (positive values) and charge depletion regions as given by PBE are shifted in the direction of the acidic hydrogen (small values of $x$) compared to HF. When the singles corrections are added to PBE, $\Delta\rho$ becomes closer to the HF result. This is consistent with the results observed for the neon dimer.Modrzejewski, Yourdkhani, and Klimes (2019)

The RPA is known to give adsorption or binding energies which are too weak.Ren et al. (2011); Eshuis and Furche (2012); Bleiziffer, Heßelmann, and Görling (2012); Olsen and Thygesen (2013b, 2014); Klimeš et al. (2015); Al-Hamdani et al. (2017) This holds also here if we consider the MP2 or RPA+GWSE values as the reference. The specific values of the GWSE corrections, that can be taken as a rough estimate of the RPA error, are then between 3.5 and 5.8 kJ/mol. These values amount to around 10% of the adsorption energy, which is consistent with typical errors of RPA binding energies.

Let us now compare our results to those presented previously in the literature. Göltl and coworkersGöltl et al. (2012) have obtained RPA and MP2 adsorption energies for methane, ethane, and propane in chabazite for four different positions of the acidic hydrogen. Moreover, Piccini et al. have obtained adsorption energies for the same molecules using the MP2:PBE-D approach.Tuma and Sauer (2006); Piccini et al. (2015) The RPA adsorption energies of methane, ethane, and propane obtained by Göltl and coworkers ($-$19.69, $-$27.65, and $-$34.96 kJ/mol, respectively) are several kJ/mol smaller in magnitude compared to the values obtained here. Note that our structures have the acidic hydrogen bound to the O(1) oxygen of chabazite and we compare to the data of Göltl et al. for the same position. We have found that part of the difference could be attributed to the PAW potentials used in Ref. Göltl et al., 2012. Göltl and coworkers used “standard” PAW potentials while we have used “hard” PAW potentials with small core radii.Klimeš, Kaltak, and Kresse (2014) We find that the “hard” potentials give binding which is 0.5 to 1 kJ/mol stronger compared to the “standard” ones for all the systems considered here. This difference is significant when comparing the quality of different methods as well as when comparison to experimental data is made. The use of hard PAWs or pseudopotentials could be especially important when trying to understand difficult systems, such as the structure of adsorbates in zeolites under high pressure conditions.Arletti et al. (2017); Fischer (2019) Turning to MP2, our adsorption energy for methane agrees well with the value of $-25.13$ kJ/mol reported by Göltl et al. for the O(1) structure. Piccini et al. gives a similar value ($-$25.34 kJ/mol), however, in their structure hydrogen is bound to the O(4) oxygen. (We note that in Ref. Piccini et al., 2015 a comparison is made to Göltl’s data for the O(3) structure.) Göltl and coworkers found that the structures where the acidic hydrogen is bound to oxygen O(1), as used here, and to oxygen O(4), as used by Piccini et al., give adsorption energies which are larger by 0.25 and 1.32 kJ/mol for ethane and propane, respectively. Our MP2 adsorption energies are larger by 0.8 kJ/mol for ethane and by 3.5 kJ/mol for propane compared to those of Piccini et al. and part of this difference could be attributed to the different adsorption site. Given the differences in computational set-up and adsorption structures, the remaining difference is reasonably small and we consider our data and the data of Piccini and co-workers to be in a good agreement. The MP2 adsorption energies of Göltl and coworkers for ethane and propane, $-40.68$ and $-55.49$ kJ/mol are then too large compared to our or Piccini’s data.

The close agreement between RPA with singles and MP2 calls for a more detailed comparison of the results which reveals clear differences, visible already in Fig 2. Specifically, alkanes (methane, ethane, propane) are bound less by MP2 than by RPA with singles (either RSE or GWSE). In contrast, the MP2 binding is stronger for ethylene and acetylene, by 2.1 and 2.8 kJ/mol, respectively, when comparing to RPA+GWSE. Comparing the different schemes, MP2 includes the second-order exchange contribution to the correlation energy, the xMP2 term, which is not accounted for in RPA. In contrast, RPA sums a class of terms up to infinite order and accounts for screening. The lack of screening leads to too large interaction energies for systems with delocalized electrons, and this would help to explain the larger MP2 adsorption energies for ethylene and acetylene. The xMP2 term is expected to contribute the most to the adsorption energy when there is a significant overlap between the fragments. Consequently, the xMP2 contributions to the adsorption energies for the hydrocarbons are small, less than 1 kJ/mol, see Table 1. Only for CO2 and water the xMP2 contributions are larger, being $-$2.9 and $-$3.8 kJ/mol, respectively. For all the systems, the difference between MP2 and RPA+GWSE adsorption energies has the same sign as the xMP2 contribution and adding terms such as the second-order screened exchangeGrüneis et al. (2009); Hummel et al. (2019) could further improve the accuracy of RPA. We reiterate that the original goal of this work was to understand how well RPA with singles performs compared to MP2 which we expected to be accurate enough to provide reference data for adsorption energies in zeolites. This is based on the small corrections between MP2 and CCSD(T) for finite clustersPiccini et al. (2015); Tuma and Sauer (2015); Rubeš et al. (2017) and also on previous results of Göltl and coworkers.Göltl et al. (2012) With the current results at hand, the question then arises, if the RPA with singles method is actually more accurate than the (more computationally demanding) MP2 scheme. Given the small energy differences involved, very precise reference data are needed and such can’t be currently obtained for the periodic material. Therefore, we used finite cluster models of the adsorption site and we discuss the results now.

IV.2 Adsorption on clusters

We now discuss the results obtained for the finite clusters, starting with the smaller 2T one. Table 2 summarizes the adsorption energies, relative differences from the CCSD(T) reference data are presented in Fig. 4. One can see that in all but one case the RPA+RSE adsorption energies are closer to the CCSD(T) reference than are the MP2 adsorption energies. Not only that, the results of RPA+RSE are rather consistent, the relative errors are within 3% for all the systems. In contrast, MP2 gives almost no error for water while the error reaches 10% for ethylene and acetylene. In agreement with previous results, the adsorption energies predicted by RPA are underestimated by 10 to 20%.

Finally, we consider adsorption on the large 4T cluster for which the results are summarized in Table 3. Comparing first MP2 with RPA+RSE, one can make similar observations as from the data for bulk and the 2T cluster. Specifically, alkanes are bound by about 1 kJ/mol less by MP2, ethylene and acetylene are bound by about 3 kJ/mol or 10% more by MP2, the binding is within 0.1 kJ/mol for CO2. MP2 predicts 1 kJ/mol stronger binding for water compared to the RPA+RSE, this is similar to the 2T cluster, but opposite to the binding in bulk. Table 3 also lists the PNO-CCSD(T)(F12*) data for reference. The data show overbinding of ethylene and acetylene by MP2. The alkanes are bound too strongly by RPA+RSE, by about 1.5 kJ/mol which is more than for the smaller 2T cluster where the difference was at most 0.3 kJ/mol for propane. A similar increase in the binding strength for the alkanes can be observed for MP2. MP2 gave about 0.5 kJ/mol too weak adsorption energies on the 2T cluster while for the 4T cluster the binding is about 0.5 kJ/mol too strong.

IV.3 Comparison to density functional theory methods

We now use the adsorption energies calculated with RPA with singles scheme to assess the predictions made by dispersion corrected DFT functionals. It is well known that adsorption energies in zeolites predicted by current DFT approximations are not very accurate and are typically substantially overestimated.Göltl et al. (2012); Thang et al. (2014); Shang et al. (2014); Piccini et al. (2015); Rubeš et al. (2018); Rocca et al. (2019); Piccini, Alessio, and Sauer (2018) One exception is the DFT/CC scheme proposed by Bludský and coworkers.Bludský et al. (2008) Here, we have obtained adsorption energies using the PBE-TS/IH method of Bučko and coworkers,Bučko et al. (2013) and the PBE-D3 and, for the finite cluster only, the PBE-D4 approaches of Grimme and coworkers.Grimme et al. (2010); Grimme, Ehrlich, and Goerigk (2011); Caldeweyher, Bannwarth, and Grimme (2017)

The adsorption energies are collected in Table 4 for bulk and in Table 5 for the 2T cluster. The PBE-D3 scheme substantially overestimates the adsorption energies of hydrocarbons, by around 16 kJ/mol for propane in bulk and by up to 3.4 kJ/mol for acetylene on the 2T cluster. Interestingly, the data for the 2T cluster show that the recent D4 correction does not bring any significant improvement compared to the D3 scheme. Moreover, while the use of hybrid functionals has been found to improve the description of structural properties of zeolites,Göltl and Hafner (2012); Cutini, Civalleri, and Ugliengo (2019) using the hybrid PBE0 functional does not improve the results substantially over those obtained with PBE. The only molecule where PBE-D3 underestimates the binding is CO2 and it’s worth mentioning that similar unexpected underbinding was also observed for CO2 crystal.Klimeš (2016) The PBE-TS/IH method was shown to improve adsorption energies for Cs-exchanged chabazite.Shang et al. (2014); Bučko et al. (2014) While PBE-TS/IH performs reasonably well for ethylene and acetylene, it considerably overestimates the adsorption energies of alkanes. The case of CO2 becomes even more puzzling as PBE-TS/IH underestimates its adsorption energy by some 6 kJ/mol for bulk and by $2.1$ kJ/mol for the 2T cluster. Overall, the results illustrate one of the benefits of RPA over dispersion corrected DFT which is the consistency of RPA errors. Finally, we note that “standard” PAW potentials give interaction energies which are artificially too strong for CO2 and H2O, by 0.9 and 1.5 kJ/mol, respectively, for alkanes the differences to the more precise “hard” PAW potentials are less than 0.2 kJ/mol.

V Discussion and Conclusions

We have shown that RPA with singles is highly accurate for predicting adsorption energies and we now briefly comment on comparison with experimental data, such as heats of adsorption which are available for some of the systems considered here, see, e.g. Refs. Hudson et al., 2012; Pham, Liu, and Lobo, 2013; Piccini et al., 2015; Hyla et al., 2019 for recent work. Our interaction energies are obtained for temperature of zero Kelvin and neglecting quantum effects of nuclei as well as relaxation of the monomers. The standard way to approximately treat the first two effects is to obtain zero-point energy and thermal corrections from vibrational frequencies obtained for molecule bound to the adsorption site. This approach was used by Piccini et al. who notably evaluated the vibrational frequencies using anharmonic potential energy surface.Piccini et al. (2015) However, at room temperature the positions of the molecules are not necessarily limited to a single adsorption site and molecular dynamics (MD) is a more appropriate tool to account for the thermal effects.Bučko et al. (2011) Given the system sizes and timescales required to converge adsorption energies in zeolites, performing MD based on RPA is not possibleRamberger, Schäfer, and Kresse (2017) and alternative strategies are needed. Rocca et al. have recently used a resampling approach to correct PBE-D adsorption enthalpies calculated from MD to RPA level.Rocca et al. (2019) Their RPA adsorption enthalpy overestimates the experimental data which is unexpected given that RPA should underestimate the binding, but we note that it is challenging to converge the resampling scheme. To overcome this issue, Chehaibou and co-workers used machine learning to first train the difference between RPA and PBE-D energies and then to use the trained model in the resampling.Chehaibou et al. (2019) This improves the convergence of the resampling significantly and leads to results which are consistent with the data presented here. In a similar spirit, Rubeš et al. resampled a PBE-based MD with PBE/CC to explain the temperature dependence of adsorption enthalpy of CO in ferrierite (H-FER).Rubeš et al. (2018) The PBE/CC was specifically reparametrized to reproduce interaction energies at CCSD(T) level for finite clusters and at RPA+RSE level for bulk. In this work the simulated and measured adsorption heats were within 1 kJ/mol. Despite these promising results, MD does not account for quantum nuclear effects which could be relevant even at room temperature for hydrogen-containing molecules. In any case, RPA with singles can be used to provide reference energies at 0 K with almost reference quality, i.e. with errors of less than 2 kJ/mol or 5 %. The values can be further used in resampling of MD or by schemes such as DFT/CC or otherFang et al. (2013) to obtain the finite temperature adsorption enthalpies as shown in Ref. Rubeš et al., 2018, possibly even with the inclusion of zero point energies, but this is beyond the scope of this paper.

While our focus here was on adsorption energies, we note that RPA with singles was shown to give very accurate results for structural properties of solids as well.Klimeš et al. (2015) Moreover, Xiao et al. found a very good agreement between RPA predictions and experiment for equilibrium volumes, bulk moduli, and relative stabilities of $\alpha$-quartz and stishovite phases of SiO2.Xiao et al. (2012) The accurate description of these properties for zeolites is also problematic for dispersion corrected DFT schemes.Fischer et al. (2019) The downside of RPA is its increased computational cost compared to standard DFT and one can ask to how large systems it can be applied. Based on our experience, systems with unit cell volumes up to approximately 3000 Å3 can be studied routinely. This means that apart from chabazite, zeolites such as ferrierite, mordenite, or LTA are accessible. As RPA does not rely on any parameters, it can be also readily applied to different cation-exchanged zeolites which show interesting adsorption properties and for which reference data would helpful.Zhang, Singh, and Webley (2008); Denayer et al. (2008); Shang et al. (2012, 2014); Coudert and Kohen (2017) Finally, RPA with singles can also provide adsorption energies and structural properties of nearly reference quality to other porous materials, such as MOFs and ZIFs.Wieme et al. (2018)

In summary, we have shown that the RPA with singles corrections scheme surpasses the MP2 approach for obtaining adsorption energies of molecules in zeolites. This is not only in terms of accuracy but also in terms of computational cost. We demonstrated this on a test set containing seven diverse molecules. The lower reliability of MP2 was the most apparent for ethylene and acetylene. These molecules possess delocalized electrons for which the description of electron screening is needed, this is present in RPA, but not in MP2. Overall, RPA with singles is a suitable scheme for obtaining highly accurate adsorption energies that can be used to compare to experimental data and to assess or even enhance the accuracy of simpler theoretical methods.

Acknowledgements.

JK acknowledges support from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 658705 and via ERC grant APES (No 759721) and support by a PRIMUS project from Charles University. DPT acknowledges support from the Max Planck Institute. Computational resources were provided by the MetaCentrum and CESNET (LM2015042), CERIT-SC (LM2015085), and by IT4Innovations (LM2015070), supported by the Czech Ministry of Education, Youth, and Sports. We thank Florian Göltl for providing us with some of the initial structures. We thank Miroslav Rubeš for useful discussions.