Interaction between water and carbon nanostructures: How good are current density functional approximations?

The interaction of molecules or liquids with nanostructured surfaces is central to many real-life applications, including catalysis, gas storage, desalination, and more. Interfaces involving water and carbon show unique and fascinating behavior, which can be employed in important applications, for instance, in water purification devices.^1–10 Topologically similar materials can have substantially different properties^11–15 emphasizing that the understanding of the nature of the interaction has to be sought at a quantum mechanical electronic structure level.

Density functional theory (DFT) is the simulation method of choice for many materials applications due to its favorable accuracy to computational cost ratio.^16–19 Modern density functional approximations (DFAs) combine semilocal expansions of the exchange-correlation with long-range corrections for missing London dispersion interactions, i.e., the attractive part of the van der Waals (vdW) forces. The approximations are physically motivated but additionally require adjustment of a small number of parameters, which are either based on exact constraints or on empirical data. Adjusting the parameters to optimally describe short, long, and middle-ranges of interactions is challenging,^20,21 and a solution that is good for interactions between small molecules is not necessarily good for the interaction of a molecule with extended surfaces or within the bulk (e.g., molecules in solutions and molecular crystals).^19,22 In either case, it is mandatory to carefully benchmark the DFT methods, especially as the “zoo” of methodologies is growing and it is often unclear what is the expected reliability of a possible DFT setup or how to pick the best DFT flavor for a specific application. In the past decade, these DFA benchmarks mainly focused on molecular properties with recent studies testing more than 200 DFAs on thousands of references including thermochemistry, kinetics, and noncovalent interactions.^23–27 Recently, some focus of DFT benchmarking moved to the description of equilibrium geometries.^28–31 Similar large-scale benchmarks for condensed phase properties are much more rare, which is mainly due to the lack of theoretical reference data. While for bulk solids, experimental lattice constants and cohesive energies have been used successfully for DFT benchmarking,^32–35 similar data for more weakly bound molecular crystals have substantially higher uncertainty.^36–38 This is due to the experimental measurement uncertainty,³⁹ indirect measurements that cannot directly be compared to simple equilibrium geometries and energies, or the challenge to do the measurement itself. The latter point holds for a single water molecule adsorbed at surfaces as water readily forms clusters.⁴⁰ 

Concerning theoretical reference calculations, exciting progress has been made in the field of high-level wavefunction methods. On the one hand, embedding techniques^41–43 and local approaches of coupled cluster theories^44–49 have made the gold standard of quantum chemistry applicable to molecular systems with a few hundred atoms and molecular crystals of small molecules.^50–52 On the other hand, new algorithmic developments in the field of diffusion Monte Carlo [DMC, a quantum Monte Carlo (QMC) technique] have made the computation of chemically accurate lattice energies of small molecular crystals feasible within reasonable computational effort.^53,54 These developments together with an increased capacity of available computational resources make the interaction energy determination of extended systems feasible. This is an important step toward the better understanding of large noncovalently interacting systems as recently highlighted.^19,55,56

Here, we capitalize on this by gathering the benchmark quality interaction energies of water with carbon nanomaterials that we have studied in the past few years. This involves the adsorption of water on graphene,⁵⁷ water on benzene and coronene as two representative aromatic hydrocarbons (abbreviated as AH),⁵⁷ and water on a carbon nanotube (CNT).⁵⁸ In practical applications, it is important to describe correctly also the interaction between water molecules; thus, we additionally analyze different phases of two-dimensional (2D) ice⁵⁹ as well as bulk (3D) ice polymorphs.⁵⁴ Therefore, we obtain a dataset of 18 configurations and associated reference interaction energies, dubbed WaC18 set. We test as many as 28 DFAs on the WaC18 set, including several recently developed vdW corrections to DFAs. Some DFT benchmarks already exist on these or related system types,^58–76 but they typically include a limited set of DFAs, the recent vdW developments are not included, and the used reference data are not equally well converged. Here, we will address all of these issues.

