Application of a Dual-Hybrid Direct Random Phase Approximation to

Aug 8, 2016 - In water clusters, there is a delicate balance of van der Waals interactions and hydrogen bonds. Although semilocal and nonlocal density...
1 downloads 12 Views 1MB Size
Subscriber access provided by University of Birmingham

Article

Application of a Dual-Hybrid Direct Random Phase Approximation to Water Clusters Pál D. Mezei, Adrienn Ruzsinszky, and Gabor I. Csonka J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00323 • Publication Date (Web): 08 Aug 2016 Downloaded from http://pubs.acs.org on August 10, 2016

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

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

Page 1 of 19

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

Journal of Chemical Theory and Computation

Application of a Dual-Hybrid Direct Random Phase Approximation to Water Clusters Pál D. Mezei,1 Adrienn Ruzsinszky,2 and Gábor I. Csonka1 1

Department of Inorganic and Analytical Chemistry, Budapest University of Technology and Economics, H-1521 Budapest, Hungary 2

Department of Physics, Temple University, Philadelphia, PA 19122, USA

Abstract In water clusters, there is a delicate balance of van der Waals interactions and hydrogen bonds. Although semi-local and non-local density functional approximations have been recently routinely applied to water in various phases, the accurate description of hydrogen bonds remains a challenge. The most popular density functional approaches fail to predict the correct ordering of the energies of water clusters. To illustrate the required accuracy, the CCSD(T) complete basis set extrapolated dissociation energy difference between the two lowest energy hexamer structures is 0.06 kcal mol-1 per monomer. In this work, we assessed interaction energies in neutral and ionic water clusters with various density functionals with or without van der Waals correction. Generally, van der Waals approximations play a significant role in clusters with increasing size, while hybrid functionals improve the description of hydrogen bonds. Despite these general trends, none of the tested density functional approximations with or without van der Waals correction and exact exchange mixing can lead to a uniform performance for neutral and ionic water clusters. The recently constructed dual-hybrid dRPA75 approximation is a successful combination of exact and semi-local exchange, and non-local correlation in its energy, while utilizing a high fraction of exact exchange. We have shown that the dRPA75 method has a systematic error, which can be efficiently compensated by the aug-cc-pVTZ basis set for small and medium-sized water clusters.

Introduction Water shows large structural and functional variety in gas-phase clusters, liquid water, and solid ice.1 2 34 In this paper, we focus on the energetics of water clusters from two well established test sets (vide infra), which show interesting features, keeping in mind that the correct energetic order can potentially lead to correct structure. Our central interest here is to perform accurate and efficient calculations for the interaction energies of various water clusters. Cohesion in water clusters is a delicate interplay of exchange-repulsion, electrostatics, induction, and dispersion interactions. The local density functionals yield a serious overbinding error for water clusters, and large error for the autoprotolysis reaction.5 The semi-local density functionals (e.g. PBE) usually give mixed results for water clusters.1 6 In PBE, there is an error compensation between the two-body and three-body exchange-repulsion interaction terms. However, the overestimation of the induction is larger than the underestimation of the dispersion; therefore, no proper error compensation can be expected.7 The recent non-empirical meta-GGA-level MGGA_MS0 functional yields somewhat better results.8 It was observed that certain hybrid functionals (e.g. PBE0) give quite accurate interaction energies for water clusters due to the delicate error compensation between the underestimated dispersion and the overestimated induction energy components.9 Note that the uncertainty of the error compensations makes the relative energy order incorrect,10 and works only for small clusters, where the dispersion energy component is relatively small.11 Adding dispersion correction to these energies leads to exaggerated cohesion.12 The decomposition of the PBE0 interaction energies shows that the major source of error is the calculated three-body exchange-repulsion interaction energy term.13 Notice that the total three-body interaction energy is sensitive to the quality of the exchange term, and the exact HF exchange yields good three-body interaction energies while many common GGA, meta-GGA, and hybrid GGA functionals fail seriously in this respect. This is why HF 1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

energies can be efficiently corrected by the D3(BJ) pairwise dispersion correction or by MP2 correlation.1 The direct random phase approximation (dRPA)14 15 16 17 18 19 is an efficient way to obtain serious improvements upon the standard semi-local density functionals; however, it requires corrections due to its well understood flaws.20 21 22 23 24 25 Recently, we have constructed a new and efficient dual-hybrid dRPA method (dRPA75), and we have shown that in organic chemistry a widely accurate description of chemical properties might be provided by a 75% exact exchange mixing both in the calculation of the reference orbitals and the final energy.26 The justification of the observed good performance is that the large fraction of exact exchange in a PBE hybrid can eliminate the consequences of the many-electron self-interaction error (i.e. too diffuse, overlapping electron densities, overestimated exchange-repulsion, induction and charge transfer).27 28 Furthermore, the non-additivity and anisotropy of the dispersion interaction energy can be seamlessly and properly taken into account by the dRPA correlation. But notice also that dRPA systematically underestimates the dispersion energy at complete basis set (CBS) limit. In dRPA75, the systematic underestimation of the dispersion effects is compensated by a systematic overestimation of the intermolecular cohesion by the aug-cc-pVTZ basis set. The significance of all these effects in the interaction energy of water clusters suggests that the recent dualhybrid dRPA75 approximation with the aug-cc-pVTZ basis set might be suitable for the accurate computation of water cluster interaction energies. In this paper, we analyze how this method performs for water clusters up to (H2O)20. We also discuss the available theoretical results.29

Figure 1. Geometries of the prism-shaped, cage-like, book-shaped and cyclic water hexamers from the BEGDB database (prism: water6PR, cage: water6CA book: water6BK1, cyclic: water6CC in Table A1) The interaction energies of the clusters from dimer to decamer were calculated accurately on the MP2/CBS and CCSD(T)/CBS levels of theory.2 The most stable water dimer is the open form while the trimer, tetramer and pentamer clusters have cyclic structure.30 The (H2O)6 cluster represents a transition from the two-dimensional cyclic to the three-dimensional all-surface structures.31 The interaction energy of the low-lying hexamers such as prism, cage, book, and cyclic shown in Figure 1 varies in a 2 kcal mol-1 range.32 Among the hexamers, the prism structure was found to be the most stable,33 34 35 although the difference between the reference energies of the two lowest energy cage and prism hexamer structures is 0.06 kcal mol-1 per monomer. After zero-point energy correction the cage hexamer becomes more stable than the prism form.36 The two most stable octamers have cubic structure with D2d or S4 symmetry.37 38 The most stable decamer is the pentagonal prism.39 The structure of 2 ACS Paragon Plus Environment

Page 2 of 19

Page 3 of 19

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

Journal of Chemical Theory and Computation

medium-sized clusters with odd number of monomers such as heptamer, nonamer, and undecamer can be derived from the structure of the lower neighboring even numbered clusters with the distortion effect of the additional water molecule.40 41 42 At lower temperatures, the compact structures are preferred. However, at higher temperatures, the configurational entropy makes the open structures more favorable.43 The lowest energy dodecamer clusters have a three-layered cuboid structure with similar partial symmetries to the octamer. An additional monomer unit distorts this cuboid geometry.44 Similarly to the dodecamer, the lowest energy (H2O)15 cluster is a three-layered pentagonal prism. The (H2O)16 isomers can be built up from two pentagonal prisms sharing one square.45 The (H2O)17 cluster shows a second transition from the all-surface to the one-water centered, internally solvated structures.46 The most stable (H2O)17 cluster is a spherical isomer with an interior water molecule. Similarly, the lowest energy (H2O)19 and (H2O)21 clusters have also an internally solvated water molecule. Contrarily, the lowest energy (H2O)18 and (H2O)20 clusters have all-surface structure.47 The low-lying (H2O)20 isomers are in the edge-sharing, face-sharing, fused cubes, and dodecahedron forms (Figure 2). Among these (H2O)20 isomers, the most stable one is the edge-sharing form, while the least stable one is the dodecahedron.48 49 In contrast to the hexamers, the zero-point vibrational energy does not change the relative energy order, but the energy spacing is reduced between the edge-sharing and dodecahedron forms.50 Above 200 K, the configurational entropy makes the dodecahedron the most common form. 51 52 The (H2O)3053 and (H2O)40 clusters54 are also discussed in the literature. Notice that in large clusters, the dispersion interactions might contribute significantly to the interaction energies.55

