Theoretical Studies of Energetics and Diffusion of Aromatic

Chemical Engineering Program, Federal University of Rio de Janeiro, Rio de Janeiro, Brazil. Atomic and molecular interactions of aromatic compounds in...
0 downloads 0 Views 128KB Size
Download PDF
Ind. Eng. Chem. Res. 2000, 39, 227-235

227

Theoretical Studies of Energetics and Diffusion of Aromatic Compounds in Supercritical Carbon Dioxide L. A. Ferreira Coelho,†,‡ A. Marchut,†,§ J. V. de Oliveira,| and P. B. Balbuena*,† Department of Chemical Engineering, University of South Carolina, Columbia, South Carolina 29208, and Chemical Engineering Program, Federal University of Rio de Janeiro, Rio de Janeiro, Brazil

Atomic and molecular interactions of aromatic compounds in carbon dioxide are studied with ab initio and molecular dynamics techniques. Ab initio calculations are used to determine the nature of the CO2-CO2, CO2-benzene, and benzene-benzene interactions. We select an explicit all-atom force field without partial atomic charges to describe the intermolecular and intramolecular pair interactions of CO2 with benzene and toluene. Molecular dynamics simulations are used to calculate diffusion coefficients for benzene and toluene at infinite dilution in CO2 along isotherms at 313.15, 323.15, and 333.15 K, in the density range from 1.35Fc to 2.10Fc (Fc ) critical CO2 density). Diffusion coefficients are also calculated with a model based on perturbation theory of simple liquids. The calculated diffusion coefficients agree fairly well with the experimental results of Suarez et al. (Chem. Eng. Sci. 1993, 48, 2419). Introduction The development of new technologies using supercritical fluids has been based almost entirely on experimental work, and one of the drawbacks is the scarcity of thermodynamic data needed for design and scale-up of processes.1 Most of the theoretical and experimental studies on supercritical fluids have focused on the description of phase equilibria, whereas the microscopic mechanisms of solute dissolution and transport in supercritical solvents have not been extensively addressed and, therefore, are not well understood. The development of theoretical methods for the description of electronic, atomic, and molecular interactions has received a large impulse due to the tremendous increase in computational power achieved in the present decade. Sophisticated solutions of the Schro¨dinger equation permit a precise study of the pair interactions for almost any molecular system.2-5 On the other hand, the theoretical study of systems consisting of many hundreds or thousands of atoms still requires the use of some radical approximations. These necessary approximations for the study of large complex systems consist of the use of effective force fields for the representation of intermolecular interactions.6 Force fields are used in molecular dynamics (MD) and Monte Carlo (MC) techniques for the prediction of thermodynamic and transport properties of new and increasingly complex systems.6-9 These simulations are routinely performed for systems containing thousands of particles.10 MD and MC techniques are also valuable tools in the development of theories of dense fluids. The most fruitful results in transport coefficients have been * To whom correspondence should be addressed. E-mail: [email protected]. † University of South Carolina. ‡ Permanent address: Chemical Engineering Program, Federal University of Rio de Janeiro. § Undergraduate participant in the NSF-REU program at University of South Carolina, Summer 1998. Current address: Department of Chemical Engineering, University of Pennsylvania, Philadelphia, PA 19104. | Federal University of Rio de Janeiro.

obtained for diffusion, because viscosity and thermal conductivity require averages over many molecules, whereas diffusion is an individual property for which greater accuracy can be obtained in the statistical averages.6-8 Force fields provide a description of the potential energy surface characterizing the interaction of a pair of atoms or molecules. They can be classified according to (1) the way the intramolecular and intermolecular potential energies are partitioned, (2) the functional form assigned to each term, and (3) the type of information used for fitting the parameters. The main intermolecular interactions are usually grouped as long range (electrostatic) and short range (van der Waals), while the intramolecular interactions may include terms such as bond stretching, valence angle bending, dihedral angle torsion, and inversion.11 Functional forms are usually parametrized to specific molecules, and parameter databases are available for a variety of systems.12 The development of new functional forms requires fitting of selected properties in a given phase space region to either experimental or ab initio values. Several force fields have been published for CO2 and aromatic molecules. However, not all of them have been developed or tested in the supercritical CO2 region. Moreover, it has been recently demonstrated9 that various sets of force-field parameters may agree on giving a good description of structural properties, but the description may not be so good for dynamic properties. This problem may be intensified in the critical region because of long-range correlation effects. Here we report the main pair interactions determined by ab initio techniques for supercritical CO2 solutions of aromatic molecules and use this information for the selection of appropriate force fields. The force fields are used in MD simulations covering three temperatures (313.15, 323.15, and 333.15 K) and densities from 1.35Fc to 2.10Fc. The benzene-CO2 and toluene-CO2 simulations were done at nine supercritical state points. The calculated diffusion coefficients were compared with those from a model based on the perturbation theory of simple liquids13 and with experimental results.14

10.1021/ie990266i CCC: $19.00 © 2000 American Chemical Society Published on Web 11/24/1999

228

Ind. Eng. Chem. Res., Vol. 39, No. 1, 2000

Table 1. Geometric Parameters of CO2 Monomer and Dimera

a

method

R (Å)

R1 (Å)

R2 (Å)

R3 (Å)

φ1 (deg)

φ2 (deg)

HF/6-31G (this work) B3PW91/6-311G** (this work) HF/6-31G*18 MP2/6-31G*19 experimental43 experimental21 experimental20

