Exploring water adsorption on isoelectronically doped graphene using alchemical derivatives

TL;DR

This paper demonstrates that alchemical derivatives can efficiently predict water adsorption energies on doped graphene with high accuracy, significantly aiding the design of novel 2D materials.

Contribution

It introduces a method using alchemical derivatives for high-throughput screening of water adsorption on doped graphene, validated against density functional theory.

Findings

Alchemical derivatives predict adsorption energies with <0.1 eV mean absolute error.

Filtering predictions based on molecular orbital analysis reduces error to ~0.02 eV.

The method enables reliable and efficient screening of doped graphene surfaces.

Abstract

The design and production of novel 2-dimensional materials has seen great progress in the last decade, prompting further exploration of the chemistry of such materials. Doping and hydrogenating graphene is an experimentally realised method of changing its surface chemistry, but there is still a great deal to be understood on how doping impacts on the adsorption of molecules. Developing this understanding is key to unlocking the potential applications of these materials. High throughput screening methods can provide particularly effective ways to explore vast chemical compositions of materials. Here, alchemical derivatives are used as a method to screen the dissociative adsorption energy of water molecules on various BN doped topologies of hydrogenated graphene. The predictions from alchemical derivatives are assessed by comparison to density functional theory. This screening method is found to predict dissociative adsorption energies that span a range of more than 2 eV, with a mean absolute error $<0.1$ eV. In addition, we show that the quality of such predictions can be readily assessed by examination of the Kohn-Sham highest occupied molecular orbital in the initial states. In this way, the root mean square error in the dissociative adsorption energies of water is reduced by almost an order of magnitude (down to $\sim0.02$ eV) after filtering out poor predictions. The findings point the way towards a reliable use of first order alchemical derivatives for efficient screening procedures.

Tables

Table 1. Table 1: Statistical analysis of data from four types of alchemical routes and resolved for mesomerically active (MA) and passive (MP) classification using threshold HOMO charges for each route. The threshold charges used for BN pair 1, BN pair 2, B2C and N2C are 0.203, 0.192, 0.140 and 0.025 e/atom, respectively. R 2 coefficient, Spearman’s rank correlation coefficient (r s s {}_{\text{s}} ), mean absolute error (MAE) in eV, and root mean square error (RMSE) in eV are listed. Numbers in parentheses specify the number of paths.

	R2	r${}_{\text{s}}$	MAE (eV)	RMSE (eV)
BN pair 1
total	0.14		0.068	0.193
MA (38)	0.10	0.46	0.296	0.504
MP (56)	0.72	0.85	0.032	0.048
BN pair 2
total	0.17		0.063	0.191
MA (42)	0.12	0.48	0.280	0.485
MP (52)	0.79	0.89	0.025	0.041
B2C
total	0.27		0.068	0.193
MA (20)	0.03	0.33	0.318	0.505
MP (74)	0.85	0.88	0.045	0.130
N2C
total	0.49		0.091	0.176
MA (20)	0.49	0.89	0.168	0.298
MP (74)	0.87	0.92	0.055	0.061

Equations

$$E(\lambda =1)=E^0+{\partial }_{\lambda }{E}^0\Delta \lambda +\frac{1}{2}{\partial }_{\lambda }^2{E}^0\Delta {\lambda }^2+\dots ,$$

$${\partial }_{\lambda }E$$

$${\partial }_{Z_I}{\mu }_I=\int d\mathbf{r}\frac{\mathrm{erf}[\sigma |{\mathbf{R}}_I-\mathbf{r}|]}{|{\mathbf{R}}_I-\mathbf{r}|}{\partial }_{Z_J}\rho (\mathbf{r}),$$

$$E_{ads}=E_{ads/sub}^{tot}-E_{sub}^{tot}-E_{ads}^{tot}$$

Full text

I Introduction

Recognising the enormous number of ways in which elements can be combined is both exciting and daunting in the search for more efficient, more sustainable, and safer materials for medical, engineering, and catalytic applications. High throughput screening in computational chemistry, otherwise known as virtual screening, is paving the way for materials discovery across academic and industrial research. There are various ways to screen through materials (see e.g. Refs. 1; 2; 3; 4; 5; 6). One particularly noteworthy example in catalysis was the study of Greeley et al. which involved the computational screening of 700 binary surface alloys to find a material with high activity for H2 evolution 7. The computational screening lead to the discovery and subsequent synthesis of BiPt which showed comparable activity to pure Pt experimentally.

