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

TL;DR

This study assesses the accuracy of 28 density functional approximations in modeling water interactions with carbon nanostructures, providing a benchmark set and identifying the most reliable methods for such systems.

Contribution

The paper introduces WaC18, a comprehensive benchmark set for water-carbon interactions, and evaluates the performance of various DFT approaches, highlighting the importance of van der Waals corrections.

Findings

Modern vdW corrections significantly improve DFT accuracy.

Best approaches identified: BLYP-D4, TPSS-D4, rev-vdW-DF2, PBE0-D4.

Mean absolute deviation varies from 135 meV to 12 meV.

Abstract

Due to their current and future technological applications, including realisation of water filters and desalination membranes, water adsorption on graphitic sp$^{2}$-bonded carbon is of overwhelming interest. However, these systems are notoriously challenging to model, even for electronic structure methods such as density functional theory (DFT), because of the crucial role played by London dispersion forces and non-covalent interactions in general. Recent efforts have established reference quality interactions of several carbon nanostructures interacting with water. Here, we compile a new benchmark set (dubbed \textbf{WaC18}), which includes a single water molecule interacting with a broad range of carbon structures, and various bulk (3D) and two dimensional (2D) ice polymorphs. The performance of 28 approaches, including semi-local exchange-correlation functionals, non-local (Fock) exchange contributions, and long-range van der Waals (vdW) treatments, are tested by computing the deviations from the reference interaction energies. The calculated mean absolute deviations on the WaC18 set depends crucially on the DFT approach, ranging from 135 meV for LDA to 12 meV for PBE0-D4. We find that modern vdW corrections to DFT significantly improve over their precursors. Within the 28 tested approaches, we identify the best performing within the functional classes of: generalized gradient approximated (GGA), meta-GGA, vdW-DF, and hybrid DF, which are BLYP-D4, TPSS-D4, rev-vdW-DF2, and PBE0-D4, respectively.

Tables

Table 1. Table 1: WaC18 reference interaction energies given per water molecule in meV. The reported error represents the uncertainty on the evaluation.

w@graphene Brandenburg et al. (2019)			w@benzene Brandenburg et al. (2019)			w@coronene Brandenburg et al. (2019)
0-leg	$-90\pm 6$		0-leg	$+43\pm 1$		0-leg	$-61\pm 3$
1-leg	$-92\pm 6$		1-leg 222Some earlier studies reported the 1-leg structure as most stable, which might be due to small differences in numerical and geometrical setups.Feller (1999); Slipchenko and Gordon (2009)	$-124\pm 3$		1-leg	$-118\pm 5$
2-leg	$-99\pm 6$		2-leg 222Some earlier studies reported the 1-leg structure as most stable, which might be due to small differences in numerical and geometrical setups.Feller (1999); Slipchenko and Gordon (2009)	$-136\pm 2$		2-leg	$-143\pm 4$
w@CNT 111Evaluation based on Ref. Al-Hamdani, Alfè, and Michaelides, 2017b, but recomputed in this work. See details in Section II.2.
external	$-85\pm 18$
internal	$-287\pm 16$
3D ice Zen et al. (2018)			2D ice Chen et al. (2016)
Ih	$-615\pm 5$		hex.	$-423\pm 3$		pent.	$-419\pm 3$
II	$-613\pm 6$		square	$-404\pm 3$		rhombic	$-389\pm 3$
VIII	$-594\pm 6$

Table 2. Table 2: Lattice energy convergence in meV of ice Ih using various numerical settings based on PBE and a Γ Γ \Gamma -centered 3 × \times 3 × \times 3 k 𝑘 k -grid. paw ( a , b ) The best estimates are highlighted in bold.

Vasp	PAW
PW cutoff [eV]	soft	normal	hard
500	-798.4	-665.3	-646.2
700	-798.1	-664.6	-636.8
1000	-798.0	-664.7	-637.1
Crystal		w.o.c.111Supramolecular approach without counterpoise (CP) correction.	CP-corrected
def2-mSVP		-1104.2	-858.8
def2-TZVPP		-721.2	-646.7
def2-QZVPP		-665.6	-638.2
CBS(TZ,QZ)222Basis set extrapolation using optimized exponents.Neese, Hansen, and Liakos (2009)		-657.9	-637.0

