Water Clusters in an Argon Matrix: Infrared Spectra from Molecular

Feb 4, 2015 - ... Density Functional-Based Tight Binding/Force-Field Potential. Aude Simon†, Christophe Iftner†, Joëlle Mascetti‡, and Fernand ...
0 downloads 0 Views 5MB Size
Article pubs.acs.org/JPCA

Water Clusters in an Argon Matrix: Infrared Spectra from Molecular Dynamics Simulations with a Self-Consistent Charge Density Functional-Based Tight Binding/Force-Field Potential Aude Simon,*,† Christophe Iftner,† Joel̈ le Mascetti,‡ and Fernand Spiegelman† †

Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, Université de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France ‡ Institut des Sciences Moléculaires, Université de Bordeaux and CNRS, 351 Cours de la Libération, 33405 Talence cedex, France S Supporting Information *

ABSTRACT: The present theoretical study aims at investigating the effects of an argon matrix on the structures, energetics, dynamics, and infrared (IR) spectra of small water clusters (H2O)n (n = 1−6). The potential energy surface is obtained from a hybrid selfconsistent charge density functional-based tight binding/force-field approach (SCCDFTB/FF) in which the water clusters are treated at the SCC-DFTB level and the matrix is modeled at the FF level by a cluster consisting of ∼340 Ar atoms with a face centered cubic (fcc) structure, namely (H2O)n/Ar. With respect to a pure FF scheme, this allows a quantum description of the molecular system embedded in the matrix, along with all-atom geometry optimization and molecular dynamics (MD) simulations of the (H2O)n/Ar system. Finite-temperature IR spectra are derived from the MD simulations. The SCC-DFTB/FF scheme is first benchmarked on (H2O)Arn clusters against correlated wave function results and DFT calculations performed in the present work, and against FF data available in the literature. Regarding (H2O)n/Ar systems, the geometries of the water clusters are found to adapt to the fcc environment, possibly leading to intermolecular distortion and matrix perturbation. Several energetical quantities are estimated to characterize the water clusters in the matrix. In the particular case of the water hexamer, substitution and insertion energies for the prism, bag, and cage are found to be lower than that for the 6-member ring isomer. Finite-temperature MD simulations show that the water monomer has a quasifree rotation motion at 13 K, in agreement with experimental data. In the case of the water dimer, the only large-amplitude motion is a distortion−rotation intermolecular motion, whereas only vibration motions around the nuclei equilibrium positions are observed for clusters with larger sizes. Regarding the IR spectra, we find that the matrix environment leads to redshifts of the stretching modes and almost no shift of the bending modes. This is in agreement with experimental data. Furthermore, in the case of the water monomer and dimer, the magnitudes of the computed shifts are in fair agreement with the experimental values. The complex case of the water hexamer, which presents several low-energy isomers, is discussed.

1. INTRODUCTION Experiments in a cryogenic environment allow the determination of the spectroscopic signatures of trapped molecular species. In particular, water clusters (H2O)n (n = 1−6) in matrices have been the object of various studies.1−10 Assigning the infrared (IR) bands to a specific structure may be delicate. First, clusters with several sizes may be embedded in the matrix, their relative abundances being likely to vary as a function of the water:matrix ratio and of the temperature. Second, several isomers for a given cluster may form. This is the case for the water hexamer, and the identification of the formed isomer(s) has received a great deal of interest lately.6,10 The interpretation of the IR spectra may thus become tedious, and theoretical support is often needed for an attempt to assign the observed bands. For instance, the experimental work by Ceponkus and co-workers has been complemented with DFT (density functional theory)6,7 or wave function8 calculations of harmonic © XXXX American Chemical Society

IR spectra. In systems such as water clusters, however, and as mentionned above, finite-temperature effects are likely to have an impact on geometries, energetics, and IR spectra. Theoretical investigations of anharmonic effects are quite advanced for the water dimer:11,12 starting from a fully flexible 12D potential, several effective 6D potentials have been used as supports for quantum dynamics calculations, leading to an accurate determination of the vibration−rotation−tunneling (VRT) levels involving the intermolecular rovibrational and tunneling states as well as the intramolecular vibrations. Such calculations were used to interpret the vibrational spectra of the water dimer5 or trimer13 in a neon matrix over a wide frequency Special Issue: Markku Räsänen Festschrift Received: August 23, 2014 Revised: February 2, 2015

A

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A range [90 to 14 000 or 11 000 cm−1, respectively]. More generally, several theoretical papers based on correlated wave function calculations have recently been dedicated to geometry optimization and vibrational frequency determination of small water clusters, unraveling anharmonic14,15 and zero-point vibrational energy effects.16 However, these advanced and accurate calculations are restricted to small systems, and the matrix environment is not taken into account. The effects of the matrix on the positions of the IR bands of a molecular system are usually small shifts7 of a few reciprocal centimeters (up to ∼ 10 cm−1). The theoretical description of these shifts is a challenge because the use of a technique allowing for the description of a large system (the impurity plus the matrix environment) is required, which would be accurate enough to describe such detailed and small spectroscopic features. It is almost impossible to describe realistic systems consisting of molecular clusters in matrices with purely ab initio calculations, unless describing a limited number of Ar atoms, for instance the 12 closest Ar neighbors.17 Several alternatives do exist, generally to treat small or simple systems, as the diatomics-in-molecules18,19 approach or equivalent schemes for chromophores such as atoms.20−23 Another possibility is simplifying the electronic structure of rare gases by substituting the whole rare gas atom by a pseudopotential, as is done for instance to treat Na clusters embedded in an Ar matrix.24−26 In the models used for these studies, the active electrons of the embedded clusters are described explicitely, while the surrounding matrix is described in terms of atomic pseudopotentials and core-polarization potentials, following an older work by Gross and Spiegelman.27 Gervais et al.24 were thus able to determine the optical absorption spectra of embedded clusters considering various possible trapping sites. The same group recently investigated the visible absorption and related luminescence of alkali atoms (Li, Na, and K) embedded in an Ar matrix, which allowed them to determine accurately the gas-to-matrix shifts for various trapping sites.28 Until now, such methods have been either restricted to rather simple systems (atoms, small molecules) or limited in the possibility of achieving exhaustive simulations (global search of the isomers, long-time molecular dynamics) when the determination of the potential energy surface (PES) is too costly (wave function schemes and even DFT approaches). In the present paper, we propose to use a hybrid self-consistent-charge density functional based tight binding29/force-field (SCC-DFTB/FF) scheme that we recently developed and benchmarked on the structures and energetics of (C6H6)0/+Arn clusters,30 with the aim of describing small water clusters (H2O)n (n = 1−6) embedded in an Ar matrix, namely (H2O)n/Ar. The water clusters are treated at the SCC-DFTB level of theory, an approximate DFT scheme that has been adapted for a correct description of structures and energetics of water clusters.31−34 Using this modified SCC-DFTB Hamiltonian, the geometries and energetics of the (H2O)6 isomers in particular were found to be in good agreement with correlated wave function calculations33 and with experimental data35 regarding the geometries. The Ar matrix is described using a FF approach, and the SCC-DFTB/FF coupling is incorporated via a scheme close to first-order degenerate perturbation theory in the linear combination of atomic orbitals (LCAO) basis supporting DFTB.30 Unlike a pure FF scheme, this hybrid SCC-DFTB/FF scheme is quantum for the embedded system and makes no assumption about electron localization in water clusters. It is also efficient enough to make possible all-atom geometry

optimizations and molecular dynamics (MD) simulations of the total (H2O)n/Ar system. The efficiency of the electronic structure calculation allows MD simulations of ∼200 ps that, repeated ∼10 times, are expected to ensure a fair statistical representation of the events, at least in the neighborhood of the trapped molecule. We thus propose in this work to theoretically investigate the effects of an Ar matrix on the structures, energetics, and mid-IR spectra of small water clusters (H2O)n (n = 1−6) with the (H2O)n/Ar system PES described within the SCC-DFTB/FF scheme. The equations of the model and the determination of the O−Ar analytical parametrized interaction functions are briefly summarized in section 2.1, followed by the description of the optimization techniques (section 2.2) and of the MD simulations from which finite-temperature IR spectra are derived (section 2.3). Benchmarks of the model for (H2O)Arn clusters are shown in section 3.1, while section 3.2 presents the structural and energetical results for (H2O)n/Ar systems (n = 1−6). Section 3.3 reports the results of MD simulations and finite-temperature IR spectra. Discussion and perspectives are proposed in section 4.

2. THEORETICAL METHODS 2.1. Electronic Structure: The SCC-DFTB/FF Model. 2.1.1. Description. We briefly summarize the SCC-DFTB/FF model introduced elsewhere.30 The total energy of the active system A consisting of atoms a surrounded by an inert rare gas (Rg) environment consisting of atoms c is calculated as occ

E

SCC‐DFTB/FF

=

A Vrep

+

A Vdisp

+

∑ ni⟨ψi|h0A + W ARg|ψi⟩ i

1 + 2

∑ (γab + γabpol)ΔqaΔqb + V RgRg ab

(1)

where the molecular orbitals (ψi) are expressed as a minimal (valence) linear combination of atomic orbitals ϕaμ (LCAO expansion) ψi =

caiμϕaμ

∑ a∈A ,μ

hA0 is the monoelectronic part of the SCC-DFTB Hamiltonian of the active system A at the reference density, ni the molecular orbital occupation number, and VArep a repulsive contribution expressed as a sum of atomic-pair terms. VAdisp is a dampeddispersion contribution for the active part, also expressed as a classical pair-additive potential36 A Vdisp =−

∑ a < b; a , b ∈ A

fdisp (Rab)

C6ab 6 Rab

The term (1/2)∑abγabΔqaΔqb covers the second-order energy correction expressed as a function of the atomic charge fluctuations (Δqa) on the active system. The atomic charges can refer either to the Mulliken charges, traditionally used in the SCC-DFTB method, or to modified charges following the class IV charge model 3 formalism37−39 (CM3). We have used here the latter charge definition as CM3 charges were determined for the O−H bond to improve the SCC-DFTB description of water clusters.32−34 In addition, CM3 charges have been shown to improve the description of the electric field originating from a molecule,40 which is mandatory for a B

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Σ (1s)2(2s)2σ1π2π2, which in terms of one-electron formulas, can be expressed as 2 +

relevant description of polarization and for obtaining IR spectra. The term (1/2)∑abγpol ab ΔqaΔqb is the polarization contribution of the host matrix atoms to the energy30 γabpol = − ∑ αcfpol (Rac)fpol (R bc) c ∈ Rg

