Ab initio potential energy surface for methane and carbon dioxide and application to vapor-liquid coexistence Sung Jin Pai and Young Chan Bae Citation: The Journal of Chemical Physics 141, 064303 (2014); doi: 10.1063/1.4891983 View online: http://dx.doi.org/10.1063/1.4891983 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/6?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Communication: A new ab initio potential energy surface for HCl–H2O, diffusion Monte Carlo calculations of D 0 and a delocalized zero-point wavefunction J. Chem. Phys. 138, 121102 (2013); 10.1063/1.4799231 Ab initio ground state phenylacetylene–argon intermolecular potential energy surface and rovibrational spectrum J. Chem. Phys. 137, 074305 (2012); 10.1063/1.4742153 Ab initio potential energy surface and intermolecular vibrations of the naphthalene-argon van der Waals complex J. Chem. Phys. 134, 064322 (2011); 10.1063/1.3555765 Spectra of Ar–CO 2 from ab initio potential energy surfaces J. Chem. Phys. 112, 5308 (2000); 10.1063/1.481120 Structural and thermodynamic properties of fluid carbon dioxide from a new ab initio potential energy surface J. Chem. Phys. 109, 3153 (1998); 10.1063/1.476922

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

THE JOURNAL OF CHEMICAL PHYSICS 141, 064303 (2014)

Ab initio potential energy surface for methane and carbon dioxide and application to vapor-liquid coexistence Sung Jin Pai and Young Chan Baea) Division of Chemical Engineering and Molecular Thermodynamics Laboratory, Hanyang University, Seoul 133-791, South Korea

(Received 16 April 2014; accepted 22 July 2014; published online 8 August 2014) A six-dimensional intermolecular potential energy surface for a rigid methane (CH4 ) and carbon dioxide (CO2 ) dimer was developed from the counterpoise-corrected supermolecular approach at the CCSD(T) level of theory. A total of 466 grid points distributed to 46 orientations were calculated from the complete basis set limit extrapolation based on up to aug-cc-pVQZ basis set. A modified site-site pair potential function was proposed for rapid representation of the high level ab initio calculations. A nonadditive three-body interaction was represented by the Axilrod-TellerMuto expression for mixtures with the polarizability and the London dispersion constant of each molecule. Second to fourth virial coefficients of CH4 and CO2 mixtures were calculated using both the Mayer sampling Monte Carlo method and the present potential functions. The virial equation of state derived from these coefficients was used to predict the pVT values and showed good agreement with experimental data below 200 bar at 300 K. The vapor-liquid coexistence curves of pure CH4 , CO2 and their mixtures were presented with the aid of Gibbs ensemble Monte Carlo simulations. The predicted tie lines agreed with the experimental data within the uncertainties up to near the critical point. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4891983] I. INTRODUCTION

Hydrocarbons are one of the most important and ubiquitous materials in all of the industry and the development of shale gas techniques has incited an increase in their usage in recent years.1 Natural or shale gases vary in their composition and the required processing for these materials changes from area to area. In the constituents of the various gases, methane (CH4 ) is the major component (80% or more); therefore, the state of the gas systems seems to be determined by the properties of pure methane. However, small amount of contaminants such as water, hydrogen sulfide (H2 S), or carbon dioxide (CO2 ) can considerably change the gas properties depending on the degree of contamination. Increasing the hydrocarbon ratio by removing those contaminants is an important operation in natural gas processing to meet pipeline quality standards. Meanwhile, it is known that as the portions of CH4 and CO2 in gas system increase (to a certain level), the condensation risk decreases.2 Therefore, an adequate amount of CO2 is needed to safely transfer the gases. To find the optimum compromise between these variables, a thorough understanding about the interactions between the components and their effects on the properties of the entire system are needed. There have been many attempts to describe experiments for the physical state of gases or liquids through thermodynamic approaches. A large amount of equations of state (EOS) have been developed by researchers, but most of these must be adjusted with several parameters in order to agree with experiments. Most of these problems are caused by a) Author to whom correspondence should be addressed. Electronic mail:

[email protected]. Tel.: +82-2-2220-0529. Fax: +82-2-2296-0568. URL: http://www.labmtl.com/.

0021-9606/2014/141(6)/064303/9/$30.00

the fact that these parameters are derived through statistical mechanics, and the magnitude of a parameter depends on the essence of a material itself. The interactions between atoms or molecules, which are integrated to be macroscopic properties, have been studied by calculating the intermolecular potentials using quantum mechanics. With dramatic improvement of hardware, highly accurate calculations such as weakly non-covalent interactions between small molecules have been made possible.3 Recently, ab initio calculations about the intermolecular potential energy surface (PES) of small molecules such as CH4 or CO2 have been carried out using very high level of theories. Hellmann et al. determined a six-dimensional PES for a CH4 dimer from the ab initio calculation using the counterpoise-corrected supermolecular approach at the couple-cluster method with single, double and noniterative triple excitations (CCSD(T)) level of theory with complete basis set extrapolation.4 Although they did use a semiempirical correction, the developed site-site potential function agreed well with the experimental second virial coefficients over the whole examined range. Ab initio calculations for CO2 have been conducted earlier due to the advantage of their simple symmetry. Bukowski et al. computed a four-dimensional PES for a CO2 dimer using manybody symmetry-adapted perturbation theory (SAPT) which is nearly identical to the second-order supermolecular manybody perturbation theory (MP2).5 Bock et al. presented another four-dimensional PES for a CO2 dimer based on the supermolecular approach of MP2 level of theory including full counterpoise correction and complete basis set extrapolation.6 With the developed site-site potential function, they were also able to derive many physical properties such as viscosity, self-diffusion, and thermal conductivity.7, 8 Bratschi et al.