1.161 1.159 1.143 1.179 1.160

3.740 4.089 3.777 3.473

1.159 1.160 1.142 1.177

1.163 1.158 1.145 1.180

179.36 179.70 179.60 179.60

126.96 125.64 123.60 120.40

3.602 3.599

122.00 121.80

R is the CO bond length; all of the other variables are dimer parameters defined in Figure 1. Table 2. CO2 Vibrational Frequencies (cm-1) vibrational mode

Figure 1. Schematic representation of the CO2 dimer. See Table 1 for calculated and experimental distances and angles.

2. Ab Initio Calculations 2.1. Methodology. The geometries of the molecular systems under study were fully optimized using Hartree-Fock (HF) and density functional theory (DFT). We obtained ground-state geometries and energies and vibrational frequencies. HF is an ab initio technique that has proved useful to reproduce molecular geometries and vibrational frequencies with good accuracy.15 However, its main deficiency is its intrinsic lack of electron correlation, resulting in an overestimation of the attractive forces between molecules. DFT is an alternative technique to solve the Schro¨dinger equation that takes into account the correlation between electrons, and it is computationally inexpensive compared to ab initio methods yielding equivalent accuracy.16 To achieve the highest possible accuracy using ab initio and DFT methods, one must carefully select the basis set used to construct the solution wave function. Calculation of systems dominated by dispersion interactions requires inclusion of electron correlation effects and the use of polarization functions. We used the 6-31G and 6-311G** basis sets for the HF and DFT calculations, respectively.15 All of these calculations were carried out with the program GAUSSIAN 94.17 We examined single molecules and systems of up to five molecules. For systems of more than two molecules, we tried several initial configurations. The values reported correspond to the most stable structures. However, we have not explored all possible spatial distributions. Therefore, some of the optimized structures may not correspond to absolute minima. 2.2. Results. Table 1 shows calculated and experimental bond lengths and angles corresponding to the CO2 molecule and its dimer. Several authors have discussed the dimer structure.18,19 We found a T-shaped conformation (C2v) and a slipped parallel structure (C2h), the latter of which is shown in Figure 1. However, experimental data20,21 and high-level ab initio calculations18,19 yielded C2h as the most stable, whereas the T-shaped conformation seems to be a transition state or a local minimum. Our results yield the slipped parallel C2h structure displayed in Figure 1, with parameters listed in Table 1, as the lowest energy structure, but the energy difference with the T-shaped structure is only 0.2 kcal/mol. Both HF and DFT methods yield geometric parameters that agree well with the experimental results. The

method

bend

symmetric stretch

asymmetric stretch

HF/6-31G B3PW91/6-311G** experimental44

656 670 667

1407 1386 1333

2374 2463 2349

Table 3. Ground-State and Binding Energies of a Pure CO2 System no. of CO2 molecules 1 2 3 4

energy (hartrees) HF DFT -187.514 95 -375.032 98 -562.550 07 -750.067 66

-188.563 82 -377.128 52 -565.692 51 -754.256 85

binding energy (kcal/mol) HF DFT 0.0 -1.93 -3.27 -4.93

0.0 -0.55 -0.66 -0.99

calculated harmonic frequencies for the monomer (Table 2) show good agreement with the experimental results. The differences between calculated and experimental frequencies may be attributed, in part, to the fact that anharmonicity effects are not included in the calculations, as has been demonstrated previously for a variety of systems.15,22 Despite the good agreement obtained for geometries and frequencies (Tables 1 and 2), HF is not an appropriate method to estimate binding energies dominated by dispersion forces. The reason is that HF does not include electron correlation effects that are considered the source of these forces. On the other hand, DFT is an accurate tool to estimate binding energies in interactive systems, and it provides the right trend regarding weak interactions, as has been shown for aromatic molecules.23 The binding energies Ebinding for complexes of two, three, and four CO2 molecules have been calculated as differences of their ground-state energies

Ebinding ) En - (En-1 + E1)

(1)

where n is the number of CO2 molecules in the complex. Note that each molecule (or complex) has an electronic distribution that results in its characteristic geometry, and its associated energy is in a global minimum (ground state). The binding energies in Table 3 are calculated using the ground-state values reported in the same table. Figure 2 shows the DFT-optimized structures for the CO2 trimer and the tetramer. Both complexes show C2v symmetry and an extended Tshaped configuration. The small values of their binding energies (Table 3) confirm the weak nature of their attractive interactions. We found T-shaped geometries as minimum energy structures for the trimer and tetramer; however, we do not discard that parallel slipped configurations can result in local minima. In close analogy to the CO2 dimer, the benzene dimer is dominated by dispersion interactions; however, the most stable structure for the benzene dimer is still

Ind. Eng. Chem. Res., Vol. 39, No. 1, 2000 229

Figure 2. Calculated (B3PW91/6-311G**) CO2 complexes. Shortest distances between intermolecular sites are as follows. (a) Trimer. C4-O3 ) C4-O8 ) 3.5 Å. (b) Tetramer: C4-O3 ) C4-O8 ) C10-O9 ) 3.5 Å. Table 4. Ground-State and Binding Energies of a Pure Benzene System energy (hartrees)

no. of benzene molecules

HF

DFT

1 2a

-230.624 47 -461.249 93

-232.213 61 -464.427 60

a

