Ab Initio Molecular Dynamics Studies of Tetrasulfur. Dynamics

12378

J. Phys. Chem. A 2010, 114, 12378–12383

Ab Initio Molecular Dynamics Studies of Tetrasulfur. Dynamics Coupling the C2W Open and D2h Closed Forms of S4 A. Ramı´rez-Solı´s,*,† Franck Jolibois,‡ and Laurent Maron‡ Departamento de Fı´sica, Facultad de Ciencias, UniVersidad Auto´noma del Estado de Morelos, AV. UniVersidad 1001, CuernaVaca, Morelos, 62209, Me´xico, and Laboratoire de Physique de Nano-Objets, INSA, IRSAMC, UniVersite´ Paul Sabatier de Toulouse, 31062 Toulouse Cedex, France ReceiVed: August 23, 2010; ReVised Manuscript ReceiVed: September 21, 2010

Ab initio molecular dynamics (AIMD) simulations were performed on the closed D2h and open C2V isomers of tetrasulfur. After a careful calibration of the electronic structure method, the calculations were done using the BPW91/aug-cc-pVTZ method. This combination of method/basis set adequately reproduces the relative benchmark CCSD(T) energy difference [Matus, M.; Dixon, D.; Peterson, K. A.; Harkless, J. A. W.; Francisco, J. S. J. Chem. Phys. 2007, 127, 174305] between these two isomers and, crucially, the fact that the D2h structure is a transition state linking two equivalent (mirror images) C2V isomers. The trajectories show that the symmetric open C2V isomers interconvert when passing through the D2h closed transition state structure and that, unlike tetraoxygen, no three-dimensional structures arise. The dynamic vibrational analysis yields peaks in good agreement with the static CCSD(T) harmonic frequencies and explains higher peaks as overtones, thus showing that unlike previous AIMD DFT-based approaches, carefully calibrated exchange-correlation functionals can produce reliable molecular dynamics results for complex PESs as the one corresponding to the lowest singlet of S4. I. Introduction From the theoretical point of view, sulfur rings and chains have been widely studied using semiempiric, DFT, and many types of ab initio methods. In particular, for S4, 11 isomers (closed rings and open chains) have been predicted. Until quite recently there was a long-standing debate on whether a “closed” cyclic D2h or the C2V cis-open chain is the lowest energy structure. Many theoretical works have been devoted since the early 1970s to solve this question using a plethora of electronic structure methods.1-25 The sophistication of these studies goes from semiempirical methods, SCF-XR-SW, HF, and DFT calculations using all types of exchange-correlation functionals, to single-reference CISD, MP2, QCISD, CASSCF, restricted multireference-CI (MRCI), coupled cluster up to CCSD(T), and even quantum Monte Carlo methods. Although many of the late studies find either the C2V cis-open chain or the D2h ring structures as the lowest-lying ones, it is interesting to note that no fewer than six different isomers had been predicted to be the lowest energy structure by early 1990. Several semiempirical studies (MNDO, SINDO1) and SCF studies predicted the D2d puckered-ring structure (with four equivalent S-S bonds as in the covalent tetraoxygen molecule) to be the lowest energy isomer.5-7 One of the earliest DFT-based studies9 predicted the planar rectangular D2h isomer to be the lowest-energy structure with the open-chain C2V species lying only 2 kcal/mol above the former. However, the accuracy of this DFT method using the local density approximation was seriously questioned, since the results obtained for the binding energy of S2 was overestimated by more than 25% with respect to the experiment. By the early 1990s there were only three ab initio studies of S4. The first one by Kao10 reported that the lowest singlet isomer * To whom correspondence should be addressed. E-mail: [email protected]. † Universidad Auto´noma del Estado de Morelos. ‡ Universite´ Paul Sabatier de Toulouse.