141, 064303-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

064303-2

S. J. Pai and Y. C. Bae

calculated the vapor-liquid coexistence of a pure CO2 system using both the site-site potential functions of Bukowski et al. and Bock et al. using Gibbs ensemble molecular dynamics (GEMD).9 They concluded that the potential function of the latter is better at describing the experimental critical point compared to the former. However, Bratschi et al. did not consider the nonadditive three-body interaction effect in their work, which is known to have contributions of about 10%20% of two-body interactions and shrink the region of vapor-liquid equilibrium (VLE).10 Had taken this into consideration, the quality of the potential function of Bukowski et al. would likely have been better. To describe macroscopic properties explicitly, a nonadditive three-body potential function and pair interactions are needed. These can be found via three-body ab initio calculations. Oakley et al. calculated ab initio calculations of CH4 and CO2 dimers and trimers, from which they developed pair and nonadditive three-body potential functions and carried out comparable Gibbs ensemble Monte Carlo (GEMC) simulations with experimental data.11, 12 However, they did not present a comparison with the experimental second virial coefficients, which is generally accepted as a necessary criterion to prove the quality of a potential function. Moreover, there were increasing deviations as the critical point was approached. These deviations were especially large for the binary mixture, and may have been caused by the low level of ab initio calculations (MP2) and the approximated accounting of the nonadditive three-body interactions. In this work, we show a quantitative description of the vapor-liquid coexistence of pure CH4 , CO2 , and their mixtures through high level ab initio calculations of dimers and explicit nonadditive interactions. Supermolecular CCSD(T) level of theory calculations for the CH4 –CO2 dimer are carried out with basis sets up to aug-cc-pVQZ quality13, 14 for 466 points distributed to 46 orientations. A relatively simple site-site pair potential function for the CH4 –CO2 dimer is proposed for rapid molecular simulations and fitted to the ab initio calculations. For pure CH4 , the results of Hellmann et al. are refitted to the new potential function, while the pair potential function of Bukowski et al. is used as is for the CO2 dimer. The Axilrod-Teller-Muto (ATM) expression15, 16 is used for the nonadditive three-body interaction using only the polarizability and London dispersion constant for each molecule. The second to fourth virial coefficients of pure CH4 , CO2 and their mixtures are calculated using the Mayer sampling Monte Carlo (MSMC) method,17, 18 and the mixture virial equation of state (VEOS) made from the coefficients is compared with experimental data. Finally, GEMC simulations for the pure CH4 , CO2 and their mixtures are carried out using the developed pair and nonadditive three-body potential functions and compared with both the experimental data and the work of Oakley et al.

II. INTERMOLECULAR POTENTIAL SURFACE A. Ab initio calculation for the CH4 –CO2 dimer

To calculate the intermolecular potential energies, the geometries of the molecules must be identified in advance. We

J. Chem. Phys. 141, 064303 (2014)

FIG. 1. CH4 –CO2 dimer geometry.

benchmarked the literatures and decided to follow the rigid molecule assumption and use the recommended values. For CH4 , the bond length of the C–H bonds were fixed to 1.099 Å (experimental zero-point vibrational average value) and the structure was set to be a regular tetrahedron.4 The linear CO2 molecule was given the bond length of the C–O bonds (1.162047 Å, from the experimental rotational constant).5 Due to the different symmetry groups belonging to each molecule, a mixed, six-dimensional polar-Euler angle coordinate system (R12 , θ 12 , ϕ 12 , ϕ 2 , θ 2 , and ψ 2 ) was selected as the calculation grid. This is the same calculation grid that was used in the work of Akin-Ojo and Szalewicz for a methanewater dimer. R12 is the carbon-to-carbon distance; θ 12 and ϕ 12 are the polar angles; and ϕ 2 , θ 2 , and ψ 2 are the Euler angles with x convention (Fig. 1).19 With the assistance of the symmetry, the calculation space was restricted: θ 12 ∈ [0, 90], ϕ 12 ∈ [0, 90], θ 2 ∈ [0, 90], and ϕ 2 = ψ 2 = 0. Starting at (R12 , 0, 0, 0, 0, 0), the polar angles θ 12 and ϕ 12 or the Euler angle θ 2 were changed by 30◦ one after another, and some coordinates near the minimum energy paths were additionally included. Such predetermined combinations of the angles in the symmetry range of 23 orientations and random combinations of 23 orientations were chosen for the ab initio calculation. For each orientation, relatively cheap trials with CCSD(T)/aug-cc-pVDZ were calculated at between 2.5 and 10 Å. Then valuable points were selected and used for higher level calculations. The total number of calculated points was 466. The selected points were calculated using the supermoecular approach with a counterpoise correction using CCSD(T) level theory13 and aug-cc-pVXZ (X = 2, 3, 4) basis sets. The calculated energies were divided into two parts: the HartreeFock (HF) energy and the correlation energy which was the difference between the CCSD(T) and HF energies. The correlation energies were extrapolated to the CBS limit using the equation of Halkier et al.,20 X CBS Vcorr = Vcorr + αX−3 ,

(1)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

064303-3

S. J. Pai and Y. C. Bae

J. Chem. Phys. 141, 064303 (2014)

