6310
J. Phys. Chem. 1996, 100, 6310-6316
Comparison of NMR Shieldings Calculated from Hartree-Fock and Density Functional Wave Functions Using Gauge-Including Atomic Orbitals Guntram Rauhut, Steve Puyear, Krzysztof Wolinski,* and Peter Pulay* Department of Chemistry and Biochemistry, The UniVersity of Arkansas, FayetteVille, Arkansas 72701 ReceiVed: October 3, 1995; In Final Form: February 7, 1996X
A numerically accurate implementation of the gauge-including atomic orbital method for the calculation of NMR shieldings in density functional theory (DFT) is presented. Results calculated by this method are compared with results of SCF and accurate coupled cluster calculations for eight small molecules. Three sets of DFT results, obtained using different exchange-correlation functionals, are further compared with each other, with the SCF data, and with experiment for a set of 10 somewhat larger organic molecules. The DFT values show a modest improvement compared to SCF theory. Potential computational savings using density functional theory are discussed.
Introduction The a priori calculation of NMR shieldings has received much interest recently.1,2 Due to the difficulty known as the gauge problem,3 such calculations are less straightforward than the calculation of energies. The most straightforward solution is the use of gauge-including atomic orbitals (GIAOs);4,5 the use of these in the SCF theory of NMR shieldings was pioneered by Ditchfield6 in the early seventies. The high computational cost of these early calculations precluded their wide use until 1980 when Kutzelnigg and co-workers7 developed the independent gauge for localized orbitals (IGLO) method. This method uses gauge factors on localized orbitals rather than on atomic basis functions, resulting in computational savings. A method similar in spirit, LORG, was devised by Hansen and Bouman.8 Subsequently, it was shown9 that the original GIAO method can be formulated very efficiently using the tools of modern analytic derivative theory;10 the latest programs generally prefer GIAOs to other formulations. GIAO SCF calculations of NMR shieldings suffer from the neglect of electron correlation; they are also more expensive than the calculation of the energy by a factor of 2-3. Several correlated methods have been developed for NMR shieldings but the most widely applicable ones are based on Møller-Plesset perturbation theory11 and coupled cluster theory.12 With large basis sets, these methods yield essentially exact results for most systems, although their computational cost is high and thus they are restricted to small systems. A cost-efficient alternative to traditional correlation theories is density functional theory13 (DFT). Earlier versions of this theory were not sufficiently accurate for most molecular problems, although they functioned well for solid-state problems for which the theory has been originally introduced. This has changed with the introduction of gradient-corrected exchangecorrelation potentials14,15 which were recently shown16,17 to be fully competitive with traditional correlation methods at substantially lower cost. Although these methods are essentially semiempirical and thus do not converge toward the exact results, they perform admirably in practice for most ground-state systems. A further advantage of DFT is that it can be significantly more efficient than even traditional SCF theory. Many workers in the DFT field have incorrectly overemphasized * Authors for correspondence. X Abstract published in AdVance ACS Abstracts, March 15, 1996.
0022-3654/96/20100-6310$12.00/0
this point, claiming that SCF theory scales with the fourth power of the molecular size while DFT scales only as O(N3) or O(N2). Although these sweeping claims have been refuted,18 it is true that replacing the Hartree-Fock exchange with a local potential can lead to important computational savings in large molecules, due to the fact that the Coulomb potential lends itself much more readily to efficient approximations than the exchange contribution. As part of our work on ab initio NMR software, we developed a GIAO-based chemical shift program for density functional wave functions, based on our program system TX90. This program is highly accurate numerically, due to the fact that no secondary fitting procedures are used, and the accuracy of the single numerical integration step is carefully controlled. At this point we can use only the usual current-independent exchangecorrelation functionals, although it is clear that the full theory should contain current-dependent exchange-correlation terms.19 Originally we intended to publish our conventional (density only) DFT results only together with results obtained from the Vignale-Rasolt formulation19 of current density functional theory, although the results in Tables 1, 2, and 8 were reported at the 1994 International Congress of Quantum Chemistry in Prague. The reason for this was that NMR shieldings calculated by conventional DFT theory turn out to be only marginally superior to SCF values, particularly for strongly correlated molecules composed of electronegative first-row atoms which constitute a stringent test. In some cases, the correlation contribution obtained from DFT calculations has even the wrong sign, compared to accurate correlated results. This contrasts with the excellent performance of gradient-corrected DFT theory for other second-order properties, for instance, force constants. Two reasons (besides slower than expected progress in implementing current density functional theory) prompted us to publish the present results. First, our recent data have shown that, in contrast to strongly correlated systems such as F2 or N2, DFT magnetic shieldings for most organic molecules represent an improvement over the SCF results. More importantly, recent advances in calculating the Coulomb interaction in DFT theory, in particular the fast multipole method20 and related techniques opened the exciting possibility of performing accurate DFT calculations with hundreds of atoms and thousands of basis functions. These approximations are difficult to apply to SCF wave functions because of the troublesome exchange term. Combining the GIAO NMR calculations with one of these © 1996 American Chemical Society
NMR Chemical Shieldings
J. Phys. Chem., Vol. 100, No. 15, 1996 6311
TABLE 1: Calculated and Experimental Isotropic Shielding Constants for Small Molecules mol CH4
nucl 13C 1
CO HCN HF H2O NH3 N2 F2
H 13C 17O 13 C 1H 15N 1H 19F 17 O 1H 15N 1 H 15N 19F
CPHF-GIAO basis Ia basis IIb 197.74 31.95 -23.88 -92.59 73.40 29.51 -50.28 29.02 413.34 332.18 31.31 266.94 32.16 -111.80 -169.96
196.60 31.88 -20.29 -80.76 73.74 29.21 -44.97 28.84 415.54 328.68 31.08 263.21 31.95 -106.14 -171.52
DFT-GIAOc basis Ia basis IIb 191.53 32.02 -7.51 -75.54 77.64 29.68 -39.47 30.79 408.41 332.57 32.21 266.24 32.56 -77.60 -265.35
187.80 31.95 -12.27 -73.60 71.74 29.40 -43.47 30.56 410.87 326.37 31.91 259.42 32.26 -80.55 -277.14
Gaussd CCSD(T) 198.9 31.6 5.6 -52.9 86.3 29.0 -13.6 29.2 418.6 337.9 30.9 270.7 31.6 -58.1 -186.5
Malkine SOS-DFPT 191.9 31.2
Zieglerf DFT-GIAO 191.2 31.4 -9.3 -68.4 91.5 8.4
29.5 410.0 325.6 31.2 257.2 31.2 -55.5 -202.8
412.5 331.5 31.2 262.0 -72.9 -282.7
exptd 198.4 ( 0.9 2.8 ( 0.9 -36.7 ( 17.2
29.2 ( 0.5 419.7 ( 6 357.6 ( 17.2 273.3 ( 0.1 -59.6 ( 1.5 -192.8
a TZP basis.34 b Pentuple-ζ + three polarization basis (see text).35 c All DFT calculations use the B-LYP (Becke-Lee-Yang-Parr) exchangecorrelation potential. d Reference 12 and references therein. e Reference 22, basis III, loc.1, PW91 exchange-correlation potential. f Reference 24.
TABLE 2: Calculated and Experimental Shielding Anisotropies for Small Moleculesa CPHF-GIAO basis basis mol nucl I II CH4
13C 1H
CO
13C 17O
HCN 13C 1H 15N HF 1H 19F H2O 17O 1H NH3 15N 1H N2 15N F2 19F
DFT-GIAO Ziegler basis basis Gauss DFTI II CCSD(T) GIAO
0.0 0.0 0.0 0.0 0.0 9.85 9.84 8.58 8.74 10.1 442.61 437.31 421.34 428.52 401.0 754.91 737.44 728.97 726.49 694.6 306.60 306.10 301.99 310.88 288.5 14.28 14.64 14.20 14.60 15.1 584.54 576.88 569.65 575.98 528.9 22.45 22.88 20.32 20.77 22.7 103.02 99.65 110.32 106.47 94.3 40.43 56.11 30.31 53.84 46.5 19.65 20.04 17.87 18.38 20.2 -41.14 -39.05 -53.66 -50.96 -43.8 15.88 15.82 14.62 14.56 16.2 675.78 667.50 626.30 631.03 596.5 987.71 990.08 1130.81 1148.63 1011.7
expt
424.1 406.1 ( 1.4 718.7 676.1 ( 26 286.1 502.9 104.2 93.8 -48.1 20.0 601.3 1057.
a See footnotes to Table 1. The anisotropy is defined by subtracting the average of the two closer shielding components from the third component; for linear or axially symmetrical molecules, this is equal to the difference of the parallel and the perpendicular components.
highly efficient DFT methods under development would result in a code which is perhaps not much more accurate but substantially more efficient (possibly by orders of magnitude for large molecules) than current SCF programs. The first chemically significant density functional calculation of NMR shielding constants were performed by Malkin, et al.21 They developed several methods, using the IGLO formalism of Kutzelnigg.7 In the complete basis set limit, their first method, which they call the uncoupled DFT method, should be equivalent to the method presented here. Both the IGLO and the GIAO methods perform well for SCF wave functions, but there are reasons to prefer GIAO, particularly for smaller basis sets.9 Malkin et al.22 call their latest method the sum-overstates (SOS) DFT perturbation theory. They argue that the DFT energy denominators (essentially the elements of the electronic Hessian) are too small as a result of the local approximation to the exchange operator and add a correction term to the orbital energy differences. The main problem with this explanation is that it should apply to other second-order properties such as polarizabilities or force constants and not just to NMR shieldings. There is, however, no sign that DFT exaggerates electronic response properties. On the contrary, force constants are very accurate in gradient-corrected DFT theories without extra correction terms in the electronic Hessian. Another formal problem is that the sum-over-states approach is not invariant
against unitary mixing of the occupied (or the virtual) orbitals; other perturbation theories can be formulated in unitary invariant manner, although at the expense of nondiagonal terms in the zeroth-order Hamiltonian. Finally, strict gauge invariance is not maintained in these calculations, as pointed out by van Wu¨llen.23 For these reasons, it is perhaps best to regard the sum-over-states formulation as an approximation to the current density functional theory, the exact form of which is unknown at present. Computationally, the sum-over-states formulation requires O(N2) numerical integration steps over molecular orbitals. Malkin et al.22 do not comment on the algorithm used to calculate these quantities. By arranging the loops so the outer loop runs over the grid points, it is probably possible to implement this part efficiently, although it appears to be more expensive than current density functional theory, contrary to the opinion expressed in ref 22. Both methods should be substantially less demanding than traditional correlation methods, e.g., Møller-Plesset perturbation theory. The results obtained using the sum-over-states theory are very promising, in spite of what appears to us a lack of solid theoretical foundation. Schreckenbach and Ziegler24 have recently presented an implementation of GIAO NMR shieldings in the Amsterdam density functional package (ADF). Although this method is identical with ours in principle, there are major differences in implementation. First, ADF uses a Slater basis while we use Gaussians. Second, our program does not use fitting functions for the Coulomb and exchange functions. Rather, we calculate the Coulomb term exactly, using traditional two-electron integral evaluation, and calculate the exchange-correlation terms by the numerical integration technique of Becke.25 Third, we use different exchange-correlation functionals than Schreckenbach and Ziegler. The comparison of the results obtained by the two techniques can thus shed light on the sensitivity of the results to numerical approximations and the choice of the exchangecorrelation potential. Theory and Computational Method The theory of DFT NMR shieldings is very similar to the corresponding SCF theory, the only substantial difference being the replacement of the exchange terms by a local exchangecorrelation potential. As pointed out by Fukui26 and Malkin et al.,21 this results in substantial computational simplification since in DFT theory it is not necessary to solve the analogues of the coupled-perturbed Hartree-Fock equations. This follows from
6312 J. Phys. Chem., Vol. 100, No. 15, 1996
Rauhut et al.
TABLE 3: B3-LYP/6-31G* Calculated and Experimentala Geometries of Selected Molecules molecule
sym coordinate calc
exp
CsH CsC CsH CdO CsH CtN CsC CsH CfdO CfsO CfsH CmsO CmsHa CmsHs CRsO CRdCβ CβsC′β CRsH CβsH CsN CsC CsC CdO CsHa CsHs CsO CsHa CsHs
1.093 1.509 1.087 1.206 1.110 1.160 1.461 1.094 1.206 1.342 1.101 1.441 1.093 1.090 1.364 1.361 1.436 1.079 1.081 1.484 1.549 1.521 1.216 1.097 1.092 1.410 1.103 1.093
1.094 1.512 1.083 1.208 1.116 1.157 1.458 1.104 1.200 1.334 1.101 1.437 1.086 1.086 1.362 1.361 1.431 1.075 1.077 1.477 1.560 1.520 1.214 1.103 1.103 1.410 1.100 1.091
CRdN CRsCβ CβdCγ CRsH CβsH CγsH
1.339 1.396 1.394 1.089 1.086 1.087
1.338 1.394 1.392 1.086 1.082 1.081
methane cyclopropane
Td D3h
formaldehyde
D2h
acetonitrile
C3V
methyl formateb Cs
furan
C2V
azetidine
Cs
acetone
C2V
dimethyl ether
C2V
pyridine
C2V
coordinate
calc
exp
HsCsH
113.9 114.0
HsCsH
115.2 116.5
direction of the external magnetic field B, and b is the direction of the nuclear magnetic moment; the subscript n denotes the nucleus in question. D is the (unperturbed) density matrix, D ) 2CoCo+, h is the one-electron Hamiltonian matrix in the atomic orbital basis; Co denotes the occupied block of the orbital coefficient matrix. Introducing the GIAO basis functions as
χp(B) ) exp(-iB‚Rp × r/2)χp(0) OdCfsO HsCfsO CfsOsCm HasCmsHs HasCmsHa
126.0 108.9 115.2 110.7 109.1
125.9 109.3 114.8 110.7 110.7
CRsOsC′R OsCRdCβ CRdCβsC′β OsCRsH CβsC′βsH CsNsC
106.8 110.5 106.1 115.6 127.3 90.5
106.6 110.7 106.1 115.9 128.0 88.0
CsCsC OsCsC HasCsHs HasCsHR CsOsC OsCsHa OsCsHs HasCsHa HasCsHs NdCRsCβ CRdNsC′R CRsCβdCγ CβsCγsC′β NdCRsH CβsCRsH CRsCβsH CγdCβsH CβdCγsH
116.5 121.8 109.6 106.8 112.3 111.8 107.2 108.1 108.9 123.8 117.1 118.4 118.5 115.9 120.3 120.3 121.3 120.7
116.0 122.0 108.4 108.4 111.7 110.8 107.2 108.7 109.6 123.8 116.9 118.5 118.4 116.0 120.2 120.1 121.4 120.8
a Experimental data taken from the following: Methane: Tarrago, G.; Dangh-Nhu, M.; Poussigue, G. J. Mol. Spectrosc. 1974, 49, 322. Cyclopropane: Butcher, R. J.; Jones, W. J. J. Mol. Spectrosc. 1973, 47, 64. Formaldehyde: Takagi, K.; Oka, T. J. Phys. Soc. Jpn. 1963, 18, 1174. Acetonitrile: Costain, C. C. J. Chem. Phys. 1958, 29, 864. Methyl formate: Curl, R. F. J. Chem. Phys. 1959, 30, 1529. Furan: No¨sberger, P.; Bauder, A.; Gu¨nthard, Hs. H. Chem. Phys. 1973, 1, 418. Azetidine: Dorofeeva, O. V.; Mastryukov, V. S.; Vilkov, L. V.; Hargittai, I. J. Chem. Soc., Chem. Commun. 1973, 772. Acetone: Iijima, T. Bull. Chem. Soc. Jpn. 1972, 45, 3526. Dimethyl ether: Blukis, U.; Kasai, P. H.; Myers, R. J. J. Chem. Phys. 1963, 38, 2753. Pyridine: Sørensen, G. O.; Mahler, L.; Rastrup-Andersen, N. J. Mol. Struct. 1974, 20, 119. b The subscript f denotes the carboxylic carbon and the subscript m the carbon of the methyl group.
the pure imaginary nature of the magnetic perturbation. If the unperturbed molecular orbitals are real (which is always possible in a nonmagnetic ground state), then the first-order perturbed molecular orbitals are pure imaginary, and their contribution to the density vanishes in first order. This computational advantage does not apply to the more accurate current density functional theories, nor to exchange-correlation functionals which retain part of the Hartree-Fock exchange. We prefer the density matrix formulation described in ref 9. This gives the ab component of the magnetic shielding tensor as a sum of matrix traces: ab a 0b σab n ) Tr(Dhn ) + Tr(D hn )
(1)
Following the convention of Ditchfield and Ellis,27 the superscripts denote differentiation with respect to the external magnetic field and the nuclear magnetic moment, in this order. A single superscript denotes differentiation with respect to the external magnetic field. Thus the superscript a stands for the
(2)
we define |p〉 as the zero-field limit of χp(B), and |pa〉 as its derivative with respect to the field:
|pa〉 ) χap ) (-iRp × r/2)χp
(3)
where Rp denotes the center of basis function |p〉. In terms of these quantities, the matrix elements are given by a 0b a (hab ˆ n |q〉 + 〈p|hˆ 0b ˆ ab n )pq ) 〈p |h n |q 〉 + 〈p|h n |q〉
(4)
(h0b ˆ 0b n )pq ) 〈p|h n |q〉
(5)
The derivative operators are27 1 2 3 hˆ ab n ) ( /2c )[δab(r - Rn)‚r - (r - Rn)arb]/|r - Rn| (6) 2 3 hˆ 0b n ) (-i/c )[(r - Rn) × ∇]b/|r - Rn|
(7)
The first-order density matrix is obtained from
Da ) -1/2DSaD + ∑iv(i - v)-1Ci+(Fa - iSa)Cv[CiC+ v -
CvCi+] (8)
where Ci and Cv denote column vectors of the occupied and virtual orbital coefficients, respectively, Sa is the derivative of the overlap matrix with respect to Ba, and Fa is the corresponding derivative of the Fock matrix. The latter is given by
Fa ) ha + Ja + Vaxc
(9)
Here, in analogy to eq 4, the first-order matrix elements of the bare-nucleus Hamiltoniann are given by
(ha)pq ) 〈pa|hˆ |q〉 + 〈p|hˆ |qa〉 + 〈p|hˆ a|q〉
(10)
hˆ a ) -(i/2)r × ∇
(11)
The first-order contributions to the Coulomb operator matrix elements originate exclusively in the gauge factors:
(Ja)pq ) ∑rs(pq|rs)aDrs
(12)
Compared to the Hartree-Fock case, the only new quantity is the first-order exchange-correlation term a a 3 (Ka)pq ) ∫(χ* p χq + χ* pχq)Vxc(F,∇F) d r
(13)
where F is the electron density. We evaluate these terms by the same Becke numerical integration scheme25 we use for the exchange matrix elements themselves. Three different exchange-correlation potentials have been implemented: the simple Dirac-Slater-Wigner28,29 (DSW) exchange, Vxc(F) ) (3/2)(3/8π)1/3F1/3; the Becke15 gradient-corrected exchange functional (BXC), and the Becke exchange combined with the correlation functional of Lee, Yang, and Parr30,31 (BLYP). These three are representative of the functionals used in modern DFT theory. The most accurate of the three is the B-LYP functional, and we will concentrate on the results obtained with it; the others are included mainly to indicate the sensitivity of the results to the choice of the exchange-correlation
NMR Chemical Shieldings
J. Phys. Chem., Vol. 100, No. 15, 1996 6313
TABLE 4: Calculated and Experimental NMR Shieldings δ(13C) of Selected Molecules absolute values
relative to CH4
SCF
BLYPa
DSWb
BXCc
BLYPa
DSWb
BXCc
C
195.68
189.04
197.60
Methane 187.95 0.00
0.00
0.00
0.00
0.00e
C
197.73
181.85
187.83
Cyclopropane 181.39 -2.05
7.19
9.77
6.55
4.2f
C
-3.78
-16.53
-16.53
Formaldehyde -12.04 199.46
205.58
214.13
199.99
196.1g
CN CH3
192.25 64.42
180.54 65.12
187.85 63.76
Acetonitrile 179.94 3.43 66.52 131.26
8.50 123.92
9.75 133.84
8.01 121.43
8.81h 123.53h
COO CH3
19.98 143.52
17.24 127.57
17.42 133.09
Methyl Formate 19.27 175.70 127.26 52.16
171.81 61.69
180.18 64.51
168.68 60.69
168.37h 57.80h
CR Cβ
39.96 76.60
35.35 68.03
36.12 69.39
155.72 119.07
153.69 121.02
161.48 128.21
150.95 118.45
149.51h 116.42h
CR Cβ
147.84 171.25
126.16 155.17
130.71 161.64
Azetidine 125.25 47.83 154.09 24.43
62.88 33.87
66.89 35.96
62.69 33.86
54.2i 28.2i
CO CH3
-20.98 164.37
-29.42 150.24
-29.40 153.85
Acetone -27.08 216.66 150.32 31.30
218.46 38.81
226.99 43.75
215.03 37.63
213.58h 37.84h
C
138.32
117.03
119.58
Dimethyl Ether 116.62 57.36
72.01
78.01
71.33
64.9j,k
CR Cβ Cγ
30.30 65.08 44.79
26.31 54.34 43.66
25.78 54.13 43.92
Pyridine 27.62 165.38 55.40 130.60 44.78 150.89
162.73 134.70 145.38
171.82 143.47 153.68
160.33 132.55 143.17
156.79h 130.63h 142.8h
atom
SCF
expd
Furan 37.00 69.49
a Becke, Lee-Yang-Parr exchange-correlation functional. b Dirac, Slater, Wigner exchange-correlation functional. c Becke’s gradient corrected exchange-correlation functional. d Experimental values corrected to CH4(g) reference by using δ(13CH4) ) -7.0 ppm relative to TMS(l). e Absolute experimental value taken as 195.1 ppm from: Jameson, A. K.; Jameson, C. J. Chem. Phys. Lett. 1987, 134, 461. f Taken from: Breitmaier, E.; Voelter, W. Carbon-13 NMR Spectroscopy; Verlag-Chemie: Weinheim, Germany, 1987. Or: Kalinowski, H.-O.; Berger, S.; Braun, S. 13C NMR Spectroscopy; Georg Thieme Verlag: Stuttgart, Germany, 1984. g Chestnut, D. B.; Foley, C. K. J. Chem. Phys. 1986, 84, 852. h Pouchert, C. J.; Behnke, J. The Aldrich Library of 13C and 1H FT NMR Spectra, 1st ed.; Aldrich Chemical Company, Inc.: Milwaukee, WI, 1983. i Lambert, J. B.; Wharry, S. M.; Block, E.; Bazzi, A. A. J. Org. Chem. 1983, 48, 3982. j Experimental value corrected from reference CS2 to CH4 by 203.1 ppm. k Dorman, D. E.; Bauer, D.; Roberts, J. D. J. Org. Chem. 1975, 40, 3729.
functional. Empirical combinations of SCF exchange and density-correlation functions, for instance, the three-parameter functional proposed by Becke32 eliminate some of the deficiencies of the B-LYP method, e.g., the overestimation of the C-H bond lengths. However, these methods also increase the computational cost by requiring a coupled-perturbed HartreeFock step. The numerical grids used in the present study are probably larger than in most similar programs, and thus the results are largely free from the numerical noise which often plagues density functional calculations. Typically, the program uses over 4000 gridpoints for each hydrogen atom, and over 6000 gridpoints for each first-row atom (C, N, O, etc.). Results Tables 1 and 2 compare the Hartree-Fock and B-LYP isotropic shieldings and their anisotropies with the accurate coupled cluster data of Gauss and with experiment for a set of small molecules. Comparison is made with the latest CCSD(T)33 values; these do not, in general, differ much from the CCSD12 ones. We used two basis sets in the DFT calculations. Set I is the triple-ζ plus polarization (TZP) basis of Scha¨fer et al.34 This set is reasonably accurate for the NMR shieldings of first-row atoms and still economical to use for large molecules. Set II is the basis recommended by Augspurger and Dykstra.35 This is an 11s8p3d set, contracted to 7s7p3d for first-row atoms, and 6s2p, contracted to 4s2p for hydrogen. The latter basis is large enough so that the main difference between our results
and the CCSD(T) ones should be the treatment of electron correlation. To make the comparison with the coupled cluster results as close as possible, we use the same experimental geometries as Gauss and Stanton.12,33 On the basis of the convergence of the results, CCSD(T) includes a very large fraction of the dynamical correlation effect and most of the effect of nondynamical correlation. As this table shows, the DFT results are, on the average, somewhat better than the SCF results. They are, however, still quite far from reaching the accuracy of the coupled cluster results. Despite the different numerical procedure and the fact that we use the Lee-Yang-Parr30,31 correlation potential and Schreckenbach and Ziegler24 use the Perdew14 potential, the two sets of results also agree reasonably well, with the exception of HCN, which we cannot explain. The geometries used by Schreckenbach and Ziegler24 are specified only as “experimental” and thus may differ from ours. Comparing the DFT and CCSD(T) values in Tables 1 and 2 it is evident that the DFT results are less shielded, i.e., the (negative) paramagnetic contributions are too large in absolute value. This is evident from both the isotropic shieldings and the anisotropies; for the mostly line molecules in these tables, the parallel (z) component of the paramagnetic shielding vanishes, and thus the overestimation of the paramagnetic contribution by DFT leads to anisotropies which are too high. For hydrogens, though, the trend is opposite: the DFT shieldings are slightly too high and the anisotropies a little too low. The general overestimation of the paramagnetic contributions by DFT lends some justification to the procedure of Malkin et al.:22
6314 J. Phys. Chem., Vol. 100, No. 15, 1996
Rauhut et al.
TABLE 5: Calculated and Experimental NMR Shieldings δ(1H) of Selected Molecules absolute values atom SCF
BLYPa
DSWb
relative to CH4 BXCc
SCF BLYPa DSWb BXCc
expd
H
Methane 31.60 31.68 31.25 31.81 0.00
H
Cyclopropane 31.91 31.63 31.00 31.84 -0.31 0.05
0.25 -0.03 0.09e
H
Formaldehyde 22.92 21.66 21.66 21.77 8.68 10.02
9.58
H
Acetonitrile 30.35 30.17 29.62 30.35 1.25 1.51
1.63
1.46 1.88g
0.00
0.00
0.00 0.00
10.04 12.3 ( 2f
Formateh
Methyl HCO 24.50 23.82 23.13 23.93 7.10 H3C 28.56 28.17 27.51 28.35 3.04
7.86 3.51
8.12 3.74
7.88 7.95g 3.46 3.62g
Furan 7.08 6.02
7.26 6.26
7.56 6.65
7.22 7.29g 6.18 6.25g
HCR 24.52 24.42 23.69 24.59 HCβ 25.58 25.42 24.60 25.63 HCR HCR HCβ HCβ HN
28.76 29.03 30.50 29.77 31.44
Azetidine 28.05 2.84 28.47 2.57 30.20 1.10 29.58 1.83 31.43 0.16
3.78 3.34 1.64 2.29 0.46
4.21 3.71 1.89 2.61 0.55
3.76 3.34 1.61 2.23 0.38
H
Acetoneh 30.13 29.83 29.15 30.02 1.47
1.85
2.10
1.79 2.04g
27.90 28.34 30.04 29.39 31.22
27.04 27.54 29.36 28.64 30.70
3.50i 2.89i 1.65i 2.13i 0.91j,k
Etherh
H
Dimethyl 29.11 28.47 27.72 28.62 2.49
3.21
3.53
3.19 3.11e
HR Hβ Hγ
Pyridine 23.12 23.09 22.21 23.27 8.48 24.78 24.70 23.83 24.92 6.82 24.07 24.27 23.41 24.47 7.53
8.59 6.98 7.41
9.04 7.42 7.84
8.54 8.47g 6.89 7.12g 7.34 7.52g
a Becke, Lee-Yang-Parr exchange-correlation functional. b Dirac, Slater, Wigner exchange-correlation functional. c Becke’s gradient corrected exchange-correlation functional. d Experimental values corrected to CH4(g) reference by using δ(C1H4) ) +0.13 ppm relative to TMS. e Emsley, J. W.; Feeney, J.; Sutcliffe, L. H. High Resolution Nuclear Magnetic Resonance Spectroscopy; Pergamon Press: Oxford, 1966; Vol. 2. f Kukolich, S. G. J. Am. Chem. Soc. 1975, 97, 5704. g Estimated from the spectrum given in: Pouchert, C. J.; Behnke, J. The Aldrich Library of 13C and 1H FT NMR Spectra, 1st ed.; Aldrich Chemical Company, Inc.: Milwaukee, WI, 1983. h Calculated values are weighted average values of the different shifts of the hydrogen atoms of a methyl group. i Dutler, R.; Rauk, A.; Sorensen, T. S. J. Am. Chem. Soc. 1987, 109, 3224. j Chauvel, J. P.; Folkendt, M. M.; True, N. S. Magn. Reson. Chem. 1987, 25, 101. k The authors of footnote i report an approximate value of 1.47 ppm (after conversion to the CH4 reference).
by increasing the energy denominators they decrease the paramagnetic contributions. The changes upon increasing the basis from I to II are modest in Tables 1 and 2, justifying our decision to use only basis I for the larger molecules in Tables 4-7. The F2 results in Tables 1 and 2 deserve comment. Although F2 is an exceptionally strongly correlated system, almost a diradical, the correlation contribution to the NMR shieldings is small, probably because of fortuitous cancellation. The DFT results are again too little shielded compared to the CCSD(T) and experimental data. Comparison with experiment is made difficult by the very strong geometry dependence of the shielding, resulting in a large (∼40 ppm) vibrational correction.37 We have calculated isotropic NMR shieldings and chemical shifts of 13C, 1H 15N, and 17O for a set of somewhat larger molecules. Their geometries have been fully optimized using the Becke three-parameter compound exchange-correlation functional37 as implemented in Gaussian92.38 This method is known to give very accurate geometry parameters for molecules composed of first-row atoms.39 The geometries are summarized, as far as practical, in Table 3. Tables 4-6 contain the NMR
quantities calculated at the SCF level, as well as with three different exchange-correlation potentials: B-LYP,15,30,31 the Becke gradient corrected exchange15 (BXC) and the simple Dirac-Slater-Wigner potential28,29 (DSW). Coupled cluster results are not available for these systems, and thus we have compared the calculated values with experimental data, although this is a less accurate comparison because of experimental errors (e.g., bulk susceptibility corrections), geometry effects, vibrational corrections, and solvent effects. All results in Tables 4-6 were obtained using the triple-ζ basis set of Scha¨fer et al. with one set of polarization functions. This is adequate for 1H and 13C but the more sensitive 15N and 17O shieldings may benefit from the addition of more polarization functions. On the average, agreement with experiment is slightly better for the DFT results than for SCF. The B-LYP and the Becke exchange results differ from each other by only a few ppm for 13C and by less than 0.1 ppm for 1H but the Dirac-Slater-Wigner results differ significantly. In general, the absolute shieldings calculated by density functional theory are lower than the SCF values for 13C and 1H. The simplest exchange-correlation potential, DSW (the Dirac-Slater-Wigner exchange), is worse than SCF. B-LYP15,30,31 and BXC (the Becke gradient-corrected exchange15) are comparable and both perform somewhat better than SCF. From Table 4, the average absolute deviation in 13C shieldings is 5.54 ppm for SCF, 4.32 for B-LYP, 2.74 for BXC, and 10.55 for DSW. However, these numbers depend quite sensitively on the absolute shieldings of the reference, in our case gaseous 13CH4; using liquid methane as reference gives somewhat different values. The situation is similar for 1H, where the average absolute deviations are 0.62 ppm (SCF), 0.28 (B-LYP), 0.32 (BXC), and 0.49 (DSW). On the whole, despite the slight improvement provided by density functional theory over SCF, both are of similar quality and probably cannot compete with MP2, unlike for other properties such as force constants. Table 7 contains the full calculated tensor values for pyridine as an example for the prediction of tensor components for larger molecules. Finally, Table 8 compares SCF, DFT (B-LYP), and MP2 results for two difficult triatomic molecules: N2O and CO2. The same large (quadruple-ζ with two sets of polarization functions) basis set was used for all calculations, except that for technical reasons we used five-component d functions for the DFT calculations and six-component ones for MP2. The MP2 calculations were carried out by the Aces-2 program system.40 SCF shieldings were determined both with five- and sixcomponent d functions; as Table 8 shows, the differences are negligible. CO2 is remarkable in that density functional theory gives a correlation contribution which is of the incorrect sign as compared to the MP2 results, for both the carbon and oxygen nuclei. The problem appears to be, as in the other cases, the overestimation of the paramagnetic contribution. The results are somewhat better for N2O but the shieldings show only a modest improvements over the SCF values. Timing data for methyl formate (TZP basis, no symmetry used) are listed in Table 9 to give the reader an idea about the computing costs involved. In our present implementation, the calculation of the exchange-correlation potential makes the SCF step significantly slower than the Hartree-Fock method, despite cutoffs being used, at least for smaller molecules. However, the GIAO-NMR step is faster, due to the absence of the coupledperturbed iteration in the DFT methods, resulting in an overall timing close to that of the Hartree-Fock method. Conclusions Density functional methods are useful for the calculation of NMR shieldings. The results improve somewhat relative to the
NMR Chemical Shieldings
J. Phys. Chem., Vol. 100, No. 15, 1996 6315
TABLE 6: Calculated and Experimental NMR Shieldings δ(15N) and δ(17O) of Selected Molecules absolute values atom
SCF
O
326.23
BLYPa 326.15
relative to H2O
DSWb
BXCc
337.57
Water 324.26
SCF
BLYPa
DSWb
BXCc
expd
0.00
0.00
0.00
0.00
0.00
769.28
780.70
762.81
719 ( 150e
O
-461.33
-443.13
-443.13
Formaldehyde -438.55 787.33
OCH3 CO
-117.40 169.53
-113.96 108.85
-123.17 109.20
Methyl Formate -112.75 443.63 108.64 156.70
440.11 217.30
460.74 228.37
437.01 215.62
400f 179f
O
52.69
14.28
7.13
Furan 13.27 274.00
311.87
330.44
310.99
276f
O
-357.64
-354.31
-379.67
Acetone -352.07 683.87
680.46
717.24
676.33
605g
O
349.86
307.16
321.72
Dimethyl Ether 302.21 -23.63
18.99
15.85
22.05
atom
SCF
BLYPa
DSWb
BXCc
BLYPa
DSWb
BXCc
264.33
263.24
274.89
Ammonia 262.10
0.00
0.00
0.00
0.00
0.00
294.09
310.46
290.53
272.6j,k
absolute values
N
16h,i
relative to NH3 SCF
exp
N
-44.85
-30.85
-35.57
Acetonitrile -28.43 309.18
N
236.53
209.54
219.59
Azetidine 206.41 27.80
53.70
55.30
55.69
25.3l
N
-107.51
-97.26
-102.28
Pyridine -95.96 371.84
360.50
377.17
358.06
316.5g
a Becke, Lee-Yang-Parr exchange-correlation functional. b Dirac, Slater, Wigner exchange-correlation functional. c Becke’s gradient corrected exchange-correlation functional. d Experimental 17O values corrected to H2O(g) reference by using δ(H217O(g)) ) -36.1 ppm relative to H2O(l). Jameson, C. J.; Jameson, A. K.; Burell, P. M. J. Chem. Phys. 1980, 73, 6013. e Neuman, D. B.; Moscowitz, J. W. J. Chem. Phys. 1969, 50, 2216. f Sugawara, T.; Kawada, Y.; Katoh, M.; Iwamura, M. Bull. Chem. Soc. Jpn. 1979, 52, 3391. g Mason, J. Multinuclear NMR; Plenum Press: New York, 1987. h Delseth, C.; Kintzinger, J.-P. HelV. Chim. Acta 1978, 61, 1327. i The authors of footnote f report a value (after conversion to H2O(g) as reference) of 6 ppm. j Jameson, C. J.; Jameson, A. K.; Oppusunggu, D.; Wille, S.; Burrell, P. M.; Mason, J. J. Chem. Phys. 1981, 74, 81. k Another gas-phase value of 242.5 ppm is reported by the author of footnote g. l Lambert, J. B.; Wharry, S. M.; Block, E.; Bazzi, A. A. J. Org. Chem. 1983, 48, 3982.
TABLE 7: Principal Components of the NMR Shielding Tensors in Pyridine at Different Computational Levels atom
SCF
BLYP
DSW
BXC
N
-426.07 -205.94 309.47 -99.07 27.81 162.18 -57.35 66.04 186.55 -84.24 26.70 191.91 19.50 22.71 27.26 21.75 24.46 28.11 20.98 23.78 27.44
-399.52 -180.65 288.40 -84.07 24.46 138.54 -58.35 53.54 167.83 -69.90 26.45 174.42 19.74 22.27 27.24 22.63 23.39 28.06 21.71 23.49 27.61
-421.60 -184.91 299.69 -85.39 20.22 142.51 -61.74 51.57 172.56 -72.16 23.10 180.83 18.68 21.41 26.52 21.78 22.39 27.30 20.74 22.63 26.84
-393.70 -179.40 285.24 -81.84 27.68 137.01 -56.41 56.60 166.01 -68.01 29.72 172.63 19.99 22.40 27.40 22.96 23.52 28.27 22.01 23.60 27.80
CR Cβ Cγ HR Hβ Hγ
Hartree-Fock method but the paramagnetic terms are overestimated without the inclusion of current-dependent terms in the exchange-correlation functional, resulting, in general, shielding values which are too low, i.e., too deshielded. DFT offers the potential for much increased efficiency for large molecules. Acknowledgment. This project has been supported by the National Science Foundation under Grant CHE-9319929 and
TABLE 8: Comparison of SCF, DFT (B-LYP) and MP2 NMR Shieldings and Shielding Anisotropies for CO2 and N2Oa isotropic shielding atom
SCF, d6
SCF, d5
DFT, d5
anisotropy MP2, d6
SCF, d6
SCF, d5
DFT, d5
53.22 53.16 51.55 67.05 344.86 344.95 352.65 CO2 CO2 223.87 223.88 213.98 246.64 286.70 286.70 301.38 NNO 64.57 64.36 89.84 134.36 416.80 417.12 381.63 NNO -29.88 -29.96 -0.24 35.94 569.72 569.85 529.91 NNO 173.77 173.70 173.43 219.06 366.61 366.74 365.29
MP2, d6 328.76 251.90 314.63 475.31 296.12
a Scha ¨ fer34 et al. quadruple-ζ basis with two polarization functions (10s6p2d/6s4p2d) for first-row atoms, (5s2p/4s2p) for hydrogen. The following geometry was used: r(CO) ) 1.15979 Å, r(NN) ) 1.1284 Å, r(NO) ) 1.1847 Å.
TABLE 9: Timing Data for Methyl Formate TZP Basis, No Symmetry (Minutes on an RS6000/390 Workstation), Traditional (Integral Nondirect) Method method
integrals
SCF
GIAO NMR
total
Hartree-Fock DFT B-LYP DFT DSW
3.0 3.0 3.0
1.9 17.9 10.3
17.2 11.7 11.1
22.0 32.6 24.4
the Air Force Office of Scientific Research under Grant No. F49620-94-1. G.R. thanks the Deutsche Forschungsgemeinschaft for financial support. We also thank the National Science Foundation and IBM Co. for an equipment grant to P.P. and J. F. Hinton. We are grateful to Prof. J. Gauss for a preprint of the CCSD(T) results, and to him and Prof. R. J. Bartlett for their permission to use the MP2 NMR shielding code implemented in the ACES-2 program package. We thank Dr. S.
6316 J. Phys. Chem., Vol. 100, No. 15, 1996 Kristyan for writing parts of the density functional wave function code and Mr. G. Magyarfalvi for calculating the values in Table 8. References and Notes (1) Kutzelnigg, W.; Fleischer, U.; Schindler, M. In NMRsBasic Principles and Progress; Springer: Berlin, 1990; Vol. 23, p 165. (2) Jameson, C. J. Nucl. Magn. Reson. 1992, 21, 36; 1993, 22, 59. (3) For example: Pulay, P.; Hinton, J. F. In Encyclopedia of NMR; Farrar, T. C., Ed. Wiley: Chichester, Vol. 2, in press. (4) London, F. J. Phys. Radium 1937, 8, 397. (5) Hameka, H. F. Mol. Phys. 1958, 1, 203. (6) Ditchfield, R. Mol. Phys. 1974, 27, 789. (7) Kutzelnigg, W. Isr. J. Chem. 1980, 19, 193. (b) Schindler, M.; Kutzelnigg, W. J. Chem. Phys. 1982, 76, 1919. (8) Hansen, A. E.; Bouman, T. D. J. Chem. Phys. 1985, 82, 5035. (9) Wolinski, K.; Hinton, J. F.; Pulay, P. J. Am. Chem. Soc. 1990, 112, 8251. (10) For example: Pulay, P. AdV. Chem. Phys. 1987, 69, 241. (11) Gauss, J. J. Chem. Phys. 1993, 99, 3629. (12) Gauss, J.; Stanton, J. F. J. Chem. Phys. 1995, 102, 251. (b) Gauss, J. Chem. Phys. Lett. 1994, 229, 198. (13) For example: Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (14) Perdew, J. Phys. ReV. B 1986, 33, 8822; 1986, 34, 7406. (15) Becke, A. D. J. Chem. Phys. 1988, 88, 2547. (16) Gill, P. M. W.; Johnson, B. G.; Pople, J. A.; Frisch, M. J. Chem. Phys. Lett. 1992, 197, 499. (17) Handy, N. C.; Maslen, P. E.; Amos, R. D.; Andrews, J. S. Chem. Phys. Lett. 1992, 197, 506. (18) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. J. Chem. Phys. 1993, 98, 5612. (19) Vignale, G.; Rasolt, M. Phys. ReV. B 1988, 37, 10685.
Rauhut et al. (20) White, C. A.; Head-Gordon, M. J. Chem. Phys. 1994, 101, 6593. (21) Malkin, V. G.; Malkina, O. M.; Salahub, D. R. Chem. Phys. Lett. 1993, 204, 80, 87. (22) Malkin, V. G.; Malkina, O. L.; Casida, M. E.; Salahub, D. R. J. Am. Chem. Soc. 1994, 116, 5898. (23) Van Wu¨llen, C. J. Chem. Phys. 1995, 102, 2806. (24) Schreckenbach, G.; Ziegler, T. J. Phys. Chem. 1995, 99, 606. (25) Becke, A. D. J. Chem. Phys. 1988, 88, 2547. (26) Fukui, H. Magn. Reson. ReV. 1987, 11, 205. (27) Ditchfield, R.; Ellis, P. D. In Topics in Carbon-13 NMR Spectroscopy; Levy, G. C., Ed.; Wiley: New York, 1974; Vol. 1, p 1. (28) Dirac, P. A. M. Proc. Cambridge Philos. Soc. 1930, 26, 376. (29) Wigner, E. P. Trans. Faraday Soc. 1938, 34, 678. (30) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (31) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Letters 1989, 157, 200. (32) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (33) Gauss, J. J. Chem. Phys., in press. (34) Scha¨fer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571. (35) Augspurger, J. D.; Dykstra, C. E. J. Phys. Chem. 1991, 95, 9230. (36) Jameson, C J.; Jameson, A. K.; Burrell, P. M. J. Chem. Phys. 1980, 73, 6013. (37) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Wong, M. W.; Foresman, J. B.; Robb, M. A.; Head-Gordon, M.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzales, C.; Martin, P. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian92/DFT; Gaussian Inc.: Pittsburgh, 1993. (39) For example: Rauhut, G.; Pulay, P. J. Phys. Chem. 1995, 99, 3093. (40) Stanton, J. F.; Gauss, J. D.; Watts, J. D.; Lauderdale, W. J.; Bartlett, R. J. Int. J. Quantum Chem. Symp. 1992, 26, 879.
JP9529127
