Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Chem Phys Lett. Author manuscript; available in PMC 2010 November 8.
Published in final edited form as:
Chem Phys Lett. 2007 June 25; 441(4-6): 187–193.
doi:  10.1016/j.cplett.2007.05.026
PMCID: PMC2975466

A comparative electron correlation treatment in H2S-benzene dimer with DFT and wavefunction-based ab initio methods


Three major conformations of H2S-benzene dimer have been located with a variety of density functional theories (DFT) and second order Møller-Plesset perturbation (MP2). In line with an experiment, MP2 results indicate that the tilted Cs-symmetry structure is a stable dimer, yet a C2v-symmetry structure is only a second-order saddle point. Although all of the examined DFT methods also predict the binding between H2S and benzene as MP2 and the coupled cluster method with single and double excitations and perturbative triples (CCSD(T)) do, they considerably underestimate binding energies as compared with CCSD(T) results. However, PW91LYP and MPWB1K reproduce the binding sequence obtained with MP2 for the dimers, and provide the best binding energies among the tested DFT methods. The method of increments with the orbitals of H2S and π orbitals of the benzene recovers 99% of the total binding from the full CCSD(T).

1. Introduction

Because of the high sensitivity of naked single-walled carbon nanotubes (SWCNTs) and the high selectivity of polymer-coated SWCNTs to gaseous molecules such as O2, NH3, NO2, etc., SWCNT-based chemical sensors, as a new generation of promising sensor at room temperature, are gaining considerable interest from both theories and experiments [1,2]. Many natural resources including petroleum, natural gas, and coal mines contain or generate hydrogen sulphide (H2S) that can cause serious effects on the nervous system at low concentrations and even may be fatal at high concentrations. Thus, in situ monitoring of H2S in technologies related to extraction or processing of these materials is also of great importance. However, to the best of our knowledge, there is no report on whether SWCNT has a significant response to hydrogen sulphide (H2S).

It is speculated that there may be only weak non-covalent interactions between H2S and SWCNTs. An accurate description of such weak interactions requires high level quantum mechanics methods for electron correlation, like MP2 and even CCSD(T). However, these methods are highly demanding in computer resources, and can be hardly applied directly to H2S-SWCNT system. As an alternative way to treat electron correlation, DFT is considerably less demanding on computational resources. Searching for a proper DFT theory to accurately describe weak interactions has been a challenge because of the difficulty to represent Pauli repulsion for such systems.[22,23] The performance of the newly developed GGA exchange functionals such as PW91 [3,4] and PBE [5] and their modifications (mPW1[6] and PBE1[7]) for the weak interaction systems has been significantly improved. Truhlar et al recently recommended several methods like PWB6K, PW6B95, MPW1B95 and MPWB1K for hydrogen bond systems such as H2O-C6H6, NH3-C6H6 systems [8] and for ••• stacking systems [9]. The accuracy of such DFT functionals to S-H••• system will be investigated.

Using the scheme of increments, in the recent years CCSD(T) has been successfully applied to extended systems like solids and polymers.[10,11] The method relies on localized orbitals or groups of localized orbitals generated in a Hartree–Fock reference calculation. The electron correlation method recently was applied to C60[12], and shows a very promising method to accurately treat electron correlation for the extended systems such as nanotubes. Whether the method of increments could be properly applied to S-H•••π system will be tested in this letter.

Dykstra et al. measured rotational spectrum of the weakly bonded C6H6-H2S dimer, and consequently proposed that H2S almost locates on the C6 axis of the C6H6 ring with a distance of 3.818 Å between S and the center of mass (c.m.) of C6H6 [13]. The angle between the C6 axis and C2v axis of H2S was determined to be 28.5°[13]. Sherrill et al has recently done the highest level calculations so far for the dimer, including a potential energy curve with CCSD(T)/aug-cc-pVDZ and ab initio limit for the C2v symmetry complex[14].

In the present work, conformations of H2S-benzene dimer are fully optimized with MP2 methods. In order to validate a relatively proper DFT for H2S-SWCNT system, a variety of DFT methods are then applied to H2S-C6H6 dimers. Higher level electron correlation treatments, full electron correlation with CCSD(T) and the method of increments have also been done to calibrate DFT/MP2 results as well as to discuss the extension of the method of increments to the extended systems of S-H•••π interaction.


The structures of the H2S-benzene dimer are firstly fully optimized with a variety of GGA and hybrid GGA types of DFT, as well as MP2 (FC) at different basis sets such as 6-311++G(d,p) and Dunning’s aug-cc-pVXZ basis sets. The binding energy is defined as the energy difference between the dimer and two individuals. We correct the binding energies for basis set superposition error (BSSE) using the counterpoise approach.[20,21] Atomic net charges were calculated using the Mulliken, CHELPG,[15] and NPA schemes[16]. All of the calculations are performed using GAUSSIAN03[17].