the two molecules. The upper part of Eq. (2) has a maximum at rmax and then goes down to –∞ as two particles approach each other. This unphysical problem was a drawback of exp-6 potential type functions, and later modified by including the lower part of Eq. (2).23 The Coulomb term was not used in order to exclude long range corrections for molecular simulations. The anisotropic dependence of each coefficient can be described with various routes. In this work, we propose a model in which there are sites with fixed centers and the coefficients are calculated by auxiliary equations such as  a kl , (3) A() = r kl k l FIG. 2. Ab initio calculations (symbols) and calculations from the pair potential function (lines) of the CH4 –CO2 dimer for selected orientations.

β() =

 k

where X is the basis set order, and α is a correlation constant. The HF energy was assumed to be saturated at the aug-ccpVQZ basis set so that extrapolation was not necessary. No more corrections were added due to the relatively small contributions. All ab initio calculations were carried out using Gaussian03.21 The calculated ab initio values are presented in the supplementary material.22 Some of the results are represented in Fig. 2 as symbols. We can see that the potential energy predominately depends on the orientation of the CO2 . Also the well depth is deepest when the symmetry axis of CO2 is orthogonal to the normal vector of the face of CH4 (orientation 1). The center-to-center distance at the minimum energy increases as the acute angle between the symmetry axis of CO2 and the normal vector of the CH4 face decreases; this is expected. For similar orientations, the energy is relatively low when the hydrogen atom of CH4 approaches the oxygen atom of CO2 . This can be explained by the opposite signs of the partial charges in the ab initio calculations. To examine the effect of the level of the ab initio method, MP2 level calculations were also conducted with the same basis sets. Comparing the two level of theories shows that when two molecules are close to one another (within 5 Å), the CBS limit energies of the CCSD(T) level are generally 2%–3% lower than those of the MP2 level. However, beyond 5 Å, this condition is reversed. B. Modified site-site pair potential function

In many works concerning the PES, site-site potential functions were fitted to their calculated values, most of which were computationally expensive. For the purpose of brevity and speediness, while maintaining anisotropic characteristics, we propose a new potential function in which the distance and orientation are separated and the orientational information is contained only in the coefficients. The form of the pair potential function is  C () C () eAij ()−βij ()Rij − 6,ij − 8,ij when Rij > rmax (2) Rij6 Rij8 , Vij = ∞ when Rij ≤ rmax (2) where A, β, C6 , and C8 are material dependent coefficients, Rij is the center-to-center distance between molecules i and j, and  represents the polar-Euler angle coordinate between

C8 () =

bkl rkl ,

(4)

l

  c8,kl k

l

rkl8

,

(5)

where akl , bkl , and c8,kl are fitting parameters for sites k and l and rkl is the distance between two sites from each molecule. The forms of Eqs. (3)–(5) were motivated by the detailed expressions of each coefficient in Eq. (2) to reflect the change in orientation and constraint the digits of the parameters to the total summation. Because CH4 and CO2 are not polar molecules, C6 was set to be isotropic in this work: C6,CH4 is 898 647 K/Å6 and C6,CH4CO2 is 978 605 K/Å6 . The former value was taken from the fitted values of the original potential function of Hellmann et al.,4 and the latter is from the geometric mean of the former and C6,CO2 (1 065 677 K/Å6 ) which was used in the original function of Bukowski et al.5 The geometry of the molecules was used for the sites: five atom positions for CH4 and three atom positions plus two virtual sites for CO2 . With different to the original potential function,4 there is no virtual site in CH4 molecule to reduce the calculation time. In practice, adding virtual sites in CH4 did not improve much the fitting quality. A remaining problem concerned the fixed distance between the two centers. Care should be taken in saying the fixed distance, which is not the real distance between two molecules but a hypothetically fixed center-to-center distance for all real orientations and distances. It is just a parameter to reflect the orientational dependency of the coefficients in Eq. (2). The recommended values for these in this work were 3 Å for the CH4 dimer and 4.5 Å for the CH4 –CO2 dimer. For the CH4 dimer, Eq. (2) was fitted to the calculated values using the original potential function for 100 random orientations, and for the CH4 – CO2 dimer, all ab initio values calculated from the present work were used to fit with the same equation. Obtained parameters using the above procedure for the CH4 dimer and the CH4 –CO2 dimer are presented in Table I. For the CO2 dimer, Eqs. (3)–(5) were not appropriate for fitting the potential energies so we needed to use another scheme. However, in this work, we decided to use the original potential function due to its relatively inexpensive computation. Fig. 2 shows the fitting results for the selected orientations of the CH4 –CO2 dimer with the ab initio values. The global minimum estimated from the fitted function is at the grid point (3.379, 87.573, 14.864,

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

064303-4

S. J. Pai and Y. C. Bae

J. Chem. Phys. 141, 064303 (2014)

TABLE I. Potential function parameters for the CH4 dimer and the CH4 –CO2 dimer. Value in parenthesis is the exponent of 10.

Species

Sites

akl ...

bkl Å−1

c8,kl KÅ−8

CH4 –CH4

C–C C–H H–H

3.7145329024(1) − 8.7452212875(0) 5.6003501771(0)

− 1.5861561410(1) 5.2082487619(0) − 1.5242632762(0)

3.5250575170(10) 2.8912515009(9) − 5.3696688452(6)

CH4 –CO2

C–C C–O C–Ma H–C H–O H–M

2.9237682704(2) 8.3504655115(1) − 2.4159497426(2) 4.0665843676(1) 7.5613544617(1) − 8.1962982441(1)

4.0611479499(2) 2.2116396338(2) − 4.3062701931(2) − 1.5993791796(2) − 8.9233087166(1) 1.7087537835(2)

7.1326882651(12) 1.2556382913(12) − 3.9497845150(12) − 1.3406524526(12) − 1.8339969482(11) 8.6317073279(11)

a

M represents the virtual site in CO2 molecule.

−21.203, 61.054, 101.383) with V12 of −510.086 K, but ab initio calculation shows that the potential energy at this point is −472.400 K. The true global minimum estimated from the geometry optimization with MP2/aug-cc-pVTZ is at the grid point (3.383, 54.251, 88.376, 276.433, 123.351, 309.236) with V12 of −491.787 K, which is located at a point slightly rotated from the minimum point of orientation 1 in Fig. 1. The fitting quality is the total root mean squared error of 121.294 K (0.284 kcal/mol), which is a fair value considering the relatively simple equations. In Fig. 3, the ab initio calculation of Oakley et al.12 (which is the only such calculation we found in the literature) is compared with the present potential function for two selected orientations. Even if they used the MP2 level of theory for their calculations, their best potential function, ISA-ad, underestimated the ab initio data of the present work by about 40% at the potential minimum. Fig. 3 also shows that the two potential calculations converged beyond 6 Å. For other orientations, the underestimation of the ISA-ad near the potential minimum is general. This should have a significant effect on the statistical properties (which are presented below in Sec. IV).