was the D2d puckered-ring structure, but the lowest energy structure of S4 was predicted to be a helical triplet (S ) 1). However, this study was not reliable since it did not include any correlation effects and, at the SCF level, the stability of triplets relative to singlets is systematically overestimated.18 Another SCF study appeared11 using a slightly larger splitvalence basis set, leading to the same results as the previous work of Kao.10 The first MP2 calculation12 predicted (using SCFoptimized geometries) that the four membered puckered-ring was not the lowest energy isomer. In a more recent ab initio study using slightly larger basis sets,18 the relative stability of the C2V cis-open, trans- and rectangular D2h isomers was addressed with the G2(MP2) method. Because correlation effects beyond second order seem to be crucial for this complex problem, these three species were optimized at the QCISD/6311G(d) level rather than at the standard MP2/6-31G(d) level. At the QCISD/6-311G(d) level, the cis-C2V conformer is predicted to be 5.3 kcal/mol more stable than the D2h rectangular cis-isomer. More importantly, whereas at this level the cis- and the trans- isomers are local minima of the PES, the D2h isomer has two imaginary frequencies (1500i cm-1 and 231i cm-1), which correspond to the two possible S-S out-of-phase stretching combinations (S1-S2 + S3-S4 and S1-S3 + S2-S4), respectively. From the DFT perspective, some newer studies17,19-21 were done to address the energy difference between several S4 isomers using either GGA (BLYP) or hybrid (B3LYP) functionals. In the most detailed DFT-based study, the authors even used the larger aug-cc-pVDZ (AVDZ) basis set and utilized the B3LYP/ AVDZ optimized geometries to perform single-point MP2 and CCSD(T)/AVDZ calculations.20 As an interesting point, note that using the somewhat smaller cc-pVDZ basis set they found that the D2d and C2V structures are minima in the B3LYP/ccpVDZ PES, but the rectangular D2h form is also a transition state at that level of theory. Their B3LYP/aug-cc-pVDZ//

10.1021/jp107994k  2010 American Chemical Society Published on Web 11/04/2010

Molecular Dynamics Studies of Tetrasulfur B3LYP/cc-pVDZ hybrid-functional data agree with those at the GGA-level of Orlova and Goddard (BLYP/6-31G)19 by predicting the C2V (1A1) cis-open form to be the lowest energy structure, followed by the D2h and D2d structures lying 2.1 and 18 kcal/ mol higher in energy, respectively. These newer DFT studies dealing with the C2V versus D2h problem are obviously based on a single Kohn-Sham determinant. Nevertheless, the study by von Niessen16 showed that from the standard MO perspective, this delicate stability problem requires the use of multireference based-methods and the use of highly correlated truly ab initio approaches. This is a key issue that we shall later address, since the accuracy and reliability of single-reference Kohn-Sham methods on this problem is crucial for our present purposes. Recent and more accurate works addressing the S4 stability problem appeared, such as the CCSD(T) single-reference studies by Matus et al.22 and by McCarthy et al.,23 the CASPT2 and limited MRCI study by Wong et al.,24 the excellent review by Bartlett and Musial,25 and the quantum Monte Carlo (QMC) study by Harkless et al.26 The QMC study was done using the single-reference CCSD(T) benchmark optimized geometries of ref 22. Nevertheless, a careful analysis of the entire body of literature on this problem reveals that the most accurate and reliable results from the quantum-theoretical point of view are, up to date, the ones reported by Matus et al.22 Therefore, in the next section we shall make a brief recall of their most important benchmark quality results. The article outline is as follows. In Section II, a detailed account of the calibration procedure is presented along with the computational details. Section III presents the details of the AIMD simulationsand we discuss the results concerning the structural and dynamical aspects of the AIMD simulations in the lowest singlet potential energy surface (PES). Finally, in Section IV, a conclusion and some perspectives are presented. II. Method A. Reliability of the CCSD(T) Results As a Reference. The authors of ref 22 used the CCSD(T) method together with systematic sequences of correlation consistent basis sets extrapolated to the complete basis set limit. With this very accurate single-reference based method, the low-lying D2h structure is a transition state lying only 1.6 kcal/mol above the C2V singlet ground state. They compared the geometry of the CCSD(T) C2V structure using the aV(T+d)Z basis set with that obtained at the CASSCF(24,16) level. At the CCSD(T) level, the short S-S bond distance is 1.915 Å and the long S-S distance is 2.162 Å while the CASSCF gives longer S-S distances of 1.937 and 2.192 Å. As an argument to support the accuracy of the singlereference CCSD(T) results they note that, at the internally contracted MRCI level of theory with a reference function consisting of the most important 50 configurations (selection CI threshold of 0.02) from a preceding CASSCF and including the multireference Davidson correction (MRCI+Q), bond distances of 1.917 and 2.168 Å are found, which are in excellent agreement with the CCSD(T) optimized geometry. In relation to the present study, it is important to highlight that their energetic results obtained with the AVTZ basis set are nearly converged with respect to those obtained with the AVQZ and AV5Z basis sets. In particular note that the AVTZ vibrational frequencies are within 2 cm-1 of those obtained with AVQZ, and the AVTZ S-S distances are within two- to six-hundredths of an angstrom.22 Both facts ensure that the large AVTZ basis set can be considered as very reliable to accurately address this problem.