Table 3. Table 3: Interaction energies and RMS deviations (meV) from the reference data for various electronic structure methods for the WaC18 benchmark set. All values have been consistently computed in the present study. Systems are as listed in Table  1 and shown in Fig.  1 . The smallest RMS in each category is highlighted in bold. Apart from reference data taken from Refs  Brandenburg  et al. , 2019 ; Zen  et al. , 2018 ; Chen  et al. , 2016 , all other data has been computed as part of this study. References to the DFAs used are also given.

	w@graphene				w@CNT			w@benzene				w@coronene				2D-ice					3D-ice
	0-leg	1-leg	2-leg		ext.	int.		0-leg	1-leg	2-leg		0-leg	1-leg	2-leg		hex.	pent.	squ.	rhom.		Ih	II	VIII		RMS
Reference	DMC				DMC			CCSD(T)				CCSD(T)				DMC					DMC
$E_{\text{int}}$111Interaction energy at fixed equilibrium geometry provided as Supporting Information.sup	-90	-92	-99		-85	-287		43	-124	-136		-61	-118	-143		-423	-419	-404	-389		-615	-613	-594		–
$\Delta$222Uncertainty estimation of reference interaction energy.	6	6	6		18	16		1	3	2		3	5	4		3	3	3	3		5	6	6		–
Local density approximation
LDA	-96	-114	-125		-121	-235		4	-172	-197		-62	-135	-155		-710	-699	-657	-637		-1016	-978	-876		193
PBE with various vdW corrections
PBEPerdew, Burke, and Ernzerhof (1996)	-9	-26	-19		-26	-82		86	-89	-81		21	-44	-50		-434	-416	-370	-354		-639	-571	-462		79
PBE-VV10Perdew, Burke, and Ernzerhof (1996); Vydrov and Van Voorhis (2010)	-123	-131	-139		-121	-375		22	-140	-149		-78	-138	-154		-498	-490	-455	-437		-755	-733	-672		63
PBE-dDsCPerdew, Burke, and Ernzerhof (1996); Steinmann and Corminboeuf (2010)	-109	-132	-143		-123	-390		26	-142	-155		-70	-142	-160		-482	-476	-447	-427		-739	-738	-688		61
PBE-TSPerdew, Burke, and Ernzerhof (1996); Tkatchenko and Scheffler (2009)	-116	-141	-160		-136	-408		28	-143	-162		-71	-145	-175		-467	-462	-429	-410		-712	-698	-621		52
PBE-MBDPerdew, Burke, and Ernzerhof (1996); Tkatchenko et al. (2012)	-93	-120	-126		-112	-310		41	-137	-146		-51	-128	-145		-478	-472	-437	-420		-721	-694	-619		41
PBE-D2Perdew, Burke, and Ernzerhof (1996); Grimme (2006)	-89	-128	-135		-119	-303		35	-147	-161		-53	-140	-154		-489	-482	-447	-432		-731	-698	-637		47
PBE-D3Perdew, Burke, and Ernzerhof (1996); Grimme et al. (2010)	-85	-117	-124		-112	-291		39	-139	-147		-49	-127	-144		-476	-467	-430	-412		-716	-679	-597		36
PBE-D4Perdew, Burke, and Ernzerhof (1996); Caldeweyher et al. (2019)	-104	-115	-118		-107	-314		30	-132	-138		-65	-122	-137		-474	-464	-426	-409		-711	-677	-593		35
GGAs with D4Caldeweyher et al. (2019) vdW correction
RPBE-D4Hammer, Hansen, and Norskov (1999)	-102	-110	-108		-100	-322		34	-115	-117		-62	-114	-124		-412	-404	-371	-355		-632	-594	-519		26
revPBE-D4Zhang and Yang (1998)	-97	-105	-105		-94	-308		43	-109	-113		-56	-110	-121		-402	-392	-357	-340		-620	-585	-515		29
BLYP-D4Becke (1988); Lee, Yang, and Parr (1988)	-109	-117	-118		-104	-332		31	-114	-119		-74	-116	-128		-439	-432	-403	-386		-659	-645	-574		22
meta-GGA (with D4Caldeweyher et al. (2019) vdW correction)
M06LZhao and Truhlar (2006)	-55	-64	-67		-58	-383		53	-93	-111		-29	-76	-95		-338	-339	-343	-321		-516	-545	-577		56
SCANSun, Ruzsinszky, and Perdew (2015)	-63	-74	-84		-78	-197		45	-123	-144		-29	-92	-116		-464	-459	-439	-421		-667	-655	-615		35
SCAN-D4Sun, Ruzsinszky, and Perdew (2015); Brandenburg et al. (2016)	-106	-116	-129		-113	-304		31	-136	-158		-62	-122	-150		-476	-473	-455	-436		-766	-694	-661		52
TPSS-D4Tao et al. (2003)	-103	-110	-113		-99	-324		40	-119	-125		-62	-114	-131		-446	-433	-387	-370		-664	-638	-549		22
1st and 2nd generation vdW-DFs
vdW-DF1Dion et al. (2004a)	-136	-134	-133		-107	-455		26	-117	-109		-99	-136	-147		-346	-346	-326	-309		-564	-557	-522		63
optB86b-vdWKlimeš, Bowler, and Michaelides (2010)	-142	-144	-150		-121	-454		23	-132	-137		-98	-150	-167		-448	-444	-413	-395		-708	-704	-661		59
vdW-DF2Lee et al. (2010)	-120	-123	-128		-110	-401		15	-124	-125		-92	-128	-140		-404	-406	-397	-378		-624	-624	-598		33
rev-vdW-DF2Hamada (2014)	-105	-110	-115		-95	-360		38	-120	-122		-67	-122	-137		-446	-439	-406	-388		-685	-666	-610		29
Hybrid functionals (with D4Caldeweyher et al. (2019) vdW correction)
HF	33	38	46		42	9		166	-41	-29		75	29	-2		-227	-226	-216	-201		-298	-271	-292		198
HF-D4	-109	-95	-106		-88	-315		68	-111	-133		-56	-93	-137		-332	-347	-356	-339		-468	-491	-568		57
revPBE0-D4Zhang and Yang (1998); Adamo and Barone (1999)	-99	-101	-105		-91	-301		51	-115	-126		-52	-106	-130		-389	-383	-353	-337		-583	-560	-516		32
B3LYP-D4Becke (1993); Stephens et al. (1994)	-111	-115	-119		-103	-328		36	-122	-132		-71	-115	-136		-442	-438	-413	-397		-641	-639	-588		18
PBE0-D4Adamo and Barone (1999)	-103	-108	-114		-100	-305		44	-131	-142		-57	-115	-140		-449	-443	-411	-395		-643	-623	-582		14
Simplified density functional approximations
sHF-3cSure and Grimme (2013); Cutini et al. (2016)	-63	-90	-108		-102	-242		65	-119	-151		-7	-123	-166		-467	-451	-425	-412		-692	-670	-563		34
HSE-3cBrandenburg, Caldeweyher, and Grimme (2016)	-114	-97	-123		-130	-260		70	-158	-194		-43	-129	-185		-511	-495	-454	-419		-735	-709	-627		54
B97-3cBrandenburg et al. (2018)	-112	-133	-137		-131	-321		40	-145	-155		-74	-138	-166		-459	-446	-412	-388		-689	-656	-590		32