To calibrate the DFT and MP2 levels, energy correlation is performed with CCSD(T) as implemented in MOLPRO [18]. The method of increments expands the correlation energy of an extended system in terms of localized orbitals or localized orbital groups numbered i, j, k.[10,11] For the H2S-benzene dimer, the groups would be the H2S molecule and the different kinds of σ and π-bonds in the benzene molecule.


Since the correlation in the dimer and two entities are almost the same for the weakly bound system, the one-body term εi has little influence on the binding. The important terms for the binding are the non-additive contributions to the correlation energy with the leading two-body increments Δεij = εij − (εi − εj) Therefore, after a Hartree-Fock calculation of the whole system and a localization of the occupied orbitals, e.g. with the Foster-Boys procedure [19], one can directly calculate the binding energy of the system. Specifically, we have localized the occupied Hartree-Fock orbitals of H2S-benzene at aug-cc-pVDZ level and selected the orbitals of the H2S as one group, the 3 orbitals of the π system of the benzene as the second group and the sigma-bonds of the benzene as the third group. The two-body increments, by including the orbitals of the H2S and π orbitals of the benzene correspond to the full CCSD(T) with H2S and benzene π system correlated. An additional advantage of the incremental method is that one has not to consider the BSSE correction, because the energies of the separate systems may not be determined. For the incremental and full CCSD(T) calculations the program package MOLPRO[18] was employed.

3. Results and Discussion

3.1 Binding energy and structures of the H2S-benzene dimer with MP2 method

Three bounded H2S-benzene dimers were located at MP2 levels. As shown in Figure 1, 1a has C2v symmetry; while other two conformations (1b and 1c) have Cs symmetry. 1a was characterized to be a second-order saddle point. Two imaginary frequencies, assigning to S-H swings in and out of H2S molecular plane for 1a at MP2/aug-cc-pVVDZ level are 92i and 76i cm−1, while the six lowest frequencies are almost zero (<2cm−1). Binding energies in Table 1 for 1a at MP2/aug-cc-pVXZ(X=D,T and Q) are rather similar to those reported by Sherrill et al.[14] Evidenced by all real frequencies for structures 1b and 1c at MP2/aug-cc-pVXZ(X=D and T), they could be considered as local minima. Structure 1b agrees well with the experimentally identified dimer structure by rotational spectrum.[13] The reported geometric parameters SX (distance between S and the center of mass of C6H6) and angle A (the angle between C6 axis of benzene and C2v axis of H2S) are 3.818Å and 28.5°; the present results at MP2/aug-cc-pVDZ shown in Table 2 are 3.625Å and 27°.

Figure 1
Three conformations for the H2S-benzene dimer (a) C2v symmetry (b) Cs symmetry with one hydrogen pointing to C in a SHC angle of approximately 180° (c) Cs symmetry with one H pointing to the center of benzene C6 ring.
Table 1
Binding energies (BE, kcal/mol) of the benzene-H2S dimers calculated with DFT, MP2, and CC methods.
Table 2
Leading structure parameters (S-X/Å: distance from S to the center of mass of C6H6; S-H-C/°: angle between donated S-H and the nearest C of C6H6; A/°: the angle between C2v axis of H2S and C6 axis of C6H6), dipole moments (D/Debye) ...

At all MP2 levels, structures 1b and 1c have higher negative binding energies than 1a by approximately 0.06–0.15kcal/mol, and 1c also displays slightly more stable than 1b by 0.02 (MP2/aug-cc-pVDZ) and 0.07 kcal/mol (MP2/aug-cc-pVTZ). The binding energy differences indicate that the potential energy surface of the dimer is quite shallow. The CHELPG charge analysis listed in Table 2 indicates that small amount of charge is transferred from benzene to H2S for both 1a and 1b. The charge carried by H2S is −0.037e and −0.047e at MP2/aug-cc-pVDZ level for 1a and 1b, respectively; while the amount of charge is decreased with NPA analysis to −0.006e for 1b.

3.2 Results from density functional methods

