© P. Redondo et al. 2020

1 Introduction

Although water ice has been found throughout the cold interstellar medium (ISM; Merrill et al. 1976; Whittet et al. 1988; Dartois et al. 1998; Gibb et al. 2004; Boogert et al. 2008), the chemical origin and longevity of the icy grain mantles in which H₂O is the dominant component is still debated. The observed abundance of water ice cannot be explained by a simple freeze-out of H₂O onto the grain surface, and the persistent survival of the mantles under an intense radiation bombardment raises a question. Because the molecule formation in the gas phase is not very efficient, on-grain processes have to play a major role in the survival of the ice mantles. These on-grain mechanisms have been proposed in chemical models in the 1980s and are still discussed today (Tielens & Hagen 1982; d’Hendecourt et al. 1985; Hasegawa & Herbst 1993; Cuppen & Herbst 2007). The three reaction pathways proposed by Tielens & Hagen (1982), that is, hydrogenation of atomic oxygen, molecular oxygen, and ozone, are not only at the origin of the chemical models but also inspired most of the laboratory experiments. The last steps of the reaction network, which are those that effectively yield H₂O, are $$\text{H+OH}\to {\text{H}}_2\text{O}$$(1) $${\text{H}}_2\text{+OH}\to {\text{H}}_2\text{O+H}$$(2) $$\text{H}+{\text{H}}_{\text{2}}{\text{O}}_{\text{2}}\to {\text{H}}_2\text{O+OH}.$$(3)

Reactions (1) and (2) appear as the final steps in the hydrogenation of the oxygen atom; they also appear in the hydrogenation of ozone after OH has been formed by $$\text{H}+{\text{O}}_{\text{3}}\to {\text{OH+O}}_{\text{2}}.$$(4)

Reaction (3) is effective in the hydrogenation of molecular oxygen after H₂O₂ has been obtained following $${\text{H+O}}_{\text{2}}\to {\text{HO}}_{\text{2}}$$(5) $${\text{H+HO}}_{\text{2}}\to {\text{H}}_2{\text{O}}_{\text{2}}.$$(6)

All these mechanisms have been investigated in the laboratory. For instance, experimental evidence for water formation through ozone hydrogenation (4) onto dust grain analogs has been obtained by Mokrane et al. (2009). The water formation by H hydrogenation of an O₂ surface has been probed by Ioppolo et al. (2010) and Cuppen et al. (2010). These studies also confirmed the implication of reaction (3) along the way (Hiraoka et al. 1998; Kouchi et al. 2008). However, although presumed to exist in the icy mantle of interstellar grains, neither O₂ nor O₃ has been identified so far.

On the other hand, the OH radical appears at one stage or another in all mechanisms of H₂O formation, and we recall that a number of these radicals exist on the ice cover or around it as the incident cosmic irradiation breaks the H₂O bonds, which produces a copious amount of radicals such as O, H, and OH. Here we focus on the $$\text{OH+OH}\to {\text{H}}_2\text{O+O}$$(7)

reaction that has not been considered properly before, probably because of an activation in the gas phase (~ 2 kcal mol⁻¹, according Karkach & Osherov 1999) and processes along a triplet surface reaction, which is not the reaction with the lowest energy.