work for most molecules (other than noble gases) due to the enormous combination of orientations. Bukowski et al. analyzed the nonadditive contributions to the three-body interactions of Ar trimers and showed that the triple-dipole contribution represents the total three-body energies (although they were caused by error cancellations).10 With the assumption that the triple-dipole term (ATM term) is representative of nonpolar molecules, this term is solely used for the nonadditive three-body interactions in this work. Midzuno and Kihara derived the ATM equation for mixture with the variation method as25 (1 + cos θi cos θj cos θk ) Vij(3)k = νij k , (6) (Rij Rj k Rik )3

C. Nonadditive three-body potential function

θ i is the inner angle at vertex i in the triangle ijk, μij is the London dispersion constant between i and j, and α i is the polarizability of i. Because Eq. (6) was derived for three atoms, it is assumed that the point molecules are at their centers of mass in the calculation of the three-body interaction. This approach also has the added advantage that we can use C6,ij (from Eq. (2)) as the London dispersion constant, μij , without any additional adjustable parameters. In other words, the nonadditive three-body interaction for trimers can be estimated from Eqs. (2) and (6) without any ab initio calculations. Although the exchange contribution is not negligible in most cases, the magnitude of the three-body interaction energy is generally 5% of the two-body interaction energy.10 We anticipate that this error will be negligible as it is in the case of molecular simulations with massive objects.

In SAPT, the nonadditive three-body interaction energy is composed of various perturbation terms.24 Strict application of the exact functions will give correct results, but the extremely long computation time limits the usefulness of this approach. Additionally, it is necessary to calculate many ab initio calculations for trimers, which is a large amount

where νij k =

2 i j k ( i + j + k ) ( i + j )( j + k )( k + i )

and 1 1 1 1 = + − i μij αk μki αj μj k αi

etc.

III. VIRIAL COEFFICIENTS AND EQUATION OF STATE A. Theory and calculation FIG. 3. Comparison with other PES for the CH4 –CO2 dimer. Selected orientations are 1 and 6 in Fig. 2.

The classical pressure virial coefficients can be obtained from statistical mechanical relations for the cluster integrals

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

064303-5

S. J. Pai and Y. C. Bae

J. Chem. Phys. 141, 064303 (2014)

with the developed potential functions. They can be expressed as linear graphs and the first three coefficients are26–29 (7)

(8)

(9) where B, C, and D are second, third, and fourth virial coefficients, respectively. Lines and shaded triangles are Mayer functions: the solid line is fij = exp (–Vij /kT) – 1, the dotted line is eij = exp (–Vij /kT), and the shaded triangle is fijk = exp (–Vijk /kT) – 1, where k is the Boltzmann constant and T is the temperature. Each vertex represents the identity of a molecule and each figure implies integration with respect to the configurational space of the molecules. For example, Eq. (7) means   dr12 d12 f12 1 1  =− (10) dr12 f12  , B=− 12 2 2 d12 where the angular bracket is the average about the subscript. For monatomic gases such as Ar, the integrations are manageable up to the fifth order coefficient.30 However, implementation of the integrations becomes arduous if the particles are diatomic, and may be impossible for polyatomic gases by normal quadratures. Singh and Kofke introduced a variant of a biased Monte Carlo quadrature, referred to as the MSMC method, for the calculation of higher virial coefficients.17 This has been successful when applied to rigid polyatomic gases such as water or nitrogen.18, 31 Because MSMC is a free energy perturbation method, there are many different versions depending on the sampling method. To calculate the second to fourth order virial coefficients for CH4 , CO2 , and their mixtures, we adapted the overlap sampling MSMC. This method considers the important region of the target and reference separately and correlates them through an overlap function for a more accurate result.18 The nth order virial coefficient is represented as a ratio of the real potential and the hard sphere potential systems as   Bnre /π re π re / BnOS /π re π re hs    ,  (11) Bn (T ) = Bn hs hs Bn /π π hs / BnOS /π hs π hs