Figure 2. Geometries of the edge-sharing, face-sharing, fused cubes and dodecahedron (H2O)20 clusters In positively charged clusters, the extra proton was found on the surface of the clusters in the hydronium ion form.56 The hydronium ion is stabilized by donating the maximal possible number of hydrogen bonds in the clusters. The pentameric cluster is the first in which a four-membered ring is closed. The lowest energy (H2O)9H+ cluster has a distorted cubic structure.57 Although above the (H2O)16H+ cluster size, the clathrate structure becomes favorable, increasing the temperature, the (H2O)17H+ cluster adopts a three-dimensional open structure. The lowest energy (H2O)20H+ cluster has a distorted dodecahedral structure,58 and that is quite different compared to the edge-sharing structure of the lowest energy (H2O)20 cluster. In negatively charged clusters, the excess electron density on the oxygen atom of the hydroxide ion is stabilized by accepting the maximal possible number of hydrogen bonds in the clusters. In the octamers, a solvation shell is formed around the hydroxide ion.59 60 3 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 4 of 19

Computational details We used two benchmarks to test the computational methods on the interactions of water clusters. The first test set contains 38 neutral water clusters from the dimer to the decamer in the BEGDB database with estimated CCSD(T)/CBS reference interaction energies calculated using RI-MP2/ADZ optimized geometries.2 We denote Dunning’s aug-cc-pVXZ basis sets by AXZ (X = D, T, Q, 5, 6) for the sake of brevity. Investigation of binding energies shows the energy difference between RI-MP2/ATZ and RIMP2/ATZ//RI-MP2/ADZ calculations is less than 0.1 kcal/mol for (H2O)n, n = 2-6. We use the average energy of the two water monomers in the dimer to estimate the average monomer energy for the larger clusters too. The second test set is the Water27 subset of the GMTKN30 database.3 5 This includes a set of fourteen neutral water clusters ((H2O)n, n = 2-6, 8, 20), five hydronium ion clusters (H3O+(H2O)n, n = 1-3, 6), seven hydroxide ion clusters (OH-(H2O)n, n = 1-6), and one autoionized water cluster (H3O+(H2O)6OH-). The reference energies are estimated CCSD(T)/CBS energies using MP2/ATZ optimized geometries except for the 20-membered clusters, where the heuristic estimated MP2/CBS reference energies are calculated using the MP2/ADZ optimized geometries.61 The dRPA gives useful results for intra- and intermolecular non-covalent interactions,62,63 for adsorption,64,65 for interlayer interactions,66 and for van der Waals bonded solids.67 We applied the dRPA method non-self-consistently inter alia with PBE68 or PBE069 orbitals (denoted as dRPA@PBE and dRPA@PBE0). The dRPA exchange-correlation energy is the sum of the exact exchange, 𝐸Xexact and dRPA correlation energies, 𝐸CdRPA: dRPA (1) 𝐸XC = 𝐸Xexact + 𝐸CdRPA The dRPA correlation energy can be obtained from the B non-antisymmetrized two-electron repulsion matrix and T double excitation amplitude matrix [eq (2)]. 1 (2) 𝐸CdRPA = tr[𝑩 𝑻] 2 The B matrix is given by the matrix elements: 𝐵𝑖𝑎,𝑗𝑏 = 〈𝑖𝑗|𝑎𝑏〉. The T matrix is given by an iterative update procedure [eq (3)] as a Hadamard product with ∆ and initialized with T(0) = 0:70 (3) 𝑻(𝑛+1) = −∆ ∘ (𝑩 + 𝑩𝑻(𝑛) + 𝑻(𝑛) 𝑩 + 𝑻(𝑛) 𝑩𝑻(𝑛) ) −1

The orbital denominator matrix is defined by the matrix elements ∆𝑖𝑎,𝑗𝑏 = (𝜀𝑎 + 𝜀𝑏 − 𝜀𝑖 − 𝜀𝑗 ) . To separate the effect of the orbitals for the exact exchange (EXx) and dRPA correlation (dRPAc), we compare the hybrid dRPAh = HF + (dRPAc@PBE) interaction energies71 with the dRPA@HF and dRPA@PBE interaction energies. We also computed the dual hybrid dRPA75 interaction energies.26 This method provides a more balanced description of non-covalent interactions. It gives better interaction energies for charge transfer complexes and hydrogen bonding than the standard non-self-consistent dRPA approaches. At the same time, it also gives excellent π-interaction energies and dipole-dipole interaction energies. Furthermore, it is good for weakly interacting complexes, reaction energies, and barrier heights. The dRPA75 energy is calculated as: dRPA75 (4) 𝐸XC = (0.75𝐸Xexact + 0.25𝐸XPBE + 𝐸CdRPA )@PBE0.75, where PBE0.75 symbolizes the self-consistent hybrid PBE functional calculated with 0.75 fraction of the exact exchange. The dRPA calculations were performed with aug-cc-pVTZ basis set (ATZ), with density fitting using the MRCC (09/07/2015) quantum chemistry software.72,73 We analyzed the effect of the basis set size on the interaction energies of water hexamers using AXZ; X = T, Q, 5 basis sets and performed a three-point CBS extrapolation using fitted power functions separately for the exchange and correlation energies similarly to our earlier work.74 We estimated the dRPA75/CBS energies of the icosameric clusters using the average CBS(3/4) extrapolation coefficients determined for the water hexamers. We also analyzed the effect of the basis set superposition error calculating the counterpoise corrections for the hexamers, and for the icosamers.75

4 ACS Paragon Plus Environment

Page 5 of 19

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

Journal of Chemical Theory and Computation

For comparison, we cite and discuss the M06L/ATZ76 results from the literature.3 Our earlier results have shown that the ATZ basis set is practically converged for DFT energy calculations. 77 We also computed the non-local VV10 dispersion corrected PBE/ATZ and revPBE/ATZ energies.78 This correction describes the pairwise instantaneous dipole – instantaneous dipole dispersion interactions, and it is based on the local polarizability model. The PBE-VV10 and revPBE-VV10 calculations were performed using the Orca 3.0 quantum chemistry software.79 The functional-dependent parameter of the correction was b = 6.2 in PBE-VV1068 80 and b = 3.7 in revPBE-VV10.81 80 Furthermore, we also computed the molecular mechanical TS82 corrected PBE and PBE0 energies. This pairwise dispersion correction uses the free interatomic C6,AA dispersion coefficients and applies the Hirshfeld partitioning of the electron density to obtain the effective C6,AA dispersion coefficients and the effective R0,A van der Waals radii to control the short-range damping. The interatomic C6,AB dispersion coefficients are calculated from the effective C6,AA and C6,BB dispersion coefficients and the atomic static dipole polarizabilities according to the London formula.83 84 The PBETS and PBE0-TS calculations were performed with the ATZ basis set using the FHI-AIMS program.85 We also included the many-body dispersion (MBD) energy correction which is computed using the coupled Quantum Harmonic Oscillator model Hamiltonian. This MBD energy contains the manybody energy terms that are missing from pairwise dispersion methods and the long-range Coulomb response contributions. The applied MBD methods were developed to work with semi-local density functionals like PBE or PBE0. In this paper, we calculate three versions of MBD: MBD@TS, MBD@SCS, and MBD@rsSCS86 with the ATZ basis set using the FHI-AIMS program,85 where the MBD@TS energy is calculated with polarizabilities without any screening based on the TS procedure; the MBD@SCS energy is calculated with polarizabilities screened with the full Coulomb potential; and the MBD@rsSCS is calculated with the range-separated self-consistent screening procedure. The MBD@rsSCS is currently the most advanced version of these MBD methods. All deviations from the reference are calculated as ΔE(method) – ΔE(reference). The accuracy and precision of the methods are characterized by the mean deviations (MD) and corrected sample standard deviations (CSSD), respectively. We also report the mean absolute deviations (MAD) for comparison with the literature data. Results and discussion Neutral water clusters First we discuss the performance of the various methods for the small and medium-sized neutral water clusters from the BEGDB database. The computed individual interaction energies are shown in Table A-1. The PBE/ATZ method overestimates the interaction energy per hydrogen bond (Table A-1). The application of the TS or VV10 dispersion corrections would increase the overestimation in the interaction energy by 0.52 or 0.85 kcal mol-1 per hydrogen bond, respectively. These corrections also have a 0.09 or 0.12 kcal mol-1 uncertainty per hydrogen bond, which would further decrease the precision of the method (CSSD = 0.16-0.18 kcal mol-1 per bond). This leads to the conclusion that the TS or VV10 corrected PBE method performs worse for the water clusters than the parent PBE functional. This is because the dispersion correction spoils the delicate error compensation as pointed out in the introduction. The hybridization with 0.25 fraction of the exact exchange shifts the average interaction energy per hydrogen bond towards the endothermic direction by 0.09 kcal mol-1 and also improve the precision (Table 1). But as the PBE0 hybrid shows an exothermic error as well, the TS or VV10 corrections would worsen the results again. For PBE0 again the dispersion correction spoils the results as pointed out earlier. The revPBE/ATZ method shows an increased deviation from the local density approximation (LDA) in its exchange part at high reduced gradient values, hence it systematically underestimates the interaction energy per hydrogen bond. This modified exchange eliminates the error compensation present in PBE and PBE0 and together with seriously underestimated dispersion effect it leads to serious underestimation of the cohesion. This can be remedied by the addition of a dispersion correction. The 5 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

