Cross Second Virial Coefficients of the H2O–H2 and H2S–H2 Systems from First-Principles

The cross second virial coefficients B12 for the interactions of water (H2O) with molecular hydrogen (H2) and of hydrogen sulfide (H2S) with H2 were obtained at temperatures in the range from 150 to 2000 K from new intermolecular potential energy surfaces (PESs) for the respective molecule pairs. The PESs are based on interaction energies determined for about 12 000 configurations of each molecule pair employing different high-level quantum-chemical ab initio methods up to coupled cluster with single, double, triple, and perturbative quadruple excitations [CCSDT(Q)]. Furthermore, the interaction energies were corrected for scalar relativistic effects. Both classical and semiclassical values for B12 were extracted from the PESs using the Mayer-sampling Monte Carlo approach. While our results for the H2O–H2 system validate the older first-principles results of Hodges et al. [J. Chem. Phys. 2004, 120, 710–720], B12 for the H2S–H2 system was, to the best of our knowledge, hitherto neither measured experimentally nor predicted from first principles.


INTRODUCTION
The prediction of thermophysical properties of a molecular fluid requires detailed information about the interactions between the individual molecules.This information is contained in the socalled potential energy surface (PES), which is the function relating the positions, spatial orientations, and internal degrees of freedom of the molecules forming the fluid to the latter's total potential energy.If we are dealing with a low-density gas and each molecule is approximated as a rigid rotor, we can formulate the PES as a sum of pair PESs, each of which is a function only of the separation and mutual orientation of the two interacting molecules.Quantum-chemical ab initio methods such as the popular coupled cluster approach with single, double, and perturbative triple excitations [CCSD(T)] 1 can be applied to obtain very accurate interaction energies for pairs of small molecules that consist only of a few atoms.Once we have sampled the whole space of pair configurations sufficiently densely by such ab initio calculations, we can fit a suitable mathematical function to the calculated interaction energies to obtain the PES in analytical form.
An important thermophysical property that can be extracted readily both classically and semiclassically from rigid-rotor pair PESs is the second virial coefficient B, which determines the lowest-order correction to the ideal gas law in the virial expansion: (1)   where p is the pressure, ρ m is the molar density, R is the molar gas constant, T is the temperature, and C is the third virial coefficient.B is a genuine two-body property and can be obtained as a function of T from statistical thermodynamics by numerically solving a multidimensional integral using a suitable integration algorithm.The virial expansion is also valid for gas mixtures.In a binary mixture with components 1 and 2, the total second virial coefficient B(T) is given through the exact relation (2)   where x 1 and x 2 are the mole fractions, B 11 and B 22 denote the corresponding second virial coefficients of the pure components, and B 12 is the cross second virial coefficient (sometimes also called interaction second virial coefficient).Cross second virial coefficients can be extracted from the respective unlike-molecule pair PESs in a manner similar to that for the pure-component second virial coefficients from the like-molecule pair PESs.
In a number of recent studies, we applied the above-described first-principles methodology to predict cross second virial coefficients B 12 for several pairs of small molecules, such as water (H 2 O) with nitrogen (N 2 ), oxygen (O 2 ), carbon monoxide (CO), carbon dioxide (CO 2 ), hydrogen sulfide (H 2 S), and sulfur dioxide (SO 2 ), see ref 2 and references therein, or H 2 S with N 2 3 and CO 2 . 4The present work extends our series of first-principles investigations on cross second virial coefficients to two systems containing molecular hydrogen (H 2 ), namely, H 2 O−H 2 and H 2 S−H 2 .Reliable thermophysical property data for these systems are important for the hydrogen economy.For example, data for the H 2 O−H 2 system are needed for modeling combustion in hydrogen gas turbines.In other applications, small amounts of H 2 O are present in H 2 gas as an impurity whose effects need to be accurately modeled.Property data for H 2 S−H 2 mixtures are also of interest since H 2 S separated from sour natural gas is a potential source for producing H 2 through decomposition.Ab initio PESs for both systems (and for H 2 O−H 2 also derived B 12 values) already exist, see, e.g., refs 5−10.However, as in our previous studies on cross second virial coefficients, we developed new pair PESs, applying higher levels of theory for the ab initio calculations.The PES models used in this work are summarized in Table 1.We stress again that the respective pure-component PESs (see, e.g., refs 11−13 for ab initio PESs from the literature) are not required for the calculation of the cross second virial coefficients.
The paper is organized as follows: In section 2, we present the new ab initio pair PESs.The details of the calculations of B 12 in the temperature range from 150 to 2000 K are provided in section 3.In section 4, we present the results for B 12 both as a data table of discrete values and as practical correlations.We also compare our results with the available experimental and firstprinciples-based values from the literature for B 12 and the related dilute gas cross isothermal Joule−Thomson coefficient ϕ 12 .Finally, our conclusions are given in section 5.

