The Ionic Product of Water in the Eye of the Quantum Cluster Equilibrium.

The theoretical description of water properties continues to be a challenge. Using quantum cluster equilibrium (QCE) theory, we combine state-of-the-art quantum chemistry and statistical thermodynamic methods with the almost historical Clausius-Clapeyron relation to study water self-dissociation and the thermodynamics of vaporization. We pay particular attention to the treatment of internal rotations and their impact on the investigated properties by employing the modified rigid-rotor-harmonic-oscillator (mRRHO) approach. We also study a novel QCE parameter-optimization procedure. Both the ionic product and the vaporization enthalpy yield an astonishing agreement with experimental reference data. A significant influence of the mRRHO approach is observed for cluster populations and, consequently, for the ionic product. Thermodynamic properties are less affected by the treatment of these low-frequency modes.

1. Introduction

Since the dawn of electrochemistry at the end of the 18th century, scientists have sought to understand the decomposition of water due to an electric current. In 1805, Theodor Grotthuss developed the first theory to describe this phenomenon [1,2,3], in which he assumed that water molecules become polarized under the influence of electric voltages of electrodes, that is, oxygen atoms become negative and hydrogen atoms become positive. The attraction and repulsion of the electrodes consequently lead to an alignment of the molecules, creating water wires that span the electric cell. Grotthuss further assumed that if the electromotive force is strong enough, the outer hydrogen and oxygen atoms of the wire, those in contact with the electrodes, will separate from their corresponding molecules, forming molecular oxygen and hydrogen at the respective electrode. The residual atoms subsequently recombine along the water wire, organizing as a new set of molecules. About 50 years later, Rudolf Clausius recognized that water molecules can also dissociate spontaneously without an applied electric field; the applied voltage just regulates and accelerates the movement of charged particles toward the corresponding electrode [4]. This self-ionization, measured in terms of pH and characterized by the ionic product, (or ), is now a cornerstone of our understanding of aqueous acid–base chemistry.

Grotthuss’s breakthrough ideas of polarization as well as continuous dissociation and reformation of water molecules were embraced by Danneel in 1905 to explain the unusually high mobility of protons in liquid water [5], a process which is nowadays referred to as the Grotthuss diffusion mechanism and has far-reaching consequences in biological systems (see Ref. [6] for a comprehensive special issue on the subject).

In technical applications, on the other hand, pH is an important stimulus. The most basic examples are probably pH indicators, whose electronic properties—and thereby their color—are altered by proton activity. More advanced examples include pH-responsive polymers, which have far-reaching applications in catalysis [7], nanomedicine [8], materials science [9], and many other fields [10,11].

Today, more than 200 years after the discoveries of Grotthuss and Clausius, theoretical chemists and physicists have an arsenal of advanced techniques to study proton transfer and transport in liquid water at high levels of theory through static calculations and dynamic simulations [12,13,14,15,16,17,18,19]. Yet, the actual computation of the ionic product is still a cumbersome task. The major challenges lie in the corresponding ionic concentrations being too small to allow direct simulations as well as water being extremely efficient in catalyzing its own dissociation, complicating static calculations. The latter can be demonstrated beautifully by comparing at standard state conditions against the value obtained for the hypothetical gas-phase dissociation of a water molecule, , at the B3LYP/6-311++G** level of theory [20]. In ab initio molecular dynamics simulations, the challenges are typically overcome through the application of artificial biases or constraints and vast amounts of computational power [21,22]. Static quantum-chemical calculations often resort to combinations with continuum theories and classical models [23,24,25,26,27,28], tainting the first principles character of the calculation.

In 2022, we remember Grotthuss’s 200th deathday and the bicentenary of Rudolf Clausius’s birth. In 2021, we also celebrated the 80th birthday of Frank Weinhold, to whom this article is dedicated. It turns out that Weinhold’s quantum cluster equilibrium (QCE) theory [29,30] establishes two important links to the works of Grotthuss and Clausius. QCE is a theory to describe associated fluids and their mixtures through representative cluster structures. The cluster weights are related to the chemical potential of the clusters, taking vibrational and rotational entropic contributions through application of the rigid-rotor–harmonic-oscillator (RRHO) approximation into account. A hallmark of QCE theory is the ability to treat phase transitions, such as in Weinhold’s seminal work, where he employed the Clausius–Clapeyron relation to describe the thermodynamics of water vaporization [30]. It turns out that QCE also provides alternative means to studying aqueous dissociation processes, void of the complications mentioned above. In the past, QCE has thus been successfully applied to compute the ionic product of water and dissociation constants of other weak acids as a function of temperature [20,31,32].