placed with molecules of the corresponding number for each virial order. For each MC trials, one translation or one rotation of a single particle was carried out randomly. A 4 × 105 equilibration times were processed prior to the sampling in order to adjust the acceptance ratio so that it achieved 50%. Then, 4 × 108 production times were carried out, of which 1% were for the reference sampling configurations. Every 105 steps were blocked to estimate the standard errors and correlations in the sampling. For a three-body interaction, a hard sphere diameter of 1 Å was applied to correct the unphysical region. The hard sphere diameter of the reference system was set to 5 Å, which was roughly estimated from the global minimum distance multiplied by 1.5. Each calculation is an average of 10 separate trials. In our previous work, we obtained analytical expressions for second to fourth virial coefficients using the semi-soft core (SSC) potential function.29 The obtained MSMC calculation points were fitted to those analytical expr essions for each virial coefficient. Used equations and the obtained fitting parameters are presented in the supplementary material.22 Using these analytical expressions, VEOS can be made for CH4 and CO2 mixtures as p = 1 + Bmix (T )ρ + Cmix (T )ρ 2 + Dmix (T )ρ 3 , (14) ρkT where p is pressure, ρ is density, and  Bmix (T ) = xi xj Bij (T ), Cmix (T ) =

i

Dmix (T ) =

i

j

j

k



xi xj xk Cij k (T ),

 i

j

k

xi xj xk xl Dij kl (T ),

l

where xi is the mole fraction of particle i. We named the EOS truncated beyond the third or fourth virial coefficients as VEOS3 and VEOS4, respectively. B. Comparisons with experimental data

To assess the quality of the ab initio calculations and the obtained potential function, the cross second virial coefficient of the CH4 –CO2 dimer was calculated using the MSMC method. This is shown in Fig. 4 for the temperature range between 280 K and 330 K and compared with the literature. It

where π is the configurational space for the importance sampling region. The superscripts re and hs refer to the real potential system and the hard sphere as a reference system, respectively. The overlap function Bn OS is defined as  hs  re Bn  |Bn | OS   Bn = , (12) α Bnhs  + |Bnre | where α is an optimization parameter which can be optimized during sampling with the criterion    OS hs  (13) Bn /π π hs = α BnOS /π re π re . Averages were obtained by the general Metropolis algorithm. Two boxes, one for the target and the other for reference, were

FIG. 4. Cross second virial coefficient for the CH4 –CO2 dimer. Symbols are experimental data: , Martin et al.;38 ●, Esper et al.;32 , Ohgaki et al.;39 , Hou et al.;33 ♦, McElroy et al.;40 , Reamer et al.;41 , Mallu et al.42

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

064303-6

S. J. Pai and Y. C. Bae

FIG. 5. Third virial coefficient of Ar. Experimental data are from Jäger et al.34 modifying the data of Gilgen et al.43 The dotted line is the calculation result using the full three-body potential function by Jäger et al. and the solid line is the calculation result of the present work.

can be seen from ure that the data of Esper et al.32 fit perfectly with the present potential function. The most recent data of Hou et al.33 also agree well with our present work. Other data appear to be overestimated by the present work, but the uncertainty ranges of them include the calculations (solid line). The calculated values of the potential function of Oakley et al. are −12.120 cm3 /mol for 300 K and −7.440 cm3 /mol for 320 K, which underestimates the experimental data (outside of the range of the figure). The shallow well depth of the potential energy used by Oakley et al. may be the cause of this underestimation. To prove the representativeness of the ATM term for the nonadditive three-body potential energy, the third virial coefficient of Ar is calculated using Eqs. (6) and (8) with the pair potential function of Jäger et al.34 The polarizability of Ar35 is 1.630 Å3 and the London dispersion constant is μij = C6,ij = 442 812.017 K/Å6 . The results are shown in Fig. 5 with the experimental data and the calculation values of Jäger et al. using the full three-body potential function from the ab initio calculation. When compared with the full potential function of Jäger et al. (dotted line), the third virial coefficients calculated from the ATM-only energies (solid line) are underestimated near the maximum point, but converge to the full potential function beyond 200 K. This shows that nonadditive terms other than the ATM term are significant near the critical point, but as temperature increases, the triple dipole term surpasses the other terms to become the major contributor or error cancellations among the other terms become true.

J. Chem. Phys. 141, 064303 (2014)

If an equation of state for a substance requires a very accurate prediction near the supercritical condition, it is necessary to carry out the high level ab initio calculations for the dimer and trimer. However, if the demand for accuracy is not severe, the prediction of the physical state could be accomplished with only ab initio calculations for dimers. With confidence in the applicability of the present work, the third and fourth virial coefficients for pure CH4 (component 1), CO2 (component 2), and their mixture were calculated using Eqs. (2) and (6). The results of the MSMC calculation for the third virial coefficients are represented in Fig. 6(a) with the experimental data. Polarizabilities are 2.593 and 2.911Å3 for CH4 and CO2 , respectively.35 As the portion of CO2 molecules increases, the maximum point moves to a higher temperature and the magnitude at the high temperature becomes larger. It can be seen from the figure that the calculation results agree well with the experimental data at high temperatures. However, unlike what was seen in the case of Ar, there are noticeable deviations near the maximum points. The deviation becomes severe from CH4 to CO2 , which may be caused by the inappropriate dimer ab initio calculation (MP2 level) or neglecting the quadrupole effect in the threebody nonadditive interactions of the CO2 molecules. The general reasons for this deviation may be multifaceted. First is about the geometry; rigid body assumption may be wrong near the critical point. Second concerns the nonadditive energy; the ATM term cannot represent the whole three-body interactions. Final candidate may be the deficiency of the quantum correction, but this is known to have little effect. The last two possibilities can be examined by additional works, but without ab initio calculation for the trimer and any adjustable parameters, the present degree of accuracy proves the feasibility of this work for applications in industry. In Fig. 6(b), the fourth virial coefficients are represented. Because there is a lack of experimental data for fourth virial coefficients and, even if they exist, the confidence in the data is low, we did not compare the calculations with the experiments. The movement of the maximum point is similar to what is observed for the third virial coefficient, but the order of the high temperature magnitudes is in reverse. It is noteworthy that the magnitude of 1112 converges to 1111 beyond 600 K, which is not seen in the case of the third virial coefficient. This may be a phenomenon caused from little difference between the CH4 –CH4 and CH4 –CO2 pair interactions in the tail region, in which the sampling is important at high temperatures. The

