The Performance of Generalized Gradient Approximation DFT Methods with Gaussian Basis Sets: Sulfur-containing Molecules

Gábor I. Csonka, N. Anh and J. Réffy

Department of Inorganic Chemistry, Technical University of Budapest

H-1521 Budapest, Hungary

E-mail: csonka@iris.inc.bme.hu

Abstract

The performance of Gradient Generalized Approximation Density Functional Theory (GGA DFT) methods with gaussian basis sets is examined by studying 5 small molecules. The geometries are optimized by HF, MP2 and DFT methods. We apply four different DFT functionals. The gradual improvement of basis sets gradually decrease the bond lengths and increase the bond angles. Accidentally the HF/6-311G(d) results are close to the experimental results and the improvement of the basis sets to 6-311G(2d,f) worsen the agreement with the experiment. The inclusion of the electron correlation effect increase the bond lengths considerably. The various GGA DFT results agree qualitatively with each other and with the MP2 results, however some functional provide exaggerated effects and poor agreement with the experiment while others yield reasonable results.

Introduction

Density functional theory (DFT) methods have received considerable research attention recently. DFT methods are attractive because the best DFT functionals include the electron correlation at smaller computational cost than the traditional correlated methods. There is a growing literature on systematic comparisons of DFT results with experiment, HF, MP and QCISD(T) results [1]. The quality of the DFT results depend on the quality of the functionals applied. It is a difficult task to construct an appropriate functional for the exact exchange-correlation hole around the electrons. A list of the physical conditions that the exact exchange-correlation hole satisfies may help the construction and testing of actual and future functionals [2]. The study of Perdew et al. [2] shows that none of the actually used functionals satisfies all the known exact conditions and it is not clear that the violations of the exact conditions which way will degrade the quality of the results. For example the correlation functional of Lee, Yang and Parr violates a large number of exact conditions, however it was successfully used for chemical problems with Becke's exchange functional (B-LYP). Moreover the form of the exact exchange is known from the HF theory and some so called hybrid functional include it in a semiempirical way. The lack of the methods which gradually improve the quality of the DFT methods make them less attractive. Actually the best way to test a functional is to perform systematic comparisons between the results calculated with various functionals and the best theoretical estimations and experimental results. Apart from the functionals another problem is the basis set dependence of the DFT results. This was recently examined for sulfur-containing molecules [3]. In this study we present the functional-dependence of the DFT results for sulfur, chlorine and fluorine-containing molecules.

Theoretical methods

The geometries of selected molecules were fully optimized at HF level of theory supplemented with the triple-split-valence plus polarization 6-311G(d), G(2d) and G(2df) basis sets [4]. This basis set involves 5 d and 7 f orbitals and the McLean-Chandler (12s,9p)[6s,5p] basis sets [5] for second row atoms.

We also combine these basis sets with four different density functional (DFT) methods. We selected the promising generalized gradient approximation (GGA) density functionals of Becke [6,7] (B and B3), Perdew [8] (P86) and Lee, Yang and Parr [9] (LYP). The combinations of the functionals were as follows:

(i) B-P86 or Becke-Perdew method, in which Becke's exchange functional is combined with Perdew's correlation functional.

(ii) B3-P86 is a hybrid method. It is a linear combination of various exchange and correlational functionals in the form:

A*Ex[S] + (1-A)*Ex[HF] + B*Ex[B]+Ec[VWN]+C*Ec[P86].

where Ex[S], Ex[HF] and Ex[B] are the Slater, HF and Becke exchange functionals; and Ec[VWN] and Ec[P86] are the Vosko, Wilk and Nussair [10] and Perdew correlation functionals respectively.

The constants A, B and C are those determined by Becke by fitting heats of formation [4].

(iii) B-LYP method, in which Becke's exchange functional is combined with the correlation functional of Lee, Yang and Parr.

(iv) B3-LYP is a hybrid method. This functional was not published before the implementation into the GAUSSIAN 92/DFT [11]. It is a logical extension of Becke's three parameter concept using different correlational functionals (e.g. LYP). The following table show the parametrization of the method. For comparison purposes we show the popular local spin density approximation (LSDA) and the B3-P86 functionals in the same table:

            Ex[S]    Ex[HF]    Ex[B]    Ec[VWN]     Ec[P86]    Ec[LYP]   
LSDA          1                            1                             
B3-P86        A       1-A        B         1           C                 
B3-LYP      0.80      0.20     0.74       0.19                  0.81     

GAUSSIAN 92/DFT uses numerical quadrature to evaluate the DFT integrals. The quadrature scheme is defined by the number of points radial and angular directions. The geometries were optimized with a normal grid (50 radial shells, 194 angular points per shell, pruned to about 3000 points per atom).

For comparison, we also calculated the molecular properties with correlated single reference many-body perturbation theory (MBPT) using Moller-Plesset partition [12] to the second order (MP2). In the MP2 calculations the atomic core was frozen and they were supplemented with the 6-311G(d) basis set.

Results and discussion