In this contribution, we build upon these foundations. Inspired by the historical background and driven by the outstanding performance of the QCE approach in describing the ionic product of water, we reconsider both the ionic product and the Clausius–Clapeyron relation for vaporization, using the so-called modified rigid-rotor–harmonic-oscillator (mRRHO) [33] model. Therein, those vibrational modes of the cluster structures, that actually represent internal rotations of molecular units and that occur at small wavenumbers, are appropriately treated by the partition function of a hindered rotor, rather than a harmonic oscillator.

The article is structured as follows. In Section 2, we revisit essential elements of QCE theory and the computational protocol, followed by the presentation and discussion of the results in Section 3. Conclusions are provided in Section 4.

2. Methods

2.1. QCE Theory

Quantum cluster equilibrium theory treats the coupled chemical equilibria between a monomer structure, , and a set of clusters, , built therefrom, according to where ℘ is a running index labelling the clusters and denotes the cluster size (number of monomers). All clusters are characterized quantum-chemically: their geometries are optimized, their ground-state energy is determined, and a vibrational frequency analysis is performed. Furthermore, cluster volumes, , are computed. QCE theory has been fully described several times in the past [29,34,35,36,37]. Here, we only review basic aspects of the theory and recent modifications [38] to the partition functions.

Assuming independent clusters with fixed populations, , and independent degrees of electronic, translational, vibrational, and rotational freedom, the canonical partition function, , at constant volume, V, and temperature, T, is given by where the factor accounts for the indistinguishability of like clusters and the terms denoted by are individual cluster partition functions corresponding to the respective degree of freedom (dof).

For the electronic partition function, only the ground-state electronic energy and a mean-field term, describing average intercluster interactions, are considered: where is called mean-field parameter and is the particle density of the system. The remaining individual partition functions are given by textbook [39] expressions for the following simple models: particle in a box (translation), harmonic oscillator (vibrations), and rigid rotor (rotations). The combination of the last two is typically referred to as rigid-rotor–harmonic-oscillator (RRHO) approximation. The accessible volume that enters into the translational partition function is reduced by an excluded volume term taking the finite size of each cluster into account: where is called excluded volume scaling parameter and is the cluster volume.

Although the theory is developed in the canonical ensemble, neither populations nor volume are input parameters in a QCE calculation. Instead, the volume is chosen to be consistent with a specified, external pressure, where is Boltzmann’s constant. The populations are chosen to satisfy the condition for chemical equilibrium, where are the chemical potentials of the clusters. These two constraints lead to a set of coupled polynomial equations, which can be solved by the Peacemaker QCE code [34,35,37].

As indicated above, QCE theory contains two adjustable parameters, and , which are typically adjusted to fix the density at ambient conditions and the temperature of phase transition, although other options are conceivable. The special case is called QCE(0) and corresponds to the description of a cluster gas.

Instead of scanning all possible combinations of and on a coarse grid as in earlier studies, the current version of Peacemaker [34,35,37] features an implementation of the downhill simplex algorithm [40] which speeds up the calculations tremendously [41] while also increasing the accuracy of the parameter optimization. Here, and were optimized within a range of 0.0 to   and 0.5 to 2.0, respectively, so as to minimize the deviations of the computed density at room temperature and of the phase transition temperature from their respective experimental values. Once the parameters are optimized, the accessible temperature and pressure ranges of the QCE approach are, in principle, limited only by the condition that Boltzmann statistics apply, i.e., the number of particles N is much smaller than the number of accessible energy states . This condition is typically fulfilled by molecular solvents at room temperature, but imposes a system-dependent lower boundary to the accessible temperature range.

Note that the computational bottleneck of a typical QCE study lies in the quantum chemical optimization of the clusters, whereas the QCE calculation itself will finish within seconds to minutes. This allows, in principle, for the inclusion of large-scale cluster sets containing hundreds of clusters [38] if the underlying quantum chemical method makes it feasible. Here, we opted for higher level quantum-chemical methods at the cost of using a smaller but representative cluster set.

2.2. Modified Partition Functions