J. Phys. Chem. A, Vol. 114, No. 47, 2010 12379 Clearly, if an ab initio molecular dynamics study is to be performed on these low-lying tetrasulfur structures, the electronic structure method must satisfy two conditions simultaneously in order to provide reliable results: (a) since hundreds of thousands of nuclear configurations must be sampled, it has be not so expensive in terms of CPU resources to yield energies and gradients in reasonably short times and, (b) it has to provide a rather good description of the potential energy surface (PES), in particular it must yield energy differences and the topological properties of the PES at the stationary points in accord with the most sophisticated description obtained so far, that is, the CCSD(T) study of Matus et al.22 Obviously, the CCSD(T) method cannot be used to perform the AIMD simulations, and only a DFT-based approach could be used. A similar situation was found for the tetraoxygen molecule, which, in spite of having some multireference character,27-29 can actually be well described by AIMD simulations using the Kohn-Sham approach if a carefully calibrated exchange-correlation functional/basis set combination is chosen.30 Therefore, a crucial first step is to find if there exists at least one exchange-correlation functional capable of yielding results in good accordance with the benchmark CCSD(T) ones. This implies a careful calibration of the DFT method in conjunction with the basis set to be used in the AIMD scheme. B. Calibration of the Electronic Structure Method: Choice of the Exchange-Correlation Functional and the Atomic Basis Set. For the reasons given in the Introduction, a proper and well balanced description of the electronic structure of the S4 molecule in its open and closed conformations is a very challenging task for ab initio quantum chemistry. As mentioned above, tetrasulfur has been shown to be a molecule with some multiconfigurational character. Therefore, before embarking on any large-scale quantum-chemical computational task, it is essential to answer the question of whether or not a DFT method, based on a single Kohn-Sham determinant, can be reliably used to perform on-the-fly electronic structure calculations for S4. For this reason we performed a thorough study of both the optimized geometries and the corresponding static vibrational spectra obtained with a set of representative exchangecorrelation functionals. Since it is known that the quality of the basis set is crucial to obtain accurate geometries (in particular the S-S-S angle and the S1-S2-S3-S4 dihedral), we have decided to use the rather large aug-cc-pVTZ (AVTZ) basis set31 for all the DFT and CCSD(T) calculations; the latter are used as the ab initio monodeterminantal references for the calibration of the DFT exchange-correlation (XC) functional. Also, in order to determine which XC functional yields the best overall geometric and topologic description of the ground state PES, these properties were compared with previous ab initio results obtained at the CCSD(T)22 level with the AVTZ basis set. As we shall see, the GGA-type BPW9132 XC functional turned out to be the most adequate to be used in the extensive ab initio molecular dynamics simulations; a discussion of this calibration is presented in Section III. C. Ab Initio Molecular Dynamics (AIMD). The BO ab initio molecular dynamics simulations were carried out with the Geraldyn2.133 code, which has been coupled to the electronic structure modules of the Gaussian98 code.34 The BO-AIMD algorithm in Geraldyn uses the velocity-Verlet integration scheme.35 The simulation times were 10 ps with a time step of 0.5 fs. A Nose´-Hoover chain of thermostats36,37 was used to control the temperature at 298 K. Electronic structure as well as the energy gradient calculations were performed at the BPW9132 level using the large aug-cc-pVTZ (AVTZ) basis set.

12380

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