Table 4. Table 4: Mean deviation (MD) and mean absolute deviation (MAD) in meV of the interaction energies from the references, see Table  1 and Fig.  1 .

	w@nano111Combination of the three adsorption sets water@graphene, water@CNT, and water@AH.			2D/3D-ice			all
Method	MD111Combination of the three adsorption sets water@graphene, water@CNT, and water@AH.	MAD		MD	MAD		MD	MAD
Local density approximation
LDA	-20	29		-302	302		-130	135
PBE with various vdW corrections
PBE	79	79		30	40		60	64
PBE-VV10	-30	30		-83	83		-51	51
PBE-dDsC	-32	32		-77	77		-49	49
PBE-TS	-40	40		-49	49		-43	43
PBE-MBD	-12	14		-55	55		-29	30
PBE-D2	-18	20		-65	65		-37	38
PBE-D3	-10	13		-46	46		-24	26
PBE-D4	-12	13		-42	43		-24	25
GGAs with D4 vdW correction
RPBE-D4	-4	14		24	29		7	20
revPBE-D4	2	12		35	37		15	21
BLYP-D4	-10	18		-12	18		-10	18
meta-GGA (with D4 vdW correction)
M06L	19	37		68	68		38	49
SCAN	22	23		-38	38		-1	29
SCAN-D4	-16	16		-72	72		-38	38
TPSS-D4	-6	12		-5	27		-6	18
1st and 2nd generation vdW-DFs
vdW-DF1	-32	39		70	70		7	51
optB86b-vdW	-44	44		-45	45		-44	44
vdW-DF2	-26	28		4	11		-14	21
rev-vdW-DF2	-11	16		-26	26		-17	20
Hybrids (with D4 vdW correction)
HF	142	142		247	247		182	182
HF-D4	1	12		79	79		32	39
revPBE0-D4	2	10		48	48		20	25
B3LYP-D4	-11	14		-14	16		-12	15
PBE0-D4	-7	9		-13	16		-9	12
Simplified density functional approximations
sHF-3c	8	20		-32	40		-8	28
HSE-3c	-16	29		-70	70		-37	45
B97-3c	-25	25		-26	28		-26	26