In the following, we briefly review the modified rigid-rotor–harmonic-oscillator (mRRHO) QCE approach, which we have recently implemented [38] in the Peacemaker [34,35,37] QCE code. In this approach, the vibrational partition function is modified in accordance to the model presented in Ref. [33], where the product is over all normal modes i with frequencies , is the traditional partition function of a harmonic oscillator, and is the partition function of a hindered rotor. The latter is given by where is the average moment of inertia of the cluster and is the moment of inertia corresponding to the normal mode. The modified partition function corresponds to that of a regular harmonic oscillator for high-frequency modes and to a hindered rotor for low-frequency modes. The switching between the two models is handled by a switching function, , for which the Chai–Head-Gordon damping function is applied [42],

Therein, is a user-defined rotor cutoff value. In this work, we applied two distinct cutoffs: 50 and 100 cm−1. The former is a sensitive choice in the recommended range of reasonable values between 20 and 100 cm−1 and the latter represents an upper limit. Choosing both allows us to investigate the influence of this parameter in a sensible way [38].

2.3. Systems Investigated

A set of clusters representative of the liquid phase of water was constructed. The set contains both regular and zwitterionic clusters and has been previously employed to compute the ionic product of water [20].

Figure 1 depicts the regular water clusters, with oxygen atoms colored blue and hydrogen atoms colored white. Their sizes range from the water monomer, W1, to a decamer cluster, W10. The trimers, W3c and W3u, as well as the clusters W5c and W6c exhibit exclusively two-fold hydrogen bond motifs, i.e., each water molecule accepts and donates one hydrogen bond. Three-fold hydrogen bond coordination can be found in clusters W5p, W8c, W8p, and W10. A four-fold hydrogen bond pattern is present in the clusters W7, W8b, and W9.

Figure 1

Ball-and-stick representation of the neutral water clusters. Oxygen: blue; hydrogen: white.

Figure 2 shows the clusters that contain separated ion pairs, i.e., net neutral clusters with zwitterionic character. Almost all of them show the motif of the Eigen-ion (H3O+ donating three hydrogen bonds), except for the cluster W10ip2, which contains the Zundel ion motif (one H+ shared between two water molecules) [16]. The ion pair clusters range in size from a pentamer, W5ip, to decamer clusters, W10ip and W10ip2. Please note that the hydroxide anion always accepts three hydrogen bonds.

Figure 2

Ball-and-stick representation of the ion pair clusters. Oxygen: blue; hydrogen: white; red shows the hydroxide anions and green shows the hydronium cations.

2.4. Quantum-Chemical Calculations

Geometry optimizations and vibrational frequency analyses were consistently performed using density functional theory with either the hybrid B3LYP [43,44] or PBE0 [45] density functional approximations and the triple-zeta def2-TZVP [46] basis set. To account for London dispersion interactions, the D3 dispersion correction with the Becke–Johnson damping function was applied [47]. Interaction energies of cluster structures are sensitive to the basis set superposition error (BSSE), which was corrected for by using a geometrical counterpoise correction (gCP) scheme [48]. This is essential, as proton transfer during cluster formation significantly alters the monomer geometries of a cluster, invalidating the assumptions of supramolecular approaches such as conventional counterpoise correction (CP) schemes [48]. As a third computational model, the composite density functional PBEh-3c was used, which includes both the gCP and D3 dispersion correction combined with a modified PBE hybrid functional and modified def2-SV(P) basis set, both optimized to maximize error cancellation [49].

Extensive studies on the importance of the underlying quantum chemical method to the quality of the QCE approach have been conducted in the past [34,50,51]. It was found that methods based on density functional theory are capable of producing results of similar quality to methods based on wave function theory. Recently, semi-empirical methods were shown to be a viable option if large-scale cluster sets are employed [52,53]. However, while some properties may be computed with reasonable accuracy for most investigated quantum chemical methods, others, such as the ionic product of water, were shown to be sensitive to the exact choice of method [20]. More details on the computational protocol can be found in our previous publication [20].

3. Results and Discussion

3.1. Variation in the QCE Parameters

In Table 1, optimized values of the QCE parameters, and , are given for different density functional approximations, using standard QCE as well as mRRHO-QCE with switching function cutoff values of 50 and 100 cm−1.

Table 1

Optimized QCE parameters in and (dimensionless) for different methods.

Method	Standard QCE		mRRHO50		mRRHO100
	a_mf	b_xv	a_mf	b_xv	a_mf	b_xv
B3LYP/D3/gCP	0.201	1.50	0.213	1.50	0.216	1.50
PBE0/D3/gCP	0.192	1.52	0.205	1.50	0.207	1.50
PBEh-3c	0.173	1.50	0.185	1.52	0.187	1.52

