Rare Gas−Benzene−Rare Gas Interactions: Structural Properties and

2266

J. Phys. Chem. A 2010, 114, 2266–2274

Rare Gas-Benzene-Rare Gas Interactions: Structural Properties and Dynamic Behavior Margarita Albertı´* IQTCUB, Departament de Quı´mica Fı´sica, UniVersitat de Barcelona, Barcelona, Spain ReceiVed: December 1, 2009; ReVised Manuscript ReceiVed: January 5, 2010

In the present work, some static and dynamic properties of trimers containing one benzene molecule and two rare gas atoms are investigated. These trimers can be formed in two different configurations, one in which the two rare gas atoms are placed in opposite sides of the benzene plane, (1|1), and the other in which the two atoms are placed on the same side, (2|0). The (1|1) configuration is more stable than the (2|0), and both minima are connected by small energy barriers. Accordingly, molecular dynamics simulations show frequent (2|0) T (1|1) interconversions, even at low temperatures. The time spent in each configuration has been related to the abundance of isomers. It has been found that at temperatures just below the dissociation, when interconversions are quite frequent, the relative abundance of (2|0) is always higher than that of (1|1), independently of the nature of the two rare gases. 1. Introduction The study of van der Waals complexes is of fundamental interest in chemistry, physics, and biology, playing a key role in the understanding of several properties of such complexes as, for instance, the mechanisms of growth and solvation and the nature of weak intermolecular interactions. In general, particular attention has been paid to molecular aggregates containing aromatic rings,1,2 mainly because of their role as models for solvation processes.3 For instance, the microsolvation of the benzene molecule (Bz) by several Ar atoms was investigated, the experimental results being consistent with the existence of different isomers.4 The existence of two general types for Bz-Ar195 aggregates allowed to understand the origin of the peak structure observed in the spectroscopic experiments of Hahn and Whetten.6 The high resolution of Fourier transform microwave (FTMW) spectroscopy allows to investigate the effect of weak intermolecular interactions, and this technique has been used to study some complexes that rare gases (Rg’s) form with aromatic molecules.7 In particular, the family of dimers containing benzene and rare gases (Bz-Rg), which can be considered to be prototype systems, have been detected by FTMW spectroscopy (except He-Bz).8 From a theoretical point of view, the modelization of the interaction and dynamic studies of heterogeneous clusters containing Bz is also a fruitful research area, both in gas and in condensed phase.3,4,9-16 The reliability of dynamical results, however, largely depends on the accuracy of the potential energy surface (PES). In general, the most powerful tools to investigate the PES characteristics are ab initio methods. Nevertheless, for some systems, the assessment of the relative role of the various components of the interaction, which are much weaker than those leading to chemical bonds, is very difficult and, accordingly, only the most relevant points (minima and saddle points) of the PES are calculated. In these cases, a proper combination of theoretical (ab initio) and experimental (spectroscopic and scattering) data often provides sufficient information to describe the PES by means of some suitable analytical functional form (semiempirical methods). (See, for instance, ref 17.) Despite * To whom correspondence should be addressed. E-mail: m.alberti@ ub.edu.

the limitations indicated above, ab initio calculations are of fundamental importance in theoretical investigations and, in the absence of experimental information, they provide the only results from which semiempirical predictions can be critically analyzed. Concerning Rg-Bz systems, a relatively large body of electronic structure calculations is available,18-21 however, because of the weakness of intermolecular interaction, an accurate construction of the related PESs is a difficult task. Intermolecular interaction in Rg-Bz clusters, because of the absence of a dipole moment in Bz and the isotropic atomic charge distribution, is dominated by weak dispersion forces including exchange-repulsion effects at short-range and induction and dispersion interaction contributions at long-range. This means that a proper modelization of the intermolecular potential energy needs to account for both short- and long-range effects. Accordingly, intermolecular interactions should be represented by functional forms capable of predicting the most important features of the system and able to provide reliable results in a wide range of configurations. Good prediction capabilities are particularly wanted when experimental investigations and ab initio calculations are not available. In molecular dynamics (MD) simulations, the interaction is usually described as a combination of a few contributions, expressed by a set of simple functions (force field). Because the most time-consuming part of almost all MD simulations is the calculation of the interaction, the simplicity of the functions is an important requirement for such simulations. Often, a model system with pairwise additive interactions is considered and, in this case, the well known Lennard-Jones function (LJ) is widely used. As it is well known, LJ exhibits a simple formulation and is able to reproduce very well both interatomic distances and energies at equilibrium. However, despite providing almost accurate descriptions at intermediate distances, LJ overestimates the strength of both the long-range attraction and the short-range repulsion. Bearing in mind that in MD studies atoms and molecules can explore different zones on the PES, the use of functional forms describing accurately the interaction, not only at intermediate distances but also in a wide range of configurations, seems to be absolutely necessary. A significant improvement of the long-range attraction description (respect

10.1021/jp9113927  2010 American Chemical Society Published on Web 01/28/2010

Rare Gas-Benzene-Rare Gas Interactions that of the LJ function) was proposed by Maitland and Smith.22 Unfortunately, the proposed modification also generates an excessive fall off of the repulsive wall. Following the basic idea given in ref 22, Pirani et al. proposed an ulterior modification of the potential function,23 which represents a further improvement on the description of the interaction at both short and long distances. The new ILJ potential function, which includes an additional parameter with respect to LJ, looks like LJ; for this reason, it has been called improved Lennard-Jones (ILJ). The additional parameter (see next section), the value of which depends on the nature of the interaction, allows us to modulate the hardness of the repulsive wall without modifying the values of the relevant parameters of the potential function (well depth and equilibrium distance), whose values can be easily predicted by scaling laws based on the use of some fundamental physical properties of the interacting partners.24 The new potential model was first applied to describe the force field (static behavior) of Rg interacting with hydrocarbon molecules (both aliphatic and aromatic)23 and then in MD simulations of Bz-Ar2 clusters (dynamic behavior).25 In the Rg2-Bz study, the ILJ function was only used to define the 12 atom bond (three-body) terms in which the Ar-Bz interaction was decomposed, whereas the atom-atom (two-body) term, describing the Ar-Ar interaction, was represented by means of an LJ function. Later, a detailed accuracy test of the new potential function, at first introduced to represent the intermolecular interaction between Rg and hydrocarbons,23 demonstrated that despite its simplicity, involving only one additional parameter with respect to LJ, the function is also able to describe a very large body of experimental data as scattering cross sections, transport properties, and spectroscopic information.26 The study of small clusters is very important to bridge the conceptual gap between single molecules and bulk materials.5 The large molecular surfaces of several aromatic compounds have many complexation sites, accepting several Rg solvent atoms. In particular, for Bz-Rgn systems, as indicated before, two general types of structures exist: one with Rg atoms distributed on both sides of the aromatic ring and one with Rg atoms placed on the same side. The clusters containing Bz and two Rg atoms are the simplest systems in which the two general types of structures can be observed. The existence of two conformers for systems containing an aromatic molecule and the same Rg’s as the prototype Bz-Ar2 have been investigated. (See, for instance, refs 27 and 28.) Clusters containing different Rg’s and organic molecules have also been detected.29 In the present study, our interest deals with the investigation of the structure and dynamic behavior of some Rg-Bz-Rg′ clusters, representing all interaction components (two- and three-body terms) by means of the ILJ function. All possible combinations of Rg are considered, except He. (The ILJ tested parameters only exist for He-Xe.) This means that the clusters that Bz forms with two identic Rg atoms (Rg-Bz-Rg) and with two different Rg atoms (Rg-Bz-Rg′) have been investigated. A molecular complex constituted of two Rg atoms attached to an aromatic molecule becomes more complicated from either a static or a dynamics point of view. However, the existence of two conformers, one with the two Rg atoms placed on the same side of the aromatic plane and the other with the two atoms placed on two different sides, adds interest to the study of these systems, allowing the investigation of isomerization processes.30 Moreover, from MD calculations, the time spent for the cluster in different conformations can be determined and related to the relative abundance of the two isomers. The consideration of the family of Rg’s also allows the investigation of size effects.