The WaC18 benchmark test will help in understanding the essential ingredients needed to describe seamlessly both the strong hydrogen bonds and the weak interaction with surfaces. Furthermore, the identification of the most accurate DFAs can be used by researchers aiming to describe these widely spread system types, complementing and updating the perspective in Ref. 77 that focused on DFT recommendations for water.

In Sec. II, we describe the benchmark systems considered, discuss the best estimates of the interaction energies including additional DMC calculations to have equally well converged reference data, and give the computational details of the DFT calculations. Following this, we report the results of a variety of DFAs and vdW corrections from several functional classes and analyze the critical aspects determining the DFA-vdW performance (Sec. III). Conclusions and a future perspective are given in Sec. IV.

In our current benchmark, we will focus on water interacting with carbon nanostructures and ice. We separate the analysis into a single water molecule interacting with graphene, a carbon nanotube (CNT), aromatic hydrocarbons (AHs), and the interaction of solid water in three-dimensional and two-dimensional ice (see Fig. 1). Water adsorption on graphene and on the AHs is tested with three different water orientations, dubbed 0-leg, 1-leg, and 2-leg. Adsorption on the CNT is considered outside the CNT in a 2-leg configuration (external) and inside the CNT in a 2-leg configuration (internal). The two-dimensional ice polymorphs have been constructed in Ref. 59 by confining a single water layer resulting in four stable polymorphs of hexagonal, pentagonal, square, and rhombic ice structures. The three bulk ice phases cover the subtle balance of the competing polymorphs Ih, II, and the high-pressure form VIII. Overall, this benchmark dataset has been designed to investigate both the water surface interaction with different adsorption motifs, different surface sizes and curvature, as well as the transferability to many-water systems at different pressures and confinements.

For all systems under consideration in our study, we use the interaction energy at fixed equilibrium geometry, where we have well converged theoretical reference energies available. We provide the reference geometries and energies as the supplementary material to make our compiled benchmark easily accessible to other researchers.

The interaction energies considered are defined as the difference between a bound and an unbound configuration. Adsorption on nanostructure is computed as

where E_w@nano is in the bound configuration and the geometries of the individual fragments E_w and E_nano for the unbound configuration are kept frozen.⁷⁸ The interaction energies for 2D and 3D ice are given per molecule as

The reference interaction energies used in the following analysis are obtained from studies published in the past four years. All the reference values are reported in Table I. The reference values for the two-dimensional ice crystals are taken from Ref. 59, while the reference values for the three bulk-ice crystals are given in Ref. 54. The reference values for water on benzene, coronene, and graphene with different orientations are taken from Ref. 57. The reference interaction energy between water and CNT where previously investigated in Ref. 58, but here we report new values that are as tightly converged as the other references.

The computational approaches used to obtain the reference values are coupled cluster theory including singles, doubles, and perturbative triples [CCSD(T)] and fixed-node DMC. In particular, the values for water on benzene and on coronene are from CCSD(T) (where DMC yields identical results within the stochastic error), and all the other values are from fixed-node DMC. Indeed, most of these systems are very challenging or even out-of-reach for CCSD(T) due to their large size and the presence of periodic boundary conditions. With DMC, it is easier to assess the binding energy in large complexes because it is straightforward to simulate periodic systems, and DMC has favorable scaling with system size. In terms of accuracy, there is generally very good agreement between CCSD(T) and DMC in the evaluation of noncovalent interactions, provided that care is taken to ensure that an accurate setup is used for each method.^53,56,79 Agreement between methods is not expected beyond a given precision. To this aim, in Table I, we report the estimated uncertainty associated with any evaluation. For the water on benzene and coronene systems, where both CCSD(T) and DMC are affordable, there is excellent agreement among them.⁵⁷ 