There is little variation in across all methods, with values consistently lying in between 1.50 and 1.52. If all other parameters were kept the same, such a shift in would lead to a change in the computed density and temperature of phase transition of about g cm−1 and  , respectively. This very small dependence of on the method is expected, as this parameter corrects the translational volume for the eigenvolume of the clusters and thus depends neither on the quantum-chemical description of the clusters nor on the precise details of the vibrational partition function.

The mean-field parameter, , in contrast, varies between 0.17 and 0.22 J m3 mol−2. Again, if all other parameters were kept the same, such a shift in would lead to an increase in the computed temperature of phase transition of about  , whereas the computed density remains almost the same, increasing only by 0.01 g cm−1. Its value depends significantly on the functional, increasing in the order PBEh-3c < PBE0/D3/gCP < B3LYP/D3/gCP. Since accounts for interactions between clusters, the parameter has in the past been observed to compensate for underbinding of the electronic structure method, in which case values become larger. Increasing the cutoffs for the mRRHO model also leads to larger values of , indicating an indirect impact of the treatment of these internal rotations on the mean-field energy. The vibrational motions at low frequencies are usually internal rotations or translations, and hence can be related to the process of cluster association. It is therefore reasonable that the mean-field parameter is sensitive to the particular form of the vibrational partition function for these modes, while is not.

3.2. Ionic Product Dependence on mRRHO

Given the populations of the ion pair clusters, it is easily possible to compute the ionic product [20], with and referring to the proton and hydroxide ion concentrations, respectively, and exploiting the fact that there is exactly one ion pair in every cluster investigated in this study.

In Figure 3, we compare the temperature dependence of between standard QCE and mRRHO-QCE for all three investigated methods (top: B3LYP, middle: PBE0, bottom: PBEh-3c, all with D3 and gCP). In general, the application of mRRHO leads to decreasing values of , which translates to an increasing ionic product, and hence, larger ion concentrations in the liquid phase. Consequently, we see an improvement in the results for B3LYP, which predicted slightly too large values for in the standard QCE model. PBE0, which produced the best data in the standard QCE approach, is now underestimating . The same is observed for PBEh-3c, which already resulted in the smallest values within the standard RRHO model.

Figure 3

Temperature dependence of the negative logarithm of the ionic product, , for all selected methods as well as experimental [54] values. The standard method is represented by solid, mRRHO50 by dashed, and mRRHO100 by dotted lines.

The temperature dependence of is also altered by the introduction of the mRRHO approach. The observed reduction in from standard QCE (solid curves) to mRRHO50 (dashed curves) is largest at high temperatures and slightly smaller at low temperatures, corresponding to a smaller slope of the curves. Changing the mRRHO cutoff value to 100 cm−1 leads to similar, but less pronounced changes, except for PBE0, where the slope of is sightly increased again.

In Table 2, we list the number of low-frequency modes per cluster that are subject to the modified partition function for the B3LYP/D3/gCP method. The full list of vibrational modes for all methods is given in Table S1 in the Supporting Information. Only 27 vibrations are affected by the mRRHO50 approach, compared to as many as 86 vibrations in the mRRHO100 model. Still, the impact on the ionic product is more pronounced in the former case, indicating that is affected most by low-frequency modes.

Table 2

Number of vibrational modes below 50 cm−1 and between 50 cm−1 and 100 cm−1, respectively, for selected clusters from Figure 1 and Figure 2 for the B3LYP/D3/gCP method.

Cluster	# Modes ν˜<50cm^−1	# Modes 50cm^−1<ν˜<100cm^−1
W5p	1	2
W5c	2	2
W6c	2	4
W7	3	4
W8p	4	5
W8b	4	4
W8c	−	5
W9	6	5
W10	−	7
W5ip	−	1
W8ip	3	5
W8cip	−	3
W10ip1	2	6
W10ip2	−	6

To substantiate our observations for the ionic product, we consider the populations of ion pair clusters in Figure 4 and regular clusters in Figure 5, using once again B3LYP/D3/gCP as a representative example (see Figures S1 and S2 for results obtained with mRRHO100). All ion pair clusters but W8ip show increased populations if treated with the mRRHO50 approach. The population of W8ip is significantly reduced. It is interesting to note that W8ip shows three modes below 50 cm−1 while the others exhibit none or only two (see Table 2). Keeping in mind that the y-axis has a logarithmic scale, it is clear that W10ip2 is the dominating ion pair cluster, causing the overall increase in .

