Basis Set and Method Dependence in Atoms in Molecules Calculations

Jan 19, 2010 - Department of Quantum Chemistry, Nicolaus Copernicus UniVersity, PL-87 100 Torun, Poland, and. Department of Crystallography and Crysta...
1 downloads 0 Views 326KB Size
Download PDF
2240

J. Phys. Chem. A 2010, 114, 2240–2244

Basis Set and Method Dependence in Atoms in Molecules Calculations Mirosław Jabłon´ski*,† and Marcin Palusiak‡ Department of Quantum Chemistry, Nicolaus Copernicus UniVersity, PL-87 100 Torun´, Poland, and Department of Crystallography and Crystal Chemistry, UniVersity of Ło´dz´, PL-90 236 Ło´dz´, Poland ReceiVed: NoVember 20, 2009; ReVised Manuscript ReceiVed: December 30, 2009

The influence of various basis sets used in HF and DFT/B3LYP calculations to the values of atoms in molecules (AIM) parameters derived from the electron density distibution for weak hydrogen-bonded systems is investigated. Using three model complexes, F3CH · · · NH3, F3CH · · · NCH, and FCCH · · · NH3, we show that values of the most important AIM parameters calculated in the bond critical point of the H · · · N hydrogen bond are almost independent of both the method and the basis set. Only the smallest Dunning-type cc-pVDZ or aug-cc-pVDZ basis sets may lead to poor results, whereas even medium-sized Pople-type basis sets can give reasonable results converting to those obtained from the use of large Dunning-type basis sets. Introduction All atomic and molecular properties are governed by the electron density distribution. Thus, the methods that deal with the analysis of the electron density distribution should have a particular appeal for chemists and help to understand the electron structure of molecules and thus interactions. Important information about the charge distribution and about its changes due to intermolecular interactions can be obtained by means of population analysis, where the total number of electrons is divided to individual atoms or even to atomic orbitals. The most popular are the population analysis introduced by Mulliken1 and Lo¨wdin.2 Still increasing weight is also attached to the natural bond orbital (NBO) method proposed by Weinhold et al.3 All these methods try to represent the molecule as a distribution of point charges. Unfortunately, their definitions are basis set dependent and specify some characteristic points that reflect the use of atom-centered basis sets. This disadvantage of traditional methods of the electron density distribution is remedied in the atoms in molecules (AIM) theory of Bader.4-9 A great interest and widespread use of this method is due to the fact that, as opposed to the wave function quantum mechanics, the AIM method gives the opportunity to have an insight into a region of a system. This opportunity meets the interest of most chemists who wish to have a theoretical tool to study a small part of a molecule only, instead of dealing with the total energy of a whole system. The goal of the AIM theory is to find a physically well-defined surface that divides a system to subsystems. In the AIM theory the surface bounding the subsystem is established by the electron b More strictly, this bounding density gradient vector field (∇F). surface is established is such a way that there must not be a b through the surface. Thus a surface satisfying this flux in ∇F condition is known as “a zero flux surface”.7 This surface passes b )b through a point called a critiacal point (CP), where ∇F 0. Of particular interest is a critical point in which the electron density is a minimum with respect to the direction of a line connecting two nuclei and a maximum with respect to all directions perpendicular to this line, thus being a saddle point * To whom correspondence should be addressed. E-mail: teojab@ chem.uni.torun.pl. Phone+48 (56) 6114695. Fax: +48 (56) 6542477. † Nicolaus Copernicus University. ‡ University of Ło´dz´.