Although reported DMC results are coming from different studies, the setup is consistent among them. DMC simulations were carried out with the CASINO code.⁸² Dirac-Fock pseudopotentials^83,84 with the localization approximation⁸⁵ (LA) are used. The trial wavefunctions were of the Slater-Jastrow type with single Slater determinants and the single particle orbitals obtained from DFT-local density approximation (LDA) calculations performed with the plane-wave self consistent field program (PWscf)⁸⁶ and re-expanded in terms of B-splines.⁸⁷ The Jastrow factor included electron-electron, electron-nucleus, and electron-electron-nucleus terms. The parameters of the Jastrow factor were carefully optimized by minimizing the variance, within a variational QMC scheme. The recently developed size-consistent DMC algorithm (ZSGMA)⁵³ was used. Finite time-step errors are carefully minimized by performing simulations with different values of the time step, until the bias appears safely smaller than the stochastic error. In periodic calculations, finite-size corrections are applied either using the model periodic Coulomb interaction^88–90 or the Kwee, Zhang, and Krakauer⁹¹ approach.

We reevaluate here the binding energy of water on the CNT because one specific aspect of the setup used in Ref. 58 is now known to be a possible issue: the optimization of the Jastrow factor for the configuration of water inside the CNT was not optimal. Since in the standard LA approach the Jastrow factor plays a major role in the pseudopotential error, we have developed a new method (to be reported elsewhere⁹²) that removes the LA bias. In repeating the evaluation, we also tuned some other aspect of the setup to be in line with the actual state-of-the-art. In particular, we used the recently developed eCEPP pseudopotentials⁹³ and the ZSGMA algorithm.⁵³ The reported binding values are obtained with time step τ = 0.01 a.u., which gives errors less than the stochastic uncertainty (standard deviation of 18 meV).

DFT test calculations based on the Perdew-Burke-Ernzerhof (PBE) functional⁹⁴ were done on selected systems to ensure convergence of all relevant numerical settings. Several different electronic structure codes and orbital basis expansions have been employed: ORCA 4⁹⁵ with large aug-cc-pVQZ, aug-cc-pV5Z,^96,97 and def2-QZVPPD⁹⁸ basis sets; CRYSTAL17^99,100 with a def2-QZVPPD(-f) basis; and Vienna Ab initio Simulation Package (VASP 5.4)^101,102 with projector-augmented plane waves (hard PAWs^103,104) with an energy cutoff of 1000 eV. In all codes, tight self-consistent field settings and large integration (and fine FFT) grids are used. The Brillouin zone sampling has been increased to converge the interaction energy to 1 meV. For the adsorption on graphene (5 × 5 supercell with 53 atoms in the unit cell) and the CNT [nonmetallic (10,0) configuration with 83 atoms in the unit cell], this reduces to a Γ-point calculation. For all PAW calculations, the nonperiodic directions use a vacuum spacing of 20 Å for the absorbed geometries and the same unit cell is used for the individual fragments, which compensates possible remaining image interactions for the binding energies. CRYSTAL and ORCA use open conditions in the nonperiodic directions consistent with the reference calculations.

For the water at the AH system, we established that the PBE interaction energy is converged within 2 meV using the def2-QZVPPD basis and the codes ORCA and CRYSTAL yield results within 1 meV. We additionally compared the PBE/def2-QZVPPD interaction energy with the PBE/hard PAW/1000 eV ones for water at graphene, water at AH, and ice Ih with a maximum deviation of 2 meV. To reach the DFT convergence for these systems, unusually tight thresholds are needed. In Table II, we list the convergence of the PBE ice Ih lattice energy for three complementary basis set expansions and software codes. For the PAW code, VASP, hard PAWs (i.e., small potential core) are needed as well as a minimum PW cutoff of 700 eV. CRYSTAL employs an all-electron basis set with atom-centered functions. Here, the interaction energy is converged employing neither a counterpoise-corrected def2-TZVPP nor a def2-QZVPP calculation. For fully converged values, a counterpoise-corrected def2-QZVPP calculation is needed and an extrapolation to the basis set limit reduced the binding by only 1.2 meV. Production level calculations for all reported DFT interactions on the full WaC18 benchmark are performed with VASP 5.4 using hard PAWs and a PW energy cutoff of 1000 eV.