revPBE-VV10 method is not particularly accurate for the small clusters; however, it works better for the medium-sized neutral clusters. The overall good performance of the revPBE-VV10 method suggests that the increased deviation from the LDA exchange and the non-local correlation are important for the proper description of these systems. Notice that the VV10 correction introduces some random error again. Table 1. Performance of different methods (in kcal mol-1 per hydrogen bond) on the neutral water clusters of BEGDB (with CCSD(T)/CBS reference interaction energies). (The positive deviations mean endothermic errors; the negative deviations mean exothermic errors. All calculations were performed with the ATZ basis set. The MP2/CBS results were taken from BEGDB.) PBE MDa MADb CSSDc a

-0.16 0.17 0.13

PBE0 -0.07 0.10 0.11

revPBEdRPA dRPA dRPA@ MP2/CBS dRPA75 VV10 @PBE @PBE0 PBE0.75 1.72 0.01 -0.02 0.29 0.01 -0.12 -0.02 1.72 0.07 0.03 0.29 0.03 0.12 0.04 0.06 0.10 0.03 0.06 0.04 0.04 0.04

revPBE

MD: mean deviation b MAD: mean absolute deviation c CSSD: corrected sample standard deviation

The MP2/CBS method shows good performance with a small exothermic shift (cf. Introduction). This suggests that the many-body dispersion terms missing from MP2 are practically negligible for the small neutral water clusters. However, for larger water clusters we observed inconsistencies between the various CBS extrapolated MP2 energies (see later). The dRPA@PBE/ATZ results show considerable average endothermic error with a quite good precision (Table 1). The endothermic error of the dRPA@PBE/ATZ method can be very well removed by replacing the PBE reference by the PBE0 reference. Comparison of the dRPA@PBE/ATZ, dRPA@PBE0/ATZ, and [email protected]/ATZ mean deviations reveals that the results are quite sensitive to the reference Slater determinants, and there is a non-linear dependence between the mean deviation and the fraction of exact exchange. Notice that the dRPA@(EXx + PBEc)/ATZ method gives smaller exothermic deviation than the [email protected] method (MD = -0.08 kcal mol-1 per bond), and the dRPA@HF/ATZ method yields again a slightly endothermic deviation (MD = 0.12 kcal mol1 per bond) (Table A-1b). Finally, the dRPA75 method also shows systematic errors. We discuss the basis set dependence of the dRPA interaction energies in the next section. The error distributions for the best performing methods discussed here can be followed in Figure 3. The dRPA75/ATZ errors are close to the normal distribution. The dRPA@PBE/ATZ, MP2/CBS, and revPBE-VV10/ATZ error distributions are more distorted.

6 ACS Paragon Plus Environment

Page 6 of 19

Page 7 of 19

dRPA75/aug-cc-pVTZ

10

10

8

8

Frequency

Frequency

dRPA@PBE0/aug-cc-pVTZ

6 4 2 0 -0.2

6 4 2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0 -0.2

0.2

-1

-0.15

MP2/CBS(3/4)

-0.05

0

0.05

0.1

0.15

0.2

revPBE-VV10/aug-cc-pVTZ 10

8

8

Frequency

10

6 4 2 0 -0.2

-0.1

Error per hydrogen bond (kcal mol-1)

Error per hydrogen bond (kcal mol )

Frequency

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

Journal of Chemical Theory and Computation

6 4 2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0 -0.2

0.2

-1

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

Error per hydrogen bond (kcal mol-1)

Error per hydrogen bond (kcal mol )

Figure 3. Histogram plots of error distributions for the dRPA@PBE0/ATZ, dRPA75/ATZ, MP2/CBS, and revPBE-VV10/ATZ methods on the interaction energies of the BEGDB water clusters. (We cut two instances from the right side of the revPBE-VV10 diagram.) Water hexamers Since water hexamers form a particularly important test for any theoretical method, we discuss here the performance of the studied methods on the interaction energies of four typical low-lying water hexamers from the BEGDB database (prism: water6PR, cage: water6CA book: water6BK1, cyclic: water6CC in Table A-1 and Figure 1). All dRPA methods predict the right energy order (prism < cage < book < cyclic) with the ATZ basis set. The dRPA75/ATZ method gives the best agreement with the reference CCSD(T)/CBS interaction energies followed by the dRPA@PBE0/ATZ method (Figure 4). It was shown earlier that the dRPA energies have a slow convergence with respect to the basis set size, and the CBS extrapolation requires particular care.74 Therefore, we have decomposed the dRPA75/ATZ errors to method error, as well as basis set incompleteness and superposition errors by complete basis set extrapolation and counterpoise correction. Table 2 shows that the dRPA75 method has a quite systematic, endothermic method error, which is comparable to the exothermic ATZ basis set error. According to Table A-2, the overbinding error of the ATZ basis set compensates the underbinding error of the dRPA75 and dRPA@PBE0 methods on the hexameric clusters. The largest deviation is 0.20 kcal mol-1 for dRPA75/ATZ and 0.43 kcal mol-1 for dRPA@PBE0/ATZ. Table 2. The dRPA75 method error (ME), as well as the ATZ basis set error (BSE), basis set incompleteness error (BSIE), and basis set superposition error (BSSE) for the low-lying water hexamers. (The errors were calculated relatively to the reference interaction energies.) prism cage book cyclic Avga CSSDb a

ME -8.8% -8.5% -7.5% -7.4% -8.0% 0.7%

BSE 8.8% 8.7% 8.0% 7.4% 8.2% 0.7%

BSIE 23.3% 23.2% 22.0% 20.5% 22.2% 1.3%

BSSE -14.5% -14.5% -14.0% -13.1% -14.0% 0.7%

Avg: average b CSSD: corrected sample standard deviation

7 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

In contrast to the dRPA75 and dRPA@PBE0 methods, note the quite poor performance of the hybrid dRPAh method with serious exothermic error (Table A-2). The application of the PBE orbitals instead of the HF orbitals in the calculation of the dRPA correlation energy results in a large exothermic shift in the interaction energies. While the same alteration in the exact exchange calculation results in a smaller endothermic shift in the interaction energies.

MP2/CBS

dRPA75

CCSD(T)/CBS

revPBE-VV10

-44.0

Interaction energy (kcal mol-1)

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

Page 8 of 19

-44.5 -45.0 -45.5 -46.0 -46.5 prism

cage

book

cyclic