Both stable structures are of singlet character and the spinunrestricted version of the Kohn-Sham equations was used. The trajectories were simulated starting from the C2V and D2h optimized equilibrium structure taken from 22 without any preferred velocity vector other than the thermal energy. Each 10 ps simulation took 28 days on a QuadCore [email protected] GHz running the Linux versions of Geraldyn2.1-G98. In order to obtain the vibrational spectra, the Fouriertransform of the velocity autocorrelation functions (FT-VACF) were built from the trajectories. The production runs were started following an initial thermalization period, which was achieved after 3 ps of simulation. For readers not acquainted with the theoretical basis of the FT-VACF method, we stress that the vibrational spectrum obtained from the velocity autocorrelation functions do not require the knowledge of the dipole moment; therefore, the calculated vibrational amplitudes yield an accurate description of the relative infrared intensities as well as of the Raman intensities. For this reason, the identification of the lowest vibrational modes relies on the previous static theoretical information. This FT-VACF approach to obtain dynamic vibrational spectra has been successfully used on other small molecules (see for instance, refs 30, 38, and 39). III. Results and Discussion A. General. First of all we should note that 22 years ago an early attempt was made to study the dynamics of tetrasulfur with density functional theory-based molecular dynamics.9 Like all MP2 and MP4 ab initio approaches, that study led to the conclusion that the closed D2h rectangular form is the most stable structure of S4. However, a recent set of Fourier transform microwave spectroscopy experiments23 showed that S4 has a planar trapezoidal shape with C2V symmetry with two similar S-S and one dissimilar S-S bond, consistent with previous matrix isolation results.40,41 The structure can be described as two S2 molecules with triplet ground states strongly coupled to form an overall singlet with two short S-S distances of 1.89 Å and one S-S bond of 2.14 Å. In the most recent microwave study,23 the authors also reported coupled cluster calculations with an approximate correction for the triples CCSD(T). Those calculations showed, in agreement with the CCSD(T) results of ref 22, that the global energy minimum is the C2V structure rather than the D2h form, consistent with the microwave experimental analysis. The conclusion that can be drawn is that the exchange-correlation functional used in the early molecular dynamics study of S4 ref 9 was not accurate enough. This fact leads us inevitably to the need of carefully calibrating the functional/basis set ensemble. B. Calibration of the XC Functional. As mentioned above, the calibration of the XC functional to be used with the large AVTZ basis set is a crucial first step before starting the AIMD simulations. Table 1 shows the optimized geometrical parameters of the closed D2h isomer versus those of the open-cis C2V structure obtained with a variety of XC functionals, along with the energy difference with respect to the most stable isomer (in all cases the C2V one) and a comparison with the corresponding ab initio CCSD(T) benchmark results of ref 22; these results were obtained with the same large AVTZ basis set. It is very interesting to analyze the relative performance of the different types of XC funcionals, that is, local-type like LDA, semilocal GGA-type, and nonlocal hybrid-type functionals. Although most functionals yield rather good SdS double bond distances for both isomers, note the rather disparate quality of the optimized single S-S bond geometries for the closed D2h isomer versus those of the open-cis C2V structure for most of

Ramı´rez-Solı´s et al. TABLE 1: Calibration of the Method: Comparison of DFT Optimized Geometries (Å and Degrees Following Notation of Figure 1) and Relative Energies for the C2W and D2h Singlet Isomers with Benchmark Ab Initio Resultsa C2V(1A1) R12(S-S) R23(SdS) LDA BLYP BP86 B3LYP B3PW91 PBE0 PBE BPW91 CCSD(T)a B3LYPc

2.261 2.284 2.252 2.139 2.117 2.094 2.230 2.235 2.188 2.239

1.899 1.942 1.927 1.924 1.912 1.907 1.923 1.925 1.924 1.939

D2h(1Ag) θ 98.14 103.06 101.83 106.96 106.17 106.34 102.21 102.48 103.50 103.60

∆E R14(S-S) R23(SdS) (kcal/mol) 2.497 2.634 2.573 2.566 2.522 2.499 2.557 2.568 2.571 2.604

1.895 1.933 1.918 1.905 1.894 1.888 1.914 1.915 1.912 1.926

0.2 0.8 0.7 3.3 3.2 3.9 0.8 0.9 1.6b 2.1

a Except where noted, all methods use the large AVTZ basis set. ∆E is the barrier from the C2V minimum to the D2h transition state. a From ref 22 with the AVTZ basis set. b From ref 22 with the CBS(Q5) scheme. c From ref 20, B3LYP/aug-cc-pVDZ//B3LYP/ cc-VDZ values.

the XC functionals. For instance, note that the LDA method largely overestimates the single S-S bond of the C2V isomer while underestimating this bond for the D2h structure. When GGA corrections are included via the BLYP functional, the overestimation of the S-S bond of the C2V isomer gets even worse, while now the long S-S bond of the D2h isomer is also dramatically overestimated. The GGA-type BP86 functional slightly improves the situation with respect to BLYP, but the prediction for the S-S bond of the C2V isomer overshoots by 0.06A the CCSD(T) benchmark value of 2.188 Å. On the other hand, the very popular B3LYP hybrid functional yields an excellent S-S bond of 2.566 Å for the D2h isomer, but underestimates the S-S bond of the C2V isomer by 0.05 Å and the optimized angle is nearly 107°, also overshooting the CCSD(T) value by more than 3°. The hybrid PBE0 functional stands out among all these functionals since it yields much too short distances for all bonds, single S-S or even double SdS bonds, as well as an overshooting angle of 106°. Although the hybrid PBE functional performs better than PBE0, yielding single S-S, double SdS bonds, and C2V angle in reasonably good agreement with the CCSD(T) values, we were somewhat surprised to find out that the BPW91 functional provides the best overall agreement with the benchmark CCSD(T) values. In particular, this GGA-type functional yields the second best possible S-S bond of the C2V isomer (2.235 for BPW91 vs 2.230 Å for PBE) while underestimating the S-S bond of the D2h isomer by only 0.003 Å. We stress that the double SdS bonds of both isomers are reproduced with stunning accuracy (a difference of 0.001 and 0.003 Å with the CCSD(T) values), as well as the S-S-S angle of the C2V isomer. These geometrical results lead to BPW91 as the natural choice, but note that PBE is (nearly) as good a choice in this respect. As found previously,36,37 the DFT descriptions yield the C2V isomer as the lowest energy structure and the D2h isomer as a transition state that links the two equivalent C2V structures. The energy difference between both isomers are also given in Table 1. Note that larger differences with the reference CCSD(T) value (1.6 kcal/mol) are found for the LDA, B3LYP B3PW91, and PBE0 functionals, but the BP86, BLYP, PBE, and BPW91 calculated barriers compare quite well (within 0.7-0.8 kcal/ mol) with that obtained at the CCSD(T)/AVTZ level, BPW91 yielding again the best agreement. Figure 1 presents the optimized geometrical parameters of both isomers with the