We selected the Cl2S, Cl2S2, ClS2H, CS2 and F2S molecules for this study. The triatomic molecules show similar tendencies: the HF theory provides the shortest bond length and the smallest bond angle for these molecules. The improper treatment of electron correlation causes bond lengths to be underestimated in HF calculations supplemented with extended basis sets (e.g. 6-311+G(3df,2p). In these type of calculations the electron-electron repulsion terms are underestimated, and the electrons tends to move too close to each other. The stronger coulombic attraction between electrons and nuclei results compact bonds. Earlier studies show a general inverse relationship between the quality of the basis set and the predicted bond lengths. In many molecules small basis sets usually overestimate the experimental bond lengths. Therefore, it is often possible to find a basis set of intermediate quality which gives a geometry very close to the experiment. At this level of theory the results are better than those obtained with much higher level of theories and this is due to opposing nature of correlation effects and basis set corrections. The basis set effects tend to decrease the correlational effects tend to increase the bond lengths. The present results for Cl2S, CS2 and F2S support this observation.

Figure 1 shows the calculated S-Cl bond length in SCl2 molecule versus the theoretical method with three different basis sets. The experimental result is a horizontal line and the MP2/6-311G(d) result is shown by an arrow. (The same notation will be used on the following figures). The figure shows clearly the opposing nature of correlation and basis set effects. The HF/6-311G(2df) S-Cl bond length agree exactly with the experimental value and this is the shortest bond length at HF level of theory. The figure also shows that the correlation effects lengthen considerably the S-Cl bond as it is expected. The electron correlation decreases the electron density in the bond centre and the nuclear repulsion leads to longer bonds. The question is that how large should be this effect? The various functional shown in the Figure increase gradually the correlation effect. The B3-P86 and the B3-LYP methods provide smaller the B-P86 and B-LYP methods provide larger effects. It is clear that the mixture of the exact HF exchange and the DFT exchange functional results better agreement with the HF results. The LYP correlational functional seem to introduce considerably stronger correlational effect than the P86 functional. This puts the functionals into an order, the smallest correlational effect method is the B3-P86, then the B3-LYP, B-P86 and B-LYP follows. The B-LYP seems to overcorrect, and it provides the worst agreement. It interesting that the HF results agree well with the experimental result. They should not, maybe the basis set is not large enough or the experimental bond length is too short. The correlation functionals do that they are supposed to do lengthen the bond length (the MP2 does the same). Consequently if the HF results agree with the experimental result the inclusion of electron correlation (MP2, CCSDT, DFT etc.) will worsen the agreement with the experiment for triatomic systems.

Figure 2 shows the calculated Cl-S-Cl bond angle in SCl2 molecule versus the theoretical method with different basis sets. The bond angle sows small basis set dependence and strong functional dependence. The bond angle seem to be converged at HF level. The HF and MP2 methods provide excellent agreement with the experimental value. The DFT functionals provide the same order as on the Figure 1 and the bond angle is too large. Again the B-LYP provide the worst agreement with the experiment.

For the four-atomic molecules the situation is somewhat more complicated, because of the dihedral bond interactions. Figure 3 shows the internal S-S bond length in S2Cl2 molecule. The HF method provide the longest S-S bond the shortest S-Cl bond (Fig. 4), smallest S-S-Cl angle (Fig. 5) and the smallest Cl-S-S-Cl dihedral angle (Fig. 6). All GGA-DFT methods perform equally for the S-S bond length, and they agree well with the MP2 results, however none of them is capable to reproduce the short experimental S-S bond length. The B3-P86 method provides the best agreement with the experiment.

Figures 4 and 5 show that for the atoms at the end of the chain behave the same way as in the triatomic molecules.

Figures 7-12 show the various bond lengths and bond angles in S2ClH molecule. The behaviour of S-S-Cl part of the molecule (Figs 9 and 11) is similar to the S-S-Cl of the S2Cl2 molecule. For the S-S-H end the B-P86 and B-LYP methods provide almost equal results (Figs 8 and 10), the main difference is in the exchange part. The B3 exchange provide considerably shorter bond than the B exchange.

Figure13 shows the C-S bond length in CS2 molecule. Figures 14 and 15 show the S-F bond length and the F-S-F bond angle in SF2 molecule. The same conclusions can be drawn as previously.

These results show that it possible to fine tune the correlation effects by selecting a series of GGA-DFT functionals. Their performance show some variation depending on the bond, however on the basis of the correlation strength there is a general order between the various functionals for the presented molecules. In increasing order: B3-P86, B3-LYP, B-P86 and B-LYP. The basis set dependence of the GGA-DFT methods is generally larger if a second d-function added to to the basis set and it is smaller if an f-function is added to the basis set compared to HF method. The B-LYP method shows the largest basis set dependence. In general the basis set dependence of the GGA-DFT methods is not larger than that of the HF method for the presented molecules and basis sets.

Acknowledgements

The authors are indebted to Professor N. C. Handy and Dr. J. A. Altmann for providing their submitted manuscript. The financial support of the Hungarian Research Foundation (OTKA I/3 No. 644 and T 014247) is acknowledged.