Figure 4. Interaction energies of four water hexamers calculated with MP2/CBS (taken from BEGDB), dRPA75/ATZ, CCSD(T)/CBS (taken from BEGDB), and revPBE-VV10/ATZ methods. The MP2/CBS method underestimates the energy difference between isomers by ~0.2 kcal mol, which results in nearly equal interaction energies for the prism and cage isomers (Figure 4). The largest MP2/CBS deviation occurs for the cyclic isomer (-0.5 kcal mol-1). These conformers are also included in the Water27 test set with slightly different optimized geometries, which difference leads to only 0.1-0.2 kcal mol1 endothermic shifts in the CCSD(T)/CBS interaction energies. However, the MP2 interaction energies obtained with the aug-def2-TZVPP and aug-def2-QZVP basis sets and extrapolated to the CBS limit for the Water27 test set5 show a considerable, 0.6 kcal mol-1 overbinding deviations compared to the other MP2/CBS interaction energies.1 3 2 These deviations show that some of the MP2/CBS energies might be uncertain due to the extrapolation method. In Table A-2, we present the effect of the higher-order many-body dispersion terms for the water hexamers of the BEGDB database. Here we compare the deviations of the TS, MBD@TS, MBD@SCS or MBD@rsSCS corrected PBE/ATZ and PBE0/ATZ interaction energies from the reference. (Comparing the PBE0/ATZ and PBE0/A5Z results, the ATZ basis set has a 0.5 kcal mol -1 error here.1) The MBD@SCS and MBD@rsSCS corrected PBE or PBE0 methods generally provide an even larger exothermic shift of the interaction energies than the TS correction, thus they worsen further the results. The difference between the TS and MBD@TS corrections shows that the higher-order manybody terms are the greatest for the three-dimensional cage isomer, and the lowest for the twodimensional ring isomer (ring < book < prism < cage). The higher-order many-body terms destabilize the cage isomer relatively to the prism isomer by ~0.2 kcal mol-1. This agrees well with our previous observation that the MP2/CBS method underestimates the interaction energy difference between the 12

8 ACS Paragon Plus Environment

Page 9 of 19

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

Journal of Chemical Theory and Computation

prism and cage isomers by ~0.2 kcal mol-1. This suggests that the higher-order many-body dispersion terms are needed to obtain correct relative interaction energies of isomers with nearly degenerate energies. The revPBE-VV10 method also underestimates the interaction energy of the prism isomer by 0.4 kcal mol-1, but the interaction energies are correct for the other isomers. According to the literature, the M06L method also gives the right energy order.3 Notice that this order corresponds to the number of hydrogen bonds (ring < book < cage < prism – 6, 7, 8, 9). 20-membered water clusters Another interesting question is whether the previously mentioned methods work well also on larger clusters. For testing the performance of these methods, we use the (H 2O)20 clusters from the Water27 test set (Table A-3) with incremental CCSD(T) | MP2/CBS(4/5) reference interaction energies.49 We discuss the accuracy of the dRPA75 method, which is more robust for non-covalent interactions than the dRPA@PBE0 method as it was shown in our earlier paper.26 Among the (H2O)20 clusters from the Water27 test set, the dRPA75/ATZ method correctly predicts the edge-sharing isomer (-212.09 kcal mol-1) for the most stable and the dodecahedron (-202.18 kcal mol-1) for the least stable but in contrast to the reference incremental CCSD(T) | MP2/CBS(4/5) method, it gives a 0.7 kcal mol-1 spacing between the face-sharing (-210.33 kcal mol-1) and fused cubes (-209.65 kcal mol-1) isomers (cf. Figure 5 and Table A-3). Furthermore, the dRPA75/ATZ interaction energies are 3.4 kcal mol-1 larger in magnitude than the reference interaction energies. Table 3. The dRPA75 method error (ME), as well as the ATZ basis set error (BSE), basis set incompleteness error (BSIE), and basis set superposition error (BSSE) for the low-lying water icosamers. (The errors were calculated relatively to the reference interaction energies.) ME BSE BSIE BSSE edge-sharing -8.4% 10.0% 25.1% -15.1% face-sharing -8.6% 10.2% 25.4% -15.2% fused cubes -8.8% 10.1% 25.1% -15.0% dodecahedron -7.4% 9.7% 24.8% -15.1% Avg -8.3% 10.0% 25.1% -15.1% CSSD 0.6% 0.2% 0.2% 0.1% a

Avg: average b CSSD: corrected sample standard deviation

We have computed the dRPA75/CBS(3/4) energies using a two-point extrapolation scheme74 with the average coefficients obtained for the hexamers from the more accurate dRPA75/CBS(3/4/5) 𝑑𝑅𝑃𝐴75𝑥 𝑑𝑅𝑃𝐴75𝑐 calculations (𝐶3/4 = 0.380 𝑜𝑟 0.398 and 𝐶3/4 = 1.035 𝑜𝑟 1.047 for the cluster or monomer energies, respectively). Table 3 shows that the dRPA75 method error is quite systematic and increases with the number of hydrogen bonds as it was observed for the hexamers; however, the opposite sign basis set error overcompensates this for the larger (H2O)20 clusters. This overestimation comes mostly from the larger basis set incompleteness error of the ATZ basis set on the larger (H2O)20 clusters. In the literature, there are different extrapolated MP2/CBS interaction energies, and these show quite large deviation up to 6 kcal mol1. Analyzing the heuristic CBS extrapolation approach applied by Fanourgakis et al.,61 we concluded that this approach might lead to considerable extrapolation errors compared to the standard focal point approach, in which the HF energy and the correlation energy are extrapolated separately. Note that the MP2/CBS(3/4) energies by Goerigk and Grimme5 still have an extrapolation error compared to the more accurate MP2/CBS(4/5) interaction energies by Anacker and Friedrich.49 Also notice that the MP2/CBS energies of the face-sharing and fused cubes isomers are close to each other. The revPBE-VV10/ATZ method results in a face-sharing isomer 1.1 kcal mol-1 more stable than the fused cubes isomer (Figure 5). The revPBE-VV10/ATZ method shows larger exothermic errors 9 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

where the contribution of the dispersion interaction tends to be larger (fused cubes < face-sharing < edge-sharing < dodecahedron). We mention that the M06L relative interaction energies (see in ref 3) simply depend on the number of hydrogen bonds in the clusters (dodecahedron < edge-sharing < facesharing < fused cubes – 30, 35, 36, 37) thus this method predicts wrong energy order for the (H2O)20 structures.

Interaction energy (kcal mol-1)

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

Page 10 of 19

dRPA75/ATZ

MP2/CBS(Grimme)

revPBE-VV10/ATZ

iCCSD(T)|MP2/CBS(4/5)

-196 -198 -200 -202 -204 -206 -208 -210 -212 -214

Figure 5. Interaction energies of four (H2O)20 clusters calculated with different methods Replacing the erroneous heuristic MP2/CBS reference interaction energies61 by the incremental CCSD(T) | MP2/CBS(4/5) binding energies for the (H2O)20 clusters, changes the overall statistics of the standard MP2/CBS(3/4) interaction energies5 for the Water27 test set from {MD = 0.21 kcal mol-1, MAD = 1.39 kcal mol-1, CSSD = 1.93 kcal mol-1} to {MD = 1.14 kcal mol-1, MAD = 1.27 kcal mol-1, CSSD = 1.06 kcal mol-1}. This change of the reference energies in the Water27 test set also alters the statistics for the density functional methods tested on the GMTKN30 database. For example, the statistics of the PBE0 binding energies improves considerably as it changes from {MD = -0.79 kcal mol1 , MAD = 2.83 kcal mol-1, CSSD = 4.79 kcal mol-1} to {MD = 0.13 kcal mol-1, MAD = 2.08 kcal mol-1, CSSD = 2.84 kcal mol-1}. Also the statistics for MGGA_MS0 improves considerably leading to {MD = -0.05 kcal mol-1, MAD = 0.68 kcal mol-1, CSSD = 1.06 kcal mol-1}. Ionic clusters After we have analyzed the performance of the different methods on neutral water clusters, the question arises how the above discussed methods work for charged clusters in comparison to their performance on neutral clusters. In Table 4, we show the deviation statistics for revPBE-VV10/ATZ, MP2/CBS, and dRPA75/ATZ compared to the CCSD(T)/CBS interaction energies of the Water27 test set without the 20-membered clusters and the autoionized water cluster (Table A-3). This latter reaction will be discussed in the next section. 10 ACS Paragon Plus Environment