HAr HAr HAr E 2HAr ≡ ucore + udisp + = us Σ

R ac R bc Rac3

+

Σ

E 3OAr Π

with αc the polarizability of the Rg atom c and f pol(Rac) a damping function (see Supporting Information for the analytical expression) for the electric field induced by Δqa on c. The Rg−Rg interactions are included through a classical pairadditive contribution VRgRg (Aziz potential41 for Ar−Ar). The perturbation term WARg is expressed as W A Rg =

∑∑





OAr OAr O Ar E 2O +Ar = ucore + udisp + 4uπOAr + uσOAr + u pol Σ

The RCCSD(T) basis set was the aug-cc-pV6Z set for H and O, and the aug-cc-pV5Z set for Ar with basis set superposition errors (BSSE) taken into account using the counterpoise method. The calculations were performed with the MOLPRO package.42 The analytical expressions of the functions uXAr λ (R) for H−Ar and O−Ar interactions are reported in the Supporting Information, along with the values of the optimized parameters (Table S1 in the Supporting Information). The computed RCCSD(T) data points and the resulting fitting curves for O0/−Ar and H0/+Ar are reported in Figures S1 and S2 in the Supporting Information, respectively. In the case of HAr, the equilibrium distance Re (6.8 bohr) and dissociation energies De (32 cm−1) are in agreement with previous RCCSD(T) results.43,44 In the case of OAr, we found De values slightly lower than those of some of the more extensive RCCSD(T) calculations published in the literature (AV6Z/bf33221 basis set)45 as we obtained 83 cm−1 versus 87.8 cm−1 for the 3Π state, 45 cm−1 versus 48 cm−1 for the 3Σ− state, and 834 cm−1 versus 852 cm−1 for the 2Σ+ state (Table S1 in the Supporting Information). 2.2. Determination of Minima on the Potential Energy Surfaces. All SCC-DFTB/FF calculations were performed with the deMonNano code.46 A two-step procedure was adopted to optimize the geometries of H2OArn clusters. First, global optimization was performed using a parallel tempering Monte Carlo (PTMC) algorithm, recently implemented in the deMonNano code46,47 and previously used for (C6H6)+/0Arn clusters.30 During the PTMC simulations, 107 configurations for each temperature were computed for 16 temperatures ranging from 10 to 100 K following a geometrical series, with structure exchange every 100 steps. The lowest-energy structures over all simulations (∼5 to 10 for each size) were further locally optimized with the gradient method, and the analysis of the harmonic vibrational frequencies by full diagonalization of the Hessian matrix was carried out to ensure the occurrence of real minima. Only the lowest-energy (expectedly global) minima obtained following this procedure are discussed in section 3.1. In the case of water clusters within an Ar matrix ((H2O)n/Ar (n = 1−6)), the Ar matrix is represented by a large cluster consisting of ∼340 Ar atoms organized in a face centered cubic (fcc) crystalline structure. As starting geometries, we substitute p neighboring Ar atoms by n water molecules (p ≤ n). A preliminary geometry optimization with conjugated gradient is performed. In the case of small water clusters (monomer to trimer), the possibility of insertion sites was tested, but after geometry optimization, the water monomers systematically moved significantly, rearranging the rare gas environment around them and layering out one Ar atom per water molecule. For these clusters, the final optimized geometries were found to

ac ac (ucore + udisp + w ac)

where the interaction of the core electrons of a with c, namely ac uac core, and the a−c dispersion potential udisp have been singled out as scalar terms. These dispersion and core terms are isotropic and insensitive to the charge redistribution in the active system. The operators wac are associated with the interactions of the valence electrons of a with Ar atoms. They take into account the anisotropy of the interactions for electrons in atomic orbitals with l ≥ 1 (l orbital quantum number) whose occupation is allowed to vary upon the interaction with Rg atoms. When the laboratory-oriented atomic orbitals ϕaμ are transformed to equivalent orbitals χaλ oriented along the ac ac pair axis via unitary transform matrices Tac λμ, the waμ,aν matrix elements can be expressed as

∑ Tλμac Tλνacuλ(Rac) λ

The one-electron interaction functions to be determined are us for s−Rg interactions, us, uσ, and uπ for sp−Rg interactions. In the following, we describe the parametrization of O−Ar and H−Ar interactions. 2.1.2. Parameterization. The ucore(R), udisp(R), and uλ(R) functions for the interactions between a ∈ A and the Ar atoms are determined so that the SCC-DFTB/FF total energies for each pair in its symmetry frame Em(R ) = ucore(R ) + udisp(R ) +

OAr OAr = ucore + udisp + 3uπOAr + uσOAr

OAr OAr OAr − = u E 3OAr + 2uσOAr core + udisp + 2uπ Σ

a ∈ A c ∈ Rg

waacμ , aν =

+

HAr H Ar E1H+Ar = ucore + u pol

R bc3

∑ nλmuλ(R) λ

fit the respective corresponding energies Em(R) of ab initio restricted CCSD(T) (RCCSD(T)) calculations for this pair. In the case of charged states (δq ≠ 0), a polarization contribution upol(R) is added to the summation (see Supporting Information for the analytical expression). The m subscripts label energies for various symmetry and/or charge states and the nλm are the associated occupation numbers. The uλ(R) functions describe the interactions of the valence orbitals with the maximum l value (1s orbital for H, 2p orbitals for O) with the Ar atoms. The interactions of the 1s and 2s orbitals of O have been included in the ucore(R) repulsive scalar term. The details regarding the procedure for the determination of the u functions parameters are detailed in the paper by Iftner et al.30 The functions ucore(R), udisp(R), and uλ(R) for H−Ar and O− Ar interactions (plus upol(R) for H+Ar and O−Ar) were fitted to reproduce the RCCSD(T) potential energy curves for HAr (2Σ+) and H+Ar (1Σ+) and for various configurations of OAr, namely 3Π (1s)2(2s)2σ1π2π1, 3Σ− (1s)2(2s)2σ2π1π1, and O−Ar C

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 1. Structural Parameters (Angles α and ϕ in Degrees, R in Angstroms (Figure 1)) and Interaction Energy (De, cm−1) for the Two Isomers of H2OAr in the Present Work (SCC-DFTB/FF, DFT, and CCSD(T)/DFT Calculations) and in Other Theoretical Works Min1 α SCC-DFTB/FF CCSD(T)//DFTB DFT CCSD(T)//DFT AW2*a PT4*b CCSD(T)*c a

ϕ

Min2 or S* R

De

α

ϕ

R

De

131

0

3.44

0

5

3.77

117

0

3.59

0

0

3.78

74.3 111 69.0

0 0 0.0

3.636 3.44 3.695

111 172 (98) 124 178 (107) 142.98 142.2 137.7

0 0 0

0 0 0

3.702 3.73 3.72

99 174 (88) 134 179 (91) 116.7 115.2 116.4

Ref 53. bRef 54. cRef 55. S* stands for saddle point. The values in parentheses are the De values including BSSE corrections.

ency upon the initial conditions and should allow a relevant sampling of the energetically favorable portion of the PES. The data points corresponding to the first 50 ps of the simulation were removed for the system to reach equilibration.

be similar for all tested initial situations, but the optimization procedure converged faster starting with initial substitution sites. Therefore, for larger water clusters, all water molecules initially occupy locations close to monosubstitution sites, with the number of Ar atoms removed from the original matrix at least equal to the number of embedded water molecules. The minimum-energy structures of the (H2O)n/Ar (n = 1−6) systems were determined using the following procedure. For each (H2O)n/Ar system, a 9 ns on-the-fly Born−Oppenheimer molecular dynamics (BOMD) simulation in the microcanonical ensemble, starting from an optimized structure, is achieved at a kinetic temperature T ∼ 15 K, which is the average temperature of Ar matrix experiments (with a time step δt = 0.9 fs in order to probe the displacement of the Ar atoms with respect to the water cluster, an accurate description of the O−H motion being less necessary). The lowest-energy structures (∼5 to 10) collected along the simulations were further locally optimized with conjugated gradient (quenching). In each case a single minimal energy structure was retained. This is not a full global optimization, which is not relevant in the present situation because it would yield clusterization of the argon part (in the absence of constraints), nor a strictly local optimization because it results from quenches within a moderately high temperature MD. We underline that all atoms were allowed to relax, including those of the Ar cluster surface. 2.3. MD Simulations and Finite-Temperature IR Spectra. On-the-fly BOMD//SCC-DFTB/FF simulations in the microcanonical ensemble were performed for the (H2O)n/ Ar (n = 1−6) systems, with T values in the [11−14 K] range depending on random initial velocities. The MD simulations were run for 0.2 ns with a time step δt = 0.1 fs. Such a small time step is mandatory to probe the stretching modes of the water clusters.34 Finite-temperature IR spectra are derived from these simulations, computing the Fourier transform of the autocorrelation function of the dipole moment: α(ω) ∝ ω 2

∫0

+∞

dt ⟨μ(0) ·μ(t )⟩e iωt

3. RESULTS AND DISCUSSION 3.1. Benchmark. We benchmarked our approach comparing the SCC-DFTB/FF structures and energetics of small

Figure 1. H2OAr geometrical parameters.

H2OArn (n = 1−3) clusters with DFT and wave function calculations performed in the present work or retrieved from the literature. For clusters with larger sizes (n ≤ 13), the SCCDFTB/FF results are compared with FF results from the literature. Geometry optimizations and frequency calculations obtained by diagonalizing the Hessian matrix were performed with the DFT and MP2 methods. For DFT calculations, we used the M062X48 functional from Truhlar’s group, which includes dispersion. Single-point calculations at the CCSD(T) level for the SCC-DFTB/FF, DFT, and MP2 optimized geometries were also carried out. All calculations were performed in conjunction with the 6-311++G(3d,3p) basis sets for O and H while the valence electrons of Ar are described with a sextuple-ζ basis set, d-orbitals for polarization, and a supplementary diffuse d function (α = 0.05) associated with a [Ne] core pseudopotential.49 All the ab initio calculations use the Gaussian09 suite of programs.50 The results for H2OAr are detailed in Table 1, referring to the geometric parameters of the scheme reported in Figure 1. The SCC-DFTB/FF, DFT, and MP2 optimized geometries of the H2OArn clusters (n = 1−3) are reported in Figure 2. The determination of an accurate PES for the H2OAr system is a theoretical challenge. It has been the object of many theoretical studies51−56 showing the flatness of the PES and the

(2)

