H 2 O 2 (s) and H 2 O 2 ·2H 2 O(s) crystals compared with ices: DFT functional assessment and D3 analysis

The H 2 O and H 2 O 2 molecules resemble each other in a multitude of ways as has been noted in the literature. Here, we present density functional theory (DFT) calculations for the H 2 O 2 (s) and H 2 O 2 ⋅ 2H 2 O(s) crystals and make selected comparisons with ice polymorphs. The performance of a number of dispersion-corrected density functionals—both self-consistent and a posteriori ones—are assessed, and we give special attention to the D3 correction and its effects. The D3 correction to the lattice energies is large: for H 2 O 2 (s) the D3 correction constitutes about 25% of the lattice energy using PBE, much more for RPBE, much less for SCAN, and it primarily arises from non-H-bonded interactions out to about 5 Å.The large D3 corrections to the lattice energies are likely a consequence of several effects: correction for missing dispersion interaction, the ability of D3 to capture and correct various other kinds of limitations built into the underlying DFT functionals, and finally some degree of cell-contraction-induced polarization enhancement. We find that the overall best-performing functionals of the twelve examined are optPBEvdW and RPBE-D3. Comparisons with DFT assessments for ices in the literature show that where the same


I. INTRODUCTION
The measured macroscopic properties of H 2 O and H 2 O 2 (hydrogen peroxide) are similar in a multitude of ways, from the molecular dipole moments in the gas-phase to the heat capacities of the liquids to the crystalline melting points.At the microscopic level, the bonding patterns of hydrogen peroxide and water are also similar: given the chance, H 2 O 2 and H 2 O each donate two and accept two hydrogen bonds and the H-bond interactions are roughly of similar strengths.A recent Raman spectroscopic study by the group of Ben-Amotz highlights such similarities by the phrase "seemless integration of hydrogen peroxide in water" (Ref.1).At the same time, there are clearly crucial differences between hydrogen peroxide and water, where the redox-active properties of the peroxide ion as a member of the ROS (reactive oxygen species) family 2,3 is just one example.Also the current paper combines hydrogen peroxide and water in condensed matter, but here for the crystalline state.The H 2 O 2 (s) and H 2 O 2 ⋅2H 2 O(s) crystals are in focus, and we will make frequent comparisons with selected ice polymorphs.First, the performance of a number of dispersion-inclusive and dispersion-corrected density functionals are assessed for the two peroxide-containing target crystals, and the result is compared with similar assessments of functionals for ice polymorphs in the literature.Then, we zoom in on the Grimme a posteriori dispersion corrections 4 for these crystals and explore to what extent the D3 correction captures not only London dispersion contributions to the lattice energies but also other energy components.The main objectives of this paper are thus to present the DFT functional assessment and to discuss the implications of the D3 London dispersion correction in relation to some common categories of intermolecular interaction: dispersion (obviously) but also polarization and H-bonding.Regarding the first objective, some commonly used DFT functionals, representative of different DFT families along Jacob's ladder of exchange-correlation functionals, 5 with and without dispersion treatment, will be assessed.An exact ranking order (1-2-3-. ..) is not the main purpose; the intent is rather to examine to what extent standard DFT methods work well for the peroxide crystals and whether functionals that have been found to perform well for ice polymorphs in the literature also work well for our peroxidecontaining crystals.Indeed, there exist several DFT benchmarking and computational method assessment papers for the ices in the literature, such as those of Santra et al., 6 Brandenburg et al., 7 Della Pia et al., 8 Gillan et al., 9 and others.In contrast, theoretical studies of H 2 O 2 -based condensed-matter systems are scarce, as evidenced by Table I, which constitutes a quite complete list of computational studies of condensed-matter systems involving H 2 O 2 .The list is short despite the fact that both force-field and electronic structure studies are included.

The Journal of Chemical Physics
Assessment efforts need access to benchmark values.For isolated molecules and small molecular aggregates, the golden standard for electronic calculations are high-level quantum-mechanical methods, such as coupled-cluster calculations with single, double, and perturbative triple excitations [CCSD(T)].For condensedmatter systems, the situation is more complicated as both short-and long-range interactions simultaneously need to be well described.For periodic molecular solids of modest complexity, diffusion quantum Monte Carlo calculations (DMC) constitute the theoretical gold standard; see, for example, the DMC study of three ice polymorphs by Santra et al. 10 and the DMC calculations by Zen et al. 11 for a set of eight small-molecule crystals, also including ice, and the assessment of fifty-one functionals against DMC-calculated lattice energies for thirteen ice polymorphs by Della Pia et al. 8 Examples of other high-level theoretical approaches are the Random Phase Approximation with exchange (RPAx) calculations for four ice polymorphs by Hellgren and Baguet, 12 the complete-basis-set (CBS) fragment The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcpCCSD(T) calculations performed for, e.g., solid CO 2 , 13 and fragment SAPT calculations for solid formamide by Bordner. 14However, for solid-state electronic and atomistic calculations, the most desirable benchmark values originate from low-temperature experiments of structure and other properties on well-defined and well-behaved crystals, whenever such data exist.Low-temperature, ambient-pressure diffraction data are available for our two hydrogen peroxide systems.However, matching experimentally and computationally derived quantities in a consistent manner can still be a challenge.This is discussed in Sec.III A (structure) and Sec.III B (lattice energy), and uncertainties in the experimental benchmark values and their consequences are discussed in Sec.III E.
In between these, Secs.III C and III D present our DFT assessment table for the peroxide-containing crystals and comparisons with DFT surveys for ices in the literature.The MARD (mean absolute relative deviation) ranges used in our color-coded grading scheme take the experimental uncertainties into account.
Section IV gives special attention to the D3 correction and whether it-in addition to the fluctuating-dipoles based London dispersion-also captures other types of intermolecular interaction, such as electrostatic polarization and hydrogen bonding and whether it captures other additional shortcomings of the underlying DFT functional.Section V provides some perspectives on which has the best H-bonding ability-water or hydrogen peroxide.