Monomer Geometries and Multipole Moments.
In this work, the monomers H 2 O, H 2 S, and H 2 were approximated as rigid rotors, meaning that all of the bond lengths and bond angles were kept at fixed values.The resulting reduction in the dimensionality of the H 2 O−H 2 and H 2 S−H 2 PESs from 9 to 5 comes at the expense of a somewhat lower (but still very good) accuracy for the derived B 12 values.The equilibrium structures (i.e., the minima of the intramolecular PESs as obtained by standard geometry optimizations) are not an optimal choice for the fixed monomer geometries.Instead, the zero-point vibrationally averaged structures are recommended to be used 14 because they at least partly account for the effects of the suppressed degrees of freedom.The respective H 2 O geometry was taken from ref 14 and is characterized by an OH bond length of 0.9716257 Å (1 Å = 10 −10 m) and an HOH angle of 104.69°.The averaged H 2 S geometry was taken from our work on the H 2 S−H 2 S system 13 and is characterized by an SH bond length of 1.3506 Å and an HSH angle of 92.219°, while for the averaged bond length of the H 2 molecule a value of 0.766638 Å 15,16 was adopted.
Electric multipole moments for the three monomers were also determined for the purpose of using them as constraints in the fits of the electrostatic parts of the analytical potential functions (see section 2.3).For H 2 , we determined the quadruple moment at the averaged bond length by interpolating the separationdependent quadrupole moment values calculated by Komasa and Thakkar, 17 yielding a value of 0.482390 au.In our recent work on the H 2 O−H 2 S and H 2 O−SO 2 systems, 2 we already determined values for the dipole moments and for the components of the quadrupole and octopole moment tensors of H 2 O and H 2 S at the CCSD(T) 1 level of theory with scalar relativistic corrections, and we do not repeat the details of these calculations here.However, because the highest level of theory used for the present ab initio calculations of H 2 O−H 2 and H 2 S− H 2 interaction energies is higher than CCSD(T), it is appropriate to also include a post-CCSD(T) correction for the multipole moments of H 2 O and H 2 S for consistency.These corrections were obtained as the differences between the respective multipole moment values obtained at the frozen-core (FC) CCSDTQ 18,19 and CCSD(T) levels of theory using the aug-cc-pVTZ 20−22 basis set with the analytical derivative approach implemented in the CFOUR 23,24 quantum chemistry code (version 2.1) with orbital relaxation included.Because the required derivatives cannot be computed at the CCSDTQ level by the employed CFOUR version, the general coupled cluster code MRCC 25,26 (2020 version) interfaced with CFOUR was used for these calculations.The relative changes in the multipole moment values due to the post-CCSD(T) correction were found to be very small, being within ±0.23% for both molecules.The details are provided in the Supporting Information.