binding energy (kcal/mol) HF DFT 0.0 -0.62

0.0 -0.23

Dimer in the C2h configuration.

controversial. Several configurations, such as T-shaped, planar sandwich, parallel displaced, and V-shaped, have been calculated using ab initio calculations at the DFT,23 Moller-Plesset MP2 and MP4,24,25 and coupled cluster levels.26 Although there is not total agreement regarding the values of the binding energies, it can be concluded that the T-shaped and parallel-displaced structures are energetically similar. Our DFT calculations yielded the parallel-displaced structure as the most stable benzene dimer, which is also the arrangement of two consecutive layers of benzene rings in the graphite structure. Absolute and binding energies for these systems are reported in Table 4. As in eq 1, the binding energy of an AB complex is defined as the difference between the ground-state energy of AB and the sum of the energies of A and B in their ground states. Parts a-d of Figure 3 display the DFT-optimized structures for the complexes of a benzene molecule with one, two, three, and four CO2 molecules. Their corresponding total and binding energies listed in Table 5 show weak interaction energies typical of van der Waals systems. The shortest distance in this complex corresponds to the O (from CO2) to H (from C6H6) distance, 3.6 Å, indicating the strongest interaction in this complex. A second CO2 molecule (Figure 3b) interacts with the first CO2 and with the benzene ring, because the like interaction is only slightly more favorable than the unlike. This can be inferred by comparison of the binding energies in Tables 3 and 5. The T-shaped configuration of the CO2 dimer remains in the (CO2)2C6H6 complex, although somehow disrupted by the presence

Figure 3. Calculated (B3PW91/6-311G**) CO2-benzene complexes. Shortest distances between intermolecular sites are as follows: (a) O15-H8 ) 3.63 Å; (b) O15-H8 ) 3.61 Å, O17-H8 ) 3.07 Å, O15-C16 ) 3.44 Å; (c) O15-H2 ) 3.49 Å, O15-C16 ) 3.44 Å, O14-C20 ) 3.44 Å, angle C20C13C17 ) 163°; (d) O18-H2 ) 3.55 Å, O15-H11 ) 3.49 Å, O14-C19 ) 3.39 Å, O21-C16 ) 3.51 Å, O23-C16 ) 3.51 Å, O17-C13 ) 3.41 Å.

230

Ind. Eng. Chem. Res., Vol. 39, No. 1, 2000

Table 5. Ground-State and Binding Energies of System CO2 + Benzene, Corresponding to Structures in Figure 3 energy (hartrees)

binding energy (kcal/mol)

type of mixture

HF

DFT

HF

DFT

benzene + 1CO2 benzene + 2CO2 benzene + 3CO2 benzene + 4CO2

-418.141 50 -605.660 67 -793.178 38 -980.698 16

-420.778 07 -609.343 04 -797.907 22 -986.473 47

-1.30 -3.95 -5.68 -8.71

-0.40 -1.12 -1.35 -2.87

of the benzene molecule. A similar situation results by addition of third and fourth CO2 molecules (Figure 3c,d). The new molecules accommodate upon the previous structures, trying to maximize the most favorable interactions. In Figure 3d we observe that all CO2 molecules locate on one side of the benzene plane; two of them interact more closely with benzene. There is a certain degree of similarity between the nature of the electronic distribution in the benzene and CO2 systems. The trimer and tetramer of benzene have been recently calculated27 using an ab initio based force field that includes manybody interaction terms, yielding configurations similar to those shown in Figures 2 and 3. The resemblance between spatial arrangements of n-mers in benzene and in CO2 shows that intermolecular forces in these systems are of a similar nature. In benzene and CO2 mixtures, there is a competition between three different sets of weak interactions (11, 22, and 12, where 1 and 2 are the pure components). The resultant configurations that correspond (at least) to local minima reflect the effect of interatomic and intermolecular forces arising from interactions among electrons and nuclei from each of the atoms in the system. Evidently, in a macroscopic system, these ordered structures can only be detected at very low densities and may be present (although short-lived) at supercritical conditions. Note also that dispersion forces arise because of fluctuations in the electronic distribution that will take place at any temperature. In the next section, we use effective force fields, which are parametrized functional forms designed to represent the true interactions, and we compare potential energy curves calculated using these force fields with the DFT results. 3. Molecular Dynamics 3.1. Force Fields. Earlier force fields for carbon dioxide include one-site,13,28-30 two-site,31 and three-site models.31,32 Some of these empirical intermolecular potentials are presented in Table 6. The critical temperature and critical density given by the one-site models are also included in Table 6. These critical properties were calculated using the Lennard-Jones (LJ) 6-12 reduced critical properties:33 Tc* ) kbT/ ) 1.35, Fc* ) Fσ3 ) 0.349. The model parameters in Table 6 have been adjusted to meet several criteria. The TD model is a LJ 6-12 potential used by Tom and Debenedetti34 to study the behavior of microstructure and solvation in supercritical mixtures. These parameters were fitted to reproduce the CO2 critical density and critical temperature. The second and third models in Table 6, LSM30 and DCO,13 were parametrized by fitting to CO2 self-diffusion measurements. However, LSM’s fitting was done in a wider temperature range (from 273 to 373 K) than DCO’s, where only supercritical CO2 state points were included. As shown in Table 6, the critical point is largely overestimated by LSM and