Table 5. Table 5: Percentage scores of selected DFAs using the scoring scheme from Ref.  Gillan, Alfè, and Michaelides, 2016 . Physical quantities scored are monomer symmetric stretch frequency f s ​ s mono superscript subscript 𝑓 𝑠 𝑠 mono f_{ss}^{\text{mono}} , dimer binding energy E b dim superscript subscript 𝐸 𝑏 dim E_{b}^{\text{dim}} , ring-hexamer binding energy per monomer E b ring superscript subscript 𝐸 𝑏 ring E_{b}^{\text{ring}} , ice Ih lattice energy E Ih subscript 𝐸 Ih E_{\text{Ih}} , difference Δ ​ E b pr-ring Δ superscript subscript 𝐸 𝑏 pr-ring \Delta E_{b}^{\text{pr-ring}} of binding energies per monomer of prism and ring isomers of the hexamer, difference Δ ​ E Ih-VIII Δ subscript 𝐸 Ih-VIII \Delta E_{\text{Ih-VIII}} of lattice energies of ice Ih and VIII, equilibrium O-O distance R OO subscript 𝑅 OO R_{\text{OO}} in dimer, and equilibrium volumes per monomer V eq Ih superscript subscript 𝑉 eq Ih V_{\text{eq}}^{\text{Ih}} , V eq VIII superscript subscript 𝑉 eq VIII V_{\text{eq}}^{\text{VIII}} of ice Ih and VIII. The total score in the final column is the average of the individual scores, unmarked values have been computed in this study.

	$f_{ss}^{\text{mono}}$	$E_b^{\text{dim}}$	$E_b^{\text{ring}}$	$E_{\text{Ih}}$	$\Delta E_b^{\text{pr-ring}}$	$\Delta E_{\text{Ih-VIII}}$	$R_{\text{OO}}$	$V_{\text{eq}}^{\text{Ih}}$	$V_{\text{eq}}^{\text{VIII}}$	Total
Reference	111References taken as gathered in Ref. Gillan, Alfè, and Michaelides, 2016.3812 cm${}^{\text{-1}}$	111References taken as gathered in Ref. Gillan, Alfè, and Michaelides, 2016.-216 meV	111References taken as gathered in Ref. Gillan, Alfè, and Michaelides, 2016.-319 meV	222Negative MD means too strongly bounded system-615 meV	111References taken as gathered in Ref. Gillan, Alfè, and Michaelides, 2016.13 meV	222For consistency within this article, references taken from Table 1.21 meV	111References taken as gathered in Ref. Gillan, Alfè, and Michaelides, 2016.2.91 Å	333Values taken from Ref. Brandenburg, Maas, and Grimme, 2015 to consistently exclude zero-point and thermal effects.30.9 Å3	333Values taken from Ref. Brandenburg, Maas, and Grimme, 2015 to consistently exclude zero-point and thermal effects.19.1 Å3	100
Tolerance	20 cm${}^{\text{-1}}$	10 meV	10 meV	10 meV	10 meV	10 meV	0.01 Å	1%	1%
LDA	60	0	0	0	100	0	0	444Values not computed as expected to be unreliable.–	444Values not computed as expected to be unreliable.–	23
HF	0	30	0	0	80	90	0	444Values not computed as expected to be unreliable.–	444Values not computed as expected to be unreliable.–	29
PBE	50	100	90	80	0	0	90	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.100	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.20	59
PBE-D4	50	90	60	0	100	10	80	666Values taken from Ref. Brandenburg, Maas, and Grimme, 2015 using the D3 dispersion correction.50	666Values taken from Ref. Brandenburg, Maas, and Grimme, 2015 using the D3 dispersion correction.100	60
BLYP-D4	30	100	90	60	100	40	100	666Values taken from Ref. Brandenburg, Maas, and Grimme, 2015 using the D3 dispersion correction.90	666Values taken from Ref. Brandenburg, Maas, and Grimme, 2015 using the D3 dispersion correction.80	77
TPSS-D4	50	100	80	50	90	10	80	666Values taken from Ref. Brandenburg, Maas, and Grimme, 2015 using the D3 dispersion correction.60	666Values taken from Ref. Brandenburg, Maas, and Grimme, 2015 using the D3 dispersion correction.70	66
optB88-vdW	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.60	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.100	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.90	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.20	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.100	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.100	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.50	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.80	555Values taken from Ref. Gillan, Alfè, and Michaelides, 2016.100	78
PBE0-D4	80	100	80	80	90	70	100	777Values taken from Ref. Gillan, Alfè, and Michaelides, 2016 using the TS dispersion correction.70	777Values taken from Ref. Gillan, Alfè, and Michaelides, 2016 using the TS dispersion correction.70	82
B97-3c	70	90	70	30	100	30	80	50	70	66