DFAs from several rungs are tested: local density approximation (LDA), generalized gradient approximation (GGA), meta-GGA, and hybrid functionals (incorporation of nonlocal Fock exchange). The semilocal DFAs are corrected for missing long-range correlation effects by means of a variety of different semiclassical and nonlocal density based corrections (D2,¹⁰⁸ D3,^109,110 D4,^111,112 TS,¹¹³ MBD,¹¹⁴ VV10,¹¹⁵ dDsC,¹¹⁶ and vdW-DF¹¹⁷). TS and MBD are used with the noniterative Hirshfeld partitioning, D3 is used in the rational damping variant including Axilrod-Teller-Muto type three-body contributions, for D4, partitioned charges are generated by the electronegativity equilibration procedure (EEQ), and many-body contributions are covered by a standard coupled fluctuating dipole expression retaining an RPA-like expression. See Refs. 118–121 for further overview on vdW interactions in the DFT framework.

Using the DMC and CCSD(T) high-quality interaction energies, we can now test the capability of standard and new DFAs and a few simplified electronic structure methods for water adsorption on nanostructured surfaces and within ice polymorphs. While discussing the performances of each method, it is important to keep in mind the numerical DFT uncertainty of about 2 meV as well as the reference uncertainty, ranging from 1 meV (w@benzene, 0-leg) to 18 meV (w@CNT, external).

Table III summarizes the individual interaction energies from several electronic structure approaches and gives a root-mean-square (rms) error over the full WaC18 set. Additional statistical information on the performance of the methods is given in Table IV.

From looking at Tables III and IV, it is clear that Hartree-Fock (HF) misses all (Coulomb) correlation effects and cannot describe any of these noncovalently bound systems appropriately.^139,140 All systems are systematically underbound by up to 342 meV. While there is some weak binding for water on benzene, this diminishes for larger substrates and becomes repulsive for the adsorption on graphene and on the CNT. Thus, exchange repulsion, induction, and electrostatics are not sufficient to lead to a net binding of water on the carbon nanostructures considered. This is consistent with our symmetry adapted perturbation theory analysis in Ref. 57 as well as many studies on molecular dimers.¹⁴¹ That the local density approximation (LDA) yields an inconsistent description of small vdW complexes has been known since the mid-1990s.¹⁴² This is confirmed here, where LDA systematically overbinds all systems, in particular, the ice polymorphs. The LDA results for water adsorption, on the other hand, seem reasonably good. However, this is a fortuitous event due to the fixed geometries. While all other DFAs give equilibrium adsorption distances within 0.1–0.2 Å (see Refs. 22, 58, and 60), the LDA minimum is at substantially smaller distance and cannot be recommended for either geometry or stability estimates of vdW bound systems.

Including semilocal exchange-correlation effects as in the popular PBE GGA functional improves the behavior although most systems are bound a bit too weakly. Clearly, the long-range correlation effects leading to vdW attraction are missing. In the past decade, several methods have been developed for including these missing interactions (see, e.g., Refs. 118, 119, and 121). We test a broad range of these vdW corrections coupled with PBE (see Tables III and IV). The error spread is still substantial, and, in particular, the older effective pair-wise schemes (PBE-VV10, PBE-dDsC, PBE-TS, and PBE-D2) do not perform well. Recent vdW developments pay off, and we can see a clear improvement of PBE-MBD¹¹⁴ and PBE-D4¹¹² over their predecessors. The many-body vdW contributions decrease the binding yielding better agreement with the references. Most of this effect is already covered at the Axilrod-Teller-Muto type three-body level^143,144 as included in the D3 method.¹⁰⁹ At the PBE-vdW level, only D4 and VV10 are able to reproduce the relative stability of the water adsorption, i.e., 0-leg vs 2-leg and 1-leg vs 2-leg stability, to good accuracy (coming within the error of the reference energy). The best PBE based method (PBE-D4) yields an excellent description of the water adsorption with mean absolute deviation (MAD) from the reference of 13 meV. The description of the ice polymorphs is less satisfactory, which can be traced back to the intrinsic overpolarization and thus overestimated induction interaction of the PBE functional.^145–149 For instance, uncorrected PBE already overbinds hexagonal ice Ih by 24 meV, which clearly is not corrected by a (mostly) attractive vdW interaction. Overall, PBE-vdW does not perform well for strong hydrogen bonded systems.⁷⁷ 

