J . Am. Chem. SOC.

1990, 112, 8251-8260 825 1

Efficient Implementation of the Gauge-Independent Atomic

Orbital Method for NMR Chemical Shift Calculations
Krzysztof Wolinski,+ James F. Hinton, and Peter Pulay*
Contributionfrom the Department of Chemistry and Biochemistry, University of Arkansas,
Fayetteville, Arkansas 72701. Received March 16, 1990

Abstract: An implementation of the gauge independent atomic orbital (CIAO) method for the ab initio self-consistent-field
(SCF) calculation of nuclear magnetic resonance (NMR) chemical shifts is described. Using modern techniques borrowed
from analytical derivative methods, we were able to improve the efficiency of the CIAO method significantly. Results with
several basis sets, some of them large, are presented for methane, methyl fluoride, cyclopropene, cyclopropane, oxirane, benzene,
carbon disulfide, the sulfate and thiosulfate anions, dimethyl sulfide, dimethyl sulfoxide, and dimethyl sulfone. Computer
timings for energy and chemical shielding calculations are given for a few large organic molecules. Comparisons are made
with the individual gauge for localized orbitals (IGLO) method of Schindler and Kutzelnigg, and with the localized orbital/lccal
origin (LORG)method of Hansen and Bouman. The GIAO method appears to converge faster than the localized techniques;
i.e., it provides the same accuracy with a smaller basis, particularly for the individual tensor components. The computational
effort for the ab initio calculation of the NMR chemical shifts is only -2.5 times of the energy calculation.