The aim of this research is to assess the possible role of the ice in the formation, or more precisely, the re-formation, of the ice mantle. In order to lead to the formation of a water molecule from two interacting OH radicals, reaction (7) has to proceed in a triplet state, regardless of the surface underneath (the end products, H₂O(¹A₁) and O(³P), are both in their ground states). The other possible reaction starts in the singlet state of the [HO...OH] complex, which implies the coupling of the two OH radicals and leads to hydrogen peroxyde H₂O₂. This is onlyone more way to start reaction (3), which is taken into account in existing chemical models. Although the oxygen atom does not participate as a reactant in the reaction OH + OH → H₂O + O and appears only as a final product when the reaction is completed, it might be wondered whether the presence of water molecules would reverse the triplet/singlet stability ordering of atomic oxygen (Molpeceres et al. 2019). In the gas phase, atomic spectroscopy has shown that the singlet state is less stable than the triplet by 1.97 eV (Kramida et al. 2019). A test of atomic oxygen on a layer of 24 H₂O molecules cut in solid water ice showed a singlet-triplet gap of ~2 eV at the MPWB1K level (Zao & Truhlar 2004). Our density functional theory (DFT) periodic calculations (often referred to as “first principle” treatments) on a solid surface of infinite dimension give ~3 eV. This great difference in energy means that a triplet-to-singlet swap is unlikely at the end of a reaction that occurs entirely on the triplet potential surface.

Here we present an extensive survey of reaction (7). This survey starts from a gas-phase situation intended to serve as reference, continues in the presence of an increasing number of H₂O molecules; then a bilayer is considered that mimics crystalline ice, and we conclude with a real solid surface of infinite dimension. We summarize our conclusions in Sect. 7.

2 Theoretical background

We investigate the possible role of ice in its own reconstruction through reaction (7). To achieve this, it is necessary to accurately describe the chemical process at the molecular level for simple species such as OH radicals or small aggregates and simultaneously at the level of the solid water that supports the reaction. These are two different domains for which elaborate and conceptually different theories have been developed to account for changing electronic correlation/dispersion effects in the course of chemical reactions. It is particularly important when weak interactions are involved. In both cluster and solid-state approaches, the reconstruction of water ice has been considered to proceed in an overall triplet state. In other words, all structures, from the pre-reactant complex [OH… OH] and transition state to the final products H₂O + O(³P), are calculated in a way that satisfies the triplet multiplicity.

2.1 Cluster approach

In ab initio post-Hartree-Fock (HF) methods designed for small-dimension systems (here simple models where the ice is represented by small H₂O aggregates), the effects of electronic correlation are treated using a systematic progression in the number of electronic excitations from the HF reference state by means of perturbative treatments for double excitations (MP2) and by high-level coupled cluster calculations for double and triple excitations (CCSD and CCSD(T); Bartlett & Shavitt 1977; Raghachavari et al. 1989). Atomic basis sets of increasing levels of flexibility, from 3–21+G(d,p) to correlation consistent cc-pVnZ of double, triple, and quadruple (n = 2–4) quality, possibly augmented with diffuse functions, were employed to perform the calculations.

Calculations for much larger systems (extended models where the ice is represented by aggregates of dozens of H₂O molecules) are the domain of the DFT. We considered the hybrid formalisms that are used most frequently, B3LYP and BHandHLYP, as well as the functionals that were recently improved to account for dispersion effects, PBE1PBE, MPW1K, and MPWB1K (Lynch et al. 2000; Zao & Truhlar 2004). The basis sets were the same in ab initio and DFT calculations.

We note a last technical point: the zero-point vibrational energies (ZPE) were included in every building of the potential energy surfaces. Because analytical frequencies cannot be calculated at the CCSD(T) level in the code we employed, we instead used the force constants that are available at this level by double numerical differentiation. Hence, harmonic vibrational frequencies were calculated at this level numerically by the standard procedure. For all other single points, the vibrational frequencies were calculated analytically at the level at which the geometries were optimized. We performed these calculations using levels of theory and basis sets as implemented in the GAUSSIAN09 package (Frish et al. 2009).

2.2 Solid-state approach

The approach presented above might suffer, a priori, from inherent dimension constraints that might introduce some distortion in the results. To avoid these limiting factors, a solution is to set up the icy environment with a proper solid surface of infinite dimension. Based on our preceding works on molecular adsorption on water ices (Lattelais et al. 2011, 2015; Bertin et al. 2017), we took a periodic representation of solid hexagonal ice I_h, composed ofbilayers of water molecules. The ice surface was described by the (0001) Miller indices. This surface was shown to be the most stable (Pan et al. 2010). This surface also satisfies the Bernal & Fowler (1933) constraint: exactly two hydrogens of neighboring water molecules form hydrogen bonds with oxygen lone pairs of the adjacent water molecules. More details can also be found in previous work on the structures of ice (Casassa & Pisani 2002; Calatayud et al. 2003; Casassa et al. 2005).