One of the current objectives is to identify a reasonable DFT method suited to describe S-H(...)π interaction for large systems. Comparing DFT results in Table 1 with those of MP2 and CCSD(T), although all of the examined DFT methods also predict the binding between H2S and benzene, they significantly underestimate the binding of H2S and benzene. However, they qualitatively predict the similar trend to that from MP2, especially C2v symmetry structure 1a being a second-order saddle point rather than a local minimum and having less negative binding energy than structures 1b and 1c. Among the DFT methods, according to Table 1 and Figure 2 GGA PW91LYP and hybrid DFT MPWB1K give the closest binding energies to those from the more sophisticated methods. At PW91LYP/6-311++G(d,p) level, the binding energies for structures 1a, 1b and 1c are −1.83 (−2.36, the data in the parentheses for without BSSE), −1.87 (−2.37) and −1.84 (−2.24) kcal/mol, respectively. Because of convergence problem of MPWB1K and MPW1B95/6-311++G(d,p) in the present system, only aug-cc-pVXZ basis sets were applied to MPWB1K and MPW1B95 methods. The binding energies yielded by MPWB1K/aug-cc-pVDZ for structures 1a, 1b and 1c are −1.76 (−2.14, without BSSE), (−2.37) and (−2.38) kcal/mol, respectively. BSSE calculations again encounter convergence problems for MPWB1K and MPW1B95 methods even with aug-cc-pVD(T)Z basis sets. The only BSSE corrected binding energy for 1a structure at MPWB1K/aug-cc-pVDZ, is slightly more negative by 0.14 kcal/mol than that from PW91LYP/aug-cc-pVDZ (−1.76 vs. −1.62 kcal/mol), which indicates that MPWB1K and PW91LYP have similar performance on S-H•••π interaction. The other DFT methods, like PW91PW91, PBE1PBE and PBEPBE, provide much lower negative binding energies than PW91LYP and MPWB1K (−1.47, −1.22 and −1.21 kcal/mol vs. −1.62 and −1.76 kcal/mol for C2v symmetry complex 1a). B3PW91 and B3LYP using conventional hybrid exchange functional B3 generate surprisingly low binding energies. Taking the binding energy of C2V symmetry structure dimer at CCSD(T)/CBS as a reference,[14] the scaling factor for the binding energy at PW91LYP/6-311++G(d,p) level is about 1.54.

Figure 2
The binding energies for H2S-C6H6 dimers as predicted by DFT methods and wave-function based methods.

According to Table 1, PW91LYP demonstrates a weak basis set dependence. The BSSE corrected binding energies predicted by PW91LYP/aug-cc-pVDZ and PW91LYP/aug-cc-pVTZ are already almost the same, and they are only lower than those from PW91LYP/6-311++G(d,p) by approximately 0.20 (1a and 1b) and 0.15 kcal/mol (1c). Other tested DFT methods such as PW91PW91, PBEPBE and PBE1PBE have similar basis set effects. A comparison for the major geometric parameters from PW91LYP, MPWB1K and MP2 is also made in Table 2. The distances (S-X) predicted from MPWB1K/aug-cc-pVDZ (3.91 and 3.83 Å for 1a and 1b) are longer than those from MP2/aug-cc-pVDZ (3.62 and 3.63 Å) by 0.2–0.3 Å, yet shorter than those with PW91LYP/aug-cc-pVDZ (4.03 and 4.15 Å). This is consistent with the above fact that the binding energies from MPWB1K are more negative than those from PW91LYP by approximately 0.1kcal/mol.

To test the accuracy of the above DFT calculations, the integration grid was increased from fine (the default grid in Gaussian 03, used throughout the paper) to ultrafine for PW91LYP/aug-cc-pVDZ. As shown in Table 1, the BSSE corrected binding energies for 1a, 1b and 1c with the ultrafine grid are −1.64, −1.69 and −1.69 kcal/mol, respectively. The binding energy for 1a is only slightly different from that yielded by the fine grid (−1.64 vs. −1.62 kcal/mol); while for both 1b and 1c they are the same as those from the above fine grid calculations. This shows that the default grid in Gaussian 03 is good enough for the current DFT calculations.

3.3 Full electron correlation with CCSD and CCSD(T)

With two different basis sets aug-cc-pV(D,T)Z, CCSD and CCSD(T) calculations are performed on 1a and 1b optimized with PW91LYP/aug-cc-pVDZ by fixing symmetry and only varying SX. As shown in Figure 3, the present binding energy for C2v-symmetry structure (1a), −2.73 kcal/mol at aug-cc-pVTZ level, is more negative than that (−2.64 kcal/mol) obtained by Sherrill et al.,[14] where rigid monomer geometries for H2S and C6H6 were used. The current higher negative binding energy is probably due to the relaxation of monomers, which were taken into account at the PW91LYP level. Upon the formation of H2S-benzene dimer, C-C bond lengths of benzene are only slightly stretched by approximately 0.001Å; however, as shown in Table 2 S-H bonds stretch by 0.001 Å in 1a and 0.003 Å in 1b, and HSH angle is compressed by as much as 1.0° Although PW91LYP and other DFT methods yield slightly higher negative binding energy for the Cs symmetry structure than for C2v structure (1a), the restricted CCSD and CCSD(T) levels with both basis sets provide the opposite results. To get the experimental results at the CC level a full geometry optimization therefore is necessary.

Figure 3
The potential energy curves of H2S-C6H6 dimers (1a and 1b) predicted by CC methods at aug-cc-pVDZ (avdz) and aug-cc-pVTZ (avtz).

3.4. Application of the method of increments