underestimated by DCO. Iwai et al.29 calculated diffusion coefficients of naphthalene and 2-naphthol in supercritical CO2 using a one-site model (YHUA). Its LJ parameters were obtained using MC simulations aimed to reproduce experimental data of solubilities of aromatic compounds in supercritical CO2. The MSB model parameters have been fitted to the second virial coefficient data and static lattice properties of carbon dioxide;35 the model reproduces qualitatively properties for CO2 in the liquid state.31 Murthy et al.31 developed new two-site models MSM1, MSM2, and MSM3, where the LJ centers, separated by a distance l (Table 6), coincide approximately with the positions of the two oxygen atoms. The third group in Table 6 includes three-site models where the third site is added to represent the carbon atom. The MSM4 model was parametrized to reproduce the structure and energetics of crystalline CO2 as well as second virial coefficient data.31 The HY1 model includes LJ terms and charges on each atom that were optimized to reproduce the molecular quadrupole moment.32 To improve the prediction of the critical properties, new parameters were developed (HY2 model; Table 6).32 The last model in Table 6 is the three-site Universal Force Field (UFF)36 that was used in this work. To estimate its critical properties, we fitted the pair interaction energy for the dimer in the parallel configuration to an equivalent one-site LJ (12-6) potential function. We found that the best fitting for the equivalent one-site model is given by /kb ) 225.5 K and σ ) 3.57 Å. The UFF critical properties in Table 6 were estimated using these one-site parameters. To illustrate the prediction of the models given in Table 6 and to compare their performance against the DFT results, we computed the pair potential energy curve for CO2 molecules (i and j) as a function of their intermolecular distance rij. For the two-site and threesite models, the potential energy curve was calculated for the T-shaped and slipped parallel configurations, with rij being the C-O distance in the first case and the C-C distance in the second case. For the three-site UFF model results shown in Figures 4-6, we did not include the intramolecular contributions, because for CO2 these forces contribute only 1-2% to the total energy. Figure 4 shows that the LSM model overestimates the potential well energy. Both the DCO and LSM parameters were fitted to CO2 self-diffusion coefficients (although at different state conditions). However, the magnitude of the  parameter for LSM (Table 6) is more than double that of the DCO model, and the σ parameter is lower for LSM than for DCO. Both differences tend to increase the interaction strength. Molecules will remain attached to each other for longer times than they would if they were interacting with weaker forces; therefore, diffusion will be slower. Thus, we expect that the diffusion coefficients will be underestimated by LSM at supercritical conditions. Figure 4 also shows the results of a two-site model (MSM1), along with those from DFT and from the UFF used in this work, all calculated for the dimer parallel configuration. The MSM1 model gives a poor representation of the potential energy curve. The UFF model yields a minimum of about 0.5 kcal/mol at a C-C distance near 4.0 Å, slightly lower than the DFT value (4.09 Å), and its corresponding minimum energy is in good agreement with that from DFT. The one-site DCO model, which contains LJ

Ind. Eng. Chem. Res., Vol. 39, No. 1, 2000 231 Table 6. Lennard-Jones (6-12) and Other Geometric Parameters for CO2 Modelsa model TD34 LSM30 DCO13 YHUA29 MSB35 MSM131 MSM231 MSM331,c MSM431 HY132 HY232 UFF36

/kb (K)

σ (Å)

225.5 500.71 195.2 236.1 129.7 136.7 130.0 90.8 29.0 (C-C) 83.1 (O-O) 28.99 (C-C) 82.99 (O-O) 28.1 (C-C) 80.5 (O-O) 52.84 (C-C) 30.20 (O-O)

3.794 3.262 3.681 3.720 2.908 2.979 3.033 3.529 3.126 (C-C) 3.383 (O-O) 2.785 (C-C) 3.064 (O-O) 2.757 (C-C) 3.033 (O-O) 3.431 (C-C) 3.118 (O-O)

Θ (DÅ)b

-3.89 -3.66 -3.69 -4.59 -3.85

l (Å)

qcc (e)

Tc (K)

Fc (g/cm3)

304.4 675.9 263.52 318.73

0.467 0.735 0.511 0.495

304.4

0.560

2.38d 2.32d 2.10d 2.32d 2.32d 1.161f

0.6645

1.149f

0.6512

1.300f

a Model critical properties were estimated based on the critical properties of the LJ fluid.33 The experimental critical properties are T c ) 304.1 K and Fc ) 0.468 g/cm3.45 b Quadrupole moment experimental value: -4.3 DÅ.46 c The charge on each oxygen is -0.5qc. d O-O separation. e Lennard-Jones 9-6. f C-O bond length.

Figure 4. Potential energy curve for formation of a CO2 dimer, calculated with one-site models LSM30 and DCO,13 the two-site model MSM1,31 and the three-site UFF model36 including only the LJ part. In the two last cases, the distance rij corresponds to R1 in Figure 1.

parameters adjusted to experimental diffusion coefficients for supercritical CO2, yields a less accurate representation of the dimer potential energy surface than the UFF model but improved with respect to the LSM model. Figure 5 displays the potential energy curve for the T-shaped dimer, for DFT and UFF, in comparison with the three-site HY1 with and without Coulombic terms. The DFT-predicted T-shaped structure for the CO2 dimer yields a minimum of magnitude similar to that of the parallel structure, located at a C-O distance of 3.5 Å. The HY1 model with charges produces a too deep potential well, which largely overestimates the attractive interactions in these systems. Again, these strong pair interaction forces will affect dynamic properties. On the other hand, very good agreement is found between the DFT and UFF potential energy curves. Similar features are observed for the slipped parallel conformation shown in Figure 6 for the same set of models of Figure 5.