A. Systems
The two compounds in focus here both form single crystals, H 2 O 2 (s) in space group P4 For the H 2 O 2 (s) crystal, no less than three single-crystal low-and ambient-temperature diffraction studies have been published [15][16][17] with consistent results.Reference 17 is a lowtemperature neutron diffraction study and will be the provider of the structural benchmark data for H 2 O 2 (s) in the current article.For the H 2 O 2 ⋅2H 2 O crystal, there exists a high-quality x-ray diffraction study at low temperature, 18 which will be the provider of benchmark  Comparisons with four ice structures (an ordered variant of ice Ih and the essentially ordered structures ice II, ice VIII, and ice IX) are made for diverse purposes in the paper.For ice Ih, we only make use of literature results, but for the other ice polymorphs, we also optimized the structures (at zero pressure) using a few of the functionals employed for the peroxide crystals.Our use of optimized zero-pressure structures for all the ice polymorphs is consistent with how we treat the peroxide crystals.Furthermore, although ice II, ice VIII, and ice IX are only stable at elevated pressures, experimental zero-pressure structural data are available for these three polymorphs: metastable ambient-pressure-recovered structures of ice II 19 and ice IX 20 and a metastable structure obtained from extrapolation of measurements at 2.4 GPa and above for ice VIII (see Ref. 21 and also Ref. 22).The same procedure of validating DFT-optimized ice structures against experimental ambient-pressure (recovered or extrapolated) ice structures is a common practice in the literature (see, e.g., Refs.7 and 8).

B. Periodic plane-wave DFT calculations using VASP
Electronic spin-unpolarized calculations within the framework of periodic plane-wave density functional theory (DFT) calculations were performed for all systems (crystals and isolated molecules), all were optimized.The computational cells used for the crystals were the crystallographic unit cells (optimized in terms of cell parameters and atomic positions).The H 2 O and H 2 O 2 monomers were optimized in 15 × 15 × 15 Å 3 cells; the resulting geometrical parameters are listed in Table S1 in the supplementary material.A 16 × 16 × 16 Å 3 cell was used for the bi-molecular complexes mentioned in Secs.IV and V.
Twelve functionals were selected.They were chosen to represent some of the quite well-established functionals, with a range of exchange-correlation descriptions, namely, GGA, meta-GGA, and hybrid functionals.With respect to dispersion correction, the functionals belong to three families: the plain GGA functionals without dispersion correction (PBE, 23 RPBE, 24 SCAN, 25 and HSE06 26 ), with a posteriori dispersion corrections (PBE-D3, 4,23 RPBE-D3, 4,24 BLYP-D3, 4,27,28 SCAN-D3, 4,25 and HSE06-D3 4,26 ), and the non-local van der Waals-inclusive functionals where the dispersion contribution is included as an integral part of the exchange-correlation functional (optB86b-vdW, 29 optPBE-vdW, 30 and vdW-DF-cx 31 ).In this way, we investigate a few steps along Jacob's ladder of DFT functionals with different exchange-correlation expressions and also include some of the lowest steps on the stairway of dispersion-correction functionals. 32he k-point sampling used a 4 × 4 × 2 Monkhorst-Pack scheme for the H 2 O 2 crystal, a 2 × 2 × 2 Γ-centered grid for the H 2 O 2 ⋅2H 2 O crystal, and just the Γ-point for the isolated molecules.4][35][36] The 1s 2 electrons were part of the O core region.The valence-core electron interactions were represented by the PAW scheme, 37 and hard pseudopotentials were used for oxygen.In all cases, the plane-wave cutoff was 1000 eV.The cell parameters and atomic coordinates were consid-ered converged when the forces acting on the atoms were less than 0.002 eV/Å.For all five of the D3-corrected methods used here, we performed structural optimizations with both the zero-damping scheme, D3(0), and the Becke-Johnson damping scheme, D3(BJ), 38 for both crystals.Lattice energies were also calculated with the two damping schemes.The resulting volumes are compared with our experimental diffraction references (Refs.17 and 18) in Fig. S1 in the supplementary material, and calculated lattice energies with the two damping schemes are shown in the same figure (for the reference values, see the supplementary material).We find that the two damping protocols generally give very similar results for our two crystals and that neither scheme performs systematically better than the other.All results presented in the text, in the figures, and in Table II were calculated with the D3(0) damping scheme.
Vibrational zero-point energies (ZPEs) were calculated for the two hydrogen peroxide crystals using the lattice dynamics protocol inherent to the VASP program.The ZPEs were then used in our calculations of sublimation enthalpies, as detailed in Sec.III B.