where μ(t) is the dipole moment vector at time t and ⟨μ(0)·μ(t)⟩ a statistical average. In practice, this statistical average is obtained using the following procedure: For a given simulation, a first average is computed over 1600 overlapping segments of 26 ps each, defined with different (random) initial conditions. This is repeated 100 times, and the results are averaged to give the IR spectrum derived from one simulation. The final IR spectra shown in the present work result from the final statistical averages of the dipole autocorrelation function over 5−10 MD runs. This procedure minimizes the dependD

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 2. Structures of (H2O)Arn clusters (n = 1−3) optimized at the SCC-DFTB/FF (left-hand column), DFT (middle column) and MP2 (right-hand column) levels of theory (see text for details). For n = 3, two views differing by a 90° rotation are represented. Interatomic distances are reported in angstroms.

Table 2. H2OArn Clustersa

Figure 3. SCC-DFTB/FF optimized structures of (H2O)Arn clusters (n = 4−13). The number (n) below each structure indicates the number of Ar atoms.

d n 1 2 3 4 5 9 11 12 13

isomer

SCC-DFTB

DFT

MP2

lit.

3.567 3.02 2.86

3.59 3.03 2.87

3.608b/3.639c 3.150b/3.039c 2.859b/2.883c 3.823c

min. min. min.

3.48 2.94 2.80 3.56 2.13 1.62 3.90 3.70 4.04

the SCC-DFTB/FF dissociation energy (De) value (111 cm−1) is underestimated by ∼20−25% with respect to the value obtained by Makarewicz (137.7 cm−1).55 It is 13 cm−1 lower than the DFT value (124 cm−1) computed in the present work and in fair agreement with the value derived from molecular beam scattering experiments (116 cm−1).57 Adjusting a parameter such as the cutoff distance of the polarization function to increase the H2O−Ar distance leads to a decrease of De. We finally considered that the couple of values (3.44 Å, 111 cm−1) was a reasonable compromise that could be obtained with our hybrid SCC-DFTB/FF scheme. Geometrical and energetical discrepancies between the SCC-DFTB/FF results and reference CCSD(T) data thus remain, despite a careful parametrization of the a−Rg interactions. This might be due to an underestimation of hybridization in the present model or the inherent use of a minimal basis set with DFTB. Regarding H2OAr, we found a second minimum (Min2 in Table 1) 12 cm−1 above Min1, characterized by α = 0°, i.e., the H atoms pointing toward the Ar atom, in a quasi C2v symmetry (ϕ = 5°). An isomer with a very similar geometry was found at the DFT level, but lying 10 cm−1 below Min1. Single points CCSD(T) calculations were performed using the SCC-DFTB/ FF and DFT geometries (see Table 1). Without taking into account BSSE, the Min1 and Min2 isomers were found to be quasidegenerate at both CCSD(T)//SCC-DFTB/FF and CCSD(T)//DFT levels of theory. With the inclusion of BSSE (counterpoise method), Min2 is found to lie 10 cm−1 (16 cm−1) above Min1 at the CCSD(T)//SCC-DFTB/FF (CCSD(T)//DFT) level of theory. These values are in agreement with the SCC-DFTB/FF results. The similarity of the De values obtained at the CCSD(T)//SCC-DFTB/FF and CCSD(T)// DFT levels of theory for Min1 and Min2 further illustrates the

a

Distance d (in angstroms) between the center of mass of the Ar cluster and that of the H2O molecule obtained with the SCC-DFTB/ FF, DFT (present work), and FF approaches. bRef 59. cRef 58.

strong dependence on the level of theory of the minimum geometry, in particular the relative orientation of the H2O molecule with respect to the Ar atom. Using the SCC-DFTB/ FF PES of H2OAr, we found a minimum-energy structure (see Figure 2 and Min1 in Table 1) with the Ar atom in the plane of the water molecule. This is in line with the DFT and MP2 optimized geometries (see Figure 2) and with previous theoretical results. However, as expected, the geometrical parameters (α and R, cf. Figure 1) are strongly dependent on the method. We found α = 131° and R = 3.44 Å within the SCC-DFTB/FF model, whereas 117° (5°) and 3.59 Å (3.68 Å) were obtained from DFT (MP2) calculations, a reference CCSD(T) calculation55 providing α = 69.0° and R = 3.695 Å. Regarding α, the SCC-DFTB/FF value is close to the DFT value (14° difference), as could be expected, and twice the value from the reference CCSD(T) calculation. The H2O−Ar distance is too short by 0.2 Å with respect to the same reference data. However, with this intermolecular bond length, E

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 3. H2OArn Clustersa Ecohes/n n

isomer

DFTB

CCSD(T)/DFTB

DFT

CCSD(T)/DFT

MP2

CCSD(T)/MP2

1

111

172

124

178

225

252

2 3 4 5 9 11

156 200 227 251 308 334 333 356 351 354 347

176 217

134 159

180 219

229 250

241 256

12 13

min. solvated min. solvated min. solvated

lit. b

131 /143c 130.2d/137.7e 176b/186c 214b/229c 240b/254c 254b/269c 321b/334c 368b/365c 345c 387c 347b/363c 347b/383c

a Cohesion energy per Ar atom Ecohes/n (cm−1) obtained with the SCC-DFTB/FF, DFT, MP2 optimized geometries, and single point CCSD(T) calculations (present work). Results from the literature: bRef 59. cRef 58 (FF calculations). dRef 56 (MP4). eRef 55 (CCSD(T)).

Figure 5. Close-up of the isomers of the water hexamer inside the matrix → inside the matrix without the surrounding Ar atoms, and isolated (right-hand side structures).

quantitatively such tiny energy differences, reproduces the energetic order and the small energy difference between the H2OAr isomers, suggesting a flat PES, as first shown by Chalasinski and co-workers52 and further confirmed by higherlevel wave function calculations.54−56 Further benchmark was performed for H2OAr2,3 clusters. When the SCC-DFTB/FF, DFT, and MP2 geometries of H2OAr2,3 clusters (Figure 2) are compared, the discrepancies in the relative orientation of H2O with respect to the Ar2 and Ar3 clusters are similar to that found for the model H2OAr system. Regarding the distance d between the center of mass of H2O and that of Ar2,3 (Table 2), the SCC-DFTB/FF value remains smaller than the DFT and MP2 values, but the agreement seems to improve as the size of the system increases (−0.08/− 0.11 Å to −0.06/−0.07 Å). The Ar−Ar distance obtained with

Figure 4. Close-up of the water molecule and clusters (H2O)n (n = 2− 4) inside the matrix (left-hand column) and isolated (right-hand column). The intermolecular bond lengths (in angstroms) refer to the smallest intermolecular O−H and O−Ar distances.

flatness of the PES. The isomer with a geometry similar to that of Min2 was found to be a transition state (quoted S*) with correlated wave function approaches (see Table 1).53−55 The energy difference between S* and the most stable isomer Min1 was found to be 21.3 cm−1 at the CCSD(T) level of theory,55 which is larger than the SCC-DFTB/FF value (12 cm−1). Therefore, our model, although not able to describe F

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 6. Optimized structures of (a) (H2O)4/Ar (the fcc elementary cell is also represented); (b) (H2O)6/Ar for cyclic (H2O)6 with a 6-Ar atom vacancy (the displaced Ar atom (upon relaxation, originally at the center of the 6-member ring) leading to the perturbation is indicated with a *); and (c) (H2O)6/Ar for cyclic (H2O)6, with a 7-Ar atom vacancy. Both (H2O)4 and cyclic (H2O)6 are located in the [111] symmetry plane.

Table 4. Energetic Data (cm−1; See Text for Definitions) for (H2O)n/Ar as a Function of Water Cluster Size (n) and Number (p) of Argon Atoms Removed n

p

Esub[(H2O)n]

Eins

Esub[nH2O]

Eint/n

Evac

1 2 3 4 6 6 6 6 6 6 6

1 2 3 4 7 6 7 6 6 7 7

158 371 1054 2261 3768 3291 3116 5160 4932 4634 4400

−551 −1047 −1072 −576 −1195 −963 −1847 +904 +678 −328 −200

158 −1 386 −4 211 −7 068 −12 144 −12 342 −12 184 −10 009 −10 009 −10 538 −10 538

−1357 −1287 −1068 −818 −817 −791 −876 −838 −876 −830 −869

1 515 2 945 4 258 5 533 8 670 8 037 8 372 10 188 10 188 9 614 9 614

prism cage bag chair/cycle boat/cycle chair/cycle boat/cycle

in better agreement with the SCC-DFTB/FF data, although lower by 17−24 cm−1. When comparing with De values obtained by wave function calculations, one should be aware that wave function results (MP2 and CCSD(T)) may be dependent on the size of the basis set. Our results also suggest an important influence of BSSE effects. An exhaustive study of basis set effects should be performed to improve the discussion, but this is beyond the scope of this work. For all cluster sizes, the SCC-DFTB/FF Ecohes/n values are compared with FF data. The former are less than 30 cm−1 (∼8% of the Ecohes/n value) smaller than the latter. The evolution of Ecohes/n as a function of n at the SCC-DFTB/FF level is also consistent with the FF results (see Table 3). The present model offers a simple quantum alternative to the use of classical FFs, some sophisticated versions being the REBO60 or the REAXFF61 reactive potentials built to describe the reactivity of hydrocarbons. In the case of water, inclusion of polarization contributions seems unavoidable. For instance, elaborate FFs providing a good description of isolated water clusters have been proposed recently.62,63 To our knowledge,

our scheme is closer to the MP2 value (−0.08/−0.11 Å; see Figure 2) than to the DFT value (longer by 0.25/0.31 Å). For larger H2OArn clusters (n = 4−13), we compared the SCCDFTB/FF structures to those obtained by force-field (FF) calculations available in the literature.58,59 For all SCC-DFTB/ FF minima for n ≤ 13, the H2O molecule is found at the surface of the cluster (Figure 3). For n ≤ 9, our results are consistent with the FF calculations.58,59 For n = 11−13, a solvated form in which the H2O molecule is embedded within the Ar cluster is also found at low energy. Such structures were actually found as the global minima in FF calculations, for n = 1358 or for n = 12,13.59 With the SCC-DFTB/FF model, they are actually not the global minima for n = 11, 12, or 13 but local minima lying +11, +55, or +91 cm−1 above the global minima, respectively. We now discuss the energetics, and in particular the values of Ecohes/n defined as the normalized values of the evaporation energy associated with the reaction H2OArn → H2O + nAr (Table 3). For H2OAr2,3 clusters, the SCC-DFTB/FF values are larger than the DFT data and lower than the MP2 data. The CCSD(T)//SCC-DFTB/FF and CCSD(T)//DFT values are G

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 7. MD simulation (T = 13 ± 2 K) of H2O in the Ar matrix. (a) Evolution of the angle θ between a chosen O−H bond and a fixed Ar−Ar direction as a function of time (first steps of the MD simulations). Finite-temperature IR spectra in the ν2 (b) and ν1−ν3 (c) regions of isolated H2O (red bands) and H2O/Ar (black bands).