We focus on an area of widespread interest that is, dissociative molecular adsorption on 2-dimensional substrates. In particular, graphene and hexagonal boron nitride (h-BN) are nearly isostructural materials with emerging applications in industry, including in catalysis8; 9; 10; 11; 12; 13; 14; 15. However, an important challenge in using graphene for catalysis, is overcoming its inertness. There are a number of ways to facilitate reactions at the surface of graphene such as using metal substrates14; 15; 16; 17; 18; 19; 20 to electronically dope graphene and in-plane doping of graphene with other elements21; 12; 13; 22. For instance, pristine graphene has been shown to be inert to the dissociative adsorption of water whereas, BN doped and hydrogenated graphene is far more likely to dissociate water21. Hydrogenating graphene breaks the large delocalized $\pi$ network of electrons in graphene, which is key to its inertness23; 24; 25. The hydrogenation of graphene has been extensively studied in experiments, with a number of methods of production (see e.g. Refs. 26; 27; 28; 29). In addition, doping graphene isoelectronically with BN atoms further facilitates the adsorption of molecules by forming stronger covalent bonds with adsorbates21; 17. The in-plane BN doping of graphene has also been realised experimentally in recent years30; 31; 32; 33, with increasing control over the doping process such that, nanometre-scale domains can be produced32; 30, as well as separated B and N atoms in the graphene surface33. Facilitating adsorption processes in such ways is vital for these materials to become energy efficient and applicable on a large scale. Here, we investigate how isoelectronically doping with BN away from the adsorption site, affects the dissociative adsorption energy of water on graphene.

Considering that computational molecular adsorption studies on graphene typically involve unit cells containing 30-50 carbon atoms, there are hundreds of ways to arrange a pair of boron and nitrogen atoms in such a unit cell (after accounting for redundancies by symmetry). However, the isoelectronic nature of doping in this study, and the proximity of boron, nitrogen, and carbon in the periodic table, can be utilized for efficient approximate screening schemes. Specifically, we can look to alchemical derivatives in density functional theory (DFT)34; 35; 36; 37. This method relies on exploiting the information encoded in the averaged electrostatic potential at each atom, which is analogous to the first order alchemical derivative, readily available after any self-consistent field (SCF) calculation. This and similar conceptual DFT has been discussed comprehensively in some contributions34; 38; 39, and later in Section II we give a brief introduction to the method employed. Note that alchemical derivatives have been used previously to predict various properties such as, intermolecular energies40, HOMO eigenvalues41, reaction energies42, doping in benzene43; 44, covalent bonds45, and binding in alkali halide crystals, 46 or transition metals 47; 48; 49.

In this study, the first order alchemical derivative is used to predict the dissociative adsorption energy of water on BN doped graphene, with doping occuring at different sites in the substrate. The predicted energies are compared to explicitly calculated energies to reveal the quality of predictions and to identify any outliers. Further, it is shown that outliers can be identified without additional calculations by simply using $\rho_{HOMO}$ of the initial state. The study begins with a description of the methods and the system setup in Section II, followed by the results of alchemical predictions in Section III. After identifying the main trends, further questions about the procedure and implications for water adsorption are discussed in Section IV before concluding in Section V.

II Methods

Let us begin with a brief background, followed later by details of the system setup and calculations. Firstly, any point in chemical compound space can be referred to as a discrete chemical thermodynamic micro-state. Within DFT, such a state is defined by the charge density, which results from solving an equivalent of Schrödinger’s equation for a given proton distribution $Z(\mathbf{r})$ and number of electrons $N_{e}$. As such, $Z(\mathbf{r})$ and $N_{e}$ can also be seen as extensive particle variables in a molecular grand-canonical ensemble34. The mutation of a chemical thermodynamic system into another can be achieved by thermodynamic integration with respect to a switching parameter $\lambda$. The parameter $\lambda$ simply tracks the change from the initial state to the final state. A converged integration would require sampling intermediate $\lambda$ and hence, several DFT calculations. Instead here, this mutation is approximated, using a Taylor expansion around the initial system and $\lambda$,

[TABLE]