Pair Geometries and Interaction
Energies.The geometries of the systems H 2 O−H 2 and H 2 S−H 2 within the rigid-rotor approximation can be described by the separation between the centers of mass of the two molecules, R, and the four Euler angles θ 1 , ψ 1 , θ 2 , and ϕ.This is illustrated exemplarily for the H 2 O−H 2 system in Figure 1, with further details given in the Supporting Information.To exhaustively sample the angular configuration space for both systems, we systematically varied θ 1 in steps of 22.5°starting with 0°, ψ 1 and θ 2 in steps of 30°s tarting with 0°, and ϕ in steps of 45°starting with 0°.The resulting set of angular configurations was checked for configurations that are related to one another through symmetry operations and are, thus, energetically equivalent.By discarding the redundant configurations, we obtained 510 unique angular configurations.For each angular configuration, we considered 24 center-of-mass separations R in the range from 1.2 to 10.0 Å, yielding 12 240 (510 × 24) pair configurations.Some configurations at very small separations were discarded immediately because of excessive overlap of the two molecules, while other configurations, also at small separations, had to be discarded because of unresolvable convergence problems in the quantum-chemical calculations.The number of remaining pair configurations was 12 097 for the H 2 O−H 2 system and 12 092 for the H 2 S−H 2 system.We used the counterpoise-corrected supermolecular approach 27 to calculate the interaction energies V for all pair configurations.An additivity scheme combining different levels of theory was employed in order to achieve the highest possible accuracy for the final interaction energies.As the lowest level of theory, for which the largest basis sets could be used, we chose the resolution of identity second-order Møller−Plesset perturbation theory (RI-MP2) 28,29 within the FC approximation, with the RI-JK approximation 30,31 applied in the Hartree− Fock (HF) step.The calculations were carried out using the augcc-pVXZ 20−22 basis sets with X = 4 (Q) and X = 5 together with the auxiliary basis sets aug-cc-pV5Z-JKFIT 32 and aug-cc-pV5Z-MP2FIT 33 for both values of X.Following Patkowski, 34 we also placed a "ghost" hydrogen atom with the same basis and auxiliary basis functions as for the real hydrogen atoms midway along the R axis of each pair configuration because such so-called bond functions usually significantly accelerate convergence toward the complete basis set (CBS) limit.The correlation contributions to the computed interaction energies, V RI−MP2 corr X , were extrapolated to the CBS limit with the widely used X −3 scheme: 35,36 (3) whereas the HF contributions were taken from the calculations at the X = 5 basis set level without extrapolation.
In the next step, we calculated the differences between the CCSD(T) 1 and the standard MP2 interaction energies within the FC approximation using the aug-cc-pVXZ 20−22 basis sets with X = 3 (T) and X = 4 (in both cases without bond functions).The differences, denoted as ΔV CCSD(T) X , were extrapolated to the CBS limit again using the X −3 scheme: (4)   We then added ΔV CCSD(T) CBS to the total RI-MP2 interaction energies, yielding values that should be very close (probably within 1% for most orientations in the well regions of the PESs) to the true CBS limits of the total CCSD(T) interaction energies within the FC approximation.
Usually, CCSD(T) or a variant thereof is the highest level of theory applied in the development of intermolecular PESs, and this is also the case for the existing PESs for the H 2 O−H 2 and H 2 S−H 2 systems. 5−7,9,10 However, because the H 2 molecule has only two electrons, the total number of electrons in the H 2 O−H 2 and H 2 S−H 2 molecule pairs is small enough to make calculations at still higher levels of theory, in particular CCSDT(Q), 37,38 with reasonably sized basis sets computationally feasible today.We calculated the CCSDT(Q) correction for all configurations by taking the differences between the interaction energies obtained at the CCSDT(Q) and CCSD(T) levels within the FC approximation employing the aug-cc-pVXZ 20−22 basis sets with X = 2 (D) and X = 3, followed again by extrapolation to the CBS limit: (5)   For both molecule pairs, ΔV CCSDT(Q) CBS is found to be always negative.Its relative influence on the total interaction energy exceeds 1%, 2%, and 4% for about 48%, 29%, and 17%, respectively, of all H 2 O−H 2 configurations and for about 61%, 37%, and 17%, respectively, of all H 2 S−H 2 configurations.
Finally, corrections were calculated for the inclusion of core− core and core−valence correlations, ΔV core , and scalar relativistic effects, ΔV rel .We determined ΔV core at the CCSD(T)/aug-cc-pwCVTZ [20][21][22]39 level by computing the differences between the interaction energies obtained with all but the sulfur 1s electrons correlated and the interaction energies resulting with the FC approximation.The relativistic corrections ΔV rel were determined at the CCSD(T) level within the FC approximation using the aug-cc-pVTZ basis set in uncontracted form by taking the differences between the interaction energies resulting with the spin-free exact two-component theory at the one-electron level 40 and the nonrelativistic interaction energies. The cmbined effect of ΔV core and ΔV rel on the total H 2 O−H 2 interaction energies is very small and essentially negligible.For half of all investigated configurations, the impact is less than 0.15%, while it is higher than 1% for only about 5% of the configurations.Not surprisingly, the two corrections are more important for the H 2 S−H 2 system; their combined impact on the interaction energy exceeds 1%, 2%, and 4% for about 42%, 21%, and 10%, respectively, of the investigated configurations.
Detailed results of the ab initio calculations for all 12 097 H 2 O−H 2 and 12 092 H 2 S−H 2 configurations are provided in the Supporting Information.The ORCA 41,42 program (version 3.0.3)was used to perform the RI-MP2 calculations, while all other calculations were carried out with CFOUR 23,24 (version 2.1).