Figure 5. Potential energy curve for formation of a CO2 dimer, calculated with three-site models HY132 with and without Coulombic terms and with the UFF model36 including only the LJ part. The distance rij corresponds to the carbon atom of one molecule and the oxygen atom of another molecule arranged in a T-shaped configuration.

To further test the effect of Coulombic forces, we used three-site LJ models with and without atomic charges for the calculation of diffusion coefficients of benzene at infinite dilution in supercritical CO2 using MD simulations. Our results (not shown) confirmed that, although both three-site models can reproduce structural results, a potential model with partial atomic charges largely underestimates diffusion coefficients in the supercritical region. The reason for this behavior is given by the weakly polar nature of the CO2 molecule. When charges are imposed on this molecule, unrealistic electrostatic forces tend to attract the solute, yielding longer lifetimes for the solute-solvent interactions. If the solute is strongly attached to the solvent, its diffusion will be slower than that in the case when only weak forces exist. We conclude that the best model for the CO2 molecule consists of three LJ sites without charges, and for the same reasons we adopted an atomistic description (LJ sites) for the aromatic compounds allowing intramolecular interactions. We employed the UFF potential function,36 which has proved to reproduce with good accuracy the conformational structure of organic mol-

232

Ind. Eng. Chem. Res., Vol. 39, No. 1, 2000

included to improve statistics. We used a equilibration stage of 100 and 200 ps for production of averages; the time step was 0.5 fs. In the computation of the diffusion coefficient we employed the Einstein relation6,39

1 D ) lim 〈|ri(t) - ri(t0)|2〉 tf∞ 6

Figure 6. Potential energy curve for formation of a CO2 dimer, calculated with three-site models HY132 with and without Coulombic terms and with the UFF model including only the LJ part. The distance rij corresponds to R1 in Figure 1. Table 7. Lennard-Jones (6-12) Parameters for Benzene and Toluene model

/kb (K)

σ (Å)

UFF36 (benzene)

52.84 (C-C) 22.14 (H-H) 531.0 52.84 (C-C) 22.14 (H-H) 575.0

3.431 (C-C) 2.571 (H-H) 5.008 3.431 (C-C) 2.571 (H-H) 5.283

DCO13 (benzene) UFF36 (toluene) DCO13 (toluene)

ecules for which charge apparently does not play an important role.37 The functional form of the UFF potential is given by

1 Eij ) Kij(r - rij)2 + Kijkl(C0 + C1 cos wijkl + 2 xij 6 xij 12 + + C2 cos 2wijkl) + Dij -2 x x m

Kijk

∑

n)0

{ [] []} m

Cn cos nθ + Kijkl2

∑ Cn cos nφijkl2

(2)

where ri(t) is the position of the aromatic compound at time t, averaged over many initial times t0. The calculated diffusion coefficients shown in Table 8 agree fairly well with the experimental values.14 Table 8 shows the average deviation (AD) found for benzene and toluene. The MD prediction is good, considering that the forcefield parameters have not been adjusted to data for these particular systems. Experimental errors were not reported by Suarez et al.;14 the authors only referred to the reproducibility of the measured diffusivities, which was found to be on the order of 1%. However, systematic errors might have been present. In fact, the authors describe several possible sources of errors, especially at low densities, but they do not quantify them. Although we have not calculated the critical point of CO2 with the three-site UFF model, our estimate shown in Table 6 (based on an effective one-site model) indicates that the solvent is in supercritical conditions at all states in these simulations. Density effects on the diffusion coefficients follow well-defined trends, as seen from Table 8. At constant temperature, diffusivities decrease when the density increases. This trend is found for both benzene and toluene in supercritical CO2. For benzene there is an exception to this trend, the value at 333.15 K and 0.863 g/cm3, which can be attributed to statistical errors, because systems dominated by weak interactions require very long simulations in order to reduce statistical uncertainties. The DCO model13,40 was also used for the calculation of diffusion coefficients (Table 8). Its LJ 6-12 parameters (Tables 6 and 7) were fitted to self-diffusion coefficients of pure components. Simple mixing rules were employed to predict mutual diffusivities in the whole range of compositions and in a wide range of temperatures.13,40 The self-diffusion coefficient can be written as41

n)0

where Kij is a bond force constant, rij is a standard or natural bond length, Kijkl is the force constant of a dihedral angle, C0, C1, and C2 are coefficients associated with a rotational barrier, wijkl is a dihedral angle, Kijk is an angle bending force constant, θ is a bond angle, Kijkl2 is an improper dihedral force constant, Cn are coefficients, and φ is the improper dihedral angle. Nonbonded interactions are represented by a LJ potential function with parameters Dij (well depth) and xij (van der Waals bond length). The rules to calculate all of the coefficients in eq 2 were described by Rappe et al.36 Because of the similar nature of the interactions between pairs CO2-CO2, and CO2-aromatic molecule, we used a similar force field, without electrostatic terms for CO2-benzene and CO2-toluene. The parameters for benzene and toluene are shown in Table 7. 3.2. MD Procedure and Results. All of the simulations were made in the NVT ensemble, using the Hoover-Nose´ thermostat which guarantees a canonical distribution in phase space.38 In all of the simulations we employed 300 molecules of CO2 and 3 molecules of the aromatic compound. More than one solute was