Table 5. Positions (cm−1) of the IR Bands (Spectra in Figure 7b,c) of the Water Molecule, Isolated H2O, and Embedded in the IR Matrix H2O/Ara mode

H2O this work

ν1 ν2 ν3

3805 3804 1565 1562 4079 4078

H2O/Ar

exp. gazb

exp. matc,d

3783 (−22) 3808 1568 (+3) 1571 4050 (−29) 4081

3657.05

3638.5, 3639.2 (−18)

1594.59

1589, 1589.1 (−5)

3755.79

3733.5, 3732.2 (−23)

other work (theor.) 3855,e 3660.28f 1621,e 1594.97f 3982,e 3757.97f

a Shifts are indicated in parentheses. The numbers in italics refer to the harmonic spectra. The experimental values are taken from the literature. bRef 67. cRef 3. dRef 2. eRef 66; MP2/aug-cc-pvtz harmonic calculations. fRef 68; post-CCSD(T) calculations with anharmonic corrections.

H

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 6. Positions (cm−1) of the IR Bands (Spectra in Figure 8a−c) of Isolated (H2O)2 and (H2O)2/Ara mode

H2O

(H2O)2 this work

ν1

3805

ν2

1565

ν3

4079

ν1Ac ν1Dc ν2Dc ν2Ac ν3Ac ν3Dc

3800 3675 1584 1563 4071 3997

(H2O)2/Ar

exp. gazb

exp. Ar matc

other work (theor.) 3838,d 3748,d 1655,d 1624,d 3955,d 3934,d

3671e 3597e 1600e 1590e 3776e 3769e

3783 3650 1589 1562 4051 3975

(−17) (−25) (+5) (−1) (−20) (−22)

3660 3601 1616 1599 3745 3735

3633 3574 1610 1593 3730 3708

(−27) (−27) (−6) (−6) (−15) (−27)

The values of the shifts (cm−1) are reported in parentheses. The experimental values are taken from the literature. bRef 1. cRef 4, as well as ab initio values. dRef 66; MP2/aug-cc-pvtz harmonic calculations. eRef 14; MP2/aug-cc-pVQZ with anharmonic corrections. a

Figure 8. Finite-temperature IR spectra (T = 13 ± 2 K) in the ν2 (a), ν1 (b), and ν3 (c) regions of isolated (H2O)2 (red bands) and (H2O)2/Ar (black bands).

no such elaborate studies describing water clusters interacting with Ar clusters and matrices are presently available; all examples are performed essentially with additive interactions

between the water centers and the rare gas atoms. The present SCC-DFTB/FF model takes into account explicitly the valence electrons of the water clusters, especially those of oxygen, and I

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 7. Positions (cm−1) of the IR Bands of Isolated (H2O)3 and (H2O)3/Ara mode

H2O

(H2O)3

(H2O)3/Ar

exp. gazb

exp. matc

ν1 ν2

3805 1565

3514 (−19) 1601

4079

3607 (−27) 1594 1576 3958 (−13)

3533

ν3

3634 1594 1576 3971

3726

3699 (−27)

layers of the total (H2O)n/Ar (n = 3−6) systems can be found in the Supporting Information (Figure S4). As mentioned in section 2.2, prior to the geometry optimization procedure, the water molecules occupy locations close to substitution sites. In (H2O)2/Ar, the two water molecules are initially centered at the nodes (0 0 0) and (0 (1/ 2) (1/2)) of the fcc lattice, i.e., the two nearest-neighbor Ar atoms are substituted by H2O molecules. This seems reasonable because the Ar nearest-neighbor interatomic distance (dAr−Ar ∼ 3.70−3.71 Å) is of the order of magnitude of the sum of the intermolecular O−O distance (2.85 Å in the water dimer) and intramolecular O−H distance (0.96 Å). For (H2O)3/Ar, the third removed Ar atom is also the nearest neighbor of the first two removed Ar atoms. For (H2O)4/Ar, the fourth removed Ar atom is the nearest neighbor of two of the previously removed Ar atoms to form a four-member ring cycle. For (H2O)6/Ar, the nearest neighbors were removed in such a way that the positions of the initially removed Ar atoms are approximately at the positions of the water molecules within the various isomers. In the following, we discuss the geometries obtained after relaxation of the whole (H2O)n/Ar system. In the case of small water clusters with 2D structures (cyclic forms, n ≤ 4), the embedding in the matrix has little influence on the water cluster structures. As can be seen in Figure 4, the water trimer and tetramer in particular maintain their cyclic structures and the intermolecular water−water distances vary by less than 0.05 Å. Reciprocally, as can be seen in Figure 6a and in Figure S4a,b in the Supporting Information, the fcc lattice structure is hardly perturbed at the vicinity of the water cluster, with the distances between the nearest-neighbor Ar atoms changed by 0.01/0.02 Å. However, the atoms close to the surface of the Ar clusters may be displaced and move slightly out of the initial [111] crystal plane (see Figure S4a,b). In the final optimized geometries of (H2O)n/Ar (n = 1−3) the water molecules remain in positions close to their initial locations (substitution sites, see Figure S4a in the Supporting Information for instance). During the optimization procedure, the intermolecular O−O distance, initially at ∼3.70 Å decreases down to 2.80 Å. In the case of H2O/Ar, the shortest H2O−Ar distances are found for the two Ar atoms in the plane of the water molecule located on the bisector plane of the HOH angle ((O−Ar)Min = 3.47 Å, (H−Ar)Min = 3.44 Å, cf. Figure 4). For (H2O)2/Ar, the shortest O−Ar (H−Ar) distance is found to be 3.39 Å (3.03 Å) with the oxygen and hydrogen atoms pertaining to the proton acceptor (PA) monomer (notation from ref 2). For (H2O)4/Ar, in contrast with (H2O)n/Ar (n = 1−3), reorientation of the water cluster and rearrangement of the neighboring Ar atoms around are observed after structural relaxation. This can be rationalized because the geometry of the cluster is not adapted to that of the matrix. The final structure can be visualized from the representative picture of the total system reported in Figure S4b in the Supporting Information. Three water molecules occupy sites close to the initial monosubstitution sites whereas the last monomer location is similar to a cramped insertion site. This corresponds to the occupation by the water cluster of a 3-Ar atom vacancy rather than the initial 4-atom vacancy. Regarding (H2O)6/Ar, all (H2O)6 isomers (which are all 3dimensional (3D) in the gas phase) have been investigated individually (Figures 5 and 6b,c and Figure S4c−g in the Supporting Information). The names of the isomers (prism, book, cage, bag, chair, boat) are those used in previous works.33,64 No minimum book structure was found in the

The values of the shifts (cm−1) are reported in parentheses. The experimental values are taken from the literature. bRef 69. cRef 8.

a

Table 8. Positions (cm−1) of the IR Bands of Isolated (H2O)4 and (H2O)4/Ara mode

H2O

(H2O)4

ν1 ν2

3805 1565

ν3

4079

3558 1599 1581 3956

a

(H2O)4/Ar 3530 1598 1582 3955

(−28) (−1) (+1) (−1)

exp. gazb

exp. matc

3416

3372 (−44)

3714

3694.8 (−19.2)

The experimental values are taken from the literature. bRef 69. cRef 6.

allows the distinction between π and σ interactions of the hybridized sp orbitals of H2O with the Ar atoms: the anisotropy of the interactions and the hybridization of the orbitals are taken into account. Another advantage of the present quantum/ FF model is that the occupation of the orbitals may vary, allowing for (i) a better account of the charges inducing polarization of the rare gas atoms and (ii) the study of charged water clusters, or water clusters with a OH− or H3O+ impurity, and therefore the determination of ionization potentials or possibly electron affinities. For instance, the influence of the size of the solvating Ar cluster on the ionization potential of benzene in C6H6Arn aggregates was the object of a previous study.30 The quantum treatment of the trapped species may also be more relevant for reactivity studies of water with other molecules. The transferability of the model has been proven from aromatic to aliphatic hydrocarbons.30 In the framework of the present work (description of O−Ar and H−Ar interactions), benchmark calculations similar to those performed for H2OArn clusters (n = 1−3) have been achieved for X−Ar clusters (X = H2, O2, H2O2,CH3OH). The results (minima geometries and dissociation energies De) are reported in Figure S3 and Table S2 of the Supporting Information. As explained at the beginning of this section, the accurate determination of the H2OAr PES is challenging. The results of the model for this system, although not in quantitative agreement with high-level correlated wave function approaches, appear to be in satisfactory agreement with DFT and experimental57 results. However, the benchmark results are quite satisfactory for the dynamic studies in large Ar clusters mimicking the matrix environment. As will be shown in section 3.3, in these simulations, the water molecule rotates in the matrix, in line with experimental data. The comparison with experiments will prove indirectly the validity of our approach, allowing us in particular to mention that the relative orientation of the H2O molecule with respect to the Ar atom is not expected to impact the IR shifts induced by the matrix. 3.2. Structures and Energetics of Water Clusters Embedded in an Ar Matrix. 3.2.1. Structures. The minimum-energy structures for (H2O)n/Ar (n = 1−6) were obtained following the geometry optimization procedure described in section 2.2. Close-ups around the water clusters are reported in Figures 4 and 5a−d. Views of representatives J

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 9. Finite-temperature IR spectra in the ν2 and ν1−ν3 regions of the prism (a and b with T = 11 ± 1 K) and bag (c and d with T = 14 ± 1 K) isomers of isolated (H2O)6 (red bands) and (H2O)6/Ar (black bands).

been removed, also leading to a stable planar structure (Figure S4d in the Supporting Information) with hardly no perturbation of the matrix (Figure 6c). In the case of the cycle isomer, all molecules occupy substitution sites in the relaxed geometries (see Figure S4c,d). The geometry of the system is indeed adapted to that of the matrix, and no significant local rearrangement is observed during structural relaxation. In the case of the prism and bag, matrix distortions are found to occur during geometry relaxation (see Figure S4e,f in the Supporting Information), which can be understood because, as in the case of the tetramer, their geometry is not adapted to that of the matrix and besides, they are 3D structures. In the prism case, five water molecules clearly occupy substitution sites (Figures S4e, parts 2 and 3). The position of the last atom is harder to describe because it is found in a distorted symmetry plane (Figure S4e, part 2). In the