J. Phys. Chem. A, Vol. 114, No. 6, 2010 2267 The article is organized as follows: the main characteristics of the interaction potential are given in Section 2, the results are discussed in Section 3, and Section 4 concludes. 2. Potential Energy Function As in our previous studies on systems containing Bz,23,31-34 this molecule has been considered to be a rigid body, and intramolecular interactions are not evaluated. Moreover, because electrostatic interactions are absent in Bz-Rgn systems, the nonelectrostatic component accounts for the total intermolecular potential energy, V, which here is decomposed as a sum of Rg-Rg and Bz-Rg interaction contributions, VRgi-Rgj and VBz-Rgi, respectively n-1

V)

n

n

∑ ∑ VRg -Rg + ∑ VBz-Rg i)1 j>i

i

j

i)1

(1)

i

Concerning the Rg-Bz-Rg′ aggregates, with only two Rg atoms, V, takes the following form

V ) VRg-Rg′ + VBz-Rg + VBz-Rg′

(2)

As has been indicated before, all terms in eq 2 are described using ILJ functions. In particular, the Rg-Rg′ interaction is described using an unique ILJ function, which depends on only one variable (the interatomic distance, r), and it is represented as

[

VRg-Rg′ ) ε

()

r0 m n-m r

n

-

( )]

r0 n n-m r

m

(3)

where ε and r0 are the relevant parameters of the potential energy function (also used in LJ functions), the well depth, and the equilibrium distance, respectively. (m and n parameters are discussed later). The interaction potential energy between Rg and Bz (generalized as VBz-Rg) is decomposed as a sum of 12 atom-bond energy contributions. Accordingly, Rg-Bz interactions are described as 6

VBz-Rg )

∑

k)1

6

VRg-(CC)k +

∑ VRg-(CH)

k)1

k

(4)

where VRg-(CC)k and VRg-(CH)k, representing the interaction between a Rg atom and a given bond (CC and CH, respectively), are also described by means of ILJ functions. However, the ILJ function describing atom-bond interactions depends on two variables (Figure 1): r, which in this case indicates the distance of Rg from the center of a given bond, and γ, the angle that each r vector forms with a particular bond. Moreover, for atom-bond interactions, the well depth, ε, and the equilibrium distance, r0, are considered to be dependent on γ, and they are modulated through a simple trigonometric formula from the corresponding perpendicular (ε⊥, r0,⊥) and parallel (ε|, r0,|) components. (See, for instance, ref 31). This formulation allows us to account for the different characteristics of both the well depth and the equilibrium distance, at different atom-bond geometrical configurations. This means that for atom-bond interactions, the ILJ function depends on the r vector.

2268

J. Phys. Chem. A, Vol. 114, No. 6, 2010

Albertı´ TABLE 1: β Parameter, Well Depth, ε, and Equilibrium Distance, r0, for Rg-Rg Interaction Rg-Rg′

β

ε/meV

r0/Å

Ne-Ne Ne-Ar Ne-Kr Ne-Xe Ar-Ar Ar-Kr Ar-Xe Kr-Kr Kr-Xe Xe-Xe

9.0 9.0 9.0 9.0 9.0 9.0 9.0 9.0 9.0 9.0

3.66 5.74 6.16 6.35 12.37 14.33 16.09 17.30 19.95 24.20

3.094 3.520 3.660 3.885 3.757 3.910 4.100 4.010 4.200 4.350

TABLE 2: Rg Bond (Rg-CC and Rg-CH) Interaction Parametersa Figure 1. Variables of ILJ function: atom-atom (panel A) and atom-bond (panel B) interactions.

In the ILJ function, the n parameter, defining the falloff of the repulsion, is expressed as

()

n ) β + 4.0

r r0

2

(5)