TABLE II. Average errors (deviations from reference values) illustrating the performance of the different functionals for the two hydrogen peroxide crystals. The mean relative deviation (MRD) is defined as MRD
and the mean absolute relative deviation (MARD) is defined as The numbers in the table are MRD values, and the color coding follows the MARD values according to the limits given in the footnotes of this table.The reference values originate from experiment in all cases (see text).The coloring of the columns marked "Method" is an attempt to summarize the individual columns pertaining to the crystals.
a Color codes for the crystals: Cell volume: green for MARD ≤ 3%, orange for MARD > 9%, yellow in between.R(O⋅ ⋅ ⋅O): green for MARD ≤ 1%, orange for MARD > 3%, and yellow in between.E latt (0 K): green for MARD ≤ 20% and orange for MARD > 20%.b MRD and MARD calculated over both crystals.c MRD and MARD calculated for H2O2(s).Our "experimental" E latt reference value after correction for thermal and ZPE effects is 71 kJ/mol (see Sec. III B).Strictly speaking, the M in the labels MRD and MARD are not fully adequate in this column as N = 1 in the error formulas here.The cell volume and O⋅ ⋅ ⋅O distance variations as a function of functional are given in Fig. 2 (and the numerical values in Table S2).Several of the methods, including optPBE-vdW and several of the D3-corrected methods, manage to reproduce the experimental results (dashed lines in the figures) quite well, while for some of the other functionals, the O⋅ ⋅ ⋅O distances are in error by close to 0.10 Å.It is seen from the figure that for the D3-corrected functionals, the D3 contribution plays a major role in determining the optimized volume, especially for RPBE but also for PBE and HSE06, but less so for the SCAN functional (consistent with the fact that SCAN already includes some amount of intermediate range dispersion interaction).The graphs in Fig. 2(a) are similar for the two crystals.

B. Lattice energies
E latt for a molecular crystal is defined as the energy needed to be added to break the intermolecular interactions between the constituent molecules at 0 K, i.e., in our case, the reaction energies for Thus, E latt = E tot (optimized crystal) − E tot (optimized isolated molecules).Lattice energy is usually defined as a 0 K property; even so, we will use the label E latt (0 K) for clarity below.Figure 3(a) presents our calculated E latt values and the numerical values are given in Table S3.It is seen that the energy vs functional variation is quite similar for the two crystals.
The melting point of H 2 O 2 (s) is 273 K, and the sublimation enthalpy has been determined experimentally at that temperature, giving ΔH sub expt (273 K) = 64.9kJ/mol. 39We would like to use this measurement as a reference to assess our DFT-calculated lattice energies of H 2 O 2 (s), but a conversion protocol is needed to connect the two types of energy.This transformation (Scheme 1) will be discussed next.S2.S3.

The Journal of Chemical Physics
Starting from the left in Scheme 1, Eq. ( 1) 40 together with the Cp expt data were used to convert the experimental ΔH sub (273 K) to a reference value at 0 K, according to Experimental values for the heat capacities Cp expt (g, T) and Cp expt (s, T) in the relevant temperature interval were taken from Refs.41 and 42, respectively.We evaluated the integral numerically; it becomes 1.5 kJ/mol, giving ΔH sub expt (0 K) = (64.9− 1.5) = 63.4 kJ/mol.
Next, the zero-point vibrational energy difference between the gas and solid phases, ΔZPE = ZPE(g) − ZPE(s), was taken into account.The ZPE values were obtained from lattice dynamics calculations for the crystal and harmonic vibrational calculations for the gas-phase molecule.It turns out that ΔZPE varies very little between the different functionals in our study, and all our ΔZPE values come out as −8±2 kJ/mol for H 2 O 2 (s), which we will use in the following equation: We give E latt expt (0 K) the superscript expt although calculated results are used in one of the terms, and it becomes E latt expt (0 K) = 63.4 − (−8) kJ/mol = 71 kJ/mol, which is the dashed line in Fig. 3(a).
Scrutiny of Fig. 3(a) shows that several of the functionals manage to land close to our experiment-based reference value for H 2 O 2 (s), in particular, the PBE-D3, HSE06-D3, SCAN, and BLYP-D3 functionals.Even for RPBE, which is the non-dispersioncorrected functional that performs the worst with respect to the experimental reference, the D3 correction does a very good job in approaching the reference and increases the lattice energy by about 50% for both crystals.This is, in fact, similar to some of the percentages obtained by Brandenburg et al. 7 for the D3 correction to RPBE in their Ice10 dataset of ten ice polymorphs.
We end this section by predicting a room-temperature sublimation energy for H 2 O 2 ⋅2H 2 O(s), ΔH sub calc (273 K), as we have not found any measurements of the sublimation enthalpy at any temperature in the literature for this crystal.From the DFT-calculated E latt calc (0 K), we will use Scheme 1 in reverse, with Eqs. ( 1) and ( 2) taken in reverse as well.
Figure 3(a) showed that for H 2 O 2 (s) the BLYP-D3 value is spot-on the experimental reference value and we, therefore, make the pragmatic decision to use this functional for the estimate of ΔH sub calc (273 K) for the dihydrate.The BLYP-D3 E latt calc (0 K) value is 196 kJ/mol.We calculated ΔZPE for the dihydrate using a number of functionals, yielding −32 ± 3 kJ/mol of H 2 O 2 ⋅2H 2 O formula units; we will use the value −32 kJ/mol.For the integral in Eq. ( 1), we use the value 3 × 1.5 kJ/mol = 4.5 kJ/mol supported The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp by our calculations of the integral for H 2 O 2 (s) and by results concerning temperature effects on ΔH sub for similar structures given by Červinka et al. 13,43 Our resulting prediction of the ΔH sub calc (273 K) value from reversing Eq. ( 1) thus becomes 196 + (−32) + 4.5 = 169 kJ per mole of formula units or 56 kJ per mol of molecules.The value 169 kJ/mol is consistent with the results of Jenkins and Glasser 44 who compiled experimental dehydration enthalpy data for a large number of crystalline hydrates at room temperature.Their data suggest a cost of ∼55 kJ/mol for every water that is released from such a hydrate, also when the dehydration occurs sequentially.This is in good agreement with our value 56 kJ/mol.