where $\lambda=0$ corresponds to the initial system, $\lambda=1$ corresponds to the target system and hence, $\Delta\lambda=1$. Indeed it is not given that the first order term in Eq. 1 is always predictive. However, it has been observed that for relative energies, such as the adsorption energy for instance, higher order terms can cancel out resulting in useful predictions of properties40; 41; 42; 43; 44; 45; 46; 47; 48; 49. Importantly, as we see below, the first order term in Eq. 1 can be evaluated from a single DFT calculation of the initial state. In general, the first order term ($\partial_{\lambda}E^{0}$) includes the variance of the energy with changes in the proton density, the nuclear positions, and the number of electrons. However, here we consider the isoelectronic doping of a graphene sheet with fixed atomic positions, and later this is shown to be a good approximation in the system considered here. Terms involving changes in atomic positions {$\mathbf{R}_{I}$} and $N_{e}$ can therefore be neglected leaving us with the electronic contribution,

[TABLE]

where the variation of the energy with respect to a small change in nuclear charge ($Z_{I}$), damped by the error-function because of the lack of intranuclear repulsion, is known as the alchemical potential $\mu_{I}$50. This is referred to as the alchemical potential, rather than the electrostatic potential, since it quantifies the first-order energy change as a result of an “alchemical” infinitesimal variation in proton number at an atomic site. When deviating from the transmutating atom’s position, the alchemical potential becomes very similar to the electrostatic potential, $\bar{V}_{ESP}(\mathbf{r})$. For practical reasons we note that the average electrostatic potential at each atom (including the nuclear contributions omitted in Eq. II) - or alchemical potential - is readily available at the end of the SCF cycle in the widely used Vienna Ab-Initio Simulation Package (VASP)51; 52; 53; 54. Hence, we can easily evaluate the first order alchemical perturbation based approximation of the energy of any doped system from the information (i.e. $\bar{V}_{ESP}(\mathbf{r})$) provided in a single DFT calculation containing all of the atoms relevant to the doping process.

Not surprisingly, however, the quality of first order based predictions can vary significantly, and it is expected that the second order derivative in Eq. 1 can improve the accuracy of predictions45 by introducing some response properties of the system. For example, the second order term includes variation of the alchemical potential with respect to nuclear charge,

[TABLE]

where $\partial_{Z_{I}}\rho({\bf r})$ corresponds to the electron density response to varying the nuclear charge at the doping atom $I$. There are various ways to calculate the electron density’s response which involve further computational effort, for this work we merely wish to estimate it in a qualitative fashion. As such, we find it useful to assume the existence of a correlation between the actual response and the Pearson’s local softness of the atom in the molecule, as measured by the local density of the highest occupied molecular orbital (HOMO) for electrophiles (such as protons), $\rho_{HOMO}$ 55.

II.1 Technical details and system setup

The dissociative adsorption of a water monomer on boron nitride doped graphene (BNDG) was calculated using DFT and VASP 5.3.251; 52; 53; 54. VASP uses plane-wave basis sets and projector augmented wave (PAW) potentials 56; 57 to model the core region of atoms. The PBE exchange-correlation functional58 is used throughout along with PBE PAW potentials and a plane-wave energy cut-off of 500 eV. Earlier work has shown that similar trends in terms of water dissociation are obtained with PBE, the hybrid B3LYP59; 60; 61; 62 functional, and the dispersion inclusive optB86b-vdW63; 64; 65 functional.21 The dissociative adsorption energy of water was found to be converged to 0.001 eV with a plane-wave energy cut-off of 500 eV when tested up to 800 eV. A $(7\times 7)$ unit cell of graphene is used, with four carbon atoms replaced by two boron and two nitrogen atoms. The dissociative adsorption energy of water is already converged with a ($5\times 5$) unit cell but using a larger cell provides more pathways for alchemical mutation of atoms. The separation between periodic images of the substrate in the z-direction is 10 Å; this achieves convergence of the adsorption energy of water to within 0.004 eV compared to a z-direction separation of 30 Å. Reciprocal space was sampled with up to $7\times 7\times 1$ k-points and the adsorption energy was found to be converged within 0.05 eV at the $\Gamma$-point. Hence, all calculations reported here were performed at the $\Gamma$-point.

The adsorption site in the substrate contains a pair of BN atoms in the surface and two adsorbed hydrogen atoms, as shown in Fig. 1. Doping and hydrogenating in this way has been shown previously to make the surface more reactive towards the dissociative adsorption of water21. Importantly, atoms other than carbon at the adsorption site remain unchanged and are not involved in any transmutations. The dissociative adsorption energy is defined as,

[TABLE]