matrix. Stable structures for the prism and bag isomers were found removing initially one extra Ar atom (removal of 7-Ar atoms instead of 6), close to or inside, respectively, the 3D structure of the water cluster. After geometry optimization, the structures of the prism and cage isomers are found to be hardly altered (fluctuations of water intermolecular distances by less than 0.05 Å, angle modifications by less than 5°). In contrast, the bag structure is quite distorted with two quasiorthogonal 4 member rings (with two shared molecules, see Figure 5). The optimizations starting from both chair and boat isomers led to the same planar monocyclic isomer, located in the [111] symmetry plane of the lattice, and designed as “cycle” hereafter. Stable structures have been obtained removing initially 6-Ar atoms (Figure S4c in the Supporting Information), but after geometry optimization, the matrix has been found to be severely perturbed (see Figure 6b). An extra Ar atom has thus K

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 10. Finite-temperature IR spectra in the ν2 and ν1−ν3 regions of the cage (a and b with T = 12 ± 1 K) and chair/cycle (c and d with T = 12 ± 1 K) isomers of isolated (H2O)6 (red bands) and (H2O)6/Ar (black bands).

matrix, the prism appears to substitute to about five Ar atoms. The bag case is similar to the prism case (Figure S4f). Regarding the cage (Figure S4g), the six water molecules occupy locations close to substitution sites (not exactly because of geometry relaxation, see Figure S4g, part 3 for instance) with one final extra Ar atom removed at the vicinity of the cluster (Figure S4g, part 2). The insertion of the cage would then rather correspond to a 7-Ar atom vacancy, and the distortion of the matrix is small. We mention here some limitations of our approach: A systematic search for all possible trapping sites has not been achieved. The deposition of the matrix layers using MD simulations in periodic boundary conditions as in the paper by Kyrychenko and Waluk65 has not been investigated in the present work. The final structures are finite-size local minima,

which has to be kept in mind when discussing the energetics and IR spectra hereafter. 3.2.2. Energetics. Several quantities are defined hereafter to describe the energetics of the water clusters embedded in the Ar matrix. We work here with a finite lattice model to describe the inert matrix, and we wish to compare various vacancies and isomer configurations for a given cluster size. Thus, as a preliminary step, we optimized the finite homogeneous matrix piece alone (without any vacancy and starting from the perfect crystal geometry) in order to provide a reference energy assuming a cancellation of the finite lattice effects (surface relaxation), assumed to be identical with or without insertion of the water molecules or clusters in the calculated energy differences. This relaxed matrix geometry will be considered as the reference, and its energy will be labeled Emat * . The energy of the host matrix (matrix with a vacancy) is labeled Ehost and L

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 9. Positions (cm−1) of the IR Bands (Spectra in Figures 9 and 10) of the Isomers of Isolated (H2O)6 (isol.) and (H2O)6/ Ara prism

a

H2O

isol.

ν1

3805

ν2

1565

3645 3603 3579 3533 3465 3293 1646 1604 1569

ν3

4079

3963 3941 3919 3900 3859 3805

bag /Ar

3632 3562 3448 3363 3328

(+102) (+32) (−82) (−167) (−202)

1645 1615 1577 1564 3976 3920 3901 3847 3806 3763

isol.

cage /Ar

3656 3609 3561 3528 3467 3412 1645 1593

3643 3602 3551 3503 3437 3345 1628 1619 1580

3974 3948 3937 3803

3971 3943 3920 3735 3720

isol.

(+113) (+72) (+21) (−27) (−93) (−185)

3615 3591 3576 3515 3498 3324 1631 1605 1576 1526 3968 3950 3923 3840 3825

chair/cycle /Ar

isol.

/Ar

3542

3514(−16)

1633 1614 1578

1596

1596

3954 3935 3923 3827 3796

3956

3936

3589 3575 3547 3490 3336

(+59) (+45) (+17) (−40) (−194)

For the ν1 bands, the numbers in parentheses indicate the shifts with respect to the ν1 band position of the tetramer (3530 cm−1).

E int = E(H2O)n /Ar − E host − E(H 2O)n

corresponds to the FF contribution to the energy of a water cluster embedded in the matrix. In the following, E(H2O)n is the SCC-DFTB energy of the optimized isolated water cluster, and E(H2O)n/Ar the energy of the total system (water cluster embedded the matrix) in its SCC-DFTB/FF geometry. We first define the substitution energy of a vacancy by the water cluster * − E(H O) Esub[(H 2O)n ] = E(H2O)n /Ar − Emat 2 n

Eint includes the interaction of the water cluster with the matrix (involving matrix polarization) in the final geometry after relaxation. This quantity can be normalized per water molecule, namely Eint/n. It must be noted that in the case of the hexamer trapping, unwanted surface reorganization of isolated surface atoms (one or two) introduces errors in the calculated differences. The values Esub, Eins, Evac, and Eint/n are reported in Table 4. As can be seen in this table, the substitution energies, Esub, for the water clusters are always positive. This is quite easy to understand in the monomer case (Esub = 158 cm−1). Indeed, in the final geometry, the interaction energy of the removed argon atom with its 12 nearest neighbors is 1079 cm−1 and cumulates up to 1515 cm−1 with the whole host system in the final geometry (the same contributions in the reference matrix geometry are 1181 and 1516 cm−1, respectively, the total vacancy energy being almost conserved because of other shell compensation). On the other hand, the water−argon energy Eint occurs to be smaller than Evac for a single monomer substitution. Indeed, the average water−argon bond interaction with the 12 nearest argon atoms is 77 cm−1 (smaller that the SCC-DFTB of the (H2O)ncomplex), and the sum is 979 cm−1. This results in a positive difference of 155 cm−1, very close to the value Eint calculated for the whole host, namely 158 cm−1. This destabilization increases with size. Obviously, the substitution energy for n independent water molecules is much lower because of incorporation of the water−water interactions in the substitution energy. Alternatively, the insertion energy, Eins, is essentially negative, obviously due to the large magnitude of the surface embedding energy of the vacancy atoms. Table 4 also provides the values of Eint/n. They range from −1347 cm−1 for one water molecule to −791 cm−1 for the hexamer (cage isomer). One may notice the decrease of the magnitude of Eint/n with the cluster size. This trend is roughly consistent with the decrease of the Ar coordination number per H2O molecule. The hexamer case offers a competitive situation. First, let us recall that the boat and chair isomers end up into a single same

(3)

This energy includes the creation cost (positive) of the relaxed vacancy (vacancy plus relaxation of the argon host), namely Evac, calculated at the final geometry, the intramolecular water relaxation (positive), and the interaction energy of the water cluster with the host, Eint. It allows to compare various vacancies and isomers for a given cluster size and does not depend on matrix size, provided convergency is reached (the finite matrix is large enough). In addition, we provide the substitution energy relative to n independent water molecules, Esub[nH2O], namely, replacing the last term of the above equation by nE(H2O). We have also estimated the insertion energy, Eins, following the expression given by Gervais et al.24 assuming that the p removed atoms of the vacancy are reembedded at the cluster surface with their bulk cohesive energies (thus keeping their interaction with the whole matrix). * − E(H O) E ins[(H 2O)n ] = E(H2O)n /Ar + pϵbulk − Emat 2 n (4)

It should be mentioned, however, that this definition conserves the number of argon atoms but implies surface change, which is not fully consistent with the present correction made for the surface energy. With the present Ar−Ar potential, the bulk energy per atom is ϵbulk = −709 cm−1. We will also use the creation energy of a vacancy at the relaxed [water cluster + matrix] geometry * Evac = E host − Emat

(6)

(5)

and the interaction energy of the water clusters with the rare gas host, namely Eint M

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 11. Finite-temperature IR spectra in the ν1 regions of (H2O)4/Ar (black lines, T = 13 ± 1 K) and the isomers of (H2O)6 (red lines)/Ar (a) prism (T = 12 ± 1 K), (b) bag (T = 14 ± 1 K), (c) cage (T = 12 ± 1 K), and (d) cycle (T = 12 ± 1 K). The intensities are normalized with respect to the ν1 band of (H2O)4/Ar.

configuration within the matrix, whatever the vacancy (planar cycle structure). As can be seen in Table 4, the lowest insertion energies are obtained for the bag, prism, and cage isomers, while the highest ones correspond to the chair, boat, and cycle isomers. The same trend is observed for the substitution energies. In the case of the hexamer, it is particularly interesting to examine the vacancy creation energies corresponding to various isomer geometries. They were found in the range [8037−10188 cm−1], the smallest value corresponding to the cage isomer, and the largest value to the cycle isomer in the 6atom vacancy. Interestingly, the chair/boat isomer is better accommodated within a 7-atom vacancy. 3.3. MD Simulations and Finite-Temperature IR Spectra of (H2O)n/Ar (n = 1−6). In this section, we describe the dynamics of a water molecule and water clusters embedded

in the matrix and the influence of the embedding on the water clusters’ IR spectra. During the BOMD//SCC-DFTB/FF simulations, the Ar atoms are not frozen but free to vibrate around their equilibrium positions. In the case of H2O/Ar, the water molecule is found to rotate, as illustrated in Figure 7a. This result is consistent with experimental data.2 In the case of (H2O)2/Ar, a rotation−torsion motion along the O−O axis, corresponding to the low-frequency torsion mode (225 cm−1), is observed, also in line with experimental features,4 except that no proton−acceptor switching is observed in our simulations. For larger clusters, apart from small-amplitude vibrational motions around equilibrium geometry, no other motion is observed, because of matrix confinement. This is also in line with the experimental data.13 In the following, we discuss the IR spectra derived from the MD simulations. N

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

