J. Phys. Chem. A 2010, 114, 7179–7186
7179
Shielding Constants and Chemical Shifts in DFT: Influence of Optimized Effective Potential and Coulomb-Attenuation Michael J. G. Peach,†,‡ John A. Kattirtzi,‡ Andrew M. Teale,§ and David J. Tozer*,‡ Institut fu¨r Theoretische Physik, Freie UniVersita¨t Berlin, Arnimallee 14, D-14195 Berlin, Germany, Department of Chemistry, UniVersity of Durham, Durham DH1 3LE, United Kingdom, and Department of Chemistry, Centre for Theoretical and Computational Chemistry, UniVersity of Oslo, P.O. Box 1033, Blindern, N-0315 Oslo, Norway ReceiVed: March 18, 2010; ReVised Manuscript ReceiVed: May 19, 2010
The influence of the optimized effective potential (OEP) and Coulomb-attenuation on shielding constants and chemical shifts is investigated for three disparate categories of molecule: main group, hydrogen bonded, and transition metal systems. Expanding the OEP in the orbital basis leads to physically sensible exchangecorrelation potentials; OEP generalized gradient approximation results provide some indication of the accuracy of the expansion. OEP uncoupled magnetic parameters from representative hybrid and Coulomb-attenuated functionals can be a dramatic improvement over conventional results; both categories yield similar accuracy. Additional flexibility is introduced by expanding the OEP in an extensive even-tempered basis set, but this leads to the well-known problem of unphysical, oscillatory potentials. Smooth potentials are recovered through the use of a smoothing norm, but deficiencies in the procedure are highlighted for transition metal complexes. The study reiterates the importance of the OEP procedure in magnetic response calculations using orbitaldependent functionals, together with the need for careful attention to ensure physically sensible potentials. It also illustrates the utility of Coulomb-attenuated functionals for computing short-range molecular properties. 1. Introduction In the absence of an external perturbation, the electronic energy in Kohn-Sham density functional theory (DFT) is written as
E ) Ts + Ene + EJ + Exc + Enn
(1)
where Ts is the kinetic energy of the Kohn-Sham noninteracting system, Ene is the nuclear-electron attraction energy, EJ is the classical Coulomb electron-electron repulsion energy, Exc is the exchange-correlation energy, and Enn is the nuclear-nuclear repulsion energy. The electronic energy is evaluated using the one-electron orbitals (or the density derived from them) that comprise the solutions to the one-electron Kohn-Sham equations
R + βerf(µr12) 1-R-βerf(µr12) 1 ≡ + r12 r12 r12
1 - ∇2 + υne(r) + υJ(r) + υxc(r) φp(r) ) εpφp(r) 2
]
[
Many approximations for Exc have been developed. The local density approximation1-3 (LDA) and generalized gradient approximations4,5 (GGAs) are explicit functionals of the densitysthe LDA depends only on the density at each point in space, whereas GGAs have an additional dependence on the gradient of the density. Meta-GGA6 and hybrid7 functionals introduce explicit occupied orbital-dependence. The former achieves this through the kinetic energy density and the latter through a fixed amount of orbital-dependent exact exchange (the orbital-dependent Hartree-Fock exchange energy expression, evaluated with the Kohn-Sham orbitals). Coulomb-attenuated (or range-separated hybrid) functionals8-10 also include orbitaldependent exchange, however the amount explicitly depends on the interelectronic distance. This dependence is typically achieved by partitioning the 1/r12 operator using the identity8-11
(4)
(2) where each of the potentials υ(r) is the functional derivative, with respect to the density F(r), of the corresponding energy in eq 1; in particular the exchange-correlation potential υxc(r) is
υxc(r) )
δExc δF(r)
(3)
* To whom correspondence should be addressed. E-mail: d.j.tozer@ durham.ac.uk. † Freie Universita¨t Berlin. ‡ University of Durham. § University of Oslo.
and evaluating the exchange energy associated with the first (long-range) and second (short-range) components using an appropriately modified orbital expression and a GGA expression, respectively. The dimensionless parameters R and β control the amount of exact exchange at short- and long-range; R is the amount at short-range and R + β is the amount at long-range. Coulomb-attenuated functionals are attracting significant interest in the literature,12-25 due to their improved description of longrange properties, particularly charge-transfer excited states. To date, there have been relatively few applications16,25 of Coulombattenuated functionals to the calculation of magnetic response parameters. The present study considers their application to shielding constants and chemical shifts.
10.1021/jp102465x 2010 American Chemical Society Published on Web 06/08/2010
7180
J. Phys. Chem. A, Vol. 114, No. 26, 2010
Peach et al.
The inclusion of explicit orbital dependence in Exc has significant implications for the calculation of second-order magnetic response parameters. Since LDA and GGA functionals are explicit functionals of the density (i.e., they have no explicit orbital dependence), the multiplicative exchange-correlation potential in eq 3 is trivial to evaluate. The magnetic Hessian matrix (which determines the imaginary component of the orbital response to a static magnetic field) is diagonal and so shielding constants (and other second-order magnetic properties) are determined in an uncoupled manner, using the Kohn-Sham orbitals and orbital energies. By contrast, meta-GGA, hybrid, and Coulomb-attenuated functionals have an explicit dependence on the orbitals, making the evaluation of the potential in eq 3 nontrivial. In conventional implementations, an alternative set of orbital equations are therefore used, obtained by minimizing the energy with respect to the orbitals, as opposed to the density. The exchange-correlation term in the equations is then nonmultiplicative, analogous to the Hartree-Fock exchange operator. In this case the magnetic Hessian matrix becomes nondiagonal and so shielding constants are evaluated in a coupled manner. In 2001, Wilson and Tozer26 used the Zhao-Morrison-Parr27 procedure to determine multiplicatiVe exchange-correlation potentials associated with a hybrid functional, for a series of main group molecules. They then used the resulting Kohn-Sham orbitals and orbital energies to determine shielding constants in an uncoupled manner. The results were on average two to three times more accurate than those determined using the conventional approach with a nonmultiplicative exchangecorrelation term and coupled equations. Subsequent studies made similar observations25,28-39 for shielding constants, chemical shifts, magnetizabilities, and rotational g tensors. To date, such studies have only considered standard hybrid functionals, with the exception of ref 25, which presented rotational g tensors and magnetizabilities for a Coulomb-attenuated functional. The rigorous methodology for determining the multiplicative exchange-correlation potential associated with an orbital-dependent functional is the optimized effective potential (OEP) approach; this procedure was used in refs 25, 35, and 37-39. This involves the determination of the multiplicative potential in eq 3, with respect to which the energy is stationary. In the present study, we use the procedure of Yang and Wu,40,41 (YW) in which the Kohn-Sham effective potential is written as
υs(r) ) υne(r) + υref(r) +
∑ btgt(r)
(5)
t
where υref(r) is a fixed reference potential, and the final term is an expansion in a Gaussian basis set {gt(r)}. The unknown parameters b are determined by directly minimizing the total electronic energy with respect to these parameters. Comparing this potential with the standard υs(r) [see eq 2] allows the exchange-correlation potential to be identified as
υxc(r) ) υref(r) - υJ(r) +
∑ btgt(r)
(6)
t
It is well-established42-47 that this minimization can lead to unphysical, oscillatory exchange-correlation potentials when the basis set used in the expansion of the potential [{gt(r)} in eq 5] is particularly different to the basis set used to represent the one-electron orbitals in eq 2. In these circumstances the basis sets are termed “unbalanced”. It is even possible for the energy
to collapse to that obtained using a nonmultiplicative operator.48,49 Although such behavior is natural in small basis sets, where the matrix representation of a nonlocal operator can be replaced to high precision by the matrix representation of a local one,50,51 in larger basis sets this requires the occupation of orbitals differing from the Aufbau choice,48,49 an option not explored in the present work. Choosing the potential basis to be the same as the orbital basis reduces the problem of oscillatory potentials, although this necessarily limits the flexibility of the expansion. Several procedures have been proposed that aim to reduce the occurrence of unphysical, oscillatory potentials. In ref 41, both truncated singular-value decomposition (TSVD) and Tikhonov regularisations were used; filter values in the range of 10-4 to 10-6 were found to yield smooth exchange-correlation potentials. A TSVD was also applied in ref 39 to the transition metal complexes of the present work. Alternatively, a basis set balancing scheme could be employed, such as that in ref 44 that requires the modification of both orbital and potential basis sets. A fourth alternative, which allows the use of arbitrary orbital and potential basis set combinations is the L-curve method of refs 45-47, which introduces a smoothing norm ||∇υb(r)|| into the energy minimization. The norm provides a measure of how oscillatory a potential is, forcing the minimization to take into account the smoothness of the potential. The functional to be minimized becomes
Ωλ(b) ) EYW(b) + λ||∇υb(r)||2
(7)
where EYW(b) is the regular OEP electronic energy. A series of calculations are performed using various λ values, and an optimal value (denoted λ*) is defined as that corresponding to the minimum slope of an “L-curve”, which is a plot of the logarithm of the smoothing norm as a function of the logarithm of the energy difference between the OEP and conventional energy evaluation. In practice, λ* is identified by locating the maximum in a plot of the reciprocal slope of the L-curve, as a function of log(λ). See refs 45-47 for full details. An alternative criterion to define λ* was explored in ref 47 and is suitable for use in exchange-only OEP calculations. The value of λ* is both system- and basis-set-dependent. It should be emphasized that each of the above procedures involves a subjective definition of “smoothness”; although the latter two are more systematic in their approach toward ensuring smooth potentials, those potentials obtained are not unique and depend on the criteria chosen to define “balanced basis sets” and “smoothness of the potential”. In one way or another all of these procedures limit the variational freedom in determining the OEP energy, relative to the use of large basis sets with minimal regularisation. In the present work we employ the L-curve analysis with the criteria of ref 46 since this allows us to maintain the same orbital basis sets as used for conventional calculations and the reciprocal slope criterion can be readily applied to different exchange-correlation functionals. A “blackbox” finite-basis OEP method capable of obtaining an appropriate smooth potential from a single calculation remains an important area for future development if OEPs are to take a place in mainstream quantum chemical applications. The aim of the present study is to quantify the influence of both the OEP approach and Coulomb-attenuation on shielding constants and chemical shifts for three categories of systems: main-group molecules, hydrogen bonded systems, and transition metal complexes. For the OEP calculations, we expand the potential in the orbital basis set and use GGA test calculations
Shielding Constants and Chemical Shifts in DFT
J. Phys. Chem. A, Vol. 114, No. 26, 2010 7181
TABLE 1: Isotropic Shielding Constants (in ppm) for Main-Group Molecules, Compared with Reference Values HF H 2O CH4 CO N2 F2 O3 PN H 2S NH3 HCN C 2H 2 C2H4 H2CO N 2O CO2 OF2 H2CN2 HCl SO2 PH3 d |d| a
F O C C O N F O O′ P N S N C N C C C O N N′ O C O O C N N′ Cl S O P
KT2
O-KT2
B3LYP
O-B3LYP
CAM-B3LYP
O-CAM-B3LYP
refa
411.5 329.5 195.1 6.8 -57.1 -60.4 -209.1 -1279.6 -810.6 52.9 -365.9 735.0 264.7 85.7 -19.7 120.1 63.0 -5.0 -379.6 102.0 11.5 177.2 62.7 221.3 -533.7 167.1 -42.1 -138.9 750.7 -148.5 -251.5 598.5 -17.1 19.8
411.4 329.5 194.9 5.7 -57.9 -61.5 -218.2 -1294.1 -817.1 48.9 -367.4 734.1 264.7 85.1 -20.8 119.6 62.4 -5.9 -382.9 101.4 10.9 176.8 62.3 220.9 -542.2 167.0 -42.9 -140.7 751.8 -154.2 -253.9 597.3 -19.3 20.8
411.1 327.7 188.7 -19.0 -81.1 -92.3 -250.4 -1693.3 -1129.1 -50.1 -431.7 700.5 259.9 69.1 -49.5 106.9 47.2 -24.4 -452.4 81.9 -11.4 173.1 48.9 213.5 -583.1 160.0 -60.1 -192.7 711.7 -262.5 -283.5 561.8 -69.0 69.0
413.7 329.9 188.3 -8.5 -53.0 -71.1 -222.2 -1245.5 -828.1 18.3 -365.7 703.8 260.7 74.5 -31.0 109.9 51.4 -20.3 -378.7 96.0 2.2 189.3 53.0 224.0 -552.1 161.5 -52.8 -148.4 719.1 -203.0 -225.4 562.6 -26.6 29.4
414.0 331.3 191.8 -20.2 -79.8 -92.8 -238.4 -1761.6 -1231.5 -66.0 -433.1 706.8 263.0 68.0 -49.9 106.4 47.4 -25.6 -461.0 80.9 -14.7 179.3 48.2 217.6 -561.0 162.2 -54.3 -206.9 709.7 -279.7 -276.6 567.6 -73.8 73.8
416.2 332.9 190.0 -7.7 -48.4 -69.1 -216.4 -1188.1 -822.1 32.6 -347.6 710.4 263.0 74.6 -28.0 109.7 51.7 -22.0 -367.8 97.7 2.2 196.1 53.0 227.8 -541.7 162.7 -51.1 -139.1 712.9 -199.7 -208.3 568.2 -20.5 28.1
419.7 357.6 198.4 2.8 -36.7 -59.6 -192.8 -1290.0 -724.0 53.0 -349.0 752.0 273.3 82.1 -20.4 117.2 64.5 -4.4 -375.0 99.5 11.3 200.5 58.8 243.4 -473.1 164.5 -43.4 -149.0 952.0 -126.0 -205.0 599.9
Reference 26 and references therein.
to provide some indication of the accuracy of the expansion. We then perform additional OEP calculations, expanding the potential in a much larger, even-tempered basis set, and highlight the unphysical potentials that result. We use the L-curve approach of Yang and co-workers45-47 to recover physically sensible potentials for representative main group and transition metal systems; to the best of our knowledge, this is the first time that this procedure has been applied to such complex systems. The molecules we consider are presented in Section 2, together with pertinent computational details. Results and conclusions are presented in Sections 3 and 4, respectively. 2. Systems Investigated and Computational Details First, we consider the 21 molecules of ref 26 in order to assess the performance for main-group systems. Following ref 26, we use the Huzinaga IGLO-IV basis set52,53 at near-experimental geometries (see refs 26, 54, 55, and references therein). Second, from the work of Kongsted et al.,56 we consider the influence of hydrogen bond complexation on formaldehyde, which is a model system providing insight into the effects of solvation on chemical shifts in biological molecules. We use the aug-ccpVTZ basis set for all atoms, and have confirmed that the results agree to within 1 ppm with those obtained using the mixed basis set of ref 56. To allow comparison with the ab initio reference values of ref 56, the CCSD geometries of that reference are used throughout. Third, we consider the first-row transition metal complexes of ref 39; the cc-pVTZ basis set is used on both the metal and the ligands and calculations are performed at the geometries of that reference.
As our representative Coulomb-attenuated functional we use the increasingly popular CAM-B3LYP functional of Yanai et al.11 For the hybrid functional we use the ubiquitous Becke-3Lee-Yang-Parr3,57-60 (B3LYP) functional. We also perform calculations using the Keal-Tozer-261 (KT2) GGA functional, which was specifically designed for the evaluation of magnetic response parameters. Results for each functional are evaluated conventionally, as well as through the OEP procedure; the latter are denoted using a prefix O-, for example, OEP B3LYP calculations are denoted O-B3LYP. We reiterate that for B3LYP and CAM-B3LYP, conventional calculations use a nonmultiplicative Kohn-Sham operator and coupled shielding constant evaluation; the OEP calculations use a multiplicative potential and uncoupled shielding constant evaluation. For KT2, both the conventional and OEP exchange-correlation potentials are multiplicative, so the shielding constants associated with each potential are evaluated in an uncoupled manner. Results evaluated using both approaches should therefore be identical. In practice they will differ due to the expansion of the potential in a finite basis set within the OEP procedure, so a comparison of the conventional and OEP results provides a measure of the accuracy of the OEP expansion. For the primary assessment, we use the same basis set for the expansion of the potential as is used for the orbitals. We subsequently consider the use of a much larger basis set for the expansion of the potential, specifically a Z-independent uncontracted even-tempered set based around 2n, consisting of s exponents with -3 E n E 9, p exponents -3 E n E 6, d exponents -3 E n E 3 and f exponents -3 E n E -1. All
7182
J. Phys. Chem. A, Vol. 114, No. 26, 2010
Peach et al.
basis sets use Cartesian basis functions, thereby facilitating comparison of the results presented here with those of previous works.26,31,37 For the spatially diffuse transition metal anions, we observed a notable dependency on the size and spatial extent of the numerical integration grid. We therefore use an extensive numerical integration grid for all systems, specifically the Lindh-Malmqvist-Gagliardi62 radial and Lebedev angular integration schemes, with a radial integration parameter of 10-16 and an angular grid capable of integrating spherical harmonics up to order 50 exactly. In addition, the thresholds for evaluating the density and atomic orbitals at the grid points were lowered to 10-12 and 10-15, respectively. For the calculation of the magnetic response we use London atomic orbitals.63-65 All calculations are performed with a development version of the Dalton program.66 For the reference potential in eq 5, we use39
(
υref(r1) ) 1-
R+β N
F(r2)
) ∫ dr |r -r | 2
1
2
(1-R - β)(3/π)1/3F1/3(r1) (8) where N is the number of electrons. This combines scaled FermiAmaldi67 and scaled LDA exchange potentials and is designed to ensure the appropriate long-range behavior of the potential for each functional (for B3LYP, R ) 0.2 and β ) 0; for CAMB3LYP, R ) 0.19 and β ) 0.46; for KT2 both R and β are 0). The reference potential is evaluated using the density obtained from a conventional DFT calculation using the respective functional. The OEP calculations use an approximate Newton scheme.41 For calculations that do not use an L-curve, we use TSVD regularisation to remove contributions to the Hessian that lie below a given threshold, corresponding to small singular values that can lead to numerical instabilities. Following ref 39, we use a TSVD filter of 10-8. For calculations using the L-curve, a TSVD filter is only used to prevent convergence problems with particularly small λ, in which case it is set to 10-10. 3. Results and Discussion We commence by comparing shieldings and chemical shifts calculated conventionally with those obtained from OEP calculations where the potential is expanded in the orbital basis set. Table 1 presents reference (experimental or rovibrationally corrected experimental values) and DFT isotropic shielding constants for the 21 main-group molecules, together with the associated mean and mean absolute errors, d and |d|. The KT2 and O-KT2 results are in good agreement, indicating that the OEP can be well-represented using the orbital basis set. The largest absolute discrepancy occurs for the ozone molecule, although the associated percentage discrepancy is small. The KT2 shieldings are significantly more accurate than the conventional B3LYP and CAM-B3LYP values. Applying the OEP approach to B3LYP and CAM-B3LYP leads to a significant improvement; errors are reduced by a factor of 2 to three, but remain slightly inferior to KT2. This is unsurprising given that KT2 was specifically designed to provide accurate magnetic parameters. The effect of using OEPs is large for challenging molecules such as ozone, but is much smaller for molecules such as water, where the conventional approach works well. Figure 1a presents the exchange-correlation potential for the representative molecule N2, determined using O-CAM-B3LYP. There is no unphysical structure, demonstrating that the smooth-
Figure 1. N2 exchange-correlation potentials υxc(r) obtained from O-CAM-B3LYP, with the nuclei at z ) (1.038 au, from TSVD and L-curve calculations. See text for details of basis sets.
ing norm is not required when the basis sets are balanced, providing there is some numerical regularization. The improved (although not ideal) -0.65/r asymptotic behavior of the potential, owing to the long-range orbital-dependent exchange contribution in the O-CAM-B3LYP functional, is also visible. Table 2 presents reference CCSD(T) and DFT shielding constants for a formaldehyde molecule complexed with two water molecules and for an isolated formaldehyde molecule. The shieldings for the isolated molecule are determined at the relaxed formaldehyde structure; we have confirmed that the effect of structural relaxation is small compared to the effect of the hydrogen bonding. The chemical shifts in the last two rows are the differences between the isolated and complexed shieldings and so quantify the effect of hydrogen bonding. There is again good agreement between KT2 and O-KT2, indicating that the potential can be well-represented using the orbital basis set. The KT2 results for the oxygen shielding constants (and hence chemical shifts) are in excellent agreement with the CCSD(T) values. Although the absolute values of the carbon shieldings are less accurately reproduced, the behavior
Shielding Constants and Chemical Shifts in DFT
J. Phys. Chem. A, Vol. 114, No. 26, 2010 7183
TABLE 2: Isotropic Shielding Constants and Chemical Shifts (in ppm) of Complexed Formaldehyde and Formaldehyde, Compared with Reference CCSD(T) Values
H2CO +2H2O H2CO
a
KT2
O-KT2
B3LYP
O C O C
-300.8 -6.0 -353.1 6.3
-309.5 -9.8 -362.7 2.8
Shielding Constants -362.1 -299.5 -28.6 -23.3 -426.2 -352.4 -15.0 -10.4
O C
52.3 -12.3
53.2 -12.5
64.0 -13.7
O-B3LYP
Chemical Shifts 52.9 -12.9
CAM-B3LYP
O-CAM-B3LYP
CCSD(T)a
-365.2 -30.0 -434.3 -16.0
-288.9 -24.9 -342.1 -12.2
-298.7 0.8 -350.6 13.0
69.2 -14.0
53.2 -12.6
52.0 -12.2
Reference 56.
TABLE 3: Isotropic Shielding Constants and Chemical Shifts (δ, in ppm) of Transition Metal Nuclei, Compared with Experimental Chemical Shifts (δexpt) TiCl4 TiCl3Me δ Cr(CO)6 CrO42δ MnO4Mn(CO)6+ δ VOCl3 VF5 δ VOF3 δ |d| a
KT2
O-KT2
B3LYP
O-B3LYP
CAM-B3LYP
O-CAM-B3LYP
-947.8 -1397.1 449.3 -644.4 -2587.1 1942.8 -3817.5 -1808.0 -2009.5 -1901.8 -1174.7 -727.1 -1205.4 -696.5 221
-955.7 -1422.6 466.8 -670.0 -2596.5 1926.5 -3822.1 -1841.9 -1980.2 -1905.6 -1174.0 -731.5 -1206.4 -699.2 207
-1000.0 -1479.1 479.0 -1058.5 -3141.1 2082.6 -4995.8 -2713.2 -2282.5 -2257.2 -1245.2 -1012.0 -1375.6 -881.6 300
-833.6 -1237.8 404.1 -571.0 -2479.2 1908.3 -3724.9 -1639.9 -2084.9 -1784.7 -1032.0 -752.7 -1088.1 -696.6 233
-975.8 -1478.3 502.5 -1160.1 -3203.5 2043.4 -5191.6 -2946.0 -2245.6 -2311.3 -1197.6 -1113.8 -1383.6 -927.8 310
-731.6 -1123.1 391.4 -506.8 -2424.7 1918.0 -3696.3 -1531.3 -2165.0 -1695.7 -963.3 -732.5 -1032.5 -663.3 264
δexpta
613 1795 -1445 -895 -757
Reference 39 and references therein.
of the shifts is similar to that observed with CCSD(T). KT2 is again more accurate than both conventional B3LYP and CAMB3LYP, for which the shielding constants are in poor agreement with CCSD(T). The B3LYP and CAM-B3LYP chemical shifts do benefit from some error cancellation, but remain inferior to KT2. Application of the OEP approach to B3LYP and CAMB3LYP improves the shieldings (particularly for oxygen) leading to chemical shifts that are in excellent agreement with CCSD(T) and KT2. Table 3 presents shielding constants and chemical shifts for the metal nuclei in a series of transition metal complexes. For the reference we use experimental values, but the absence of solvent and vibrational effects in our calculations makes it difficult to judge the DFT accuracy. The absolute agreement between KT2 and O-KT2 is inferior to that observed with the main group molecules, although the discrepancy in all cases is less than 34 ppm, which is small in percentage terms. The results suggest thatsas in Tables 1 and 2sKT2 is on average an improvement over conventional B3LYP and CAM-B3LYP and that application of the OEP procedure improves the results of the orbital-dependent functionals. Discrepancies compared to experiment are largest for the manganese-containing systems, however there is some question over the accuracy of the experimental chemical shift in this case.68 Figure 2a presents the exchange-correlation potential for the representative molecule, TiCl4 determined using O-CAM-B3LYP; again, there is no unphysical structure. In summary, Tables 1-3 highlight the significant influence of the OEP approach on B3LYP and CAM-B3LYP shielding constants and chemical shifts; improvements can be dramatic. By expanding the potential in the orbital basis set, smooth potentials are obtained and the comparison of KT2 and O-KT2
quantifies the accuracy of the potential expansion. Despite the different long-range properties of the two orbital-dependent functionals, they yield similar accuracy for shieldings and chemical shifts, reflecting the fact that the short-range behavior is most relevant for these properties. (This is also evident in the observation that a long-range asymptotic correction has minimal effect on shieldings.69) 3.1. Use of a Larger Potential Expansion Basis. Expanding the OEP in the orbital basis necessarily limits the flexibility. To investigate this, additional flexibility was introduced through the use of a much larger basis set expansion. For example, O-KT2 calculations on N2 and TiCl4 using the even-tempered basis set of Section 2 and a TSVD filter of 10-8 give respective shielding constants of -60.5 and -951.0 ppm, compared to values of -61.5 and -955.7 ppm obtained using the expansion in the orbital basis. The new values are closer to the conventional KT2 values of -60.4 and -947.8 ppm, indicating that the potential is better reproduced. However, applying this procedure to B3LYP and CAM-B3LYP leads to a fundamental problems the potentials become unphysical, as illustrated in Figures 1b and 2b, where the O-CAM-B3LYP potentials are presented. We now illustrate the use of the L-curve approach for alleviating the problem of unphysical potentials by applying it to O-CAMB3LYP calculations on the representative systems N2 and TiCl4/ TiCl3Me. Figure 3a presents the L-curve for N2, and Figure 3b presents the associated reciprocal slope plot used to identify λ*. There is a well-defined maximum in the latter, with an associated λ* value of 2.5 × 10-4, which is indicated in Figure 3a. Figure 3c plots the variation in the shielding constant with λ. There is a strong variation, but at the value of λ* the shielding constant is -71.0 ppm, which is in good agreement with the value in Table
7184
J. Phys. Chem. A, Vol. 114, No. 26, 2010
Peach et al.
Figure 3. L-curve and reciprocal slope plots for N2, as illustrations of the procedure used to determine the optimal value of λ, together with the variation of shielding constant with λ. Note the well-defined maximum in panel b and the strong variation in panel c.
Figure 2. TiCl4 exchange-correlation potentials υxc(r) obtained from O-CAM-B3LYP, with the Ti nucleus at z ) 0 and the Cl nucleus at z ) 4.101 au, from TSVD and L-curve calculations. See text for details of basis sets.
1. The corresponding potential is plotted in Figure 1c. It is physically sensible and closely resembles the potential in Figure 1a that was obtained using the orbital basis set. The variation of the shielding constant with λ in Figure 3c is dominated by the paramagnetic component of the shielding, the
changes in which are large compared to those observed in the diamagnetic contribution. As λ increases the shielding initially becomes more negative, reflecting the closing of the HOMOLUMO gap, on which the negative paramagnetic contribution has an inverse dependence. As λ becomes larger the potential is oversmoothed, strongly affecting the orbitals obtained, leading to both a more positive paramagnetic contribution and overall shielding constant. Figures 4a-c repeat the analysis for TiCl4. For this system, the reciprocal slope plot is much broader and actually exhibits two maxima. The absolute maximum value is associated with λ* ) 2.5 × 10-4, for which the shielding constant is -761.8 ppm. This differs by 30 ppm from the value in Table 3. The O-CAM-B3LYP potential associated with this value of λ* is plotted in Figure 2c. It is smooth, as required, and closely resembles the potential in Figure 2a. Calculations on TiCl3Me (not shown) exhibit the same featuressa broad reciprocal slope plot with two maxima. The absolute maximum value gives λ* ) 5 × 10-4, for which the shielding is -1244.0 ppm, a change
Shielding Constants and Chemical Shifts in DFT
J. Phys. Chem. A, Vol. 114, No. 26, 2010 7185 4. Conclusions
Figure 4. L-curve and reciprocal slope plots for TiCl4, as illustrations of the procedure used to determine the optimal value(s) of λ, together with the variation of shielding constant with λ. Note the broad feature in panel b and the strong variation in panel c.
of over 100 ppm from the value in Table 3. The chemical shift associated with these shieldings is 482.2 ppm, at variance with the value in Table 3 of 391.4 ppm. The broad nature of the plot in Figure 4b indicates that there is a range of λ values that yields essentially the same reciprocal slope and so the precise choice of λ* is somewhat ambiguous. This is undesirable because, as shown in Figure 4c, the shielding constant is very sensitive to this choice. For example, the second maximum in Figure 4b is associated with λ* ) 10-3, for which the shielding constant is -816.2 ppm. For TiCl3Me, the corresponding values are λ* ) 10-3 and -1302.1 ppm. Despite the changes in the shielding constants from the values determined using the absolute maxima, the chemical shift is relatively unchanged, at 485.9 ppm. The O-CAM-B3LYP potential associated with this λ* for TiCl4 is plotted in Figure 2d; the shell structure observed in Figure 2c is marginally less evident, due to the increased influence of the smoothing norm.
Both the OEP approach and Coulomb-attenuated functionals have attracted significant recent interest in the DFT literature. In the present study, we considered their influence on the calculation of shielding constants and chemical shifts in three disparate categories of moleculesmain group, hydrogen bonded, and transition metal systems. OEP calculations expanding the potential in the orbital basis set yield physically sensible potentials, when using minimal TSVD regularisation to avoid numerical instabilities. OEP GGA results provide some indication of the accuracy of the expansion. For the main-group and hydrogen-bonded systems, the use of OEP (with associated uncoupled shielding evaluation) yields dramatic improvements for the representative hybrid (B3LYP) and Coulomb-attenuated (CAM-B3LYP) approximations. The two functionals yield results of similar quality. Improvements are also observed for the transition metal systems. For representative systems, additional flexibility was introduced by expanding the OEP in a much larger basis set. Simply using an unaltered TSVD filter value results in the wellestablished problem of unphysical, oscillatory potentials, as illustrated for N2 and TiCl4. Application of the L-curve approach successfully recovers smooth potentials in both cases. For N2, a well-defined maximum is observed in the reciprocal slope plot and the associated shielding is in good agreement with that obtained using the orbital basis. For TiCl4, however, the reciprocal slope plot is broad, highlighting an ambiguity in the choice of λ*. Despite large changes in the shielding constants, the chemical shifts based on the two maxima are similar but are in rather poor agreement with that obtained using the orbital basis set. This study reiterates the importance of the OEP procedure in magnetic response calculations using orbital dependent functionals, together with the need for careful attention to ensure physically sensible potentials are obtained. It also illustrates the utility of Coulomb-attenuated functionals for computing shortrange molecular properties. For the molecules considered, O-B3LYP and O-CAM-B3LYP results remain inferior to those of KT2, which is unsurprising given that the latter was specifically designed to provide accurate magnetic response. Given the nontrivial nature of the OEP calculations, the use of functionals with well-defined potentials such as KT2 is therefore to be preferred at present. As orbital-dependent functionals and OEP methodologies improve, however, these calculations will become increasingly important. Acknowledgment. This work was supported through an EPSRC studentship (MJGP) and by the Norwegian Research Council through Grant No. 171185/V30 and 197446/V30 (AMT) and the CoE Centre for Theoretical and Computational Chemistry through Grant No. 179568/V30. References and Notes (1) Dirac, P. A. M. Proc. Cam. Phil. Soc. 1930, 26, 376–385. (2) Slater, J. C. Phys. ReV. 1951, 81, 385–390. (3) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200– 1211. (4) Perdew, J. P.; Yue, W. Phys. ReV. B 1986, 33, 8800–8802. (5) Becke, A. D. J. Chem. Phys. 1986, 85, 7184–7187. (6) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. ReV. Lett. 2003, 91, 146401. (7) Becke, A. D. J. Chem. Phys. 1993, 98, 1372–1377. (8) Gill, P. M. W.; Adamson, R. D.; Pople, J. A. Mol. Phys. 1996, 88, 1005–1010. (9) Leininger, T.; Stoll, H.; Werner, H.-J.; Savin, A. Chem. Phys. Lett. 1997, 275, 151–160.
7186
J. Phys. Chem. A, Vol. 114, No. 26, 2010
(10) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. J. Chem. Phys. 2001, 115, 3540–3544. (11) Yanai, T.; Tew, D. P.; Handy, N. C. Chem. Phys. Lett. 2004, 393, 51–57. (12) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. J. Chem. Phys. 2004, 120, 8425–8433. (13) Baer, R.; Neuhauser, D. Phys. ReV. Lett. 2005, 94, 043002. ´ ngya´n, J. G. Chem. Phys. Lett. 2005, 415, 100– (14) Gerber, I. C.; A 105. (15) Vydrov, O. A.; Scuseria, G. E. J. Chem. Phys. 2006, 125, 234109. (16) Peach, M. J. G.; Helgaker, T.; Sałek, P.; Keal, T. W.; Lutnæs, O. B.; Tozer, D. J.; Handy, N. C. Phys. Chem. Chem. Phys. 2006, 8, 558–562. (17) Fromager, E.; Toulouse, J.; Jensen, H. J. A. J. Chem. Phys. 2007, 126, 074111. (18) Chai, J.-D.; Head-Gordon, M. J. Chem. Phys. 2008, 128, 084106. (19) Peach, M. J. G.; Benfield, P.; Helgaker, T.; Tozer, D. J. J. Chem. Phys. 2008, 128, 044118. (20) Peach, M. J. G.; Le Sueur, C. R.; Ruud, K.; Guillaume, M.; Tozer, D. J. Phys. Chem. Chem. Phys. 2009, 11, 4465–4470. (21) Borini, S.; Limacher, P. A.; Lu¨thi, H. P. J. Chem. Phys. 2009, 131, 124105. (22) Jacquemin, D.; Wathelet, V.; Perpete, E. A.; Adamo, C. J. Chem. Theory Comput. 2009, 5, 2420–2435. (23) Rohrdanz, M. A.; Martins, K. M.; Herbert, J. M. J. Chem. Phys. 2009, 130, 054112. (24) Wiggins, P.; Williams, J. A. G.; Tozer, D. J. J. Chem. Phys. 2009, 131, 091101. (25) Lutnæs, O. B.; Teale, A. M.; Helgaker, T.; Tozer, D. J.; Ruud, K.; Gauss, J. J. Chem. Phys. 2009, 131, 144104. (26) Wilson, P. J.; Tozer, D. J. Chem. Phys. Lett. 2001, 337, 341–348. (27) Zhao, Q.; Morrison, R. C.; Parr, R. G. Phys. ReV. A 1994, 50, 2138– 2142. (28) Wilson, P. J.; Tozer, D. J. J. Mol. Struct. 2002, 602-603, 191– 197. (29) Allen, M. J.; Keal, T. W.; Tozer, D. J. Chem. Phys. Lett. 2003, 380, 70–77. (30) Patchkovskii, S.; Autschbach, J.; Ziegler, T. J. Chem. Phys. 2001, 115, 26–42. (31) Teale, A. M.; Tozer, D. J. Chem. Phys. Lett. 2004, 383, 109–114. (32) Hieringer, W.; Della Sala, F.; Go¨rling, A. Chem. Phys. Lett. 2004, 383, 115–121. (33) Arbuznikov, A. V.; Kaupp, M. Chem. Phys. Lett. 2004, 386, 8–16. (34) Arbuznikov, A. V.; Kaupp, M. Chem. Phys. Lett. 2004, 391, 16– 21. (35) Cohen, A. J.; Wu, Q.; Yang, W. Chem. Phys. Lett. 2004, 399, 84– 88. (36) Arbuznikov, A. V.; Kaupp, M. Int. J. Quantum Chem. 2005, 104, 261–271. (37) Teale, A. M.; Tozer, D. J. Phys. Chem. Chem. Phys. 2005, 7, 2991– 2998. (38) Lutnæs, O. B.; Teale, A. M.; Helgaker, T.; Tozer, D. J. J. Chem. Theory Comput. 2006, 2, 827–834.
Peach et al. (39) Teale, A. M.; Cohen, A. J.; Tozer, D. J. J. Chem. Phys. 2007, 126, 074101. (40) Yang, W.; Wu, Q. Phys. ReV. Lett. 2002, 89, 143002. (41) Wu, Q.; Yang, W. J. Theor. Comput. Chem. 2003, 2, 627–638. (42) Hirata, S.; Ivanov, S.; Grabowski, I.; Bartlett, R. J.; Burke, K.; Talman, J. D. J. Chem. Phys. 2001, 115, 1635–1649. (43) Go¨rling, A. Phys. ReV. Lett. 1999, 83, 5459–5462. (44) Heβelmann, A.; Go¨tz, A. W.; Della Sala, F.; Go¨rling, A. J. Chem. Phys. 2007, 127, 054102. (45) Heaton-Burgess, T.; Bulat, F. A.; Yang, W. Phys. ReV. Lett. 2007, 98, 256401. (46) Bulat, F. A.; Heaton-Burgess, T.; Cohen, A. J.; Yang, W. J. Chem. Phys. 2007, 127, 174101. (47) Heaton-Burgess, T.; Yang, W. J. Chem. Phys. 2008, 129, 194102. (48) Staroverov, V. N.; Scuseria, G. E.; Davidson, E. R. J. Chem. Phys. 2006, 124, 141103. (49) Staroverov, V. N.; Scuseria, G. E.; Davidson, E. R. J. Chem. Phys. 2006, 125, 081104. (50) Harriman, J. E. Phys. ReV. A 1986, 34, 29–39. (51) Harriman, J. E.; Hoch, D. E. Int. J. Quantum Chem. Chem. 1997, 63, 111–119. (52) Huzinaga, S. Approximate Atomic Functions; University of Alberta: Edmonton, 1971. (53) Kutzelnigg, W.; Fleischer, U.; Schlindler, M. NMRsBasic Principles and Progress; Springer: Heidelberg, 1990; Vol. 23. (54) Wilson, P. J.; Amos, R. D.; Handy, N. C. Chem. Phys. Lett. 1999, 312, 475–484. (55) Handy, N. C.; Tozer, D. J. Mol. Phys. 1998, 94, 707–715. (56) Kongsted, J.; Aidas, K.; Mikkelsen, K. V.; Sauer, S. P. A. J. Chem. Theory Comput. 2008, 4, 267–277. (57) Becke, A. D. Phys. ReV. A 1988, 38, 3098–3100. (58) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (59) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623–11627. (60) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785–789. (61) Keal, T. W.; Tozer, D. J. J. Chem. Phys. 2003, 119, 3015–3024. (62) Lindh, R.; Malmqvist, P.-Å.; Gagliardi, L. Theor. Chem. Acc. 2001, 106, 178–187. (63) London, F. J. Phys. Radium 1937, 8, 397–409. (64) Hameka, H. F. Mol. Phys. 1958, 1, 203–215. (65) Pople, J. A. Mol. Phys. 1958, 1, 175–180. (66) Dalton, A Molecular Electronic Structure Program, Release 2.0, 2005. http://www.kjemi.uio.no/software/dalton/dalton.html, see http:/ /www.kjemi.uio.no/software/dalton/dalton.html. (67) Fermi, E.; Amaldi, E. Mem. Cl. Sci. Fis. Mat. Nat. Acad. Italia 1934, 6, 117–149. (68) Bu¨hl, M. Theor. Chem. Acc. 2002, 107, 336–342. (69) Wilson, P. J.; Amos, R. D.; Handy, N. C. Mol. Phys. 1999, 97, 757–768.
JP102465X