Molecular Dynamics Studies of Tetrasulfur

J. Phys. Chem. A, Vol. 114, No. 47, 2010 12381 TABLE 2: Calibration of the Method: Comparison of DFT Vibrational Frequencies (cm-1) for the C2W and D2h Isomers with Previous DFT and ab Initio Resultsa C2V (1A1) isomer

Figure 1. S4 optimized geometries. (a) C2V isomer, (b) D2h isomer. Reference CCSD(T) values in parentheses.

BPW91 functional and its comparison with the reference ab initio CCSD(T)/AVTZ values. Table 2 presents the harmonic frequencies of both structures with a variety of XC functionals and the comparison with previous DFT and ab initio CCSD(T)/AVTZ results. Again note that the BPW91 frequencies for the C2V isomer are in good agreement with those reported at the CCSD(T)/AVTZ level. For the D2h isomer the difference between the BPW91 and the CCSD(T) values is only 4 cm-1 for the smallest (imaginary mode) frequency and only 17 wavenumbers for the largest frequency. The only significant difference is found for the fifth vibrational mode. We stress that all DFT approaches predict this mode above 600 cm-1, while the ab initio CCSD(T) result is 320 cm-1. Note that a previous B3LYP result similar to ours was reported a few years ago (see ref 37); however, no explanation for the difference between the DFT and ab initio frequency for this ag mode of the D2h isomer has been found yet. As a conclusion of the calibration stage, the BPW91 exchange-correlation functional combined with the AVTZ basis set yields the best agreement with the reference CCSD(T) benchmark results of ref 22. C. AIMD Structural and Dynamic Results. Some very interesting insights concerning the dynamical behavior of these isomers can be obtained by analyzing the time evolution of the S-S distances of both isomers. Figures 3 and 4 show the time evolution at 298 K of the S-S bonds for the C2V and D2h structures, respectively. Let us analyze first the behavior of the trajectory starting from the optimized C2V isomer. Note (Figure

ω1

ω2

ω3

ω4

ω5

ω6

method

a1

a2

b2

a1

b2

a1

LDA BP86 B3LYP B3PW91 PBE0 PBE BLYP BPW91 CCSD(T)a CASSCFb B3LYP/6-31G(2df)c

84.0 88.1 105.0 110.8 117.4 91.7 78.2 88.6 118.8 110.8 104

240.0 214.2 206.4 215.3 219.7 216.0 200.7 212.8 210.5 181.1 207

314.7 296.1 326.9 335.4 341.1 308.5 272.4 302.8 325.3 287.3 329

336.1 317.6 358.4 381.2 406.0 320.4 304.8 318.4 328.7 313.5 373

674.9 640.7 640.4 663.2 673.5 648.4 615.6 644.9 640.2 627.9 649

713.3 671.3 666.1 690.9 702.6 678.4 643.6 674.5 681.2 642.3 674

D2h(1Ag) isomer

LDA BP86 B3LYP B3PW91 PBE0 PBE BLYP BPW91 CCSD(T)a B3LYP/6-31G(2df)c

ω1

ω2

ω3

ω4

ω5

ω6

ag

au

b3g

b2u

ag

b1u

58.0i 82.1i 137.3i 143.5i 158.6i 87.3i 83.4i 87.4i 91.2i 145 i

251.5 234.5 242.3 251.8 258.5 238.0 221.9 235.5 233.4 245

307.6 274.7 283.4 301.0 312.1 279.7 251.2 275.7 233.5 284

338.5 321.0 331.2 341.9 348.9 324.5 307.3 322.5 278.1 334

680.5 650.2 662.5 684.1 696.2 658.4 625.5 654.7 319.6 671

727.1 695.1 717.2 740.4 755.7 703.8 668.7 699.9 716.7 724

a All methods use the large AVTZ basis set (see text). a From ref 22. b From ref 36. c From ref 37.

3) that two distances remain small (around 1.94A) during the whole simulation up to 10 ps, and these correspond to the pair of stable double bonds mentioned in Section A. These distances oscillate with small amplitudes (of around 0.08 Å) at a rather high frequency of one cycle every 45 fs, thus revealing the nature of the double bonds. On the other hand, one readily recognizes that the two larger distances oscillate differently in three regimes (0-3.3 ps, 3.3-5.6 ps, and 5.6-10 ps) with much larger amplitudes (maximum amplitude of 1.4 Å). Note that the S-S distance which starts at 2.25 Å describes the single

