Prediction of cross second virial coefficients for dimer H2-O2 from ab initio calculations of intermolecular potentials

120 Journal of Chemistry, Vol. 46 (1), P. 120 - 126, 2008 Prediction of cross second virial coefficients for dimer H2-O2 from ab initio calculations of intermolecular potentials Received May 23, 2007 Pham Van Tat1, U. K. Deiters2 1Department of Chemistry, University of Dalat 2Institute of Physical Chemistry, University of Cologne, Germany Summary The intermolecular interaction potentials of the dimer H2-O2 were calculated from quantum mechanics, using coupled-cluster theory CCSD(T) and correlation-consistent basis sets aug-cc- pVmZ (m = 2, 3); the results were extrapolated to the basis set limit aug-cc-pV23Z. The quantum mechanical results were used to construct 5-site pair potential functions. The cross second virial coefficients of the dimer hydrogen-oxygen were obtained by integration; in these cases corrections for quantum effects were included. The results agree well with experimental data and empirical correlations. I - Introduction Computer simulation techniques, Monte Carlo as well as Molecular Dynamics, cannot work without some input, however: It is necessary to know the interaction potentials of the systems under study. The usual procedure is to assume a simple model potential. A system is the fluid mixture (H2-O2). Its thermodynamic properties are important for the design of efficient rocket engines [1], but there are remarkably few publications of experimental results only -for evident reasons. Recently an alternative approach has become feasible, for which the name “global simulation” has been coined [2]. One of the first attempts in such global simulations was that of Deiters, Hloucha and Leonhard [3] for neon to predict the vapour-liquid phase equilibria without recourse to experimental data. Further global simulation attempts for noble gases were published by the group of Huber [4]. Using a functional form for the dispersion potentials of argon and krypton proposed by Korona et al. [5]. Leonhard and Deiters constructed a 5-site Morse potential to represent the pair potential of nitrogen [8] and were able to predict vapour pressures and orthobaric densities. Nasrabad and Deiters even predicted phase high-pressure vapour-liquid phase equilibria of noblegas mixtures [7]. Bock et al. also used a 5-site pair potential for carbon dioxide [9]. Naicker et al. used SAPT (symmetry-adapted perturbation theory) to develop a 3-site pair potential for hydrogen chloride [11]; they then successfully predicted the vapour-liquid equilibria of hydrogen chloride with Gibbs ensemble Monte Carlo simulations. In this work we use quantum mechanical calculations at a sufficiently high level of approximation to obtain pair potential data of the dimer H2-O2. These data are then represented by analytical pair potential functions. These in turn are used to calculate 121 cross 2nd virial coefficients. The virial coefficients can then be compared to experimental data - as far as such data are available. II - Computational details 1. Molecular Orientation Hydrogen and oxygen molecules are represented as 5-site models, with two sites placed on the atoms (H or O), one site in the center of gravity (M), and two sites halfways between the atoms and the center (N). The molecules are treated as rigid; the interatomic distances are set to 0.74130 Å for hydrogen and 1.20741  for oxygen [12]. As hydrogen and oxygen are linear molecules, the intermolecular pair potential is a function of distance r (distance between the centers of gravity) and three angular coordinates, α, β and φ. Interaction energies were calculated for all values of r from 2.6 to 15 Å with increment 0.2 Å; the anglesα, β and φ; were varied from 0 to 180o with increment 45o. 2. Quantum chemical calculations The method CCSD(T) and the correlation- consistent basis sets of Dunning et al. [13]: aug- cc-pVDZ (for oxygen: 10s5p2d/4s3p2d, for hydrogen: 5s2p/3s2p), aug-cc-pVTZ (for oxygen: 12s6p3d2f/5s4p3d2f, for hydrogen: 6s3p2d/4s3p2d) were used in this work. The ab initio energy results were corrected for the basis set superposition error (BSSE) [14]: ΔEint = EAB - (EAb + EaB) (1) where EAB denotes the total electronic energy of a dimer AB, EAb the energy of a dimer consisting of an A atom and a B ghost atom (an atom without nucleus and electrons, but having its orbital), and EaB vice versa. The electronic energies are then extrapolated to the complete basis set limit [15]: ΔE(m) = ΔE(∞) + cm-3 (2) with m = 2 (for the aug-cc-pVDZ basis set) or 3 (for the aug-cc-pVTZ basis set). If results for two basis sets are available, it is possible to calculate the energy value for an infinite basis set from Eq.2; this result is referred to as aug- cc-pV23Z below. 3. Potential function Two our new 5-site pair potential functions were developed for dimer H2-O2: 5 5 1 2 1 1 6,8,10 0 ( ) ( ) 4 ij ij ij r i jij n e a ij ijn i j n ij ij q qCu D e f r f r r r α πε − = = = ⎡ ⎤= + +⎢ ⎥⎢ ⎥⎣ ⎦∑ ∑ ∑ (3) with 15 )2(2 1 )1()( −−−+= ijijrija erf δ and ijijrij erf β−−= 1)(2 5 5 1 2 1 1 6,8,10,12 0 ( ) ( ) 4 ij ij ij r i jij n e b ij ijn i j n ij ij q qCu D e f r f r r r α πε − = = = ⎡ ⎤= + +⎢ ⎥⎢ ⎥⎣ ⎦∑ ∑ ∑ (4) with ∑ = −−= 10 0 1 ! )( 1)( k k ijijr ijb k r erf ijij δδ and ijij r ij erf β−−= 1)(2 Here the rij site-site distances, the qi, qj electric charges of sites, and the nijC dispersion coefficients; the leading dispersion term is always proportional to r. The two models differ mostly in the choice of the damping functions f1a(rij) and f1b(rij). 4. Fitting potential function The adjustable parameters of the ab initio pair potential functions can be estimated by nonlinear least-square fitting to the ab initio interaction energy values. The fit process has to be carried out by two steps. The global minima are coarsely located by means of the Genetic algorithm (GA), and 122 the parameters resulting from the Genetic algorithm are used as initial values for the Marquardt-Levenberg algorithm. The values of root mean-square deviations (rms), multiple correlation coefficients (R2), and residuals of the fitted analytical potential functions are given in table 1. The statistical estimates here are important for assessing the fitting quality. The fitting results in Table 1 show that the residual discrepancies between the pair potentials Eq. 3 and Eq. 4 resulting from the least-square fit are insignificant. Table 1: The statistical results for fitting the intermolecular potentials Eq. 3 and Eq. 4. The values are in μEH Residual Potential R2 rms min max Eq. 3 0.9981 6.3268 -7.903 7.347 Eq. 4 0.9978 6.5214 -8.742 5.106 III - Second virial coefficients 1. Classical virial coefficients The classical viral coefficients 0clB of hydrogen resulting from the formula Eq. 7 using the ab initio 5-site pair potentials Eq. 3 and Eq. 4 are shown in table 2. Table 2: Cross second viral coefficient, B2(T) (given in cm 3/mol) of the mixture hydrogen-oxygen; correl.: empirical correlation [24]; exp.: experimental data. T/K Method 0clB 1 rB 1 aIB 1 aB m B2(T) Ref. Eq.3 -138.338 0.267 0.199 0.015 -137.857 Eq.4 -136.921 0.057 0.037 0.028 -136.798 correl. -142.097 49.8 exp. -110.0 [20] Eq.3 -61.875 0.283 0.228 0.057 -61.307 Eq.4 -60.906 0.085 0.037 0.032 -60.751 correl. -63.887 80.0 exp. -72.9 [19] Eq.3 -55.745 0.405 0.388 0.078 -54.874 Eq.4 -54.678 0.325 0.073 0.057 -54.223 correl. -56.806 85.0 exp. -54.0 [19] Eq.3 -54.073 0.314 0.177 0.047 -53.535 Eq.4 -52.976 0.128 0.073 0.014 -52.761 correl. -54.881 86.5 exp. -58.1 Eq.11 Eq.3 -50.428 0.352 0.344 0.070 -49.662 Eq.4 -49.261 0.191 0.054 0.046 -48.971 correl. -50.664 90.0 exp. -32.5 [19] 123 2. Quantum corrections In the case of hydrogen the matter is more complicated because of quantum effects. The first order corrections to the second virial coefficient of linear molecules have been worked out by Pack [10] and Wang [17]. Following the latter, the virial coefficient up to first order can be written as: ∫ ∫ ∫ ∫∫∫ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎦ ⎤⎢⎣ ⎡ +−−= 212102 21 2 )(12 11)/exp(1 2 ΩΩΩΩ dddrdruHTkTkuddu NB B B A (5) Here NA is Avogradro’s constant, kB Boltzmann’s constant, T the temperature, and u(r, α, β, φ) the pair potential; its parameters, the center-center distance and the relative orientation angles must be calculated from the center vectors ri and the absolute orientations i. H0 is the translation-rotation Hamiltonian for a pair of molecules. This expression can be broken down into a classical term and first-order quantum corrections (radial part, angular part proportional to I-1 (moment of inertia), angular part proportional to μ-1 (reduced mass)): 0 1 1 1 cl a a( ) ( ) ( ) ( ) ( )r IB T B T B T B T B Tm= + + + (6) The classical part is given by 2 0 2 cl 0 0 0 0 sin sin exp 1 4 A B N uB d d d r dr k T π π π ϕ β β α α ∞ ⎛ ⎞⎛ ⎞= − − −⎜ ⎟⎜ ⎟⎜ ⎟⎝ ⎠⎝ ⎠∫ ∫ ∫ ∫ (7) All these integrals were evaluated numerically with a 4D Gauss quadrature method [18]. Experimental values for the cross 2nd virial coefficients of the hydrogen-oxygen interaction are difficult to find in the literature. There are some experiments, however, from which these virial coefficients can be calculated: - Van Itterbeek and van Doninck measured the speed of sound in (hydrogen + oxygen) mixtures at low temperatures and pressures [19]. The pressure dependence of this property is related to the virial coefficients. The values of the cross virial coefficient obtained by these authors lie reasonably close to our predictions (see Fig. 1); it should be noted, however, that their evaluation method involved several simplifications (linearizations, neglect of temperature derivatives of the virial coefficient), and that their results exhibit an uncertainty of about 20%. - McKinley at al. measured solid-fluid equilibria of the (hydrogen + oxygen) system [20]. With the usual assumptions and simplifications (no hydrogen dissolved in the solid oxygen, neglect of higher virial coefficients) it is possible to estimate cross second virial coefficients from these data. The result agrees reasonably well with the ab initio prediction. The sublimation pressure of g- oxygen, which is required for the equilibrium calculation, was taken from the work of Roder [21]. In an earlier publication on high-pressure phase equilibria of the (hydrogen + oxygen) system it had been suggested that the parameters of the hydrogen-oxygen interaction could be interpolated from those of the the systems (hydrogen + nitrogen), (neon + nitrogen), and (neon + oxygen) [22]. This idea can be extended to second virial coefficients as follows: Eq. 7 can be simplified although with some loss of accuracy - by performing the integrations of the orientation variables: 124 Figure 1: Cross second virial coefficient of the hydrogen-oxygen system. ⎯: ab initio prediction (this work) based on Eq. 3; : ab initio prediction based on Eq. 4 (this work); V : empirical correlation [24];  : interpolation from Eq. 11; other symbols: experimental data. 0 2 0 ( , )( ) 2 exp 1A B u r TB T N r dr k T π ∞ ⎡ ⎤⎛ ⎞= − − −⎢ ⎥⎜ ⎟⎝ ⎠⎣ ⎦∫ (8) Here ( , )u r T denotes an angle-averaged pair potential. For small molecules like the ones studied here the assumption of conformal pair potentials is usually acceptable, i.e., pair potentials can be written as ( , ) ( , )u r T u r Tε= %% % with rr σ=% and Bk TT ε=% (9) where ( , )u r T%% % is a universal (reduced) pair potential function. Then Eq. 8 becomes 0 3 2 0 ( , )( ) 2 exp 1A u r TB T N r dr T π σ ∞ ⎡ ⎤⎛ ⎞= − − −⎢ ⎥⎜ ⎟⎝ ⎠⎣ ⎦∫ %% % % %% (10) We use this equation for the cross virial coefficient of the hydrogen-oxygen system, assuming 2 2 2 2 2 2; ; ;H O H N O N ε ε ε= + Δ , where the last term is supposed to be small. Taylor expansion of the Boltzmann factor of this term, truncation after the linear term, and rearrangement yield 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 0 , , , , 2 , , , , ,0 ( ) ( ) ( , ) 2 ( , ) exp 1H O H N O N H NH N A H O A H O H N H N B T B T u r T u r T r dr N N T T επσ σ ∞ ⎡ ⎤⎛ ⎞Δ≈ − − −⎢ ⎥⎜ ⎟⎜ ⎟⎢ ⎥⎝ ⎠⎣ ⎦∫ %% %%% % % %% % (11) Where the integral is a function of the reduced temperature T only. A similar equation holds for the cross virial coefficients of the (neon + oxygen) and (neon + nitrogen) systems, for which experimental data are available [23]. Therefore 2 2O ;N εΔ can be determined from the neon data and then substituted into Eq. 11 to give the second virial coefficient of (hydrogen + oxygen) at the same reduced temperature as the (neon + oxygen) system (tacitly assuming that this value also holds for the hydrogen systems). For the parameters ε and σ usual Lennard-Jones parameters [16] and Berthelot-Lorentz combining rules were used. It turns out that the hydrogen-oxygen cross virial 50 100 150 200 250 300 350 400 450 -140 -120 -100 -80 -60 -40 -20 0 20 T/K B 2 (T )/c m 3 m ol -1 125 coefficient obtained from this interpolation (-58 cm3/mol) agrees reasonably well with the ab initio predictions as well as with the experimental values (see Fig. 1). Recently Estela-Uribe and Jaramillo [24] published empirical correlation equations for second virial coefficients which are based on the corresponding-states approach of Lee and Kesler [25]. In their work, binary interactions are characterized by so-called pseudocritical parameters, which are interpolations of the pure-fluid critical temperatures and densities [24]: 1/3 1/3 1/3 , , , 1 ( ) 2 ij c ij c i c j dρ ρ ρ− − −+= + and 1/2 , , , (1 )( ) 1 / ( ) ij c i c j c ij ij k T T T c M T −= + (12) with 1 1 1 1 ( ) 2ij i j M M M− − −= + and , 1/2 , , 1 ( ) ij c ij ij c i c j a k ρ ρ ρ= − Here the Mi denote molar masses of the pure components, Mij an “interaction molar mass”, and c is a constant (21.8 K g/mol). It turns out that the adjustable parameters aij and dij are very close to zero for a large number of chemical compounds. We have set these two correlation parameters to zero for the hydrogen-oxygen interaction and used the correlations of Estela-Uribe and Jaramillo to predict cross second virial coefficients. The results show a remarkably good agreement with the predictions from quantum mechanics. IV - Conclusion We conclude that our ab initio pair potentials for the hydrogen-oxygen interaction are reliable, and that the calculation of thermodynamic properties from quantum mechanical results can be useful, if experimental data are scarce. Acknowledgments: The Regional Computer Center of Cologne (RRZK) contributed to this project by a generous allowance of computer time as well as by efficient software support; we wish to thank Dr. L. Packschies for technical help with the Gaussian03 software. References 1. V. P. Dawson. A History of the Rocket Engine Test Facility at the NASA Glenn Research Center, NASA, Glen Research Center, Cleveland (2004). 2. H. Popkie, H. Kistenmacher, and E. Clementi. J. Chem. Phys., 59, 1325 - 1336 (1973). 3. U. K. Deiters, M. Hloucha, and K. Leonhard. Experiments?. No, thank you!, in T. M. Letcher, editor, Chemistry for the 21st Century: Chemical Thermodynamics, IUPAC Monograph Series, pages 187-195. Blackwell Science, Oxford (1999). 4. M. Venkatraj, C. Bratschi, H. Huber. Fluid Phase Equilib., 218, 285 - 289 (2004). 5. T. Korona, H. L. Williams, R. Bukowski, B. Jeziorski, and K. Szalewicz. J. Chem. Phys., 106, 5109 - 5122 (1997). 6. A. E. Nasrabad and U. K. Deiters. J. Chem. Phys., 119, 947 - 952 (2003). 7. A. E. Nasrabad, R. Laghaei, and U. K. Deiters. J. Chem. Phys., 121, 6423 - 6434 (2004). 8. K. Leonhard and U. K. Deiters. Mol. Phys., 100, 2571 - 2585 (2002). 9. S. Bock, E. Bich, and E. Vogel. Chem. Phys., 257, 147 - 156 (2000). 10. R. T. Pack. J. Chem. Phys., 78, 7217 - 7222 (1983). 126 11. P. K. Naicker, A. K. Sum, and S. I. Sandler. J. Chem. Phys., 118, 4086 - 4093 (2003). 12. L. E. Sutton, editor. Table of Interatomic Distances and Configurations in Molecules and Ions: supplement 1956-1959. Number 18 in Special publication. Chemical Society, London (1965). 13. R. A. Kendall, T. H. Dunning, Jr., and R. J. Harrison. J. Chem. Phys., 96, 6796 - 6806 (1992). 14. S. F. Boys and F. Bernardi. Mol. Phys., 19, 553 - 566 (1970). 15. S. Y. Park and L. J. Sh. J. Chem. Phys., 116, 5389 - 5394 (2002). 16. J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird. Molecular Theory of Gases and Liquids. John Wiley, New York (1954). 17. W. F. Wang, J. Quant. Spectrosc. Radiat. Transfer, 76, 23 - 30 (2003). 18. W. Squire. Integration for Engineers and Scientists. Elsevier, New York (1970). 19. A. van Itterbeek and W. van Doninck, Proc. Phys. Soc. London B, 62, 62 - 69 (1949). 20. C. McKinley, J. Brewer, and W. E. S. J., Adv. Cryogen. Eng., 7, 114 - 124 (1961). 21. H. M. Roder. The Thermodynamic Properties of Slush Hydrogen and Oxygen. NBSIR 77-859. National Bureau of Standards, Boulder (1977). 22. U. K. Deiters, M. Neichel, E. U. Franck, Ber. Bunsenges. Phys. Chem., 97, 649 - 657 (1993). 23. J. H. Dymond and E. B. Smith, The Virial Coefficients of Pure Gases and Mixtures. Clarendon Press, Oxford (1980). 24. J. F. Estela-Uribe and J. Jaramillo, Fluid Phase Equilib., 231, 84 - 98 (2005). 25. B. I. Lee and M. G. Kesler, AIChE J., 21, 510 - 527 (1975).