where $E^{tot}_{ads/sub}$ is the total energy of the adsorption system, $E^{tot}_{sub}$ is the total energy of the substrate (with two hydrogen atoms adsorbed), and $E^{tot}_{ads}$ is the energy of the intact water molecule in the gas phase.

Four types of alchemical mutation routes between carbon, boron, and nitrogen are considered here, illustrated in Fig. 2. There are a set of paths associated with each route, where a path defines the starting and final states for a given transmutation. The initial state in each path contains a pair of BN atoms near the edge of the unit cell which can be involved in transmutation. Note that in all four alchemical routes, the graphene sheet is also hydrogenated and contains a second pair of BN atoms at the dissociation site, but these particular dopants are excluded from alchemical mutation. In two types of routes, referred to as BN pair 1 and BN pair 2, a pair of BN atoms are transmutated to different sites across the graphene sheet, as illustrated with examples in Fig. 2. These two routes are distinguishable due to the existence of two sublattices within graphene. In BN pair 1, the transmutating BN atoms occupy the same sublattice in graphene as the unchanging BN atoms at the dissociation site. Whereas in BN pair 2, the transmutating BN atoms occupy the other sublattice. The third type of route, B2C, refers to alchemical changes involving only the boron atom. Similarly N2C refers to the swapping of carbon atoms with nitrogen while keeping the boron atom fixed. In each type of route there are 94 possible paths for this unit cell size, such that we have validated a total of 376 paths for this study. Note that only two single point DFT calculations are needed to make alchemical predictions for a set of 94 paths.

Thanks to the geometrical similarity of graphene and hexagonal boron nitride, doping graphene with BN atoms has a small impact on the structure. The largest change in bond lengths upon relaxation of target systems was seen for boron-carbon bonds, which changed by up to 0.06 Å. The energy of relaxation gained from this is up to $\sim 0.3$ eV and does not alter the trends observed. This makes fixing the geometry in all calculations a reasonable approximation to begin with.

III Results

The PBE energy of water dissociation has been calculated for each transmutation path without relaxing the positions of the atoms and compared with the alchemically predicted dissociation energy. Fig. 3 shows scatter plots comparing these energies, for each alchemical route. It can be seen that the PBE adsorption energies range from $-0.3$ to $-2.8$ eV, revealing that the precise location of the dopants has a significant impact on the reactivity of the active site. The large range of adsorption energies for seemingly similar surfaces can be understood in terms of two chemical effects from the doping boron and nitrogen atoms, namely, resonance and induction. In the former, the non-bonding valence electrons of nitrogen partake in $\pi$ conjugation with p-states on carbon atoms. This has a long-range impact on the electron density of the surface and therefore the reactivity of the adsorption site. Second, the difference in electronegativity between boron, carbon, and nitrogen atoms leads to local inductive effects and this is likely to have a particularly large impact when the doping atoms are near the active site. Upon considering how well the alchemical derivatives capture this behaviour, it can be seen that the majority of predictions is good. There are, however, a number of outliers resulting in a poor R2 correlation coefficient of 0.14 for the BN pair 1 route. The R2 coefficients for the other alchemical routes are similarly unimpressive between 0.17-0.49, and in all cases there are clear outliers. In addition, the few outliers correspond to configurations with either the most or least favorable adsorption energies – and the predictive power of the first order alchemical derivatives is worse for the outliers with the less favorable adsorption energies. These potentially interesting configurations are considered in more detail in Section IV, but first it is important to avoid predicting misleading trends for the outliers. It follows that for an effective screening process, it would be better to identify outliers without further computational cost. In the following section it is demonstrated how that is possible using the HOMO in the initial states.

III.1 Filtering outliers using highest occupied molecular orbitals

Let us first consider doped graphene in which the substituent atoms and dopants have a mesomeric effect on the electronic structure of the surface, i.e. they have either an electron withdrawing or electron releasing impact. This effect resonates across the surface giving rise to mesomerically active and passive sites. The mesomeric role of atoms can be probed using a Bader charge density partition66 per atom of the HOMO charge density, which indicates the prominence of the HOMO at a given atom site. Atoms with charge density above a chosen cut-off value in the HOMO are considered mesomerically active, and those under are mesomerically passive. For a given path, the charge at the sites of mutation in the initial state can be summed to obtain a measure of the extent of mesomeric activity. This combined Bader charge and the corresponding relative absolute error (RAE) for each path is shown in Fig. 4. It can be seen that most paths have a RAE less than 0.01, whilst those which have substantial errors also have large HOMO charges associated with them. As a result, the partitioned HOMO charge can be used to eliminate the outliers. Note that the correlation is not direct, there are some paths with a high associated HOMO charge but small errors.