DFAs with a nonlocal kernel to describe vdW interactions have been pioneered by Dion et al. (vdW-DF1).^117,150 This first-generation nonlocal functional gives unsatisfactory results on our benchmarks, the adsorption strength is overestimated, and the ice lattice energy is underestimated. The revised version with optimized semilocal exchange-correlation optB86b-vdW improves upon this, but the overall MAD is at 44 meV still rather high. Binding to graphene and the CNT seems to be extremely challenging for the nonlocal functionals with maximum deviation of 168 meV for water inside the CNT, as already noted in Ref. 58. Consistent with previous studies,^130,131 the second generation of nonlocal functionals significantly improves the results on all benchmark systems and the revised variant rev-vdW-DF2 more than halves the overall MAD to 20 meV. Only water inside the CNT remains challenging being overestimated by 72 meV, which is worse than all other vdW corrected semilocal DFAs (with the exception of PBE-TS).

For a better visual comparison of the performance of the DFAs for the different WaC18 subsets, we show a graphical representation of the individual rms errors separated into the different functional classes in Fig. 2. LDA is not reported as it yields very bad results. In the GGA panel, we show the results for the three most accurate GGA functionals, after inclusion of D4 for long-range vdW interactions. PBE-D4 is not included as some of the various modifications of the PBE exchange enhancement factors prove better, most notably RPBE and revPBE that are both known to give more reasonable hydrogen bond strengths.^38,148 While the water adsorption computed by RPBE-D4 and revPBE-D4 is very similar to PBE-D4, we see a clear improvement for the 2D/3D ice polymorphs. However, the MADs at 29 and 37 meV are still unexpectedly high. Especially, the denser ice structures (rhombic 2D-ice and high-pressure ice VIII) are systematically underbound. Note that the use of normal PAWs or smaller basis set expansions results in a systematic shift toward more strongly bound systems (see Table II), removing most of the bias for RPBE-D4 and revPBE-D4 and giving artificially better results.^{151,164–166} The most successful GGA tested by us is BLYP-D4 giving a very consistent performance with MADs below 20 meV for all considered benchmark sets.

Including higher derivatives (like the kinetic energy density τ) in the exchange-correlation enhancement factors give rise to the meta-GGA class. Formally, their computational cost scales with system size as the GGAs, but a stable self-consistent field convergence can be numerically more involved and typically requires larger integration grids.^20,152 TPSS is based on the PBE GGA and has a rather weak τ-dependency but still improves over PBE for many physical properties.²³ This also holds for our benchmark systems, TPSS-D4 has a rather balanced description of very accurate adsorption energies and decently good ice lattice energies, and its overall MAD of 18 meV is identical to that of BLYP-D4. SCAN and M06L are modern meta-GGAs that can cover part of the medium-range correlation and have been shown to yield good structures and energies for diversely bonded systems.^126,153 Still, both underbind all water adsorption systems, which can be partially compensated by correction schemes (see SCAN-D4). However, since SCAN already includes some vdW forces, it is nontrivial to combine SCAN with correction schemes.^20,21 This is particularly relevant for 2D/3D-ice, where plain SCAN already overbinds all systems.