For atom-atom interactions, n depends only on r (the distance between the two atoms). However, for atom-bond interactions, because of the angular dependence of r0, n depends not only on r (the distance between an atom and the center of a bond) but also on γ (the angle that r forms with a bond). The main advantage of this formulation is that n is no longer a constant and, accordingly, the repulsive wall can be modified depending on the interaction characteristics, by changing conveniently the value of β. Moreover, for a given value of β, n is modulated throught its dependence on r. These particular characteristics of the model add flexibility to the potential function with respect to LJ, allowing us to obtain good descriptions for a variety of interactions, for instance H bonds (small clusters and liquid water),35-37 neutral systems with partners interacting weakly,23,25 and compounds containing ions (positive, negative).26,31-34,38-40 Recently, the model has been extended to investigate the ternary cation-Bz-anion systems by considering the important induction effects originated by the presence of a second ion.41 Moreover, the model has also been extended to formulate bond-bond interactions.42 A detailed study performed to investigate the capability of ILJ describing different experiments26 suggested the use of lower values of β than those initially proposed for Rg-Bz systems.23,25 Accordingly, in the present study, the value of β has been lowered from 10, the value used in previous studies involving Bz and Rg,23,31-34 to 9. The value of m (see eq 3) has been taken to be equal to 6 (Rg-Rg and Rg-bond interactions), the typical value for neutral-neutral interactions. The remaining potential parameters are given in Tables 1 (Rg-Rg interactions)26 and 2 (Rg bond contributions). 2.1. Static Properties of Rg-Bz-Rg′ Clusters. Before the study of the static properties of Rg-Bz-Rg′ clusters, the results for Rg-Bz dimers are presented. Rg-Bz aggregates were widely investigated at experimental and theoretical levels.23,43,44 In these studies, the geometries for different approaches of Rg to Bz were analyzed, and theoretical predictions were compared with experimental and ab initio results. Because of the similarity of the geometry and the energy results obtained for Rg-Bz

Rg bond

β

ε⊥/meV

ε|/meV

r0,⊥/Å

r0,|/Å

Ne-CC Ne-CH Ar-CC Ar-CH Kr-CC Kr-CH Xe-CC Xe-CH

9.0 9.0 9.0 9.0 9.0 9.0 9.0 9.0

1.706 2.544 3.895 4.814 4.824 5.667 5.654 6.233

1.862 1.960 4.910 3.981 6.365 4.781 7.849 5.384

3.629 3.316 3.879 3.641 3.997 3.782 4.163 3.976

4.015 3.549 4.189 3.851 4.284 3.987 4.425 5.384

a Perpendicular and parallel components for the well depth (ε⊥,ε|) and equilibrium distance(r0,⊥,r0,|).

TABLE 3: Well Depth, De, and Equilibrium Distance of Rg from the Center of Mass of Bz, Re, for Three Relevant Cuts of the Rg-Bz PES out-of-plane

in-plane vertex out-of-plane vertex

Rg-Bz Re/Å De/meV Re/Å De/meV

Re/Å

De/meV

Ne-Bz

5.250 5.280 5.495 5.520 5.607 5.650 5.760 5.800

6.41 6.30 14.95 14.60 18.94 18.50 22.97 22.30

Ar-Bz Kr-Bz Xe-Bz

3.285 3.300 3.558 3.570 3.684 3.690 3.862 3.880

19.71 19.50 44.37 44.10 54.88 54.60 64.01 63.70

4.844 4.860 5.128 5.150 5.255 5.270 5.429 5.460

10.34 10.30 22.52 22.30 27.78 27.40 32.55 32.10

present ref 23 present ref 23 present ref 23 present ref 23

using β ) 10 and 9 (eq 4), in the present article (β ) 9), only selected geometries related to three relevant cuts of the PES are presented and compared with the previous ones (β ) 10). Results are shown in Table 3, where it can be seen that the most stable geometry is the one with Rg placed along the C6 symmetry axis of Bz. Moreover, the Table shows that results obtained by using the two values of β are quite similar, and here the theoretical results are not compared with the experimental ones because a detailed discussion about Rg-Bz dimers can be found in the bibliography.23 The Rg-Bz-Rg′ total intermolecular interaction, as shown in eq 2, can be described by means of two different interaction contributions, Rg-Rg′ and Rg-Bz (and Rg′-Bz). Both contributions are appreciable when Rg and Rg′ are placed on the same side of the Bz plane (isomers (2|0)). However, when Rg atoms are placed in opposite sides of the Bz plane (isomers (1|1)), the Rg-Rg′ interaction is almost negligible. Despite this, the (1|1) isomer is the most stable one because Rg and Rg′ are placed along the C6 symmetry axis of Bz (the most favorable geometry for Rg-Bz interactions).23 Geometries and energies, at the (1|1) equilibrium configuration, for all Rg-Bz-Rg′ clusters investigated are given in Table 4, whereas the geometries (in Cartesian coordinates) and energies, associated with the (2|0) equilibrium configuration, are given in Table 5. Calculations show that the (2|0) equilibrium structures depend on the Rg mass. It has been found that all Rg-Bz-Rg′ clusters

Rare Gas-Benzene-Rare Gas Interactions

J. Phys. Chem. A, Vol. 114, No. 6, 2010 2269

TABLE 4: Equilibrium Distances from the Atoms to Bz (R1, R2) Equilibrium Energy (Given with Respect to the Rg+Bz+Rg′ dissociation Limit: De), and Potential Energy Contribution of the Bond-Bond Interaction, VRg-Rg′ ((1|1) Equilibrium Configuration)a Rg-Bz-Rg′

R1/Å

R2/Å

De/meV

VRg-Rg′

Ne-Bz-Ne Ar-Bz-Ar Kr-Bz-Kr Xe-Bz-Xe Ne-Bz-Ar Ne-Bz-Kr Ne-Bz-Xe Ar-Bz-Kr Ar-Bz-Xe Kr-Bz-Xe

3.285 3.556 3.682 3.859 3.284 3.283 3.283 3.556 3.555 3.681

3.285 3.556 3.682 3.859 3.557 3.684 3.861 3.682 3.860 3.860

-39.45 -89.10 -110.37 -129.09 -64.23 -74.76 -83.94 -99.73 -109.01 -119.70

-0.05 -0.36 -0.62 -1.0 -0.17 -0.17 -0.22 -0.48 -0.63 -0.82

a For all Rg-Bz-Rg′ clusters, R1 and R2 represent the distances from Rg and Rg′ to the CM of benzene, respectively.

containing Ne show equilibrium configurations in which the heavier Rg (independently of its mass) is placed along the symmetry axis. Also, in the (2|0) equilibrium configuration of Ar-Bz-Ar, one of the atoms is placed along the symmetry axis of Bz. On the contrary, the remaining clusters have equilibrium structures with both atoms placed out of the symmetry axis. This can be explained from the increase in the Rg-Rg′ equilibrium distance with the mass (size) of Rg, which originates an increase in the equilibrium Rg-Bz center of mass (CM)-Rg′ angle. In particular, this angle is equal to approximately 50, 55, and 60° for Ne-Bz-Ar, Ar-Bz-Kr, and Kr-Bz-Xe, respectively. By increasing the mass of Rg, if the heavier atom is placed along the C6 axis, the lightest one is forced to occupy more repulsive zones of the PES. Consequently, when Rg atoms are heavy enough, the Rg-Bz-Rg′