The use of HOMO charges can be demonstrated by comparing the quality of predictions for two sets of paths, defined by a cut-off in their combined HOMO charges. More specifically, paths with a combined HOMO charge higher than a given cut-off charge are referred to as mesomerically active, and those with a lower charge are referred to as mesomerically passive. Here, the cut-off charge is chosen as the lowest combined HOMO charge found in paths with a RAE $>0.1$. In this way, we knowingly class all paths with RAE $>0.1$ as mesomerically active, and paths with RAE less than 0.1 are classed as mesomerically passive. Using this hindsight classification, the cut-off charges for the four routes are 0.203, 0.192, 0.140, and 0.025 e/atom for BN pair 1, BN pair 2, B2C, and N2C, respectively. Later we discuss how the cut-off charge can be chosen a priori without a threshold RAE, but its usefulness is first demonstrated in Fig. 5. It can be seen that all outliers belong to the mesomerically active paths (see filled circles in Fig. 5). In addition, the mesomerically passive paths deviate less from the PBE calculated energies and are therefore better predicted than mesomerically active paths.

The effectiveness of this procedure is more clearly seen in Table 1 where the R2, Spearman’s rank coefficient (r${}_{\text{s}}$), mean absolute errors (MAE) and root mean square errors (RMSE) are reported for each route. The MAE and RMSE are an order of magnitude larger for paths involving mesomerically active sites compared to passive sites. The MAE for mesomerically passive sites is $\sim$0.03 eV for the BN pair routes and thus within the so called chemical accuracy ($\sim$0.04 eV) of the PBE adsorption energies. Similarly, the MAE for mesomerically passive paths in the B2C and N2C routes is only slightly larger ($\sim$0.05 eV). Interestingly, the errors are generally larger for N2C (see in Fig. 4(b) the comparison with B2C) and as a result a smaller charge cut-off was used based on the threshold RAE of 0.1. The larger errors for N2C may seem at odds with the very good r${}_{\text{s}}$ coefficient for both mesomerically active (0.89) and passive (0.92) sites. Indeed from Fig. 5(d), it can be seen that there is only one obvious outlier in the N2C route. However, it has been shown previously that predictions for right-to-left transformations in the periodic table are not equivalent to the reverse and entail larger errors 45. We see this in the N2C route in which a nitrogen atom takes the place of different carbon atoms across the surface. Encouragingly, a strong correlation is still present between alchemically predicted and PBE calculated adsorption energies in the N2C route, despite a general shift away from the calculated energies.

IV Discussion

Partitioning the HOMO charge density for the initial states is shown to be an effective means of filtering out particularly weak predictions. However, at least two important questions need to be addressed with regards to this process and let us also draw some chemical insights.

First, how should the initial threshold value for the HOMO charge density be chosen without performing further calculations? This is somewhat of an arbitrary choice but some guidelines can be used. For example, the threshold charge can be chosen by considering the distribution of combined HOMO charges of all paths, and finding the point at which the combined HOMO charge begins to deviate from the majority of paths. For example, without considering the RAE and focusing only on the spread of values in the combined Bader charge in Fig. 4, most combined charges are below 0.15 e/atom. Importantly, this choice does not rely on a knowledge of the direct PBE results and interestingly, it is comparable to the informed choice of threshold values for the BN pair and B2C routes in Table 1. Indeed, according to Fig. 4 a threshold value of 0.15 e/atom would still correspond to small errors for all routes considered.

Second, how can the filtered mesomerically active paths be salvaged? In the current context that would be very useful because the most negative dissociation energies arise from doping at mesomerically active sites (see Fig. 4). Two particular solutions can be pursued. One is to simply perform DFT calculations for the mesomerically active paths - this is somewhat unimaginative but straightforward. The second possibility is to go beyond the first order alchemical derivative, and improve the prediction by including second order terms. Recently Chang et al. compared three approximations to the second order term namely, the coupled perturbed (CP) approach, the independent particle approximation (IPA) and the finite difference method, for the density response to alchemical coupling 45. The CP approach is shown to be superior to IPA for horizontal isoelectronic transformations in many-electron systems. However, all higher order alchemical derivative terms require additional computational cost. As such, it depends on the implementation of second order derivative approaches whether they would be more efficient than to directly calculate the DFT energies for mesomerically active paths.