Equations

$$E_{\text{int}}=E_{\text{w@nano}}-(E_{\text{w}}+E_{\text{nano}}),$$

$$E_{\text{int}}=E_{\text{solid}}/N_{\text{mol}}-E_{\text{w}}.$$

Full text

Keywords: Density functional theory, water adsorption, carbon nanomaterials, graphene, ice, benchmark, van der Waals

I Introduction

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.Fumagalli et al. (2018); Abraham et al. (2017); Joshi et al. (2014); Nair et al. (2012); Algara-Siller et al. (2015); Secchi et al. (2016); Radha et al. (2016); Hummer, Rasaiah, and Noworyta (2001); Strogatz et al. (2005); Holt (2006); Secchi et al. (2016) Topologically similar materials can have substantially different properties Tocci, Joly, and Michaelides (2014); Michaelides (2016); Esfandiar et al. (2017); Siria, Bocquet, and Bocquet (2017); Yoshida et al. (2018) 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.Burke (2012); Becke (2014); Yu, Li, and Truhlar (2016); Maurer et al. (2019) Modern density functional approximations (DFAs) combine semi-local 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,Brandenburg et al. (2016); Hermann and Tkatchenko (2018) 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, molecular crystals).Maurer et al. (2019); Al-Hamdani et al. (2017) 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, so 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 non-covalent interactions.Goerigk et al. (2017); Mardirossian and Head-Gordon (2016a); Mardirossian et al. (2017); Mardirossian and Head-Gordon (2017); Wang et al. (2017) Recently, some focus of DFT benchmarking moved to the description of equilibrium geometries.Grimme et al. (2015); Witte et al. (2015); Piccardo et al. (2015); Grimme and Steinmetz (2013) 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,Staroverov et al. (2004); Pernot et al. (2015); Peintinger, Oliveira, and Bredow (2013); Zhang et al. (2018) similar data for more weakly bound molecular crystals have substantially higher uncertainty.Otero-de-la-Roza and Johnson (2012); Reilly and Tkatchenko (2013); Brandenburg, Maas, and Grimme (2015) This is due to the experimental measurement uncertainty,Chickos (2003) 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.Björneholm et al. (2016)