equilibrium configuration is the one with both atoms placed out of the symmetry axis. In Figure 2, as a representative example of the different arrangements at the (2|0) configuration, the snapshots of the equilibrium geometry for Ne-Bz-Ar (on the left hide side (lhs)), Ar-Bz-Kr (on the center), and Kr-BzXe (on the right-hand size (rhs)) are shown. As can be seen in the Figure, when the mass of the cluster increases, the heavier atom is no longer placed on the symmetry axis, whereas the lightest atom tends to approach it. In Tables 4 and 5, I report the equilibrium energies of the (1|1) and (2|0) configurations, respectively. The two sets are very similar for all investigated clusters. Moreover, both stable isomers are connected by means of low energy barriers. The presence of low barriers interconnecting the minima has been observed for all clusters investigated. As an example, the equipotential contours for Ne-Bz-Ar are represented in Figure 3, where Bz and Ar (in red) have also been represented. Blue color has been used to identify the lowest energy contour (-62 meV). Consecutive contours are spaced by 20 meV. The existence of the (1|1) and (2|0) equilibrium configurations is apparent from the Figure, where it can also be seen that a low barrier connects both equilibrium configurations. 2.2. Molecular Dynamics Calculations. The DL_POLY code45 has been used to perform MD calculations on Rg-Bz-Rg′ clusters. The simulations, as indicated in Section 2, have been performed considering Bz as a rigid body. Taking into account the range of temperatures investigated in the present study, the rigidity of Bz should not affect dynamic results. To carry out the calculations, a Fortran routine, returning the value of the adopted PES and the values of the related derivatives with respect to the internuclear distances, was assembled. The initial conditions of the system were set to correspond to a microcanonical ensemble (NVE) in which the number of particles (N), volume (V), and total energy (Etotal) are conserved, and various

TABLE 5: Equilibrium Positions of the Atoms (xi, yi, zi), Equilibrium Geometry (De), and Potential Energy Contribution of the Bond-Bond Interaction, VRg-Rg′a Rg-Bz-Rg′

x1/Å

y1/Å

z1/Å

x2/Å

y2/Å

z2/Å

De/meV

VRg-Rg′′ /meV

Ne-Bz-Ne Ar-Bz-Ar Kr-Bz-Kr Xe-Bz-Xe Ne-Bz-Ar Ne-Bz-Kr Ne-Bz-Xe Ar-Bz-Kr Ar-Bz-Xe Kr-Bz-Xe

2.955 3.786 1.414 -0.586 3.294 3.371 3.445 3.447 3.447 1.414

0.000 0.000 2.449 0.339 0.000 0.000 0.000 0.000 0.000 2.449

2.955 2.651 3.371 3.841 2.763 2.828 2.219 2.893 2.893 3.371

0.000 0.000 0.000 3.677 0.000 0.000 0.000 -0.330 -0.519 0.245

0.000 0.000 -1.334 0.000 0.000 0.000 0.000 0.571 0.435 -1.381

3.291 4.622 3.665 3.085 3.566 3.696 3.855 3.742 3.841 3.853

-36.10 -82.85 -102.14 -125.48 -62.27 -72.25 -81.20 -94.31 -105.12 -115.07

-3.38 -11.88 -17.24 -24.20 -5.35 -5.36 -5.48 -14.33 -16.09 -19.95

a Benzene molecule is placed on the xy plane, with its CM on the origin of coordinates. The equilibrium energy is given with respect to the Rg+Bz+Rg′ dissociation limit ((2|0) equilibrium configuration). For all Rg-Bz-Rg′ clusters, indexes 1 and 2 are used to indicate the coordinates of Rg and Rg′, respectively.

Figure 2. Ne-Bz-Ar (lhs), Ar-Bz-Kr (center), and Kr-Bz-Xe (rhs) at the (2|0) configuration equilibrium. Red color is used to identify the heavier atom of each cluster, Ar (lhs), Kr (center), and Xe (rhs).

2270

J. Phys. Chem. A, Vol. 114, No. 6, 2010

Albertı´ TABLE 6: Mean Values of Temperature, T, Total Energy, Etotal, potential energy, Ecfg, Rg-Rg Interaction Energy, ERg-Rg, and Bz-Rg Interaction Energy, EBz-Rg, for Ar-Bz-Ar

Figure 3. Equipotential energy contours for the Ne-Bz-Ar system. The benzene molecule is placed on the xy plane, and the Ar atom (in red) is fixed at z ) 3.56 Å, whereas Ne moves on the xz plane. Blue color indicates the lowest energy contour of -62 meV, and consecutive contours are spaced by 20 meV.

values of Etotal were considered. Calculations are performed without imposing any boundary conditions. At each new energy, calculations were started by taking as the initial configuration of velocities and forces the final values of the previous run. A time step of 1 fs was adopted for the integration of the motion equations. At each step, the total potential energy is calculated (eq 2), allowing the calculation of the kinetic energy, Ek, i. Then, the instantaneous temperature, Ti, was calculated at each time step from Ek, i by means of the relationship Ti ) 2Ek, i/kBf (where kB is the Boltzmann constant and f is the number of degrees of freedom in the system). This means that each step has assigned a value of potential energy, Ti and Ek, i. The total integration time for the simulations was generally set at 50 ns. At the end of the simulation, mean values of temperature, T, kinetic energy, Ekin, and potential energy, Ecfg, are obtained. Accordingly, with eq 2, Ecfg has contributions of the Rg-Rg′, Bz-Rg, and Bz-Rg′ interactions (ERg-Rg′, EBz-Rg, and EBz-Rg′). 3. Results and Discussion In general, it has been observed that all clusters perform frequent configuration inconversions, even at low temperatures. In particular, the lighter clusters, governed by weaker interactions than the heavier ones, isomerizate at lower temperatures. By increasing T (below the temperature value from which interconversions are possible), the atoms explore different zones of the PES, as far from equilibrium configurations as high as the temperature is. This means that the aggregate, at sufficient high temperatures, can adopt rather different instantaneous (at each step) geometries having quite different values of potential energy. Consequently, in the context of NVE, ensemble of particles, in which the total energy (kinetic plus potential energy) is conserved, different instantaneous values of the kinetic energy are also obtained. This originates noticeable fluctuations on Ti along the dynamic simulation. These fluctuations should indicate that the simulation conditions are close to those from which isomerization processes can occur. Accordingly, a subsequent small increase in Etotal, causing a small increase in T, will originate the beginning of the isomerization process, which for Rg-Bz-Rg′ aggregates is also characterized by appreciable variations of ERg-Rg′. When an (1|1) initial configuration is considered, the isomerization process to (2|0) configurations originates a noticeable enhancement of the interaction between the two atoms and, consequently, a decrease in ERg-Rg′ (more negative). Contrariously, for (2|0) initial configurations of a given aggregate, an increase in ERg-Rg′ can be related to the corresponding isomerization process. For a given cluster, the number of interconversions depends on T, and, in general, at the lower values of T at which the