Analytical Potential Functions.
The two new pair PESs are represented analytically by site−site functions.We chose an arrangement with nine sites for both the H 2 O and H 2 S molecules.All of the sites are located in the molecular planes, and in both molecules, three of the nine sites are located directly on the C 2v symmetry axis.For the smaller H 2 molecule, five interaction sites were deemed sufficient.With these site arrangements, due to symmetry, H 2 O and H 2 S each only have six distinct types of sites, while H 2 only has three distinct types of sites.The number of distinct types of site−site interactions and the total number of site−site interactions in each molecule pair are thus 18 (6 × 3) and 45 (9 × 5), respectively.We represent each of the site−site interactions by the following function: (6)   where R ij is the separation between site i in H 2 O or H 2 S and site j in H 2 , and the damping functions f 6 are those of Tang and Toennies: 43 (7)   The total pair potential is obtained by summing over all site−site contributions: (8)   We determined the parameters A, α, b, and C 6 for the distinct types of site−site combinations (72 parameters in total per molecule pair), the precise locations of the sites in the molecules, and the site charges q (fixed at zero for one of the H 2 O sites) by performing nonlinear least-squares fits to the 12 097 H 2 O−H 2 and 12 092 H 2 S−H 2 ab initio interaction energies.The fits were constrained such that the molecules have zero charge and that the ab initio multipole moments of the monomers (see section 2.1) are reproduced.Furthermore, to ensure that the fit quality is appropriately balanced over the entire PESs, the ab initio interaction energies were weighted with the following function: (9)   where a = 350 for the H 2 O−H 2 PES and a = 220 for the H 2 S−H 2 PES.The denominator in eq 9 increases the weight of configurations as the interaction energy decreases toward the most negative values (V > −350 K and V > −220 K for all investigated H 2 O−H 2 and H 2 S−H 2 configurations, respectively), whereas the numerator ensures that the fit quality does not seriously deteriorate at large R values.We used similar weighting functions previously (e.g., in ref 2).Note that throughout this work, we quote energies in units of kelvin; i.e., we divided the energies by the Boltzmann constant k B .However, we omit k B from the notation for brevity.
Figure 2 visualizes the positions of the interaction sites relative to the atoms in the molecules for both potential functions.
Although barely discernible, the positions of the five sites in H 2 are different for the H 2 O−H 2 and H 2 S−H 2 potential functions.
Figure 3 depicts the deviations of the interaction energies obtained with the fitted potential functions from the corresponding ab initio calculated ones.With the relative deviations being mostly within ±2%, the fit quality can be considered as sufficient for accurate calculations of thermophysical properties because the effects on the calculated property values by fitting errors with positive and negative sign are expected to mostly cancel out.We restricted Figure 3 to interaction energies lower than 8000 K because higher energies are of quickly diminishing relevance for the present applications of the potential functions.For both PESs, the full range of ab initio interaction energies extends to more than 160 000 K. A visual inspection of many one-dimensional cuts through the potential functions along R and the four angles showed that the functions are very smooth; i.e., there are no apparent unphysical "wiggles," which can occur when a function with too many adjustable parameters is fitted to too few ab initio points.While the potential functions, due to the chosen form and applied constraints, behave correctly at R values above the investigated range, the behavior of the potential functions at the very small R values below the investigated range is partly unphysical.This is manifested by the occurrence of spurious saddle points for some angular configurations, so that V actually decreases as R is further decreased, often to the point that V becomes negative.Such a "hole" is problematic if the saddle point separating it from the rest of the PES is thermally accessible at the temperatures at which we calculate the thermophysical properties.Fortunately, even at the highest temperature investigated in this work (2000 K), the lowest of these spurious saddle points with interaction energies of about 95 000 K (H 2 O−H 2 ) and 73 000 K (H 2 S−H 2 ) are still completely inaccessible.As described in section 3, the holes can be relatively easily addressed in this case.   .Deviations of interaction energies obtained using the new potential functions from the corresponding ab initio calculated interaction energies as a function of the latter for positive interaction energies up to 8000 K (upper panels) and negative interaction energies (lower panels).The dashed red lines are a guide to the eye indicating relative deviations of ±2%. the minima on the H 2 O−H 2 and H 2 S−H 2 PESs of the present work, Figure 6 shows contour plots of V as a function of the angles θ 1 and θ 2 , with the remaining three variables chosen such that V is minimized.It can be seen that for both PESs the saddle point separating the global minimum and the secondary minimum lies energetically only slightly above the secondary minimum.The full details regarding the minima of the potential functions are given in the Supporting Information.Fortran 90 programs that compute the potential functions are also provided there.