FIG. 6. (a) Third and (b) fourth virial coefficients of CH4 (1), CO2 (2) and their mixtures. Experimental data are:  and ◦ by Holleran;44  and  by Hou et al.33 Lines are the calculated results of the present work.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

064303-7

S. J. Pai and Y. C. Bae

FIG. 7. Phase diagram of CH4 (1), CO2 (2) and their mixtures at 300 K. Experimental data are:  and from NIST;45 ◦ by Esper et al.32 Solid lines are the calculated results of VEOS4 and dashed lines are those of VEOS3. High density points for pure CO2 ( > 0.01 mol/cm3 ) are the liquid phases.

reason that this phenomenon occurs only for the fourth virial coefficient may be that the number of CH4 –CH4 pair interactions is larger for the fourth than for the third virial coefficient. In light of this inference, convergence may occur in the third virial coefficient at an even higher temperature. Another point of interest concerns the magnitude of the maximum points, of which a trend comes from the two-body energies. The reason for the highest value at 1222 may be that small amounts of CH4 break the stabilized structure of the CO2 clusters and direct them to high energy orientations. However, the accuracy is not reliable enough near these points to prove this claim; it would be better not to mention this further in this work. The predicted phase diagram for the CH4 and CO2 mixtures at 300 K developed with the VEOSs are represented and compared with the experimental data in Fig. 7. For pure CH4 and the mixture, the prediction agrees well with the experimental data below 200 bar, and VEOS4 (solid lines) shows better results than VEOS3 (dashed lines) for the whole region. It may be inferred from this result that higher VEOSs would enlarge the region to even higher pressures. In the case of pure CO2 , however, the prediction agrees well only within the vapor region. At this temperature, the high density region (above 0.01 mol/cm3 ) shows the liquid density for pure CO2 , in which the deviation grows significantly from VEOS3 to VEOS4 according with the divergent character of VEOS in this region. Nevertheless, it can be seen from the result that the developed VEOS4 can be used as an accurate guide for predicting the phase diagrams for natural or shale gas systems.

J. Chem. Phys. 141, 064303 (2014)

FIG. 8. Vapor-liquid coexistence of pure CH4 . Solid line represents experimental data.46 Symbols are GEMC simulation points. ✩ is the predicted critical point from the present work.

A 50% acceptance ratio was maintained through the production steps. A cut-off radius, rc , of a minimum value between 15 Å and half the box length was used for both pair and nonadditive interactions. The closest countable distance for nonadditive interactions was set to be 1 Å. Long-range correction was applied only to van der Waals type pair interactions since a CO2 molecule has no charge or dipole moment and the nonadditive three-body tail correction energies are insignificant enough to ignore. The critical temperature and density were derived from a combination of the scaling raw ρl − ρg = B(T − Tc )β

(15)

and the law of rectilinear diameters ρl + ρg = ρc + A(T − Tc ), (16) 2 where the subscripts l, g, and c represent the liquid, gas, and critical point, respectively. β is the critical exponent set to 0.32 and A and B are proportionality constants. For mixtures, NpT GEMC simulations were performed because the pressures were fixed in advance in order to reduce calculation costs. For mixture simulations, we placed 200 molecules in each box using the experimental total molar ratios because the particle number of each species is important to draw the result to the experimental condition. Using the GEMC simulation method and the present pair and nonadditive three-body potential functions, vapor-liquid coexistence diagrams of pure CH4 , CO2 , and their mixtures were calculated and are presented in Figs. 8–10. All simulation results were compared with both the experimental

IV. VAPOR-LIQUID COEXISTENCE

The vapor-liquid coexistence of simple molecules can be described via many EOS or simulation techniques. Among the various approaches, GEMC simulation enables direct reduction of the phase coexistence properties from the present potential functions.36 For pure CH4 and CO2 , NVT GEMC simulations with constant volume were performed over the experimental range. Four hundred molecules, allotted in two boxes, were initially placed as the center atoms in lattice points with random orientations. 1 × 105 steps for equilibration and 2 × 105 steps for production were processed in each simulation. In a trial step, a random composite of 200 translations, 200 rotations, 40 exchanges, and 2 volume changes was conducted.

FIG. 9. Vapor-liquid coexistence of CO2 .46 All specifications are the same as in Fig. 8.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

064303-8

S. J. Pai and Y. C. Bae

J. Chem. Phys. 141, 064303 (2014)

FIG. 10. CH4 /CO2 pressure-composition diagram at (a) 230 K, (b) 250 and 270 K. Solid lines are experimental data37 and symbols are GEMC simulation points.

data37 and the GEMC simulations results of Oakley et al.12 For pure CH4 (Fig. 8), the results of the present work, including the nonadditivity three-body interaction (◦), perfectly agree with the experimental data. This is especially true near the critical point. When considering only pair interaction energies (), it can be seen that the predicted coexistence region is enlarged as expected. The estimated critical point (✩) showed ρ c = 0.01012 mol/cm3 and Tc = 189.713 K, which is very close to what was observed in the experimental data (ρ c exp = 0.01015 mol/cm3 and Tc exp = 190.6 K). The results of Oakley et al. deviate in the vapor region close to the critical point. Looking at their results, it is clear that the deviation caused by the incorrect well depth of their pair potential function becomes remarkable near the critical region. For pure CO2 (Fig. 9), the results were similar to the case of pure CH4 . Although the accuracy decreases somewhat near the critical point, the quality of the prediction is still better than the work of Oakley et al. The present work estimated the critical point (✩) of CO2 as ρ c = 0.01059 mol/cm3 and Tc = 305.579 K, which is comparable to the experimental data (ρ c exp = 0.01064 mol/cm3 and Tc exp = 304.19 K).