T/K

Etotal/meV

Ecfg/meV

ERg-Rg/meV

EBz-Rg/meV

5.0 ( 0.2 15.1 ( 0.3 25.0 ( 0.4 30.5 ( 12.8 39.8 ( 16.2

-86.41 -80.80 -75.10 -71.09 -63.35

-87.70 -84.80 -82.33 -79.00 -73.64

-0.36 -0.35 -0.34 -0.31 -7.75

-87.36 -84.45 -81.97 -78.69 -65.89

cluster isomerizates, only a few interconversions are observed. However, once the isomerization process begins, enhancements on the interconversion frequency are observed, even for small increases in T. Further increments of T, however, favor the dissociation of aggregates (Bz-Rg+Rg′). This means that Rg-Bz-Rg′ clusters can isomerize only in a short range of temperatures. Out of this range, at lower T values, the isomerization is not possible because there is not enough energy to surmount the barrier connecting (1|1) and (2|0) structures, whereas at higher values of T, the aggregate tends to dissociate. The values of T at which isomerization occurs are the expected ones from the expansion of a supersonic jet, originating the formation of clusters.The isomer population (at the end of the colliding process) results from those clusters formed at intermediate temperatures (T ≈ 40 K), for which interconversions are quite likely.27,46 The criterion adopted to decide whether an interconversion between the two isomers has taken place is the same as that used in other dynamical studies of the system:25,27 the isomer is required to remain around a particular geometry for at least 5 ps. Moreover, a particular Rg-Bz-Rg′ geometry is associated with either the (2|0) or the (1|1) isomer only when the distances of the Rg atoms from Bz CM are >1.5 Å. 3.1. Rg-Bz-Rg Clusters. Interconversions between (1|1) and (2|0) configurations have been observed at temperatures of about 18, 35, 50, and 65 K for Ne-Bz-Ne, Ar-Bz-Ar, Kr-Bz-Kr, and Xe-Bz-Xe, respectively. By increasing T, the number of interconversions increases. However, for all investigated clusters, an increase in T higher than 10 K reverts on a noticeable enhancement of the dissociation probability. This means that isomerizations can be observed only in a short range of temperatures. Some results of dynamic simulations for the Ar-Bz-Ar aggregate, performed at different values of T, are shown in Table 6, where it can be seen that for temperatures ranging from 5 to 30 K the Rg-Rg interaction is very small (values of ERg-Rg close to zero). However, a remarkable decrease in ERg-Rg is observed at T ) 39.8 K. The selected initial configuration for this simulation has been the (1|1) configuration, and the decrease in ERg-Rg observed at T ) 39.8 K indicates a diminution of the mean Rg-Rg′ interatomic distance, which can be associated with the formation of (2|0) configurations. Moreover, if one bears in mind that the minimum potential energy interaction between Ar atoms is (in absolute value) equal to 12.37 meV (see Table 1), the result of ERg-Rg ) -7.75 meV obtained at T ) 39.8 K should indicate that interconversions between (1|1) and (2|0) configurations occur. Table 6 also shows that at the mean temperature of 30.5 K, just before the beginning of isomerizations, the temperature associated with each step fluctuates noticeably along the trajectory (high fluctuations in kinetic energy and, consequently, in potential energy). By further analyzing the Rg-Rg interaction, it has been observed that an initial increase in T originates an enhancement

Rare Gas-Benzene-Rare Gas Interactions

J. Phys. Chem. A, Vol. 114, No. 6, 2010 2271

Figure 4. Time evolution of the Xe-Xe distance on the Xe-Bz-Xe aggregate under the T conditions indicated in the Figure. Top panel: the initial configuration is the one with the two atoms placed on the same side of the benzene plane. Lower panel: the initial configuration is the one with the two atoms placed on opposite sides of the benzene plane.

of the interaction energy. This can be explained by considering the fact that the initial configuration of the cluster corresponds to molecular structures close to equilibrium and that small increases of T allow the system to explore more repulsive zones of the PES. In particular, further increases in T (just before isomerization) allow the atoms to explore zones close to the Bz plane. In this case, often, one of the atoms remains in the preferred geometry (along the symmetry axis of Bz), and the distance between Rg atoms increases. This originates a diminution of the interatomic interaction and, consequently, a small increase in ERg-Rg, which can be observed before its important decrease associated with the isomerization process. Time evolution of the Rg-Rg distance is analyzed in Figure 4 for Xe-Bz-Xe at two similar temperatures. MD calculations performed at constant total energies of Etotal ) -90.60 meV (T ) 66. 0 ( 27.0 K) (top panel, (2|0) initial configuration) and of Etotal ) -91.60 meV (T ) 64.4 ( 26.1 K) (lower panel, (1|1)initial configuration) show that the time spent for Xe-Bz-Xe around a particular geometry, (2|0) or (1|1), is very similar in both cases. MD calculations, at two similar temperatures but different initial configurations, have also been performed for the remaining Rg-Bz-Rg aggregates. The results allow us to conclude that the time spent for a given aggregate around some particular geometry depends on T but does not depend on the initial configuration. In Figure 5, the effect of increasing T is analyzed for the Kr-Bz-Kr system. At the lower temperature (top panel), the system, initially at a (1|1) configuration, remains around this configuration (without isomerization) for ∼5 ns; then, (1|1) Kr-Bz-Kr isomerizates. Other isomerizations have not been observed along the trajectory. The ratio of the time spent around the (2|0) and the (1|1) configuration, t(2|0)/t(1|1), at T ) 49.6 ( 19.6 K is equal to 8.16. At higher temperatures, the frequency of the interconversions increases, as can be seen in the medium and lower panel of the Figure. Moreover, the t(2|0)/t(1|1) ratio diminishes with T, and at temperatures of 56.4 ( 23.14 and 64.4 ( 27.5 K, t(2|0)/t(1|1) ratios of 5.23 and 3.96 have been obtained, respectively. For all clusters, it has been found that interconversions privilege the (2|0) isomer over the (1|1) isomer, probably because the (2|0) isomer explores a large configuration space. The values

Figure 5. Time evolution of the Kr-Kr distance on the Kr-Bz-Kr at three different temperatures. Top and medium panel: the initial configuration is the one with the two atoms placed on opposite sides of the benzene plane. Lower panel the initial configuration is the one with the two atoms placed on the same side of the benzene plane.

of t(2|0)/t(1|1) can be taken to be proportional to the relative isomeric abundance. However, as it has been pointed out elsewhere,46 converged estimates can be obtained only for long simulation runs. In the previous work on Ar clusters, it was found that 40 ns of simulation suffices to reach converged estimates. Nevertheless, to ensure that converged results are obtained, a simulation time of 50 ns is long enough; longer simulations of 150 ns have been performed for some systems under some conditions. No differences have been found accumulating results during 50 or 150 ns. 3.2. Rg-Bz-Rg′ Clusters. For Rg-Bz-Rg′ aggregates, the enhancement of kinetic energy, originated by an increase in T, mainly affects the mobility of the lightest Rg, which, when the isomerization process takes place, moves from one side to the other of the Bz plane. Contrariously, the heavier atom remains placed close to its equilibrium position. This means that the Rg-Bz-Rg′ dynamic behavior is mainly governed by the movement of the lightest Rg; consequently, those aggregates containing the same light atom but different heavy Rg should begin the isomerization process at the same temperature. This, in fact, has been observed from MD calculations, whose results show that the temperature at which the isomerization process begin depends on the nature of the lightest atom on the cluster but is almost independent of the heavier one. Accordingly, interconversions between the (1|1) and (2|0) configurations for Ne-Bz-Ar, Ne-Bz-Kr, and Ne-Bz-Xe have been observed at a temperature of ∼18 K (as for Ne-Bz-Ne). However, some differences on the mobility of the heavier atoms are observed when comparing the dynamic behavior of those clusters having (2|0) equilibrium geometries with the heavier atom placed along the symmetry axis, as, for instance, the Ne-Bz-Rg series, with those with the heavier atom placed out of the symmetry axis, as, for instance, Ar-Bz-Kr, Ar-Bz-Xe, and Kr-Bz-Xe. (See

2272

J. Phys. Chem. A, Vol. 114, No. 6, 2010

Albertı´

Figure 6. Ne-Bz-Ar caloric curve.

Figure 8. Time evolution of the distances of Xe (upper panel) and Ne (lower) from the center of mass of benzene before the isomerization process. The initial configuration is the one with Xe and Ne placed on opposite sides of the benzene plane.

Figure 7. Time evolution of the distances between Ar and Ne for different values of mean temperature, T. Top and medium panels: in the initial configuration, Ar and Ne are placed on opposite sides of the benzene plane. Lower panel: the initial configuration is the one with Ar and Ne placed on the same side of the benzene plane.

Table 5.) It has been observed that the mobility of the heavier atom is lower when it is placed along the symmetry axis. This can be explained from the higher stability of the Rg-Bz dimers when the atom is placed along the symmetry axis. The dynamic study for the Ne-Bz-Rg family shows that for Ne-Bz-Ar, the beginning of the isomerization processes occurs at Etotal of about -53 meV, obtaining after 50 ns of simulation a mean temperature of 17.642 K. From this value, small increases in total energy originate an enhancement of the isomerization processes, which take place at almost constant kinetic energy, Ekin (and mean temperature). This behavior is shown in the caloric curve represented in Figure 6, which shows a plateau at temperatures of ∼18 K. After a very short interval of T, further increases in Etotal revert on an increase of Ekin and to an enhancement of the isomerization process, as shown in Figure 7, where the evolution of the distance from Ne to Ar is represented for three different values of T. In general, because of the dynamic behavior of those clusters having the same lightest atom (Ne-Bz-Kr and Ne-Bz-Xe), no noticeable differences with respect to Ne-Bz-Ar have been observed. To study the size effect, the Rg-Bz-Xe family (Ne-Bz-Xe, Ar-Bz-Xe, and Kr-Bz-Xe) has also been investigated. Taking into account the previous results, it can be expected that the T at which the isomerization of Ne-Bz-Xe begins is very

Figure 9. Time evolution of the distances of Xe (upper panel) and Ne (medium panel) from the center of mass of benzene. The evolution of the interatomic distance between Rg atoms is shown in the lower panel. The initial configuration is the one with Xe and Ne placed on the same side of the benzene plane.

similar to that of Ne-Bz-Rg (Rg ) Ne, Ar, Kr). Similarly, the lowest T for the isomerization of Ar-Bz-Xe should be very similar to that of Ar-Bz-Rg (Rg ) Ar, Kr) and so on. In fact, it has been observed that the temperature at which the isomerization begins is 18, 35, and 50 K for Ne-Bz-Xe, Ar-Bz-Xe, and Kr-Bz-Xe, respectively. The higher mobility of the lightest atom (even before the isomerization process) in comparison with the heavier one is shown in Figure 8 for Ne-Bz-Xe at a mean temperature value of 12.945 K. For the same aggregate at a mean temperature of 17.268 K (the lowest T at which isomerizations are observed), Figure 9 shows the evolution of the distances of Xe and Ne from the benzene CM (top and medium panels, respectively) as well as the interatomic distance between Xe and Ne (lower panel). The difference of the mobility of Xe and Ne can be observed by comparing the distance evolutions shown on the higher and medium panels of

Rare Gas-Benzene-Rare Gas Interactions

J. Phys. Chem. A, Vol. 114, No. 6, 2010 2273 ferent results. To obtain good results, independent of the simulation, a small increase in Etotal is required. Then, the relative abundance tends to a certain value that is independent of the initial configuration of the cluster ((1|1) or (2|0)) and of the number of interconversions. In all cases, it has been found, as before, that interconversions privilege the (2|0) isomer over the (1|1) one, mainly for the heavier clusters, accordingly with the higher stabilization of the Rg-Rg′ interaction when the mass increases. MD simulations give (2|0)/(1|1) ratios of about 1.2, 2.3, and 3.0 for Ne-Bz-Xe, Ar-Bz-Xe, and Kr-Bz-Xe trimers, respectively. 4. Conclusions

Figure 10. Time evolution of the distances of Xe (upper panel) and Ar (lower panel) from the center of mass of benzene for the Ar-Bz-Xe cluster. In the initial configurations, both Rg atoms are placed on the same side of the benzene plane.

Figure 9. Moreover, the lower panel shows that Ne-Bz-Xe, at the beginning of the trajectory, is in the (2|0) configuration and that it isomerizates after ∼8 ns, remaining in (1|1) configurations for ∼10 ns. After that, the cluster isomerizates again to (2|0) configurations. Dynamics calculations for the Ar-Bz-Xe and Kr-Bz-Xe clusters show that when the heavier atom in the (2|0) equilibrium configuration of the cluster is not placed along the symmetry axis, the mobilities of the light and the heavier Rg are not so different. However, also for Rg-Bz-Xe family of clusters, the mobility of the lightest Rg is higher, for instance, at Etotal ) -81.62 meV (44.833 K of mean temperature), the distances of Ar from the CM of benzene in the Ar-Bz-Xe cluster oscillate from 3.22 to 6.50 Å, whereas those of Xe from the CM of benzene oscillate from 3.55 to 5.60 Å. This can be observed in Figure 10, where the evolution of the distances of Xe and Ar to the CM of benzene are shown. In regards to the temperature range in which isomerizations can be observed, molecular dynamic simulations have shown that it depends on the cluster mass. For instance, for Ne-Bz-Xe clusters, interconversions are only observed in a range of ∼10 K, whereas the range of T in which isomerizations are observed is of about 30 and 40 K for Ar-Bz-Xe and Kr-Bz-Xe, respectively. This can be explained because of the higher stability of the heavier clusters. The stability of the Rg-Bz-Rg′ clusters increases with the Rg’s mass. This means that a same increase in Etotal (causing similar increases in T) cannot affect the stabilization of the heavier clusters but can increase the dissociation probability of the lightest ones. In the context of the NVE ensemble, for which total energy is conserved, an increase in Etotal allows the cluster to explore more repulsive zones of the PES. This means that small increases in total energy (to a maximum of a few millielectronvolts) originate enhancements of the interconversion number. Moreover, when the cluster explores repulsive zones of the PES, the dissociation probability increases, mainly for the lightest clusters that have lower equilibrium energies than the heavier ones. The relative isomeric abundance has also been calculated for Ne-Bz-Xe, Ar-Bz-Xe, and Kr-Bz-Xe. It has been observed that at the lower isomerization temperature with a very short number of interconversions, different simulations can predict different relative isomeric abundances. Sometimes, at very low values of Etotal (and of T), only one or two interconversions are observed. Accordingly, similar simulations can predict very dif-