CALCULATION OF THE CROSS SECOND VIRIAL COEFFICIENTS
Classical statistical thermodynamics provides a relatively simple expression for the cross second virial coefficient of a pair of rigid molecules: (10)   where N A is Avogadro's constant, R is the separation vector between the two molecules, and the angle brackets indicate proper averaging over the spatial orientations of the two molecules, which we represent by the variables Ω 1 and Ω 2 .
Particularly due to the H 2 molecule's very small mass and moment of inertia, it is essential to account for nuclear quantum effects, which in this work we did with two well-established semiclassical schemes.The quadratic Feynman−Hibbs (QFH) approach 44 adds a temperature-dependent term to the potential V:   Here, ℏ is Planck's constant h divided by 2π; μ denotes the reduced mass of the molecule pair; x, y, and z are the Cartesian components of R; I 1,a , I 1,b , and I 1,c are the three principal moments of inertia of molecule 1; I 2 is the moment of inertia of molecule 2 perpendicular to its bond axis; and the angles ψ 1,a , ψ 1,b , ψ 1,c , ψ 2,a , and ψ 2,b describe rotations around the principal axes of the molecules with the exception of the bond axis of molecule 2. The second semiclassical approach that we considered differs from the QFH approach in that it does not modify the potential V but instead modifies the classical expression by introducing its first-order quantum correction: (13)   where ΔV(T) has the same form as for the QFH approach, although it is possible to simplify it in this case through integration by parts such that it involves only first derivatives.
Both the QFH approach and the first-order quantum correction approach (which we refer to as 1qc approach) are correct only up to order ℏ 2 with respect to the complete semiclassical expansion of B 12 in powers of ℏ 2 .While apparently B 12 1qc does not contain any contributions of orders ℏ 4 , ℏ 6 , and so forth, B 12 QFH does due to ΔV(T) appearing inside the exponential.These higher-order contributions can be viewed as extremely crude estimates of the true higher-order quantum corrections, so that in most cases, B 12 QFH can be expected to be closer to the fully quantum-mechanical result than B 12 1qc .To evaluate B 12 numerically for both molecule pairs at the classical and the two semiclassical levels of theory, we used the Mayer-sampling Monte Carlo (MSMC) approach devised by Singh and Kofke. 45The calculations were performed with the same in-house code as in all of our previous studies on cross second virial coefficients.During the MSMC runs, B 12 averages were accumulated simultaneously for 99 temperatures from 150 to 2000 K by applying the multitemperature variant of the MSMC approach, 46 with the temperature governing the sampling distribution being 300 K. To avoid sampling of the holes in the PESs at very small separations, we placed hard spheres with a diameter of 1 Å on all interaction sites in the molecules and checked that the spheres are large enough to cover the holes completely but not so large that their presence introduces systematic errors into the computed B 12 values.Each attempted MSMC move consisted of a random displacement and rotation of one of the two molecules, which was also chosen randomly.The maximum step sizes were iteratively adjusted in short equilibration runs until an acceptance rate of close to 50% was achieved.The derivatives of the pair PESs needed for the semiclassical schemes were computed analytically.As the reference system required for the MSMC approach, we used simple hard spheres of 4.5 Å diameter placed at the centers of mass of the molecules.For each molecule pair and level of theory, 16 independent simulation runs, each comprising 5 × 10 10 attempted moves, were carried out, and the resulting B 12 values for each temperature were averaged.The statistical standard uncertainties of the final B 12 values due to the integration are always smaller than 0.001 cm 3 •mol −1 and are thus so small that they are irrelevant to the overall uncertainty assessment.

