- Email: [email protected]
Conformational dynamics and finite-temperature infrared spectra of the water octamer adsorbed on coronene
Computational and Theoretical Chemistry 1021 (2013) 54–61
Contents lists available at SciVerse ScienceDirect
Computational and Theoretical Chemistry journal homepage: www.elsevier.com/locate/comptc
Conformational dynamics and finite-temperature infrared spectra of the water octamer adsorbed on coronene Aude Simon ⇑, Fernand Spiegelman Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, Université de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France
a r t i c l e
i n f o
Article history: Received 17 May 2013 Received in revised form 18 June 2013 Accepted 18 June 2013 Available online 29 June 2013 Keywords: Water clusters PAH Molecular dynamics DFTB Infrared
a b s t r a c t We present in this work a detailed study of all-atoms conformational dynamics and finite temperature infrared (IR) spectra of the water octamer (H2O)8 adsorbed on coronene (C24H12), a compact Polycyclic Aromatic Hydrocarbon (PAH). The potential energy surface is calculated on-the-fly within the self-consistent-charge density-functional based tight-binding (SCC–DFTB) approach, with modifications insuring the correct description of water–water and water–PAH interactions. On-the-fly Born–Oppenheimer molecular dynamics simulations are performed in the temperature T range 10–300 K starting from the most stable cubic conformer of (H2O)8 adsorbed on C24H12. As T increases, diffusion of (H2O)8 starts, followed by isomerisation processes and ‘‘melting’’ of (H2O)8 for T in the range 150–200 K. These dynamical processes are quantified through the temperature evolution of chosen indicators. We show that the adsorption of (H2O)8 on C24H12 leads to lower isomerisation temperatures than for the bare cluster, in line with an increase of the density of isomers. A comparative analysis of the low-T (harmonic) spectra of bare and adsorbed (H2O)8 is provided. Finally, the evolution of the spectra as a function of T, illustrating the conformational dynamics, is discussed. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Water clusters have raised a growing and fundamental interest [1] in several aspects (i) as finite size systems with specific properties [2–6], (ii) as precursors to understand the bottom-up building properties of liquid and ice water [7], and (iii) to understand nanoscale solvation [8], a feature essential in biological systems. The interaction of water molecules and clusters impurities is also of wide interest in atmospheric chemistry, for instance in nucleative processes initiating cloud formation, and various investigations have considered the interaction between water droplets or clusters and impurities such as an extra proton [9,10], or carbon-containing molecular impurities [11]. Among the latter, Polyaromatic Aromatic Hydrocarbon (PAH) molecules show significant abundance in the atmosphere, due to their efficient formation as by-products of natural processes (biomass burning) or of human activities (combustion of fossil fuels) [12]. PAH molecules and derived-PAH species have also raised a strong and lasting interest since they were proposed in the mid-eighties [13,14] as the carriers of the aromatic infrared (IR) bands, that are a set of mid-IR emission bands observed in many regions of the interstellar medium. Despite many experimental and theoretical studies, a specific molecule has not been identified yet [15]. These interstellar PAH ⇑ Corresponding author. Tel.: +33 561556033. E-mail address: [email protected] (A. Simon). 2210-271X/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.comptc.2013.06.023
molecules would form at the border of molecular clouds, in Photodissociation Regions (PDRs). The importance of PAH species, under various forms, in interstellar carbonaceous (very small) grains and their role in the interstellar chemistry is still in debate. Inside molecular clouds, temperatures are expected to be in the range 20–100 K and in PDRs, temperatures can reach between 100 K and 1000 K. Regarding atmospherical conditions, the temperatures in the troposphere, that contains an important quantity of water, are expected to be in the range 200–300 K. It seems therefore relevant to investigate the influence of temperature on the structures, energetics, and spectral properties of water–PAH systems of interstellar and atmospherical interest. Such properties are difficult to obtain experimentally as forming such species is a challenge, as well as controlling their internal temperature. Many theoretical studies have been published on pure and protonated water clusters [2–4,16–20,16,21,7,22,23]. Some special sizes such as for instance the octamer [2,16–18,24] or the 20mer[19] motivated very detailed investigations. From the structural point of view, for instance, the octamer was exhaustively studied by Jensen et al. [17] who could identify 29 conformers within 45 kJ mol1. The energetically most stable isomers are cubes, differing from one another via the various arrangements of the hydrogen bonds. Obviously finite temperature simulations are needed in such degenerate cases, and a number of Monte Carlo and/or molecular dynamics (MD) simulations have been achieved, either using force fields [9] or, for the smaller species, ab initio or
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
Car–Parinello MD [10,24], sometimes with quantum vibrations [25]. In between force fields and ab initio methods to describe the ground state potential energy surfaces, approximate (albeit still quantal) approaches such as Density Functional Tight Binding (DFTB) [26–29] methods offer an appealing route to finite temperature simulations for reasonably large systems, or whenever intensive simulations are aimed, beyond single trajectory pictures not necessarily representative. Several DFTB studies of water clusters (neutral or protonated) have been performed recently [20,16,21,7,22]. Yu and Qui [22] in particular simulated the selfconsistent-charge (SCC–) DFTB [27,28] dynamical spectra of protonated water clusters and showed that the obtained spectral trends are satisfactory. Finite temperature investigations of the interaction between carbonaceous systems such as PAHs and water molecules or clusters have been much less considered in the literature, and essentially until very recently mainly from a structural point of view (see Ref. [30] and references therein). Recently, Hernandez-Rojas and coworkers [11] reported a rather exhaustive investigation of thermal properties for carbonaceous systems including PAHs using polarizable force fields in the rigid approximation for the carbonaceous parts and the water molecules, focusing on the solid-to-liquid changes and determining the caloric curves from Monte Carlo simulations. The determination of the IR spectra from molecular dynamics for intramolecular modes however obviously makes compulsory to depart from the rigid molecule approximation. Very recently, we achieved a dynamical SCC–DFTB study of small water clusters (n 6 7) interacting with coronene C24H12, a compact PAH, emphasizing the conformational changes and isomerization processes with increasing temperature, demonstrating that the changes could be evidenced via spectral intra-molecular signatures [31]. The present paper reports a dynamical DFTB investigation of finite-temperature structural, energetical and spectral properties of the water octamer (H2O)8 adsorbed on coronene. The water octamer is indeed interesting because of its cubic shape at equilibrium and is a relevant candidate to understand the 3-dimensional (3D) to 2-dimensional (2D) changes induced by temperature in the presence of the interaction with a PAH molecule. Section 2 provides a brief presentation of the DFTB version with the details of the electronic structure calculations and the conditions of the dynamical simulations. Section 3 is devoted to the analysis and discussion of the results, emphasizing in particular on the link between the conformational changes and the evolution of the IR spectra. Finally the conclusion summarizes the results and brings some perspectives.
2. Computational methods In the present paper, we apply the method benchmarked and used for (C24H12)(H2O)n (n = 1–7) complexes [30,31] to the particular case of (C24H12)(H2O)8. Briefly, on-the-fly Born– Oppenheimer molecular dynamics (BOMD) simulations are performed in the microcanonical ensemble. Given the size of the clusters and their floppiness, reasonably long simulation times are needed for a good statistics in probing the overall phase space. This requires an efficient computational method to calculate on-the-fly energies and gradients possibly in a quantum framework. The SCC– DFT method [26–29] is well suited as it is an approximated DFTbased approach significantly faster than DFT due to the use of parametrised integrals and a reduced basis set. Briefly, (1) The Kohn–Sham energy is developed up to the second order around a reference density q0 leading to:
ESCC—DFTB ¼
55
occ X X ^ q j w i þ Erep ðq Þ þ 1 ni hwi j h½ c ðRab Þdqa dqb i 0 0 2 ab ab i
ð1Þ where ^ q is the (i) wi are the Kohn–Sham molecular orbitals, h½ 0 monoelectronic Hamiltonian at the reference density and ni is the molecular orbital occupation number, (ii) Erep(q0) is a repulsive contribution expressed as a sum of atomic-pair terms:
Erep ðq0 Þ ¼
1X U½qa0 ; qb0 2 ab
ð2Þ
(iii) the last term covers the Coulomb interactions between the fluctuating charges located on atoms a and b respectively. The second-order charge density fluctuations are approximated by the atomic charge fluctuations dqa (Mulliken in the initial SCC–DFTB scheme). The cab terms are functions dependent on the inter-atomic distance and the values of the chemical hardness ga of the atom involved. (2) The molecular orbitals are expressed as a minimal (valence) linear combination of atomic orbitals /l functions:
wi ¼
X cil /l l
(3) All three-center contributions are neglected in the Kohn– Sham matrix. We use the SCC–DFTB implementation in the deMonNano code [32] that we modified to treat PAH–water clusters [30]. As explained in these papers, tackling PAH–water clusters requires a delicate description of the balance between long range interactions, namely hydrogen bonds, Coulomb interactions between distant charge fluctuations, London dispersion, polarization and Pauli repulsion. This was achieved via (i) charge corrections: we showed that the long range electrostatic interactions in these systems was improved by replacing the Mulliken definition of charges in the potential by the Class IV – Charge Model 3 (hereafter CM3 charges) developed by the group of Truhlar [33], (ii) inclusion of dispersion contributions adding an empirical term in order to describe dispersion interactions [34–36]:
X C6 0 Edisp ¼ fdamp ðRkk0 Þ 6kk Rkk0 k–k0
ð3Þ
where k and k0 are atomic indices. The C 6ab dispersion coefficient values were taken from the paper by Goursot et al. [37]. These corrections, and particularly the replacement of Mulliken charges by CM3 charges in the initial SCC–DFTB hamiltonian, led to relevant geometries and thermodynamics for (H2O)2,6 and (C6H6)(H2O)1,3 systems in regards of correlated wavefunction calculations and experimental results [30,31]. In the following, we will refer to this modified SCC–DFTB hamiltonian as simply DFTB. The kinetic temperature T of the BOMD/DFTB simulation is defined by:
T¼2
hEk i ð3n 6Þkb
ð4Þ
where n is the number of atoms and kb the Boltzman constant. T is increased dT = 50 K in the range 10–300 K. For a given T, 10 MD simulations of 0.1 ns each, with a timestep of 0.1 fs, are performed. We
56
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
found out that such a small timestep was mandatory to correctly describe the OAH stretching motion of H2O. The starting-point geometry is that of the lowest-energy isomer, with kinetic energies distributed randomly over atomic velocities. Finite-temperature IR spectra are derived from these simulations, computing the Fourier transform of the auto-correlation function of the dipole moment:
aðxÞ / x2
Z
þ1
dthlð0Þ lðtÞi eixt
ð5Þ
0
where l(t) is the dipole moment vector at time t and hl(0)l(t)i a statistical average over 100 16,000 overlapping segments of 32,958 fs each, defined with different (random) initial conditions. For a spectrum at a given T, the statistical averages of the dipole autocorrelation function are collected over 10 100 16,000 segments to calculate the spectra according to Eq. (5). This procedure minimizes the dependency on the initial conditions and should allow a relevant sampling of the energetically favorable portion of the PES. In order to benchmark the BOMD/DFTB algorithm at low temperature vs harmonic hessian diagonalisation, we first checked that the MD spectra obtained at 10 K and the harmonic spectra were similar. 3. Results and discussion 3.1. Dynamical evolution of (C24H12)(H2O)8 In this section, we describe the evolution of (C24H12)(H2O)8 during the dynamics for T between 10 K and 300 K. Comparison with bare (H2O)8 is made. All simulations start from (C24H12)(H2O)8 (1) reported in Fig. 1. This structure results from the optimization
Fig. 1. DFTB optimized structures of low-energy isomers of (C24H12)(H2O)8. The distances (in Å) between the centers of mass of (H2O)8 and C24H12 are reported, as well as the computed relative enthalpies at 0 K (in kJ mol1) with respect to isomer (1).
of (C24H12)(H2O)8 with the most stable (cubic) structure of (H2O)8 (Fig. 4) deposited on C24H12. Local minima optimized from conformers formed during 100 K and 150 K dynamics (quenching) are reported in Figs. 1 and 2, those from 200 K simulations are reported in Fig. 3. The relative enthalpies at 0 K with respect to the most stable conformer (DH(0 K)) are also reported. Isomers of bare (H2O)8 found as conformers adsorbed on C24H12 during the simulations, were also optimized. The structures of the obtained local minima along with DH(0 K) values, are reported in Fig. 4. Regarding our exhaustive studies on (i) the isomers of (H2O)6 for which DFTB DH(0 K) values were compared to CCSD (T) data and (ii) the C24H12(H2O)n (n = 1–10) complexes for which the DFTB structures and intermolecular (H2O)nAC24H12 binding energies were compared to DFT and force-field data [31], these DFTB structures and DH(0 K) values are expected to be reliable. To support our discussion, two parameters are followed: (i) the average distance hdi between the centers of mass of (H2O)8 and C24H12 (cf. Fig. 5), (ii) the Lindeman criterion, namely the fluctuations of the intermolecular O—H distances within (H2O)8 in (C24H12)(H2O)8, quantified by d defined as:
X 2 d¼ nðn 1Þ i
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hR2ij i hRij i2 hRij i
ð6Þ
n is the number of oxygen atoms and Rij the interatomic distances. For infinite systems and clusters, the d parameter can be used as an indicator of phase change. For instance, a sudden increase of d may be indicator of melting (although the signature of phase change in
Fig. 2. DFTB optimized structures of low-energy isomers of (C24H12)(H2O)8. The distances (in Å) between the centers of mass of (H2O)8 and C24H12 are reported, as well as the computed relative enthalpies at 0 K (in kJ mol1) with respect to isomer (1).
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
57
Fig. 5. Average distance (in Å) between the centers of mass of (H2O)8 and C24H12 within (C24H12)(H2O)8 as a function of T. Error bars correspond to the root mean square values.
Fig. 3. DFTB optimized structures of low-energy isomers of (C24H12)(H2O)8. The distances (in Å) between the centers of mass of (H2O)8 and C24H12 are reported, as well as the computed relative enthalpies at 0 K (in kJ mol1) with respect to isomer (1).
Fig. 6. Evolution of d (Lindeman criterion, see text) for (H2O)8 in the (C24H12)(H2O)8 complex (red squares) and in the bare clusters (empty circles) as a function of T. The d values derived from individual simulations are reported, as well as the average hdi values, joint by plain ((C24H12)(H2O)8) or dashed ((H2O)8) lines.
Fig. 4. DFTB optimized structures for (H2O)8. Representative intermolecular distances (in Å) are reported, as well as the computed relative enthalpies at 0 K (in kJ mol1) with respect to isomer (10 ).
finite system is not unique). The values of d for (H2O)8 in (C24H12)(H2O)8 and in bare clusters computed in the present work
are reported in Fig. 6. Furthermore, the distribution of kinetic energy between the intramolecular vibrational modes of (H2O)8 and C24H12, and the intermolecular modes (C24H12)A(H2O)8, as a function of T, is shown in Fig. 7. The reasonable linear fits of the data points shows a correct kinetic energy repartition in (H2O)8 and C24H12. At 10 K, vibration of intramolecular OAH bond not involved in H—OH intermolecular bonds is essentially observed. No conformer change (hdi = 0.0088) and no displacement of (H2O)8 with respect to C24H12 (hdi = 9.6 Å, with rms = 0.007) is observed. At 50 K, similarly, vibration of the free OH bonds is observed, along with CH vibration (mainly out-of-plane bending). No conformer change (hdi = 0.022) occurs, and motions of small amplitude of (H2O)8 with respect to C24H12 are shown (hdi = 9.75 Å, with rms = 0.022).
58
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
Fig. 7. Computed kinetic energy repartition between intramolecular vibrational modes of C24H12 (hEkiC24H12, full circles) and (H2O)8 (hEki(H2O)8, crosses), and intermolecular modes (C24H12)A(H2O)8 (hEkiinter, triangles) for (C24H12)(H2O)8. The data are fitted by linear fits (plain lines).
At 100 K (hEki(H2O)8 32 kJ mol1), the relative motion of (H2O)8 with respect to C24H12 is enhanced, and diffusion starts (hdi = 10.7 Å, rms = 0.1). (H2O)8 mostly maintains its initial cubic structure (hdi = 0.04) and migrates to the edge of C24H12. A local optimization of the resulting observed structure leads to isomer (7) (Fig. 3), located 4.4 kJ mol1 above isomer (1). The diffusion to the edge is sometimes followed by the opening of an intermolecular OAH bond, leading to the formation of isomer (2) (Fig. 1, +0.4 kJ mol1 with respect to (1)). Isomer (6), with a pentagonal ring (+3.5 kJ mol1 with respect to (1)), is also formed. At 150 K (hEki(H2O)8 48 kJ mol1), isomerisation processes are current. This is shown as the range of the d values increases (0.045–0.11), with an average of 0.07 (cf. Fig. 6). Besides, as can be seen in Fig. 5 (hdi = 11.7 Å, rms = 0.1), edge coordination prevails. In particular, the formation of isomer (3), quasi degenerate with isomer (1), is frequently observed. In isomer (3), an alternative cubic form of (H2O)8 (that we will call C2, C1 being the cubic conformer in isomer (1)), resulting from the rotation of the water monomers with respect to C1, is coordinated to the edge of coronene. In the case of bare (H2O)8 cluster, C2 was found to be 14 kJ mol1 higher in energy than C1 (Fig. 4). The edge-coordination to C24H12 stabilises C2 with respect to bare water octamer clusters. Other low-energy isomers for (H2O)8, most of them coordinated to the edge of coronene, also form at 150 K. Several are reported in Figs. 1 and 2. In addition to isomers (2), (6) and (3), isomers (4) and (5) also form. In isomer (4), (H2O)8 can be seen as a prism structure plus two water monomers insuring the coordination with two hydrogen atoms of coronene, whereas in isomer (5), (H2O)8 can also be seen as a prism plus one water monomer coordinated on the edge of coronene, plus one water molecule forming a four-member ring. Interestingly, the corresponding optimized bare (H2O)8 clusters are found 14–27 kJ mol1 higher in energy than the C1 form (Fig. 4). The (mostly edge) coordination to coronene stabilizes the (H2O)8 opened structures with respect to cubic forms. This results in an increase of the density of structures. This also corresponds to a step towards a 3D–2D morphology change. At 200 K (hEki(H2O)8 57 kJ mol1), fast isomerisation into a variety of conformers is observed. This is shown as d values are found
in the range [0.13–0.25] with an average of 0.21. Beside those reported in Figs. 1 and 2, some higher energy isomers form. Two of them (optimized locally) are reported in Fig. 3: Isomer (8) (+6.1 kJ mol1 above (1)), where (H2O)8 adopts a distorted Prism structure with a joint four-membered ring, and isomer (9) (+6.6 kJ mol1 above (1)), where one water molecule departs from the cubic form of (H2O)8 and prefers to coordinate to the edge of C24H12. Consistently with a large d value, the notion of permanent structure is questioned at 200 K. At 250 K (hEki(H2O)8 70 kJ mol1), the notion of permanent structure disappears and the loss of one water molecule is observed (2 over 10 simulations). It is interesting to note that the temperature ranges for isomerisation, ‘‘melting’’ behavior, and dissociation starting processes are similar to those found for (C24H12)(H2O)6,7 [31]. Besides, when comparing the d values of bare and adsorbed (H2O)8, the temperatures for isomerisation starting and phase-change like processes for (H2O)8 are decreased by 50 K upon adsorption (Fig. 6). Overall, the melting temperature for (H2O)8 adsorbed on C24H12 should be in the 150–200 K range and that of bare (H2O)8 in the 200–250 K range. This is consistent with the previously mentioned increase of the density of structures for the lowest-energy isomers of (H2O)8 when they are adsorbed on a PAH such as coronene. This is also in line with the results by Hernandez-Rojas et al. [11] who found that, when adsorbed on Cþ 60 , the melting point for (H2O)8 is lowered by 80 K, from 200 K in bare cluster to 120 K when adsorbed on Cþ 60 . Their results are obtained from heat capacity calculations with parallel tempering Monte Carlo simulations in the rigid approximation using a polarizable force-field potential. One should be aware that the temperature mentioned in the present microcanonical simulations is the kinetic temperature. Moreover since the Lindeman criterion is only a qualitative criterion not univoquely related with a maximum in the heat capacity, it is difficult to achieve fully quantitative comparison with the former work which moreover concerns ions. One may also recall that all the internal degrees of freedom are considered in the present work and thus trap a part of the energy. This might also contribute to a shift in the values of the precursor isomerization and melting threshold temperatures with respect to rigid approximation. 3.2. Finite-temperature IR spectra The finite-temperature IR spectra of bare and adsorbed (H2O)8 were computed with the method described in Section 2. The resulting mid-IR spectra in the regions of the intramolecular water stretching modes are reported in Fig. 8(a) and (b). The positions of the bands discussed here are reported in Table 1. In the following, we first describe the harmonic spectra (similar to the 10 K spectra) providing the description of the main modes. We discuss in particular the influence of the adsorption on C24H12. Secondly, we describe the dynamical evolution as a function of T. 3.2.1. IR spectrum of (H2O)8 at low T: influence of the adsorption on C24H12 To refer to the modes, we use the notation of our previous paper [30]. m1 and m3 refer to the symmetric and antisymmetric stretching modes of H2O, and m2 to the bending modes. In a cluster, softer modes corresponding to the inter-monomer vibrations are also activated, this was shown for instance for (H2O)2 [30]. These are out-of-plane HOH bending (cHOH) and in-plane HOH bending (dHOH) modes, refering to the plane of the water molecule. Besides, we will refer to the well known modes of a PAH molecule that are the CAH stretching (mCH), CAH out-of-plane bending (cCH), CAC stretching (mCC) and CAH in-plane bending (dCH) modes. We will call free H the hydrogen atoms of the water cluster that are not involved in the intermolecular OAH bonds.
59
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
Table 1 Positions (in cm1) and intensities (in parenthesis, in km mol1) of relevant water octamer bands in the DFTB harmonic spectra of (H2O)8 and (C24H12)(H2O)8. Mode
(H2O)8
m3
3925 3925 3925 3925 3805 3784 3784 3756
(0) (165) (165) (27) (877) (495) (495) (0)
3928 3921 3914 3870 3809 3795 3779 3744
(C24H12)(H2O)8 (68) (89) (152) (56) (489) (525) (146) (96)
m1
3552 3543 3543 3535 3362 3362 3323 3291
(0) (300) (300) (653) (1988) (1988) (5) (0)
3560 3549 3546 3518 3409 3371 3323 3264
(51) (331) (253) (206) (1136) (1145) (149) (760)
m2
1619 1611 1598 1598 1565 1564 1563 1563
(0) (278) (103) (103) (70) (0) (74) (74)
1620 1611 1607 1593 1572 1565 1563 1533
(33) (86) (101) (106) (36) (49) (26) (41)
cHOH
897 891 796 796 724 684 684
(223) (0) (168) (168) (197) (452) (452)
897 874 790 776 725 693 668
(98) (52) (95) (58) (152) (364) (255)
dHOH
519 519
(14) (14)
522 511
(29) (11)
mOAH
370 370
(41) (41)
375 367
(45) (28)
(C24H12)(H2O)8 (367 and 375 cm1). Overall, adsorption on a PAH leads to complex spectroscopic patterns as:
Fig. 8. Finite-temperature BOMD/DFTB spectra in the m1 and m3 regions for (a) (H2O)8 starting from isomer (10 ) and (b) (C24H12)(H2O)8 starting from isomer (1). The 10 K spectra are reproduced in dashed lines to compare with higher T spectra.
As can be seen in Table 1, in bare (H2O)8, the m3 modes are found in the 3756–3925 cm1 range, the m1 modes in the 3291– 3552 cm1 range and the m2 modes in the 1563–1619 cm1 range. When (H2O)8 is adsorbed on C24H12, the ranges are slightly broader i.e. they become 3744–3928 cm1 (m3), 3264–3560 cm1 (m1) and 1533–1620 cm1 (m2). For each mode type, the lower energy bands are redshifted while the higher energy bands are blueshifted with respect to bare (H2O)8. This is due to the fact that the bonds that are equivalent in bare (H2O)8 become inequivalent in (C24H12)(H2O)8, leading to degeneracy quenching. For instance, the higher energy bands in the m3 region are assigned to free H. These modes are degenerate in bare (H2O)8 but are split upon adsorption on the PAH. The bands of higher energy (3920 and 3928 cm1) refer to the free H that are not in interaction with C24H12, and those of lower energy (3870 and 3914 cm1) are assigned to the free H that are in interaction with C24H12. In the lower energy part of the spectrum, the OAH intermolecular bonds stretching modes (mOAH in Table 1) are degenerate in the case of bare (H2O)8 (370 cm1) and are split into two bands in
(i) symmetry decreases, resulting in degeneracy quenching as just been explained, (ii) some modes of (H2O)8 couple with the modes of C24H12. For instance, the dHOH modes perturb the cCH mode, leading to the occurrence of two bands at 824 and 826 cm1 instead of one at 824 cm1 in bare C24H12, (iii) at very low energy, new modes of low intensity assigned to the relative motions of (H2O)8 against C24H12 (translations and rotations) are found at 13, 17, 22, 29, 41 and 49 cm1. The band at 49 cm1 in particular is due to the [(H2O)8]A[C24H12] ‘‘stretching’’ motion. Interestingly, only the mCH mode of C24H12 does not couple to any mode of (H2O)8. As mentioned in our previous paper [30], the IR absolute band positions of the water modes computed with the DFTB method have to be scaled to fit experimental data, the scaling factor depending on the mode. This is expected in part due to the inherent errors of DFTB and also as no low-T anharmonic effect due to nuclear quantum effect likely to be crucial when hydrogen atoms are involved, are not taken into account. Regarding water clusters, scaling factors were obtained for (H2O)3, for which experimental data in cryogenic matrix are available [5], in the m1 (0.97), m2 (1.02), and m3 (0.94) regions. Scaling factors are very commonly used for DFT harmonic spectra when comparison with
60
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
experimental data is achieved [38]. We showed that the scaling factors determined for (H2O)3 could be reasonably applied to larger water clusters, we checked on (H2O)6 for which experimental data are available. Besides, the trends in the shifts of the band positions of (H2O) upon adsorption on C6H6 were shown to be correctly described with the DFTB method against experimental data [30]. The spectral trends induced by the adsorption of water clusters on PAHs are therefore expected to be correctly described with our DFTB method. 3.2.2. Evolution of the IR spectra as a function of T The dynamical IR spectra of bare and adsorbed (H2O)8 at several temperatures, in the m1 and m3 regions, are reported in Fig. 8(a) and (b) respectively. As mentionned previously, the dynamical IR spectra at 10 K are similar to the harmonic spectra, with a maximum discrepancy of 2% on the positions. Indeed, no low-T anharmonic effect due to the quantum nature of nuclei, is taken into account in our approach. In our experience [39,40], the expected evolution of uncoupled vibrational modes when T increases is a redshift of the position, along with a broadening of the band. This is the case for instance of the cCH, mCC and mCH modes of C24H12. We showed that the temperature evolution of the redshift for these modes could be fit by linear laws for C24H12 neutral, cationic and complexed with Si+ or Fe+ [39,40], the slope of the linear fit being the so-called anharmonicity factor v. Obviously, in case of isomerization or melting, such perturbative definition of anharmonicity breaks down. Our computed v value for the cCH mode of bare C24H12 (quoted hereafter vcCH ) was found to be in agreement with experimental data [39]. We showed that the adsorption of Si and Fe on C24 Hþ 12 leads to a decrease of vcCH with respect to bare C24 Hþ 12 [39,40]. We also showed that although the adsorption of one H2O molecule on C24H12 hardly affects vcCH , the adsorption of (H2O)2 leads to an increase of vcCH with respect to bare C24H12, and we proposed that this could be regarded as a signature for edge-coordination [30]. In the case of (H2O)8 (bare and adsorbed), a simple description of the evolution of the band positions in the m1 and m3 regions, and thus the determination of a simple linear law to quantify it, is impossible. This is expected as, due to the high density of bands, modes are likely to couple. Besides, at 150 K for (C24H12)(H2O)8 and 200 K for (H2O)8, isomerisations processes start, while a liquid-like behavior is observed at 200 K and 300 K for (C24H12)(H2O)8 and (H2O)8 respectively. As can be seen in Fig. 8(a) and (b), broadening, merge and blueshifts are observed. From 10 K to 50 K, broadenings are observed. The most energetic m3 band is sligthly redshifted (-2 cm1 for bare (H2O)8). The other bands are blueshifted, the m3 lowest energy band within bare (H2O)8 by + 9 cm1, the m1 highest (resp. lowest) energy band by +5 cm1 (resp. +13 cm1). From 50 K to 100 K, enhanced broadenings and blueshifts occur. The m3 lowest energy band for bare (H2O)8 is indeed shifted by +20 cm1 while the m1 highest (resp. lowest) energy bands is shifted by +11 cm1 (resp. +17 cm1). Interestingly, the most shifted modes are assigned to OH bonds with H atoms involved in intermolecular bonds. The evolution of the IR spectrum of (C24H12)(H2O)8 is similar and the phenomena (broadenings, merge and blueshifts) are enhanced. The spectra at 200 K for bare (H2O)8, and at 150 K for (C24H12)(H2O)8, can be regarded as ‘‘transitory’’ as isomerisation processes start. In the case of bare (H2O)8, two groups of bands can be distinguished for each mode type (m1 and m3), that could be assigned to free OH on the one hand (higher energy m3 band, lower energy m1 band), and to OH with H involved in intermolecular OAH bonds on the other hand (lower energy m3 band, higher energy m1 band). In the case of (C24H12)(H2O)8, the situation is more complicated as other OH bond types are involved (H coordinated to the PAH surface). Finally, at 300 K for bare (H2O)8, and at
200 K for (C24H12)(H2O)8, only one broad m1 band (with maxima at 3523 and 3539 cm1 respectively) and one broad m3 band (with maxima at 3877 and 3883 cm1 respectively) are found: all hydrogen atoms are equivalent, in line with a liquid-phase like regime. When comparing the finite-temperature IR spectra of (C24H12)0,1(H2O)6 published in Ref. [31] with those of (C24H12)0,1(H2O)8 (this work), It may be interesting to mention that, at any temperature, even when the ‘‘melting’’ point is reached, the spectrum of (C24H12)0,1(H2O)6 [31] is clearly different from that of (C24H12)0,1(H2O)8. This is obviously related to the finite-size of the system.
4. Conclusions The present paper provides a detailed study of finite-temperature BOMD/DFTB simulations in the range 10–300 K starting from the lowest energy structure of the complex (H2O)8C24H12. For T 6 50 K, only small amplitude vibrations are observed. At 100 K, diffusion of (H2O)8 on coronene starts, with sporadic isomerisation of (H2O)8. Above 150 K, isomerization processes occur, and at slightly higher T, a melting-like behavior of (H2O)8 is observed. At 250 K, loss of one water molecule is observed. This evolution is quantified by the determination of a chosen set of indicators, and visible on the finite-temperature IR spectra derived from the dynamics. A comparison between the dynamical behavior of bare and adsorbed (H2O)8 is discussed. We show in particular that (H2O)8, when adsorbed on C24H12 isomerizes at a lower T than under its bare form, and that it is consistent with a larger density of structures induced by the adsorption on the PAH. From this work, we infer that the melting temperature for (H2O)8 adsorbed on C24H12 should be in the 150–200 K range and that of bare (H2O)8 higher by 50 K, which is somewhat consistent with the results obtained by Hernandez-Rojas et al. [11] for cationic water systems from parallel tempering Monte Carlo simulations using the rigid molecule approximation. Thanks to the inclusion of all degrees of freedom in the simulations, we also show the influence of the adsorption of (H2O)8 (in its most stable cubic conformation) on the IR spectrum at low T. As expected, the IR spectrum gets more complex as (i) degeneracy is quenched leading to band splittings, (ii) new modes are activated, and (iii) some modes of the water cluster couple to the modes of coronene. The analysis of finite-T IR spectra shows that at all T (prior to the evaporation of the cluster), the spectra of bare and adsorbed (H2O)8 can be distinguished. It is also interesting to quote that, at any temperature, even when the ‘‘melting’’ threshold is reached, the spectrum of (C24H12)0,1(H2O)8 is clearly different from that of (C24H12)0,1(H2O)6, and thus a signature of the structure of the (finite)-size cluster. The present work certainly calls for several extensions. The most important one should concern the quantum character of the vibration which plays a non-negligible role in the case of water clusters [25]. Isomer interconversion will be affected energetically and entropically by quantum delocalization, account of vibrational zero point energy (ZPE) and tunneling effect. Differences in the ZPEs of various isomers will certainly influence the conformational changes. Wang et al. [25] have recently investigated the isomer populations in the isolated water hexamer using Path Integral Molecular Dynamics (PIMD). They found that while in the classical simulations, only a prism-type isomer was populated at low temperature (T = 30 K), decaying into a cage-type isomer for T close to 100 K, and that alternatively in the quantum simulation both isomers already coexist with a 60/40% ratio at T = 30 K. Quantum effects will contribute to a decrease in the melting temperature. Incorporating quantum effects could be targeted by combining the present DFTB code with a PIMD approach. The quantum
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
character influence on the IR spectra could also certainly be examined using techniques extending the use of the autocorrelation functions of the dipole moments, as done recently on isolated PAH molecules, either following the dynamics of the replica (or polymers), or that of the centroids [41]. Account of quantum vibration is expected to shift and broaden the bands, and possibly change some of the peak intensities. This may be influential in the case of the water octamer, where numerous low energy isomers contribute to the IR spectra even at low temperature. Among other possible relevant theoretical extensions, it would also be interesting to determine the heat capacities for the present system, as was done by Choi [16] for pure water clusters. It would also certainly be relevant to extend the size of the carbonaceous cluster, in particular to examine the size effects (number of vibrations involved) in the coupling between the intramolecular and intermolecular degrees of freedom. In this direction, another development would also be to use thermostatted MD in order to approach the evolution of a surface with a given canonical temperature. Such perspectives would require increased computing resources. References [1] R. Ludwig, Water: from clusters to the bulk, Angew. Chem. Int. Ed. 40 (2001) 1808–1827. [2] C.J. Gruenloh, J.R. Carney, C.A. Arrington, T.S. Zwier, S.Y. Fredericks, K.D. Jordan, Infrared spectrum of a molecular ice cube: the S-4 and D-2d water octamers in benzene–(water)8, Science 276 (1997) 1678–1681. [3] A.D. Kulkarni, R.K. Pathak, L.J. Bartolotti, Structures, energetics, and vibrational spectra of H2O2 (H2O)n, n = 1–6 clusters: Ab initio quantum chemical investigations, J. Phys. Chem. A 109 (2005) 4583–4590. [4] A.D. Kulkarni, S.R. Gadre, S. Nagase, Quantum chemical and electrostatic studies of anionic water clusters H2 O n , J. Mol. Struct. Theochem 851 (2008) 213–219. [5] B. Tremblay, B. Madebène, M. Alikhani, J. Perchard, The vibrational spectrum of the water trimer: Comparison between anharmonic ab initio calculations and neon matrix infrared data between 11,000 and 90 cm1, Chem. Phys. 378 (2010) 27–36. [6] J. Ceponkus, P. Uvdal, B. Nelander, Water tetramer, pentamer, and hexamer in inert matrices, J. Phys. Chem. A 116 (2012) 4842–4850. [7] P. Miro, C.J. Cramer, Water clusters to nanodrops: a tight-binding density functional study, Phys. Chem. Chem. Phys. 15 (2013) 1837–1843. [8] A.D. Kulkarni, K. Babu, S.R. Gadre, L.J. Bartolotti, Exploring hydration patterns of aldehydes and amides: ab initio investigations, J. Phys. Chem. A 108 (2004) 2492–2498. [9] J. Douady, F. Calvo, F. Spiegelman, Effect of an ionic impurity on the caloric curves of water clusters, Eur. Phys. J. D 52 (2009) 47–50. [10] N.J. Singh, M. Park, S.K. Min, S.B. Suh, K.S. Kim, Magic and antimagic protonated water clusters: Exotic structures with unusual dynamic effects, Angew. Chem. Int. Ed. 45 (2006) 3795–3800. [11] J. Hernandez-Rojas, F. Calvo, F. Rabilloud, J. Breton, J.M. Gomez Llorente, Modeling water clusters on cationic carbonaceous seeds, J. Phys. Chem. A 114 (2010) 7267–7274. [12] B. Finlayson-Pitts, J. Pitts (Eds.), Atmospheric Chemistry, Fundamental and Experimental Techniques, WILEY-VCH, Verlag Berlin GmbH, 1986. [13] L.J. Allamandola, A.G.G.M. Tielens, J.R. Barker, Polycyclic aromatichydrocarbons and the unidentified infrared-emission bands – auto exhaust along the milky-way, Astrophys. J. 290 (1985) L25–L28. [14] A. Léger, J.L. Puget, Identification of the unidentified ir emission features of interstellar dust, Astron. Astrophys. 137 (1984) L5–L8. [15] A.G.G.M. Tielens, Interstellar polycyclic aromatic hydrocarbon molecules, An. Rev. Astron. Astrophys. 46 (2008) 289–337. [16] T.H. Choi, Simulation of the (H2O)8 cluster with the SCC–DFTB electronic structure method, Chem. Phys. Lett. 543 (2012) 45–49. [17] J.O. Jensen, P.N. Krishnan, L.A. Burke, Theoretical-study of water clusters – octamer, Chem. Phys. Lett. 246 (1995) 13–19. [18] J.S. Kim, B.J. Mhin, S.J. Lee, K.S. Kim, Entropy-driven structures of the water octamer, Chem. Phys. Lett. 219 (1994) 243–246. [19] G.S. Fanourgakis, E. Apra, W.A. de Jong, S.S. Xantheas, High-level ab initio calculations for the four low-lying families of minima of (H2O)20. II.
[20] [21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39]
[40]
[41]
61
Spectroscopic signatures of the dodecahedron, fused cubes, face-sharing pentagonal prisms, and edge-sharing pentagonal prisms hydrogen bonding networks, J. Chem. Phys. 122 (2005) 134304. T.H. Choi, K.D. Jordan, Application of the SCC–DFTB Method to H+(H2O)6, H+(H2O)21, and H+(H2O)22, J. Phys. Chem. B 114 (2010) 6932–6936. P. Goyal, M. Elstner, Q. Cui, Application of the SCC–DFTB method to neutral and protonated water clusters and bulk water, J. Phys. Chem. B 115 (2011) 6790– 6805. H. Yu, Q. Cui, The vibrational spectra of protonated water clusters: a benchmark for self-consistent-charge density-functional tight binding, J. Chem. Phys. 127 (2007) 234504. S. Maheshwary, N. Patel, N. Sathyamurthy, A.D. Kulkarni, S.R. Gadre, Structure and stability of water clusters (H2O)n, n = 8–20: an ab initio investigation, J. Phys. Chem. A 105 (2001) 10525–10537. S. Karthikeyan, M. Park, I. Shin, K.S. Kim, Structure, stability, thermodynamic properties, and infrared spectra of the protonated water octamer H+(H2O)8, J. Phys. Chem. A 112 (2008) 10120–10124. Y. Wang, V. Babin, J.M. Bowman, F. Paesani, The water hexamer: cage, prism, or both. Full dimensional quantum simulations say both, J. Am. Chem. Soc. 134 (2012) 11116–11119. D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, R. Kaschner, Construction of tight-binding-like potentials on the basis of density-functional theory – application to carbon, Phys. Rev. B 51 (1995) 12947–12957. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, G. Seifert, Self-consistent-charge density-functional tight-binding method for simulations of complex material properties, Phys. Rev. B 58 (1998) 7260–7268. T. Frauenheim, G. Seifert, M. Elsterner, Z. Hajnal, G. Jungnickel, D. Porezag, S. Suhai, R. Scholz, A self-consistent charge density-functional based tightbinding method for predictive materials simulations in physics, chemistry and biology, Phys. Status Solidi (b) 217 (2000) 41–62. A. Oliveira, G. Seifert, T. Heine, H. duarte, Density-functional based tightbinding: an approximate DFT method, J. Braz. Chem. Soc. 20 (2009) 1193– 1205. A. Simon, M. Rapacioli, J. Mascetti, F. Spiegelman, Vibrational spectroscopy and molecular dynamics of water monomers and dimers adsorbed on polycyclic aromatic hydrocarbons, Phys. Chem. Chem. Phys. 14 (2012) 6771–6786. A. Simon, F. Spiegelman, Water clusters adsorbed on polycyclic aromatic hydrocarbons: energetics and conformational dynamics, J. Chem. Phys. 138 (2013) 194309. T. Heine, M. Rapacioli, S. Patchkovskii, J. Frenzel, A. Koster, P. Calaminici, H.A. Duarte, S. Escalante, R. Flores-Moreno, A. Goursot, J. Reveles, D. Salahub, A. Vela, demon-nano experiment (2009).
Contents lists available at SciVerse ScienceDirect
Computational and Theoretical Chemistry journal homepage: www.elsevier.com/locate/comptc
Conformational dynamics and finite-temperature infrared spectra of the water octamer adsorbed on coronene Aude Simon ⇑, Fernand Spiegelman Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, Université de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France
a r t i c l e
i n f o
Article history: Received 17 May 2013 Received in revised form 18 June 2013 Accepted 18 June 2013 Available online 29 June 2013 Keywords: Water clusters PAH Molecular dynamics DFTB Infrared
a b s t r a c t We present in this work a detailed study of all-atoms conformational dynamics and finite temperature infrared (IR) spectra of the water octamer (H2O)8 adsorbed on coronene (C24H12), a compact Polycyclic Aromatic Hydrocarbon (PAH). The potential energy surface is calculated on-the-fly within the self-consistent-charge density-functional based tight-binding (SCC–DFTB) approach, with modifications insuring the correct description of water–water and water–PAH interactions. On-the-fly Born–Oppenheimer molecular dynamics simulations are performed in the temperature T range 10–300 K starting from the most stable cubic conformer of (H2O)8 adsorbed on C24H12. As T increases, diffusion of (H2O)8 starts, followed by isomerisation processes and ‘‘melting’’ of (H2O)8 for T in the range 150–200 K. These dynamical processes are quantified through the temperature evolution of chosen indicators. We show that the adsorption of (H2O)8 on C24H12 leads to lower isomerisation temperatures than for the bare cluster, in line with an increase of the density of isomers. A comparative analysis of the low-T (harmonic) spectra of bare and adsorbed (H2O)8 is provided. Finally, the evolution of the spectra as a function of T, illustrating the conformational dynamics, is discussed. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Water clusters have raised a growing and fundamental interest [1] in several aspects (i) as finite size systems with specific properties [2–6], (ii) as precursors to understand the bottom-up building properties of liquid and ice water [7], and (iii) to understand nanoscale solvation [8], a feature essential in biological systems. The interaction of water molecules and clusters impurities is also of wide interest in atmospheric chemistry, for instance in nucleative processes initiating cloud formation, and various investigations have considered the interaction between water droplets or clusters and impurities such as an extra proton [9,10], or carbon-containing molecular impurities [11]. Among the latter, Polyaromatic Aromatic Hydrocarbon (PAH) molecules show significant abundance in the atmosphere, due to their efficient formation as by-products of natural processes (biomass burning) or of human activities (combustion of fossil fuels) [12]. PAH molecules and derived-PAH species have also raised a strong and lasting interest since they were proposed in the mid-eighties [13,14] as the carriers of the aromatic infrared (IR) bands, that are a set of mid-IR emission bands observed in many regions of the interstellar medium. Despite many experimental and theoretical studies, a specific molecule has not been identified yet [15]. These interstellar PAH ⇑ Corresponding author. Tel.: +33 561556033. E-mail address: [email protected] (A. Simon). 2210-271X/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.comptc.2013.06.023
molecules would form at the border of molecular clouds, in Photodissociation Regions (PDRs). The importance of PAH species, under various forms, in interstellar carbonaceous (very small) grains and their role in the interstellar chemistry is still in debate. Inside molecular clouds, temperatures are expected to be in the range 20–100 K and in PDRs, temperatures can reach between 100 K and 1000 K. Regarding atmospherical conditions, the temperatures in the troposphere, that contains an important quantity of water, are expected to be in the range 200–300 K. It seems therefore relevant to investigate the influence of temperature on the structures, energetics, and spectral properties of water–PAH systems of interstellar and atmospherical interest. Such properties are difficult to obtain experimentally as forming such species is a challenge, as well as controlling their internal temperature. Many theoretical studies have been published on pure and protonated water clusters [2–4,16–20,16,21,7,22,23]. Some special sizes such as for instance the octamer [2,16–18,24] or the 20mer[19] motivated very detailed investigations. From the structural point of view, for instance, the octamer was exhaustively studied by Jensen et al. [17] who could identify 29 conformers within 45 kJ mol1. The energetically most stable isomers are cubes, differing from one another via the various arrangements of the hydrogen bonds. Obviously finite temperature simulations are needed in such degenerate cases, and a number of Monte Carlo and/or molecular dynamics (MD) simulations have been achieved, either using force fields [9] or, for the smaller species, ab initio or
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
Car–Parinello MD [10,24], sometimes with quantum vibrations [25]. In between force fields and ab initio methods to describe the ground state potential energy surfaces, approximate (albeit still quantal) approaches such as Density Functional Tight Binding (DFTB) [26–29] methods offer an appealing route to finite temperature simulations for reasonably large systems, or whenever intensive simulations are aimed, beyond single trajectory pictures not necessarily representative. Several DFTB studies of water clusters (neutral or protonated) have been performed recently [20,16,21,7,22]. Yu and Qui [22] in particular simulated the selfconsistent-charge (SCC–) DFTB [27,28] dynamical spectra of protonated water clusters and showed that the obtained spectral trends are satisfactory. Finite temperature investigations of the interaction between carbonaceous systems such as PAHs and water molecules or clusters have been much less considered in the literature, and essentially until very recently mainly from a structural point of view (see Ref. [30] and references therein). Recently, Hernandez-Rojas and coworkers [11] reported a rather exhaustive investigation of thermal properties for carbonaceous systems including PAHs using polarizable force fields in the rigid approximation for the carbonaceous parts and the water molecules, focusing on the solid-to-liquid changes and determining the caloric curves from Monte Carlo simulations. The determination of the IR spectra from molecular dynamics for intramolecular modes however obviously makes compulsory to depart from the rigid molecule approximation. Very recently, we achieved a dynamical SCC–DFTB study of small water clusters (n 6 7) interacting with coronene C24H12, a compact PAH, emphasizing the conformational changes and isomerization processes with increasing temperature, demonstrating that the changes could be evidenced via spectral intra-molecular signatures [31]. The present paper reports a dynamical DFTB investigation of finite-temperature structural, energetical and spectral properties of the water octamer (H2O)8 adsorbed on coronene. The water octamer is indeed interesting because of its cubic shape at equilibrium and is a relevant candidate to understand the 3-dimensional (3D) to 2-dimensional (2D) changes induced by temperature in the presence of the interaction with a PAH molecule. Section 2 provides a brief presentation of the DFTB version with the details of the electronic structure calculations and the conditions of the dynamical simulations. Section 3 is devoted to the analysis and discussion of the results, emphasizing in particular on the link between the conformational changes and the evolution of the IR spectra. Finally the conclusion summarizes the results and brings some perspectives.
2. Computational methods In the present paper, we apply the method benchmarked and used for (C24H12)(H2O)n (n = 1–7) complexes [30,31] to the particular case of (C24H12)(H2O)8. Briefly, on-the-fly Born– Oppenheimer molecular dynamics (BOMD) simulations are performed in the microcanonical ensemble. Given the size of the clusters and their floppiness, reasonably long simulation times are needed for a good statistics in probing the overall phase space. This requires an efficient computational method to calculate on-the-fly energies and gradients possibly in a quantum framework. The SCC– DFT method [26–29] is well suited as it is an approximated DFTbased approach significantly faster than DFT due to the use of parametrised integrals and a reduced basis set. Briefly, (1) The Kohn–Sham energy is developed up to the second order around a reference density q0 leading to:
ESCC—DFTB ¼
55
occ X X ^ q j w i þ Erep ðq Þ þ 1 ni hwi j h½ c ðRab Þdqa dqb i 0 0 2 ab ab i
ð1Þ where ^ q is the (i) wi are the Kohn–Sham molecular orbitals, h½ 0 monoelectronic Hamiltonian at the reference density and ni is the molecular orbital occupation number, (ii) Erep(q0) is a repulsive contribution expressed as a sum of atomic-pair terms:
Erep ðq0 Þ ¼
1X U½qa0 ; qb0 2 ab
ð2Þ
(iii) the last term covers the Coulomb interactions between the fluctuating charges located on atoms a and b respectively. The second-order charge density fluctuations are approximated by the atomic charge fluctuations dqa (Mulliken in the initial SCC–DFTB scheme). The cab terms are functions dependent on the inter-atomic distance and the values of the chemical hardness ga of the atom involved. (2) The molecular orbitals are expressed as a minimal (valence) linear combination of atomic orbitals /l functions:
wi ¼
X cil /l l
(3) All three-center contributions are neglected in the Kohn– Sham matrix. We use the SCC–DFTB implementation in the deMonNano code [32] that we modified to treat PAH–water clusters [30]. As explained in these papers, tackling PAH–water clusters requires a delicate description of the balance between long range interactions, namely hydrogen bonds, Coulomb interactions between distant charge fluctuations, London dispersion, polarization and Pauli repulsion. This was achieved via (i) charge corrections: we showed that the long range electrostatic interactions in these systems was improved by replacing the Mulliken definition of charges in the potential by the Class IV – Charge Model 3 (hereafter CM3 charges) developed by the group of Truhlar [33], (ii) inclusion of dispersion contributions adding an empirical term in order to describe dispersion interactions [34–36]:
X C6 0 Edisp ¼ fdamp ðRkk0 Þ 6kk Rkk0 k–k0
ð3Þ
where k and k0 are atomic indices. The C 6ab dispersion coefficient values were taken from the paper by Goursot et al. [37]. These corrections, and particularly the replacement of Mulliken charges by CM3 charges in the initial SCC–DFTB hamiltonian, led to relevant geometries and thermodynamics for (H2O)2,6 and (C6H6)(H2O)1,3 systems in regards of correlated wavefunction calculations and experimental results [30,31]. In the following, we will refer to this modified SCC–DFTB hamiltonian as simply DFTB. The kinetic temperature T of the BOMD/DFTB simulation is defined by:
T¼2
hEk i ð3n 6Þkb
ð4Þ
where n is the number of atoms and kb the Boltzman constant. T is increased dT = 50 K in the range 10–300 K. For a given T, 10 MD simulations of 0.1 ns each, with a timestep of 0.1 fs, are performed. We
56
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
found out that such a small timestep was mandatory to correctly describe the OAH stretching motion of H2O. The starting-point geometry is that of the lowest-energy isomer, with kinetic energies distributed randomly over atomic velocities. Finite-temperature IR spectra are derived from these simulations, computing the Fourier transform of the auto-correlation function of the dipole moment:
aðxÞ / x2
Z
þ1
dthlð0Þ lðtÞi eixt
ð5Þ
0
where l(t) is the dipole moment vector at time t and hl(0)l(t)i a statistical average over 100 16,000 overlapping segments of 32,958 fs each, defined with different (random) initial conditions. For a spectrum at a given T, the statistical averages of the dipole autocorrelation function are collected over 10 100 16,000 segments to calculate the spectra according to Eq. (5). This procedure minimizes the dependency on the initial conditions and should allow a relevant sampling of the energetically favorable portion of the PES. In order to benchmark the BOMD/DFTB algorithm at low temperature vs harmonic hessian diagonalisation, we first checked that the MD spectra obtained at 10 K and the harmonic spectra were similar. 3. Results and discussion 3.1. Dynamical evolution of (C24H12)(H2O)8 In this section, we describe the evolution of (C24H12)(H2O)8 during the dynamics for T between 10 K and 300 K. Comparison with bare (H2O)8 is made. All simulations start from (C24H12)(H2O)8 (1) reported in Fig. 1. This structure results from the optimization
Fig. 1. DFTB optimized structures of low-energy isomers of (C24H12)(H2O)8. The distances (in Å) between the centers of mass of (H2O)8 and C24H12 are reported, as well as the computed relative enthalpies at 0 K (in kJ mol1) with respect to isomer (1).
of (C24H12)(H2O)8 with the most stable (cubic) structure of (H2O)8 (Fig. 4) deposited on C24H12. Local minima optimized from conformers formed during 100 K and 150 K dynamics (quenching) are reported in Figs. 1 and 2, those from 200 K simulations are reported in Fig. 3. The relative enthalpies at 0 K with respect to the most stable conformer (DH(0 K)) are also reported. Isomers of bare (H2O)8 found as conformers adsorbed on C24H12 during the simulations, were also optimized. The structures of the obtained local minima along with DH(0 K) values, are reported in Fig. 4. Regarding our exhaustive studies on (i) the isomers of (H2O)6 for which DFTB DH(0 K) values were compared to CCSD (T) data and (ii) the C24H12(H2O)n (n = 1–10) complexes for which the DFTB structures and intermolecular (H2O)nAC24H12 binding energies were compared to DFT and force-field data [31], these DFTB structures and DH(0 K) values are expected to be reliable. To support our discussion, two parameters are followed: (i) the average distance hdi between the centers of mass of (H2O)8 and C24H12 (cf. Fig. 5), (ii) the Lindeman criterion, namely the fluctuations of the intermolecular O—H distances within (H2O)8 in (C24H12)(H2O)8, quantified by d defined as:
X 2 d¼ nðn 1Þ i
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hR2ij i hRij i2 hRij i
ð6Þ
n is the number of oxygen atoms and Rij the interatomic distances. For infinite systems and clusters, the d parameter can be used as an indicator of phase change. For instance, a sudden increase of d may be indicator of melting (although the signature of phase change in
Fig. 2. DFTB optimized structures of low-energy isomers of (C24H12)(H2O)8. The distances (in Å) between the centers of mass of (H2O)8 and C24H12 are reported, as well as the computed relative enthalpies at 0 K (in kJ mol1) with respect to isomer (1).
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
57
Fig. 5. Average distance (in Å) between the centers of mass of (H2O)8 and C24H12 within (C24H12)(H2O)8 as a function of T. Error bars correspond to the root mean square values.
Fig. 3. DFTB optimized structures of low-energy isomers of (C24H12)(H2O)8. The distances (in Å) between the centers of mass of (H2O)8 and C24H12 are reported, as well as the computed relative enthalpies at 0 K (in kJ mol1) with respect to isomer (1).
Fig. 6. Evolution of d (Lindeman criterion, see text) for (H2O)8 in the (C24H12)(H2O)8 complex (red squares) and in the bare clusters (empty circles) as a function of T. The d values derived from individual simulations are reported, as well as the average hdi values, joint by plain ((C24H12)(H2O)8) or dashed ((H2O)8) lines.
Fig. 4. DFTB optimized structures for (H2O)8. Representative intermolecular distances (in Å) are reported, as well as the computed relative enthalpies at 0 K (in kJ mol1) with respect to isomer (10 ).
finite system is not unique). The values of d for (H2O)8 in (C24H12)(H2O)8 and in bare clusters computed in the present work
are reported in Fig. 6. Furthermore, the distribution of kinetic energy between the intramolecular vibrational modes of (H2O)8 and C24H12, and the intermolecular modes (C24H12)A(H2O)8, as a function of T, is shown in Fig. 7. The reasonable linear fits of the data points shows a correct kinetic energy repartition in (H2O)8 and C24H12. At 10 K, vibration of intramolecular OAH bond not involved in H—OH intermolecular bonds is essentially observed. No conformer change (hdi = 0.0088) and no displacement of (H2O)8 with respect to C24H12 (hdi = 9.6 Å, with rms = 0.007) is observed. At 50 K, similarly, vibration of the free OH bonds is observed, along with CH vibration (mainly out-of-plane bending). No conformer change (hdi = 0.022) occurs, and motions of small amplitude of (H2O)8 with respect to C24H12 are shown (hdi = 9.75 Å, with rms = 0.022).
58
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
Fig. 7. Computed kinetic energy repartition between intramolecular vibrational modes of C24H12 (hEkiC24H12, full circles) and (H2O)8 (hEki(H2O)8, crosses), and intermolecular modes (C24H12)A(H2O)8 (hEkiinter, triangles) for (C24H12)(H2O)8. The data are fitted by linear fits (plain lines).
At 100 K (hEki(H2O)8 32 kJ mol1), the relative motion of (H2O)8 with respect to C24H12 is enhanced, and diffusion starts (hdi = 10.7 Å, rms = 0.1). (H2O)8 mostly maintains its initial cubic structure (hdi = 0.04) and migrates to the edge of C24H12. A local optimization of the resulting observed structure leads to isomer (7) (Fig. 3), located 4.4 kJ mol1 above isomer (1). The diffusion to the edge is sometimes followed by the opening of an intermolecular OAH bond, leading to the formation of isomer (2) (Fig. 1, +0.4 kJ mol1 with respect to (1)). Isomer (6), with a pentagonal ring (+3.5 kJ mol1 with respect to (1)), is also formed. At 150 K (hEki(H2O)8 48 kJ mol1), isomerisation processes are current. This is shown as the range of the d values increases (0.045–0.11), with an average of 0.07 (cf. Fig. 6). Besides, as can be seen in Fig. 5 (hdi = 11.7 Å, rms = 0.1), edge coordination prevails. In particular, the formation of isomer (3), quasi degenerate with isomer (1), is frequently observed. In isomer (3), an alternative cubic form of (H2O)8 (that we will call C2, C1 being the cubic conformer in isomer (1)), resulting from the rotation of the water monomers with respect to C1, is coordinated to the edge of coronene. In the case of bare (H2O)8 cluster, C2 was found to be 14 kJ mol1 higher in energy than C1 (Fig. 4). The edge-coordination to C24H12 stabilises C2 with respect to bare water octamer clusters. Other low-energy isomers for (H2O)8, most of them coordinated to the edge of coronene, also form at 150 K. Several are reported in Figs. 1 and 2. In addition to isomers (2), (6) and (3), isomers (4) and (5) also form. In isomer (4), (H2O)8 can be seen as a prism structure plus two water monomers insuring the coordination with two hydrogen atoms of coronene, whereas in isomer (5), (H2O)8 can also be seen as a prism plus one water monomer coordinated on the edge of coronene, plus one water molecule forming a four-member ring. Interestingly, the corresponding optimized bare (H2O)8 clusters are found 14–27 kJ mol1 higher in energy than the C1 form (Fig. 4). The (mostly edge) coordination to coronene stabilizes the (H2O)8 opened structures with respect to cubic forms. This results in an increase of the density of structures. This also corresponds to a step towards a 3D–2D morphology change. At 200 K (hEki(H2O)8 57 kJ mol1), fast isomerisation into a variety of conformers is observed. This is shown as d values are found
in the range [0.13–0.25] with an average of 0.21. Beside those reported in Figs. 1 and 2, some higher energy isomers form. Two of them (optimized locally) are reported in Fig. 3: Isomer (8) (+6.1 kJ mol1 above (1)), where (H2O)8 adopts a distorted Prism structure with a joint four-membered ring, and isomer (9) (+6.6 kJ mol1 above (1)), where one water molecule departs from the cubic form of (H2O)8 and prefers to coordinate to the edge of C24H12. Consistently with a large d value, the notion of permanent structure is questioned at 200 K. At 250 K (hEki(H2O)8 70 kJ mol1), the notion of permanent structure disappears and the loss of one water molecule is observed (2 over 10 simulations). It is interesting to note that the temperature ranges for isomerisation, ‘‘melting’’ behavior, and dissociation starting processes are similar to those found for (C24H12)(H2O)6,7 [31]. Besides, when comparing the d values of bare and adsorbed (H2O)8, the temperatures for isomerisation starting and phase-change like processes for (H2O)8 are decreased by 50 K upon adsorption (Fig. 6). Overall, the melting temperature for (H2O)8 adsorbed on C24H12 should be in the 150–200 K range and that of bare (H2O)8 in the 200–250 K range. This is consistent with the previously mentioned increase of the density of structures for the lowest-energy isomers of (H2O)8 when they are adsorbed on a PAH such as coronene. This is also in line with the results by Hernandez-Rojas et al. [11] who found that, when adsorbed on Cþ 60 , the melting point for (H2O)8 is lowered by 80 K, from 200 K in bare cluster to 120 K when adsorbed on Cþ 60 . Their results are obtained from heat capacity calculations with parallel tempering Monte Carlo simulations in the rigid approximation using a polarizable force-field potential. One should be aware that the temperature mentioned in the present microcanonical simulations is the kinetic temperature. Moreover since the Lindeman criterion is only a qualitative criterion not univoquely related with a maximum in the heat capacity, it is difficult to achieve fully quantitative comparison with the former work which moreover concerns ions. One may also recall that all the internal degrees of freedom are considered in the present work and thus trap a part of the energy. This might also contribute to a shift in the values of the precursor isomerization and melting threshold temperatures with respect to rigid approximation. 3.2. Finite-temperature IR spectra The finite-temperature IR spectra of bare and adsorbed (H2O)8 were computed with the method described in Section 2. The resulting mid-IR spectra in the regions of the intramolecular water stretching modes are reported in Fig. 8(a) and (b). The positions of the bands discussed here are reported in Table 1. In the following, we first describe the harmonic spectra (similar to the 10 K spectra) providing the description of the main modes. We discuss in particular the influence of the adsorption on C24H12. Secondly, we describe the dynamical evolution as a function of T. 3.2.1. IR spectrum of (H2O)8 at low T: influence of the adsorption on C24H12 To refer to the modes, we use the notation of our previous paper [30]. m1 and m3 refer to the symmetric and antisymmetric stretching modes of H2O, and m2 to the bending modes. In a cluster, softer modes corresponding to the inter-monomer vibrations are also activated, this was shown for instance for (H2O)2 [30]. These are out-of-plane HOH bending (cHOH) and in-plane HOH bending (dHOH) modes, refering to the plane of the water molecule. Besides, we will refer to the well known modes of a PAH molecule that are the CAH stretching (mCH), CAH out-of-plane bending (cCH), CAC stretching (mCC) and CAH in-plane bending (dCH) modes. We will call free H the hydrogen atoms of the water cluster that are not involved in the intermolecular OAH bonds.
59
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
Table 1 Positions (in cm1) and intensities (in parenthesis, in km mol1) of relevant water octamer bands in the DFTB harmonic spectra of (H2O)8 and (C24H12)(H2O)8. Mode
(H2O)8
m3
3925 3925 3925 3925 3805 3784 3784 3756
(0) (165) (165) (27) (877) (495) (495) (0)
3928 3921 3914 3870 3809 3795 3779 3744
(C24H12)(H2O)8 (68) (89) (152) (56) (489) (525) (146) (96)
m1
3552 3543 3543 3535 3362 3362 3323 3291
(0) (300) (300) (653) (1988) (1988) (5) (0)
3560 3549 3546 3518 3409 3371 3323 3264
(51) (331) (253) (206) (1136) (1145) (149) (760)
m2
1619 1611 1598 1598 1565 1564 1563 1563
(0) (278) (103) (103) (70) (0) (74) (74)
1620 1611 1607 1593 1572 1565 1563 1533
(33) (86) (101) (106) (36) (49) (26) (41)
cHOH
897 891 796 796 724 684 684
(223) (0) (168) (168) (197) (452) (452)
897 874 790 776 725 693 668
(98) (52) (95) (58) (152) (364) (255)
dHOH
519 519
(14) (14)
522 511
(29) (11)
mOAH
370 370
(41) (41)
375 367
(45) (28)
(C24H12)(H2O)8 (367 and 375 cm1). Overall, adsorption on a PAH leads to complex spectroscopic patterns as:
Fig. 8. Finite-temperature BOMD/DFTB spectra in the m1 and m3 regions for (a) (H2O)8 starting from isomer (10 ) and (b) (C24H12)(H2O)8 starting from isomer (1). The 10 K spectra are reproduced in dashed lines to compare with higher T spectra.
As can be seen in Table 1, in bare (H2O)8, the m3 modes are found in the 3756–3925 cm1 range, the m1 modes in the 3291– 3552 cm1 range and the m2 modes in the 1563–1619 cm1 range. When (H2O)8 is adsorbed on C24H12, the ranges are slightly broader i.e. they become 3744–3928 cm1 (m3), 3264–3560 cm1 (m1) and 1533–1620 cm1 (m2). For each mode type, the lower energy bands are redshifted while the higher energy bands are blueshifted with respect to bare (H2O)8. This is due to the fact that the bonds that are equivalent in bare (H2O)8 become inequivalent in (C24H12)(H2O)8, leading to degeneracy quenching. For instance, the higher energy bands in the m3 region are assigned to free H. These modes are degenerate in bare (H2O)8 but are split upon adsorption on the PAH. The bands of higher energy (3920 and 3928 cm1) refer to the free H that are not in interaction with C24H12, and those of lower energy (3870 and 3914 cm1) are assigned to the free H that are in interaction with C24H12. In the lower energy part of the spectrum, the OAH intermolecular bonds stretching modes (mOAH in Table 1) are degenerate in the case of bare (H2O)8 (370 cm1) and are split into two bands in
(i) symmetry decreases, resulting in degeneracy quenching as just been explained, (ii) some modes of (H2O)8 couple with the modes of C24H12. For instance, the dHOH modes perturb the cCH mode, leading to the occurrence of two bands at 824 and 826 cm1 instead of one at 824 cm1 in bare C24H12, (iii) at very low energy, new modes of low intensity assigned to the relative motions of (H2O)8 against C24H12 (translations and rotations) are found at 13, 17, 22, 29, 41 and 49 cm1. The band at 49 cm1 in particular is due to the [(H2O)8]A[C24H12] ‘‘stretching’’ motion. Interestingly, only the mCH mode of C24H12 does not couple to any mode of (H2O)8. As mentioned in our previous paper [30], the IR absolute band positions of the water modes computed with the DFTB method have to be scaled to fit experimental data, the scaling factor depending on the mode. This is expected in part due to the inherent errors of DFTB and also as no low-T anharmonic effect due to nuclear quantum effect likely to be crucial when hydrogen atoms are involved, are not taken into account. Regarding water clusters, scaling factors were obtained for (H2O)3, for which experimental data in cryogenic matrix are available [5], in the m1 (0.97), m2 (1.02), and m3 (0.94) regions. Scaling factors are very commonly used for DFT harmonic spectra when comparison with
60
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
experimental data is achieved [38]. We showed that the scaling factors determined for (H2O)3 could be reasonably applied to larger water clusters, we checked on (H2O)6 for which experimental data are available. Besides, the trends in the shifts of the band positions of (H2O) upon adsorption on C6H6 were shown to be correctly described with the DFTB method against experimental data [30]. The spectral trends induced by the adsorption of water clusters on PAHs are therefore expected to be correctly described with our DFTB method. 3.2.2. Evolution of the IR spectra as a function of T The dynamical IR spectra of bare and adsorbed (H2O)8 at several temperatures, in the m1 and m3 regions, are reported in Fig. 8(a) and (b) respectively. As mentionned previously, the dynamical IR spectra at 10 K are similar to the harmonic spectra, with a maximum discrepancy of 2% on the positions. Indeed, no low-T anharmonic effect due to the quantum nature of nuclei, is taken into account in our approach. In our experience [39,40], the expected evolution of uncoupled vibrational modes when T increases is a redshift of the position, along with a broadening of the band. This is the case for instance of the cCH, mCC and mCH modes of C24H12. We showed that the temperature evolution of the redshift for these modes could be fit by linear laws for C24H12 neutral, cationic and complexed with Si+ or Fe+ [39,40], the slope of the linear fit being the so-called anharmonicity factor v. Obviously, in case of isomerization or melting, such perturbative definition of anharmonicity breaks down. Our computed v value for the cCH mode of bare C24H12 (quoted hereafter vcCH ) was found to be in agreement with experimental data [39]. We showed that the adsorption of Si and Fe on C24 Hþ 12 leads to a decrease of vcCH with respect to bare C24 Hþ 12 [39,40]. We also showed that although the adsorption of one H2O molecule on C24H12 hardly affects vcCH , the adsorption of (H2O)2 leads to an increase of vcCH with respect to bare C24H12, and we proposed that this could be regarded as a signature for edge-coordination [30]. In the case of (H2O)8 (bare and adsorbed), a simple description of the evolution of the band positions in the m1 and m3 regions, and thus the determination of a simple linear law to quantify it, is impossible. This is expected as, due to the high density of bands, modes are likely to couple. Besides, at 150 K for (C24H12)(H2O)8 and 200 K for (H2O)8, isomerisations processes start, while a liquid-like behavior is observed at 200 K and 300 K for (C24H12)(H2O)8 and (H2O)8 respectively. As can be seen in Fig. 8(a) and (b), broadening, merge and blueshifts are observed. From 10 K to 50 K, broadenings are observed. The most energetic m3 band is sligthly redshifted (-2 cm1 for bare (H2O)8). The other bands are blueshifted, the m3 lowest energy band within bare (H2O)8 by + 9 cm1, the m1 highest (resp. lowest) energy band by +5 cm1 (resp. +13 cm1). From 50 K to 100 K, enhanced broadenings and blueshifts occur. The m3 lowest energy band for bare (H2O)8 is indeed shifted by +20 cm1 while the m1 highest (resp. lowest) energy bands is shifted by +11 cm1 (resp. +17 cm1). Interestingly, the most shifted modes are assigned to OH bonds with H atoms involved in intermolecular bonds. The evolution of the IR spectrum of (C24H12)(H2O)8 is similar and the phenomena (broadenings, merge and blueshifts) are enhanced. The spectra at 200 K for bare (H2O)8, and at 150 K for (C24H12)(H2O)8, can be regarded as ‘‘transitory’’ as isomerisation processes start. In the case of bare (H2O)8, two groups of bands can be distinguished for each mode type (m1 and m3), that could be assigned to free OH on the one hand (higher energy m3 band, lower energy m1 band), and to OH with H involved in intermolecular OAH bonds on the other hand (lower energy m3 band, higher energy m1 band). In the case of (C24H12)(H2O)8, the situation is more complicated as other OH bond types are involved (H coordinated to the PAH surface). Finally, at 300 K for bare (H2O)8, and at
200 K for (C24H12)(H2O)8, only one broad m1 band (with maxima at 3523 and 3539 cm1 respectively) and one broad m3 band (with maxima at 3877 and 3883 cm1 respectively) are found: all hydrogen atoms are equivalent, in line with a liquid-phase like regime. When comparing the finite-temperature IR spectra of (C24H12)0,1(H2O)6 published in Ref. [31] with those of (C24H12)0,1(H2O)8 (this work), It may be interesting to mention that, at any temperature, even when the ‘‘melting’’ point is reached, the spectrum of (C24H12)0,1(H2O)6 [31] is clearly different from that of (C24H12)0,1(H2O)8. This is obviously related to the finite-size of the system.
4. Conclusions The present paper provides a detailed study of finite-temperature BOMD/DFTB simulations in the range 10–300 K starting from the lowest energy structure of the complex (H2O)8C24H12. For T 6 50 K, only small amplitude vibrations are observed. At 100 K, diffusion of (H2O)8 on coronene starts, with sporadic isomerisation of (H2O)8. Above 150 K, isomerization processes occur, and at slightly higher T, a melting-like behavior of (H2O)8 is observed. At 250 K, loss of one water molecule is observed. This evolution is quantified by the determination of a chosen set of indicators, and visible on the finite-temperature IR spectra derived from the dynamics. A comparison between the dynamical behavior of bare and adsorbed (H2O)8 is discussed. We show in particular that (H2O)8, when adsorbed on C24H12 isomerizes at a lower T than under its bare form, and that it is consistent with a larger density of structures induced by the adsorption on the PAH. From this work, we infer that the melting temperature for (H2O)8 adsorbed on C24H12 should be in the 150–200 K range and that of bare (H2O)8 higher by 50 K, which is somewhat consistent with the results obtained by Hernandez-Rojas et al. [11] for cationic water systems from parallel tempering Monte Carlo simulations using the rigid molecule approximation. Thanks to the inclusion of all degrees of freedom in the simulations, we also show the influence of the adsorption of (H2O)8 (in its most stable cubic conformation) on the IR spectrum at low T. As expected, the IR spectrum gets more complex as (i) degeneracy is quenched leading to band splittings, (ii) new modes are activated, and (iii) some modes of the water cluster couple to the modes of coronene. The analysis of finite-T IR spectra shows that at all T (prior to the evaporation of the cluster), the spectra of bare and adsorbed (H2O)8 can be distinguished. It is also interesting to quote that, at any temperature, even when the ‘‘melting’’ threshold is reached, the spectrum of (C24H12)0,1(H2O)8 is clearly different from that of (C24H12)0,1(H2O)6, and thus a signature of the structure of the (finite)-size cluster. The present work certainly calls for several extensions. The most important one should concern the quantum character of the vibration which plays a non-negligible role in the case of water clusters [25]. Isomer interconversion will be affected energetically and entropically by quantum delocalization, account of vibrational zero point energy (ZPE) and tunneling effect. Differences in the ZPEs of various isomers will certainly influence the conformational changes. Wang et al. [25] have recently investigated the isomer populations in the isolated water hexamer using Path Integral Molecular Dynamics (PIMD). They found that while in the classical simulations, only a prism-type isomer was populated at low temperature (T = 30 K), decaying into a cage-type isomer for T close to 100 K, and that alternatively in the quantum simulation both isomers already coexist with a 60/40% ratio at T = 30 K. Quantum effects will contribute to a decrease in the melting temperature. Incorporating quantum effects could be targeted by combining the present DFTB code with a PIMD approach. The quantum
A. Simon, F. Spiegelman / Computational and Theoretical Chemistry 1021 (2013) 54–61
character influence on the IR spectra could also certainly be examined using techniques extending the use of the autocorrelation functions of the dipole moments, as done recently on isolated PAH molecules, either following the dynamics of the replica (or polymers), or that of the centroids [41]. Account of quantum vibration is expected to shift and broaden the bands, and possibly change some of the peak intensities. This may be influential in the case of the water octamer, where numerous low energy isomers contribute to the IR spectra even at low temperature. Among other possible relevant theoretical extensions, it would also be interesting to determine the heat capacities for the present system, as was done by Choi [16] for pure water clusters. It would also certainly be relevant to extend the size of the carbonaceous cluster, in particular to examine the size effects (number of vibrations involved) in the coupling between the intramolecular and intermolecular degrees of freedom. In this direction, another development would also be to use thermostatted MD in order to approach the evolution of a surface with a given canonical temperature. Such perspectives would require increased computing resources. References [1] R. Ludwig, Water: from clusters to the bulk, Angew. Chem. Int. Ed. 40 (2001) 1808–1827. [2] C.J. Gruenloh, J.R. Carney, C.A. Arrington, T.S. Zwier, S.Y. Fredericks, K.D. Jordan, Infrared spectrum of a molecular ice cube: the S-4 and D-2d water octamers in benzene–(water)8, Science 276 (1997) 1678–1681. [3] A.D. Kulkarni, R.K. Pathak, L.J. Bartolotti, Structures, energetics, and vibrational spectra of H2O2 (H2O)n, n = 1–6 clusters: Ab initio quantum chemical investigations, J. Phys. Chem. A 109 (2005) 4583–4590. [4] A.D. Kulkarni, S.R. Gadre, S. Nagase, Quantum chemical and electrostatic studies of anionic water clusters H2 O n , J. Mol. Struct. Theochem 851 (2008) 213–219. [5] B. Tremblay, B. Madebène, M. Alikhani, J. Perchard, The vibrational spectrum of the water trimer: Comparison between anharmonic ab initio calculations and neon matrix infrared data between 11,000 and 90 cm1, Chem. Phys. 378 (2010) 27–36. [6] J. Ceponkus, P. Uvdal, B. Nelander, Water tetramer, pentamer, and hexamer in inert matrices, J. Phys. Chem. A 116 (2012) 4842–4850. [7] P. Miro, C.J. Cramer, Water clusters to nanodrops: a tight-binding density functional study, Phys. Chem. Chem. Phys. 15 (2013) 1837–1843. [8] A.D. Kulkarni, K. Babu, S.R. Gadre, L.J. Bartolotti, Exploring hydration patterns of aldehydes and amides: ab initio investigations, J. Phys. Chem. A 108 (2004) 2492–2498. [9] J. Douady, F. Calvo, F. Spiegelman, Effect of an ionic impurity on the caloric curves of water clusters, Eur. Phys. J. D 52 (2009) 47–50. [10] N.J. Singh, M. Park, S.K. Min, S.B. Suh, K.S. Kim, Magic and antimagic protonated water clusters: Exotic structures with unusual dynamic effects, Angew. Chem. Int. Ed. 45 (2006) 3795–3800. [11] J. Hernandez-Rojas, F. Calvo, F. Rabilloud, J. Breton, J.M. Gomez Llorente, Modeling water clusters on cationic carbonaceous seeds, J. Phys. Chem. A 114 (2010) 7267–7274. [12] B. Finlayson-Pitts, J. Pitts (Eds.), Atmospheric Chemistry, Fundamental and Experimental Techniques, WILEY-VCH, Verlag Berlin GmbH, 1986. [13] L.J. Allamandola, A.G.G.M. Tielens, J.R. Barker, Polycyclic aromatichydrocarbons and the unidentified infrared-emission bands – auto exhaust along the milky-way, Astrophys. J. 290 (1985) L25–L28. [14] A. Léger, J.L. Puget, Identification of the unidentified ir emission features of interstellar dust, Astron. Astrophys. 137 (1984) L5–L8. [15] A.G.G.M. Tielens, Interstellar polycyclic aromatic hydrocarbon molecules, An. Rev. Astron. Astrophys. 46 (2008) 289–337. [16] T.H. Choi, Simulation of the (H2O)8 cluster with the SCC–DFTB electronic structure method, Chem. Phys. Lett. 543 (2012) 45–49. [17] J.O. Jensen, P.N. Krishnan, L.A. Burke, Theoretical-study of water clusters – octamer, Chem. Phys. Lett. 246 (1995) 13–19. [18] J.S. Kim, B.J. Mhin, S.J. Lee, K.S. Kim, Entropy-driven structures of the water octamer, Chem. Phys. Lett. 219 (1994) 243–246. [19] G.S. Fanourgakis, E. Apra, W.A. de Jong, S.S. Xantheas, High-level ab initio calculations for the four low-lying families of minima of (H2O)20. II.
[20] [21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39]
[40]
[41]
61
Spectroscopic signatures of the dodecahedron, fused cubes, face-sharing pentagonal prisms, and edge-sharing pentagonal prisms hydrogen bonding networks, J. Chem. Phys. 122 (2005) 134304. T.H. Choi, K.D. Jordan, Application of the SCC–DFTB Method to H+(H2O)6, H+(H2O)21, and H+(H2O)22, J. Phys. Chem. B 114 (2010) 6932–6936. P. Goyal, M. Elstner, Q. Cui, Application of the SCC–DFTB method to neutral and protonated water clusters and bulk water, J. Phys. Chem. B 115 (2011) 6790– 6805. H. Yu, Q. Cui, The vibrational spectra of protonated water clusters: a benchmark for self-consistent-charge density-functional tight binding, J. Chem. Phys. 127 (2007) 234504. S. Maheshwary, N. Patel, N. Sathyamurthy, A.D. Kulkarni, S.R. Gadre, Structure and stability of water clusters (H2O)n, n = 8–20: an ab initio investigation, J. Phys. Chem. A 105 (2001) 10525–10537. S. Karthikeyan, M. Park, I. Shin, K.S. Kim, Structure, stability, thermodynamic properties, and infrared spectra of the protonated water octamer H+(H2O)8, J. Phys. Chem. A 112 (2008) 10120–10124. Y. Wang, V. Babin, J.M. Bowman, F. Paesani, The water hexamer: cage, prism, or both. Full dimensional quantum simulations say both, J. Am. Chem. Soc. 134 (2012) 11116–11119. D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, R. Kaschner, Construction of tight-binding-like potentials on the basis of density-functional theory – application to carbon, Phys. Rev. B 51 (1995) 12947–12957. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, G. Seifert, Self-consistent-charge density-functional tight-binding method for simulations of complex material properties, Phys. Rev. B 58 (1998) 7260–7268. T. Frauenheim, G. Seifert, M. Elsterner, Z. Hajnal, G. Jungnickel, D. Porezag, S. Suhai, R. Scholz, A self-consistent charge density-functional based tightbinding method for predictive materials simulations in physics, chemistry and biology, Phys. Status Solidi (b) 217 (2000) 41–62. A. Oliveira, G. Seifert, T. Heine, H. duarte, Density-functional based tightbinding: an approximate DFT method, J. Braz. Chem. Soc. 20 (2009) 1193– 1205. A. Simon, M. Rapacioli, J. Mascetti, F. Spiegelman, Vibrational spectroscopy and molecular dynamics of water monomers and dimers adsorbed on polycyclic aromatic hydrocarbons, Phys. Chem. Chem. Phys. 14 (2012) 6771–6786. A. Simon, F. Spiegelman, Water clusters adsorbed on polycyclic aromatic hydrocarbons: energetics and conformational dynamics, J. Chem. Phys. 138 (2013) 194309. T. Heine, M. Rapacioli, S. Patchkovskii, J. Frenzel, A. Koster, P. Calaminici, H.A. Duarte, S. Escalante, R. Flores-Moreno, A. Goursot, J. Reveles, D. Salahub, A. Vela, demon-nano experiment (2009).