with the (3,-1) curvature. A critical point with such characteristics is called a bond critical point (BCP) and indicates the accumulation of the electron density. Thus BCP can be used in the recognition of chemical bonds between atoms and then in studies of the character of this interatomic interaction. Large attention is paid to the proper values of the electron density, FBCP, and its Laplacian, ∇2FBCP, in BCP of the XH · · · Y interaction, since they may indicate a hydrogen bond if these parameters have proper values.10 The Laplacian of FBCP has a dominant role in the AIM analysis, since it determines where the electron density is locally concentrated (∇2F < 0) and depleted (∇2F > 0) within a molecular system.5 Moreover, ∇2FBCP is related6 to the local potential and kinetic energy densities (V(r b) and G(r b), respectively) that are components of H, the total energy density. Most of AIM calculations are performed with the use of rather small basis sets; thus one may ask about the dependence of values of AIM parameters on both the basis set and the method used in AIM calculations. In other words, one may ask about the reliability of such calculations, i.e., about the stability of values of AIM parameters with respect to the change of both the basis set and the method used in AIM calculations. The aim of this paper is to study this issue. Methodology To study the stability of AIM calculations with respect to the basis set and the method, we chose three model complexes: F3CH · · · NH3, F3CH · · · NCH, and FCCH · · · NH3. The choice of weak hydrogen-bonded systems was thought to find differences between values of AIM parameters (if any) on various levels of approximation, since all variations due to the use of a given basis set and method should be a substantial part of investigated AIM parameters. Various AIM parameters commonly used in the description of inter- and intramolecular interactions were then computed using a large set of medium and large basis sets of Pople11 and Dunning12,13 type. Small basis sets of the Pople type were not used in this paper since they give unsatisfactory results in calculations concerning intermolecular interactions. HF11 and DFT/B3LYP14-17 approximations were used. Geometry optimizations were performed at the level of B3LYP/aug-cc-pVTZ for collinear arrangement of individual monomers to use symmetry relations and thus to simplify

10.1021/jp911047s  2010 American Chemical Society Published on Web 01/19/2010

Basis Set and Method Dependence in AIM Calculations

J. Phys. Chem. A, Vol. 114, No. 5, 2010 2241

Figure 1. Values of chosen AIM parameters depending on a basis set (horizontal axis) and method (HF on the left-hand side, DFT/B3LYP on the right-hand side) for weak H-bonded complexes: F3CH · · · NH3 (upper row), F3CH · · · NCH (central row), and FCCH · · · NH3 (lower row). Basis sets are given numerical values: (0) 6-31++G(d,p), (1) 6-311++G(d,p), (2) 6-311++G(2df,2pd), (3) 6-311++G(3df,3pd), (4) cc-pVDZ, (5) aug-ccpVDZ, (6) cc-pVTZ, (7) aug-cc-pVTZ, (8) d-aug-cc-pVTZ, (9) cc-pVQZ with g-functions removed, (10) aug-cc-pVQZ with g-functions removed. Points have been connected for a better visualization of results.

calculations. It should be added that in fact performing a fully unconstrainted geometry optimization should not be so important here, since the aim was to compare values of AIM parameters obtained by means of different approximations for systems with the same geometries, whatever the model of the geometry construction, and to ensure the existence of only one hydrogen bond. Results following from this part of calculations are discussed in the first part of the next section. To study the influence of both the method and the basis set used in AIM calculations to the relation between either the hydrogen bond energy or the H · · · N distance and commonly the most often used AIM parameters, we also performed AIM calculations for a wide range of H · · · N distances in the F3CH · · · NH3 complex. In this case, for monomers, geometry parameters were used as from their isolated states to ensure the proper value of the interaction energy in the dissociation limit. Calculations with the H · · · N distance scan were performed for

two basis sets: 6-311++G** and aug-cc-pVTZ, together with HF or DFT/B3LYP approximations. Results of these calculations are discussed in detail in the second part of the next section. Calculations were performed with the use of the Gaussian03 set of codes.18 A detailed analysis of the electron distribution function was made according to the concept of the “atoms in molecules” theory proposed by Bader,4-9 using the AIM2000 program.19,20 Properties of the electron density calculated at bond critical points (BCP) of H-bonds under consideration were characterized. The electron density at BCP, FBCP, its Laplacian, ∇2FBCP, the density of the total electron energy, HBCP, and its two components (densities of potential, VBCP, and kinetic, GBCP, electron energies) were discussed. Stability of AIM Calculations Basis Set and Method Dependence. Most of AIM calculations are performed using rather small Pople type basis sets.

2242

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