Page 11 of 19

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

Journal of Chemical Theory and Computation

The remaining 22 water clusters can be sorted into three subsets composed of seven negatively charged OH-(H2O)n, n = {1-6}; ten neutral (H2O)n, n = {2-6, 8}; and five positively charged H3O+(H2O)n, n = {1-3, 6} clusters. Inspection of the results in Table 4 reveals that the revPBE-VV10/ATZ method shows a quite uneven performance, and fails seriously for positively and negatively charged clusters, despite the reasonable performance for the neutral clusters. The MP2/CBS is considerably more consistent, and it systematically overestimates the interaction energies by about 0.75-1.06 kcal mol-1 for these three subsets. The results also show that the MP2/CBS performs worse for the ionic clusters than for the neutral clusters. The dRPA75/ATZ method performs the best, and it gives quite balanced results for all the three subsets. For comparison, we mention that the M06L/ATZ method is significantly less precise than the MP2/CBS method.3 The MGGA_MS0 performs better than MP2/CBS for charged water clusters and it yields close to dRPA75/ATZ quality results for water clusters. Table 4. Performance of different methods (in kcal mol-1) on the seven negatively charged, ten neutral and five positively charged water clusters from the Water27 test set (with CCSD(T)/CBS reference interaction energies). The 20-membered clusters and the autoprotolysis reaction are not included. (The positive deviations mean endothermic deviations; the negative deviations mean exothermic deviations.) Negatively charged

MDa MADb CSSDc a

Neutral

Positively charged

revPBErevPBErevPBEMP2 dRPA75 MP2 dRPA75 VV10 VV10 VV10 /CBS /ATZ /CBS /ATZ /ATZ /ATZ /ATZ -1.12 -1.04 -0.45 -0.19 -0.75 -0.25 -2.08 1.18 1.04 0.47 0.42 0.75 0.30 2.08 0.97 0.33 0.29 0.44 0.39 0.27 0.42

MP2 /CBS -1.06 1.06 0.20

dRPA75 /ATZ -0.76 0.76 0.15

MD: mean deviation b MAD: mean absolute deviation c CSSD: corrected sample standard deviation

Autoionization reaction For the exothermic reaction of the H3O+(H2O)6OH- leading to neutral (H2O)8 cluster, the reference CCSD(T)/CBS reaction energy is -28.5 kcal mol-1 as shown in the Table A-3. The revPBEVV10/ATZ and MP2/CBS methods yield large 36% and 6% endothermic deviations, respectively. Contrarily, the dRPA75/ATZ method shows a much better performance giving only a negligible 0.15 kcal mol-1 (0.5%) endothermic deviation. For comparison, the M06L and M06L-D3 methods give about 13% exothermic deviations.5 The MGGA_MS0 yields 3% endothermic deviation. Conclusion For the 38 neutral water clusters from the dimer to the decamer in BEGDB database, we observed that the PBE/ATZ method overestimates the interaction energy per hydrogen bond, and the application of the TS or VV10 dispersion corrections increases this overestimation. We also observed that the TS or VV10 corrections introduce some random error, thus they worsen the precision of the calculated reaction energies. The hybridization of PBE with 0.25 fraction of exact exchange leads to somewhat decreased overestimation of the interaction energy per hydrogen bond. But the application of the TS or VV10 corrections worsens the results again. The best performing PBE variant is revPBE-VV10 because the underestimation in the revPBE interaction energy is compensated by the VV10 correction. The revPBEVV10 method is not particularly accurate for the small clusters; however, it works better for the mediumsized neutral clusters. But it does not give sufficiently accurate results for water hexamers or icosamers, for ionic clusters, and for the autoionization reaction in the Water27 test set. For comparison, we also discussed the M06L and D3 corrected results, which showed mixed performance. Our results are in agreement with the earlier observations that many common density functionals fail for water clusters due to the improper error compensation between the exchange and dispersion errors. In functionals 11 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 12 of 19

where these opposite errors compensate each other a proper dispersion correction destroys the error compensation and leads to serious overbinding error. For the Water27 test set, quite different MP2/CBS interaction energies are available in the literature. We have shown that the earlier heuristic MP2/CBS extrapolation leads to significant errors especially for the lowest energy (H2O)20 clusters. We have used the incremental CCSD(T) | MP2/CBS(4/5) interaction energies as reference. We have shown that, if the more accurate reference energies are used for the (H2O)20 clusters, many functionals e.g. MGGA_MS0 show even better performance than the published results suggest. The standard MP2/CBS method shows good performance with a small exothermic shift for the cluster formation. As mentioned the total three-body interaction energy is sensitive to the quality of the exchange term and the exact HF exchange yields good results. The MP2 correlation yields quite good two-body interaction energy for water clusters, thus HF exchange combined with MP2 correlation leads to good results. This suggests that the many-body dispersion terms missing from MP2 are practically negligible for the formation of small neutral water clusters. However, the standard MP2/CBS method underestimates the energy differences among the water hexamers leading to nearly equal interaction energies for the prism and cage isomers and considerable exothermic error for the cyclic isomer. The analysis of the differences between the MBD@TS and TS corrected interaction energies suggests that the missing many-body dispersion terms can cause this inaccuracy in MP2 for the relative energies. Our results also show that the standard MP2/CBS results are worse for ionic clusters than for neutral clusters, and this method yields significant error (6%) for the autoionization reaction. We have also tested several dRPA variants as dRPA@PBE, dRPA@PBE0, [email protected], dRPA@(EXx + PBEc), dRPA@HF, the hybrid dRPAh, and the new dual-hybrid dRPA75. We used ATZ basis set, and also performed a CBS(3/4/5) extrapolation for the water hexamers. These results show that the dRPAh method gives very poor results, and all the other dRPA methods yield significant endothermic errors at the CBS limit for the formation energies of water clusters. The dRPA75 method shows quite systematic method errors (~8%) for small, medium-sized and large water clusters. For small and medium-sized clusters, the endothermic errors of the dRPA75 method (and similarly those of the dRPA@PBE0 method) can be efficiently compensated by the systematic exothermic error of the ATZ basis set. However, for larger clusters the ATZ basis set error tends to overcompensate the method error. In this paper, we have found that the dRPA75/ATZ results are in good agreement with the CCSD(T)/CBS reference for the water clusters of the BEGDB and Water27 test sets. The dRPA75/ATZ method also shows good performance for water hexamers, for charged clusters, and for the autoionization reaction of the water octamer cluster. The large, 0.75 fraction of the exact exchange in dRPA75 gives reasonable total three-body interaction energy. The origin of the good performance of the dRPA75/ATZ is that it also provides the proper balance between the exchange energy (the total three-body interaction energy is sensitive to the quality of the exchange term) and the many-body dispersion energy terms. Acknowledgement This work was supported by the Department of Energy under Grant No. DE-SC0010499. The authors thank the Hungarian National Information Infrastructure Development Institute for the computer time.

12 ACS Paragon Plus Environment

Page 13 of 19

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

Journal of Chemical Theory and Computation

Appendix Table A-1a Calculated interaction energies (kcal mol-1) for the water clusters of the BEGDB database System code water2Cs water3UUD water3UUU water4S4 water4Ci water4Py water5CYC water5CAA water5CAB water5CAC water5FRA water5FRB water5FRC water6BAG water6BK1 water6BK2 water6CA water6CB1 water6CB2 water6CC water6PR water7PR3 water7BI2 water7BI1 water7PR1 water7CH3 water7CH2 water7CA2 water7HM1 water7CA1 water7PR2 water7CH1 water8D2d water8S4 water9D2dDD water9S4DA water10PP1 water10PP2

H-bonds

CCSD(T) /CBS

1 3 3 4 4 5 5 7 7 7 6 6 6 7 7 7 8 6 6 6 9 11 8 8 11 8 8 9 8 9 10 8 12 12 13 13 15 15