The next DFA rung requires the inclusion of nonlocal (Fock) exchange resulting in hybrid functionals. The most widely used hybrid DFAs are PBE0 and B3LYP, and while they are typically only of medium quality for many chemical properties,²³ PBE0-D4 and B3LYP-D4 consistently improve over their GGA parents. In particular, the improvement for the 2D/3D ice polymorphs is significant. Of all tested methods, PBE0-D4 has the closest agreement with the reference interaction energies for all considered systems with an overall smallest MAD of 12 meV. Importantly, all tested relative stability sequences (ice polymorphs and the adsorption motifs) of PBE0-D4 are correct. The relative adsorption energies of the different binding motifs on graphene and the CNT are actually within the stochastic uncertainty of the DMC references.

The simplified DFAs (sHF-3c,^135,136 HSE-3c,¹³⁷ B97-3c;¹³⁸ see Ref. 154 for an overview) give mixed results. Overall, their accuracy is similar to the average vdW corrected DFA. Especially, HSE-3c has problems at describing the strong hydrogen bonds in 2D/3D-ice, most likely due to remaining basis set superposition errors that cannot be fully compensated. B97-3c, on the other hand, employs a slightly larger basis set expansion and gives reasonably balanced results. As those methods have been designed for increased computational speed (speedup of up to 2 orders of magnitude compared to converged basis set DFT¹⁵⁴), they might still be useful for screening applications.

We find the following points essential for a well-balanced description of both water adsorbed on nanostructures as well as within ice polymorphs:

• Correlation effects beyond the local density approximation (avoid HF and LDA).

• Use of a modern vdW correction (D3/D4, MBD, or vdW-DF2).

• Converged numerical settings with hard cores for PAWs or expansions beyond triple-ζ quality for atom-centered basis sets.

• GGA, meta-GGA, and hybrid DFAs are similarly good for adsorption.

• Fock exchange improves strong H-bonds in ice polymorphs.

As shown in Fig. 2, the best GGA, meta-GGA, vdW-DF, and simplified DFT methods (BLYP-D4, TPSS-D4, rev-vdW-DF2, and B97-3c, respectively) fulfill the first three points and perform rather similarly. Significantly more accurate are the best hybrid functionals mainly due to their improved description of 2D/3D ice. In terms of computational cost, GGAs seem to have the best accuracy vs effort ratio, while hybrids should be used when aiming for the highest accuracy. The most successful DFAs have errors consistently well below chemical accuracy (1 kcal/mol = 43 meV), challenging experimental errors of sublimation enthalpies.^39,155

In a previous effort to judge “How good is DFT for water?”⁷⁷ a scoring scheme has been devised to judge the performance of approximated methods for the properties of the water monomer, the dimer, the hexamer, and ice structures. Physical quantities scored are monomer symmetric stretch frequency $fssmono$⁠, dimer binding energy $Ebdim$⁠, ring-hexamer binding energy per monomer $Ebring$⁠, ice Ih lattice energy E_Ih, difference $ΔEbpr-ring$ of binding energies per monomer of prism and ring isomers of the hexamer, difference ΔE_Ih-VIII of lattice energies of ice Ih and VIII, equilibrium O-O distance R_OO in dimer, and equilibrium volumes per monomer $VeqIh$⁠, $VeqVIII$ of ice Ih and VIII.¹⁵⁶ For more details on the scoring system, see Ref. 77. For some of the DFAs examined in our benchmark study, we report their performance for this scoring system in Table V. As expected, LDA and HF are not reliable to describe water though some individual scores like the relative stability of the prism and ring hexamer are fortuitously good. The same holds for uncorrected as well as dispersion corrected PBE as already recognized in the original study.⁷⁷ On the other hand, BLYP-D4, TPSS-D4, optB88-vdW are performing reasonably well. Especially energetic and geometric properties are well reproduced although the lattice energy of ice Ih seems to be problematic. Symmetric stretch frequencies of the water monomer are as usual underestimated by all (meta-)GGA functionals. Here, it is worth pointing out that out of the nonhybrid functionals, the low-cost method B97-3c yields the best frequencies and overall performs competitively to the dispersion corrected DFAs. In agreement with our current benchmark, PBE0-D4 is the best performing method yielding an overall score of 82, which is indeed higher than any other DFT method tested on this set so far.