RESULTS AND DISCUSSION
Table 2 lists the values for the cross second virial coefficients obtained with the new pair PESs at the classical and the semiclassical 1qc and QFH levels of theory for the H 2 O−H 2 and H 2 S−H 2 systems at 47 of the 99 investigated temperatures.In addition, the table provides values for a "tuned" QFH level for the H 2 O−H 2 system: (14)   where the factor 0.6 is the result of approximately matching the differences between B 12 QFH−T and B 12 1qc at 150, 200, and 250 K with the differences between the fully quantum-mechanically calculated B 12 values of Hodges et al. 5 and their respective 1qc values at these temperatures.This tuning approach is particularly justified because the first-order quantum corrections obtained in this work and by Hodges et al. agree almost perfectly.
Finally, Table 2 lists the estimated combined expanded uncertainties (coverage factor k = 2, corresponding approximately to a 95% confidence level) of the B 12 QFH−T values for the H 2 O−H 2 pair and of the B 12 QFH values for the H 2 S−H 2 pair.The uncertainty due to the imperfections of the pair PESs stems mostly from their reduced-dimensionality treatment, i.e., the rigid-rotor approximation.Both the level of theory for the individual PES points and the quality of the analytical fits should not be significant error sources for the final B 12 values.However, the semiclassical evaluation of B 12 contributes to some uncertainty at the lowest temperatures.Taking these considerations and experience from previous works into account, we conservatively estimated the combined expanded uncertainties to be (15)   for the H 2 O−H 2 system and (16)   for the H 2 S−H 2 system.
To provide simple correlations for practical applications, we fitted polynomials in T −1/2 to the 99 calculated B 12 QFH−T and B 12 QFH values for the H 2 O−H 2 and H 2 S−H 2 systems, respectively.The polynomial structures were optimized with the Eureqa software (version 1.24.0). 47The resulting correlations, in which T* = T/(100 K), are given by (17)   for the H 2 O−H 2 system and (18)   for the H 2 S−H 2 system.The correlations reproduce the calculated B 12 values within 0.04 cm 3 •mol −1 and extrapolate to temperatures below 150 K and above 2000 K in a physically reasonable manner.Figure 8 shows the dilute gas cross isothermal Joule− Thomson coefficient ϕ 12 of the H 2 O−H 2 system derived from the B 12 correlations of this work and of Hodges et al. 5 and Scribano et al. 8 using the relation (19)   The figure also displays the available experimental ϕ 12 data, 53,55 with those for ref 55 having been obtained by Hodges et al. 5 from the measured enthalpies of mixing given in that reference.Wormald and Lancaster 53 used their ϕ 12 data to derive the above-mentioned B 12 values that we omitted from Figure 7 because these B 12 values cannot provide more information than the ϕ 12 data they are based on.It is thus also not surprising that the pattern of the deviations for B 12 mirrors that for ϕ 12 .As in the case of B 12 , the agreement between the present ϕ 12 values and those of Hodges et al. is excellent, while the values for the correlation of Scribano et al. are noticeably higher and the quality of the experimental data is again insufficient for using them to assess the quality of the calculated values.
For the H 2 S−H 2 system, to the best of our knowledge, neither B 12 nor ϕ 12 were previously calculated from first-principles or measured experimentally.