Structural and dynamic properties of the symmetric RgBz-Rg and asymmetric Rg-Bz-Rg′ trimers have been analyzed. The ILJ modified version of the LJ potential has been used to calculate the interaction potential. The potential energy has been calculated by considering two- and three-body (atom-bond) interactions. All trimers have two equilibrium structures, the most stable being the one with the two Rg atoms placed on opposite sides of the benzene plane, (1|1). In the second structure, the two atoms are placed on the same side of the benzene plane, (2|0). The corresponding (2|0) equilibrium structure presents some differences as a function of the cluster mass. The two lighter symmetric trimers, Ne-Bz-Ne and Ar-Bz-Ar, have (2|0) equilibrium structures in which one Rg atom is placed along the symmetry axis of benzene. Contrariously, the heaviers Kr-Bz-Kr and Xe-Bz-Xe have (2|0) equilibrium structures in which Rg atoms are not placed on the symmetry axis of benzene. In regards to the asymmetric Rg-Bz-Rg′, in the (2|0) equilibrium structures, all trimers containing Ne have the heavier Rg placed along the symmetry axis. The increase in the mass of the lightest atom affects the equilibrium position of the heavier atom on the trimer, and all of the remaining trimers have (2|0) equilibrium structures without atoms placed along the benzene symmetry axis. MD calculations show that before the beginning of the isomerization, the cluster remains in the initial configuration ((2|0) or (1|1)) without changes. This behavior can be observed for a small range of total energy values (and T). This range depends on the mass of the cluster being smaller for the lightest clusters. The Rg mass affects their mobility in the isomerization processes because on one hand the atoms placed on the symmetry axis hardly interact with benzene, and on the other hand, the heavier Rg’s interact with each other stronger than the lighter ones. This originates that the relative abundance of the isomers depends on the mass combination of Rg’s on the clusters. For a given isomer, the relative abundance of the (2|0) one is always higher than that of (1|1), according to the existence of more equivalent isomers for the (2|0) conformer (twelve, six on each side of the Bz plane) than for the (1|1) (only two equivalent isomers). Acknowledgment. The author acknowledges financial support from the Ministerio de Educacio´n y Ciencia (Spain, project CTQ2007-61109), the Ministerio de Ciencia e Innovacio´n (Spain, Mobility Program, project PR2008-0251), and the Generalitat de Catalunya-DURSI (project 2005 PEIR 0051/69). Thanks are also due to the Centre de Supercomputacio´ de Catalunya CESCA-C4 and Fundacio´ Catalana per a la Recerca for the allocated supercomputing time. References and Notes (1) Mu¨ller-Dethelfs, K.; Hobza, P. Chem. ReV. 2000, 100, 143. (2) Alonso, J. L.; Antolı´nez, S.; Bianco, S.; Lesarri, A.; Lon´pez, J. C.; Caminati, W. J. Am. Chem. Soc. 2004, 126, 3244.