Figure 4

Temperature dependence of the monomer-normalized populations for the ion pair clusters, showing B3LYP/D3/gCP data as a representative example. Solid lines: conventional QCE, dashed lines: mRRHO50.

Figure 5

Temperature dependence of the monomer-normalized populations for the regular clusters, showing B3LYP/D3/gCP data as a representative example. Solid lines: conventional QCE, dashed lines: mRRHO50. Note the semilogarithmic scale.

While the ion pair clusters are essential for the calculation of the ionic product, they hardly contribute to macroscopic thermodynamic quantities of liquid water, which are determined by the regular clusters, instead. In Figure 5, we show the populations of those clusters, grouped by the structural motifs (monomer, dimer, cyclic clusters, spiro clusters, cubelike, sandwich-shaped). Since the monomer and dimer do not exhibit any vibrational modes lower than 100 cm−1, the populations of these clusters are not influenced by the modified vibrational partition function. The populations of clusters that are dominated by ring structures, be it cyclic clusters or condensed rings in spiroclusters, are reduced if the mRRHO50 model is employed.

Considering the number of low-frequency modes (Table 2), W6c, W7, and W9 contain (many) modes below 50 cm−1 and their individual population is reduced. However, W5c and W8b also show these modes and obtain slightly higher or unaltered populations. For the sandwich-type decamer, the population is significantly increased if the hindered rotor model is used to treat internal rotations. Please note that this cluster shows no modes below 50 cm−1, see Table 2. While it is still not the dominant cluster in the liquid phase, its population is enhanced by a factor of around four so that within the mRRHO50 model up to 20% of all water molecules are arranged in this shape. Interestingly, the population of the cubic W8c cluster is decreased at low temperatures and increased at higher temperatures.

3.3. Clausius–Clapeyron Analysis

Motivated by the capability of the QCE approach to treat phase transitions and the option to apply different pressures, we performed a Clausius–Clapeyron analysis following the original work of Weinhold [30]. QCE calculations with the parameters listed in Table 1 have been conducted at different pressures ranging from 0.095898 bar to 4.7572 bar and the observed temperatures of phase transition have been noted, see Table S2 of the Supporting Information.

Overall, we observe that the temperatures of phase transition deviate only slightly between standard QCE and either of the mRRHO approaches. The maximum deviation between the approaches is 3 in the investigated pressure range. A similar observation is made for the different density functional approximations, with phase transition temperature variations of up to only 5 . Of course, since the parameters have been optimized at the standard pressure of bar, deviations at very small and very large pressures are larger than at the reference pressure.

Once again, we show B3LYP/D3/gCP data as an example (Figure 6). The agreement between calculated and experimental data is astonishing, especially recalling the finite cluster set and the fact that only two parameters are optimized to reproduce one of these data points. There is some deviation in the high-temperature/high-pressure regime, which can, however, still be considered to be reasonable agreement.

Figure 6

Pressure–temperature graph for B3LYP/D3/gCP (green line) and compared to experimental [55] data (red line).

The integrated form of the Clausius–Clapeyron relation states that the natural logarithm of the pressure is proportional to the inverse temperature of the phase transition, which is why we show vs. in Figure 7. In this plot, rather small differences in the linear graph in Figure 6 at small temperatures and pressures are amplified, now occurring in the bottom right of the plot.

Figure 7

Clausius–Clapeyron plot: Logarithmic plot of pressure P versus inverse temperature T for B3LYP/D3/gCP (green line) and experiment [55] (red line).

Naturally, and due to the parameter-optimization procedure, a perfect agreement at ambient pressure () is observed, while the slight deviations from the experiment have a comparable magnitude below and above that point. With the Clausius–Clapeyron relation, an estimate of the enthalpy of vaporization can be deduced from the slope of these curves. The resulting calculated enthalpies are listed and compared to the experimental value in Table 3. As already visible from Figure 7, the agreement of the calculated values with the experimental reference is excellent. While all calculations overestimate the enthalpy of vaporization, a deviation of only 2 kJ mol−1 is observed for PBE0/D3/gCP. All mRRHO calculations lead to an overestimation of this quantity, whereas the B3LYP data are least sensitive to the mRRHO approach and yield the best overall data, which is in agreement with the observations made for the ionic product. With the exception of PBEh-3c, all methods reproduce the experimental enthalpy of vaporization to within 4 kJ mol−1, a threshold often referred to as a targeted chemical accuracy.