The importance of the correct pair potential energy was more striking in the case of mixtures. In Fig. 10, CH4 /CO2 pressure-composition diagrams are shown at 230, 250, and 270 K. At 230 K (Fig. 10(a)), it is clear that the GEMC results of Oakley et al. deviate from the experimental data in the high pressure liquid region and the low pressure vapor region. The estimated critical pressure from the work of Oakley et al. seems to be far below the experimental data. Although the critical points could not be calculated due to the limitation of the GEMC method, the calculated points by the present work agree well with the experimental data in the whole range and the estimated critical points also seem to be proximate to the experimental points. Fig. 10(b) shows that the accuracy of the present work does not decrease, even at high temperatures. The phase properties obtained from the GEMC simulations for CH4 /CO2 mixtures are presented in Table II. The values in parentheses are the estimated uncertainties, all of which are included in the symbols in Fig. 10. In the last column of Table II, the densities of the gas phases calculated from the VEOS4 are inserted for comparison. The densities obtained from the VEOS4 coincide well with simulations below 50 bar, but are underestimated at higher

TABLE II. The phase properties of CH4 /CO2 mixtures from the GEMC simulation of the present work. Values in parentheses are the uncertainties of the last digits. Temperature K 230

250

270

Pressure bar

xCH4 -. . .

yCH4 ...

20.265 32.423 40 45.6 51 55.728 60 30.397 40.529 50.662 60.794 50.631 58.575 70.207

0.043(1) 0.102(1) 0.154(2) 0.197(3) 0.251(5) 0.298(4) 0.352(4) 0.048(1) 0.088(2) 0.143(2) 0.221(2) 0.074(1) 0.107(3) 0.164(3)

0.518(2) 0.672(2) 0.719(2) 0.728(2) 0.740(2) 0.741(2) 0.755(2) 0.372(3) 0.479(2) 0.542(1) 0.590(2) 0.290(1) 0.346(3) 0.398(3)

cm3

ρl mol−1

0.02521(3) 0.02437(4) 0.02365(5) 0.02521(7) 0.02217(9) 0.02176(8) 0.02023(10) 0.02337(4) 0.02252(6) 0.02169(5) 0.02017(10) 0.02018(7) 0.01945(10) 0.01802(17)

cm3

ρg mol−1

0.001257(2) 0.002163(4) 0.002796(6) 0.001257(11) 0.004114(26) 0.004810(25) 0.005356(34) 0.001912(4) 0.002676(7) 0.003569(15) 0.004688(34) 0.003456(17) 0.004267(29) 0.005214(86)

ρ g VEOS4 cm3 mol−1 0.00125 0.00216 0.00282 0.00344 0.00407 0.00469 0.00557 0.00188 0.00265 0.00352 0.00466 0.00348 0.00421 0.00572

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01

064303-9

S. J. Pai and Y. C. Bae

pressures. It can be inferred from these results that more virial terms are needed to precisely predict the saturated vapor pressure above 50 bar. V. CONCLUSION

In order to improve our ability to predict the properties of natural or shale gas systems, we studied the macroscopic behaviors of CH4 , CO2 , and their mixtures using quantum mechanics and molecular simulations. Because a quantitative description of the phase properties requires accurate potential energies, we conducted CCSD(T) level of ab initio calculations for the CH4 –CO2 dimer. The global minimum was found to be near the point where the side of CO2 parallels to the perpendicular of one of the faces of the CH4 molecule. There was a 2%–3% lowering of the well depth by elevating the calculation level from MP2 to CCSD(T), which will likely produce only a small effect in practical problems. To increase the simulation speed, a relatively simple but accurate pair potential function was proposed and only the ATM expression was used as the nonadditive three-body potential function. For pure CH4 , the lap time for a molecular simulation was shortened by about 2.5 times compared to the original potential function of Hellmann et al. The reliability of the proposed potential functions was verified by comparing them with experimental virial coefficients. For CO2 , we used the original pair potential function of Bukowski et al. due to problems in fitting, but the advanced schemes in Eqs. (3)–(5) may facilitate the application of the present potential function. The overestimation near the critical point in the third virial coefficient showed the need for more terms in nonadditive three-body interactions (such as induction by quadrupoles). A VEOS4 with the derived virial coefficients agreed well with the pVT data of pure CH4 , CO2 , and their mixtures up to 200 bar at 300 K. The vapor-liquid coexistence curve calculated by the present potential functions and the GEMC simulations agreed well with the experimental data. The results were particularly accurate for binary systems, even near the critical point. Our results differed from the results of Oakley et al., from which it was proven that quantitative prediction for phase behaviors of simple and nonpolar molecules is possible if there are accurate two- and three-body potential functions. To go a step further, more quantum mechanical studies are needed in order to expand the proposed work to complex molecules such as ethane or larger hydrocarbons that have internal degrees of freedom. ACKNOWLEDGMENTS

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science, and Technology (Grant No. 2011-0024330). 1 K.

J. Chew, Philos. Trans. R. Soc. A 372, 20120324 (2014). E. Voulgaris, C. J. Peters, and J. de Swaan Arons, Int. J. Thermophys. 16, 629 (1995).