Figure 2. Temporal evolution of the S-S distances (Å) for the AIMD simulation starting from the C2V isomer. Horizontal axis in femtoseconds.

12382

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

Ramı´rez-Solı´s et al.

Figure 3. Temporal evolution of the S-S distances (Å) for the AIMD simulation starting from the D2h isomer. Horizontal axis in femtoseconds.

Figure 4. BO-AIMD vibrational spectrum of the C2V isomer at 298 K. Vertical axis in arbitrary units, horizontal axis in wavenumbers.

bond bridging the S2 pairs and remains approximately constant from the start of the simulation up to 3.3 ps, when this single bond gets broken and the other long distance (3.05 Å, which initially is related to the nonbonded S-S pair) gets shorter and becomes that of the other singly bonded pair. This transition of the initial C2V structure to the other takes around 450 fs and the new C2V structure lasts 2.3 ps, when another transition occurs, bringing back the molecule to its original conformation 5.6 ps after the start of the simulation. Clearly, only much longer simulation times could reveal a statistically reliable lifetime for one of the two equivalent C2V isomers. A completely different situation arises when the simulation starts from the D2h isomer. In Figure 4 it can be seen that the two shortest distances (corresponding to a couple of double S-S bonds) oscillate with small amplitudes throughout the simulation. However, only 250 fs after the start of the simulation the molecule acquires C2V symmetry, since one of the two initially identical longer distances (both starting at 2.5 Å) becomes much longer (around 3.6 Å) while the other shortens to around 2.2 Å. This situation reverses many times, revealing a molecule that goes from one C2V form to the other with a period of around 400 fs until 9 ps. It is quite interesting to find that the molecule changes its dynamical regime after only 20 C2V f D2h f C2V transitions and gets locked into one of the C2V isomers. After 9 ps the dynamical behavior of the molecule becomes similar (oscillation amplitudes for the short and long S-S bonds as well as the oscillation periods) to that found of the

Figure 5. BO-AIMD vibrational spectrum of the D2h isomer at 298 K. Vertical axis in arbitrary units, horizontal axis in wavenumbers.

simulation that started from the C2V isomer. We note that both structures remain mostly planar throughout the 10 ps simulations. D. AIMD Vibrational Spectra and Nonharmonic Effects. After full thermalization was achieved (around 3 ps after the start of the simulations), the FT-VACF vibrational spectra at 298 K are given in Figures 4 (starting from the C2V isomer) and 5 (starting from the D2h isomer) up to 2000 wavenumbers since, after that cutoff energy, statistical noise dominates the spectra. In continuous black line we show the Gaussian-convoluted curve, and we have marked with the letters (a, b, ...., i) the peaks that have been identified as arising from the fundamental modes. The nonharmonic (NH) effects were determined according to the usual definition (see, for instance, ref 39). Thus, for the six fundamental modes of the stable structure the NH were determined by the difference of the dynamically derived frequency ωAIMD(ν1, ν2, ν3, ν4, ν5, ν6) minus the corresponding static purely harmonic BPW91/AVTZ frequencies (ω1, ω2, ..., ω6) given in Table 2, that is, NH(ν1, ν2, ν3, ν4, ν5, ν6) ) ωAIMD(ν1, ν2, ν3, ν4, ν5, ν6) ν1ω1 - ν2ω2 - ν3ω3 - ν4ω4 - ν5ω5 - ν6ω6 The vibrational peaks appear in Table 3 for both isomers with the corresponding NH effects for the C2V structure. It can be seen that the NH effects are small and negative for the a1 and b2 modes, whereas they are positive and slightly larger for the remaining vibrational modes. It is interesting to note that the g and h peaks might either arise as a combination of overtones of the first four fundamental modes or as simple overtones of the fifth and sixth modes, respectively.

Molecular Dynamics Studies of Tetrasulfur

J. Phys. Chem. A, Vol. 114, No. 47, 2010 12383