-5.03 -15.70 -15.09 -27.43 -26.58 -23.88 -36.01 -34.54 -33.82 -34.69 -33.13 -34.88 -32.44 -44.59 -45.51 -45.14 -45.93 -43.57 -43.51 -44.60 -46.14 -56.89 -53.33 -53.40 -57.37 -53.08 -53.87 -54.88 -52.35 -55.69 -57.10 -54.38 -72.55 -72.56 -81.70 -81.46 -92.89 -92.89

PBE /ATZ -4.98 -15.81 -14.88 -28.65 -27.73 -23.57 -37.96 -34.33 -33.94 -34.72 -34.05 -35.55 -33.35 -46.11 -47.11 -46.89 -46.43 -45.77 -45.57 -46.79 -46.23 -57.53 -55.80 -55.83 -58.12 -54.95 -56.18 -56.18 -54.13 -57.03 -57.86 -56.75 -73.92 -73.93 -84.05 -83.75 -95.36 -95.46

PBEPBETS VV10 /ATZ /ATZ -5.34 -5.60 -16.76 -17.95 -15.82 -16.99 -30.56 -32.02 -29.64 -31.07 -25.78 -27.33 -40.23 -42.24 -37.69 -39.72 -37.31 -39.25 -38.12 -40.12 -36.92 -38.90 -38.38 -40.36 -36.17 -38.09 -50.12 -52.45 -50.79 -53.18 -50.64 -53.05 -51.18 -53.56 -48.45 -50.93 -48.29 -50.74 -49.43 -51.89 -50.89 -53.62 -63.25 -66.31 -59.85 -62.90 -59.96 -62.97 -63.94 -66.93 -59.10 -61.97 -60.38 -63.26 -61.63 -64.43 -57.60 -60.68 -62.42 -65.33 -63.72 -66.82 -60.87 -63.74 -81.82 -85.36 -81.81 -85.35 -92.28 -96.50 -92.05 -96.27 -105.17 -109.79 -105.11 -109.71

13 ACS Paragon Plus Environment

PBE0 /ATZ -5.05 -15.64 -14.83 -28.24 -27.35 -23.33 -37.44 -33.94 -33.41 -34.21 -33.36 -35.14 -32.71 -45.30 -46.42 -46.10 -45.69 -45.25 -45.14 -46.35 -45.49 -56.68 -54.65 -54.70 -57.22 -54.26 -55.28 -55.26 -53.57 -56.08 -56.94 -55.94 -72.66 -72.66 -82.44 -82.12 -93.80 -93.74

revPBE /ATZ -3.42 -10.67 -9.88 -20.55 -19.72 -15.44 -27.58 -22.81 -22.61 -23.19 -23.34 -24.67 -22.79 -32.27 -33.30 -33.10 -31.53 -33.29 -33.10 -34.20 -30.93 -39.22 -39.69 -39.70 -39.68 -38.97 -40.13 -38.78 -38.73 -39.50 -39.41 -40.63 -50.69 -50.70 -58.44 -58.13 -66.12 -66.19

revPBEVV10 /ATZ -4.68 -15.14 -14.23 -27.46 -26.56 -23.09 -36.32 -33.72 -33.32 -34.07 -33.14 -34.44 -32.38 -44.92 -45.53 -45.47 -45.88 -43.78 -43.60 -44.63 -45.76 -56.68 -54.01 -54.08 -57.21 -53.09 -54.32 -55.19 -52.03 -56.00 -57.29 -54.71 -73.45 -73.40 -83.05 -82.90 -94.53 -94.23

Journal of Chemical Theory and Computation

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

Table A-1b Calculated interaction energies (kcal mol-1) for the water clusters of the BEGDB database System code water2Cs water3UUD water3UUU water4S4 water4Ci water4Py water5CYC water5CAA water5CAB water5CAC water5FRA water5FRB water5FRC water6BAG water6BK1 water6BK2 water6CA water6CB1 water6CB2 water6CC water6PR water7PR3 water7BI2 water7BI1 water7PR1 water7CH3 water7CH2 water7CA2 water7HM1 water7CA1 water7PR2 water7CH1 water8D2d water8S4 water9D2dDD water9S4DA water10PP1 water10PP2 a

H-bonds 1 3 3 4 4 5 5 7 7 7 6 6 6 7 7 7 8 6 6 6 9 11 8 8 11 8 8 9 8 9 10 8 12 12 13 13 15 15

MP2 /CBS -5.00 -15.72 -15.17 -27.64 -26.80 -23.81 -36.38 -34.46 -33.65 -34.59 -33.11 -34.95 -32.42 -44.85 -45.77 -45.40 -45.98 -44.06 -43.99 -45.09 -46.02 -56.95 -53.63 -53.68 -57.41 -53.39 -54.23 -55.04 -52.70 -55.86 -57.09 -54.76 -72.60 -72.64 -81.92 -81.68 -93.06 -93.01

dRPA dRPA dRPA dRPA dRPA dRPA75 @PBE @PBE0 @PBE0.75 @PBE1.0 @HF /ATZ /ATZ /ATZ /ATZ /ATZa /ATZ -4.77 -5.07 -5.28 -5.29 -5.14 -5.09 -14.59 -15.41 -15.83 -15.77 -15.23 -15.46 -13.97 -14.77 -15.21 -15.17 -14.68 -14.80 -25.94 -27.22 -27.88 -27.74 -26.91 -27.46 -25.13 -26.40 -27.05 -26.92 -26.10 -26.66 -22.50 -23.66 -24.19 -24.08 -23.17 -23.64 -34.07 -35.76 -36.68 -36.52 -35.51 -36.14 -32.73 -34.32 -35.02 -34.83 -33.53 -34.32 -32.17 -33.75 -34.42 -34.21 -32.90 -33.80 -32.97 -34.55 -35.24 -35.04 -33.72 -34.59 -31.52 -33.07 -33.71 -33.47 -32.28 -33.24 -33.03 -34.64 -35.43 -35.26 -34.10 -34.81 -30.85 -32.39 -33.04 -32.82 -31.65 -32.58 -42.59 -44.57 -45.46 -45.16 -43.73 -44.90 -43.34 -45.39 -46.38 -46.12 -44.71 -45.71 -43.08 -45.12 -46.07 -45.79 -44.36 -45.46 -43.84 -45.87 -46.75 -46.47 -44.81 -46.01 -41.16 -43.22 -44.39 -44.23 -43.04 -43.67 -41.07 -43.15 -44.34 -44.20 -43.01 -43.58 -42.07 -44.17 -45.37 -45.23 -44.03 -44.60 -44.11 -46.13 -46.95 -46.66 -44.93 -46.15 -54.52 -56.97 -58.03 -57.67 -55.69 -57.11 -50.91 -53.30 -54.40 -54.05 -52.36 -53.73 -50.99 -53.37 -54.44 -54.08 -52.39 -53.79 -55.02 -57.48 -58.51 -58.12 -56.11 -57.62 -50.49 -52.91 -54.14 -53.87 -52.26 -53.31 -51.32 -53.74 -54.91 -54.60 -52.97 -54.19 -52.59 -55.00 -56.05 -55.68 -53.82 -55.30 -49.51 -51.93 -53.24 -53.03 -51.50 -52.30 -53.32 -55.76 -56.86 -56.51 -54.66 -56.05 -54.76 -57.22 -58.23 -57.84 -55.79 -57.37 -51.79 -54.25 -55.47 -55.18 -53.56 -54.71 -70.01 -73.00 -74.13 -73.57 -71.00 -73.19 -69.99 -72.97 -74.12 -73.56 -71.00 -73.18 -78.89 -82.27 -83.61 -82.99 -80.26 -82.59 -78.67 -82.02 -83.31 -82.68 -79.93 -82.35 -90.09 -93.87 -95.32 -94.58 -91.46 -94.22 -89.99 -93.79 -95.26 -94.54 -91.42 -94.15

We use the shorter [email protected] notation here for the dRPA@(EXx + PBEc) method.

14 ACS Paragon Plus Environment

Page 14 of 19

Page 15 of 19

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

Journal of Chemical Theory and Computation

