Role of Hartree-Fock exchange in density functional theory

18 November 1994

CHEMICAL PHYSICS LETTERS

Chemical Physics Letters 230 (1994) 189-195

ELSEVIER

Role of Hartree-Fock exchange in density functional theory. Some aspects of the conformational potential energy surface of glycine in the gas phase Francesco Lelj a,*,Carlo Adamo a, Vincenzo Barone b aDipartimento di Chimica, Universitd della Basilicata, via N. Sauro 85. 85100 Potenza, Italy b Dipartimento di Chimica. Universitir “Federico II”, via Mezzocannone 4, 80134 Napoli, Italy Received 28 June 1994; in final form 7 September 1994

Abstract A density functional study of the glycine molecule using different basis sets in the local spin density approximation with gradient corrections is reported. The results compare well with previous calculations which take into account correlation by different post Hartree-Fock procedures and with microwave and electron diffraction studies. Energetic results show that the OH-N containing conformation is more stable than the conformation of C, symmetry containing NH2-OH bifurcated H-bonds. In the case of this family of conformations the C, symmetry conformation is a saddle point. Furthermore the results show that inclusion of a fraction of the exact Hartree-Fock exchange improves both the molecular geometries and energetics.

1. Introduction

Glycine is the smallest amino acid. Apart from its intrinsic interest as an interstellar amino acid, which may provide insight into the role of second-row atoms in prebiotic events, it represents a nice benchmark for assessing the quality of different computational schemes in determining structural and electronic properties and intramolecular interactions in biological molecules. Moreover this molecule is an interesting example of the relevance of the dialectics between high-quality computational results [ l-31 and experimental investigations. In fact the finding of a new conformer in the gas phase by microwave spectroscopy [ 21 was triggered by the increasing accuracy of ab initio computations [ 4 ] and recently the * Corresponding author. 0009-2614/94/$07.00

gas phase structure has been determined by electron diffraction [ 5 1. The glycine molecule contains three internal rotations around H,N-C, C-COOH and C-OH bonds (Fig. 1). Different experiments and theoretical studies [ 2,3,5,6] have pointed out that in the gas phase among the many conformations available two of the conformers are mainly populated. The former is characterized by the H-bond interaction between the hydrogen atoms of the NH2 moiety and the oxygen of the carbonylic group, and the latter by the interaction between the hydrogen atom of the OH moiety with nitrogen. The most recent experimental studies in the gas phase have been performed by microwave [ 1 ] and electron diffraction [ 5 ] techniques. The former, further supported by Sellers and Schafer’s calculations [ 31 suggested that the most stable conformer was I whereas the latter indicated II as the most stable.

0 1994 Elsevier Science B.V. All rights reserved

SSD10009-2614(94)01156-7

190

F. Lelj et al. /Chemical

Physics Letters 230 (1994) 189-195

III

(GC) [ 12,13-l 5 1, which include a large amount of the correlation energy, appear to be a promising tool for studying systems in the gas phase [ 161 and model solutions [ 17,18 1. It has been shown by different authors [ 7,16,19,20 ] that GC, despite showing, at least in model Hamiltonians, some unexpected behavior [ 2 11, are essential for describing energies and geometries close to the values requested for chemical utility even though some problems still remain in reproducing heavy atom hydrogen bond lengths [ 221. Recently, the suggestion [ 23-25 ] of taking into account non-local exchange contributions by means of a term proportional to the Hartree-Fock exchange has greatly improved the agreement with experiment at least in the case of non-transition metal containing compounds [ 25 ] and there are indications [ 26 ] that in the case of 3d transition metal element monocarbonyls the improvement is of the same quality.

(Cl)

2. Methods Fig. 1. Plot of the three most stable conformations of glycine as computed by B3LYP. II is the transition state between I and III.

Conformer II, with C, symmetry, has been shown to be a first-order saddle point (SP) on the BomOppenheimer potential surface at the Hartree-Fock (HF) and second-order Moller-Plesset (MP2) levels, being the SP for the interconversion between I and III [ 71. Further theoretical calculations with an improved treatment of correlation have confirmed that II is less stable than I and III, but no check of the SP nature of II has been undertaken [ 81. After many HF studies [ 3,4,9, lo] it has been shown [ 68] that a correlated level of theory is necessary for reproducing the glycine structure and for describing its conformational behavior. These results also appear to be relevant in suggesting geometrical parameters [ 51 essential in some investigations where the model, describing the molecule for quantitative interpretation of the experimental data, contains a number of parameters larger than the number of independent observations [ 51. Among the different theoretical schemes for computing molecular properties, density functional (DF) theories [ 111 with different gradient corrections

The density functional computations have been carried out using several different local density gradient corrected approaches: the standard Becke [ 14,15 ] exchange functional with either Lee-YangParr [ 141 or Perdew [ 131 correlation functionals (BLYP and BP respectively), and a mixed approach which incorporates some contributions from the Hartree-Fock exchange as suggested by Becke [23,24]. Different basis sets (Slater and Gaussian type) have been used and all the conformers have been characterized by computing the harmonic frequencies, which have previously been computed only at the Hartree-Fock [ 31 and MP2 [ 7 ] levels. The ADF ’ and GAUSSIAN 92/DFI [28] set of programs have been used for computing geometrical structures and Hessian matrices by the number and analytical algorithms of the two programs respectively. In the case of ADF, a double-zeta ST0 basis set [ 29 ] was used plus a set of polarization functions. Fit functions including up to g-STOs were used for all ’ A modified version of Amsterdam Density Functional System (ADF), Department of Theoretical Chemistry, Vrije Universiteit, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands 1271.

F. I_& et al. /Chemical Physics Letters 230 (1994) 189-195

second row atoms and up to d-STOs for hydrogen atoms [ 29 1. In the case of second row atoms all 1s cores were kept frozen and one ST0 was added in order to ensure proper orthogonalization. All integrals were computed with accuracy parameter set to 5 according to the te Velde and Baerends integration scheme [ 30 1. Standard default criteria for density of integration points were chosen. Local spin density (LSD) contributions have been computed using Slater exchange and Vosko, Wilk and Nusair (VWN) correlation functionals [ 3 11. Gradient corrections have been introduced using Becke exchange [ 151 and Perdew correlation [ 13 ] (BP). In the computations carried out by the GAUSSIAN 92/DFT package, gradient corrections have been introduced using Becke exchange [ 151 and LeeYang-Parr [ 141 correlation (BLYP) or the Perdew contribution [ 13 ] (BP). Becke has argued [24] that “exact exchange energy must play a role in highly accurate density functional theories”. Starting from the adiabatic connection formula [ 23 ] Eq. ( 1) represents the simplest possible mixture of local density gradient corrections and exact exchange that recovers the uniform electron gas limit: E,, =EE’” +A,(E,HF -E,“s”) +A, &,+A,

MC 3

(1)

where the AE terms represent gradient corrections. We will use here Becke’s gradient correction to exchange, but the Perdew and Wang (PW) correlation potential, used in Ref. [ 241, has been replaced either by the Perdew (leading to the B3P functional) or by the LYP one (leading to the B3LYP functional ). Since the LYP functional contains both a local part and a gradient correction, only the latter contribution should be used to obtain a coherent implementation. As a consequence we approximate the last term in Eq. (1) by AE,, zEkyp - E,VWN. A number of tests showed that the values of the three semi-empirical coefficients appearing in Eq. ( 1) near 0.8 provide the best results, irrespective of a particular form of different functionals. We will use here the values (A,=0.20, A,=0.72 and ApO.81)

191

determined by Becke from a best fitting of the heats of formation of a standard set of molecules [ 241. Note that in the original application by Becke [ 241, HF exchange and GC were only added in the computation of energies using converged LSD electron densities. Furthermore single point energy computations were performed at experimental geometries. Here we use instead a fully coherent implementation in which the SCF process, geometry optimization and computation of second derivatives are performed with the complete density functionals including GC and possibly part of the HF exchange. Our reference basis set (hereafter referred to as DZP) is the Dunning double-c contraction [ 321 of the Huzinaga (9,5,1; 4,1) primitive set. We have also performed some computations with a larger (10,6,2;5,2)/ [4,3,2;3,2] basis set [33] (hereafter referred to simply as TZ2P), whose polarization functions are taken from the recent correlation consistent sets [ 341.

3. Results The whole set of accessible conformations has been studied [ 351 whereas in this preliminary account only the results on the three most stable conformations suggested by previous experimental and theoretical results are reported. Here we compare also the effect of using Gaussian or Slater type basis functions on geometries and relative energies of these three conformers. 3.1. Relative stability Concerning the relative stability of different conformers LSD-GC does not lead to the energy order for conformation II and III with respect to I as has been found using different basis sets and different correlation treatments (Table 1). This behaviour does not seem to be related to the particular exchange-correlation functional used (BP or LYP) nor to the basis set functional dependence (Slater versus Gaussian). In fact the energy of II with respect to I is - 1.335 and -0.229 kcal/mol in thecaseofBP/STO and BLYP/DZP respectively. Detailed study of the stability of conformations I, II and III shows that only the B3LYP functional is

192

F. Lelj et al. /Chemical Physics Letters 230 (1994) 189-195

Table 1 Relative energies (kcal/mol) of structures II and III with respect to glycine structure I. All the values refer to the corresponding optimized structures. In parentheses are reported the number of imaginary frequencies MP2J TZP b

HF/ DZP a I II III

CISD/ CCSD/ DZP’ DZP

O.OOO(Oi) O.OOO(Oi) 0.000 0.000 3.109(li) OS72(li) 1.901 1.431 2.9OO(Oi) 0.514(Oi) 1.769 1.330

“Ref. [8].

bRef. [7].

CCSD(T)/ DZP

MP2/ BP/ best estimate b ST0

0.000 1.060 1.009

0.000 0.560 0.492

BLYP / DZP

B3LYP J DZP

B3LYPJ TZZP

Exp. ’

0.000 O.OOO(Oi) O.OOO(Oi) O.OOO(Oi) 0.000 -1.335 -0.229(li) 0.217(li) 0.262(li) l.4~04 - 1.023 -0.267(Oi) O.l68(Oi) 0.246(01) ’

‘Ref. [9].

able to restore I as the most stable conformation and II as the saddle point between I and III though the absolute value of the energy remains lower than that computed by CISD ( 1.901 kcal/mol) and CCSD ( 1.4 13 kcal/mol) using a DZP basis set. On the other hand a larger basis set, including polarization and diffuse functions, reduces the energy differences among the three conformations to 0.572 and 0.514 kcal/mol (see Table 1). Our B3LYP data compare well with these values, being 0.262 and 0.246 kcal/ mol respectively, showing a slight increase in passing from the DZP to TZZP basis set. Most interestingly II is a SP between I and III, confirming that the second stable conformer has C, and not C, symmetry. The experimental energy stability of (II/III) with respect to I which amounts to 1.4 k 0.4 kcal/mol [ 61 appears large. It should be pointed out that the experimental value has been determined by some assumption about the partition function of the two states responsible for the intensity of the MW absorption. This could have effected the value of the error which was assumed to be within the value of 0.4 kcal/mol obtained as an estimation of the instrument error. Therefore it appears that even if the B3LYP result is in good agreement with the experimental value the energy difference between the conformations is probably underestimated and the experimental value is probably an upper bound. 3.2. Structure The use of LSD-GC gives, as already known, too long H-X bonds. This effect is present irrespective of the exchange-correlation functional expression: BP or BLYP. On the other hand bond lengths involving heavy atoms appear sensitive to the chosen basis set

and functional expression. In particular, STOs give geometries closer to the experimental values and to correlated calculation geometries as reported in Table 2. NC, CC and C=O bond lengths computed at BP/ST0 are closer to MP2, CISD or CCSD values than the BLYP/DZP ones. The introduction of the B3LYP functional gives better geometries in good agreement with CCSD ones. As a general trend the H-X bond lengths computed with B3LYP are closer to the values found from geometry optimization including correlation. Geometrical parameters obtained by ED study are biased by some assumptions of the structural parameters used for fitting the electron diffraction data obtained from older HF results [ 4, lo]. Furthermore, as already observed the NC bond length and NCC valence angle ( 1.467 8, and 112.1’ respectively) are strongly correlated in the ED data and their disagreement with previous CCSD/DZP (1.458 8, and 115.4” ) and our DFT/TZ2P results ( 1.45 1 A and 115.5 ’ ) suggests that the experimental values suffer from this. Of particular interest is the comparison between parameters computed and directly measured by microwave spectroscopy. Our B3LYP/TZ2P rotation constants (Table 3) are the closest till now obtained for the study molecules confirming that I is the most stable conformer. The closer values to the experimental data shown by III with respect to II suggest that the equilibrium structure has C, and not C, symmetry. Also the B3LYP/TZ2P dipole moments, in fairly good agreement with the experimental ones and with correlated wavefunction results, further support this hypothesis. Table 4 reports the DFT computed and unscaled harmonic frequencies and intensities compared with those previously computed. Our DFT results con-

F. L&j et al. / Chemical Physics Letters 230 (1994) 189-I 95 Table 2 Optimized

parameters

(8, and deg) for the most stable conformer

HF/ DZP a NH NC CH cc c-o OH C-0 HNH HNC NCC HCH HCC cc-o COH cc=0 C=G...N

1.002 1.440 1.086 1.517 1.331 0.95 1 1.191 106.4 110.5 115.3 106.1 107.3 111.7 108.6 125.6 2.794

CISD/ DZP =

CCSD/ DZP

BP/ STO-DZP

BLYP/ DZP

B3LYP/ DZP

B3LYP / TZ2P

Exp. c

1.014 1.447 1.094 1.519 1.356 0.968 1.209 106.2 110.2 115.6 106.1 107.4 110.9 106.6 125.7

1.011 1.445 1.090 1.518 1.344 0.962 1.204 105.3 109.5 115.2 106.2 107.4 111.5 107.1 125.6 2.771

1.021 1.458 1.098 1.525 1.359 0.972 1.216 104.7 108.8 115.4 106.3 107.5 111.5 105.9 125.7 2.782

1.034 1.454 1.107 1.526 1.367 0.999 1.214 104.2 108.6 115.0 106.1 108.2 111 105.3 125.7 2.848

1.029 1.471 1.105 1.543 1.379 0.985 1.228 104.7 108.6 115.3 105.9 107.7 111.2 105.9 125.8 2.911

1.019 1.455 1.097 1.529 1.360 0.974 1.215 105.4 109.3 115.2 105.8 107.6 111.5 106.7 125.7 2.877

1.012 1.451 1.091 1.522 1.356 0.968 1.204 105.7 109.9 115.5 105.6 107.7 111.5 107.1 125.7 2.867

(1.001) 1.467 (1.081) 1.526 1.355 (0.966) 1.205 (110.3) (113.3) 112.1 (107.0)

rotational

constants

HF/DZP

a

[8].

firm the intensities.

were assumed

(MHz)

MP2/TZP

from low level theoretical

and dipole moments b

CCSD/TZP

calculations

(D) for structures ’

(HF/4-2

I-III obtained

from optimized

geometries

BLYP/DZP

B3LYP/DZP

B3LYP/TZZP

Exp. ’

10342 3876 2912 % 1.00

10279 3877 2908 1.20

10229 3841 2884 1.24

9991 3760 2821 1.185

10227 3867 2881 1.239

10371 3855 2902 1.169

10297 4092 3028 5.98

10175 4076 3011

10069 4036 2980 5.79

9970 3973 2939 5.882

10156 4040 2990 5.953

10193 4059 3002 5.701

10196 4110 3054

10128 4085 3024

10015 4050 3001

9920 3978 2946

10105 4049 3001

10167 4065 3009

“Ref.

[7].

previous

5.59 “Refs.

5.61

111 (112.3) 125.1

1G, Ref. [ 41) .

10605 3917 2955 1.34

5.79 ‘Ref.

of glycine

MP2/ TZZP b

“Ref. [8]. bRef. [7]. ’ Ref. [ 5 1,the values in parentheses Table 3 Computed

193

5.779

5.844

5.644

10129 4071 3007 33

[1,2].

computed

frequencies

and

4. Conclusions Inclusion of contributions from HF exchange seems to correct to a good extent deficiencies shown by

‘simple’ LSD-GC functionals particularly in those cases where tiny energy differences are relevant as in the case of most stable glycine conformers in the gas phase. Our results further confirm the previous calculations including correlation and the experimental results (MW and ED) on the stability order of these

F. Lelj et al. / Chemical Physics Letters 230 (1994) 189-195

194 Table 4 Harmonic

i

frequencies Sym.

(cm-‘)

and intensities

MP2/ TZP B

(km/mol)

for the most stable conformer

BLYP/ DZ

of glycine

BP/ STO-DZP

B3LYP/ DZP

B3LYP/ TZZP

Vi

int.

*i

int.

v,

Vi

int.

Vi

int.

I 2 3 4 5 6 7 8 9 10 11 12 13 14 15

A’

3808 3562 3106 1825 1680 1474 1422 1318 1191 1147 949 845 642 472 260

75.9 2.8 13.5 264.7 22.5 14.4 22.7 19.1 56.2 264.1 124.5 76.8 5.9 31.5 9.4

3602 3396 2988 1747 1633 1416 1343 1263 1119 1072 905 790 604 442 247

29.8 0.1 17.3 253.5 15.8 16.8 13.4 6.6 62.7 201.0 162.9 50.7 6.1 29.5 9.0

3853 3553 3036 1803 1694 1407 1391 1285 1243 1061 998 848 672 493 278

3767 3517 3068 1833 1677 1456 1405 1308 1169 1135 933 826 631 461 255

54.4 1.2 17.2 295.8 20.2 17.4 26.2 9.6 128.5 155.1 150.1 78.3 6.9 31.0 10.0

3751 3518 3060 1807 1685 1466 1398 1312 1159 1118 919 818 635 463 258

73.1 7.1 11.2 293.6 20.1 18.3 17.0 10.6 121.4 182.4 126.9 83.3 5.7 29.0 10.4

16 17 18 19 20 21 22 23 24

A”

3652 3158 1401 1192 927 626 497 237 54

9.3 4.8 0.3 0.9 1.0 85.5 53.7 49.7 3.4

3471 3396 1341 1145 892 629 486 224 69

0.7 0.1 0.2 0.3 3.2 88.2 44.0 47.0 4.8

3616 3103 1308 1109 808 585 422 213 44

3596 3112 1377 1179 920 648 503 221 70

3.2 9.8 0.2 0.6 4.5 99.3 43.9 50.7 5.6

3585 3094 1392 1192 919 659 516 211 75

11.3 1.5 0.0 1.5 4.4 94.8 33.2 50.3 4.1

“Ref.

[lo].

conformers and indicate that I is the most stable conformation and II is a saddle point between I and III. Analysis of the geometrical parameters and the good agreement with the rotational constants indicate that the ED data for the NC bond length and the NCC valence angle are incorrect and should be shorter and wider respectively.

Acknowledgement This work has been financially supported by Italian CNR, CT12 and by MURST. One of us (FL) thanks NSERC of Canada for a grant received during a summer stay at the Department of Chemistry of University of Calgary (Alberta) and Professor T. Ziegler for hospitality and stimulating discussions. We thank Professor E.J. Baerends for allowing us to use

the latest version of the ADF program and CISIT of University of Basilicata for technical support.

References [I] R.D.Brown,P.D.Godfrey,J.W.V. StoreyandM.P.Bassez, J. Chem. Sot. Chem. Commun. (1978) 547. [2] R.D. Suenram and F. Lovas, J. Mol. Spectry. 72 ( 1978) 372. [3] H. Sellers and L. Schafer, J. Am. Chem. Sot. 100 (1978) 7728. [4] L. Schafer, H. Seller, F.J. Lovas and R.D. Suenram, J. Am. Chem. Sot. 102 (1980) 6656. [ 5 ] K. Ijima, K. Tanaka and S. Onuma, J. Mol. Struct. 246 (1991) 257. [6] M. Ramek, VKW. Cheng, R.F. Frey, S.Q. Newton and L. Schafer, J. Mol. Struct. 254 ( 199 1) 1. [7] A.G. Csaszar, J. Am. Chem. Sot. 114 (1992) 9568. [ 81 C.H. Hu, M.Z. Shen and H.F. Schaefer III, J. Am. Chem. Sot. 115 (1993) 2923.

F. Lefi et al. / Chemical Physics Letters 230 (1994) 189-195 [9] R.D. Suenram and F.J. Lovas, J. Am. Chem. Sot. 100 (1980) 7180. [lo] S. Vishveshwara and J.A. Pople, J. Am. Chem. Sot. 99 (1977) 2422. [ 111 E.S. Kryachko and E.V. Ludena, Density functional theory of many electron systems (Kluwer, Dordrecht, 199 1). [ 121 F. Herman, J.P. van Dyke and I.P. Ortenburger, Phys. Rev. Letters 22 (1969) 807. [ 131J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [ 141 C. Lee, W. Yang and R.G. Parr, Phys. Rev. B 37 (1988) 785. [ 151 A.D. Becke, Phys. Rev. A 38 (1988) 3098. [16]T.Ziegler,Chem.Rev.91 (1991)651. [ 171 C. Adamo and F. Lelj, Chem. Phys. Letters 223 (1994) 54. [ 181 F. Lelj and C. Adamo, submitted for publication. [ 191 F. Sim, A. St-Amant, I. Papai and D.R. Salahub, J. Am. Chem. Sot. 114 (1992) 4391. [20] B.G. Johnson, P.M.W. Gill and J.A. Pople, J. Chem. Phys. 98 (1993) 5612. [ 211 C. Filippi, C.J. Umrigar and M. Taut, J. Chem. Phys. 100 (1994) 1290. [ 221 A. Berces and T. Ziegler, J. Can. Chem., in press. [23] A.D. Becke, J. Chem. Phys. 98 (1993) 1372. [24] A.D. Becke, J. Chem. Phys. 98 (1993) 5648.

195

[25] V. Barone, Chem. Phys. Letters 226 (1994) 392; J. Chem. Phys. 101 (1994), in press. [ 261 C. Adamo, V. Barone and F. Lelj, unpublished results. [27] E.J. Baerends and P. Ros, Intern. J. Quantum Chem. S 12 (1978) 169. [28] M.J. Frisch, G.W. Trucks, M. Head-Gordon, P.M.W. Gill, M.W. Wong, J.B. Foresman, B.G. Johnson, H.B. Schlegel, M.A. Robb, E.S. Replogle, R. Gomperts, J.L. Andres, K. Raghavachari, J.S. Binkley, C. Gonzalez, R.L. Martin, D.J. Fox, D.J. DeFrees, J. Baker, J.J.P. Stewart and J.A. Pople, GAUSSIAN 92/DFf (Gaussian Inc., Pittsburgh, 1993). ~29I J.G. Snijders, P. Vemooijs and E.J. Baerends, At. Data Nucl. Data Tables 26 ( 1982) 483; P. Vemooijs, J.G. Snijders and E.J. Baerends, Internal Report ( Vrije Universiteit, Amsterdam, 198 1). [30 G. te Velde and E.J. Baerends, J. Comput. Phys. 99 (1992) 84. ]31 S.H. Vosko, L. Wilk and M. Nusair, J. Can. Phys. 58 ( 1980) 55. [32] T.H. Dunning and J.R. Hay, in: Modem theoretical chemistry, Vol. 3 (Plenum Press, New York, 1977) pp. l[33] PH. Dunning Jr., J. Chem. Phys. 72 (1980) 5639. [34] T.H. Dunning Jr., J. Chem. Phys. 90 (1989) 1007. [ 351 V. Barone, C. Adamo and F. Lelj, J. Chem. Phys., in press.