(3)

D)

3 kT 1/2 SP P (F*) 8Fσ2 πm

( )

(4)

where F is the density, m is the mass of the substance, T is the temperature, σ is the hard-sphere diameter, and PSP(F*) is a correction used to more accurately account for the density dependence of D41

(

PSP(F*) ) 1 -

F* [1 + F*2(0.40 - 0.83F*2)] (5) 1.09

)

where F* (≡Fσ3) is the reduced density. An effective hard-sphere diameter d(F*,T*) that depends on both LJ parameters42 can be used to substitute the size parameter σ in eq 4. For a binary system, the mutual diffusion coefficient is given by13

D12 )

(

3 kT 2 2πm 8Fmd12 12

)

1/2

PSP(F/m)

(6)

where d12 is the cross hard-sphere diameter and m12 is the reduced mass of solute 1 and solvent 2. The parameter d12 is calculated as the arithmetic mean of

Ind. Eng. Chem. Res., Vol. 39, No. 1, 2000 233 Table 8. Calculated (This Work) and Experimental14 Diffusion Coefficient D12 × 109 (in m2/s) for Benzene and Toluene in Supercritical CO2 benzene

a

toluene

T (K)

F (g/cm3)

exp

MD

DCO model

exp

MD

DCO model

313.15 313.15 313.15 323.15 323.15 323.15 333.15 333.15 333.15 % ADa

0.7810 0.8807 0.9361 0.7008 0.8350 0.9003 0.6071 0.7812 0.8630

12.99 10.27 9.01 15.58 11.72 10.58 18.18 13.60 11.57

13.01 9.30 10.01 17.73 10.22 8.77 16.51 12.40 15.98 13.36

12.37 9.69 8.42 15.32 11.09 9.44 19.65 12.89 10.59 6.28

12.23 10.39 9.32 15.26 11.93 10.32 17.18 13.42 11.41

11.65 9.50 9.00 13.81 17.57 9.17 17.79 13.57 8.10 13.14

11.23 8.80 7.65 13.91 10.07 8.57 17.84 11.70 9.61 12.80

NP % AD ) 100/NP∑i)1 |(Diexp - Dicalc)/Diexp|.

the pure-component hard-sphere diameters, d1 and d2, at the temperature of the mixture and density of the pure substances. The critical temperature of the DCO model (263.51 K) is largely underestimated with respect to the experimental model (304.1 K), and its critical density (0.511 g/cm3) is slightly above the experimental critical density (0.468 g/cm3). Therefore, we do not expect a good performance of this set of LJ parameters for the prediction of densities or other thermodynamic properties in the near-critical region. However, eq 6 provides an excellent model to predict diffusivities of aromatic compounds in the supercritical CO2 region (Table 8). Its functional form gives a simple explanation for the experimentally found density and temperature dependence. The success of this model can be explained in terms of the weak attractive nature of the pair interactions among solvent and solute molecules. To test the effect of a purely repulsive model, we have calculated the diffusion coefficients using the plain hard-sphere diameter σ instead of d(F*,T*) for all of the conditions in Table 8. The results of the model using σ (not shown) improve the mean deviation from 6.28 to 4.83% for benzene and from 12.80 to 11.03% for toluene. Clearly, repulsive forces are dominant under these conditions. We emphasize that this last remark is valid for the set of conditions that we have tested, i.e., at pressures and temperatures sufficiently far from the critical point, whereas close to the critical point, the role of attractive forces may become important. 4. Conclusions DFT was used to determine the most stable conformations for systems composed of one to four CO2 molecules with and without the presence of a benzene molecule. A slipped parallel structure found for the CO2 dimer and T-shaped conformations obtained for its trimer and tetramer resemble similar interactions between benzene and other aromatic molecules, dominated by weak dispersion interactions. MD simulations were used to examine previous models for CO2 based on simple Lennard-Jones potential functions with and without Coulombic terms. It was found that models including charges largely underestimate diffusion coefficients in the supercritical region. Another test of the different force fields was done by computing the potential energy curve for dissociation of the CO2 dimer. On the basis of these tests, we selected an all-atom force field without partial atomic charges to represent the interactions involving CO2, benzene, and toluene. Diffusion coefficients of benzene and toluene in supercriti-