C. Assessment table
Table II (given in Sec.II B) attempts to summarize the performance of the various DFT functionals in terms of mean relative deviations (MRDs) and mean absolute relative deviations (MARDs) from the experimental reference values for three properties: cell volume, H-bond distance, and lattice energy.The MRD is used for the numbers in Table II and is defined in the caption, and the MARD is used for the color coding, which is defined in a table footnote.We have chosen the color code limits based on our experience of what should be considered high-quality experimental structure or energy determinations and based on our assessment of experimental uncertainties (reported in Sec.III E).Admittedly, the limits are somewhat subjective, and this is also why we do not make any cells flaming red, but rather orange.
Overall, the functionals behave in a systematic manner [see, e.g., the cell volumes in Fig. 2(a) or the lattice energies in Fig. 3(a)], which is useful in discussions of trends even in cases when predictive accuracy is not achieved.Incidentally, plotting the optimized cell volumes in Fig. 2(a) against the corresponding lattice energy values in Fig. 3(a) yields a rather good correlation (not shown here).
We find that the overall best-performing ("the most green and yellow" in Table II) functionals in terms of reproducing the experimental benchmark values are optPBE-vdW and RPBE-D3, but several of the others perform well too, such as vdW-DF-cx, HSE06-D3, BLYP-D3, PBE-D3, PBE, and SCAN.We also note that for the pure PBE, RPBE, and HSE06 functionals, the D3 correction achieves marked improvements (not so for SCAN).

D. Comparison with assessment of ices in the literature
Brandenburg et al. 7 reported cell volumes and lattice energies for ten ice polymorphs using a number of functionals, five of which are in common with our study (PBE, PBE-D3, RPBE, RPBE-D3, and BLYP-D3).Their lattice energy values for the ordered ice II and IX polymorphs and an ordered form of ice Ih are compared with our H 2 O 2 (s) results in Fig. 3(b), and the variation with respect to functional is seen to be similar for all four crystals.In addition, the absolute values are reasonably close to each other, within 10 kJ/mol for a given method.This again supports the notion of significant chemical similarity between H 2 O and H 2 O 2 .
The extensive benchmarking study by Della Pia et al. 8 scrutinized the performance of 45 variants of DFT functionals for 13 ice polymorphs using diffusion Monte Carlo (DMC)-calculated lattice energies as the benchmark references.The authors found that it was difficult to find one single functional that achieved reliable performance for all ice phases and with respect to both absolute E latt energies and the relative stability order among the polymorphs.It should be noted, however, that the criterion for "reliable performance" was set very tightly in that study, partly necessitated by the fact that the DMC lattice energies for the 13 structures are very similar, all lying within 57 ± 2.5 kJ/mol.
The resulting assessment table in Ref. 8 shows that the performance varies considerably between the functionals tested, from a MAE (mean absolute error, i.e., the same as we call MAD in this paper) value of 0.9 kJ/mol (for RSCAN and revPBE-D3) to 36 kJ/mol (for LDA).For our hydrogen peroxide crystals, there is also significant E latt variations with the functional [Fig.3(a)].Many of the functional variants in Ref. 8 yield a MAD value below 4 kJ/mol ("chemical accuracy"), while some others, such as the well-known PBE-D3 and optPBE-vdW functionals, both gave MAD values of about 9 kJ/mol, which translates to a relative error (MRD) of +16%.This is similar in sign and size to our values in Table II; our calculated PBE-D3 lattice energy is 10% larger in magnitude than the reference value, and for optPBE-vdW, it is 17% larger.
Altogether we have five functionals in common with the ice study of Ref. 8. For those functionals, our MRD values for the H 2 O 2 crystal, and theirs (averaged over the thirteen ice crystals), compare as shown in the following list, where a plus sign means that the DFT-calculated energy is larger in magnitude than the reference energy.The list is as follows: optPBE-vdW (this work: +17%, Ref. 8: +9%), PBE-D3 (+10%, +9%), optB86b-vdW (+22%, +12%), and SCAN (+6%, +8%).Finally, for the PBE method, our value is −17%, while for the ices in Ref. 8 the variation is large between the different polymorphs: from +5% larger than the DMC reference for some polymorph to −20% for some other.In summary, the large deviations from the experimental reference that we find for the E latt values for many of the functionals in Table II, and which cause(d) us some concern, on the whole appear to be consistent with the ice results of Ref. 8.