Concerning theoretical reference calculations, exciting progress has been made in the field of high-level wavefunction methods. On the one hand, embedding techniques Masur et al. (2016); Bygrave, Allan, and Manby (2012); Beran et al. (2014) and local approaches of coupled cluster theories Riplinger et al. (2013, 2016); Maurer et al. (2013); Schütz and Werner (2001); Nagy, Samu, and Kállay (2018); Paulus (2006) have made the gold standard of quantum chemistry applicable to molecular systems with a few hundred atoms and molecular crystals of small molecules.Beran (2016); Yang et al. (2014); Gruber et al. (2018) On the other hand, new algorithmic developments in the field of diffusion Monte Carlo (DMC, a quantum Monte Carlo technique) have made the computation of chemically accurate lattice energies of small molecular crystals feasible within reasonable computational effort.Zen et al. (2016, 2018) These developments together with an increased capacity of available computational resources makes the interaction energy determination of extended systems feasible. This is an important step towards the better understanding of large non-covalently interacting systems as recently highlighted.Michaelides et al. (2015); Al-Hamdani and Tkatchenko (2019a); Maurer et al. (2019)

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,Brandenburg et al. (2019) water on benzene and coronene as two representative aromatic hydrocarbons (abbreviated as AH),Brandenburg et al. (2019) and water on a carbon nanotube (CNT).Al-Hamdani, Alfè, and Michaelides (2017a) 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) iceChen et al. (2016) as well as bulk (3D) ice polymorphs.Zen et al. (2018) Therefore, we obtain a dataset of eighteen 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,Ma et al. (2011); Ajala et al. (2019); Chen et al. (2016); Al-Hamdani, Alfè, and Michaelides (2017a); Santra et al. (2013); He et al. (2012); Hirata et al. (2017); Hamada and Otani (2010); Voloshina et al. (2011); Lei et al. (2016); Hamada (2012); Silvestrelli and Ambrosetti (2014); Ambrosetti and Silvestrelli (2011); Lorenz et al. (2014); Leenaerts, Partoens, and Peeters (2009); Heßelmann (2013); Rubeš et al. (2009); Jenness, Karalti, and Jordan (2010); Jordan and Heßelmann (2019) but they typically include a limited set of DFAs, the recent vdW developments are not included, and the used reference data is 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. Gillan, Alfè, and Michaelides, 2016 that focused on DFT recommendations for water.

In Section 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 (III). Conclusions and a future perspective are given in section IV.

II Benchmark Setup

II.1 Systems under consideration

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. Chen et al., 2016 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 supporting information files to make our compiled benchmark easily accessible to other researchers.sup

II.2 Reference interaction energies

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

[TABLE]

where $E_{\text{w@nano}}$ is in the bound configuration and the geometries of the individual fragments $E_{\text{w}}$ and $E_{\text{nano}}$ for the unbound configuration are kept frozen.geo The interaction energies for 2D and 3D ice are given per molecule as

[TABLE]

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 1. The reference values for the 2-dimensional ice crystals are taken from Ref. Chen et al., 2016, while for the three bulk-ice crystals are given in Ref. Zen et al., 2018. The reference values for water on benzene, coronene, and graphene with different orientations are taken from Ref. Brandenburg et al., 2019. The reference interaction energy between water and CNT where previously investigated in Ref. Al-Hamdani, Alfè, and Michaelides, 2017b, 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 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 non-covalent interactions, provided that care is taken to ensure that an accurate setup is used for each method.Zen et al. (2016); Dubecký, Mitas, and Jurecka (2016); Al-Hamdani and Tkatchenko (2019b) Agreement between methods is not expected beyond a given precision. To this aim, in Table 1 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.Brandenburg et al. (2019)

Although reported DMC results are coming from different studies, the setup is consistent among them. DMC simulations were carried out with the casino code.Needs et al. (2010) Dirac-Fock pseudopotentials Trail and Needs (2005a, b) with the localization approximation Mitas, Shirley, and Ceperley (1991) (LA) are used. The trial wavefunctions were of the Slater-Jastrow type with single Slater determinants and the single particle orbitals obtained from DFT-LDA plane-wave calculations performed with pwscf Giannozzi et al. (2009) and re-expanded in terms of B-splines.Alfè and Gillan (2004) The Jastrow factor included electron-electron, electron-nucleus and electron-electron-nucleus terms. The parameters of the Jastrow were carefully optimized by minimizing the variance, within a variational QMC scheme. The recently developed size-consistent DMC algorithm (ZSGMA)Zen et al. (2016) was used. Finite time-step errors are carefully minimized by performing simulations with different values of the time-step, untill the bias appears safely smaller than the stochastic error. In periodic calculations, finite-size corrections are applied either using the model periodic Coulomb interactionFraser et al. (1996); Williamson et al. (1997); Kent et al. (1999) or the Kwee, Zhang, and Krakauer (2008) approach.