The aim of this paper is to analyze the reliability of AIM calculations with the use of medium and large Pople type basis sets by comparing with results obtained with the use of larger Dunning type basis sets, i.e., of cc-pVXZ or aug-cc-pVXZ type, where X ) D or T. These calculations have been widened to d-aug-cc-pVTZ, as well as to cc-pVQZ and aug-cc-pVQZ basis sets with g-functions removed in the case of the QZ quality basis sets, because of technical limitations. For the purposes of our analysis we chose three weak hydrogen-bonded complexes, namely, F3CH · · · NH3, F3CH · · · NCH, and FCCH · · · NH3. Parameters of the electron density function (F, its Laplacian, H, and its components: V and G) were measured at critical points corresponding to the proper HBs linking H-bonded complexes under investigation. The graphical representation of the interrelation between values of these parameters and basis sets used in these calculations is shown in Figure 1. The first conclusion that flows from Figure 1 is most likely that the pattern of point distribution is almost the same for all six figures (see Figure 1), irrespective of whether HF or DFT/ B3LYP calculations were performed for a given complex and irrespective of which complex these AIM calculations concern. The electron density in a bond critical point, FBCP, is surely the most important AIM parameter used in the field of intra- and intermolecular interactions. This is because it is often correlated to values of the interaction energy.21-27 For all but one basis sets investigated in this paper we obtained practically the same values of electron density in the bond critical point of the H · · · N interaction for a given system. Only a poor cc-pVDZ basis set leads to noticeably higher values; however, the augmentation of this basis set to aug-cc-pVDZ corrects this value, leading to the result being closer to that obtained also with the use of much larger Dunning type basis sets (see Figure 1).

Jabłon´ski and Palusiak Much larger differences in absolute values, on comparing various basis sets, were obtained for the Laplacian of the electron density in BCP, ∇2FBCP. For ∇2FBCP both cc-pVDZ and augcc-pVDZ (and in a smaller amount 6-31++G(d,p) basis set as well) basis sets lead to remarkably smaller values than in the case of larger Pople or Dunning type basis sets. If the cc-pVDZ basis set is enlarged to cc-pVTZ, one obtains a much higher value of ∇2FBCP for both HF and DFT/B3LYP calculations; nevertheless, enlarging this basis set to aug-cc-pVTZ, d-augcc-pVTZ, or QZ quality basis sets leads to a fast convergence. It is noteworthy that the medium and the large Dunning type basis sets used in this paper lead to values of ∇2FBCP similar to those obtained from Pople type basis sets (not valid for the smallest 6-31++G(d,p) basis set). The total electron energy density, HBCP, is a sum of two components having opposite signs: potential electron energy density, VBCP, which is negative, and kinetic electron energy density, GBCP, which is positive. Thus the value of HBCP, for a given method and a basis set, depends on both the value of VBCP and the value of GBCP and should be close to zero for weak interactions. Again, it is seen from Figure 1 that patterns of distribution of points indicating values of HBCP, VBCP, and GBCP are very much similar for all six figures, irrespective of the method and the complex. Somewhat smaller values of HBCP are obtained for the smallest basis sets within each group of basis sets used here, i.e., of Pople and Dunning type. Summarizing, one can say that values of AIM parameters are practically independent of the basis set used in HF or DFT/ B3LYP calculations. This conclusion is also valid for the Laplacian of the electron density in BCP, which is in general more sensitive to changes of basis sets. Moreover, in this case, values similar to those of cc-pVDZ and aug-cc-pVDZ are also

Figure 2. Dependence between AIM parameters (in au) calculated at the BCP of the H · · · N hydrogen bond of the F3CH · · · NH3 complex with fixed d(H · · · N) distances and the hydrogen bond energy. Point distibutions were fited by functions of the f(E) ) A(e-aE - 1) type for better visualization of results.