matrix size on the IR spectrum of H2O/Ar increasing the matrix by a supplementary Ar atoms’ layer (1330 atoms). The position of the ν2 was found not to be affected. The ν1 and ν3 bands were found to be shifted by +5 cm−1 and +8 cm−1 with respect to the results from the small matrix (343 atoms). Working with this larger matrix thus led to a slight improvement of the comparison with experimental shifts regarding the ν3 band, −21 cm−1 (large matrix) instead of −29 cm−1 (small matrix) with an experimental shift of −23 cm−1. For the ν1 band, a shift of −16 cm−1 (large matrix) instead of −22 cm−1 (small matrix) was obtained (experimental value, −18 cm−1). As explained in section 2.2, the IR spectra are obtained from MD simulations computing the Fourier transform of the autocorrelation function of the dipole moment, i.e., within a classical description of the nuclei. The study of the rotational structure of the vibrational spectrum, probably hidden in the band profiles, would likely require a dedicated quantum approach which is beyond the scope of the present paper. The IR spectra of (H2O)2, isolated and embedded in the Ar matrix, are reported in Figure 8a−c and, the band positions are shown in Table 6. The embedding in the matrix leads to a splitting of the bands. This can be understood as, similarly to the bare water dimer, there are several quasi-isoenergetic local minima on the PES, differentiated by the relative orientation of one monomer with respect to the other. In the matrix these minima have different Ar environments that can be distinguished because there is more confinement than in the case of the monomer, leading to sensibly different IR positions. In order to compare with experimental data, the band positions for the (H2O)2/Ar system reported in Table 6 are the results of an average over those of the maxima of the band structures (the same treatment was applied to the IR spectra of (H2O)3,4/Ar discussed hereafter). As in the case of the water molecule discussed previously, we find with our model that the ν1 and ν3 bands are red-shifted with respect to isolated (H2O)2. This result is in agreement with experimental data, and the order of magnitudes of the shifts are also satisfactory with respect to experimental values (see Table 6). Also in line with experimental trends, the positions of the ν2 bands are hardly modified. The positions of the IR bands of the cyclic (H2O)3,4 clusters and the (H2O)3,4/Ar systems are reported in Tables 7 and 8. The influence of the matrix is similar in both cases: the position of the ν2 mode is not affected, while redshifts are observed for the ν1 and ν3 modes. The calculated values of the shifts for the ν1 bands are in correct agreement with the experimental data, whereas those for the ν3 bands are found to be too small, with hardly no shift (−1 cm−1) in the case of the tetramer. The case of the water hexamer is more complex. The IR spectra of the prism and bag structures in the intramolecular water mode regions are reported in Figure 9a,b and Figure 9c,d, respectively. Those of the cage and cycle isomers are reported in Figure 10a,b and Figure 10c,d, respectively. The positions of the bands are reported in Table 9. As for H2O/Ar, we investigated the effect of the size of the matrix on the IR spectrum of the cycle isomer adding a layer of Ar atoms (1325 Ar atoms instead of 337). The position of the ν2 band is not affected; those of the ν1 and ν3 bands are red-shifted by 7 and 5 cm−1, respectively. From both test calculations performed in the present work, the ν2 band position remains unaffected when increasing the size of the matrix, whereas small variations (lower than ±8 cm−1 according to present calculations) may be expected for the ν1 and ν3 band positions.

The accurate determination of the positions of the isolated water clusters’ IR bands with correlated wave function approaches has been the object of many studies (see for instance refs 14−16 and 66, and references therein). The SCCDFTB absolute values of the positions of the IR bands for the ν2 bending mode differ from the best ab initio values by less than 30 cm−1 for (H2O)1,2 (see Tables 5 and 6). Regarding the symmetric and asymmetric stretching modes ν1 and ν3, the SCC-DFTB values are larger by about 100 cm−1 for (H2O)2 and 200 cm−1 for (H2O)3−6 with respect to recent CCSD(T) calculations taking into account anharmonic effects.15 These differences are due to the approximate SCC-DFTB approach (truncation of the DFT Hamiltonian, minimal valence basis set, parametrization). The set of chosen charges for the selfconsistency also affects the potential energy surface and thus the positions of the IR bands. Besides, nuclear quantum effects, likely to affect the band positions (and intensities), are not included in our scheme. Therefore, and as mentioned in a previous paper,34 the absolute band positions are not reproduced with our SCC-DFTB approach. However, this scheme has been shown to reproduce correctly the trends regarding the influence of the adsorption on a surface (benzene) on the positions of the IR bands of a water molecule with respect to experimental data.34 Therefore, in terms of relative positions, SCC-DFTB calculations are likely to provide reliable results; however, a benchmark is necessary. Regarding the relatives intensities of the IR bands for a given spectrum, they are dependent on the set of charges chosen for the self-consistency. Despite the CM3 charges (see section 2.1.1) were determined to reproduce the intermolecular binding energy of the water dimer (H2O)2,34 they do not allow the reproduction of the relative intensities of the water monomer.5 Besides, despite our efforts to reach ergodicity, the relative intensities may still be influenced by the initial conditions. Therefore, in the following, the relative intensities of the bands of a given spectrum are not discussed. In all the reported spectra, they are normalized with respect to the most intense band. In the following, we discuss the influence of the matrix on the positions of the IR bands of the water molecule and of water clusters comparing the results of our SCC-DFTB/ FF scheme with experimental data. For all water clusters embedded in the matrix, the positions of the IR bands are shifted with respect to those of the harmonic spectra, revealing significant low-temperature anharmonic effects. Interestingly, such effects are not observed for the H2O molecule interacting with small Ar clusters (see Table S3 of the Supporting Information). The finite-temperature IR spectra of H2O both isolated and embedded in the Ar matrix are reported in Figure 7b,c. For the three water modes (symmetric and asymmetric stretching modes ν1 and ν3, respectively, and bending mode ν2), the matrix induces a broadening of the bands, certainly caused by rotation. The redshifts of the ν1 and ν3 bands are likely due to the flattening of the PES induced by the interaction with the Ar atoms of the matrix. These shifts were not observed in the case of the interaction with small Ar clusters. The values of the band positions are reported in Table 5, along with the experimental data. The values of the shifts induced by the interaction with the matrix are reported in parentheses. The match between the shifts computed in the present work and the experimental data is satisfactory, showing the validity of the approach. As explained in section 2.2, all simulations are performed with an Ar cluster of ∼340 atoms. We tested the influence of the O

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

sufficient to reproduce accurately the relative intensities of the bands of a given spectrum, it reproduces the splitting of the oscillator strength into the various mode components. On the basis of energetical criteria (see section 3.2.2), the concentration of the cycle should be lower than that of the quasidegenerate cage, bag, and prism isomers and the intensity of its ν1 band in the experiment should then be smaller than that reported in Figure 11d. However, the population ratio of the various isomers in the matrix also depends on their formation mechanism during deposition and on their diffusion in the matrix. Therefore, using only energetical data does not allow us to conclude definitively on the experimental relative abundances and IR band intensities. The computed spectrum of the cage isomer in the ν1 region is compatible with the two bands shifted by −47 and −161 cm−1 with respect to the ν1 band of (H2O)4 in the experiment: the cage spectrum presents two bands shifted by −40 cm−1 and −194 cm−1 with respect to the computed (H2O)4 ν1 band (see Figure 11c). Therefore, the presence of the cage could contribute to the two bands assigned to the cycle and cage isomer, respectively,6 all the more our results indicate it should be the most stable isomer in the matrix. The computed higher-energy band in the cage spectrum might contribute to the blue tail of the (H2O)4 band observed experimentally. The presence of the cycle isomer is also possible even though the shift with respect to the (H2O)4 ν1 band is found to be smaller than the experimental shift (−16 versus −47 cm−1). Despite its stability, the presence of the bag may be excluded because of the band at 3437 cm−1 (−93 cm−1 with respect to the ν1 band of (H2O)4) that does not seem to match the experimental spectrum. Regarding the theoretical spectrum of the prism, which was also found to be stable in the matrix, many bands of low intensity are found in the ν1 region (see Figure 11a). Its presence cannot be ruled out as its lowintensity bands might be hidden in the bands of the other water clusters in the experimental spectrum. In this comparison with the experimental spectrum,6 the results from our model seem to be in line with the formation of the cage and possibly the cycle isomers. These conclusions should be taken with caution as our model contains many approximations, as mentioned above.

In the following, we discuss the influence of the matrix on the ν1 and ν3 bands, starting from the prism to the cycle. In the case of the cycle, only the most stable isomer (7-Ar atom vacancy) is discussed, the IR spectrum of the isomer with 6-Ar atom vacancy can be found in Figure S5 of the Supporting Information. The ν2 mode has not been investigated experimentally to our knowledge. We may only mention that the influence of the matrix is minor on this mode (Figures 9a,c and 10a,c and Table 9). Regarding the prism structure, the six bands of the ν1 modes in the free cluster are found in the [3293−3645 cm−1] range. In the matrix environment, eight bands are found in a narrower range of energy, [3328−3632 cm−1]: the lowest-energy band is blueshifted, whereas most bands are red-shifted (see Figure 9b). The ν3 modes are found in the [3805−3963] cm−1 range in the free cluster and in the [3763−3976 cm−1] range in the matrix environment; the lowest-energy band appears as the most red-shifted. Regarding the bag isomer, the six bands of the ν1 mode are found in the [3412−3656 cm−1] range for the free cluster and in the [3345− 3643 cm−1] range in the matrix (see Figure 9c); all the bands are red-shifted, and the larger shift concerns the lowest-energy band (−67 cm−1). The ν3 modes are similarly red-shifted ([3803−3974 cm−1] for the free cluster and [3720−3971 cm−1] for the embedded cluster), with the larger shift for the lower-energy band. In the case of the cage structure (Figure 10b), the spectrum is hardly modified; all bands are red-shifted except the lowest-energy ν1 band, similarly to the prism spectrum. The ν1 bands are found in the [3324−3615 cm−1] range for the isolated cluster and in the [3336−3589 cm−1] range for the embedded cluster. The ν3 bands are found in the [3825−3968 cm−1] range for the isolated cluster and in the [3796−3954 cm−1] range for the embedded cluster. As in the case of the prism structure, the lowest-energy band is the most red-shifted. Finally, regarding the chair/cycle isomer (Figure 10c), the ν1 and ν3 bands of the embedded cluster are broadened and red-shifted with respect to those of the isolated chair isomer. The ν1 band is found at 3514 cm−1, and this corresponds to a redshift of 28 cm−1 with respect to the ν1 band of the isolated chair isomer (3542 cm−1). The ν3 band is redshifted by 20 cm−1 after embedding in the matrix (3936 cm−1 versus 3956 cm−1). In IR spectra of water clusters obtained from Ar matrix experiments,6 the band at 3372 cm−1 was assigned to (H2O)4 and those at 3325 cm−1 and 3211 cm−1, of lower intensity, were assigned to the cycle (c-(H2O)6) and to the cage isomers of (H2O)6, respectively. Thus, in the IR spectrum, two bands of lower intensity are found shifted by −47 and −161 cm−1 with respect to the ν1 band of (H2O)4. In order to compare our results with these experimental data, the IR spectra of the prism, bag, cage, and cycle isomers along with that of (H2O)4, in the ν1 region, are reported in Figure 11. Intensities are normalized with respect to the ν1 band of (H2O)4. The computed shifts of the bands of the (H2O)6 isomers with respect to that of (H2O)4 are reported in parentheses in Table 9. Assuming the same abundance for all clusters, the relative intensities in the ν1 region are low for the isomers of (H2O)6 with respect to (H2O)4 (see Figure 11a−c), except for the cycle isomer (Figure 11d), which is in agreement with the experimental spectrum.6 This was expected because for the prism, bag, and cage isomers, the oscillator strength is split into many components, which is not the case for cyclic (H2O)4 and c-(H2O)6, whose ν1 mode intensities are concentrated into one or two bands (Figure 11d). Although our approach may not be