cal CO2 were calculated using MD simulations and a simple model based on perturbation theory. Both simulation and perturbation theory results agree fairly well with experimental results in the entire range of supercritical conditions examined in this study. We conclude that the weak nature of the interactions in these systems is well represented by a LJ model without partial atomic charges. Further, we demonstrate that the choice of a force field is a crucial step that requires a detailed understanding of the systems involved. In summary, we have used three levels of theoretical analysis to examine the behavior of aromatic molecules in supercritical CO2. DFT calculations revealed the nature of intramolecular and intermolecular interactions, and effective force fields representing such interactions were used in MD simulations. Finally, a simple model based on the addition of attractive forces as a first-order perturbation to a hard-sphere model was used to compute diffusion coefficients. The diffusion coefficients calculated by the simple model and from the MD simulations agree with the experimental results. The results of the more sophisticated theoretical methods permitted us to interpret the success of the simple model and to become aware of its limitations. Acknowledgment A.M.’s participation in this work was supported by the NSF/REU program Grant DMR-9732227 at the University of South Carolina. We thank Andre´s Ma´rquez, who performed the DFT calculations for benzenebenzene interactions. L.A.F.C. is indebted to CNPqBrazilian Government for a fellowship and to Dr. Pedro Derosa and Andre´s Ma´rquez for helpful discussions. Partial support for this work from NSF Grant CTS9810053 to P.B.B. is gratefully acknowledged. Nomenclature l: distance between two Lennard-Jones centers for twosite models (MSB, MSM1, MSM2, MSM3, MSM4) or the distance between carbon-oxygen in the models HY1, HY2, and in this work, Å qc: charge on the carbon atom used in models HY1 and HY2 kb: Boltzmann constant, 1.380 662 × 10-23 J/K Kij: force constant, (kcal/mol)/Å2 rij: standard or natural bond length, Å Kijkl: force constant of the dihedral angle, kcal/mol wijkl: angle between the il axis and the ijk plane, deg d: effective hard-sphere diameter, Å D: diffusion coefficient, m2/s m: mass of the molecule, kg

234

Ind. Eng. Chem. Res., Vol. 39, No. 1, 2000

P: density correction r(t): position vector in Cartesian coordinates at a given time, Å r: distance between two molecules, Å t: time, s Dij: well depth, kcal/mol E: energy, kcal/mol or Hartrees n: number of CO2 molecules in a complex xij: van der Waals bond length, Å Kijk: angle bending force constant, kcal/mol Kijkl2: improper dihedral force constant, kcal/mol Cn: coefficients chosen to satisfy appropriate boundary conditions such as that the function has a minimum at the natural bond angle θ0 in the case of angle bend or in the case of torsion they are determined by the rotational barrier, dimensionless R1: distance between carbon centers (Figure 1), Å R2: distance between oxygen and carbon (Figure 1), Å R3: distance between oxygen and carbon (Figure 1), Å T: temperature, K Greek Symbols : Lennard-Jones parameter, kcal/mol σ: Lennard-Jones parameter, Å Θ: quadrupole moment, DÅ θ: bond angle, deg F: density of CO2, g/cm3 φ: improper dihedral angle, deg φ1: defined in Figure 1, deg φ2: defined in Figure 1, deg Subscripts 1: solvent 2: solute c: critical state 0: initial condition m: mixture Superscripts p: perturbation ref: reference *: indicative of reduced units SP: density correction proposed by Speedy41

Literature Cited (1) Brennecke, J. F.; Chateauneuf, J. E. Homogeneous organic reactions as mechanistic probes in supercritical fluids. Chem. Rev. 1999, 99, 433. (2) Jensen, F. Introduction to Computational Chemistry; Wiley: Chichester, U.K., 1999. (3) Molecular Dynamics: from classical to quantum methods; Balbuena, P. B., Seminario, J. M., Eds.; Elsevier Science: Amsterdam, The Netherlands, 1999; Vol. 7. (4) Simons, J.; Nichols, J. Quantum mechanics in chemistry; Oxford University Press: New York, 1997. (5) Foresman, J. B.; Frisch, A. Exploring chemistry with electronic structure methods, 2nd ed.; Gaussian, Inc.: Pittsburgh, 1996. (6) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1990. (7) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: New York, 1996. (8) Haile, J. M. Molecular Dynamics Simulation; John Wiley & Sons: New York, 1992. (9) Balbuena, P. B.; Wang, L.; Li, T.; Derosa, P. A. Ab initio and molecular dynamics studies of cation-water interactions. In Molecular dynamics: from classical to quantum methods; Balbuena, P. B., Seminario, J. M., Eds.; Elsevier Science: Amsterdam, The Netherlands, 1999; Vol. 7; pp 431-469. (10) Yin, D.; Mackerell, A. D. Combined ab initio/empirical approach for optimization of Lennard-Jones parameters. J. Comput. Chem. 1998, 19, 334.