CONCLUSIONS
The cross second virial coefficients B 12 of the H 2 O−H 2 and H 2 S−H 2 systems were determined employing a proven firstprinciples methodology.In the first step, new H 2 O−H 2 and H 2 S−H 2 pair PESs were developed using the counterpoisecorrected supermolecular approach 27 with high-level ab initio methods up to CCSDT(Q) 37,38 and with consideration of scalar relativistic effects.The PESs are represented analytically by site− site functions, which we provide in the form of Fortran 90 codes as part of the Supporting Information.
In the second step, we extracted B 12 values for the temperature range from 150 to 2000 K from the analytical potential functions classically and with two semiclassical approaches using the MSMC technique. 45Our final recommended results are provided together with conservative uncertainty estimates both as a data table of discrete values and in the form of simple correlations, eqs 17 and 18, which not only represent the calculated values very accurately but also extrapolate in a reasonable manner to temperatures far below and above the investigated range.
Our B 12 values for the H 2 O−H 2 system are in excellent agreement with the first-principles values of Hodges et al., 5 whereas the more recent first-principles values of Scribano et al. 8 are significantly higher.The available experimental B 12 data for  The values of this work were calculated using eq 17 and those of Hodges et al. 5 and Scribano et al. 8 using the correlations provided in these papers.The experimental ϕ 12 data of Wormald and Lancaster 53 were given in the original paper, while those of Lancaster and Wormald 55 were derived by Hodges et al. 5

Figure 1 .
Figure 1.Visualization of the variables R, θ 1 , ψ 1 , θ 2 , and ϕ used to describe H 2 O−H 2 and, in an analogous manner, H 2 S−H 2 configurations.The center of mass of the H 2 O molecule is located at the origin of a Cartesian coordinate system defined by the axes x, y, and z, while the H 2 molecule's center of mass is at (x, y, z) = (0, 0, R).The moleculefixed axes z 1 and z 2 are the symmetry axis of H 2 O and the bond axis of H 2 , respectively.The molecule-fixed axis y 1 lies in the molecular plane of H 2 O and is perpendicular to z 1 ; it crosses z 1 in the center of mass of the H 2 O molecule, i.e., at (x, y, z) = (0, 0, 0).
Figures 4 and 5 show cuts along R through the H 2 O−H 2 and H 2 S−H 2 potential functions, respectively, for 12 selected angular configurations.The figures demonstrate that the anisotropy of the pair interactions is especially strong with respect to the well depth.

Figure 2 .
Figure 2. Optimized positions of the nine H 2 O, nine H 2 S, and five H 2 interaction sites for the two potential functions.The sites are all located in the molecular planes.

Figure 3
Figure3.Deviations of interaction energies obtained using the new potential functions from the corresponding ab initio calculated interaction energies as a function of the latter for positive interaction energies up to 8000 K (upper panels) and negative interaction energies (lower panels).The dashed red lines are a guide to the eye indicating relative deviations of ±2%.

Figure 4 .
Figure 4. H 2 O−H 2 interaction energies as a function of the center-of-mass separation R for 12 of the 510 investigated angular configurations.Symbols indicate the discrete ab initio interaction energies, while solid lines indicate the interaction energies obtained for the respective angular orientations with the fitted potential function.

Figure 5 .
Figure 5. H 2 S−H 2 interaction energies as a function of the center-of-mass separation R for 12 of the 510 investigated angular orientations.The meaning of the symbols and lines is the same as that in Figure 4.

Figure 6 .
Figure 6.Contour plots of V (in kelvin) for the analytical H 2 O−H 2 and H 2 S−H 2 PESs of this work as a function of the angles θ 1 and θ 2 , with the centerof-mass separation R and the angles ψ 1 and ϕ chosen such that V is minimized for each combination of θ 1 and θ 2 .The minima of the PESs (apart from the dubious third minimum of the H 2 S−H 2 PES) are marked in the plots.