Table 3

First line: enthalpy of vaporization in kJ mol−1 calculated through a Clausius–Clapeyron analysis. The experimental value is based on the data shown in the SI from Ref. [56]. Next lines: entropy of vaporization in J mol−1 K−1 at 298 and at the temperatures of phase transition (373 ). See Section 3.4 for details. Experimental values are from Refs. [57,58,59].

	B3LYP/D3/gCP			PBE0/D3/gCP			PBEh-3c			Exp
mRRHO	−	50	100	−	50	100	−	50	100
ΔH_vap	43.74	45.89	45.93	43.57	45.61	46.26	46.31	48.20	49.03	41.58
ΔS_vap(298 K)	128.69	131.21	130.68	123.45	129.70	129.92	131.47	132.82	133.90	118.76
ΔS_vap(trs)	115.78	121.46	121.73	115.64	121.21	122.37	122.97	128.62	130.51	109.54

3.4. Entropy

Finally, the entropy of liquid water and water vapor in the liquid range and beyond the phase transition point are presented for the different methods and mRRHO treatments. Calculated and experimental values for the entropy are shown in Figure 8 and phase change data are listed in Table 3.

Figure 8

Temperature dependence of the entropy for all methods (colored lines) as well as experimental values (red marks). Standard QCE results are presented as solid, mRRHO50 results as dashed, and mRRHO100 results as dotted lines. The grey line is a linear extrapolation of the gas-phase curve. Experimental values are from Refs. [57,58,59].

It can be seen that the gas-phase entropy is hardly affected by any variation of either the quantum-chemical methods or the models for the vibrational partition function. This is plausible, since the gas phase is mainly composed of monomer units and a small number of dimers, both of which are unaffected by the mRRHO approach since they do not exhibit low-frequency vibrations. Furthermore, at such large phase volumes, the translational contributions become dominant which leads to a negligible influence of the electronic structure method.

In the liquid phase and at lower temperatures, the contributions of the individual clusters and the treatment of internal rotations become more important, leading to differences between the density functionals and between the mRRHO treatments. All three functionals underestimate the entropy at ambient temperature (298 ), and therefore, overestimate the entropy of vaporization at that temperature. This is even more pronounced if mRRHO is applied, which leads to even smaller liquid phase entropies in all cases. The effect of the cutoff value has a minor influence here, with mRRHO100 values being only slightly smaller than for mRRHO50.

4. Conclusions

In this contribution, we combined Frank Weinhold’s QCE theory with a treatment of low-frequency vibrational modes as hindered rotations rather than harmonic vibrations (mRRHO). The approach has been applied to the characterization of water, for which we computed a diverse set of properties, ranging from the microscopic effect of self-dissociation to macroscopic thermodynamic data over a range of pressures and temperatures.

We observed a significant influence of the mRRHO model on the cluster populations and derived properties such as the ionic product. The effect on thermodynamic quantities was present, but less pronounced. Since the relationship between the individual cluster partition functions and the final thermodynamic quantities of the associated phase is very complex, no clear predictions can be made in terms of which properties are affected most, how the properties are affected, or what the best cutoff value is. The influence of the mRRHO approach does indicate, however, that improving the computational chemistry tool set beyond the RRHO approximation is key to the successful description of thermodynamic data and should continue to be addressed in future studies.

The performance of QCE calculations at different pressures, and thereby, the prediction of phase-change data using the Clausius–Clapeyron relation is compelling. With the small cluster set presented in this study and by optimizing only two empirical parameters, we were able to reproduce the enthalpy of vaporization with chemical accuracy.

The QCE approach has been applied to various systems in the past (see Ref. [60]), and its areas of application are continuously being extended. Recent examples include protic ionic liquids [61], aqueous solutions of weak acids [31,32], and various organic solvents [38]. The mRRHO model has led to significant improvements in the treatment of weakly interacting aprotic compounds such as methane and chloroform. While its applicability is currently restricted to fluid phases, investigations of critical phase systems or of solids dissolved in liquids are conceivable and are possible subjects of future studies.

Once more, we conclude that the QCE approach—despite its simplicity and limited cluster set—is capable of providing a surprisingly accurate description of a complicated liquid such as water. All investigated quantities could be reproduced with astonishing agreement with the experimental reference data, rendering the QCE approach a powerful tool for future investigations of associated liquids in various contexts.