4. CONCLUSION We applied the SCC-DFTB/FF model developed to describe the interactions of molecular systems with rare gas atoms30 to the determination of the structural, energetical, dynamical, and spectral properties of water clusters embedded in an Ar matrix. This modeling first required the parametrization of O−Ar and H−Ar interactions and benchmark studies on H2OArn clusters. Using this benchmarked SCC-DFTB/FF potential, the structures and energetics of small water clusters (H2O)n (n = 1−6) embedded in an Ar matrix represented by a fcc crystalline cluster consisting of ∼340 Ar atoms were determined. The structures of 2D clusters (n = 1−4) are hardly modified by the embedding in the matrix. In the case of 3D structures (n = 6), noticeable distortions may occur: the chair and boat become planar while the two 4-member rings of the bag structure become planar and orthogonal (one ring in a reticular plane of the fcc lattice). Insertion and substitution energies of water clusters in the matrix were determined for all cluster sizes. Regarding the hexamer, insertion and substitution energies for the prism, cage, and bag were found to be lower than that for the 6-member ring isomer. Using this SCC-DFTB/FF potential, BOMD simulations were performed. The behavior P

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

through the grant NEXT ANR-10-LABX-0037 in the framework of the “Programme des Investissements d’Avenir”.

of the studied systems at a kinetic temperature within 11−14 K is in line with existing experimental data: the rotation of the water monomer is observed, as well as the rotation−torsion structure of the dimer, whereas no global motion of the larger clusters is found to occur. Finite-temperature IR spectra are derived from these simulations. The IR spectra of (H2O)n/Ar systems present broader and more structured bands than those of bare (H2O)n clusters. Some new bands also appear because of the different environment seen by the H atoms of the water molecules during their internal rotation motions. For all clusters, the position of the ν2 band is hardly modified, whereas the ν1 and ν3 bands are red-shifted. In the case of the water monomer and dimers, the shift values are in fair agreement with the experimental data. The case of the water hexamer, which has received a great deal of interest lately, is the most complex as many isomers are likely to coexist. The comparison of the computed spectra with the experimental spectra by Ceponkus and co-workers6 tends to confirm the presence of the cage (and probably the cycle) isomer, but given the approximations within our scheme, no unambiguous assignment can be achieved. We may indeed underline that the present theoretical work suffers from the following weaknesses: (i) The exploration of the (H2O)n/Ar possible structures is certainly not exhaustive. (ii) The matrix is represented by a cluster, and periodic conditions, although probably not too important for IR spectra, may modify internal shell constraints. (iii) In the SCC-DFTB/FF scheme, the sp hybridization of the O atoms under the interaction with the Ar is not fully taken into account. Finally, (iv) nuclear quantum effects are likely to play a role in the structures, energetics, and finite-temperature IR spectra of molecular systems containing such light atoms as hydrogen. Despite these limitations, we hope that the present work will contribute to the interpretation of the influence of a rare gas matrix on the IR spectra of molecular clusters.





(1) Buck, U.; Huisken, F. Infrared Spectroscopy of Size-Selected Water and Methanol Clusters. Chem. Rev. (Washington, DC, U.S.) 2000, 100, 3863−3890. (2) Perchard, J. P. Anharmonicity and Hydrogen Bonding. III. Analysis of the Near Infrared Spectrum of Water Trapped in Argon Matrix. Chem. Phys. 2001, 273, 217−233. (3) Michaut, X.; Vasserot, A.-M.; Abouaf-Marguin, L. Temperature and Time Effects on the Rovibrational Structure of Fundamentals of H2O Trapped in Solid Argon: Hindered Rotation and RTC Satellite. Vib. Spectrosc. 2004, 34, 83−93. (4) Ceponkus, J.; Uvdal, P.; Nelander, B. Acceptor Switching and Axial Rotation of the Water Dimer in Matrices, Observed by Infrared Spectroscopy. J. Chem. Phys. 2010, 133, 074301. (5) Bouteiller, Y.; Tremblay, B.; Perchard, J. The Vibrational Spectrum of the Water Dimer: Comparison between Anharmonic ab Initio Calculations and Neon Matrix Infrared Data between 14,000 and 90 cm−1. Chem. Phys. 2011, 386, 29−40. (6) Ceponkus, J.; Uvdal, P.; Nelander, B. Water Tetramer, Pentamer, and Hexamer in Inert Matrices. J. Phys. Chem. A 2012, 116, 4842− 4850. (7) Ceponkus, J.; Engdahl, A.; Uvdal, P.; Nelander, B. Structure and Dynamics of Small Water Clusters, Trapped in Inert Matrices. Chem. Phys. Lett. 2013, 581, 1−9. (8) Ceponkus, J.; Karlstro, G.; Nelander, B. Intermolecular Vibrations of the Water Trimer, a Matrix Isolation Study. J. Chem. Phys. A 2005, 109, 7859−7864. (9) Hirabayashi, S.; Yamada, K. M. Infrared Spectra and Structure of Water Clusters Trapped in Argon and Krypton Matrices. J. Mol. Struct. 2006, 795, 78−83. (10) Hirabayashi, S.; Yamada, K. M. The Monocyclic Water Hexamer Detected in Neon Matrices by Infrared Spectroscopy. Chem. Phys. Lett. 2007, 435, 74−78. (11) Leforestier, C.; van Harrevelt, R.; van der Avoird, A. VibrationRotation-Tunneling Levels of the Water Dimer from an ab Initio Potential Surface with Flexible Monomers. J. Phys. Chem. A 2009, 113, 12285−12294. (12) Babin, V.; Leforestier, C.; Paesani, F. Development of a “First Principles” Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient. J. Chem. Theory Comput. 2013, 9, 5395−5403. (13) Tremblay, B.; Madebène, B.; Alikhani, M.; Perchard, J. P. The Vibrational Spectrum of the Water Trimer: Comparison between Anharmonic Ab Initio Calculations and Neon Matrix infrared data between 11,000 and 90 cm−1. Chem. Phys. 2010, 378, 27−36. (14) Dunn, M. E.; Evans, T. M.; Kirschner, K. N.; Shields, G. C. Prediction of Accurate Anharmonic Experimental Vibrational Frequencies for Water Clusters, (H2O)n, n = 2−5. J. Phys. Chem. A 2006, 110, 303−309. (15) Miliordos, E.; Apra, E.; Xantheas, S. S. Optimal Geometries and Harmonic Vibrational Frequencies of the Global Minima of Water Clusters (H2O)n, n = 2−6, and several hexamer local minima at the CCSD(T) level of theory. J. Chem. Phys. 2013, 139, 114302. (16) Temelso, B.; Shields, G. C. The Role of Anharmonicity in Hydrogen-Bonded Systems: The Case of Water Clusters. J. Chem. Theory Comput. 2011, 7, 2804−2817. (17) Trakhtenberg, L. I.; Fokeyev, A. A.; Zyubin, A. S.; Mebel, A. M.; Lin, S. H. Matrix Reorganization with Intramolecular Tunneling of H Atom: Formic Acid in Ar Matrix. J. Chem. Phys. 2009, 130, 144502. (18) Galindez, J.; Calvo, F.; Paska, P.; Hrivnak, D.; Kalus, R.; Gadea, F. DIM Modelings for Realistic Simulations of Ionic Rare-Gas Clusters, Test on Structures and Photoabsorption Spectra of Arn+ (n = 3−8). Comput. Phys. Commun. 2002, 145, 126−140. (19) Kuntz, P.; Valldorf, J. A DIM Model for Homogeneous NobleGas Ionic Clusters. Z. Phys. D: At., Mol. Clusters 1988, 8, 195−208.

ASSOCIATED CONTENT

S Supporting Information *

Additional information regarding (i) the parametrization (section 2.1.2); (ii) the benchmark (transferability of the SCC-DFTB/FF model, section 3.1); (iii) the SCC-DFTB/FF optimized structures of (H2O)n/Ar (section 3.2); (iv) the differences in the harmonic and finite-temperature IR spectra of (H2O)Arn clusters (n = 1−13); and (v) finite-temperature IR spectra. This material is available free of charge via the Internet at http://pubs.acs.org.



REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +33561556033. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors acknowledge Mathias Rapacioli for technical help, the computing facility CALMIP for allocation of computer resources, and CNRS and Region Midi-Pyrénées for PhD grant to C.I. This study has been supported by ANR 13-BS08-0005 PARCS “Polycyclic Aromatic hydrocarbons Reactivity in Cryogenic Solids”, GDR Groupement de Recherche 3533 EMIE Edifices Moléculaires Isolés et Environnés, and partially supported Q

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (20) Gaveau, M.; Briant, M.; Fournier, P.; Mestdagh, J.; Visticot, J.; Calvo, F.; Baudrand, S.; Spiegelman, F. Spectroscopy of Calcium Deposited on Large Argon Clusters. Eur. Phys. J. D 2002, 21, 153− 161. (21) Jungwirth, P.; Gerber, R. Quantum Dynamics Simulations of Nonadiabatic Processes in Many-Atom Systems: Photoexcited BaAr10 and BaAr20 Clusters. J. Chem. Phys. 1996, 104, 5803−5814. (22) Krylov, A.; Gerber, R.; Gaveau, M.; Mestdagh, J.; Schilling, B.; Visticot, J. Spectroscopy, Polarization and Nonadiabatic Dynamics of Electronically Excited BaArn Clusters: Theory and Experiment. J. Chem. Phys. 1996, 104, 3651−3663. (23) Plata, J. J.; Heitz, M.-C.; Spiegelman, F. Effect of Structure and Size on the Excited States Dynamics of CaArn Clusters. Eur. Phys. J. D 2013, 67, 17. (24) Gervais, B.; Giglio, E.; Jacquet, E.; Ipatov, A.; Reinhard, P.; Suraud, E. Simple DFT Model of Clusters Embedded in Rare Gas Matrix: Trapping Sites and Spectroscopic Properties of Na Embedded in Ar. J. Chem. Phys. 2004, 121, 8466−8480. (25) Gervais, B.; Giglio, E.; Jacquet, E.; Ipatov, A.; Reinhard, P.; Fehrer, F.; Suraud, E. Spectroscopic Properties of Na Clusters Embedded in a Rare-Gas Matrix. Phys. Rev. A: At., Mol., Opt. Phys. 2005, 71, 015201. (26) Douady, J.; Gervais, B.; Giglio, E.; Ipatov, A.; Jacquet, E. Spectroscopic Properties of Na3+ in Bulk and on the Surface of Ar Droplets. J. Mol. Struct. 2006, 786, 118−122. (27) Gross, M.; Spiegelmann, F. Theoretical Study of Spectroscopical Properties of Na and Na2 in Argon Clusters and Matrices. J. Chem. Phys. 1998, 108, 4148−4158. (28) Jacquet, E.; Zanuttini, D.; Douady, J.; Giglio, E.; Gervais, B. Spectroscopic Properties of Alkali Atoms Embedded in Ar Matrix. J. Chem. Phys. 2011, 135, 174503. (29) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Material Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260−7268. (30) Iftner, C.; Simon, A.; Korchagina, K.; Rapacioli, M.; Spiegelman, F. A Density Functional Tight Binding/Force Field Approach to the Interaction of Molecules with Rare Gas Clusters: Application to (C6H6)+/0Arn Clusters. J. Chem. Phys. 2014, 140, 034301. (31) Rapacioli, M.; Simon, A.; Dontot, L.; Spiegelman, F. Extensions of DFTB to Investigate Molecular Complexes and Clusters. Phys. Status Solidi B 2012, 249, 245−258. (32) Simon, A.; Spiegelman, F. Conformational Dynamics and FiniteTemperature Infrared Spectra of the Water Octamer Adsorbed on Coronene. Comput. Theor. Chem. 2013, 1021, 54−61. (33) Simon, A.; Spiegelman, F. Water Clusters Adsorbed on Polycyclic Aromatic Hydrocarbons: Energetics and Conformational Dynamics. J. Chem. Phys. 2013, 138, 194309. (34) Simon, A.; Rapacioli, M.; Mascetti, J.; Spiegelman, F. Vibrational Spectroscopy and Molecular Dynamics of Water Monomers and Dimers adsorbed on Polycyclic Aromatic Hydrocarbons. Phys. Chem. Chem. Phys. 2012, 14, 6771−6786. (35) Perez, C.; Muckle, M. T.; Zaleski, D. P.; Seifert, N. A.; Temelso, B.; Shields, G. C.; Kisiel, Z.; Pate, B. H. Structures of Cage, Prism, and Book Isomers of Water Hexamer from Broadband Rotational Spectroscopy. Science 2012, 336, 897−901. (36) Elstner, M. The SCC-DFTB Method and its Application to Biological Systems. Theor. Chem. Acc. 2006, 116, 316−325. (37) Li, J.; Zhu, T.; Cramer, C.; Truhlar, D. New Class IV Charge Model for Extracting Accurate Partial Charges from Wave Functions. J. Phys. Chem. A 1998, 102, 1820−1831. (38) Thompson, J. D.; Cramer, C. J.; Truhlar, D. G. Parameterization of Charge Model 3 for AM1, PM3, BLYP, and B3LYP. J. Comput. Chem. 2003, 24, 1291−1304. (39) Winget, P.; Thompson, J. A.; Xidos, J. A.; Cramer, C. J.; Truhlar, D. G. Charge Model 3: A Class IV Charge Model based on Hybrid Density Functional Theory with Variable Exchange. J. Phys. Chem. A 2002, 106, 10707−10717.

(40) Rapacioli, M.; Spiegelman, F.; Talbi, D.; Mineva, T.; Goursot, A.; Heine, T.; Seifert, G. Correction for Dispersion and Coulombic Interactions in Molecular Clusters with Density Functional derived Methods: Application to Polycyclic Aromatic Hydrocarbon Clusters. J. Chem. Phys. 2009, 130, 244304. (41) Aziz, R. A. A Highly Accurate Interatomic Potential for Argon. J. Chem. Phys. 1993, 99, 4518−4525. (42) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G. et al. MOLPRO, version 2012.1, a Package of Ab Initio Programs, http://www.molpro.net. 2012. (43) Partridge, H.; Schwenke, D. W.; Bauschlicher, C. W. Theoretical Study of the Ground States of the Rare-Gas Hydrides, HeH, NeH, and ArH. J. Chem. Phys. 1993, 99, 9776−9782. (44) Partridge, H.; Bauschlicher, C. W. The Dissociation Energies of He2, HeH, and ArH: A Bond Function Study. Mol. Phys. 1999, 96, 705−710. (45) Garand, E.; Buchachenko, A. A.; Yacovitch, T. I.; Szcześniak, M. M.; Chałasiński, G.; Neumark, D. M. Study of ArO and ArO via Slow Photoelectron Velocity-Map Imaging Spectroscopy and Ab Initio Calculations. J. Phys. Chem. A 2009, 113, 4631−4638. (46) Heine, T.; Rapacioli, M.; Patchkovskii, S.; Frenzel, J.; Koster, A.; Calaminici, P.; Duarte, H. A.; Escalante, S.; Flores-Moreno, R.; Goursot, A. et al. deMon-Nano Experiment http://physics.jacobsuniversity.de/theine/research/deMon/. 2009. (47) Dontot, L. Theoretical Investigation of Structural and Spectral Properties of Cationic PAH Clusters. Ph.D. thesis, University of Toulouse III - Paul Sabatier, Toulouse, France, 2014. (48) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157−167. (49) Castex, M. C.; Morlais, M.; Spiegelman, F.; Malrieu, J. P. Comparison between Experimentally and Theoretically Determined Potential Curves of the Ar2* Lowest States. J. Chem. Phys. 1981, 75, 5006−5018. (50) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian09, revision A.1. 2009; Gaussian Inc.: Wallingford CT, 2009. (51) Bulski, M.; Wormer, P. E. S.; van der Avoird, A. Ab initio potential energy surfaces of Ar−H2O and Ar−D2O. J. Chem. Phys. 1991, 94, 8096. (52) Chalasińsḱ i, G.; Szcześniak, M. M.; Scheiner, S. Ab Initio Study of the Intermolecular Potential of Ar−H2O. J. Chem. Phys. 1991, 94, 2807. (53) Cohen, R. C.; Saykally, R. J. Determination of an Improved Intermolecular Global Potential Energy Surface for Ar-H20 from Vibration-Rotation-Tunneling Spectroscopy. J. Chem. Phys. 1993, 8, 6007−6030. (54) Hodges, M. P.; Wheatley, R. J.; Harvey, A. H. Intermolecular Potentials and Second Virial Coefficients of the Water−Neon and Water−Argon Complexes. J. Chem. Phys. 2002, 117, 7169. (55) Makarewicz, J. Ab Initio Intermolecular Potential Energy Surfaces of the Water-Rare Gas Atom Complexes. J. Chem. Phys. 2008, 129, 184310. (56) Tao, F.-M.; Klemperer, W. Accurate Ab Initio Potential Energy Surfaces of Ar−HF, Ar−H2O, and Ar−NH3. J. Chem. Phys. 1994, 101, 1129. (57) Aquilanti, V.; Cornicchi, E.; Moix Teixidor, M.; Saendig, N.; Pirani, F.; Cappelletti, D. Glory-Scattering Measurement of Water− Noble-Gas Interactions: The Birth of the Hydrogen Bond. Angew. Chem., Int. Ed. 2005, 44, 2356−2360. (58) Liu, S.; Bacic, Z.; Moskowitz, J. W.; Schmidt, K. E. ArnH2O (n = 1−14) van der Waals Clusters: Size Evolution of Equilibrium Structures. J. Chem. Phys. 1994, 101, 8310. (59) Borges, E.; Ferreira, G. G.; Braga, J. P. Structures and Energies of ArnH2O (n = 1−26) Clusters Using a Nonrigid Potential Surface: A Molecular Dynamics Simulation. Int. J. Quantum Chem. 2008, 108, 2523−2529. R

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (60) Brenner, D. W. Empirical Potential for Hydrocarbons for use in Simulating the Chemical Vapor Deposition of Diamond Films. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 42, 9458−9471. (61) van Duin, A.; Dasgupta, S.; Lorant, F.; Goddard, W. A., III. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 15, 9396−9409. (62) Babin, V.; Medders, G. R.; Paesani, F. Toward a Universal Water Model: First Principles Simulations from the Dimer to the Liquid Phase. J. Phys. Chem. Lett. 2012, 3, 3765−3769. (63) Real, F.; Vallet, V.; Flament, J.-P.; Masella, M. Revisiting a Many-Body Model for Water Based on a Single Polarizable Site: From Gas Phase Clusters to Liquid and Air/Liquid Water Systems. J. Chem. Phys. 2013, 139, 114502. (64) Bates, D. M.; Tschumper, G. S. CCSD(T) Complete Basis Set Limit Relative Energies for Low-Lying Water Hexamer Structures. J. Phys. Chem. A 2009, 113, 3555−3559. (65) Kyrychenko, A.; Waluk, J. Molecular Dynamics Simulations of Matrix Deposition. I. Site Structure Analysis for Porphyrin in Argon and Xenon. J. Chem. Phys. 2003, 119, 7318−7327. (66) Xantheas, S. S.; Dunning, T. H. Ab Initio Studies of Cyclic Water Clusters (H2O)n, n = 1−6 I. Optimal Structures and Vibrational Spectra. J. Chem. Phys. 1993, 99, 8774−8792. (67) Benedict, W. S.; Gailar, N.; Plyler, E. K. RotationVibration Spectra of Deuterated Water Vapor. J. Chem. Phys. 1956, 24, 1139− 1165. (68) Kahn, K.; Kirtman, B.; Noga, J.; Ten-no, S. Anharmonic Vibrational Analysis of Water with Traditional and Explicitly Correlated Coupled Cluster Methods. J. Chem. Phys. 2010, 133, 1−12. (69) Huisken, F.; Kaloudis, M.; Kulcke, A. Infrared Spectroscopy of Small Size - Selected Water Clusters. J. Chem. Phys. 1996, 104, 17−25.

S

DOI: 10.1021/jp508533k J. Phys. Chem. A XXXX, XXX, XXX−XXX