We reevaluate here the binding energy of water on the CNT because one specific aspect of the setup used in Ref. Al-Hamdani, Alfè, and Michaelides, 2017b 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 elsewhereZen et al. (2019)) 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 pseudopotentialsTrail and Needs (2017) and the ZSGMA algorithm.Zen et al. (2016) The reported binding values are obtained with time step $\tau=0.01$ a.u., which gives errors less than the stochastic uncertainty (standard deviation of 18 meV).

II.3 DFT Computational Details

DFT test calculations based on the PBE functional Perdew, Burke, and Ernzerhof (1996) 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 Neese (2012) with large aug-cc-pVQZ, aug-cc-pV5Z Dunning (1989); Kendall, Dunning, and Harrison (1992), and def2-QZVPPD Weigend, Furche, and Ahlrichs (2003) basis sets, Crystal17 Erba et al. (2017); Dovesi et al. (2018) with a def2-QZVPPD(-f) basis; and Vasp 5.4 Kresse and Furthmüller (1996a, b) with projector-augmented plane waves (hard PAWs Blöchl (1994); Kresse and Joubert (1999)) 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$\times$5 supercell with 53 atoms in the unit cell) and the CNT (non-metallic (10,0) configuration with 83 atoms in the unit cell), this reduces to a $\Gamma$-point calculation. For all PAW calculations the non-periodic 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 non-periodic directions consistent with the reference calculations.

For the water@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@graphene, water@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 2 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 neither converged employing 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 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 semi-local DFAs are corrected for missing long-range correlation effects by means of a variety of different semi-classical and nonlocal density based corrections (D2 Grimme (2006), D3 Grimme et al. (2010); Grimme, Ehrlich, and Goerigk (2011), D4 Caldeweyher, Bannwarth, and Grimme (2017); Caldeweyher et al. (2019), TS Tkatchenko and Scheffler (2009), MBD Tkatchenko et al. (2012), VV10 Vydrov and Van Voorhis (2010), dDsC Steinmann and Corminboeuf (2010), vdW-DF Dion et al. (2004a)). TS and MBD are used with the non-iterative 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. Klimeš and Michaelides, 2012; Grimme et al., 2016; Hermann, DiStasio, and Tkatchenko, 2017; Berland et al., 2015 for further overview on vdW interactions in the DFT framework.

III Results and Discussion

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 3 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 4.

III.1 Impact of vdW interactions

From looking at Tables 3 and 4, it is clear that Hartree-Fock (HF) misses all (Coulomb) correlation effects and cannot describe any of these non-covalently bound systems appropriately.London (1937); Stone (1997) 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. Brandenburg et al., 2019 as well as many studies on molecular dimers.Jurečka et al. (2006) That the local density approximation (LDA) yields an inconsistent description of small vdW complexes has been known since the mid-90s.Kristyán and Pulay (1994) 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 to 0.2 Å (see Ref. Al-Hamdani, Alfè, and Michaelides, 2017a; Al-Hamdani et al., 2017; Ma et al., 2011), the LDA minimum is at substantially smaller distance and cannot be recommended for either geometry or stability estimates of vdW bound systems.

Including semi-local 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. Klimeš and Michaelides, 2012; Grimme et al., 2016; Berland et al., 2015). We test a broad range of these vdW corrections coupled with PBE (see Table 3, 4). The error spread is still substantial and in particular the older effective pair-wise schemes (PBE-VV10, PBE-dDsC, PBE-TS, PBE-D2) do not perform well. Recent vdW developments pay off and we can see a clear improvement of PBE-MBD Tkatchenko et al. (2012) and PBE-D4 Caldeweyher et al. (2019) 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 levelAxilrod and Teller (1943); Muto (1944) as included in the D3 method.Grimme et al. (2010) 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.Peach, Williamson, and Tozer (2011); Johnson, Otero-de-la Roza, and Dale (2013); Bryantsev et al. (2009); Goerigk, Kruse, and Grimme (2011); Řezáč et al. (2015) 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.Gillan, Alfè, and Michaelides (2016)

DFAs with a nonlocal kernel to describe vdW interactions have been pioneered by Dion et al. (vdW-DF1).Dion et al. (2004b, 2005) 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 semi-local 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. Al-Hamdani, Alfè, and Michaelides, 2017b. Consistent with previous studies,Lee et al. (2010); Hamada (2014) 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 semi-local DFAs (with the exception of PBE-TS).