E. Uncertainties in the experimental benchmark values
The accuracy of the experimental benchmark values, i.e., the dashed lines drawn in Figs.2(a) and 3 and used in Table II, will be discussed here.Clearly, the numbers and details presented here pertain only to our systems, but the kind of uncertainty analysis that we present could hopefully be of more general relevance.We start with the statistical errors related to the experimentally determined crystal structures.The experimentally reported standard deviations of the cell parameters derived from the least-squares diffraction refinements are ≤0.001Å for a, b, and c of H 2 O 2 (s), 17 and for H 2 O 2 ⋅2H 2 O(s), they are 0.002 Å for a and b, 0.01 Å for c, and 0.03 ○ for the angle. 18As a test for the less accurately determined H 2 O 2 ⋅2H 2 O(s) structure, we increased the values by two experimental standard deviations for all three experimental cell parameters, The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp and the result was a volume increase by 1.0%.Our chosen upper MARD limit for the "green class" (cf.Table II) in the volume column is 3%, which we deem reasonable also with respect to the uncertainty of the experimental cell parameter values.The same is true for R(O⋅ ⋅ ⋅O): the published standard deviations on all experimental R(O⋅ ⋅ ⋅O) distances are ≤0.005Å for the two crystals, so 2× the experimentally given standard deviations gives 0.3% of the distance, which is well below our 1% MARD limit for the green R(O⋅ ⋅ ⋅O) class.
Regarding systematic experimental errors, the difficulty in determining H positions experimentally (actually also using neutron diffraction) was mentioned earlier, and we do not use such distances as benchmarks.Furthermore, temperature effects should be considered.The structural reference data here originate from diffraction studies at 110 and 83 K, respectively.It has been shown in the literature that temperature effects on cell parameters and atomic coordinates can be appreciable for molecular crystals when going from room temperature to low-temperature measurements.However, below 100 K, say, the changes are modest, as demonstrated in a small literature survey 45 performed for diffraction-determined crystalline hydrate structures that had each been measured at three different temperatures.The survey showed that between 300 and 100 K, the differences are appreciable, but when going from 100 to 20 K or below, the O⋅ ⋅ ⋅O hydrogenbond distances decrease by less than 0.01 Å and the cell volumes by less than 1%.Another example is ice IX, where diffraction studies have shown that the cell volume decrease is only 0.3% between 110 and 30 K. 46 As long as low-temperature experimental structures are available, the temperature and zeropoint energy effects on the cell parameters and H-bond distances are much smaller than the variations with the selection of functional.
For the lattice energies (third data column in Table II), the reference value is the E latt expt (0 K) value of H 2 O 2 (s) derived via the equations in Sec.III B and based on measurements of ΔH sub (273 K).Our E latt expt (0 K) value is 71 kJ/mol with an estimated uncertainty of ±6 kJ/mol.
Below, we report how this uncertainty is derived.
We have considered three sources of uncertainties.(i) ΔH sub measurements are typically published with an uncertainty of 1 or a couple of kJ/mol; 13 let us use 2 kJ/mol.(ii) As for the integral in Eq. ( 1), its value is very small (1.5 kJ/mol), and the tables in Ref. 13 for the elevated temperature measurements suggest that also the uncertainty is small, just one or a couple of kJ/mol.Here, we set it tentatively to 2 kJ/mol.(iii) For the uncertainty of ΔZPE, we will tentatively use the spread between the different functionals (±2 kJ/mol) although this might hide some errors, e.g., the error due to the harmonic approach.Altogether, then, our reference ΔH sub expt (273 K) value may be fraught with an uncertainty of the order of 2 + 2 + 2 kJ/mol, which out of a total value of 71 kJ/mol is 8%.To be on the safe side in our assessment of the quality of the predictions by the different functionals, we therefore used rather liberal limits for the lattice energy color code: the "good and green" cells in Table II are those with MARD ≤ 10%, the most deviating ones with MARD > 20% are orange, and yellow is used for those in between.