Basis Set and Method Dependence in AIM Calculations obtained with the use of the smallest 6-31++G(d,p) basis set, which is also of split-valence double-ζ type. It is worth mentioning here a very recent report by Matta.28 In that report the comparison between AIM-derived data obtained with the use of different methods (in conjunction with the 6-311++G(d,p) basis set) was performed for biologically relevant set of molecules. It was found that also the level of theory does not influence essentially the reliability of AIM analysis. AIM Parameters vs EHB. Although in the literature one can find some correlations between the hydrogen bond energy and the electron density in the bond critical point of this hydrogen bonding,21-27 similar correlations for other AIM parameters are hardly found. Thus it seems attractive enough to study correlation patterns between some AIM parameters, the most often met in the field of intra- and intermolecular interactions, and the hydrogen bond energy. To study the influence of the level of approximation used in AIM calculations on these correlations, we have performed scan calculations for the F3CH · · · NH3 complex.BothHFandDFT/B3LYPcalculations,with6-311++G** and aug-cc-pVTZ basis sets, each belonging to the group of triple-ζ split-valence basis sets with both polarization and diffuse functions have been performed. Both these basis sets appear to be the smallest reasonable in each group (i.e., of Pople and Dunning type). For each fixed distance between interacting monomers we have kept geometry parameters of individual monomers frozen from their isolated fully optimized states. In this way one obtains a proper value in the dissociation limit. Both the decrease and the increase of the equilibrium distance leads to the increase of the potential energy; thus, to avoid doubling of similar values of EHB coming from both sides of the minimum on the potential energy curve, we have only investigated the right part of it. Results showing the influence of both the method and the basis set used in AIM calculations

J. Phys. Chem. A, Vol. 114, No. 5, 2010 2243 to relations between various AIM parameters and the hydrogen bond energy are presented in Figure 2. It is seen from Figure 2 that point distributions corresponding to various levels of approximation used in calculations of AIM parameters have similar shapes, although they are shifted with respect to each other as a result of different equilibrium energies. It means that the correlation pattern does not depend much on the method and the basis set used in AIM calculations, assuming that at least a basis set of split-valence triple-ζ quality with polarization and diffuse functions is used. It may be particularly important in calculations for larger systems, where the use of more time-consuming larger basis sets is unwilling or even impossible in a reasonable time. It is in agreement with chemist’s intuition and literature21-27 that absolute values of all AIM parameters investigated in this paper grow monotonically with the increase of the hydrogen bond energy (see Figure 2). Thus, a larger value (for VBCP, absolute) of AIM parameter can in general indicate a stronger hydrogen bond. This conclusion should be valid at least for weak hydrogen-bonded complexes as the F3CH · · · NH3 complex investigated herein. It is noteworthy that the correlation between the electron density in BCP of the H · · · N interaction is not linear with respect to the energy of this interaction. It is clearly seen from Figure 2. Nonlinear dependence between FBCP and EHB has been already obtained earlier by Grabowski26 for a group of malonaldehyde derivatives, where the quadratic regression

EHB ) 1805.59FH · · · O2 - 302.1FH · · · O - 0.7161

(1)

with R2 ) 0.928 was proposed.26 In the case of small range of FBCP values this relation can be expressed with linear regression,27 as has been shown for various complexes with meth-

Figure 3. Dependence between AIM parameters (in au) calculated at the BCP of the H · · · N hydrogen bond and the H · · · N distance (in Å) in the F3CH · · · NH3 complex shown for a few levels of approximation.

2244

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

oxybenzene, where the H-bond energy was found on the basis of the EHB ) 1/2VBCP relation proposed by Espinosa et al.23 Since, for a given Y, there is a relation between the hydrogen bond energy and the H · · · Y distance, there should be also relations (and rather smooth) between values of AIM parameters in the H · · · Y BCP and the H · · · Y distance. The influence of both the basis set and method used in AIM calculations on these relations for the F3CH · · · NH3 complex is depicted in Figure 3. Relations between investigated AIM parameters and the H · · · N distance in the F3CH · · · NH3 complex are of similar character irrespective of both the method and the basis set. Since the smaller distance between interacting subsystems in general means the larger interaction energy, the smaller value of d(H · · · N) leads to larger values of AIM parameters (for VBCP, the absolute value). Interestingly enough, it is not valid for HBCP. In this case one obtains a significant decrease of HBCP values for smaller H · · · N distances and negative values can even be obtained for very large shortening of the H · · · N distance (see DFT results in Figure 3). This is because the density of the electron potential energy, which is the negative component of HBCP, becomes the overwhelming factor. Such Morse-type dependence of HBCP with the intermolecular separation distance has already been found by Dominiak at al.29 for DMAN complexes and Schiff bases in crystals. A similar Morse-type curve should be also obtained for ∇2FBCP in the case of shorter separation H · · · N distances.29,30 Conclusions AIM calculations are often used in investigations of intraand intermolecular interactions, and thus they are very helpful in studies of hydrogen bonds. It flows from the fact that many correlations between various AIM parameters and the hydrogen bond energy or some other parameters that indirectly describe the strength of the hydrogen bond were found.21-27,29 These calculations, however, were performed using small or medium basis sets of Pople type. In this paper we show that such an approach is rather justified, since values of AIM parameters seem to be independent of both the method (HF or DFT/B3LYP) and the basis set used in the AIM calculations. Merely poor cc-pVDZ or aug-cc-pVDZ basis sets of Dunning type should rather be omitted, since they may lead to wrong values of AIM parameters compared to those obtained from much more reliable basis sets of split-valence triple-ζ or even better quality. The fact that AIM calculations seem to be rather method/basis set independent may be very helpful in much more time-consuming calculations for large systems. The generalization of the method/

