J. Phys. Chem. B 2004, 108, 7415-7423
7415
Cyclohexane-Benzene Mixtures: Thermodynamics and Structure from Atomistic Simulations Giuseppe Milano*,†,‡ and Florian Mu1 ller-Plathe‡ Dipartimento di Chimica, UniVersita di Salerno, I-84081 Baronissi (SA), Italy, and International UniVersity of Bremen, School of Engineering and Science, D-28759 Bremen, Germany ReceiVed: February 7, 2004; In Final Form: March 17, 2004
Mixtures of cyclohexane and benzene are investigated by means of all-atom molecular dynamics simulations. Thermodynamic properties such as density and enthalpy of mixing and diffusion coefficients were calculated for five different compositions and compared with experimental results. In agreement with experimental results, the calculated enthalpy of mixing indicates a slight nonideal behavior. Angle-dependent radial distribution functions have been analyzed in order to obtain detailed information about ring mutual orientations in pure liquids and in the mixtures. The most relevant configuration for benzene-benzene dimers in pure liquid and in a mixture with cyclohexane is the T-shape. Parallel configurations are distributed mainly in a broad range of shifting extent. In liquid cyclohexane the T-shape and sandwich configurations are equally relevant. Among the possible sandwich configurations, the nonshifted sandwich is the most relevant. In the cyclohexanebenzene dimers the nonshifted sandwich configuration dominates.
1. Introduction Atomistic simulations are a powerful tool in the analysis and interpretation of the thermodynamics and microscopic structure of liquids.1-4 Methods based on semiempirical force fields adapted to every single substance on the basis of experimentally observed quantities are accurate and reliable. The approach has been widely applied to various single component liquids, and some all-atom simulations of organic-organic or organicinorganic mixtures have been reported.5-10 Since it is wellknown that binary mixtures show a rich variety of interesting thermodynamic effects,11 it is promising to investigate them theoretically. Aromatic and aliphatic rings are two of the most elementary structures in organic chemistry. The cyclohexane ring represents the backbone of many sugars. Industrially, benzene and cyclohexane and their derivatives are used as solvents in the production and manufacture of several products. Cyclohexanebenzene mixtures are good model systems to test force-field parameters for polymer/solvent or polymer/polymer interactions. Furthermore, such force fields can be applied later in the study of polymer blends. In proteins the side chains of several amino acids consist of aromatic rings. Noncovalent interactions, especially those involving aromatic rings, are pivotal to protein-ligand recognition and therefore to drug design. Indeed, the vast majority of X-ray crystal structures of protein complexes with small molecules reveals bonding interaction involving aromatic amino acid side chains of the receptor and/or aromatic rings and heteroaromatic rings of the ligand.12,13 Several statistical analyses of PDB structures have elucidated the preferential arrangement of specific aromatic amino acid pairs in proteins or R-helices.14,15 The resulting general view is a competition between stacked and T-shaped complexes. * Corresponding author. E-mail:
[email protected]. † Universita di Salerno. ‡ International University of Bremen.
The benzene dimer has often been chosen as a model system for the investigation of these π-π interactions. High level ab initio calculations of benzene-benzene gas-phase dimers identified two main configurations, a T-shape and a shifted parallel one, as nearly isoenergetic.16,17 In the liquid state, the more favorable gas-phase interactions can be considerably damped by solvent and entropic effects. As revealed by 1H NMR investigations,18 the T-shape motif is preferred for liquid and crystalline benzene.19 Free energy profiles for benzene dimer association have been calculated using all-atom Monte Carlo simulations in liquid benzene, water, and chloroform.20 Additionally, some molecular dynamics simulations were used to better understand the local structure and mutual orientation of benzene rings in water.21 Furthermore, molecular dynamics simulations together with neutron diffraction studies on liquid benzene and mixtures with its fluorinated derivatives have been reported.22,23 With this precedent, cyclohexane-benzene mixtures can be a good model system to validate all-atom potential models for aromatic-aliphatic interactions. Furthermore, comparison between benzene and cyclohexane behavior and analysis of like (benzene-benzene, cyclohexane-cyclohexane) and unlike (benzene-cyclohexane) dimer configurations in neat liquid and in the mixtures are very useful for a better understanding of the role of the underlying interactions. We simulated five different systems. Thermodynamic properties such as the density and enthalpy of mixing are used to check the validity of the models employed. Previous all-atom models are available for the pure fluids, but not all atom models were developed for the properties of mixtures. United-atom simulations of cyclohexane-benzene mixtures have already been reported in the literature.24,25 However, several limitations of united-atom models have been previously demonstrated.26 This paper is organized as follows: In section 2.1 the technical details of the simulations are presented. Section 2.2 refines the definition of ring orientations. Section 3.1 delineates
10.1021/jp0494382 CCC: $27.50 © 2004 American Chemical Society Published on Web 05/05/2004
7416 J. Phys. Chem. B, Vol. 108, No. 22, 2004
Milano and Mu¨ller-Plathe
TABLE 1: Simulated Systems 1 2 3 4 5
x(C6H6)
C6H12
C6H6
N
0.00 0.25 0.50 0.75 1.00
250 324 215 72 0
0 108 215 215 215
250 432 430 287 215
t
sim
[ns]
6.62 5.65 6.82 9.0 10.0
TABLE 2: Geometry and Force Field Parameters of Cyclohexane and Benzenea parameter H C C-C C-H
the results for isolated dimer structure and stability. Calculated thermodynamic properties and liquid structure descriptions including radial distribution functions and angular correlations are reported in sections 3.2 and 3.3, respectively.
C-C-C C-C-H H-C-H
2. Simulation and Data Analysis 2.1. Details of Simulation. Different mixtures of cyclohexane (C6H12) and benzene (C6H6) were simulated at ambient conditions (T ) 300 K, p )1013 hPa). We investigated five systems (see Table 1) with different composition going from pure cyclohexane to pure benzene. The molecular dynamics package YASP27 was used to run the all-atom simulations under constant pressure (Berendsen manostat with coupling time τp ) 2 ps) and constant temperature (Berendsen thermostat28 with τT ) 0.2 ps at a time step of 2 fs). All bond lengths were kept rigid by the SHAKE procedure. The cutoff for nonbonded interactions was rc ) 1.1 nm with a Verlet neighbor list29 cutoff of 1.2 nm. No charges were taken into account for cyclohexane. The benzene model has small partial charges on the carbon and hydrogen atoms in order to reproduce the electric quadrupole moment of the benzene molecule. The geometric structure parameters are well-known, and from several experimental techniques bond lengths, angles, and dihedrals are available. The force field parameters were taken from former simulations;30,31 see Table 2. The nonbonded Lennard-Jones parameters have to be estimated from similar substances and then reoptimized individually. This is usually a lengthy procedure, which has to be done manually. For benzene we adjusted the nonbonded parameters of Jorgensen and Severance20 to reproduce the experimental density (at 1 atm) and heat of vaporization of liquid benzene. Cyclohexane nonbonded parameters have been obtained automatically by adapting the Lennard-Jones parameters to the experimental values of density and enthalpy of vaporization by the simplex method.30 This approach gives more accurate results than the traditional trial-and-error approach. With the benzene force field, different properties such as diffusion coefficients and the reorientational motion of benzene in the neat liquid and in polystyrene-benzene solutions have been measured and calculated.31,32 The cyclohexane force field has been used to calculate different properties such as thermodynamics properties, diffusion coefficients, and reorientation times for pure cyclohexane and cyclohexane-cyclohexene mixtures.8,30 The system was prepared as follows: The appropriate number of molecules was put at random positions and orientations in a large periodic box (about half the real density). This system was simulated for short intervals at a time step that started from 0.001fs and was increased to 2 fs. With the final value, the simulation was run until total energy and density had equilibrated, which took several hundred picoseconds. The density has been calculated (eq 1) from the box size and the heat of vaporization from the intermolecular nonbonded energy ∆Hvap ) - 〈Einter pot 〉 + RT. From the ∆Hvap values
C6H12
C6H6
Atomic Masses [u] M 1.00787 M 12.01
1.00787 12.01
Bond Length Constraints [nm] 0.1526 0.1090 Angle Potential φ0 [deg], kφ [kJ/(mol rad2)] φ0 109.5 Kφ 335 φ0 109.5 Kφ 420 φ0 109.5 Kφ 290
0.1390 0.1080 120 376.6 120 418.8
Torsions (Periodicity p )3), τ0 [deg], kτ [kJ/(mol rad2)] C-C-C-C τ0 180 kτ 10 Harmonic Dihedrals δ0 [deg], kδ [kJ(mol rad2)] 0 C-C-C-C δ0 kδ 167.4 C2-C3-C1-H [on C2] 0 δ0 kδ 167.4 C H
Lennard-Jones [kJ/mol], σ [nm] 0.299 σ 0.328 0.189 σ 0.258 Charges [q/e] 0.000b 0.000b
C H
0.294 0.335 0.126 0.242 -0.115 +0.115
a All intramolecular nonbonded interaction have been excluded. For definition of the mathematical form of the potential terms see ref 27. b In force field II the charges for C6H12 are -0.12 on carbons and +0.06 on hydrogens.
obtained from simulations for pure liquids A and B and their mixtures it is possible to calculate the mixing enthalpy:
- ∆HMIX ) ∆Hvap(mixture) - xA∆Hvap(A) (1 - xA)∆Hvap(B) (1) The ∆G (B f A) for the replacement of molecule A with a molecule B in the coordination environment of A is calculated from g(R) in the following way:
ln
gBA(R)
)-
gAA(R)
∆G(B f A) RT
(2)
Tracer diffusion coefficients are calculated from the meansquared displacements ∆r2 of the centers of mass of the molecules via the Einstein relation:
D≈
1d 〈∆r2〉 6 dt
(3)
Theoretically, eq 3 becomes exact in the limit t f ∞. Practically, the slope of the mean squared displacement is calculated in the linear regime. In our simulations of a NpT ensemble, the calculation of ∆r2 is nontrivial because the manostat rescales all coordinates in every step with a factor s depending on the difference between actual pressure and target pressure. This means an unphysical displacement of the atoms, which depends on their absolute positions, thus destroying translational invariance. A method, derived by one us,8 which minimizes the artifacts introduced
Cyclohexane-Benzene Mixtures
J. Phys. Chem. B, Vol. 108, No. 22, 2004 7417
TABLE 3: Calculated Tracer Coefficients in 10-5 cm2/s (T ) 300 K)a
a
x(C6H6)
C6H12
0.00 0.25 0.50 0.75 1.00
0.48 ( 0.02 (1.48) 0.84 ( 0.05 (1.78) 1.03 ( 0.08 (2.02) 1.29 ( 0.04 (2.13)
C 6H 6 1.04 ( 0.07 (2.17) 1.22 ( 0.06 (2.32) 1.43 ( 0.08 (2.36) 1.35 ( 0.04 (2.25)
In parentheses the experimental values from ref 33 are reported.
SCHEME 1: Ring Numbering Figure 1. Definition of T-shape and parallel shifted and nonshifted dimer configurations. θ is the angle between two ring normals. For a T-shape dimer the value of θ is 90°, and θ is 0° for a sandwich one. The angle between a vector joining the center of mass of two rings and one of the ring normals is φ. Nonshifted configurations present φ ) 0°.
by the rescaling of the periodic box, has also been applied in the present simulations. The diffusion coefficients of Table 3 were calculated for all the mixtures. Calculated and experimental33 diffusion coefficients of both components show an increase from cyclohexane to pure benzene, with a maximum located at a benzene molar fraction of 0.75. The calculated behavior of diffusion coefficients for both cyclohexane and benzene is in good qualitative agreement with the experimental one,33 with an underestimation of the diffusion coefficient by a factor of about 2. The geometries reported in section 3.1 have been obtained keeping the cyclohexane and benzene bond distances, angles, and dihedrals fixed to the values of the corresponding constraints or the minimum value of the potential reported in Table 2. Only the distance between the centers of mass, R, was varied. To locate the minimum energy geometries the potential energy surface was scanned by varying the distance, R, in steps of 0.01 nm. 2.2. Ring Orientation. The mutual orientation of two rings can be described by the angle θ between the plane normal nˆ 1 and nˆ 2 of two rings. For a six-membered ring, according to the labeling of Scheme 1, the plane normal to the ring is calculated as the mean of the normal of the two planes defined by carbon atoms 1-3-5 and 2-4-6. To define the T-shape and parallel shifted and nonshifted (or sandwich) arrangements of two rings the criteria explained below have been adopted. As depicted in Figure 1, the value of θ is 90° for the T-shape dimer and 0° for the sandwich one. Shifted and nonshifted arrangements can be discriminated considering the angle φ between the plane normal and the vector joining the centers of mass of two rings. 3. Result and Discussion 3.1. Isolated Dimer Structure and Interaction Energies. The structure of an isolated benzene dimer has been addressed in many experimental and theoretical studies.17 To our knowledge no studies have been reported about cyclohexane as well as cyclohexane-benzene dimers. In Figures 2, 3, and 4 the optimized geometries and the corresponding interaction energies for benzene, cyclohexane, and cyclohexane-benzene dimers are reported. In agreement with previous results,17 for the benzene dimer the T-shape and parallel displaced arrangements are favored with respect to the sandwich. On the contrary, for cyclohexane dimer,
Figure 2. Optimized geometries and interaction energies for like benzene-benzene and cyclohexane-cyclohexane dimers. Distances between rings center of mass are in nm. Energies are reported in kJ/ mol.
the sandwich structure is more stable by 5.5 kJ/mol than the T-shape arrangement. As for the unlike cyclohexane-benzene dimers, the sandwich arrangements indicated as S1 and S2 in Figure 3 are the most stable ones. To our knowledge no experimental data are available for cyclohexane-benzene dimers. However, these results can be compared with those obtained by rotational coherence spectroscopy for similar sandwich dimers formed by aromatic and aliphatic ring complexes such as fluorenecyclohexane and perylene-cyclohexane.34 The results reported in Table 5, obtained using a different force field (force field II) including atomic charges also for cyclohexane carbon and hydrogen atoms,35 are very close to those described above. This shows that, as a result of the low value for the quadrupole moment of cyclohexane (1040θ/Cm2 ) 3.0)36 with respect to that for benzene (1040θ/Cm2 ) -28.3),37 the electrostatic interactions are too small to determine molecular orientations. They can be easily included in an averaged way in the Lennard-Jones potentials. The results show that the stabilization for benzene-benzene and cyclohexane-cyclohexane dimers comes from different types of interactions. Although for the sandwich arrangement the Lennard-Jones stabilization
7418 J. Phys. Chem. B, Vol. 108, No. 22, 2004
Milano and Mu¨ller-Plathe TABLE 5: Interaction Energies for Isolated Dimers (Force Field II) force field II R [nm] electrostatic Lennard-Jones total [kJ/mol] T-shape P-displaced sandwich
0.519 0.450 0.380
Benzene-Benzene -2.14 -6.68 +1.42 -10.57 +5.73 -13.11
-8.82 -9.15 -7.38
T-shape sandwich
Cyclohexane-Cyclohexane 0.570 +0.02 -7.07 0.460 -0.06 -12.53
-7.05 -12.59
S1 S2 T1 T2 T3
0.430 0.420 0.520 0.510 0.559
Benzene-Cyclohexane +0.50 -11.56 +0.53 -11.78 -0.41 -7.59 -0.10 -8.80 -0.20 -6.04
-11.06 -11.25 -8.00 -8.90 -6.24
TABLE 6: Thermodynamic Properties of Mixtures (T ) 300 K) Figure 3. Optimized geometries and interaction energies for unlike benzene-cyclohexane dimers. Distances between rings center of mass are in nm. Energies are reported in kJ/mol.
F [g/cm3] x(C6H6)
sim
exptla
0.00 0.25 0.50 0.75 1.00
0.7742 0.7938 0.8160 0.8448 0.8752
0.7694 0.7864 0.8086 0.8357 0.8680
a
Figure 4. Calculated (b) and experimental (0) enthalpies of mixing (kJ/mol). The value of enthalpy of mixing obtained using a different force field (2) for the 1:1 mixture is also reported.
TABLE 4: Interaction Energies for Isolated Dimers (Force Field I) force field I
R [nm]
elecrostatic
Lennard-Jones
total [kJ/mol]
T-shape P-displaced sandwich
0.519 0.450 0.380
Benzene-Benzene -2.14 -6.68 +1.42 -10.57 +5.73 -13.11
-8.82 -9.15 -7.38
T-shape sandwich
0.590 0.460
Cyclohexane-Cyclohexane -7.16 -12.70
-7.16 -12.70
S1 S2 T1 T2 T3
0.430 0.430 0.530 0.51 0.579
Benzene-Cyclohexane -11.80 -11.91 -7.80 -9.17 -6.06
-11.80 -11.91 -7.80 -9.17 -6.06
term is about 6 kJ/mol larger than that for the T-shape one, the electrostatic interactions disfavor it. The benzene-benzene T-shape arrangement, being stabilized from both terms, results in the best compromise between Lennard-Jones and electrostatic interactions. For the cyclohexane dimer, the electrostatic interactions have no relevant role in determining the stabilization energy, and in this case the sandwich arrangement, maximizing the contact surface between the two molecules, is the most
∆Hmix [kJ/mol] sim
exptla
0.00 0.578 0.701 0.466 0.00
0.00 0.541 0.718 0.490 0.00
(T ) 303 K); from ref 38.
favored. In the same way the sandwich cyclohexane-benzene unlike dimer shows a stabilization energy larger than that of the T-shape one. It is worth noting that the stabilization energies calculated in this section include, to some extent, many-body effects inherent with the parametrization against bulk thermodynamic quantities. This contribution, however, can be considered small and practically constant for different molecular orientations. 3.2. Thermodynamic Properties. In Figure 4 the calculated and experimental enthalpies of mixing for different compositions are reported. In agreement with experiment38,39 the mixture behavior is not ideal, showing a positive excess enthalpy of mixing. The quantitative agreement between calculated and experimental mixing enthalpies is surprisingly good for all compositions (see Table 6). In fact, the force fields have not been designed to describe mixtures, and the common LorentzBerthelot mixing rules were used to treat unlike Lennard-Jones interactions. For the equimolar mixture, the enthalpy of mixing, as shown in Figure 4, has also been calculated utilizing a force field including charges on cyclohexane atoms.20 The calculated enthalpy of mixing is very close to the one obtained, not including point charges on cyclohexane atoms. These results are consistent with those obtained in the previous section for ∆E of stabilization of dimers. Also the reproduction of mixing enthalpy atomic charges on cyclohexane molecule are not necessary, and the small effect of electrostatic interactions can be included in averaged way in the Lennard-Jones potentials. The density increases monotonically from 0.7742 to 0.8752 g/cm3 (Table 6). For the pure systems, the values match the experimental ones, because the force field was optimized against them. A small offset between the experimental and simulated behavior was observed because the experiments were performed at a higher temperature (T ) 303 K).39 As shown in Figure 5, the increase for both experimental and calculated data is not linear. Fitting the density behavior with a quadratic polynomial, as indicated in Figure 5, gives
Cyclohexane-Benzene Mixtures
J. Phys. Chem. B, Vol. 108, No. 22, 2004 7419 TABLE 7: Coordination Numbers in First Neighbor Shella,b x(C6H6)
C6H12-C6H12
C6H12-C6H6
C6H6-C6H12
C6H6-C6H6
0.00 0.25 0.50 0.75 1.00
12.78 9.83 6.73 3.40
9.28 6.45 3.26
3.09 6.45 9.74
3.06 6.16 9.47 12.95
a Coordination numbers have been obtained by integration of the radial distribution of center of mass function to its first minimum. b The table head means that the first component has n next neighbors of the second component.
Figure 5. Calculated (b) and experimental (0) density (g/cm3). The value of A1 and A2 coefficients obtained from fitting experimental and calculated data by a second power polynomial are also reported.
Figure 7. Free energy (kJ/mol) for the exchange of a benzene with a cyclohexane molecule (solid line) in the benzene coordination environment and for the exchange of a cyclohexane with a benzene ring in the cyclohexane coordination environment (dashed line) as function of intermolecular distance.
Figure 6. Radial distribution functions of the centers of mass for (a) pure liquids and (b) equimolar mixture.
similar values for both the linear and quadratic coefficient. This good reproduction of thermodynamics properties can be ascribed to the individually reoptimized force fields, to yield heat of vaporization and densities of the pure liquids.30 3.3. Structure of Mixture. Radial Distribution Functions. In Figure 6a and b the radial distribution functions g(R) of the centers of mass are shown for the pure liquids and the equimolar mixture. The calculated coordination numbers, obtained by integrating the g(R) until its first minimum, are reported in Table
7. The system results fully mixed. So the coordination of cyclohexane and benzene is almost equal in the liquids. At low cyclohexane concentration both cyclohexane and benzene have slightly more cyclohexane neighbors and the same holds at high cyclohexane concentrations. From Figure 6b it is clear that for distances between 0.4 and 0.5 nm the cyclohexane-cyclohexane g(R) value is almost zero, whereas for the benzene-benzene and the benzene-cyclohexane pairs it is close to 1. This means that at small intermolecular distances the cyclohexane coordination shell is made up of benzene molecules. In Figure 7 are reported the calculated free energy values for the exchange of a benzene with a cyclohexane molecule in the benzene coordination environment (continuous line) and of a cyclohexane with a benzene molecule in the cyclohexane coordination environment (dashed line) as function of intermolecular distance. At small intermolecular distances replacement of a cyclohexane by a benzene molecule in the cyclohexane coordination environment leads to a lower free energy. In contrast, positive values of free energy correspond to replacement of cyclohexane with a benzene molecule in the benzene coordination environment. This feature can be easily explained in terms of the different size and shape of the two species. At close intermolecular distances the flatter benzene ring is able to intercalate between cyclohexane molecules. At larger intermolecular distances, the curves present some structure favoring slightly benzene or cylclohexane replacement until 1.3 nm. After that the mixing becomes practically random. In Figures 8 and 9, partial radial distribution functions are plotted for carbon and hydrogen atoms, respectively. For both cases, the intramolecular part was excluded for clarity. The carbon atoms show the anisotropy and local structures of the molecule. There is a shoulder at smaller distance, so some
7420 J. Phys. Chem. B, Vol. 108, No. 22, 2004
Figure 8. Intermolecular C-C radial distribution function of pure liquids and 1:1 mixture.
Figure 9. Intermolecular H-H radial distribution function of pure liquids and 1:1 mixture.
carbons can approach each other more closely. This feature results in a more pronounced change on going from pure cyclohexane to pure benzene. The hydrogen radial distribution functions are practically indistinguishable from the pure liquid and the mixture, and the interspecies RDF can be viewed as the average between the pure ones. Angular Correlations. For a detailed description of the structural features it is worthwhile to consider and compare the mutual orientations between the cyclohexane and benzene rings in the pure liquids and the mixtures. To this aim, as explained in the section 2.2, we refer to the absolute value of the scalar product |nˆ 1‚nˆ 2| between normal planes of two different rings. This is 1 for two coplanar rings, 0 for a T-shaped arrangement, and 1/2 for a random distribution of orientations. In Figure 10, this orientational distribution function (ODF), as a function of the average distance between a pair of benzene (continuous line) and cyclohexane (dashed line) rings in the pure liquids, is shown. As a result of the different thicknesses of benzene and cyclohexane the two corresponding curves appear shifted by about 0.1 nm. For liquid benzene, as already reported,31 at a close distance the rings are predominantly coplanar. However, as is evident from the RDF, there are only very few pairs at this distance. Around the typical nearest-neighbor distance of 0.5 nm, there is a minimum (