Although it is not necessary for such small dimers to apply the method of increments, because a full CCSD(T) treatment is feasible, it is a good example to validate the method of increments for S-H•••π interaction system. For the C2v-symmetry structure, we have localized the occupied Hartree-Fock orbitals of the systems at aug-cc-pVDZ level and selected the orbitals of the H2S as one group, the 3 orbitals of the π system of the benzene as another group and the sigma-bonds of the benzene as the third group. As shown in Figure 4, the binding energy by only including the second-order increment, i.e., the first two correlated groups, is rather comparable to the overall binding energy, and this contribution yields 99% of the total binding, which supports the conclusion from Tauer et al that the attraction in H2S-benzene mainly come from the partially positive hydrogens in H2S and the negatively charged π-cloud of the benzene [14]. Moreover, the optimized distances from S to the mass center of the benzene with the two-body increment and full CCSD(T) treatment are almost identical (3.88Å).

Figure 4
The potential energy curves of H2S-C6H6 C2v-symmetry dimer (1a) predicted by a full CC calculation and the method of increments, both with a basis set of aug-cc-pVDZ (avdz).

4. Conclusions

H2S-benzene dimer was investigated with a variety of DFT and wave-function based quantum mechanics methods such as MP2, CCSD(T) and the method of increments. In line with experimental findings in the rotational spectrum, MP2 and PW91LYP predicted that a Cs-symmetry structure (an angle of approximately 30° between C2v axis of H2S and C6 axis of C6H6) is a stable dimer, yet the previously reported C2v-symmetry structure is a second-order saddle point instead of a minimum. PW91LYP and MPWB1K reproduce the MP2 sequence of binding strength for the three located dimer conformations, and provide the best binding energies among the tested DFT methods. In agreement with Tauer et al,[14] the method of increments with the orbitals of H2S and π orbitals of the benzene recovers 99% of the total binding from the full CCSD(T) calculation.


The work was funded by a Faculty Development Award (Y. Wang) from Albany State University through grants P20 MD001085 and G11 HD049635 from the National Institutes of Health. Y. Wang greatly appreciates this support. Visiting MPIPKS as a guest scientist (Y. Wang) in the summer of 2006 is also acknowledged.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


1. Kong J, Franklin NR, Zhou C, Chapline MG, Peng S, Cho K, Cho K, Dai H. Science. 2000;287:622. [PubMed]
2. Collins PG, Bradley K, Ishigami M, Zettl A. Science. 2000;287:1801. [PubMed]
3. Burke K, Perdew JP, Wang Y. In: Electronic Density Functional Theory: Recent Progress and New Directions. Dobson JF, Vignale G, Das MP, editors. Plenum: 1998.
4. Perdew JP, Burke K, Wang Y. Phys. Rev. B. 1996;54:16533. [PubMed]
5. Perdew JP, Burke K, Ernzerhof M. Phys. Rev. Lett. 1996;77:3865. [PubMed]
6. Adamo C, Barone V. J. Chem. Phys. 1998;108:664.
7. Adamo C, Barone V. J. Chem. Phys. 2002;116:5933.
8. Zhao Y, Tishchenko O, Truhlar DG. J. Phys. Chem. B. 2005;109:19046. [PubMed]
9. Zhao Y, Truhlar DG. J. Phys. Chem. A. 2005;109:4209. [PubMed]
10. Stoll H. Phys. Rev. B. 1992;46:6700. [PubMed]
11. Paulus B. Phys. Rep. 2006;428:1.
12. Paulus B. Int. J. Quant. Chem. 2004;100:1026.
13. Arunan E, Emilsson T, Gutowsky HS, Fraser GT, Oliveira GD, Dykstra CE. J. Chem. Phys. 2002;117:9766.
14. Tauer TP, Derrick ME, Sherrill CD. J. Phys. Chem. A. 2005;109:191. [PubMed]
15. Breneman CM, Wiberg KB. J. Comp. Chem. 1990;11:361.
16. Reed AE, Curtiss LA, Weinhold F. Chem. Rev. 1988;88:899.
17. Frisch MJ, et al. Gaussian 03, Revision A. D02. Wallingford CT: Gaussian, Inc.; 2004.
18. Werner H-J, Knowles PJ, Lindh R, Schütz M, et al. MOLPRO, Version 2002.6, A Package of Ab initio Programs. 2003.
19. Foster JM, Boys SF. Rev. Mod. Phys. 1960;32:296.
20. Simon S, Duran M, Dannenberg JJ. J. Chem. Phys. 1996;105:11024.
21. Boys SF, Bernardi F. Mol. Phys. 1970;19:533.
22. Tsuzuki S, Lüthi HP. J. Chem. Phys. 2001;114:3949.
23. Zhang Y, Pan W, Yang W. J. Chem. Phys. 1997;107:7921.