Jabłon´ski and Palusiak basis set independency of AIM calculations should be studied also for systems with stronger hydrogen bonds, although any variations due to the change of the method/basis set would be not so substantial a part of the value of the AIM parameter in this case. Acknowledgment. Calculations have been carried out on the multiprocessor cluster at the Information and Communication Technology Center of Nicolaus Copernicus University, Torun´ (http://www.ucu.umk.pl). References and Notes (1) Mulliken, R. J. Chem. Phys. 1955, 23, 1833. (2) Lo¨wdin, P. O. J. Chem. Phys. 1950, 18, 365. (3) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. ReV. 1988, 88, 899. (4) Bader, R. F. W. Acc. Chem. Res. 1975, 8, 34. (5) Bader, R. F. W.; McDougall, P. J.; Lau, C. D. H. J. Am. Chem. Soc. 1984, 106, 1594. (6) Bader, R. F. W.; Esse´n, H. J. Chem. Phys. 1984, 80, 1943. (7) Bader, R. F. W. Acc. Chem. Res. 1985, 18, 9. (8) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: New York, 1990. (9) Popelier, P. L. A. Atoms in Molecules. An Introduction; Longman: Singapore, 2000. (10) Koch, U.; Popelier, P. L. A. J. Chem. Phys. A 1995, 99, 9747. (11) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to AdVanced Electronic Structure Theory; Macmillan Publishing Co., Inc.: New York, 1982. (12) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 9747. (13) Davidson, E. R. Chem. Phys. Lett. 1996, 260, 514. (14) Hohenberg, P.; Kohn, W. Phys. ReV. 1964, 136, B864. (15) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133. (16) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (17) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (18) Frisch, M. J.; et al. Gaussian 03; Gaussian, Inc.: Wallingford, CT, 2004. (19) Biegler-Ko¨nig, F. W.; Bader, R. F. W.; Tang, T. H. J. Comput. Chem. 1982, 3, 317. (20) Biegler-Ko¨nig, F. W. AIM2000 program; 2000. (21) Boyd, R. J.; Cheng Choi, S. Chem. Phys. Lett. 1985, 120, 85. (22) Boyd, R. J.; Cheng Choi, S. Chem. Phys. Lett. 1986, 129, 62. (23) Espinosa, E.; Molins, E.; Lecomte, C. Chem. Phys. Lett. 1998, 285, 170. (24) Grabowski, S. J. J. Phys. Chem. A 2000, 104, 5551. (25) Grabowski, S. J. J. Phys. Chem. A 2001, 105, 10739. (26) Grabowski, S. J. J. Mol. Struct. 2001, 562, 137. (27) Palusiak, M.; Grabowski, S. J. J. Mol. Struct. 2002, 642, 97. (28) Matta, C. F. J. Comput. Chem., DOI: 10.1002/jcc.21417. (29) Dominiak, P. M.; Makal, A.; Mallinson, P. R.; Trzcinska, K.; Eilmes, J.; Grech, E.; Chruszcz, M.; Minor, W.; Woz´niak, K. Chem.sEur. J. 2006, 12, 1941. (30) Grabowski, S. J.; Małecka, M. J. Phys. Chem. A 2006, 110, 11847.

JP911047S