III.2 Performance by Jacob’s ladder classification

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.Goerigk, Kruse, and Grimme (2011); Brandenburg, Maas, and Grimme (2015) 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 towards more strongly bound systems (see Table 2), removing most of the bias for RPBE-D4 and revPBE-D4 and giving artificially better results.ice 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 $\tau$) 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.Mardirossian and Head-Gordon (2016b); Brandenburg et al. (2016) TPSS is based on the PBE GGA and has a rather weak $\tau$-dependency, but still improves over PBE for many physical properties.Goerigk et al. (2017) 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, its overall MAD of 18 meV is identical 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.Sun et al. (2016); Zhao and Truhlar (2006) 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 non-trivial to combine SCAN with correction schemes.Brandenburg et al. (2016); Hermann and Tkatchenko (2018) This is particularly relevant for 2D/3D-ice, where plain SCAN already overbinds all systems.

The next DFA rung requires the inclusion of non-local (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,Goerigk et al. (2017) 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 Sure and Grimme (2013); Cutini et al. (2016), HSE-3c Brandenburg, Caldeweyher, and Grimme (2016), B97-3c Brandenburg et al. (2018), see Ref. Caldeweyher and Brandenburg, 2018 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 DFTCaldeweyher and Brandenburg (2018)), they might still be useful for screening applications.

III.3 Essential ingredients for well-balanced DFA

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-$\zeta$ 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), fulfil 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 successfull DFAs have errors consistently well below chemical accuracy (1 kcal/mol = 43 meV), challenging experimental errors of sublimation enthalpies.Whalley (1984); Chickos (2003)

III.4 Comparison to water scoring scheme

In a previous effort to judge *”How good is DFT for water?”*Gillan, Alfè, and Michaelides (2016) 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 $f_{ss}^{\text{mono}}$, dimer binding energy $E_{b}^{\text{dim}}$, ring-hexamer binding energy per monomer $E_{b}^{\text{ring}}$, ice Ih lattice energy $E_{\text{Ih}}$, difference $\Delta E_{b}^{\text{pr-ring}}$ of binding energies per monomer of prism and ring isomers of the hexamer, difference $\Delta E_{\text{Ih-VIII}}$ of lattice energies of ice Ih and VIII, equilibrium O-O distance $R_{\text{OO}}$ in dimer, and equilibrium volumes per monomer $V_{\text{eq}}^{\text{Ih}}$, $V_{\text{eq}}^{\text{VIII}}$ of ice Ih and VIII.gil For more details on the scoring system see Ref. Gillan, Alfè, and Michaelides, 2016. For some of the DFAs examined in our benchmark study, we report their performance for this scoring system in Table 5. 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.Gillan, Alfè, and Michaelides (2016) On the other hand, BLYP-D4, TPSS-D4, optB88-vdW are performing reasonably well. Especially energetic and geometric properties are well reproduced, though 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 non-hybrid 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.

IV Conclusions

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 S22Jurečka et al. (2006) and S66.Řezáč, Riley, and Hobza (2011) The 0D, 1D, 2D, and 3D periodicity covered here, also includes new aspects compared to benchmark sets of large supramolecular complexes like L7Sedlak et al. (2013) and S12l.Grimme (2012) In contrast to benchmarks using back-corrected experimental references, our high-level theoretical interaction energies make the comparison much more straight-forward without the complication of thermal and zero-point energy contributions. We span a broad range of interaction energies from non-binding (water@benzene, 0-leg motif) to moderately strong binding (ice Ih with lattice energy of by $-615$ meV). Fig. 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 the 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 geometriesGrimme et al. (2015) and organo-metallic reaction energies.Dohm et al. (2018) In the large GMTKN55 benchmarkGoerigk et al. (2017) 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 other DFAs, classical force fields, machine learning potentials, tight-binding Hamiltonians, and to test other many body electronic structure methods like Random Phase Approximation (RPA) or Møller-Plesset perturbation theory (MP).

Acknowledgements.

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, US Air Force, under Grant 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 (FP/2007-2013)/ERC Grant Agreement 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 DE-AC05-00OR22725. We are also grateful, for computational resources, to ARCHER UK National Supercomputing Service, United Kingdom Car–Parrinello (UKCP) consortium (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 (EP/P020194/1).

J.G.B and A.Z. contributed equally to this work.