We have gathered and computed well converged reference interaction energies of water with carbon nanostructures and of water within two- and three-dimensional ice polymorphs compiled in the new WaC18 benchmark set. Combined, this gives a challenging set of large and complex systems, ideally suited to benchmark approximated methods. Importantly, those systems are larger than standard noncovalent interaction benchmark sets based on molecular dimers of small to medium sized molecules like S22¹⁴¹ and S66.¹⁵⁷ The 0D, 1D, 2D, and 3D periodicity covered here also includes new aspects compared to benchmark sets of large supramolecular complexes like L7¹⁵⁸ and S12l.¹⁵⁹ In contrast to benchmarks using back-corrected experimental references, our high-level theoretical interaction energies make the comparison much more straightforward without the complication of thermal and zero-point energy contributions. We span a broad range of interaction energies from nonbinding (water@benzene, 0-leg motif) to moderately strong binding (ice Ih with lattice energy of by −615 meV). Figure 3 shows an overview of the different binding strengths by correlating the reference energies with a few considered DFAs. The correlation highlights that we roughly follow Jacob’s ladder classification of DFAs with HF and LDA being unreliable, and PBE, PBE-D4, and PBE0-D4 step by step increasing the accuracy. Of all methods considered, BLYP-D4, TPSS-D4, rev-vdW-DF2, and PBE0-D4 are the most accurate within their respective functional class. Replacing D4 with the older D3 or the MBD vdW correction leads to minor deterioration and can be used when D4 is not available. Our present benchmark focuses on specific noncovelent interactions only; from previous studies, it is known that TPSS-D3 and PBE0-D3 yield very good equilibrium geometries²⁸ and organometallic reaction energies.¹⁶⁰ In the large GMTKN55 benchmark,²³ BLYP-D3 and TPSS-D3 are among the best performing GGAs and meta-GGAs, respectively, consistent with our findings.

We see our benchmark results as a guideline for future simulation studies of water in the condensed phase (liquid or solid) and water–carbon nanostructure interfaces. Additionally, the provided WaC18 dataset can help as a challenging cross check for other DFAs, classical force fields, machine learning potentials, tight-binding Hamiltonians, and to test other many body electronic structure methods such as the Random Phase Approximation (RPA) or Møller-Plesset (MP) perturbation theory.

See supplementary material for geometries and reference energies of all considered structures of the WaC18 benchmark set, including water@graphene, water@CNT, water@AH, 2D-ice, and 3D-ice. Geometries are provided in VASP POSCAR and XYZ format.

We thank Stefan Grimme and Eike Caldeweyher for providing access to the dftd4 code. J.G.B. acknowledges support by the Alexander von Humboldt Foundation. A.Z. and D.A.’s work is sponsored by the Air Force Office of Scientific Research, Air Force Material Command, U.S. Air Force, under Grant No. FA9550-19-1-7007. A.Z. and A.M. were supported by the European Research Council (ERC) under the European Union’s Seventh Framework Program (No. FP/2007-2013)/ERC Grant Agreement No. 616121 (HeteroIce project). This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC05-00OR22725. We are also grateful, for computational resources, to the ARCHER UK National Supercomputing Service, United Kingdom Car–Parrinello (UKCP) consortium (No. EP/F036884/1), the London Centre for Nanotechnology, University College London (UCL) Research Computing, and the UK Materials and Molecular Modelling Hub, which is partially funded by EPSRC (No. EP/P020194/1).