Figure 7
depicts the B 12 values for the H 2 O−H 2 system obtained in this work at the tuned QFH level and those obtained by Hodges et al.5 and Scribano et al.,8 who each provided their results in the form of a correlation.Hodges et al.'s correlation is based on fully quantum-mechanical calculations of B 12 at temperatures up to 250 K and semiclassical calculations at higher temperatures using their own ab initio PES, while Scribano et al.'s correlation is based solely on semiclassical B 12 values at the 1qc level calculated using the ab initio PES of Valiron et al. 7 The figure shows an excellent agreement between the present results and those of Hodges et al., whereas Scribano et al.'s values are somewhat higher and below ambient temperature also outside the estimated expanded uncertainty range of our values.More positive virial coefficient values indicate that the PES is on average less attractive, which is in line with the discussion in section 2.3 regarding the depths of the minima on the PESs of Valiron et al. and of this work.In contrast, the maximum well depths of the present PES and of that of Hodges et al. are very close, so that assuming that this is not a fluke, it is unsurprising that the resulting B 12 values are also in excellent agreement.Apart from Wormald and Lancaster, 53 whose B 12 values we briefly discuss below in connection with Figure 8, Hodges et al. 5 provided the only experimentally based B 12 values available in the literature, which they derived from other measured quantities given in refs 48−52 and which are also shown in Figure 7. Since the scatter and the uncertainties of the data exceed the uncertainties of the calculated values substantially, the data are not suitable for verifying the calculated values.Note that we repeated Hodges et al.'s analysis of the data of Seward et al. 52 for the second virial coefficients of different H 2 O−H 2 mixtures using newer values for the second virial coefficient of pure H 2 , 54 but the changes in the B 12 values due to the reanalysis do not exceed 0.1 cm 3 •mol −1 , which would not even be visible in the figure.

Figure 7 .
Figure 7. Cross second virial coefficient B 12 of the H 2 O−H 2 system as a function of temperature T. The calculated values of this work, for which the gray shading indicates the expanded (k = 2) uncertainty, correspond to the tuned QFH level of theory, while those of Hodges et al. 5 (barely distinguishable from our values) correspond to semiclassical and, at the lowest temperatures, fully quantum-mechanical levels of theory.Scribano et al. 8 calculated their B 12 values at the 1qc level.All three sets of calculated values are based on different ab initio PESs.The experimental B 12 data shown were derived by Hodges et al. 5 from other measured quantities provided in refs 48−52.

Figure 8 .
Figure 8. Dilute gas cross isothermal Joule−Thomson coefficient, ϕ 12 = B 12 − T(dB 12 /dT), of the H 2 O−H 2 system as a function of temperature T.The values of this work were calculated using eq 17 and those of Hodges et al.5 and Scribano et al.8 using the correlations provided in these papers.The experimental ϕ 12 data of Wormald and Lancaster53 were given in the original paper, while those of Lancaster and Wormald55 were derived by Hodges et al.5 from the measured enthalpies of mixing given in ref55.

Table 1 .
PES Models Used in This Work

Table 2 .
Cross Second Virial Coefficients B 12 (in cm 3 •mol −1 ) of the H 2 O−H 2 and H 2 S−H 2 Molecule Pairs at Selected Temperatures T at the Classical, the Semiclassical 1qc, the Semiclassical QFH, and, for H 2 O−H 2 , the "Tuned" Semiclassical QFH from the measured enthalpies of mixing given in ref55.theH 2 O−H 2 system are far less accurate than any of the firstprinciples sets of B 12 values.In the case of the H 2 S−H 2 system, we are not aware of any previous first-principles or experimental determinations of B 12 .sıSupportingInformationThe Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jced.3c00300.Further details on the internal coordinates of the H 2 O− H 2 and H 2 S−H 2 molecule pairs, detailed results of the ab initio calculations of the pair interaction energies and the multipole moments of H 2 O and H 2 S, the details of the distinct minimum structures obtained with the new analytical pair PESs, and Fortran 90 implementations of the PESs (ZIP) Institut fur Thermodynamik, Helmut-Schmidt-Universität/Universität der Bundeswehr Hamburg, 22043 Hamburg, Germany; orcid.org/0000-0003-3121-6827;Email: robert.hellmann@hsu-hh.deComplete contact information is available at: https://pubs.acs.org/10.1021/acs.jced.3c00300 *