TABLE 3: Ab Initio Molecular Dynamics Vibrational Peaks Obtained When Starting Geometries Are the C2W and D2h Singlet Isomers of S4a BO-AIMD C2V isomer a b c d e f g h

a b c d e

BPW91/AVTZ NHb

ν1

ν2

ν3

ν4

ν5

ν6

1 0 0 0 0 0 1 0 1 0

0 1 0 0 0 0 1 0 1 0

0 0 1 0 0 0 2 0 1 0

0 0 0 1 0 0 1 0 2 0

0 0 0 0 1 0 0 2 0 0

0 0 0 0 0 1 0 0 0 2

82 215 317 331 664 687 1260

0 0 0 0 0

0 0 0 0 0

0 1 0 0 0

D2h isomer 0 0 0 0 1 0 0 1 0 0

0 0 0 0 1

206 291 319 636 690

-6 2 -14 13 19 13

1293

a

Vibrational mode assignments for the g and h peaks are given with tentative harmonic quantum numbers (see text). Energies in wavenumbers. b For definition of the nonharmonicities see text.

IV. Conclusions Ab initio molecular dynamics simulations at room temperature were performed on the closed D2h and open C2V isomers of tetrasulfur, which have been predicted to be lowest energy ones. A careful calibration of the DFT electronic structure method and basis set was mandatory in order to ensure that it adequately reproduces the relative benchmark CCSD(T) energy difference between these two isomers and the fact the D2h structure is a transition state linking two equivalent C2V isomers. Thus, the Kohn-Sham calculations were done using the BPW91 GGA-type exchange-corelation functional with the large AVTZ basis The AIMD derived vibrational spectrum of the C2V isomer yields peaks in good agreement with the static CCSD(T) harmonic frequencies. The simulation starting from the C2V isomer clearly reveals the transition of the initial C2V structure to ist mirror image after 450 fs and the new C2V structure lasts 2.3 ps, when another transition occurs bringing back the molecule to its original conformation 5.6 ps after the start of the simulation. Clearly, only much longer simulation times could reveal a statistically reliable lifetime for one of the two equivalent C2V isomers. The simulation starting from the D2h transition state shows a completely different behavior, where after only 20 C2V f D2h f C2V transitions the molecule gets locked into one of the C2V isomers; after 9 ps the dynamical behavior becomes similar to that found of the simulation which started from the C2V isomer. The nonharmonic effects were found to be a small percent of the vibrational energies of all modes for the C2V structure. As an important corollary, we stress that unlike previous AIMD DFTbased approaches,9 carefully calibrated exchange-correlation functionals can actually produce reliable molecular dynamics results for complex PES as the one corresponding to the lowest singlet of S4. Acknowledgment. A. R. S. thanks support from CONACYT project number 45986-E and SESIC-SEP/FOMES2000 project “Co´mputo cientı´fico” for unlimited CPU time on the IBM-p690 supercomputer at UAEM. References and Notes (1) Meyer, B.; Spitzer, K. J. Phys. Chem. 1972, 76, 2214. (2) Kao, J.; Allinger, N. L. Inorg. Chem. 1977, 16, 35.