(11) Israelachvilii, J. Intermolecular and surface forces; Academic Press: London, 1992. (12) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187-217. (13) Dariva, C.; Coelho, L. A. F.; Oliveira, J. V. A kinetic approach for predicting diffusivities in dense fluid mixtures. Fluid Phase Equilib. 1999, in press. (14) Suarez, J. J.; Bueno, J. L.; Medina, I. Determination of Binary Diffusion Coefficients of Benzenes and Derivatives in Supercritical Carbon Dioxide. Chem. Eng. Sci. 1993, 48, 2419. (15) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; John Wiley & Sons: New York, 1986. (16) Recent Developments and Applications of Modern Density Functional Theory; Seminario, J. M., Ed.; Elsevier Science Publishers: Amsterdam, The Netherlands, 1996; Vol. 4. (17) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A. GAUSSIAN 94, Revision E.1; Gaussian Inc.: Pittsburgh, 1997. (18) Illies, A. J.; McKee, M. I.; Schlegel, H. B. Ab initio studies of the CO2 dimer and the CO2 ion complexes (CO2)2+ and (CO2)3+. J. Phys. Chem. 1987, 91, 3489. (19) Eggenberger, R.; Gerber, S.; Huber, H. The carbon dioxide dimer. Mol. Phys. 1991, 72, 433. (20) Jucks, K. W.; Huang, Z. S.; Miller, R. E.; Fraser, G. T.; Pine, A. S.; Lafferty, W. J. Structure and vibrational dynamics of the CO2 dimer from the sub-Doppler infrared spectrum of the 2.7 µm Fermi diad. J. Chem. Phys. 1988, 88, 2125. (21) Walsh, M. A.; England, T. H.; Dyke, T. R.; Howard, B. J. Pulsed molecular beam infrared absorption spectroscopy of CO2 dimer. Chem. Phys. Lett. 1987, 142, 265. (22) Wong, M. W. Vibrational frequency prediction using density functional theory. Chem. Phys. Lett. 1996, 256, 391. (23) Seminario, J. M.; Zacarias, A. G.; Tour, J. M. Theoretical interpretation of conductivity measurements of thiolane sandwiches. A molecular scale electronic controller. J. Am. Chem. Soc. 1998, 120, 3970-3974. (24) Smith, G. D.; Jaffe, R. L. Comparative study of force fields for benzene. J. Phys. Chem. 1996, 100, 9624. (25) Chipot, C.; Jaffe, R.; Maigret, B.; Pearlman, D. A.; Kollman, P. A. Benzene Dimer: A good model for π-π interactions in proteinssA comparison between the benzene and the toluene dimers in the gas phase and in aqueous solution. J. Am. Chem. Soc. 1996, 118, 11217-11224. (26) Hobza, P.; Selzle, H.; Schlag, E. W. Potential energy surface for the benzene dimer. Results of ab initio CCSD(T) calculations show two nearly isoenergetic structures: T-shaped and parallel displaced. J. Phys. Chem. 1996, 100, 18790-18794. (27) Hobza, P.; Engkvist, O.; Selzle, H. L.; Schlag, E. W. Benzene trimer and benzene tetramer and properties determined by the nonempirical model (NEMO) potential calibrated from the CCSD(T) benzene dimer energies. J. Chem. Phys. 1999, 110, 5758. (28) Chialvo, A. A.; Debenedetti, P. G. Molecular dynamics study of solute-solute microstructure in attractive and repulsive supercritical mixtures. Ind. Eng. Chem. Res. 1992, 31, 1391. (29) Iwai, Y.; Higashi, H.; Uchida, H.; Arai, Y. Molecular dynamics simulation of diffusion coefficeints of naphthalene and 2-naphthol in supercritical carbon dioxide. Fluid Phase Equilib. 1997, 127, 251. (30) Liu, H.; Silva, C. M.; Macedo, E. A. Unified approach to the self-diffusion coefficients of dense fluids over wide ranges of temperature and pressure-hard sphere, square well, Lennard Jones and real substances. Chem. Eng. Sci. 1998, 53, 2403. (31) Murthy, C. S.; Singer, K.; McDonald, I. R. Interaction site models for carbon dioxide. Mol. Phys. 1981, 44, 135-143. (32) Harris, J. G.; Yung, K. H. Carbon Dioxide’s liquid-vapor coexistence curve and critical properties as predicted by a simple molecular model. J. Phys. Chem. 1995, 99, 12021-12024. (33) Powles, J. G. The liquid-vapor coexistence line for Lennard-Jones type fluids. Physica 1984, A126, 289-299.

Ind. Eng. Chem. Res., Vol. 39, No. 1, 2000 235 (34) Tom, J. W.; Debenedetti, P. G. Integral equation study of micostructure and solvation in model attractive and repulsive supercritical mixtures. Ind. Eng. Chem. Res. 1993, 32, 2118. (35) McRury, T. B.; Steele, W. A.; Berne, B. J. Intermolecular potential models for anisotropic molecules with applications to N2, CO2 and benzene. J. Phys. Chem. 1976, 64, 1288. (36) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. J. Am. Chem. Soc. 1992, 114, 10024-10035. (37) Casewit, C. J.; Colwell, K. S.; Rappe, A. K. Application of a universal force field to organic molecules. J. Am. Chem. Soc. 1992, 114, 10035. (38) Hoover, W. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695. (39) Einstein, A. Investigations on the theory of the brownian movement; Dover: New York, 1926. (40) Dariva, C.; Coelho, L. A. F.; Oliveira, J. V. Predicting diffusivities in dense fluid mixtures. Braz. J. Chem. Eng. 1999. (41) Speedy, R. J. Diffusion in the hard-sphere fluid. Mol. Phys. 1987, 62, 509.

(42) Souza, L. E. S.; Ben-Amotz, D. Optimized perturbed hard sphere expressions for the structure and thermodynamics of Lennard-Jones fluids. Mol. Phys. 1993, 78, 137. (43) Harmony, M. D.; Laurie, V. W.; Kuczkowski, R. L.; Schwendeman, R. H.; Ramasay, D. A.; Lovas, F. J.; Lafferty, W. J.; Maki, A. G. Molecular structures of gas-phase polyatomic molecules determined by spectroscopic methods. J. Phys. Chem. Ref. Data 1979, 8, 619. (44) Handbook of Chemistry and Physics, 77th ed.; Lide, D. R., Ed.; CRC Press: Boca Raton, FL, 1997. (45) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The properties of gases and liquids, 4th ed.; McGraw-Hill: Boston, 1987. (46) Buckingham, A. D.; Disch, R. L. Proc. Royal Soc. London, A 1963, 272, 275.

Received for review April 13, 1999 Revised manuscript received September 27, 1999 Accepted October 6, 1999 IE990266I