IV. DISCUSSION OF D3 EFFECTS AND INTERMOLECULAR INTERACTIONS
In Sec.III, we have seen that the D3 London dispersion corrections overall do a very good job of correcting the underlying non-dispersion-corrected DFT volumes and lattice energies of our target systems.This section discusses "the contents" of the D3 correction.
A. Is the D3 correction energy for our crystals mostly London dispersion?
London dispersion interactions are fundamentally a quantummechanical phenomenon related to the correlated motion between dynamically fluctuating dipolar electron densities of (neutral) atoms or molecules.They give rise to a 1/r 6 "long-range" distance behavior (see, e.g., Stone 47 and references therein for a comprehensive review).
University textbooks in chemistry teach us that oxygen atoms, which are more valence-electron-rich than H, display stronger London dispersion interaction than H, such that the O-O dispersion interaction will be stronger than O-H and H-H.However, contrary to such expectations, the O-O, O-H, and H-H  Another indication of non-dispersion effects being handled by the D3 correction is the well known large variation of the D3 contribution with choice of DFT functional, both regarding energy magnitude and regarding the "contraction force" that it invokes.As the crystals' electron densities with the different GGA functionals are quite similar, the variation with functional suggests that the D3 correction not only adds London dispersion but also manages to correct for various limitations inherent in the underlying functionals.This is highlighted in Fig. 5, where we performed an isotropic volume scan for the H 2 O 2 (s) crystal, contracting and expanding the crystal cell rigidly around the optimized volume while keeping the internal molecular geometries fixed.The D3 correction energy was calculated directly from the stand-alone DFT-D3 program 4,48 for different functionals.The PBE functional, for example, gives the blue curve (third from the top) in Fig. 5(a), and this curve is equivalent to the blue ΔE D3 curve in the right figure, which is a "pedagogical illustration" of how the D3 correction shifts the energy minimum at the PBE(-D3) level.For the D3-corrected functionals used in Fig. 2(a), the slopes in Fig. 5(a) are seen to vary in the following order from smallest to largest magnitudes: in good agreement with the cell volume contractions caused by the D3 correction after full optimizations as displayed in Fig. 2(a) and Table S2: At the same time, there are (of course) ample reasons to believe that London dispersion is still appreciable in our two peroxide crystals as well as in the ices, and that dispersion makes a significant contribution to the D3 correction.(One such straightforward indication is that calculations of the "intermolecular electron correlation energy" from high-level quantum-chemical calculations for small water clusters compared to converged basis-set Hartree-Fock binding energies yield appreciable values for the dispersion energy already for the water dimer; cf., e.g., the results in Refs.50 and 51.It should be noted, however, that the "intermolecular electron correlation energy" is not solely dispersion energy but also includes the effect of changes of the static molecular dipole moment and the polarizability when the Hartree-Fock Hamiltonian is changed into a correlated Hamiltonian.).
Figures 6(a In summary, we believe that the large D3 contributions to E latt for all five systems in Figs.6(c) and 6(d) are the result of both the many close molecule-molecule neighbor interactions and their accompanying London ("fluctuating dipolar") dispersion and the fact that the D3 energies correct for various other shortcomings of the underlying DFT functionals.The effect of the increased polarization induced by the D3-invoked cell contraction will be discussed next.

B. Increased polarization as a consequence of D3 contraction
There exist many different schemes to classify intermolecular interactions, e.g., in terms of energy components such as electrostatic, induction (polarization), dispersion, and exchange, and variants of these.The dispersion contribution to the D3 correction was discussed in Sec.IV A. We will not discuss exchange, which accompanies any electron density redistribution invoked by other attractive or repulsive energy components, but we will examine the electrostatic and polarization interactions in the context of the D3 correction, namely to what extent the addition of the D3 correction indirectly alters the molecular polarization by means of the cell contraction.The answer will be "modest," and we use Fig. 7 to justify it.
The experimental permanent electric dipole moments (μ total ) of the gas-phase water and peroxide molecules are both appreciable and they are quite similar   S2).The μ local value a targeted OH group was calculated from the electron distributions within the two AIM basins pertaining to the OH group and then it was projected along the O-H bond.μ local is free of the angle dependence plaguing the total dipole moment for H 2 O 2 and is thus a more appropriate dipole moment descriptor to probe of the external perturbation caused by intermolecular interactions.From Fig. 7(b) we note that the polarization effect is substantial both in the bimolecular complexes and in the crystals, and varies according to μ local, H2O2 (crystals) ≫ μ local, H2O2 (bimolecular complexes) ≫ μ local, H2O2 [H 2 O 2 (g)].Also the water molecules in the dihydrate are appreciably polarized, by ∼35% compared to the gasphase value.All calculated values in Fig. 7, plus the water dipole moments in the different systems, are given in Tables S4 and S5.
The in-crystal enhancement of the molecular dipole moments will enhance the electrostatic part of the lattice energies, but the question we want to discuss is, more specifically, if the D3-induced cell contraction is expected to enhance the lattice energy appreciably.The quantity on the x-axis in Fig. 7(b) is the optimized R(O⋅ ⋅ ⋅O) H-bond distance.A clear "μ vs R(O⋅ ⋅ ⋅O)" correlation is seen, not only for different systems calculated with one method but also the five different functionals approximately fall on the same curve (piece-wise lines), and in the same order.This is a consequence of (i) the resemblance between electron densities calculated with different DFT functionals and (ii) the dominance of the H-bonds in steering the polarization effect.The D3-induced foreshortenings of the R(O⋅ ⋅ ⋅O) distances in H 2 O 2 (s) are 0.02-0.06Å depending on the functional (Table S2) and 0.03-0.14Å for the dihydrate.According to the Δμ/ΔR O⋅ ⋅ ⋅O slopes in Fig. 7(b), such contractions enhance the in-crystal dipole moments by some 5%-15%, which will result in a significant, yet modest extra contribution to the lattice energies.We conclude that the large increase of the lattice energies when adding D3 for some of the functionals originates to some extent also from increased polarization effects.