(3) (a) Foti, A. E.; Smith, V. H., Jr.; Salahub, D. R. Chem. Phys. Lett. 1978, 57, 33. (b) Salahub, D. R.; Foti, A. E.; Smith, V. H., Jr. J. Am. Chem. Soc. 1978, 100, 7847. (4) Dewar, M. J.; McKee, M. L. J. Comp. Chem. 1983, 4, 84. (5) Baird, N. C. J. Comp. Chem. 1984, 5, 35. (6) Bolotin, A. B.; Pipiraite, P. P.; Ruzhene, A. Y.; Gorbunov, B. N.; Tarkhov, G. V.; Simanenkova, L. B. J. Siruci. Chem. 1988, 29, 182. (7) (a) Jug, K.; Iffert, R. J. Comput. Chem. 1987, 8, 1004 and references therein. (b) Jug, K.; Iffert, R. J. Mol. Struct. (Theochem) 1989, 186, 347. (8) (a) Stillinger, F. H.; Weber, T. A.; LaViolette, R. A. J. Chem. Phys. 1986, 85, 6460. (b) Stillinger, F. H.; Weber, T. A. J. Chem. Phys. 1987, 91, 4899. (9) Hohl, D.; Jones, R. O.; Car, R.; Parrinello, M. J. Chem. Phys. 1988, 89, 6823. (10) Kao, J. Inorg. Chem. 1977, 16, 2085. (11) Feng, W. L.; Novaro, O. Int. J. Quantum Chem. 1984, 26, 521. (12) (a) Dixon, D. A.; Wasserman, E. J. Phys. Chem. 1990, 94, 5772. (b) McLean, A. D.; Chandler, G. S. J. Chem. Phys. 1980, 72, 5639. (13) Quelch, G. E.; Schaeffer, H. F., III; Marsden, C. J. J. Am. Chem. Soc. 1990, 112, 8719. (14) Knowles, P. J.; Werner, H.-J. Chem. Phys. Lett. 1985, 115, 259. Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1985, 82, 5053. (15) Ragavachari, K.; Rohlfing, C. M.; Binkley, J. S. J. Chem. Phys. 1990, 93, 5862. (16) von Niessen, W. J. Chem. Phys. 1991, 95, 8301. (17) Hunsicker, S.; Jones, R. O.; Gantefo¨r, G. J. Chem. Phys. 1995, 102, 5917. (18) Abboud, J. L. M.; Esseffar, M.; Herreros, M.; Mo´, O.; Molina, M. T.; Notario, R.; Ya´n˜ez, M. J. Phys. Chem. A 1998, 102, 7996. (19) Orlova, G.; Goddard, J. D. J. Phys. Chem. A 1999, 103, 6825. (20) Millefiori, S.; Alparone, A. J. Phys. Chem. A 2001, 105, 9489. (21) Jones, R. O.; Ballone, P. J. Chem. Phys. 2003, 118, 9257 and references therein. (22) Matus, M.; Dixon, D.; Peterson, K. A.; Harkless, J. A. W.; Francisco, J. S. J. Chem. Phys. 2007, 127, 174305. (23) McCarthy, M.; Thorwirth, C. S.; Gottlieb, C. A.; Thaddeus, P. J. Am. Chem. Soc. 2004, 126, 4096. McCarthy, M. C.; Thorwirth, S.; Gottlieb, C. A.; Thaddeus, P. J. Phys. Chem. 2004, 121, 632. Thorwirth, S.; McCarthy, M. C.; Gottlieb, C. A.; Thaddeus, P.; Gupta, H.; Stanton, J. F. J. Chem. Phys. 2005, 123, 054326. (24) Wong, M. W.; Steudel, R. Chem. Phys. Lett. 2003, 379, 162. (25) Bartlett, R. J.; Musial, M. ReV. Mod. Phys. 2007, 79, 291. (26) Harkless, J. A. W.; Francisco, J. J. Phys. Chem A 2008, 112, 2088. (27) Herna´ndez-Lamoneda, R.; Ramı´rez-Solı´s, A. J. Chem. Phys. 2000, 113, 4139. Herna´ndez-Lamoneda, R.; Ramı´rez-Solı´s, A. Chem. Phys. Lett. 2000, 321, 191. (28) Herna´ndez-Lamoneda, R.; Ramı´rez-Solı´s, A. J. Chem. Phys. 2004, 120, 10084. (29) Caffarel, M.; Scemama, A.; Herna´ndez-Lamoneda, R.; Ramı´rezSolı´s, A. Phys. ReV. Lett. 2007, 99, 153001. (30) Ramı´rez-Solı´s, A.; Jolibois, F.; Maron, L. Chem. Phys. Lett. 2010, 485, 16. (31) Peterson, K. A.; Wilson, A. K.; Woon, D. E.; Dunning, T. H. Theor. Chem. Acc. 1997, 97, 251. (32) Becke, A. D. Phys. ReV. A 1988, 38, 3098. Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244. Burke, K.; Perdew, J. P.; Yang, W. Electronic Density Functional Theory: Recent Progress and New Directions, Dobson, J. F. eds; Vignale, G.; Das, M. P. 1998. (33) Raynaud, C.; Maron, L.; Daudey, J. P.; Jolibois, F. Phys. Chem. Chem. Phys. 2004, 6, 4226. (34) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Rega, N.; Salvador, P.; Dannenberg, J. J.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andre´s, J. L.; Gonza´lez, C.; Head-Gordon, M.; Replogle, E. S. Pople, J. A., Gaussian 98 A.11.3; Pittsburg,USA, 2002. (35) Verlet, L. Phys. ReV. 1967, 159, 98. (36) Nose´, S. J. Chem. Phys. 1981, 81, 511. (37) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (38) Maron, L.; Ramı´rez-Solı´s, A. J. Phys. Chem. A 2007, 111, 3173. (39) Jolibois, F.; Maron, L.; Ramı´rez-Solı´s, A. J. Molec. Struct. THEOCHEM 2009, 899, 9. (40) Brabson, D.; Mielke, Z.; Andrews, L. J. Phys. Chem. 1991, 95, 79. (41) Hassazadeh, P.; Andrews, L. J. Phys. Chem. 1992, 96, 6579.

JP107994K