Table A-2 Deviations from the reference CCSD(T)/CBS interaction energies (kcal mol-1) for the water hexamers of the BEGDB database Method CCSD(T)/CBS PBE-TS/ATZ PBE-MBD@TS/ATZ PBE-MBD@SCS/ATZ PBE-MBD@rsSCS/ATZ PBE0-TS/ATZ PBE0-MBD@TS/ATZ PBE0-MBD@SCS/ATZ PBE0-MBD@rsSCS/ATZ dRPA@HF/ATZ dRPA@HF/AQZ dRPA@HF/A5Z dRPA@HF/CBS a dRPAh/ATZ dRPAh/AQZ dRPAh/A5Z dRPAh/CBS a dRPA@PBE/ATZ dRPA@PBE/AQZ dRPA@PBE/A5Z dRPA@PBE/CBS a dRPA@PBE0/ATZ dRPA@PBE0/AQZ dRPA@PBE0/A5Z dRPA@PBE0/CBS a [email protected]/ATZ [email protected]/AQZ [email protected]/A5Z [email protected]/CBS a dRPA75/ATZ dRPA75/AQZ dRPA75/A5Z dRPA75/CBS a

prism -46.14 -4.75 -4.47 -4.99 -5.53 -3.27 -3.27 -3.77 -4.01 1.21 2.39 3.14 4.67 -8.19 -6.54 -5.60 -3.88 2.03 3.73 4.56 6.09 0.01 1.56 2.42 4.12 -0.81 0.63 1.54 3.41 -0.01 1.28 2.17 4.04

cage -45.93 -5.25 -4.78 -5.25 -5.79 -3.75 -3.55 -4.00 -4.25 1.12 2.24 2.97 4.53 -7.88 -6.31 -5.36 -3.55 2.09 3.72 4.54 6.17 0.06 1.54 2.45 4.41 -0.82 0.54 1.44 3.34 -0.08 1.15 2.02 3.91

a

book -45.51 -5.28 -5.21 -5.48 -6.06 -3.88 -3.98 -4.25 -4.57 0.80 1.77 2.43 3.87 -6.99 -5.57 -4.71 -3.06 2.17 3.64 4.40 5.88 0.12 1.44 2.28 4.15 -0.87 0.33 1.14 2.91 -0.20 0.87 1.66 3.42

cyclic -44.60 -4.83 -5.02 -5.15 -5.73 -3.79 -4.08 -4.22 -4.58 0.57 1.23 1.83 3.38 -5.55 -4.50 -3.70 -1.85 2.53 3.61 4.32 6.02 0.43 1.40 2.15 4.05 -0.77 0.39 0.82 1.08 0.00 0.75 1.45 3.03

The CBS extrapolation is based on the ATZ, AQZ, and A5Z energies. We used fitted power functions to estimate the convergence of the exchange and correlation energies.

15 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 16 of 19

Table A-3 Calculated interaction energies (kcal mol-1) for the Water27 test set est. dRPA75 MP2/CBSb CCSD(T)/CBSa /ATZ

# Cluster code 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

H2O2 H2O3 H2O4 H2O5 H2O6 H2O6c H2O6b H2O6c2 H2O8d2d H2O8s4 H2O20 H2O20fc H2O20fs H2O20es H3OpH2O H3OpH2O2 H3OpH2O3 H3OpH2O63d H3OpH2O62d OHmH2O OHmH2O2 OHmH2O3 OHmH2O4c4 OHmH2O4cs OHmH2O5 OHmH2O6 H3OpH2O6OHm

-5.01 -15.8 -27.4 -35.9 -46.0 -45.8 -45.3 -44.3 -72.6 -72.6 -197.7 -207.1 -207.1 -208.7 -33.5 -56.9 -76.5 -117.8 -114.9 -26.6 -48.4 -67.6 -84.8 -84.8 -100.7 -115.7 -28.5

-5.04 -15.99 -28.02 -36.79 -46.63 -46.54 -46.26 -45.44 -73.75 -73.72 -198.9 -206.3 -206.7 -208.7 -34.30 -57.86 -77.58 -118.93 -116.25 -27.67 -49.72 -68.58 -85.19 -86.24 -101.68 -116.79 -26.74

a

-5.06 -15.57 -27.50 -36.17 -46.15 -46.01 -45.70 -44.58 -73.29 -73.23 -202.18 -209.65 -210.33 -212.09 -34.29 -57.48 -77.12 -118.70 -115.80 -27.27 -49.21 -68.10 -84.75 -85.28 -100.93 -116.23 -28.35

revPBEVV10 /ATZ -4.71 -15.34 -27.62 -36.58 -45.58 -45.87 -45.77 -44.91 -73.19 -72.99 -204.07 -206.93 -208.05 -210.67 -36.19 -59.14 -78.20 -119.46 -116.99 -29.03 -50.43 -69.23 -84.61 -85.91 -100.84 -116.41 -18.35

incremental CCSD(T) | MP2/CBS(4/5) reference energies for the (H 2O)20 clusters (ref 49) and CCSD(T)/CBS reference energies for the rest (ref 61) b MP2/CBS(4/5) energies for the (H2O)20 clusters (ref 49) and MP2/CBS(3/4) energies by Goerigk and Grimme (ref 5) for the rest

References

16 ACS Paragon Plus Environment

Page 17 of 19

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

Journal of Chemical Theory and Computation

1