V. WHICH MOLECULE FORMS THE STRONGEST H-BONDS: H 2 O OR H 2 O 2 ?
We end by returning to the question raised in the Introduction concerning which of the two molecules has the best H-bonding capacity.This is not entirely simple to answer although it has often been stated in the literature that the H 2 O molecule is a better H-bond acceptor than H 2 O 2 and that H 2 O 2 is a better H-bond donor than H 2 O, with or without convincing argumentation.We performed calculations for the H 2 O⋅ ⋅ ⋅H 2 O 2 and H 2 O⋅ ⋅ ⋅H 2 O bimolecular complexes at the CCSD(T) level for geometries optimized at the MP2 level, in both cases with a close to converged basis-set.The interaction energy E tot (complex) − E tot (H 2 O 2 ) − E tot (H 2 O) was calculated, with all terms corresponding to optimized isolated gas-phase species.The H-bond donor is H 2 O in both complexes and the acceptor is either H 2 O 2 or H 2 O.It turns out that in both cases only one Hbond forms, which simplifies our possibility to interpret the interaction energy as a hydrogen-bond energy.The H 2 O⋅ ⋅ ⋅H 2 O makes the strongest complex, which is a strong indicator that water is the better acceptor.As for the H-bond donating ability, we refer to the results from an infrared spectroscopic matrix isolation study performed by Goebel et al. 56 In that study, measurements were performed for HOH⋅ ⋅ ⋅O(CH 3 ) 2 and HOOH⋅ ⋅ ⋅O(CH 3 ) 2 gas-phase complexes.As these complexes are structurally similar, the larger OH vibrational frequency downshift observed for the second complex constitutes a strong indication that H 2 O 2 is a better H-bond donor than H 2 O. Vener et al. 57 calculated vibrational frequencies to reach a similar conclusion for crystals that contained either water⋅ ⋅ ⋅amino acid or H 2 O 2 ⋅ ⋅ ⋅amino acid hydrogen bonds, i.e., again a case of the same acceptor but different donors.They derived "H-bond energies" based on energy-frequency correlation curves from the literature and found the H 2 O 2 ⋅ ⋅ ⋅amino acid energy to be larger.Finally, it can also be mentioned that the H 2 O 2 ⋅2H 2 O(s) crystal in our study contains H 2 O 2 ⋅ ⋅ ⋅H 2 O, H 2 O⋅ ⋅ ⋅H 2 O, and H 2 O⋅ ⋅ ⋅H 2 O 2 hydrogen bonds.Their respective R(H⋅ ⋅ ⋅O) hydrogen-bond distances increase in that order, which at least is consistent with the rule just established.
This order of the binding strength is also supported by radial distribution functions for hydrogen peroxide in aqueous solution generated in Refs.58 and 59 by QM/MM-MD and QM/MM-Monte Carlo simulations, respectively.However, given the small QM clusters often used in such simulations combined with the sensitivity to the MM model quality and the subtle differences between H 2 O 2 and H 2 O hydrogen bonding, some caution should likely be exercised in the interpretation of such simulation results.

VI. CONCLUDING REMARKS
We are aware of only a small number of previous periodic DFT studies of the H 2 O 2 (s) and H 2 O 2 ⋅2H 2 O(s) crystals (see Table I).Here we have focused on these crystalline systems, which contain two-in many ways kindred-molecules.We have presented computational results that highlight similarities and differences between H 2 O 2 (s) and the dihydrate, and we also made selected comparisons with ice polymorphs.The main results are as follows.
(1) For our crystals, the overall best-performing functionals in terms of reproducing the experimental benchmark values are optPBE-vdW and RPBE-D3 (the two best ones), as well as vdW-DF-cx, HSE06-D3, BLYP-D3, PBE-D3, and SCAN.(2) It is the many non-H-bonded neighbors out to, say, 5 Å (Fig. 6) that make the D3 corrections to the lattice energies very sizable.For H 2 O 2 (s) and PBE, the D3 correction constitutes about 25% of the lattice energy, much more for RPBE, and much less for SCAN (Fig. 3).The D3 contribution related to the H-bonded interactions is minor as only the nearest-neighbor O⋅ ⋅ ⋅O and H⋅ ⋅ ⋅O contacts are involved.(3) The large contributions of the D3 correction to the lattice energies is likely a consequence of several effects: correction for missing dispersion interaction, the ability of D3 to capture and correct various other kinds of limitations built into the underlying DFT functionals, and finally some degree of contraction-induced polarization enhancement.On the whole, the gas → crystal molecular polarization is large (Fig. 7) but the extra contribution from D3 to the polarization (via the structure contraction) is modest.The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp(5) For the functionals that we have in common with the extensive ice polymorph studies of Brandenburg et al. 7 and Della Pia et al., 8 we find that their assessments, and ours for the two H 2 O 2 -containing crystals, largely agree, both in magniand trends with respect to the references chosen.For technical/computational reasons, this is gratifying, and it also highlights the kinship between ices and our crystals.

SUPPLEMENTARY MATERIAL
The supplementary material contains the following: 1 2 1 2 and H 2 O 2 ⋅2H 2 O(s) in C2/c.The standard crystallographic unit cell contains four formula units for H 2 O 2 (s) and four for H 2 O 2 ⋅2H 2 O(s).The structures are shown in Fig. 1.All H 2 O 2 molecules in H 2 O 2 (s) are structurally equivalent.The H 2 O 2 molecules in H 2 O 2 ⋅2H 2 O(s) are also structurally equivalent as are the H 2 O molecules.The structures were fully optimized with a range of functionals, as detailed below.

