Monte carlo simulation of vapor - Liquid equilibria of hydrogen using ab initio intermolecular potentials

529 Journal of Chemistry, Vol. 47 (5), P. 529 - 534, 2009 MONTE CARLO SIMULATION OF VAPOR - LIQUID EQUILIBRIA OF HYDROGEN USING AB INITIO INTERMOLECULAR POTENTIALS Received 6 September 2007 PHAM VAN TAT1, U. K. DEITERS2 1Department of Chemistry, University of Dalat 2University of Cologne, Germany abstract The vapor-liquid equilibria of pure fluid hydrogen were predicted by Gibbs ensemble Monte Carlo simulation techniques using our two different ab initio intermolecular pair potentials. The ab initio intermolecular interaction potentials were obtained from coupled-cluster calculations, using the CCSD(T) level of theory and Dunning's correlation consistent basis sets aug-cc-pVmZ (m =2, 3) [9]. The phase diagram, critical properties, thermodynamic properties, vapor pressures and orthobaric densities based on them are found to be in very good agreement with experimental data. Keywords: Vapor-liquid equilibria, Gibbs ensemble Monte Carlo simulation, ab initio potentials. I - INTRODUCTION Hydrogen and the mixture hydrogen-oxygen are used in several industrial applications. It could become the most important energy carrier of tomorrow. Liquid hydrogen, oxygen are the usual liquid fuels for rocket engines [1]. The National Aeronautics and Space Administration (NASA) is the largest user of liquid hydrogen in the world [2]. Computer simulations have become indispensable tools for studying pure fluids and fluid mixtures. One of the first attempts Nasrabad and Deiters predicted phase high- pressure vapour-liquid phase equilibria of noble- gas mixtures [3,4] from the global simulations using the intermolecular potentials. Leonhard and Deiters used a 5-site Morse potential to represent the pair potential of nitrogen [5] and were able to predict vapour pressures and orthobaric densities. Naicker et al. developed the 3-site pair potentials for hydrogen chloride [6]; they predicted successfully the vapour- liquid phase equilibria of hydrogen chloride with GEMC (Gibbs Ensemble Monte Carlo simulations [7]). In this work we report the simulation results of the vapor-liquid equilibria for the pure fluid hydrogen using Gibbs Ensemble Monte Carlo (GEMC) simulation techniques with our new ab initio intermolecular pair potentials resulting from quantum mechanical calculations at a sufficiently high level of approximation of dimer H2-H2 [8]. The phase equilibrium results density, vapour pressure and enthalpy of vaporization obtained from GEMC simulation are compared with experimental data and with those from literature data. II - COMPUTATIONAL TECHNIQUE 1. Simulation details The Gibbs ensemble Monte Carlo (GEMC) simulation techniques were used for calculating the thermodynamic properties of the pure fluid 530 hydrogen. The GEMC-NPT simulation was used to calculate the density, and the internal energy of the fluid hydrogen to examine the accuracy of the pair potentials. This simulation was investigated on isobars at 1.0 MPa and 5.0 MPa and for temperatures from 26.0 K to 250 K, respectively. GEMC-NVT simulations were performed to obtain coexisting liquid and vapor densities, and vapor pressures. They were in the temperature range 18.0 K to 32.0 K with an increment 2.0 K. The two our new 5-site pair potential functions for hydrogen resulting from ab initio calculations [8] were used for both simulation cases: ∑∑ ∑ = = = − ++= 5 1 5 1 10,8,6 0 211 ]4 )()([ i i n ij ji ijn ij ij n ija rij e r qq rf r CrfeDu ijij πε α (1) ∑∑ ∑ = = = − ++= 5 1 5 1 12,10,8,6 0 212 ]4 )()([ i i n ij ji ijn ij ij n ijb rij e r qq rf r CrfeDu ijij πε α (2) With: 15 )2(2 1 )1()( −−−+= ijij rija erf δ , ∑ = −−= 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 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). The optimal adjustable parameters of these functions were estimated by nonlinear least-square fitting to the ab initio interaction energy values, as shown in [8]. Total number of particles N = 512 were used in both GEMC simulations with the standard periodic boundary conditions and the minimum image convention. For GEMC-NVT simulation runs the equilibration between two phase required 1-2 x 106 cycles. The simulation parameters were set for 50% acceptance ratios for translations and volume fluctuations. All movements were performed randomly with defined probabilities. The accumulative averages of desired quantities were established within 103 cycles, after initial equilibration had been reached within 5.0 x 104 cycles. The simulation data were exported using block averages with 1000 cycles per block. The statistical errors in the simulation runs were estimated by dividing each run into 100 blocks and taking the largest deviation of a block mean from the total mean as error. The simulations were started with equal densities in two phases. The simulation systems were equilibrated for about 1.0 x 106 cycles. The cut-off radius rc was set to 7.5  for hydrogen. Corrections for long- range interactions for the internal energy were computed by the standard relations [10]. 2. Structural properties The structural properties of the fluid hydrogen were studied for the liquid phase at different temperatures with the GEMC-NVT and -NPT simulations, respectively; in both cases the temperature dependence is shown by site-site pair distribution function gr. Because of the 5-site model of dimer H2-H2 was constructed with two sites placed on the atoms H, one site in the center of gravity M, and two sites halfways between the atoms and the center N [8]. Consequently, the 5-site pair correlation functions also consisted of the interactions H-H, N-N, M-M, N-M, H-N and H-M for the fluid hydrogen. These were achieved by simulations. The structural properties of fluid were compared with experimental data and with data from literature, if available. 3. Phase coexistence properties The critical temperature Tc/K, density ρc/g.cm-3 and volume Vc/cm3.mol-1 of the pure 531 fluid hydrogen were derived from least-squares fits to the densities of coexisting phase using the relations of the rectilinear diameter law [10]: βρρ ρρρ )( )( 2 21 21 c cc TTB TTA −=− −+=− (3) where ρl and ρv are the coexistence liquid density and vapor density, β is the critical exponent (β = 0.325). A and B are adjustable constants. The critical pressure Pc/ MPa was calculated with the Antoine equation CT BAP −−=ln (4) where A, B and C are Antoine constants. The relation between vapor pressure, heat of vaporization ΔvapH and temperature is given by the Clausius-Clapeyron equation: ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +Δ−= 212 1 11ln TTR H P P vap (5) For the standard state P0 = 0.101 MPa this relation is rewritten as R S TR H P P vapvap Δ+Δ−= 1ln 0 (6) Here T1 and T2 are the temperatures at the pressure P1 and P2. The slope and the intercept of lnP with respect to 1/T are proportional to ΔvapH and ΔvapS. III - RESULTS AND DISCUSSION 1. Structural properties This section describes the features of the site-site pair distribution functions resulting from two GEMC-NVT and -NPT simulation techniques for the pure fluid hydrogen. The ab initio pair potentials Eqs. 1 and 2 of hydrogen, respectively, were used for those simulations. The temperature dependence of the radial distribution functions at two different pressures 1.0 MPa and 5.0 MPa is obtained from the global simulations GEMC-NPT and GEMC- NVT. Figures 1 and 2 show simulation results for hydrogen. As 5-site potential models Eq.1 and Eq. 2 were used for these simulations, they consist of the interactions of ghost sites N, M on the molecules, so these were also represented here by 5-site pair correlation functions [8]. The height of site-site pair correlation functions decreased with increasing temperature from 26.0 K to 250.0 K, as shown in Table 1. In general the peaks for the interaction of sites including an atomic nucleus were higher than those without a nucleus. a) b) Figure 1: Temperature dependence of gH-H for hydrogen at P = 1.0 MPa a) simulation GEMC-NPT and b) simulation GEMC-NVT, in both cases using pair potential Eq.1 2 3 4 5 6 7 0 1 2 3 r/ Å g ( H -H ) T=26K T=30K T=60K T=90K T=120K T=250K 2 3 4 5 6 7 0 1 2 3 4 r/ Å g ( H -H ) T=18K T=20K T=22K T=24K T=26K T=28K T=30K T=32K 532 All first peaks of the site-site correlation functions for hydrogen are located between 2.893 Å and 3.205 Å. The second peaks are located between 6.081 Å and 6.234 Å. Table 1: The height of first peaks of the site-site distribution functions derived with GEMC-NPT simulation using the ab initio pair potentials Eq. 1 and Eq. 2 at 1.0 MPa at 5.0 MPa Equation T/ K gH-H gN-N gM-M gH-M gH-H gN-N gM-M gH-M 26.0 2.92 1.97 2.53 2.69 2.79 1.85 2.40 2.56 30.0 2.69 1.87 2.37 2.51 2.55 1.76 2.24 2.39 60.0 1.75 1.51 1.67 1.71 1.67 1.44 1.59 1.62 Eq. 1 90.0 1.49 1.35 1.44 1.47 1.42 1.28 1.37 1.40 26.0 2.88 1.94 2.51 2.67 2.74 1.84 2.37 2.53 30.0 2.79 1.90 2.44 2.59 2.67 1.80 2.32 2.46 60.0 2.41 1.74 2.15 2.27 2.27 1.65 2.04 2.14 Eq. 2 90.0 2.20 1.66 2.01 2.09 2.11 1.57 1.91 1.99 2. Phase coexistence properties The simulation results are shown in tables 2 and 3. The vapor-liquid coexisting phase curve of the fluid hydrogen is illustrated in figure 2. Experimental data [12, 13], values from the modified empirical equation of state [11] as well as from Johnson simulation results using the Goldman (SG) potential [14] are also included. The vapor pressures and enthalpies of vaporization derived from the same simulations are depicted in figures 3a and 3b. These vapor pressures differ on absolute average from the experimental data typically by about 3.49% and 9.14%. These differences are small within statistical uncertainties of experimental resources and a few previous publications [11, 14]. Figure 2: Vapor-liquid coexistence diagram of hydrogen. •, experimental data [12, 13]; o, modified BWR equation [11]; ◊, simulated by Johnson [14]; ―, --- : ab initio pair potentials Eq. 1 and Eq. 2. a) b) Figure 3: a) Vapor pressure, b) Vaporization enthalpy; for an explanation see Fig. 2 18 20 22 24 26 28 30 32 34 0.0 0.2 0.4 0.6 0.8 1.0 1.2 T/K p/ M Pa 18 20 22 24 26 28 30 32 34 0.0 0.2 0.4 0.6 0.8 1.0 T/K ΔH va p/K J m ol -1 0.00 0.02 0.04 0.06 0.08 18 21 24 27 30 33 T/ K ρ/g cm-3 533 Table 2: Critical properties of hydrogen resulting from the GEMC-NVT simulation using equations Eq. 1 and 2; EOS: equation of state [11]; Exp.: experimental values Method Tc, K ρc, g.cm-3 Pc, MPa Vc, cm3mol-1 Ref. Eq. 1 33.216 0.0313 1.1258 64.3806 this work Eq. 2 33.024 0.0311 1.0990 64.7907 this work EOS 32.972 0.0312 1.2837 64.1539 [11] Exp. 33.190 0.0312 1.2928 64.1026 [12] Exp. 33.00 0.0310 1.2930 64.5677 [13] The critical properties of the pure fluid hydrogen could not be calculated directly from the simulation, but they could be obtained from the orthobaric densities of vapor-liquid equilibria by the least-square fit to the relations (3), as shown in table 2. The critical pressure of hydrogen was calculated from the Antoine equation Eq. 4. The results agreed reasonable well with experimental data. The thermodynamic properties of this fluid are also shown in table 3. Table 3: Enthalpy of vaporization ΔvapH, entropy of vaporization ΔvapS and boiling temperature Tb at P = 101.3 kPa predicted from simulation vapor pressures Method ΔvapH, kJ mol-1 ΔvapS, kJ.mol-1.K-1 Tb, K ref. Eq. 1 1.17148 0.05608 20.8911 this work Eq. 2 1.21621 0.05717 21.2740 this work EOS 1.07399 0.05305 20.2457 [11] Exp. 1.07786 0.05299 20.3900 [12] Exp. 1.07752 0.05314 20.2754 [13] The discrepancies between predicted results and experimental data are insignificant. IV - CONCLUSION The vapour-liquid phase equilibria and thermodynamic properties of the pure fluid hydrogen were calculated successfully with our developed computer simulation programs GEMC-NVT and GEMC-NPT using ab initio intermolecular pair potentials. The simulation results turn out to be in very good agreement with experimental data and with those from literature data. This also confirmed that the two our new 5-site analytical potential functions resulting from ab initio calculations in the work [8] are of high quality, accurate and reliable. Acknowledgments: The Regional Computer Center of Cologne (RRZK) contributed to this project by a generous allowance of computer time. We would like to thank Dr. Naicker, Prof. A. K. Sum and Prof. S. I. Sandler (University of Delaware, USA) for making available their computer programs. References 1. L. J. Krzycki, build and test small liquid- fuel rocket engines, Oregon (1967). 2. V. P. Dawson. A History of the Rocket Engine Test Facility at the NASA Glenn Research Center, NASA, Glen Research Center, Cleveland (2004). 3. E. Nasrabad and U. K. Deiters. J. Chem. Phys., 119, 947 - 952 (2003). 4. A. E. Nasrabad, R. Laghaei, and U. K. Deiters. J. Chem. Phys., 121, 6423 (2004). 5. K. Leonhard and U. K. Deiters. Mol. Phys., 534 100, 2571 - 2585 (2002). 6. P. K. Naicker, A. K. Sum, and S. I. Sandler. J. Chem. Phys., 118, 4086 - 4093 (2003). 7. Z. Panagiotopoulos. Mol. Phys., 61, 813 - 826 (1987). 7. Pham Van Tat, U. K. Deiters. Journal of Chemistry, Vol. 5(5), 656 - 660 (2007). 8. R. A. Kendall, T. H. Dunning, Jr., J. Chem. Phys., 96, 6796 - 6806 (1992). 9. M. P. Allen and D. J. Tildesley. Computer Simulation of Liquids, Clarendon Press, Oxford, (1991). 10. A. Younglove. J. Phys. Chem. Ref. Data, Vol. 11(1), 1 - 161 (1980). 11. R. D. McCarty and L. A. Weber. Thermophysical properties of parahydrogen, Government Printing Office, Washington, D.C, 1972. 12. R. Lide. Handbook of Chemistry and Physics, CRC Press, 82nd Edition, Boca Raton (2002). 13. Q. Wang and J. K. Johnson and J. Q. Broughton. Mol. Phys., Vol. 89, 1105 (1996).