1. Introduction ner-Prettre et aL8 and by Fukui et aL9 The CIAO method has
Modern high-field, multipulse N M R spectroscopy has proven been successfully used in calculations for small- and medium-size
to be an exceptionally powerful technique for solution of many systems.I0
types of problems in chemistry and biochemistry. However, the In spite of these successes, it appears that the existing imple-
problem of correct signal assignment as well as the understanding mentations of the CIAO method require too much in computing
of the relationship between the chemical shift and molecular resources to be routinely applicable to large molecules. Chesnut
structure can be quite difficult. Ab initio calculations are now et a1.I’ have recently suggested a way to improve the efficiency
becoming affordable and accurate enough to be useful in the of CIAO calculations by using “attenuated” basis sets, which
solution of some of these problems. A comparison of the ex- describe well the nuclei under consideration and less accurately
perimental and theoretical spectra can be very useful in making the others. In our opinion this method may be inefficient when
correct assignments and understanding the basic chemical chemical shifts for a number of nuclei are needed, e.g., in C16H24
shift-molecular structure relationship. discussed later in this paper. Other, more promising methods
Calculations of magnetic shielding tensors have been performed improve the efficiency of the CIAO calculations by applying the
at various theoretical levels. Semiempirical methods’ provide a gauge factors to localized molecular orbitals instead of every
correct qualitative understanding but are not accurate quantita- atomic orbital. Two such methods have been developed recently:
tively. Almost all a b initio results are at the S C F level, and are the individual gauge localized orbital (IGLO) method by Schindler
thus based on the coupled Hartree-Fock perturbation theory and KutzelniggIz and the localized orbital/local origin (LORG)
(CHF).2 A common difficulty in the calculation of magnetic. method by Hansen and B o ~ m a n . ’ In ~ practice, both methods
properties is that the usual wave functions do not guarantee gauge
invariance? i.e., in the simplest case, the results may depend on ( I ) Fukui, H. Magn. Reson. Reo. 1987, 11, 205.
the position of the molecule in the Cartesian frame. The physical (2) Stevens, R. M.; Pitzer, R. M.; Lipscomb, W. N. J. Chem. Phys. 1963,
reason for this is that a magnetic field, say in the z direction, leads 38, 550.
to a perturbation in the momenta, which twists the molecular (3) Hameka, H. F. Aduanced Quantum Chemistry; Addision-Wesley:
New York, 1963; p 162.
orbitals around the z axis in the imaginary xy plane. If an atom (4) Epstein, S.T. J. Chem. Phys. 1965, 42, 2897.
is situated on the z axis, its basis set is closed under this rotation, (5) (a) Holler, R.; Lischka, H. Mol. Phys. 1980, 41, 1017, 1041. (b)
and it can be described equally well both unperturbed and in the Galasso, V. Theor. Chim. Acta 1983,63, 35. (c) Galasso, V. J . Mol. Struct.
presence of the perturbation. However, for an atom far from the 1983. 93, 201. (d) Jaszunski, M.: Adamowicz, L. Chem. Phys.Lett. 1981, 79,
133. (e) Lazzeretti, P.; Zanasi, R. J . Chem. Phys. 1980, 72, 6768. (0
z axis, the rotation of the orbitals can only be described by using Lazzeretti, P.; Zanasi, R. J . Chem. Phys. 1981, 75, 5019. (9) Lazzeretti, P.;
high angular momentum basis functions, which are not normally Zanasi, R. J. Am. Chem. SOC.1983, 105, 12. (h) Tossell, J. A,; Lazzeretti,
included in the basis. P.; Vaughan, D. J. J . Magn. Reson. 1987, 73, 334. (i) Lazzeretti, P.; Tossell,
There are two general approaches to solve the gauge problem. J . A. J. Phys. Chem. 1987, 91. 800.
(6) London, F. J. Phys. Radium 1937, 8, 397.
Very large basis sets lead to approximate gauge invariance since (7) Ditchfield, R. Mol. Phys. 1974, 27, 789.
in the limit of complete basis sets the results should be gauge (8) Ribas Prado, F.; Giessner-Prettre, C.; Daudey, J.-P.; Pullman, A,;
independent! A number of calculations have been published by Young, F.; Hinton, J. F.; Harpool, D. J. Magn. Reson. 1980, 37, 431.
this method for small molecule^.^-^ In our opinion, however, this (9) Fukui, H.; Miura, K.; Yamazaki, H.; Nosaka, T. J . Chem. Phys. 1985,
82, 1410.
is not a viable technique for larger systems. The second approach (IO) (a) Ditchfield, R. J. Chem. Phys. 1976.65, 3123. (b) Ditchfield, R.
introduces local gauge origins to define the vector potential of the Chem. Phys. Left.19’16.40, 53. (c) Ditchfield, R.; McKinney, R. E. Chem.
external magnetic field. This idea was first suggested and used Phys. 1976,13. 187. (d) Ditchfield, R. frog. Theor. Org. Chem. 1977.2.503.
by London in the study of molecular diamagnetism over 50 years (e) Hinton, J. F.; Bennett, D. L. Chem. Phys. Lett. 1985, 116, 292. (f)
McMichael Rohlfing, C.; Allen, L. C.; Ditchfield, R. Chem. Phys. Lett. 1982,
ago.6 For magnetic shielding calculations, it was first adapted 86, 380. (9) Chesnut, D. E. Chem. Phys. 1986, 110, 415. (h) Chesnut, D.
by Ditchfield in the gauge invariant atomic orbitals (CIAO) B.; Foley, C. K. J. Chem. Phys. 1986, 84, 852.
method.’ In Ditchfield’s method each atomic orbital has its own (11) (a) Chesnut, D. B.; Foley, C. K. Chem. Phys. Lett. 1985, 118, 316.
local gauge origin placed on its center. Following Ditchfield’s (b) Chesnut, D. B.; Moore, K. D. J . Comput. Chem. 1989, 10, 648.
(12) (a) Kutzelnigg, W. Isr. J . Chem. 1980, 19, 193. (b) Schindler, M.;
work, the CIAO method has been also implemented by Giess- Kutzelnigg, W. J . Chem. Phys. 1982,76, 1919. (c) Schindler, M.;Kutzelnigg,
W. J . Am. Chem. SOC.1983, 105, 1360. (d) Schindler, M.; Kutzelnigg. W.
Mol. Phys. 1983, 48. 781.
’
Permanent address: Institute of Chemistry, Maria Curie-Sklodowska
University, 20031 Lublin, Poland.
(13) (a) Hansen, A. E.; Bouman, T. D. J. Chem. Phys. 1985,82, 5035.
(b) Hansen, A. E.; Bouman, T. D. J . Chem. Phys. 1989, 91, 3552.

0002-78631901 I5 12-825 I SOZ.SO/O 0 1990 American Chemical Society

8252 J . Am. Chem. Soc., Vol. 112, No. 23, 1990 Wolinski et al.
have been remarkably successful. However, they may encounter derivation of the CIAO CHF7 method uses the S C F orbitals
difficulties in molecules with delocalized electron structure and, explicitly, we found it convenient to formulate it in terms of the
consequently, strong ring currents. Moreover, like all methods Fock-Dirac density matrix,I5 which is invariant with respect to
based on localized MOs, the results are slightly ambiguous when unitary transformation of the orbitals. In this formalism with a
there is no unique localization. One of the purposes of the present nonorthogonal and fielddependent basis set, the CHF perturbation
paper is to compare a modern CIAO implementation with the equations can be written in matrix form as
IGLO and LORG methods. F'DOSO + F'D'SO + F'DOS' = S'DOF' + SOD'F' + SODOF'
In spite of the good performance of the methods based on (6)
localized orbitals, an efficient implementation of the CIAO method
would be, in our opinion, preferable for the following reasons:
(i) Compared with the CIAO method, the localized methods D'SoDo + D W D 0 + DoSoD1= 2Dl (7)
introduce a further approximation in the form of the closure Equation 7 arises from orthonormality constraints, while eq 6
relation, which holds only for complete basis sets. represents the Hartree-Fock equations in first-order perturbation
(ii) The wave function is more flexible in the CIAO method theory. In eq 6, and in the following, we do not explicitly denote
than in the localized ones. Both (i) and (ii) would lead to the the Cartesian components, only the order of the perturbation due
expectation that the localized methods are more sensitive to the to the external magnetic field. F, D, and S are the Fock, density,
quality of the basis set than the GIAO technique. and overlap matrices, respectively. Symbols without superscript
(iii) The localized theories are not expected to perform well denote zeroth-order (unperturbed) quantities. The first-order Fock
for systems with inherently delocalized orbitals. The electronic matrix is given by
structure of most simple organic molecules (even that of conju-
gated polyenes) is well localizable. Nevertheless, a number of F' = h' + C(DI,$) + G(Do,gl) (8)
electron-deficient compounds, inorganic compounds, and unusual where
structures exist where localizability is questionable or'poor.
(iv) Finally, generalization of the CIAO formalism to correlated G(D',$),, = EfsD'fsgOpq,rs
wave functions is straightforward. Such generalization appears
to be more difficult for the localized theories.
11. Theory and
Below we recapitulate the theory of magnetic ~ h i e l d i n g ~in? ~ ~ ' ~ h',, = (p'lhlq) + (plh'lq) + (plhlq')
order to establish the notation. The magnetic shielding constant,
B,, is the second derivative of the molecular energy with respect The first derivative of the bare nucleus Hamiltonian with respect
to the external magnetic field, H,and the nuclear magnetic to the external magnetic field is h' = (-i/2c)r X V. The two-
moment of a given nucleus n. This second-order quantity is electron integrals are denoted by gfqfs= $pqfs = (pqlrs) - 0.5-
represented by a nonsymmetric tensor with nine independent (pslrq); their first derivatives are g . The solution of the C H F
components. In the coupled Hartree-Fock theory, using atomic eq 6 provides the first-order density matrices D", a = x, y , z .
orbital throughout, the tensor component ab (a, b = x, y, z) for Equation 6 can be conveniently rewritten as
nucleus n is given by
D' = 0.5 X DOS'DO + Eia(ei - ea)-' X Cit(F' - eiS1)C, X
+
Rb= T r ( P X kbl Tr(DBoX h",} (1) (CiC,t - C,C,+) (9)
where D is the first-order reduced density matrix, and h is the where Ci and C, denote column vectors of the unperturbed oc-
matrix of the one-electron part of the Hamiltonian. The super- cupied and virtual orbitals, respectively, and the corresponding
scripts represent differentiation with respect to the external orbital energies are ei and e,. Equation 9 must be solved iteratively
magnetic field, and the nuclear magnetic moment, in this order. because F' depends on D'.
The first term in eq 1 is the diamagnetic component of the
shielding and the second term represents the paramagnetic part; 111. Efficiency of the CIAO Method
the latter depends on the first-order densit matrix Po.The
formulas for the derivative matrix elements h: Y depend on the type The amount of computer time and memory required by CIAO
of basis set employed: fixed or field dependent. In the GIAO programs was a major obstacle for the routine application of this
method with explicitly field-dependent basis functions method for NMR chemical shift calculations. In the past 20 years,
much effort has been spent on improving the efficiency of ab initio
x P ( W = exp[(-i/2c)(H X Rp)~r1xp(O) (2) analytical derivative theory. It was pointed out recentlyI6 that
these advances can be used to improve significantly the efficiency
these are the derivatives of the matrix elements of the bare nucleus of NMR chemical shift calculations. Although the physical
Hamiltonian h: perturbation caused by changing the nuclear coordinates and by
(h",), = ( p n K b l q )+ (plhib14")+ (plhiblq) (3) applying an external magnetic field is very different, mathe-
matically they have much in common. In both cases, the most
where pa denotes first derivative of the basis function, eq 2, with time-consuming part of the calculation arises from perturbations
respect to the a component of the magnetic field. The derivative in the basis sets. As pointed out above, a uniform magnetic field
one-electron operators in the above equations are given as in the z direction rotates the CIAO functions, eq 2, in the im-
aginary xy plane around the z axis. A compelling analogy with
htb = (1 /2c2)[r.(r - R) - r, (r -R)bl/lr - R13 (4) the gradient method is obtained by absorbing the gauge factor
in the center coordinates of Gaussians, which now move in the
h:b = ( - i / c ) [ ( r - R) X V],/lr - R13 (5) complex space by an amount proportional to the external magnetic
where r is the electron position vector, R is the position vector of field.I6 The possibility of representing gauge factors by single
nucleus n, c is the speed of light, and the symbols for dot and cross Gaussians with complex coordinates was first pointed out by Hall."
products have their usual significance. The presence of imaginary quantities, instead of the real per-
The first-order density matrices needed for paramagnetic
shielding are solutions of the appropriate C H F equations with the
external magnetic field as a pert~rbation."~While the original ( I S ) (a) McWeeny, R. Phys. Reu. 1962, 126, 1028. (b) D d s , J . L.;
McWeeny, R.; Sadlej, A . J . Mol. Phys. 1977, 34, 1779. (c) Wolinski, K.;
Sadlej, A . J . Mol. Phys. 1980, 41, 1419.
(14) Ditchfield, R. In Topics in Carbon-13 N M R Spectroscopy; Levy, G. (16) Pulay, P. Ado. Chem. Phys. 1987, 69, 241.
C..Ed.; Wiley: New York, 1974; Vol. I, p I . (17) Hall, G. G.In?. J . Quantum Chem. 1973, 7 , 15.
GIAO Method Jor NMR Chemical Shift Calculations J . Am. Chem. SOC.,Vol. 112. No. 23, 1990 8253
Table I. Convergence of the Calculated NMR Shielding with the Table Ill. Convergence of the Calculated Isotropic NMR Shielding
Basis Set in the CIAO Method for Methane and Methyl Fluoride with the Basis Set in the CIAO Method for Three-Membered Rings
(PPm) CvcloDroDene
Methane ~ ~~

isotropic shielding isotropic shielding

basis set energy, au '3C 1H basis set energy. au I3CH7 "C= 'Ha
6-31G -40.180513 207.36 32.93 6-3 I G -1 15.765 846 209.56 88.30 32.07
6-31G (d, p) -40.201 662 202.62 31.75 6-31G (d,p) -115.829 123 201.20 89.24 31.17
6-31 IG -40. I88 1 13 203.59 32.83 6-31 IG -I 15.788724 203.59 73.88 31.77
6-31 IG (d, p) -40.209003 197.09 31.88 6-31 IG (d,p) p) -115.850880 195.10 73.62 31.19
6-31 IG+(d,p) -1 15.852080 194.98 73.32 31.08
6-31 IG+(d,p) -40.209 080 196.97 3 1.87
6-21 1 IG+(d,p) -40.210 176 194.73 31.62 6-31 IG+(2d,2p) -1 15.857 193 193.99 73.08 30.89
6-31 IG+(2d,2p) -40.212095 196.33 31.67
Cyclopropane
Methyl Fluoride" isotropic
13C shielding 19F shielding shielding
basis set is0 par perp is0 par perp basis set energy, au "C 1H'
6-31G 142.28 198.94 113.95 489.35 439.93 514.06 6-3 1 G -1 17.007492 208.98 33.08
6-31G (d,p) 137.00 192.69 109.15 485.48 438.71 508.86 6-31G (d,p) -1 17.068 328 203.45 3 1.90
6-3 I IG 136.82 194.61 107.93 487.18 437.71 511.92 6-31 IG -1 17.029754 204.09 32.99
6-31 IG (d,p) -1 17.087969 197.97 32.07
6-31 IG (d,p) 128.36 186.95 99.08 483.06 437.59 505.79
6-31 IG+(d,p) 127.39 187.53 97.32 484.44 438.95 507.18 6-31 IG+(d,p) -1 17.088 500 197.84 32.04
6-31 I G + ( 2 d , 2 ~ ) 126.23 187.09 95.79 483.34 438.88 505.57 6-31 IG+(2d,2p) -1 17.093995 197.09 31.84
'par and perp denote the component of the shielding parallel and Ethylene Oxide
perpendicular to the C,axis; is0 is the isotropic average. In the calcu- isotropic
lations with diffuse functions, the latter were placed only on heavy at- shielding
oms. basis
set enerev. au "0 I3C 'H'
Table II. Comparison of the IGLO,' LORG,' and CIAO Results for 6-3 I G -152.783 107 375.08 172.26 30.94
I3C and I9F Absolute Shieldings and Anisotropies in CH3F and CH4 6-31G (d,p) -152.871 980 377.0 165.76 29.77
CH3F 6-31 IG -152.823643 373.75 167.15 30.86
13C shielding I9F shielding 6-31 IG (d,p) -152.905999 375.76 158.74 29.89
6-31 lG+(d,p) -152.909907 375.67 157.70 29.88
method. basis set is0 aniso is0 aniso CHI
6-31 IG+(2d,2p) -152.916 178 375.36 156.41 29.65
IGLOb
A 164 72 465 +26 219 'H shieldings for the hydrogen in the methylene groups.
B 158 81 542 +7 222
C 147 488 215 method^.'^ in particular t h e Obara-Saika recursive method a s
D 124 90 450 195 implemented by the Pople group,20have lower operation counts
11 127 87 446 -60 than the original Rys polynomial method. However, the latest
Ill 122 92 474 -66 implementations of the Rys polynomial method are much more
LORG
4-31G 165 77 403 +30 22 1 efficient than the original one and they vectorize excellently on
4-31G (d) 148 450 212 vector computers.*' For this reason, we chose the Rys method.I*
4-3 1G (d,p) 145 450 207 However, we are currently exploring t h e use of the Obara-Sai-
6-311G (d,p) 132 85 454 197 ka-type methods. No complications arise from t h e presence of
CIAO imaginary terms.
4-31G 142 85 489 -75 206 (2) Elimination of the storage for perturbed two-electron in-
4-31G (d) 135 82 484 -71 200 tegrals. This step is a bottleneck in existing implementations of
4-31G (d,p) 137 83 484 -69 20 1 the G l A O methods as well as in the first implementations of
6-311G (d,p) 128 88 483 -68 197
analytical force constant calculation. T h e bottleneck is eliminated
6-311G+(d,p) 127 90 484 -68 197 by storing only the three contributions to the derivative Fock
6-311G+(2d,tp) 126 91 483 -67 196 matrices, C(Do,g') in eq 8, instead of the integrals.22
exptlC 132 47 1 197 (3) Formulation of t h e coupled-perturbed Hartree-Fock
exptld 116.8 90 471 -90 4 * equations fully in atomic orbital ( A O ) basis, thus avoiding the
IGLO (A-D) and LORG results from ref 35, IGLO (I1 and 111) four-index transformation of the integrals that is present in other
results from ref 36. The anisotropy is u,, - ox, where z is the symme- implementations? These equations have traditionally been for-
try axis. bHuzinaga's basis sets: A, C.F(73/42), H(3/2). B, C,F- mulated in a molecular orbital ( M O ) , or mixed MO-A0 basis.23
(95/42), H(3/2). C, C,F(951/531), H(51/31). D, C,F(I 172/762), T h e full A 0 formulation was recommended by Pulay22 and Os-
H(72/62). II, C,F(951/541), H(51/31). 111, C,F(I 172/762), H- amura et al.24 For analytical force cosntant calculations, the A 0
(62/42). CReference35. dChestnut. D. B.; Phung, C.G. J . Chem. method may not be more advantageous than the MO one, par-
Phys. 1989, 9/. 6238.
ticularly on vector computers, because a large number of CHF

turbation caused by moving the nuclei in gradient theory, is only

(19) (a) McMurchie, L. E., Davidson, E. R.J. Conipur. Phys. 1978, 26,
a minor complication. 218. (b) Obara, S.;Saika, A. J . Chem. Phys. 1986. 84, 3963.
T h e following techniques have been used to implement a n (20) (a) Head-Gordon, M.;Pople, J. A. J . Chem. Phys. 1988,89, 5777.
efficient a b initio GlAO chemical shift program: (b) Gill. P. M. W.;Head-Gordon, M.; Pople, J. A. Int. J . Quanrum Chem.
( I ) Calculation of the perturbed two-electron integral by t h e Symp. 1989, 23, 269.
(21) Komornicki, A,, private communication, 1989.
quadrature technique of Dupuis, Rys, and King.I8 This is one (22) Pulay, P. J . Chem. Phys. 1983, 78, 5043.
of the modern methods for integral evaluation, although competing (23) Stevens, R. M.; Pitzer, M.; Lipscomb, W. N. J . Chem. Phys. 1963,
38, 413.
(24) Osamura. Y.;Yamaguchi, Y.; Saxe, P.; Fox, D.; Vincent, M. A.;
(18) Dupuis, M.; Rys, J.; King, H. F. J . Chem. Phys. 1976, 65. 1 1 1. Schaefer, H. F., 111 J . Mol. Srruer. 1983, /03, 183.
8254 J . Am. Chem. Soc., Vol. 112, No. 23, 1990 Wolinski et al.
Table 1V. "C NMR Shift Tensors in the Methylene Groups of Table V. Calculated NMR Shielding and Chemical Shift Tensors
Three-Membered-Ring Compounds Relative to Methane (ppm)" for Benzene
principal componentsC (a) Convergence of the Calculated N M R Shielding with the Basis
method is0 aniso AA BB CC Sets in the GIAO Method for Benzene (ppm)
Cyclopropene isotopic
IGLO (DZ) 4 89 54 13 -55 shielding
LORG (TZP) -2 107 -736 basis set enerw. au "C IH
G IAO
6-31G -2 105 43 22 -72 6-31G -230.623649 74.20 2554
6-31G (d,p) 1 103 42 30 -67 6-3 IG (d,p) -230.712713 73.21 24.49
6-31 IG 0 108 44 28 -72 DZ (9s5p/4s) -230.641136 66.17 25.14
6-31 IG (d,p) 2 109 43 34 -71 6-31 IG -230.661957 62.70 25.53
6-3 1 I G+(d,p) 2 109 43 34 -71 6-3 I IG (d,p) -230.7531 19 58.61 24.57
6-31 IG+(2d,2p) 2 108 43 34 -70 6-31 IG+(d,p) -230.755905 58.54 24.54
exptl 3 94 40 29 -59 6-31 IG+(2d,2p) -230.766561 57.73 24.20
exptl gas" 57.2
Cyclopropane
IGLO (DZ) IO 31 42 -I -10 (b) Chemical Shift Tensors in Benzene Relative to CHI (ppm)b
I l l (TZP) -3 55.5 25 6 -40 I'C
IV (QZP) -1.3 58 27 9 -40 principal
LORG (TZP) 0 53 -356 components 'H
CIAO
6-31G -2 49 26 4 -35 method is0 aniso A A BB CC is0 aniso
6-31G (d,p) -1 52 25 7 -35 IGLO (DZ) 131 211 8.69 -5.1
6-31 IG -I 53 27 8 -36 LORG (DZ) 135 197 260 141 4
6-31 IG (d,p) -I 56 26 IO -38 CHF 140 201 252 162 6
6-31 lG+(d,p) -1 57 26 IO -38 GIAO:
6-31 IG+(2d,2p) -1 56 26 IO -38 6-31G 133 193 247 147 5 7.39 -4.13
exptl -4 48 22 2 -36 6-31G (d,p) 129 184 240 141 7 7.26 -5.53
6-31 IG 141 199 263 151 8 7.30 -3.95
Ethylene Oxide 6-311G (d,p) 138 194 259 148 9 7.31 -5.33
IGLO (DZ) 31 34 74 IO 8
6-311G+(d,p) 139 196 260 147 8 7.23 -5.02
I l l (TZP) 39.3
45 48 92 6-311G+(2d,2p) 139 195 260 148 8 7.47 -5.35
V (QZW 29 13
LORG (TZP) 38 45 ab
exptl (gas) 137.2
CIAO exptl (solid) 132 181 236 148 I 1 7.24 -3.9
6-3 1G 35 43 79 20 6 Reference 38. IGLO results from ref 12c; LORG,CHF, and ex-
6-31G (d,p) 37 45 81 23 7 perimental results from ref 13, except for the gas-phase data.38 The
6-31 IG 36 45 80 22 7 A A direction is in the CCC plane parallel to the C-H bond, BB is in
6-31 IG (d,p) 38 47 84 24 7 the CCC plane perpendicular to the C-H bond, and CC is out of plane.
6-31 lG+(d,p) 39 46 85 24 9
6-31 IG+(2d,2p) 40 45 86 24 IO derivatives are estimated by the uncoupled Hartree-Fock pro-
exptl 44 37 93 19 I9 (??) cedure.28 In this approximation the terms containing the first-
bCalculated from isotropic and anisotropy shift. C A Ais the CCC or order density matrix, D', are neglected in eq 8, resulting in a
COC bisector, BB is in the CCC or COC plane, and CC is out of noniterative formula for the first-order density. It is now accepted
plane. that the uncoupled Hartree-Fock theory28 usually overestimates
second-order properties and may be in error by a factor of 2 or
3. We believe, however, that it is accurate enough to estimate
equations must be solved. Once the transformation has been made the order of magnitude of the contributions. A convenient property
the solution is significantly faster in MO basis (see, e.g., ref 16). of the present formulation is that the estimated contributions to
This argument does not apply, however to the present case where the chemical shielding can be written in a simple effectioe density
there are only three equations. matrix form. In the uncoupled approximation, the two-electron
(4)Solution of the CHF perturbation eq 9 by a conjugate part of the first-order density matrix in the external magnetic field,
gradient type method.25 This is a modification of the convergence PI,is given by
acceleration introduced by Pople et a1.26for the solution of the
C H F equations. Subsequently, it was pointed out by Wormer P' = Ci,,(ei - ea)-' X C:G(D,gl)C, x (CiC: - CaC,') (10)
et that this method is equivalent to the conjugate gradient
technique and can make use of a three-term recursion formula. Substituting this into the formula for the paramagnetic part of
The advantage of using the conjugate gradient formulation is the chemical shielding
reduced storage requirement (only three densities and the cor- on = Tr(Q'h:l]
responding residuals need to be stored) and simpler program
organization (the matrix dimensions do n ow as the iteration and extracting from the resulting formula the coefficient of the
progresses). In spite of its clear advantages, the method of derivatives of the integral with indexes p , u, p, and u yields
Wormer et aL2' has not become widely recognized. The con-
vergence rate of the modified method is identical with the original ( P * ~ I P * ~ ) ~ ( ~ ~ , , ~-Q D,,Q,,
,~ + 2D,,Q,, + DcuQyp) +
one: convergence to an accuracy of IO4 in the first-order density ( ~ * ~ I P * ~ ) ' ( ~ D-, D,,QVu
, ~ Q , ,- 2DPuQ,,+ DvaQHp)
( 1 1)
matrix is usually achieved in IO or fewer iterations. All three ( x ,
y . z ) C H F equations are solved simultaneously. where the matrix Q = Q" for the nth nucleus is defined as
(5) Screening of the perturbed two-electron integrals and the
calculation of only those that contribute significantly to the
magnetic shieldings. The contributions of the two-electron integral We precalculate the maximum absolute value of each matrix
element of Q, QF = maxQ;,, where n runs over all nuclei for
which the shielding is calculated. Using techniques similar to those
(25) Stiefel, E. L. Nail. Bur. Stand., Appl. Math. Ser. 1958, 49, I . of SchlegeP for gradient calculations, we can determine a nonstrict
(26) Pople, J. A.; Krishnan, R.;Schlegel, H. B.: Binkley, J. S. Int. J .
Quantum Chem. Symp. 1979, 13, 225. -~
(27) Wormer. P. E. S.;Visser, F.;Paldus, J. J . Compur. fhys. 1982,48, (28) Dalgarno, A. Ado. PhyJ. 1962, 1 1 , 281.
23. (29) Schlegel, H.B. J . Chem. fhys. 1982, 77, 3676
CIAO Method for N M R Chemical Shift Calculations J . Am. Chem. SOC..Vol. 112, No. 23, 1990 8255
Table VI. Convergence of the Calculated NMR Shielding with the Basis Sets in the GIAO Method for Hydrogen Sulfide and Carbon Disulfide
Carbon Disulfide
”C shielding shielding
basis set is0 paral PerP is0 paral perp
6-31G -56.46 29 1.64 -230.5 I 530.09 1062.73 263.77
6-31G (0,d) -55.32 290.58 -228.28 560.70 1062.40 309.85
6-31G (d,O) -36.89 290.14 -200.41 553.59 1062.28 299.05
6-31G (d,d) -36.94 289.67 -200.25 575.45 1062.45 331.96
6-31G+(d,d) -35.78 289.67 -198.50 578.35 1062.36 336.35
6-3 IG+(2d,d) -32.93 289.85 -194.32 582.16 1062.49 342.00
6-31 IG -53.62 291.61 -226.24 475.33 1062.56 181.71
6-31 IG (d,d) -50.72 289.90 -22 1.05 551.19 1062.35 295.61
6-31 IG (2d,d) -49.23 289.99 -218.93 560.78 1062.44 309.96
6-3 1 I G (3d,2d) -5 1.45 289.90 -222.1 2 564.99 1062.43 316.26
6-2 I I 1G (3d,2p) -53.30 289.82 -224.86 532.64 1062.41 267.76b
6-31 lG+(d,d) -54.35 289.84 -226.45 551.52 1062.34 296.1 1
6-31 IG+(2d,d) -5 I .69 289.92 -222.50 559.70 1062.39 308.36
6-31 IG+(2d,2d) -52.49 289.89 -223.67 566.94 1062.38 319.21
6-31 IG+(3d,2d) -52.44 289.86 -223.59 567.06 1062.38 3 19.40
6-21 1 IG(3d.2d) -53.30 289.82 -224.86 532.64 1062.41 267.76
6-21 1 IG+(3d,2d) -53.89 289.80 -225.74 534.20 1062.36 270.1 1
lGLO
II -45.7 5 16.0
111 -56.1 287.0 -227.7 536.1 1047.4 280.5
IV -57.1 287.1 -229.3 549.9 1068.3 290.8
Hydrogen Sulfide
‘H shielding I3S shielding
basis set is0 A B C is0 A B C
6-31 IG 32.8 26.0 29.9 42.3 767.9 648.9 738.8 916.0
6-21 I IGC 32.8 26.0 29.9 42.3 745.1 6 14.2 713.0 908.3
6-21 1 IG+ 32.8 26.0 30.0 42.3 750.5 617.4 714.4 919.7
6-21 1 l G + ( 2 ~ ) ~ 32.8 26.0 30.0 42.3 750.6 617.6 714.3 919.8
6-21 I IG+(d,p) 31.0 24.4 26.9 41.8 730.7 573.5 690.2 928.4
6-21 1 IG+(3d,p) 30.9 24.5 26.8 41.6 713.2 542.8 670.6 926.2
6-21 1 IG+(3d,2p) 30.8 24.3 26.6 41.5 715.7 546.1 672.6 928.3
6-21 1 IG+(5d,3p) 30.8 24.3 26.5 41.5 7 12.6 541.7 668.7 927.2
a (nd,md) denotes n d orbitals on sulfur and tn d orbitals on carbon IGLO results from ref 42; the basis sets 11, 111, and IV are all large; 111 is an
S(12,8,3), C(I 1,7,2) basis, comparable to our best basis set. *This is a calculation with a decontracted 6-31 IG (3d,2d) basis, made to check the
effect of contraction on the results. ‘The hydrogen s-type basis set was kept at the 6-31 IG level since preliminary calculations have shown that
further increase in the flexibility of the s basis set on hydrogen had negligible effect. dBasis set augmented by two s-type functions on the sulfur
atom.

upper limit for the contribution of a shell of integral derivatives reason, we have calculated chemical shifts for some representative
to the N M R shieldings. This estimate includes the elements of sulfur compounds.
D and Qmax, as well as estimates for the integral derivatives. (3) To explore the limits of the GlAO chemical shift calculations by
Examples for the efficiency of this scheme are given in section presenting timings for large molecules and systems with large basis sets.
Calculated shieldings are listed for the following molecules: methane,
V. The only significant disadvantage of this procedure is that it methyl fluoride. cyclopropane, cyclopropene, oxirane (ethylene oxide),
requires the storage (or the recalculation) of the matrix elements and benzene, as well as the sulfur compounds hydrogen sulfide, carbon
of h:‘. For a large molecule with all shieldings calculated, this disulfide, the sulfate and thiosulfate ions, dimethyl sulfide, dimethyl
is a substantial amount of storage. If the integral magnitude sulfoxide, and dimethyl sulfone. Representative timings are given for
estimation is not used, then an alternative program organization some of these molecules, and also for some larger systems, like styrene,
is preferable. In the latter, the first-order density in the external (E)-stilbene, the adamantyl cation, oxOcane (C70Hll). octahydrocyclo-
field is determined first, and the matrix elements of b:’ are con- paraphane (CI6HZl),and the water cluster (H20),,.
tracted with it as they are evaluated. Unless noted otherwise, the calculations have been performed using
the experimental geometries;” in most cases this means r, geometries
The above methods have been incorporated in our TEXAS ab when available. For the first-row elements Pople’s basis sets 4-31G
initio program package.30 through 6-31 IG+(2d,2~)’~ have been employed. For sulfur the McLe-
an-Chandler ( I 2s,9p) basis setJ3was used, contracted to [6s5p], and
IV. Calculations augmented by one, two, or three sets of d functions (with exponents
A large number of NMR chemical shielding calculations have been 0.515, or 0.2575 and 1.030, or 0.12875, 0.515, and 2.06) and by a set of
performed with the CIAO program; however, only a few will be dis- diffuse s and p functions (exponents 0.05 and 0.03) has been used.
cussed. The examples have been selected for the following reasons: Decontracting the 6-31 IG basis set further to 6-21 1 IG had only a very
( I ) To compare the accuracy and the efficiency of the GIAO imple- small effect, less than 2 ppm in the absolute shieldings and much less in
mentation with other methods, with special focus on the individual tensor th relative shifts. However, the M~Lean-ChandIer~~ basis set is appar-
components, since there is less experience available with them than with ently too strongly contracted in the 2p region, and decontraction to
the isotropic average. [6s,6p] (separating the exponents, which are 5.5045 and 2.2433 for
(2) To examine the convergence of the calculated shielding constants sulfur) results in significant changes in the absolute shieldings and even
with respect to the basis set quality. We are not aware of any systematic in the shifts. Essentially all of this change is due to the paramagnetic
study of this extent. Establishing the convergence behavior is particularly
important for second-row atoms, as the results for these appear to be (3 I ) Landolt-Bornstein. Zahlenwerre and Funktionen; New Series 1117;
more sensitive to basis set truncation than first-row atoms. For this Hellwege, K. H., Ed.: Springer: Berlin, 1976.
(32) Frisch, M. J.; Pople, J. A.; Binkley, J. S . J . Chem. Phys. 1984, 80,
3265.
(30) Pulay, P. Theor. Chim. Acta 1979, 50. 299. (33) McLean, A. D.; Chandler, G. S . J . Chem. Phys. 1980, 72, 5639.
8256 J . Am. Chem. Soc., Vol. 112, No. 23, 1990 Wolinski et al.
Table VII. Convergence of the Calculated NMR Shielding with the Table VIII. Isotropic N M R Chemical Shifts (ppm) in the
Basis Set in the CIAO Method for the Sulfate and Thiosulfate Ions" Thiosulfate Ion Relative to the Sulfate Ion
33s
Sulfate
3
3s 170
method thio central 170

basis set is0 is0 aniso energy, au GIAO

6-3 I G (d,d) -300.19 65.75 68.73
6-31G (d.d) 272.91 202.05 120.49 -696.779259 6-31G+(d,d) -275.06 43.94 67.14
6-31G+(d,d) 249.21 185.95 123.37 -696.823301
6-31 IG (d,d) 262.76 182.16 125.48 -696.868280 6-31 IG (d,d) -282.97 38.41 75.98
6-31 IG (2d.d) 220.57 169.38 137.28 -696.905919 6-3 I IG (2d,d) -304.29 39.49 78.82
6-3 1 1 G (3d,2d) 219.12 168.33 125.44 -696.948668 6-31 1'3 (3d,2d) -294.27 37.46 77.75
6-31 lG+(d,d) 246.72 167.65 125.38 -696.903067 6-31 lG+(d,d) -249.79 37.41 70.97
6-31 IG+(2d,d) 2 10.37 159.63 132.89 -696.9 39406 6-31 IG+(2d,d) -275.1 1 36.98 74.29
6-31 IG+(3d,2d) 213.59 162.21 135.53 -696.976374 6-31 lG+(3d,2d) -280.22 40.01 73.65
6-21 1 IG+(3d,2d) 170.88 163.76 134.46 -696.978567 6-21 I IG+(3d,2d) -295.21 47.89 74.12
Thiosulfate exDtl" 36.3 61
central 33Sthio 170 Reference 45.
basis set is0 aniso is0 aniso is0 aniso
-
6-31G (d,d) 207.16 -35.82 573.10 626.66 133.32 149.58 Table I X . I3S Shifts in the Sulfate and Thiosulfate Ions Relative to
6-31G+(d.d) 205.27 -49.34 524.27 723.94 118.81 150.22 Carbon Disulfide
6-31 1'3 (d,d) 224.35 -52.00 545.73 656.06 106.18 163.35 thiosulfate ion
6-31 IG (2d,d) 181.08 -62.40 524.86 672.93 90.56 170.06
6-31 IC (3d,2d) 181.66 -66.76 513.39 701.54 90.58 163.78 sulfate ion central thio
6-3 I IG+(d,d) 209.31 -39.66 496.51 742.23 96.68 160.00 method is0 is0 is0
6-3 I I G+(2d,d) 173.39 -48.31 485.48 755.03 85.34 164.01 CIAO
6-3 I IG+(3d,2d) 173.58 -51.65 493.81 746.98 88.56 165.10 6-31'3 (d,d) 302.54 368.29 2.35
6-21 1 1G+(3d,2d) 123.99 -53.14 466.09 782.57 89.64 164.05 6-31G+(d.d) 329.14 373.08 54.08
E (nd,md) denotes n d orbitals on sulfur and m d orbitals on oxygen. 6-311'3 (d,d) 288.43 326.84 5.46
The anisotropy is the difference of the parallel and perpendicular 6-31 IG (2d,d) 340.21 379.70 35.92
components of the shielding. 6-3 1 I G (3d,2d) 345.87 383.33 51.60
6-31IG+(d,d) 304.80 342.21 55.01
contribution. This is somewhat surprising in that the 2p function is 6-3 I 1 G+(2d,d) 349.33 386.31 74.22
expected to be largely atomic in nature and therefore should have no 6-31 IG+(3d,2d) 353.47 393.48 73.25
paramagnetic contribution. When comparing different calculations, we 6-21 1 IG+(3d,2d) 363.32 410.2 68.1
make use of the double-f or triple-f (DZor TZ) and singly polarized exptl
(+P). doubly polarized (+2P) etc. notations. For instance, the 6-31G** (a) 350.3" 386.6"
basis in this notation is DZ + P. The addition of diffuse functions is (b) 332b 368.3b
denoted by a + sign, as usual.
I n most cases, the calculated values are compared with both Har- Calculated from data in ref 45. *Reference 46.
tree-Fock limit values (obtained by large basis sets) and with experi-
mental chemical shifts. However, the Hartree-Fock limit values are with an analogous T Z + 2P basis set. However, the convergence
more relevant in comparing different techniques for the following reasons. of the GIAO results is much faster than both IGLO and LORG.
First, all methods considered here are based on the Hartree-Fock wave Even with the small 4-31G basis, the GIAO results are satisfactory,
function, and therefore, any agreement better than the Hartree-Fock and closer to the Hartree-Fock limit than the IGLO and LORG
limit must be fortuitous. Second, comparison with experimental values results with the DZ + P basis set. This faster convergence is
suffers from uncertainties in molecular geometries, from condensed phase
effects, and, in the case of the individual tensor components, from ex- particularly evident for the 19F shielding tensor. For instance,
perimental difficulties in some instances. the isotropic 19Fshielding varies for different basis sets by 92 and
Calculated N M R chemical shifts for two structures of the nor- 51 ppm in the IGLO and LORG methods, respectively; with the
bornadienyl cation and for norbornadiene, obtained with an earlier ver- GIAO method the maximum difference is only 6 ppm. The
sion of this program, have been published previously." anisotropy of the I9F shielding is even more striking. Without
polarization functions, the IGLO and LORG values are of the
V. Results and Discussion wrong sign, while all of the CIAO values are in the correct range
Table 1 shows the convergence of the calculated shielding tensors with correct signs, compared to both experiment3' and the Har-
with the basis set quality in methane and methyl fluoride. tree-Fock limit. However, the IGLO (and presumably the
Convergence is essentially achieved at the triple-{ plus polarization LORG) values improve very much when polarization functions
basis set level for 'H, I3C, and I9Fisotropic shieldings as well as are added.36
for principal components of the carbon and fluorine shielding Table I11 shows the GIAO isotropic shieldings for the three-
tensors in methyl fluoride. The influence of the diffuse orbitals membered rings cyclopropene, cyclopropane, and ethylene oxide.
is insignificant. Convergence is obtained at the T Z + P 6-31 lG(d,p) level for all
Recently, IGLO and LORG results were published for a series nuclei. The principal components of the methylene 13Cchemical
of fluor~methanes~~ with a focus on the influence of the basis set shift tensors (relative to methane; this is close to the generally
on the calculated shielding tensors. These results, and the CH3F used tetramethylsilane standard) are compared with IGLO,
data from an extensive set of calculations on fluoro compounds LORG, and experimental values in Table IV. (As the absolute
, ~ ~compared with ours in Table
by Fleischer and S ~ h i n d l e rare shielding in methane is probably larger than in tetramethylsilaneP*
I1 for CHI and CH3F. In the case of the I3C shielding the IGLO it is likely that the experimental values quoted should be increased
and CIAO results converge to essentially the same values; pre- by a few ppm. Following other workers in the field, we did not
sumably, the LORG predictions would converge to a similar value include this correction, as its magnitude is somewhat uncertain,
and it is of the same order as condensed-phase effects, which
preclude very accurate absolute comparisons of chemical shifts.)
(34) Bremer, M.; SchBtz, K.; Schleyer, P. v. R.; Fleischer, U.; Schindler,
M.;Kutzelnigg, W.;Koch, W.; Pulay, P. Angew. Chem., In?. Ed. Eng. 1989, With large TZ + P basis sets the GIAO and LORG results agree
28, 1042.
(35) Facelli, J. C.; Grant, D. M.; Bouman, T. D.; Hansen, A. E.J. Com- (37) Jameson, C . J.; Jameson, A. K.:Burrell. P. M. J. Chem. Phys. 1980,
pur. Chem. 1990. 1 1 , 32. 73, 6013.
(36) Fleischer, U.; Schindler, M . Chem. Phys. 1988, 120, 103. ( 3 8 ) Jameson, A. K.; Jameson, C. J. Chem. Phys. Leu. 1987, 134, 461.
C I A O Method for NIMR Chemical Shift Calculations J . Am. Chem. SOC.,Vol. 112, No. 23, 1990 8251

Table X. Convergence of the Calculated NMR Shielding (ppm) with Respect to the Basis Sets in the CIAO Method. 33S,"0,
and I3C
Shieldings in Dimethyl Sulfide, Dimethyl Sulfoxide, and Dimethyl Sulfone"
"S 'IC 170

basis set is0 aniso is0 aniso is0 aniso

Dimethyl Sulfide
6-3 IG (d,O,O) 673.16 389.52 191.77 23.84
6-3 1G (d,d,p) 686.08 281.58 184.99 25.40
6-3 IG+(d,d,p) 680.09 273.53 183.93 27.14
6-3 I 1G (d,d.p) 675.72 285.19 175.78 27.47
6-3 I 1G (2d,d,p) 657.22 287.53 175.41 27.52
6-31 lG+(d,d,p) 674.63 266.56 175.07 27.93
6-31 IG+(2d,d,p) 654.23 279.33 175.02 27.72
6-31 IG+(3d,2d,p) 650.52 285.62 173.77 26.95
6-21 1 IG+(3d,2d,p) 619.21 302.00 174.43 27.06
Dimethyl Sulfoxide
6-31G (d,O,O) 327.13 335.00 173.10 5 1.09 328.74 242.69
6-31G (d,d,p) 345.41 324.23 167.83 49.73 34 1.28 229.53
6-3 IG+(d,d,p) 346.84 329.50 166.23 52.28 334.37 209.79
6-3 I IG (d,d,p) 358.46 318.56 158.13 52.35 326.85 236.99
6-3 I IG (2d,d,p) 3 15.98 356.81 157.48 52.32 3 12.84 240.23
6-31 IG+(d,d,p) 355.57 314.90 157.58 52.69 327.41 226.63
6-31 IG+(Zd,d,p) 316.83 348.47 157.13 52.49 317.92 229.16
6-31 IG+(3d,2d,p) 299.48 354.92 155.61 51.60 315.60 224.51
6-21 I IG+(3d,2d,p) 250.88 376.38 156.08 51.88 317.46 226.92
IGLO I1 (2d,d,p) 273.5 147.4 318.0
IGLO 111 253.7 368.8 142.3 73.9 318.2 285.1
Dimethyl Sulfone
6-3 1G (d,O,O) 272.48 226.49 168.1 1 55.60 180.00 1 10.08
6-31G (d,d,p) 300.36 197.72 163.08 52.22 198.01 106.33
6-31G+(d,d,p) 299.26 192.14 161.68 52.94 193.28 108.93
6-31 IG (d,d,p) 3 13.40 205.15 152.96 54.99 174.89 115.71
6-31 IG (2d,d,p) 264.93 227.51 152.98 53.94 153.98 119.03
6-31 IG+(d,d,p) 307.77 203.58 152.44 55.30 177.21 1 15.29
6-31 IG+(2d,d.p) 265.87 2 19.00 152.70 54.17 161.70 122.50
6-31 IG+(3d,2d,p) 261.94 209.48 151.48 52.74 162.04 116.31
6-21 1 IG+(3d,2d,p)
.. 215.35 221.33 151.92 52.8 1 164.21 118.58
+
"The anisotropy is defined as A - (B C)/2 where A, B, and C are absolute shieldings and A > B > C. IGLO results from ref 42. Basis set I11
used in the IGLO calculations is a verv larae uncontracted basis. The basis sets designated by the plus (+) sign (diffuse functions) contain diffuse
functions only on the heavy atoms (S,'C, 6).
very well. Unfortunately, smaller (say DZ or DZ+P) LORG TZ+P basis set level. However, for relative shifts reasonable values
results are not available for these molecules; such results are a are obtained with the 6-31G basis set, as the errors cancel one
much more critical test of the method than large basis set results, another. The DZ IGLO and LORG chemical shifts are also in
and they are also more important for the projected routine ap- good agreement with experiment, presumably for the same reason.
plications for large systems. The IGLO results at the DZ level, An interesting observation is that in all our results for first-row
particularly the individual tensor components, are in much poorer elements, the CIAO isotropic shielding decreases monotonically
agreement with the large basis set values than the CIAO DZ with increasing basis set quality. No such pattern is observed for
(6-31G) data. Again, the addition of polarization functions im- IGLO or LORG (see, e&, the I9F results in Table 11). This
proves the IGLO result very m ~ c h , " ~although
,~* at considerable monotonic trend probably makes the CIAO chemical shifts more
cost. We believe that our largest basis sets (TZ+P) approach the accurate than the IGLO or LORG ones. The final absolute
Hartree-Fock limit for these molecules. Agreement between the shieldings are in excellent agreement with the gas-phase values
large basis (near-Hartree-Fock) and experimental tensor com- of Jameson and J a m e ~ o n . ~ *
p o n e n t ~is~only
~ fair. It is not clear at this time whether the Relatively few accurate N MR chemical shift calculations have
remaining discrepancy between theory and experiment is due to been published for second-row atoms beyond hydrides. Recent
the neglect of electron correlation or is mainly experimental in examples include the large-basis IGLO calculations by S ~ h i n d l e r ~ ~
its origin. For instance, the experimental assignment of the tensor and the GlAO calculations on 29Siby Van Wazer at al.43 To
components in ethylene oxide39 was based on the DZ E L 0 test the CIAO method with second-row nuclei, 33Shas been chosen
calculations,"0 which predict approximately equal BB and CC as the representative nucleus. For 33S,carbon disulfide is used
components. More accurate calculations, however, predict that as a common reference; theoretical calculations are sometimes
the out-of-plane shielding is significantly smaller than the other referred to hydrogen sulfide. Our basis set tests for H2S and CS2
components. are presented in Table VI. For H2S, decontracting the 2p region
The CIAO results for benzene are shown in Table V. As has a large effect on the 33Schemical shift. Adding diffuse
before, convergence for the absolute shielding is reached at the functions results in a change of -5 ppm, while enlarging the
McLean-Chandler basis further with s functions, or diffuse p
functions, leads only to minor changes. Polarization functions
(39) Zilm, K. W.; Conlin, R. T.; Grant, D.; Michl, J. J . Am. Chem. SOC.
1980, 102,6672. are, as expected, important, and three sets of d functions are
(40) Facelli, J. C.; Orendt, A. M.; Beeler, A. J.; Solum, M. S.;Depke, G.; needed to obtain convergence to a few ppm. Similar results are
Malsch, K. D.; Downing, J. W.; Murthy, P. S.;Grant, D. M.; Michl, J . J . Am,
Chem. Sot. 1985, 107, 6749.
(41) Orendt, A. M.;Facelli, J. C.; Grant, D. M.; Michl. J.; Walker, F. H.; (42) kichindler, M . J. Chem. Phys. 1988, 88, 7638.
Dailey, W.P.; Waddell, S.T.; Wiberg, K. B.;Schindler, M.; Kutzelnigg, W. (43) Van Wazer, J. R.; Ewig, C. S.;Ditchfield. R. J. Phys. Chem. 1989,
Theor. Chim. Acta 1985, 68, 421. 93, 2222.
8258 J. Am. Chem. Soc.. Vol. 112, No. 23, 1990 Wolinski et ai.

Table XI. 33SIsotropic Shifts (ppm) Relative to CS2 in Dimethyl Table XIV. Timings for Large Basis Set NMR Chemical Shift
Sulfide, Dimethyl Sulfoxide, and Dimethyl Sulfone Calculations on Dimethyl Sulfone, Benzene, and Phenylacetylene'
method MeSMe MeSOMe MeSOIMe energy shift ratio
basis set dimension* time time Ts/TE
GlAO
6-31G (d,O,O) - I 19.57 226.46 281.11
6-31G (d,d,p) -1 19.63 230.04 275.09 6-31 IG 28.3 2.23
6-3 lG+(d,d,p) -101.74 231.51 279.09 6-31 IG (d,d,p) 134 38.3 95.0 2.48
6-31 IG (2d,d,p) 139 44.0 108.9 2.48
6-3 I IG (d,d,p) -124.53 192.73 245.79 6-31 IG+(d,d,p) 154 60.6 153.5 2.53
6-3 I 1 G+(d,d,p) -1 23.1 1 195.95 243.75 159 69.0 183.5
6-31 IG+(Zd,d.p) 2.66
6-31 IG (2d,d,p) -96.44 244.80 295.85 6-31 IG+(3d,2d,p) 184 117.7 308.5 2.62
6-31 IG+(2d,d,p) -94.53 242.87 301.07
6-31 IG+(3d,2d,p) -83.46 267.58 305.12 6-31 IG 14.1 1.89
6-21 I IG+(3d,2d,p) -85.01 283.32 318.85 6-31 IG (d,p) 144 22.9 51.8 2.27
exptl 233 f 20' 315E 6-31 IG+(d,p) 168 36.3 93.6 2.58
312b 6-31 IG+(2d,2p) 216 89.7 225.1 2.51
GlAO (9s5p/4s2p) 72 5.5 12.4 2.26
'Reference 48. *Reference 49. CAnnunziata, R.; Barbarrella, G. IGLO (9s5p/4s2p) 72 total 20.0
Org. Magn. Reson. 1984, 22, 250.
Phenylacetylene (C,)
Table XII. I3C NMR Shift Tensor Relative to CH, (ppm) in GlAO (9s5p2d/5slp) 212 135.9 391.2 2.87
Dimethvl Sulfide and Dimethvl Sulfone' IGLO (9s5p2d/5slp) 212 231 31 1 1.34
principal components 'Timings in minutes on the Apollo DN10000 workstation, for all
method is0 aniso AA BB CC nonequivalent nuclei. The symmetry group shown may be a subgroup
of the molecular point group as TEXAS can use only idempotent sym-
MeSMe mctry operations (C2,u, i). and only those symmetry operations that
IGLO (DZ) 12 34 31 15 -1 1 transform different atoms into each other can be used efficiently.
G IAO IGLO results from refs 35 and 51, scaled by factors of 1/4035 and 451
6-31G 16.71 27.11 30.90 20.99 -1.76 to account for the speed of the computers used (see text).
6-31G (d.d,p) 17.63 25.40 30.75 21.45 0.70
6-3 I IG 20.62 28.89 35.59 24.34 1.36 not surprizing in a molecule with low-lying T orbitals.
6-3 I IG (d,d,p) 21.31 27.47 35.71 25.23 3.00 The thiosulfate ion is an interesting case. Only one 33Sreso-
6-31 lG+(d.d,p) 21.90 27.93 36.08 26.33 3.28 nance singal was observed, although the two sulfur atoms are
6-21 1 IG+(3d,2d,p) 20.30 27.06 34.44 24.20 2.26 chemically different.45 The missing signal is the result of the fact
exptl 22 28 43 20 4 that one of the sulfur nuclei has a resonance line too broad to be
observed.45 The observed 33Ssignal was first assigned to the thio
MeS0,Me sulfur atom.45 However, it was shown later that this signal belongs
IGLO (DZ) 41 85 77 66 -14 to the central sulfur atom.46 This is also supported by ab initio
GlAO electric field gradient calculation^.^^ Our calculated results for
6-31G 42.91 67.49 72.21 58.59 -2.09
6-3 I G (d,d,p) 39.54 52.22 60.96 52.94 4.73 the sulfate and thiosulfate ions are shown in Table VII. Even
with the largest basis sets employed, convergence is not yet fully
6-31 IG 46.88 65.29 75.10 62.15 3.33 reached, although the final results for both 33Sand 170seem to
6-31 IG (d,d,p) 44.13 54.99 66.54 58.39 7.47
44.53 55.30 66.83 59. IO
be quite close to the convergence limit. The relative values shown
6-3 I I G+(d,d,p) 7.67
6-21 I IG+(3d,2d,p) 42.81 52.81 63.47 57.36 7.60 in Table VI11 agree very well with the experimental values. The
signal of the central sulfur in S2032-is shifted downfield relative
exotl 44 53 62 62 9 to S042- by 40 ppm (experiment, 36 ppm). The calculated
"IGLO and experimental results from Solun, M. S.;Facelli, J.; downfield shift of 170nuclei, 74 ppm, agrees reasonably with the
Michl, J.; Grant, D. M. J . Am. Chem. SOC.1986, 108, 6464. experimental 61 ppm.& The calculated 33Schemical shifts relative
to CS2, shown in Table IX, agree well with experiment, in view
Table XIII. Effect of the Neglect of the Perturbed Two-Electron of the difficulties of the latter. In these negative ions, diffuse s
Integrals in Dimethyl Sulfone, with the 6-31 IC (d,p) Basis Set and p functions are, understandably, important. Their inclusion
threshold no. of time, "S shielding, significantly decreases the calculated shielding.
I 0- 2e int. S vpm Calculated shieldings for a different class of sulfur-containing
6 I 247 495 1598 359.49 molecules, dimethyl sulfide, dimethyl sulfoxide, and dimethyl
7 1 494 012 1722 358.42 sulfone, are presented in Table X. As one would expect, con-
8 1 712 058 1835 358.46 vergence with the basis set is fastest for 13Cand the slowest for
9 1 902 199 1936 358.46 shielding. Diffuse functions appear important only for I7O
IO 2077 310 2027 358.46 shielding, probably because of the partial negative charge on the
I5 2 750960 2349 358.46 oxygen atom. Table XI contains the calculated and experimental
33Schemical shifts in these molecules, relative to CS1. In the case
obtained for CS2. For this molecule the paramagnetic contribution of two contradicting experimental values for the s u l f o ~ i d e , 4the
~*~~
to the parallel component of the shielding vanishes by symmetry. theoretical results are intermediate between the results of Ret-
The remaining diamagnetic component, a s a simple expectation cofsky et al.48 and Belton et al.49 Table XI1 compares the cal-
value of i ' ,is insensitive to the basis set. The I3C shielding is culated I3C shifts (relative to methane) with IGLO and experi-
essentially converged with the 6-31 IG+(Zd,d) basis set to -52.4 mental results. For the sulfide and the sulfone, the agreement
ppm, close to Schindler's estimate4*of the Hartree-Fock limit. between the experimental results and the GIAO calculated values
However, the sulfur shielding needs a 6-21 11G+(3d,2d) basis set
(44) Wasylishen, R. E.; Connor, C.; Friedrich, J. 0. Con. J . Chem. 1984,
for convergence. Additional d functions, and also f functions on 62, 981.
the sulfur (not shown in Table VI), cause little further change. (45) Lutz, 0.;Nepple, W.; Nolle, A. 2.Nuturforsch. A 1976, 3 / . 978.
It is somewhat disappointing that the final value for the 33S (46) Hinton, J. F.; Buster, D. J . Mugn. Reson. 1984, 58, 324.
isotropic shielding is still 47 ppm too low compared to the ex- (47) Hinton, J. F. In Annual Reports on NMR Spectroscopy: Weeb. G.
A,, Ed.; Academic: London, 1987; Vol. 19. p 191.
perimental absolute shielding, 581 ppm.# The theoretical absolute (48) Retcofsky, H. L.; Friedel, R. A. J . Am. Chem. Soc. 1972, 94,6579.
shielding for I3Cshift is also 46 ppm lower than the experimental (49) Belton, P. S.;Cox, I . J.; Harris, R. K. J . Chem. Soc., Faraday Trans.
gas-phase value?8 -8 ppm. This appears to be a correlation effect, 2 1985, 81. 63.
CIAO Method for N M R Chemical Shijt Calculations J . Am. Chem. SOC.,Vol. 112, No. 23, 1990 8259
Table XV. Timings of the NMR Chemical Shielding Calculation for Larger Molecules"
molecule, formula, symmetry usedb geometry, total energy basis set, dimensione Tint TSCF Tshieldin8

styrene, C8H8,CI 4-21G adjusted! -307.147 34 4-31G, 88 14.6 7.5 57.0

styrene, CBH8,C, 4-21G adjusted! -307.582 12 6-31G(d), 128 48.7 19.2 204.2
oxocane, C7HI40,C, boat-chair I 4-31G, -347.563 18 4-31G, 100 13.8 8.5 50.9
oxocane, C7HI40,C, boat-chair 3 4-31G. -347.56385 4-31G, 100 26.8 13.3 97.6
adamantyl cation, CIoHl5+ (C,) 6-31Gt,' -387.178 89 6-31G(d), 170 71.1 46.1 285.1
stilbene, CI4Hl2,Cz (trans) 4-21G adjusted! -536.37603 4-31G, 150 38.1 20.3 128.4
stilbene, C,,HIz, C2 (trans) 4-21G adjusted," -537.13626 6-31G(d), 220 130.1 65.0 523.4
C16H2J CI MMX, -618.89678 4-31G, 192 233.3 174.2 1164.9
(H20)17 cs X-ray diffr, -1290.638 21 4-31G, 221 288.4 94.0 859.4
"Times in CPU minutes on the Apollo DN loo00 workstation, energies in atomic units. *The symmetry used is sometimes a subgroup of the full
molecular point group; see the footnote in Table XIV. CNumberof contracted basis functions. dHargitai, R.; Szalay, P. G.; Fogarasi, G., to be
published. For the adjustment, see: Pulay, P.; Fogarasi, G.; Pongor, G.; Boggs, J. E.; Vargha, A. J . A m . Chem. Soc. 1983, 105, 7037. 'Professor
P. v. R. Schleyer, private communication. /Ipso-octahydroparacyclophane,geometry optimized by the MMX molecular mechanics method.
dimethyl sulfone, benzene, and phenylacetylene, while Table XV
shows computer times for larger systems: the adamantyl cation,
two conformers of oxocane, styrene, stilbene, an isomer of the
strained molecule octahydroparacyclophane, and a microscopic
ice cluster, (H20)17. Table XIV also includes a comparison with
0 5 1 2 . 1 1 ' calculations were made on different
IGLO t i m i n g ~ ~ ~As, ~the
I 31 2 37
0.19 2 . 2 3 computers, the comparison is approximate; we made the (con-
servative) assumption that the VAX 11/750 computer used in
(cl
0 . 7
" 7
ref 35 is 40 times slower overall than our Apollo DN10000, and
that the CYBER 205, with two vector pipes, is 4 times faster than
the DN 10000 for vectorized code. Based on the two comparisons
(benzene with DZ basis, C, symmetry, and phenylacetylene with
a 9sSpld basis), our implementation of CIAO and the IGLO
program of Kutzelnigg and Schindler shows comparable effi-
1 . 5
I 1
-0 7
- 1 . 2
ciencies. IGLO ought to be more efficient conceptually, and we
surmise that the Gaussian lobe integral program compromises the
Figure 1. I3C NMR chemical shifts in styrene and (E)-stilbene. All efficiency of the IGLO implementation.I2
values in ppm relative to the para position; positive values correspond to The timings relative to the energy (integrals + SCF) calcula-
a downfield shift. (a) and (b): calculated 4-31G (upper number) and tions are also important, because they are less influenced by details
6-31G(d) (lower number) chemical shifts. (c) Comparison of the theo-
retical 6-3 IG(d) (upper number) and experimental (lower number) of the coding. The results in Table XIV and XV show that GIAO
chemical shifts in styrene. Note the averaging of the ortho and meta chemical shielding calculations take, on the average, 2.5 times
positions, due to internal rotation. See: Stothers, J. B. Curbon-13 N M R longer than the calculation of the energy. Therefore, if the energy
Spectroscopy; Academic: New York, 1972; pp 71 and 197. A somewhat calculation is feasible, then the NMR chemical shielding can also
newer reference (Hamer, G. K.;Peat, 1. R.; Reynolds, W. F. Can. J . be calculated. Nevertheless, work is in progress to increase the
Chrm. 1973,5/, 897). on which most recent compilations are based, gives efficiency of the chemical shielding program further. The
slightly different values for the a and @chemicalshifts (+9.5 and -14.3 shielding/energy timing ratio is close to 1 for IGL0,51suggesting
ppm, respectively). The overall agreement with the theoretical results that ultimately an IGLO calculation may be 1.8 times faster than
remains about the same. (d) Comparison of the theoretical 6-31G(d) a CIAO one. LORG, in its present implementation" includes
(upper number) and experimental (lower number) "C chemical shifts an N5 integral transformation step (in our opinion, this could be
in (E)-stilbene. Note the averaging in the orth and meta positions, due
to internal rotation. The experimental shifts were measured at -50 OC eliminated), and thus it is probably not competitive for large
in chloroform-d. See: Meic, Z.; Vikic-Topic, D.; Gusten, H. Org. Mugn. systems.
Reson. 1984, 22, 239. The chemical shift of the a-carbon relative to the Our styrene and stilbene results are summarized in Figure I .
para ring carbon changes decreases slightly with decreasing temperature, It is interesting to note that even minor differences in chemical
indicating that at least part of the observed sign discrepancy (-0.7 vs +0.6 shifts, like the 0.7 ppm difference between the para and meta
ppm) is due to the effect of internal motion. See Table XV for the carbon in styrene, are correctly reproduced by the calculations,
geometries. as shown in Figure I .
is almost perfect for both isotropic shift and anisotropy. At the VI. Conclusions
double-!: level (6-31G in our calculations), the CIAO results
generally agree better with experimental than the IGLO ones, The comparison of the CIAO and localized (IGLOIZ and
in agreement with a trend noted earlier. The IGLO data in Table LORG1j) methods shows that the latter are more sensitive to the
XI1, particularly those on dimethyl sulfone, illustrate another point. quality of the basis set employed. In other words, convergence
It is possible to obtain fortuitously good agreement for the isotropic of calculated chemical shieldings is faster with the GIAO method.
chemical shift and have significant errors in the individual tensor In particular, calculations with double-!: basis sets provide quite
components. good results for organic molecules with the CIAO method. The
We conclude this section with some timing results obtained on localized method need polarization functions to achieve comparable
an Apollo DN10000 workstation. Table XI11 illustrates our accuracy. This is what one would expect, since the CIAO ap-
procedure of estimating and selectively neglecting two-electron proach internally extends the basis set with higher angular mo-
integrals, described in section Ill. For a molecule of modest size, mentum orbitals, which are necessary for the correct description
like dimethyl sufoxide, raising the threshold in the integral de- of the perturbed system. In the localized methods, the flexibility
rivative evaluation from IO-15 to IC7produces a savings of almost of this extension is reduced, as all atomic orbitals participating
a factor of 2 with little loss of accuracy. Although this is sig- in a localized molecular orbital share the same gauge factor,
nificant, more impressive savings have been reported for gradientz9
and SCF50 calculations. The reason for this is not quite under- (50)Haser. M.; Ahlrichs, K. J. C ~ m p uCheni.
~. 1989, 10, 104.
stood, and work i s continuing on the integral neglect part of the ( 5 1 ) Kutzelnigg, W.; Fleischer, U.;Schindler, M. Ab Initio Calculation
and Interpretation of N M R Chemical Shifts and Magnetic Susceptibilities
program. by Means of the IGLO Method. In NMR-Basic Principles and Progress;
Table XIV contains timings for large basis set calculations on Springer: Heidelberg, in press.
8260 J . Am. Chem. SOC.1990, 112, 8260-8265

roughly the average of the atomic gauge factors. The greater chemical shift calculations, are available for a nominal fee from
flexibility of the CIAO wave function results in faster convergence. the University of Arkansas. The whole TEXAS system, including
As the CIAO method can thus employ smaller basis sets, it may many other features, some unique, such as a very compact integrals
be ultimately the more efficient method, in spite of the savings storage, UNO-CAS gradients, and a very efficient geometry
offered by the localized techniques. optimization using automatically generated internal coordinates,
Another conclusion emerging from our calculations is that, will also be made generally available in the fall of 1990 from the
unless large basis sets are used, the accuracy of individual tensor University of Arkansas.
components is often lower with the localized methods than with
the CIAO. This is true even in the cases where the isotropic
average is correctly predicted by IGLO and LORG. By contrast, Acknowledgment. We gratefully acknowledge partial support
the tensor components are usually fairly well predicted by the for this project by the National Science Foundation under Grant
CIAO method. CHE-8500487. We thank IBM Corp. for financial support and
Finally, we would like to emphasize that our implementation for the loan of an PC-RT workstation. The Apollo DNlOOOO used
of the GIAO method can be routinely used for N M R chemical for the large calculations was purchased with funds provided by
shift calculations on fairly large molecules, say in the Cl0-C,, N S F under Grant CHE-8814143. We thank Prof. P. R. v.
range, on departmental workstations and minicomputers. Indeed, Schleyer for the ab initio geometry of the adamantyl cation, Prof.
all results shown in this paper have been obtained on minicom- G. Fogarasi for the ab initio optimized geometry of stilbene, Profs.
puters. Our program currently runs on the IBM 4381 under CMS N . Allison and W. L. Meyer for the MMX optimized geometry
and VSI, on the IBM RT PC workstation under AIX, and on of octahydroparacyclophane, and Mr. P. Taylor for the ab initio
the Celerity 1200 and the Apollo DN10000 workstations under geometries of oxocane. We are also grateful to one of our re-
Unix. These versions, including integral, closed-shell SCF, and viewers for careful and constructive criticism.

Structure and Energetics of C2H4Br+: Ethylenebromonium Ion

vs Bromoethyl Cations
Tracy P. Hamilton and Henry F. Schaefer, III*
Contribution No. 91 from the Center for Computational Quantum Chemistry, University of
Georgia, Athens, Georgia 30602. Received March 8, 1990

Abstract: The cyclic and acyclic isomers of C2H4Br+are compared by using ab initio quantum mechanical techniques, including
the use of electron correlation and polarization functions. The cyclic bromonium ion is found to be more stable than the acyclic
1 -bromoethyl cation by 1.5 kcal/mol, in very good agreement with experiment. A transition state for the interconversion of
thcse two forms is reported, the energy barrier being -25 kcal/mol. The relative energies of the cyclic and acyclic minima
are remarkably insensitive to basis set and electron correlation effects, validating the results of previous low-level studies. The
2-bromoethyl cation does not exist as a minimum on the potential energy surface and spontaneously collapses to the bromonium
ion upon torsion.

Introduction The open 2-bromoethyl cation is expected to be unfavorable

The addition of Br2 to olefins is usually highly stereospecific, because of the positive charge residing on a primary carbon:
and the postulation in 1937 of a cyclic bromonium ion as an methyl substitution should stabilize this carbonium center relative
intermediate by Roberts and Kimball’ was a radical suggestion. to the bridged isomer. The positive charge on the carbon in the
The detection of these intermediates was to come 30 years later 1-bromoethyl cation is stabilized by Br, which is able to back-
in the pioneering NMR studies in superacid media from Olah’s donate electron density. The global energy minimum should be
group.* In the interim, the mechanism for addition of Br2 to either the open 1-halo form or the onium ion. The trend among
double bonds via a bridged intermediate was widely accepted halogens from experiment favors the open I-halo isomer for F and
because it neatly accounted for why addition was anti: Br+ is CI and bridging for Br and I.4,s There may not even be a cyclic
added to the double bond, forming a 3-membered ring, and attack ethylenefluoronium ion.
by the remaining Br- must come from the other side. The positive Ion cyclotron resonance experiments form Beauchamp’s group4
charge in the ethylenebromonium ion does not reside on the Br found two non-interconverting isomers of C2H4Br+and attributed
atom at all but primarily on the carbon atoms. Thus there is no the spectra to the 3-membered ring being 1.4 kcal/mol lower than
attraction between the Br atom in the ring and Br-, and steric the I-bromoethyl cation which was made from a different source.
factors are free to exert their full effect. A recent review by This procedure has been criticized because the two different modes
Ruasse) gives an account of the bromination of olefins, with of product production may lead to differing levels of internal
consideration of whether bromonium ions or @-bromocarbocations energy in the two isomers6 Because of its stereoselectivity, it
are intermediates for various alkenes. has recently been reported that C2H4Br+ has been used as a

( I ) Roberts, 1.; Kimball, G.E.J . Am. Chem. SOC.1937, 59, 947. (4) Berman, D. W.; Anicich, V.;Beauchamp, J. L. J . Am. Chem. SOC.
(2) (a) Olah. G.A.; Bollinger, J. M. J . Am. Chem. Soc. 1968.90.947. (b) 1979, 101, 1239.
Olah. G. A.: Bollinger, J. M.; Brinich, J. M. J. Am. Chem. SOC.1968. 90, ( 5 ) Olah, G . A. Halonium Ions; Wiley: New York, 1975.
2587. (6) For example, see: Angelini, G.; Speranza, M. J . Am. Chem. Soc. 1981,
(3) Ruasse. M. F. Acc. Chem. Res. 1990, 23. 87. 103, 3792 and ref 14 therein.

0002-7863/90/ I5 12-8260$02.50/0 0 1990 American Chemical Society