FIG. 1 .
FIG. 1. Crystal structures of the crystallographic cells and local coordinations around the constituent molecules in (a) H 2 O 2 (s) and (b) H 2 O 2 ⋅2H 2 O (s).The structures shown were optimized with the BLYP-D3 method and are very similar to the experimental structures (see the supplementary material for tables of structural parameters with different methods).The blue lines mark hydrogen bonds in the structures.The unit cell content of H 2 O 2 ⋅2H 2 O(s) displayed is not complete in order to better highlight the bonding situation.The H 2 O 2 molecule is symmetric with a two-fold axis in both crystals.
ARTICLE pubs.aip.org/aip/jcpIII.RESULTS: DFT TRENDS AND ASSESSMENT OF THE FUNCTIONALS A. Optimization of the crystal structures The 110 K neutron diffraction results from Ref. 17 is our reference structure for H 2 O 2 (s), as mentioned.For H 2 O 2 ⋅2H 2 O(s), the x-ray diffraction study by Olovsson and Templeton at 83 K 18 was used as the reference.It is of an adequate quality, but this does not change the fact that x-ray diffraction refinements lead to severely foreshortened covalent O-H bond distances and a corresponding elongation of the H⋅ ⋅ ⋅O hydrogen-bond distance.Therefore, we will not use the experimental H⋅ ⋅ ⋅O distances as benchmark values for H 2 O 2 ⋅2H 2 O(s) but only the O⋅ ⋅ ⋅O distances.

FIG. 3 .
FIG. 3. Resulting optimized lattice energies, E latt =−[E(optimized crystal) − ∑ E(optimized molecules)], as a function of computational method.Positive values mean stronger binding.(a) E latt for the DFT-optimized crystals are given; the experimental reference value for H 2 O 2 (s) is marked by a dashed line (71 kJ/mol).(b) This figure compares our E latt results for H 2 O 2 (s) with literature results for three ice polymorphs, namely ice Ih, ice II, and ice IX from Ref. 7, using five functionals.The numerical values are listed in TableS3.

FIG. 4 .
FIG. 4. The pie-diagrams display the relative D3 energy contributions from OO, OH, and HH pairs in the 0-5 Å range (upper row) and the relative number of these atom-atom distances (lower row).RPBE-(D3) values were used (see text for details).The ice polymorph results included in this figure originate from structural optimizations performed by us with similar technical settings as for the peroxide crystals.
Now altogether Figs. 4 and 5(a) have provided reasons to conclude that the D3 term includes other effects than merely London dispersion (see also the discussion by Shahbaz and Szalewicz 49 ).

FIG. 5 .FIG. 6 .
FIG. 5. Results from volume scans performed for the H 2 O 2 (s) crystal, where all the cell parameters were expanded and contracted uniformly (in percent) from the optimized PBE structure, keeping the internal geometry of the molecules fixed, and the total and D3 energies were calculated using different DFT-D3 functionals.(a) Using the D3 program of Grimme 4,48 the D3 energy contribution was recorded along the scan for different functionals.The slopes suggest the relative magnitudes of the structural shift that the D3 correction will induce.(b) Illustration of the "perturbation" caused by the D3 energy on the structure and energy of the H 2 O 2 (s) crystal.The potential energy curves from the PBE and PBE-D3 scans illustrate how the D3 energy (i.e., the difference ΔE D3 between the two curves) can be interpreted as the perturbation that shifts the optimized cell minimum.In both (a) and (b), the energies are given per crystallographic unit cell, i.e., per four H 2 O 2 molecules.

2 FIG. 7 .
FIG. 7. "In situ" dipole moments of the H 2 O 2 molecule in different systems, calculated using different functionals and the Bader formalism (see the text for details).The gas-phase H 2 O molecule and H 2 O 2 molecule dipole moments are also presented.(a) Total molecular dipole moments of the optimized isolated H 2 O and H 2 O 2 molecules.The experimental values are marked as dashed lines, and references are given in the text.(b) The local OH dipole moment of H 2 O 2 , μ local , calculated over the O and H atomic AIM basins and projected onto the intramolecular OH bond axis.The strip at the bottom refers to μ local for the isolated H 2 O 2 molecule with the same five DFT variants.

Figure S1 :
Comparison of optimized structures and resulting lattice energies using the D3(0) and the D3(BJ) damping schemes for H 2 O 2 (s) and H 2 O 2 ⋅2H 2 O(s).

TABLE I .
Published computational studies of condensed-matter systems with H 2 O 2 .Both quantum-mechanical and force-field publications are included.Solids and liquids are included but not ad-molecules on surfaces.

ARTICLE pubs.aip.org/aip/jcp data
for that compound.Neither of the two compounds appear to suffer from the partial disorder that is present to various degrees in many of the ice polymorphs.
molecule because the variation of the total molecular dipole moment when the molecule is bound is largely dominated by the variation of the molecule's dihedral angle, which is a soft angle (the angle values are given in Table

Table S1 :
Optimized geometries from VASP calculations for H 2 O 2 (g) and H 2 O(g), including some results from the literature.Table S2: Distances, angles, and cell parameters for the optimized structures, and experimental references for H 2 O 2 (s) and H 2 O 2 ⋅2H 2 O(s).Table S3: Calculated lattice energies, zero-point vibrational energies (ZPE, harmonic approximation), and sublimation enthalpies for H 2 O 2 (s) and H 2 O 2 ⋅2H 2 O(s).Table S4: Calculated Bader dipole moments at optimized equilibrium geometries.Table S5: Bader dipole moments for H 2 O and H 2 O 2 in various surroundings.