The electronic treatment relies on the generalized gradient approximation (GGA) in the form of the PBE functional. Long-range van der Waals and hydrogen-bonding interactions are better described within the corrective scheme (PBE+D2) of Grimme et al. (2010). The core electrons were kept frozen and were replaced by pseudo-potentials generated by the plane-augmented wave (PAW) method.

Cluster calculations show that the functional must include the exact HF exchange. We performed solid-state calculations of the molecular two-body OH + OH → H₂O + O case system using a cell large enough to prevent any interaction between the molecules in adjacent cells. To make the solid-state description comparable to that of the MPWB1K hybrid DFT approaches we employed for cluster calculations, a 50% percentage of the exact HF exchange was added to the PBE+D2 functional.

These solid-state periodic DFT treatments, also known as first-principle calculations, were carried out using the Vienna ab initio simulation package (VASP; Kresse & Hafner 1994; Kresse & Furthmüller 1996; Kresse & Joubert 1999). More details on the computational procedure using VASP are given in the appendix.

Fig. 1 Reaction profile for the gas-phase reaction OH(2Π) + OH(2 Π) → H2O(1A1) + O(3 P). Energies in kcal mol−1 are from CCSD(T)/aug-cc-pVTZ (including ZPE). Less expansive calculations (CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ are given in parentheses.

3 Two-body OH + OH reaction as a case study

An extensive study of the reaction profile of reaction (7) was carried out. We anticipated calculations on much larger systems andtherefore used multiple levels of wave functions, from highly correlated post-HF to DFT treatments using a large panel of functionals.

The reaction pathway illustrated in Fig. 1 from the values listed in Table 1 shows the same stationary points at all levels of theory. A shallow complex [OH...OH] forms in the entrance channel, followed by a transition state whose position with respectto the starting species is largely dependent on the level of theory, and another shallow complex in the exit channel [O...H₂O] that is close in energy to the isolated products H₂O and O(³P).

We found that the entire reaction proceeds in the same plane, in agreement with a previous study (Karkach & Osherov 1999) at the QCISD(T)/6-311(d,p) level. The reaction is exothermic with a value of approximately −19 kcal mol⁻¹ at the CCSD(T) and QCISD(T) levels. This value is reduced to approximately −16 kcal mol⁻¹ for coupled cluster calculations when zero-point energies are taken into account. At the same level of theory, the entrance channel complex is stabilized by approximately −2.4 kcal mol⁻¹ and the exit channel complex is found to be approximately −1.2 kcal mol⁻¹ below the separated products; the transition state is ~1.6 kcal mol⁻¹ above the two separated OH radicals, that is, ~4 kcal mol⁻¹ above the [OH...OH] complex.

Compared to the highest level of theory (CCSD(T)/aug-cc-pVTZ), the exothermicity of the reaction is overestimated at MP2 and underestimated at the DFT level, except for the MPW1K and MPWB1K functionals. The complexes on the reaction pathway are also best represented by the MPW1K and MPWB1K functionals. The energies given by the MP2 calculations, even with augmented basis sets, are too high. This suggests that higher excitations are needed to determine reliable activation energies. However, weexpect that the geometry obtained at this level is correct; when we used the MP2 geometries, we found about the same activation energy in the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ as in the fully optimized CCSD(T)/aug-cc-pVTZ treatment¹.

Although the list of functionals we tested is far from being complete, the DFT approach to be used in forthcoming studies appears essentially limited to MPW1K and MPWB1K functionals (augmented with diffuse functions). We verified and confirmed this point inthe following developments by comparing it to an increasing complexity of the reactive environment. It is worth mentioning that B3LYP, PBE1PBE, and MPWB1K without diffuse functions are inappropriate enough to lead to artificial hidden barriers (i.e., below the reactants).

4 Test evolution when H₂O molecules are added

After reconsidering the gas-phase reaction to determine what level of theory is appropriate for the calculations on larger systems, we tested the potential energy surface of reaction (7) in the presence of an increasing number of H₂O molecules atthe CCSD(T) and DFT levels. We study the influence of the icy support here, therefore we use as starting point reference the situation where one OH is attached to the support (from one H₂O molecule to the full solid). In this context of reactivity on the surface, we consider as final point the complex in which the oxygen is still attached to the support.

Fig. 2 Reaction profile for OH(2Π) + OH(2 Π) → H2O(1A1) + O(3 P) at the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ level in presence of one H2O, ZPVE included (relative energies are in kcal mol−1).

4.1 Potential surface at CCSD(T) level

Starting with one H₂O support, we obtained the reaction profile illustrated in Fig. 2 from the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ energies reported in Table 2. As a consequence of the presence of a witness H₂O, two different complexes are now possible on the entrance channel. These complexes [H₂O...HO] and [HO...H₂O] correspond to the two possibilities of hydrogen bonding, or between the H atom of the OH radical and the oxygen lone pair of the H₂O witness, giving the most stable of the two complexes (which is used as energetic reference), or between one of the H atoms of the H₂O witness and the oxygen lone pair of the OH radical, 1.65 kcal mol⁻¹ above in energy. Then, a possible three-body complex is found before the activation barrier, 2.6 kcal mol⁻¹ below the reference starting point. It is followed by a transition state whose position on the energy scale appears to be 1.7 and 0.1 kcal mol⁻¹ below the starting [HO...H₂O] and [H₂O...HO] complexes, respectively. The final product in this model reaction consists of an oxygen atom linked to the H₂O dimer. This atom can be removed from the support later.

Taking a support of two H₂O in the form of the dimer obtained above, we obtained the reaction profile illustrated in Fig. 3 from the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ energies reported in Table 3. The starting point of the reaction was a triangular complex, [H₂Odimer...OH + OH], in which one OH radical makes two hydrogen bonds with the two H₂O molecules forming the water dimer. The reaction progresses by accreting the second OH radical in the form of a quadrangular complex that is better represented by the OH...OH complex found for the gas phase that interacts at both ends with the water dimer. The next stationary point, that is, the transition state, follows from the latter complex with a structure similar to that of the gas-phase transition state (Sect. 3) that interacts with the supporting water dimer. It is found 1.2 kcal mol⁻¹ below the starting point and 3.7 kcal mol⁻¹ above the intermediate complex. We confirmed that this former triangular complex also correlates with another transition state that is much higher in energy (~13.2 kcal mol⁻¹) above the reference energy, which makes it a highly unfavorable pathway (not shown). The final product is a complex between the oxygen atom and the water trimer that is formed as before.

With a support of three H₂O added in the form of the optimized trimer obtained above, we found a reaction profile different from those previously obtained, that is, the trend that the activation barrier is lowered while the number of water molecules increases is no longer true for n = 3. We found a real activation barrier of 1.7 kcal mol⁻¹ on the way to a cubic structure that is the most stable form. This agrees with the result obtained by Wales & Hodges (1998) for n = 4 in a systematic study of the structure of the ice cluster (n = 2−21) H₂O molecules. This structure does not reproduce even a partial geometry of solid ice, however. We therefore discontinued this type of modeling for ice, that is, optimized clusters. However, the indication given by the former very small clusters remains valid as a trend.

Fig. 3 Reaction profile for OH(2Π) + OH(2 Π) → H2O(1A1) + O(3 P) at the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ level in presence of H2O dimer, ZPVE included (relative energies in kcal mol−1).

4.2 Potential surface at the DFT level

Knowing that CCSD(T) is not applicable to larger models, we reconsidered the above potential surfaces at the MPW1K and MPBW1K that provided the best approaches to the reference in the gas phase (Table 1). With one H₂O support (see Table 2), the MPW1K/ MPWB1K functionals provide a fair approach to the CCSD(T) exothermicity and activation energy, with a slight advantage for MPWB1K.

In the case of two H₂O support (seeTable 3), as in the studies we referred to above, the results obtained with MPW1K/ MPWB1K are similar to those of the CCSD(T) calculations for both the exothermicity and the activation energy. The trend is identical: the result for MPWB1K is slightly better.

The same drawbacks as in CCSD(T) appear when the H₂O trimer is the starting complex that reacts with the OH radicals.

4.3 Modeling ice with small clusters

We found roughly the same energy profile with similar structures for each of the stationary points. The two approaches we used, CCSD(T) and DFT (especially MPWB1K) give similar values for the exothermicities and activation barriers.

The parameter that controls the facility for the OH + OH reaction to proceed is the activation energy. It decreases from its gas phase value with the addition of one or two molecules of water. However, the fact that this tendency, which is valid for CCSD(T) and DFT, is not followed when more H₂O is added to the support is consistent with the results of Wales & Hodges (1998), who reported a systematic study of the structure of ice clusters with increasing dimension (up to 21 H₂O molecules). From the H₂O trimer on, none of them can be taken as a substructure of the most stable I_h compact ice. We add that numerically extrapolating the interaction energy between one single H₂O and a free-flyer molecule to an adsorption energy on ice can lead to erroneous conclusions, especially when more than one adsorption site is possible (Vakelam et al. 2017). Consequently, the approach based on the successive growth of simple water clusters can only give an indication of trends that has to be confirmed by a more elaborate model.

5 Reaction on water layers

To approach the ice situation while keeping a tractable model support, an H₂O aggregate of reasonable dimension can be obtained by an appropriate cut of 24 H₂O molecules (24–1L) assembled in the form of six fused cyclohexane-like structures that might represent a significant part of the surface bilayer of the crystalline ice (I_h). An even more realistic description can be obtained by superimposing two of these bilayers (48 H₂O molecules) to model the role of the underlying bulk. Both clusters satisfy the Bernal & Fowler (1933) constraint.

With the MPWB1K functional we previously selected, we can study the reconstruction of the ice surface following the reaction path in additon to the ice model. This is not possible at the highly correlated level we used previously for the treatment of the small aggregates. We anticipated calculations on (48–2L) clusters and therefore employed two basis sets of different dimensions for the (24–1L) cluster: aug-cc-pVTZ as in the small cluster calculations, and 6–311++G(3df, 2p) to investigate alternative solutions of lower cost. Because the reaction profiles and stationary points are quite similar for the 24(1–L) and 48(2–L)ice supports, we present only one figure as illustration with snapshots to show the structural details of each critical point (Fig. 4).

(i) The first test was calculating OH + OH on the single bilayer (24–1L). We were aware of the problems of an eventual flexibility of the structure (as in the 3 H₂O case) and therefore kept the bilayer frozen in order to keep an ice-like structure. Only the reactants were optimized along the reaction path.

The reaction starts from an OH radical on the ice surface in the form of a bridged complex in which OH is attached to the 1, 3 positions of the hexagonal surface in a chair conformation and was taken as energy reference for the process we considered. The next stationary point is also a bridged complex in which the incoming OH forms an H₂O dimer with the OH that was previously adsorbed on the surface. This complex is stabilized by ~12 kcal mol⁻¹ below the energy of the starting products. The next step on the reaction path is a transition state that is ~10 kcal mol⁻¹ above the bridge complex, but still ~2 kcal mol⁻¹ below the reactants energy. This confirms the trend we previously found with the clusters, that is, the presence of a hidden barrier. This means that there is no activation energy to oppose the continuation of the reaction. The final product is the complex composed of the original bilayer increased by one H₂O molecule with the oxygen atom attached to the surface. In this position, the molecule is well placed to react with any diffusing hydrogen atom to form a new OH radical and reinitiate the whole mechanism. The two basis sets lead to practically the same values within the uncertainty margin (Table 4). This justifies the use of 6–311++G(3df, 2p) for the next larger calculations.

(ii) A second calculation was performed with the two-bilayer structure (48–2L). The subsurface frozen bilayer is intended to maintain the structure that would be imposed by the bulk of real ice. By imposing this constraint, we leave the possibility for relaxation in the upper bilayer by optimizing the six nearest H₂O molecules that form the cyclohexane-like environment that supports the reaction without distorting the whole structure dramatically. The reaction profile of OH + OH supported by a two-bilayer of water is very similar to the profile that is obtained with one bilayer. The energy values reported in Fig. 4 are those where the ice support is optimized.

The reaction starts with an incoming second OH radical that leads to the formation of a bridged complex (~7.5 kcal mol⁻¹ below the reference energy) in which the OH...OH dimer is attached to the 1, 3 positions of the hexagonal surface in a chair conformation. The next stationary point is a transition state whose structure is shown in Fig. 4, situated between 1 and 2 kcal mol⁻¹ below the reactants level. It shows again that we are in presence of a hidden barrier on the reaction path and that there is no activation energy to oppose the continuation of the reaction. As with one bilayer, the final product is a structure composed of the original layer augmented by a H₂O molecule to which the oxygen atom is attached at the surface. This system can either react with any diffusing hydrogen atom to form a new OH radical and then reinitiate the reconstruction mechanism of the ice mantle, or it can detach from the ice.

Fig. 4 Schematic top view of the reaction profile for OH(2Π) + OH(2 Π) → H2O(1A1) + O(3 P) in presence of an ice double bilayer at the MPWB1K/6-311++G(3df,2p). Energies are given in kcal mol−1. Details of the evolution of the reaction sites during the reaction are given on snapshots of the topmost bilayer for each step of the energy profile.

6 Reaction on solid ice

We finally configured the icy environment with a solid surface of infinite dimension. In this periodic approach, the model unit cell that is cut into apolar hexagonal ice I_h contains twoH₂O bilayers. It is an elementary volume whose replication in three spatial directions creates a series of ice surfaces of infinite dimensions. This volume is defined by three structural parameters a, b, and c. Here, the basal plane dimensions (corresponding to the surface on which the reaction proceeds) are a = 7.145 Å and b = 8.73 Å. The third dimension of the cell was determined so that there is no interaction between molecules generated by a translation perpendicular to the ice surface. The same reaction path as for clusters (see Fig. 4) was considered, starting from one OH radical adsorbed on the ice model surface of infinite dimensions, taken as reference. As in the cluster case, the second OH is a free flyer, and we follow the bridge conformation throughout, as for the 48(2–L) simpler model.

The trend for stabilization shown by the interaction energy of the transition state with the ice support when n(H₂O) increases is even emphasized. It shows a hidden barrier of −3.8 kcal mol⁻¹ below the reactants, that is, deeper than that found in the cluster case. This confirms the hypothesis of a possible ice reconstruction through this mechanism.

Snapshots of the energy profile of the solid-state reaction mechanism are given in Fig. 5. The critical energies are reported in Table 5 together with the results of the preceding calculations for comparison.

Fig. 5 Views of the reactive sites of OH + OH → H2O + O on the topmost ice bilayers of periodic ice Ih. Reactants (left), transition state (middle), and products (right).

7 Summary and concluding remarks

The irradiation of the interstellar ices is known to break the water molecules into three fragments: atomic hydrogens, oxygens, and OH radicals. The main goal of this work was to increase the chemical information on the reconstruction mechanisms of ice mantles of the interstellar grains at the molecular level. The direct recombination of H and OH fragments (reaction (1)) was originally proposed for the formation of the ice. However, the hydrogen atoms produced by the diffusing irradiation easily give mostly molecular H₂. These hydrogen molecules could also react with the OH radicals (reaction (2)), and it has been shown that this reaction could proceed through the tunneling effect (Oba et al. 2012; Meisner et al. 2017).

The reconstruction of the ice under the hypothesis that the OH radicals produced by destructive irradiations might recombine through reaction (7) at the top of the remaining ice was considered. Two types of numerical simulations based on conceptually different approaches were employed. The first approach relies on individualizing the ice support down to the molecular level using high-level ab initio calculations (CCSD(T) and recently developed DFT functionals. The second approach relies on globalizing the ice support as a solid of infinite dimensions using a first-principle periodic calculation. The two approaches give a similar qualitative answer: the results for large H₂O molecular clusters are similar to the results we obtain for the solid surface that has infinite dimensions.

We found that in both approaches, no activation energy prevents the pursuit of the H₂O synthesis at the surface of water ice. In the present case, where we constrained the top bilayer to respect the structure of the underlying bulk, the new H₂O molecule that formed is well positioned to help in the reconstruction of the next bilayer. The oxygen atom released in the process is also accessible to diffusing H radicals or to other H atoms that are present in the environment. The first-principle periodic calculation carried out finally agreed with the large cluster treatments: the reconstruction of the ice mantles by reaction OH + OH → H₂O + O is a barrierless process when it is catalyzed by ice support.

Acknowledgements

This work was supported by CNRS national program PCMI (Physics and Chemistry of the Interstellar Medium) and COST Action CM 805, (The Chemical Cosmos). P.R. acknowledges the funding support from the Spanish Ministerio de Economia Industria y Competitividad (Grant AYA2017-87515-P) and the Junta de Castilla y León (Spain; Grant No. VA010G18).

Appendix A VASP keywords

Keywords for using the PBE+D2 functional with 50% of exact HF exchange in VASP:

GGA = PE

IVDW = 1

LHFCALC =.TRUE.

AEXX = 0.5

AGGAX = 0.5

AGGAC = 1.0

ALDAC = 1.0

PREFOCK = Normal

ALGO = All

TIME = 0.4

For more insight into the Vienna Ab initio Simulation Package, see https://cms.mpi.univie.ac.at/vasp/vasp/Hartree_Fock_HF_type_hybrid_functional_ calculations.html

References

All Tables

All Figures

	Fig. 1 Reaction profile for the gas-phase reaction OH(2Π) + OH(2 Π) → H2O(1A1) + O(3 P). Energies in kcal mol−1 are from CCSD(T)/aug-cc-pVTZ (including ZPE). Less expansive calculations (CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ are given in parentheses.
In the text

	Fig. 2 Reaction profile for OH(2Π) + OH(2 Π) → H2O(1A1) + O(3 P) at the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ level in presence of one H2O, ZPVE included (relative energies are in kcal mol−1).
In the text

	Fig. 3 Reaction profile for OH(2Π) + OH(2 Π) → H2O(1A1) + O(3 P) at the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVDZ level in presence of H2O dimer, ZPVE included (relative energies in kcal mol−1).
In the text

	Fig. 4 Schematic top view of the reaction profile for OH(2Π) + OH(2 Π) → H2O(1A1) + O(3 P) in presence of an ice double bilayer at the MPWB1K/6-311++G(3df,2p). Energies are given in kcal mol−1. Details of the evolution of the reaction sites during the reaction are given on snapshots of the topmost bilayer for each step of the energy profile.
In the text

	Fig. 5 Views of the reactive sites of OH + OH → H2O + O on the topmost ice bilayers of periodic ice Ih. Reactants (left), transition state (middle), and products (right).
In the text