Beyond the implications of efficiently screening isoelectronically doped configurations of graphene, one can take a closer look at the resulting favorable dissociative adsorption configurations to gain some chemical insight. Fig. 6 shows the configuration with the most favorable water adsorption energies obtained from alchemical predictions as well as direct PBE calculations, for each route, which range from $-2.1$ to $-2.8$ eV. Despite different starting sublattices for BN pair 1 and BN pair 2, the same configuration is identified as the most favorable for water dissociation. In this state, the hydrogen atom of water adsorbs on a carbon atom between two nitrogen atoms. This is not surprising given that the central carbon atom becomes more positive as a result of the electronegative nitrogen atoms, and is stabilized by bonding to a hydrogen atom. This is in agreement with the patterns identified in a previous DFT study21. Similarly, in the configuration found for N2C, the boron atom is between two nitrogen atoms and thus forms a stronger bond with the OH fragment of water. More interestingly, the favorable configuration from B2C is less intuitive, with a boron atom that is not directly bonded to a surface atom at the active site (see Fig. 6) and yet it corresponds to an adsorption energy of $-2.1$ eV. The reason behind the large negative adsorption energy for this peculiar configuration is still not fully understood and highlights the usefulness of screening through various topological possibilities. Whilst these adsorption energies are large and exothermic, indicating that they are likely to be observed in experiments, it is nonetheless important to also compute activation barriers. Note that although this is outside the scope of this study, alchemical derivatives can be used to predict activation barriers as previously shown42. Let us also note that the most unfavorable adsorption configurations found through the alchemical screening involve nitrogen-nitrogen single covalent bonds in the surface, suggesting that water adsorption on such sites is particularly unfavorable.

The exothermic adsorption energies of a water molecule on BN doped graphene found in this study span a remarkably large range ($>1$ eV) as a result of BN doping at various sites in the surrounding graphene sheet. Given that BN of doping graphene has been achieved experimentally33, and in-plane mixtures of graphene and h-BN have also been produced67; 68, it would be particularly interesting to verify our findings with experiments. In addition, hydrogenation of graphene has also been experimentally achieved (see e.g. Ref. 69), and thus it is timely to explore the thermochemistry at BN doped and hydrogenated graphene surfaces with adsorption measurements and surface studies.

V Conclusion

It has been shown that predictions using alchemical derivatives in DFT can be used to explore the impact of isoelectronic doping in activated graphene on the dissociative adsorption of water. Doping at different sites around the adsorption site in the substrate leads to a spread of $\sim$2 eV in the adsorption energy. Such a wide spread of adsorption energies shows that doping away from the dissociative adsorption site on graphene can have a significant impact on the adsorption energy of water. This suggests that BN doping of graphene could be a potential method for tuning surface reactions. Importantly, it has been demonstrated that poor alchemical predictions can be filtered out by identifying mesomerically active and passive sites using a Bader analysis of the HOMO in the initial state. In this way, one can efficiently screen through the majority of configurations with very good accuracy. For instance, in this study the MAE is as low as 0.025 eV in the dissociative adsorption energy of water. This corresponds to less than $1\%$ error for hundreds of PBE dissociative adsorption energies using minimal computational effort (eight self-consistent field DFT calculations). The use of alchemical derivatives for screening can also provide an efficient way to study the adsorption of various industrially important molecules such as hydrogen and methane on doped graphene surfaces. More broadly, there is scope for going beyond the first step in this study and screening materials and adsorbates in complex catalytic processes with alchemical derivatives. Such pre-screening could significantly reduce the number of DFT calculations that would need to be performed whilst providing useful chemical insight.

Acknowledgements.

We are grateful for support from University College London and Argonne National Laboratory (ANL) through the Thomas Young Centre-ANL initiative. O.A.v.L. acknowledges funding from the Swiss National Science foundation (No. PP00P2$\_$138932). This research was partly supported by the NCCR MARVEL, funded by the Swiss National Science Foundation. Some of the research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement number 616121 (HeteroIce project). A.M. is supported by the Royal Society through a Wolfson Research Merit Award. In addition, we are grateful for computing resources provided by the London Centre for Nanotechnology and Research Computing at University College London.