2274

J. Phys. Chem. A, Vol. 114, No. 6, 2010

(3) Brupbacher, Th.; Makarewicz, J.; Bauder, A. J. Chem. Phys. 1994, 101, 9736. (4) Schmidt, M.; Mons, M.; Le Calve´, J.; Millie´, P. Chem. Phys. Lett. 1991, 177, 371. (5) Adams, J. E.; Stratt, R. M. J. Chem. Phys. 1990, 93, 1358. (6) Hahn, M. Y.; Whetten, R. Phys. ReV. Lett. 1988, 61, 1190. (7) Caminati, W.; Grabow, J. U. Microwave Spectroscopy: Molecular Systems. In Frontiers of Molecular Spectroscopy; Laane, J., Ed.; Elsevier: New York, 2008; Chapt.15, p 508. (8) Evangelisti, L.; Favero, L. B.; Giuliano, B. M.; Tang, S.; Melandri, S.; Caminati, W. J. Phys. Chem. A 2009, 113, 14227. (9) Ondrechen, M. J.; Berkovitch-Yellin, Z.; Jortner, J. J. Am. Chem. Soc. 1981, 103, 6586. (10) Fried, L. E.; Mukamel, S. J. Chem. Phys. 1991, 96, 116. (11) Schmidt, M.; Le Calve`, J.; Mons, M. J. Chem. Phys. 1993, 98, 6102, and references therein. (12) Hobza, P.; Bludsky, O.; Selzle, H. L.; Schlag, E. W. Chem. Phys. Lett. 1996, 250, 402. (13) Lenzer, T.; Luther, K. J. Chem. Phys. 1996, 105, 10944. (14) Dullweber, A.; Hodges, M. P.; Wales, D. J. J. Chem. Phys. 1998, 106, 1530. (15) Mons, M.; Courty, A.; Schmidt, M.; Le Calve`, J.; Piuzzi, F.; Dimicoli, I. J. Chem. Phys. 1997, 106, 1676. (16) Bernshtein, V.; Oref, I. J. Chem. Phys. 2000, 112, 686. (17) Roncero, O.; Villarreal, P.; Delgado-Barrio, G.; Gonza´lez-Platas, J.; Breto´n, J. J. Chem. Phys. 1998, 109, 9288. (18) Hobza, P.; Selzle, H.; Schlag, E. W. Chem. ReV. 1994, 94, 1767. (19) Klopper, W.; Leuthi, H. P.; Brupbacher, Th.; Bauder, A. J. Chem. Phys. 1994, 101, 9747. (20) Koch, H.; Ferna´ndez, B.; Makariewicz, J. J. Chem. Phys. 1998, 108, 2784. (21) Koch, H.; Ferna´ndez, B.; Makariewicz, J. J. Chem. Phys. 1999, 111, 198. (22) Maitland, G. C.; Smith, E. B. Chem. Phys. Lett. 1973, 22, 443. (23) Pirani, F.; Albertı´, M.; Castro, A.; Moix Teixidor, M.; Cappelletti, D. Chem. Phys. Lett. 2004, 394, 37. (24) Pirani, F.; Maciel, G. S.; Cappelletti, D.; Aquilanti, V. Int. ReV. Phys. Chem. 2006, 25, 165. (25) Albertı´, M.; Castro, A.; Lagana, A.; Pirani, F.; Porrini, M.; Cappelletti, D. Chem. Phys. Lett. 2004, 392, 514. (26) Pirani, P.; Brizi, S.; Roncaratti, L. F.; Casavecchia, P.; Cappelletti, D.; Vecchiocattivi, F. Phys. Chem. Chem. Phys. 2008, 10, 5489.