Santra, B.; Michaelides, A.; Fuchs, M.; Tkatchenko, A.; Filippi, C.; Scheffler, M. J. Chem. Phys. 2008, 129, 194111. 2 Temelso, B.; Archer, K. A.; Shields, G. C. J. Phys. Chem. A 2011, 115, 12034-12046. 3 Shields, R. M.; Temelso, B.; Archer, K. A.; Morrell, T. E.; Shields, G. C. J. Phys. Chem. A 2010, 114, 11725– 11737. 4 Bryantsev, V. S.; Diallo, M. S.; Van Duin, A. C. T.; Goddard, W. A. J. Chem. Theory Comput. 2009, 5, 1016– 1026. 5 Goerigk, L.; Grimme, S. Phys. Chem. Chem. Phys. 2011, 13, 6670–6688. 6 Morawietz, T.; Behler, J. J. Phys. Chem. A 2013, 117, 7356–7366. 7 Gillan, M. J. J. Chem. Phys. 2014, 141, 224106. 8 Hao, P.; Sun, J.; Xiao, B.; Ruzsinszky, A.; Csonka, G. I.; Tao, J.; Glindmeyer, S.; Perdew, J. P. J. Chem. Theory Comput. 2013, 9 (1), 355–363. 9 Svozil, D.; Jungwirth, P. J. Phys. Chem. A 2006, 110, 9194–9199. 10 Su, J. T.; Xu, X.; Goddard, W. A., III J. Phys. Chem. A 2004, 108, 10518. 11 Dahlke, E. E.; Truhlar, D. G. J. Phys. Chem. B 2005, 109, 15677–15683. 12 Klimeš, J.; Michaelides. J. Chem. Phys. 2012, 137 (12), 120901. 13 Wang, F.-F.; Jenness, G.; Al-Saidi, W. A.; Jordan, K. D. J. Chem. Phys. 2010, 132, 134303. 14 Langreth, D. C.; Perdew, J. P. Solid State Commun. 1975, 17, 1425–1429. 15 Langreth, D. C.; Perdew, J. P. Phys. Rev. B 1977, 15, 2884. 16 Langreth, D. C.; Perdew, J. P. Phys. Rev. B 1980, 21, 5469–5493. 17 Furche, F. J. Chem. Phys. 2008, 129, 114105. 18 Kubo, R. Rep. Prog. Phys. 2002, 29, 255–284. 19 Scuseria, G. E.; Henderson, T. M.; Sorensen, D. C. J. Chem. Phys. 2008, 129, 231101. 20 Mori-Sánchez, P.; Cohen, A. J.; Yang, W. Phys. Rev. Lett. 2008, 100, 146401. 21 Paier, J.; Janesko, B. G.; Henderson, T. M.; Scuseria, G. E.; Grüneis, A.; Kresse, G. J. Chem. Phys. 2010, 132, 094103. 22 Paier, J.; Janesko, B. G.; Henderson, T. M.; Scuseria, G. E.; Grüneis, A.; Kresse, G. J. Chem. Phys. 2010, 133, 179902. 23 Jiang, H.; Engel, E. J. Chem. Phys. 2007, 127, 1–11. 24 Furche, F. Phys. Rev. B 2001, 64, 195120. 25 Harl, J.; Schimka, L.; Kresse, G. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 115126. 26 Mezei, P. D.; Csonka, G. I.; Ruzsinszky, A.; Kállay, M. J. Chem. Theory Comput. 2015, 11, 4615–4626. 27 Sai, N.; Barbara, P. F.; Leung, K. Phys. Rev. Lett. 2011, 106, 226403. 28 Atalla, V.; Yoon, M.; Caruso, F.; Rinke, P.; Scheffler, M. Phys. Rev. B 2013, 88, 165122. 29 Gadre, S. R.; Yeole, S. D.; Sahu, N. Chem. Rev. 2014, 114, 12132–12173. 30 Kistenmacher, H.; Lie, G. C.; Popkie, H.; Clementi, E. J. Chem. Phys. 1974, 61, 546. 31 Saykally, R. J.; Wales, D. J. Science 2012, 336, 814. 32 Xantheas, S. S.; Burnham, C. J.; Harrison, R. J. J. Chem. Phys. 2002, 116, 1493. 33 Xantheas, S. S. Chem. Phys. 2000, 258, 225. 34 Kim, K.; Jordan, K. D.; Zwier, T. S. J. Am. Chem. Soc. 1994, 116, 11568. 35 Krishnan, P. N.; Jensen, J. O.; Burke, L. A. Chem. Phys. Lett. 1994, 217, 311. 36 Liu, K.; Brown, M. G.; Saykally, R. J. J. Phys. Chem. A 1997, 101, 8995–9010. 37 Xantheas, S. S.; Aprà, E. J. Chem. Phys. 2004, 120, 823. 38 Maeda, S.; Ohno, K. J. Phys. Chem. A 2007, 111, 4527. 39 Shields, R. M.; Temelso, B.; Archer, K. A.; Morrel, T. E.; Shields, G. C. J. Phys. Chem. A 2010, 114, 11725. 40 Kim, J.; Majumdar, D.; Lee, H. M.; Kim, K. S. J. Chem. Phys. 1999, 110, 9128. 41 Jensen, J. O.; Krishnan, P. N.; Burke, L. A. Chem. Phys. Lett. 1996, 260, 499. 42 Lee, H. M.; Suh, S. B.; Kim, K. S. J. Chem. Phys. 2001, 114, 10749. 43 Acelas, N.; Hincapié, G.; Guerra, D.; David, J.; Restrepo, A. J. Chem. Phys. 2013, 139, 044310. 44 Bulusu, S.; Yoo, S.; Aprà, E.; Xantheas, S. S.; Zeng, X. C. J. Phys. Chem. A 2006, 110, 11781. 45 Leverentz, H. R.; Qi, H. W.; Truhlar, D. G. J. Chem. Theory Comput. 2013, 9, 995. 46 Lagutschenkov, A.; Fanourgakis, G. S.; Niedner-Schatteburg, G.; Xantheas, S. S. J. Chem. Phys. 2005, 122, 194310. 47 Yoo, S.; Aprà, E.; Zeng, X. C.; Xantheas, S. S. J. Phys. Chem. Lett. 2010, 1, 3122. 48 Wales, D. J.; Hodges, M. P. Chem. Phys. Lett. 1998, 286, 65. 49 Anacker, T.; Friedrich, J. J. Comput. Chem. 2014, 35, 634–643. 50 Fanourgakis, G. S.; Aprà, E.; de Jong, W. A.; Xantheas, S. S. J. Chem. Phys. 2005, 122, 134304. 51 Lenz, A.; Ojamäe, L. J. Chem. Phys. 2009, 131, 134302. 52 Parkkinen, P.; Riikonen, S.; Halonen, L. J. Phys. Chem. A 2013, 117, 9985–9998.

17 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

53

Page 18 of 19

Qian, P.; Song, W.; Lu, L.; Yang, Z. Int. J. Quantum Chem. 2010, 110, 1923. Liu, J.; Wang, L.; Zhao, J. J. Comput. Theor. Nanosci. 2009, 6, 454. 55 Iwata, S.; Bandyopadhyay, P.; Xantheas, S. S. J. Phys. Chem. A 2013, 117, 6641. 56 Searcy, J. Q.; Fenn, J. B. J. Chem. Phys. 1974, 61, 5282. 57 Bankura, A.; Chandra, A. Chem. Phys. 2011, 387, 92. 58 Parkkinen, P.; Riikonen, S.; Halonen, L. J. Phys. Chem. A 2012, 116, 10826. 59 Novakovskaya, Y. V.; Stepanov, N. F. J. Phys. Chem. A 1999, 103, 3285. 60 Novakovskaya, Y. V.; Stepanov, N. F. J. Phys. Chem. A 1999, 103, 10975. 61 Fanourgakis, G. S.; Aprà, E.; Xantheas, S. S. J. Chem. Phys. 2004, 121, 2655. 62 Eshuis, H.; Furche, F. J. Phys. Chem. Lett. 2011, 2, 983–989. 63 Zhu, W.; Toulouse, J.; Savin, A.; Angyán, J. G. J. Chem. Phys. 2010, 132, 244108. 64 Ren, X.; Rinke, P.; Scheffler, M. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 1–8. 65 Schimka, L.; Harl, J.; Stroppa, A.; Grüneis, A.; Marsman, M.; Mittendorfer, F.; Kresse, G. Nat. Mater. 2010, 9, 741–744. 66 Lebègue, S.; Harl, J.; Gould, T.; Ángyán, J. G.; Kresse, G.; Dobson, J. F. Phys. Rev. Lett. 2010, 105, 196401. 67 Harl, J.; Kresse, G. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 1–8. 68 Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. 69 Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. 70 Heßelmann, A. Phys. Rev. A 2012, 85, 012517. 71 Ren, X.; Tkatchenko, A.; Rinke, P.; Scheffler, M. Phys. Rev. Lett. 2011, 106, 153003. 72 MRCC, a quantum chemical program suite written by M. Kállay, Z. Rolik, I. Ladjánszki, L. Szegedy, B. Ladóczki, J. Csontos, and B. Kornis. See also Z. Rolik, L. Szegedy, I. Ladjánszki, B. Ladóczki, and M. Kállay, J. Chem. Phys. 139, 094105 (2013), as well as: www.mrcc.hu. 73 Kállay, M. J. Chem. Phys. 2015, 142 (20), 204105. 74 Mezei, P. D.; Csonka, G. I.; Ruzsinszky, A. J. Chem. Theory Comput. 2015, 11, 3961–3967. 75 Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553–566. 76 Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215–241. 77 Csonka, G. I.; Ruzsinszky, A.; Perdew, J. P. J. Phys. Chem. B 2005, 109, 21471–21475. 78 Vydrov, O. A.; Van Voorhis, T. J. Chem. Phys. 2010, 133, 244103. 79 Neese, F. (2012) The ORCA program system, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2, 73-78. 80 Aragó, J.; Ortí, E.; Sancho-García, J. C. J. Chem. Theory Comput. 2013, 9, 3437–3443. 81 Zhang, Y.; Yang, W. Phys. Rev. Lett. 1998, 80, 890−890. 82 Tkatchenko, A.; Scheffler, M. Phys. Rev. Lett. 2009, 102, 073005. 83 F. London Z. Physik. Chem. 1930, 11, 222. 84 K. T. Tang Phys. Rev. 1969, 177, 108. 85 Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab Initio Molecular Simulations with Numeric Atom-Centered Orbitals. Comp. Phys. Comm. 2009, 180, 2175-2196 86 Ambrosetti, A.; Reilly, A. M.; DiStasio, R. A.; Tkatchenko, A. J. Chem. Phys. 2014, 140, 18A508. 54

18 ACS Paragon Plus Environment

Page 19 of 19

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

Journal of Chemical Theory and Computation

TOC graphic 85x47mm (300 x 300 DPI)

ACS Paragon Plus Environment