2 M.

J. Chem. Phys. 141, 064303 (2014) 3 P.

Hobza, R. Zahradník, and K. Müller-Dethlefs, Collect. Czech. Chem. Commun. 71, 443 (2006). 4 R. Hellmann, E. Bich, and E. Vogel, J. Chem. Phys. 128, 214303 (2008). 5 R. Bukowski, J. Sadlej, B. Jeziorski, P. Jankowski, K. Szalewicz, S. A. Kucharski, H. L. Williams, and B. M. Rice, J. Chem. Phys. 110, 3785 (1999). 6 S. Bock, E. Bich, and E. Vogel, Chem. Phys. 257, 147 (2000). 7 S. Bock, E. Bich, E. Vogel, A. S. Dickinson, and V. Vesovic, J. Chem. Phys. 117, 2151 (2002). 8 S. Bock, E. Bich, E. Vogel, A. S. Dickinson, and V. Vesovic, J. Chem. Phys. 120, 7987 (2004). 9 C. Bratschi, H. Huber, and D. J. Searles, J. Chem. Phys. 126, 164105 (2007). 10 R. Bukowski and K. Szalewicz, J. Chem. Phys. 114, 9518 (2001). 11 M. T. Oakley and R. J. Wheatley, J. Chem. Phys. 130, 034110 (2009). 12 M. T. Oakley, H. Do, and R. J. Wheatley, Fluid Phase Equilib. 290, 48 (2010). 13 S. F. Boys and F. Bernardi, Mol. Phys. 19, 553 (1970). 14 K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chem. Phys. Lett. 157, 479 (1989). 15 B. M. Axilrod and E. Teller, J. Chem. Phys. 11, 299 (1943). 16 Y. Muto, Proc. Phys. Math. Soc. Jpn. 17, 629 (1943). 17 J. K. Singh and D. A. Kofke, Phys. Rev. Lett. 92, 220601 (2004). 18 K. M. Benjamin, J. K. Singh, A. J. Schultz, and D. A. Kofke, J. Phys. Chem. B 111, 11463 (2007). 19 O. Akin-Ojo and K. Szalewicz, J. Chem. Phys. 123, 134311 (2005). 20 A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, J. Olsen, and A. K. Wilson, Chem. Phys. Lett. 286, 243 (1998). 21 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., GAUSSIAN 03, Revision D.01, Gaussian, Inc., Wallingford, CT, 2004. 22 See supplementary material at http://dx.doi.org/10.1063/1.4891983 for used geometries and ab initio calculation values of CH4 and CO2 dimer and fitting constants in virial coefficients. 23 W. E. Rice and J. O. Hirschfelder, J. Chem. Phys. 22, 187 (1954). 24 V. F. Lotrich and K. Szalewicz, J. Chem. Phys. 106, 9668 (1997). 25 Y. Midzuno and T. Kihara, J. Phys. Soc. Jpn. 11, 1045 (1956). 26 E. A. Mason and T. H. Spurling, The Virial Equation of State (Pergamon Press, Oxford, 1969). 27 K. O. Monago, Chem. Phys. 316, 9 (2005). 28 J. Wiebke, Chem. Phys. 411, 43 (2013). 29 S. J. Pai and Y. C. Bae, Fluid Phase Equilib. 338, 245 (2013). 30 J. A. Barker, P. J. Leonard, and A. Pompe, J. Chem. Phys. 44, 4206 (1966). 31 R. Hellmann, Mol. Phys. 111, 387 (2013). 32 G. J. Esper, D. M. Bailey, J. C. Holste, and K. R. Hall, Fluid Phase Equilib. 49, 35 (1989). 33 H. Hou, J. C. Holste, K. R. Hall, K. N. Marsh, and B. E. Gammon, J. Chem. Eng. Data 41, 344 (1996). 34 B. Jäger, R. Hellmann, E. Bich, and E. Vogel, J. Chem. Phys. 135, 084308 (2011). 35 D. R. Lide, CRC Handbook (CRC Press LLC, 2003). 36 A. Z. Panagiotopoulos, J. Phys.: Condens. Matter 12, R25 (2000). 37 H. Knapp, R. Döring, L. Oellrich, U. Plöcker, and J. M. Prausnitz, VaporLiquid Equilibria for Mixtures of Low Boiling Substances (DECHEMA, Frankfurt/Main, 1982). 38 M. L. Martin, R. D. Trengrove, K. R. Harris, and P. J. Dunlop, Int. Data Ser. Sel. Data Mixtures, Ser. A 1987, 57. 39 K. Ohgaki, T. Mizuhaya, and T. Katayama, J. Chem. Eng. Jpn. 14, 71 (1981). 40 P. J. McElroy, L. L. Kee, and C. A. Renner, J. Chem. Eng. Data 35, 314 (1990). 41 H. H. Reamer, R. H. Olds, B. H. Sage, and W. N. Lacey, Indust. Eng. Chem. 36, 88 (1944). 42 B. V. Mallu and D. S. Viswanath, J. Chem. Thermodyn. 22, 997 (1990). 43 R. Gilgen, R. Kleinrahm, and W. Wagner, J. Chem. Thermodyn. 26, 383 (1994). 44 E. Holleran, Fluid Phase Equilib. 251, 29 (2007). 45 See http://webbook.nist.gov/chemistry/fluid/ for NIST homepage. 46 B. D. Smith and R. Srivastava, Thermodynamic Data for Pure Compounds (Elsevier Science Publishers, 1986).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Tue, 16 Sep 2014 08:46:01