Albertı´ (27) Vacek, J.; Konvicka, K.; Hobza, P. Chem. Phys. Lett. 1994, 220, 85. (28) Schmidt, M.; Mons, M.; Le Calve´, J.; Millie´, P.; Cossart-Magos, C. Chem. Phys. Lett. 1991, 183, 69. (29) Melandri, S.; Giuliano, B. M.; Maris, A.; Evangelisti, L.; Velino, B.; Caminati, W. ChemPhysChem 2009, 2503. (30) Ben-Horin, N.; Even, U.; Jortner, J. Chem. Phys. Lett. 1992, 188, 73. (31) Albertı´, M.; Castro, A.; Lagana`, A.; Moix, M.; Pirani, F.; Cappelletti, D.; Liuti, G. J. Phys. Chem. A 2005, 110, 9002. (32) Albertı´, M.; Aguilar, A.; Lucas, J. M.; Lagana`, A.; Pirani, F. J. Phys. Chem. A 2007, 111, 1780. (33) Huarte-Larragan˜a, F.; Aguilar, A.; Lucas, J. M.; Albertı´, M. J. Phys. Chem. A 2007, 111, 8072. (34) Albertı´, M.; Aguilar, A.; Lucas, J. M.; Cappelletti, D.; Lagana`, A.; Pirani, F. Chem. Phys. 2006, 328, 221. (35) Albertı´, M.; Aguilar, A.; Cappelletti, D.; Lagana`, A.; Pirani, F. Int. J. Mass. Spectrom. 2009, 280, 50. (36) Albertı´, M.; Aguilar, A.; Bartolomei, M.; Cappelletti, D.; Lagana`, A.; Lucas, J. M.; Pirani, F. Phys. Script. 2008, 78, 058108. (37) Faginas Lago, N.; Huarte-Larragan˜a, F.; Albertı´, M. Eur. Phys. J. D 2009, 55, 75. (38) Albertı´, M.; Castro, A.; Lagana`, A.; Moix, M,; Pirani, F.; Cappelletti, D. Eur. Phys. J. D 2006, 38, 185. (39) Albertı´, M.; Aguilar, A.; Lucas, J. M.; Pirani, F. Theor. Chem. Acc. 2009, 123, 21. (40) Cappelletti, D.; Bartolomei, D.; Pirani, F.; Aquilanti, V. J. Phys. Chem. A 2002, 106, 10764. (41) Albertı´, M.; Aguilar, A.; Pirani, F. J. Phys. Chem. A 2009, 113, 14741. (42) Thibault, F.; Cappelletti, D.; Pirani, F.; Bartolomei, M. J. Phys. Chem. A 2009, 113, 14867. (43) Albertı´, M.; Aguilar, A.; Lucas, J. M.; Pirani, F.; Coletti, C.; Re, N. J. Phys. Chem. A 2009, 113, 14606. (44) Pirani, F.; Porrini, M.; Cavalli, S.; Bartolomei, M.; Cappelletti, D. Chem. Phys. Lett. 2003, 367, 405. (45) The DL_POLY Molecular Simulation Package. http://www. cse.clrc.ac.uk/ccg/software/DL_POLY/index.shtml. (46) Vacek, J.; HoBza, P. J. Phys. Chem. 1994, 98, 11034.

JP9113927