Spin–spin coupling constants and triplet instabilities in Kohn–Sham theory

(1)

Vol. 108, Nos. 19–20, 10–20 October 2010, 2579–2590

INVITED ARTICLE

Spin–spin coupling constants and triplet instabilities in Kohn–Sham theory

Ola B. Lutnæs^a, Trygve Helgaker^a* and MichalJaszun´ski^b

aDepartment of Chemistry, Centre for Theoretical and Computational Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway;^bInstitute of Organic Chemistry,

Polish Academy of Sciences, Kasprzaka 44, Warsaw, 01224, Poland (Received 31 May 2010; final version received 24 July 2010)

Indirect nuclear spin–spin coupling constants calculated using restricted Hartree–Fock theory are unreliable since the usually dominant Fermi-contact (FC) contribution arises from a triplet perturbation of the electronic system, poorly described in the Hartree–Fock theory – in particular, at geometries close to the onset of triplet instabilities. These problems are usually but not invariably overcome in Kohn–Sham theory, which typically provides good spin–spin coupling constants. We here examine the sensitivity of spin–spin coupling constants to triplet instabilities in Kohn–Sham and Hartree–Fock theories by correlating the quality of the spin–spin coupling constants and the quality of the lowest triplet excitation energy for a number of small molecules. In general, the FC contributions are most stable in the local density approximation (LDA) and slightly less stable in the generalised gradient approximation (GGA); on the other hand, GGA coupling constants are usually more accurate than LDA constants. Importantly, although hybrid theory often gives better results than the GGA theory, it is also more susceptible to triplet instabilities (inheriting some of the problems of the Hartree–Fock theory) and therefore less reliable than the GGA theory for spin–spin coupling constants. For calculations of spin–spin coupling constants, we recommend the Perdew–Burke–Ernzerhof GGA exchange-correlation functional, which provides a good compromise of accuracy and robustness.

Keywords: spin–spin coupling constants; triplet excitations; triplet instabilities; Hartree–Fock theory

1. Introduction

Many molecular properties describing the response of a molecular system to perturbations, such as a nuclear magnetic moment or an externally applied magnetic field can be calculated as derivatives of the electronic energy with respect to these perturbations. In such calculations, the quality of the theoretical values strongly depends on the approximation applied in describing the electronic system. In particular, the underlying computational model must have the flexi- bility not only to describe the electronic ground state accurately; it must also be sufficiently flexible to describe the changes induced in the ground state by the perturbations. Often, these changes are conve- niently discussed in terms of the excited states of the system. Viewed from this perspective, the electronic- structure model that is applied to the study of molecular properties must not only provide an accurate representation of the electronic ground state but also provide a balanced description of the excited states that contribute to the property of interest.

Thus, any computational model that provides an unbalanced description of the ground state and the most important (lowest) excited states of a given symmetry will also provide a poor description of perturbations of that symmetry.

The calculation of indirect spin–spin coupling constants in nuclear magnetic resonance (NMR) spec- troscopy at the restricted Hartree–Fock (RHF) level of theory provides a good illustration of this point. At the non-relativistic level of theory, four separate terms contribute to the spin–spin coupling constants: the diamagnetic spin–orbit (DSO) term, the paramagnetic spin–orbit (PSO) term, the spin–dipole (SD) term, and the Fermi-contact (FC) term [1]. Whereas the DSO and PSO terms are usually well-described by the RHF theory, the SD term and the dominant FC contributions are poorly described and may sometimes be in error by several orders of magnitude. The reason for this behaviour is the inability of the RHF theory to provide a balanced representation of the ground state and the lowest triplet state. Indeed, it is often observed

*Corresponding author. Email: [email protected]

ISSN 0026–8976 print/ISSN 1362–3028 online

!2010 Taylor & Francis DOI: 10.1080/00268976.2010.513344 http://www.informaworld.com

Downloaded By: [Aarhus Universitets Biblioteker] At: 14:45 5 November 2010

(2)

that, in the RHF theory, the lowest triplet state occurs below the singlet ground state – in particular, for molecules with stretched bonds. As a consequence of such triplet instabilities or near instabilities, molecular properties that involve triplet perturbations are typically poorly described at the RHF level of theory.

RHF triplet instabilities and their consequences have been thoroughly discussed in the literature – see, for example [2–5]. We here note that such triplet instabilities affect not only the RHF results but also the results obtained using correlated models that employ the RHF state as the zero-order model. For example, coupled- cluster (CC) calculations of indirect spin–spin coupling constants are reliable if carried out only in a special manner, without orbital relaxation [6].

In contrast, computational methods based on density-functional theory (DFT) are less prone to triplet instabilities and offer a substantial improvement over the RHF theory for the calculation of spin–spin coupling constants [7]. Indeed, the results obtained with restricted Kohn–Sham (RKS) theory are often close to those obtained with highly accurate but more expensive correlated methods. However, although the RKS theory normally produces fairly accurate predic- tions of spin–spin coupling constants, there are cases where RKS theory also shows symptoms of triplet instabilities, as recently demonstrated for the near biradicalo-benzyne molecule [8]. In this study, two key observations were made concerning triplet-instability problems in calculations of spin–spin coupling constants: first, they strongly depend on which exchange–

correlation functional is used; second, they are much less pronounced at optimised geometries than at experimental geometries. Therefore, to avoid triplet- instability problems, indirect spin–spin coupling constants should be evaluated at the optimised geometry and the exchange–correlation functional should be chosen with care.

Although a low-lying triplet state or a multiple bond may signal triplet-instability problems, it is in general not possible to know in advance for which molecules triplet instabilities will give rise to poor spin–spin coupling constants in RKS calculations. It is therefore desirable to gain more insight into the problem of RKS triplet instabilities and, in particular, how these affect the calculation of spin–spin coupling constants. This is the focus of this article, where we explore this topic by calculating FC contributions for five molecules (carbon monoxide, nitrogen, acetylene, ethylene and ethane) at various non-equilibrium geometries using a representative set of RKS exchange–correlation functionals; for comparison, we also apply RHF theory.

2. Theory and background

2.1. Triplet instabilities in RHF theory

In the RHF theory, the wave function is constrained to a certain spin symmetry – in the closed-shell case, where all orbitals are doubly occupied, the RHF wave function has singlet spin symmetry. Apart from pro- viding computational simplification, the purpose of this constraint is to enforce properties of the exact wave function on the approximate Hartree–Fock wave function. However, as a result of enforcing spin symmetry in the RHF theory, situations may arise when there exists an unrestricted Hartree–Fock (UHF) wave function with a lower energy than the RHF singlet wave function. The RHF wave function is then said to be triplet unstable. Triplet instability is a well- known phenomenon that has been investigated in the literature on numerous occasions – see, for example Ref. [2].

As an elementary illustration, we consider the dissociation of H2, for which the potential-energy curves at different levels of theory are displayed in Figure 1 (with the nuclear–nuclear repulsion removed).

The calculations have been carried out in an atomic- orbital (AO) basis consisting of one Slater orbital on each atom A and B: 1sA=BðrÞ ¼exp $rA=B

! "

= ffiffiffi p!

, where r_Aand r_B are the distances between the electron and nuclei A and B, respectively. The corresponding bonding and antibonding molecular orbitals (MOs) are given by 1"_g/u(r)¼N_g/u[1s_A(r)%1s_B(r)], whereN_g/u is the normalisation constant. From these MOs, we have constructed the closed-shell RHF wave functions

1!^þ_gðg²Þ

$$

$ E

¼ 1"g1""g

$$ %

, ð1Þ

1!^þ_gðu²Þ

$$

$ E

¼j1"u1""ui, ð2Þ the open-shell RHF wave functions

$$³!^þ_uðguÞ%

¼&$$1"g1"u% , 1

ffiffiffi

p2!$$1"_g1""u%

þ$$1""g1"u%"

,$$1""g1""ui '

, ð3Þ

1!^þ_uðguÞ

$$ %

¼ 1 ffiffiffi p2 1"g""u

$$ %

$ 1""g"u

$$ %

! "

, ð4Þ

and the ground-state full configuration–interaction (FCI) wave function and its orthogonal counterpart of the form

1!^þ_gðFCIÞ

$$

$ E

¼Cg 1"g1""g

$$ %

þCuj1"u1""ui: ð5Þ In the following, we review the behaviour of these wave functions as the bond distance increases from R¼0

(3)

(helium atom) through R¼1.4a0 (hydrogen molecule at equilibrium) to R¼ 1 (two separate hydrogen atoms).

As is well-known, the RHF ground-state configu- rationj¹!^þ_gðg²Þiprovides a good representation of the electronic system at short bond lengths, correlating with the ¹S(1s²) ground-state of helium. However, in the limit of complete dissociation,j¹!^þ_gðg²Þiprovides a mixed ionic and covalent description ð1s²_Aþ1s²_BÞþ

1sA1sBþ1sB1sA

ð Þ. Thus, it gives an increasingly poor representation of the ground-state wave function upon dissociation, as seen from the following FCI wave functions:

1!^þ_gðFCIÞ

$$

$ E

He¼1:00 1"g1""g

$$ %

$0:03 1"j u1""ui, ð6Þ

1!^þ_gðFCIÞ

$$

$ E

H2

¼0:99 1"g1""g

$$ %

$0:11 1"j u1""ui, ð7Þ

1!^þ_gðFCIÞ

$$

$ E

HþH¼0:71 1"g1""g

$$ %

$0:71 1"j u1""ui: ð8Þ

Thus, in the dissociation limit, covalent dissociation 1sA 1sBþ1sB1sA requires a linear combination of bonding and antibonding configurations. By contrast, the lowest triplet state of the hydrogen molecule is well-described by$$³!^þ_uðguÞ%

at all internuclear separations correlating with the ³P(1s2p) state of helium in the limit of a united atom and with the covalent ground state in the limit of complete dissociation, with spatial wave function 1s_A1s_B$1s_B1s_A.

As a result of the poor ground-state dissociation but correct triplet-state dissociation in the RHF theory, the lowest triplet excitation energy of the hydrogen molecule is poorly described, even becoming zero at large bond distances, where the RHF wave function becomes triplet unstable. In our simple minimal-basis calculations in Figure 1, the onset of the triplet instability occurs at a bond distance of about 3.13a0, almost twice the equilibrium bond distance of 1.67a0in this basis. For other systems and other basis sets, this onset may occur closer to the equilibrium bond distance and even at bond distances shorter than the equilibrium distance. Since the quality of calculated triplet properties depend critically on the quality of the Figure 1. The electronic energy without nuclear repulsion (hartree) of H₂plotted against bond distance (bohr) at different levels of theory in a minimal STO basis. Black lines represent closed-shell RHF wave functions; dashed lines represent open-shell RHF wave functions and grey lines represent FCI wave functions.

(4)

triplet excitation energies and, in particular, the lowest triplet excitation energy, we conclude that the RHF theory is ill-suited to the study of such properties.

For some systems, where the onset of the triplet instability is sufficiently far removed from the equilibrium bond distance to cause problems, it may provide reasonable accuracy. For other systems, the onset may occur before the equilibrium geometry, leading to physically meaningless results. In short, the RHF theory is unsuitable for the prediction of triplet properties, such as indirect spin–spin coupling constants.

For an accurate description of triplet properties and, in particular, indirect spin–spin coupling constants, we must admit more determinants into the description, as done, for example, in configuration–

interaction (CI) theory [9], in multiconfigurational self- consistent field (MCSCF) theory [10], in CC theory [11,12] and in the second-order polarisation propaga- tor approximation (SOPPA) [13,14]. In this manner, a balanced treatment of singlet and triplet states may be achieved, removing problems associated with triplet instabilities. Alternatively, we may resort to DFT.

Indeed, the main focus of this article is the description of spin–spin coupling constant in RKS theory. For a recent review of the application of DFT as well as wave-function based methods to the study of spin–spin coupling constants, see Ref. [15].

2.2. Triplet instabilities in Kohn–Sham theory

In RKS theory, no attempt is made at representing the correlated wave function of the interacting electrons.

Rather, we construct a determinantal wave function

#¼detj#1,#2,. . .j for a fictitious system of non-

interacting electrons of the same density $ as the interacting system. In principle, an exact description of the electronic system can be achieved in this manner.

In practice, the quality of the description depends critically on the quality of the exchange–correlation functional Exc[$], from which the effects of exchange and correlation are extracted directly from the electron density. Over the years, many approximate exchange–

correlation functionals have been developed. In the local-density approximation (LDA), the functional is obtained by applying locally relations that are valid globally for the uniform electron gas. The LDA model provides a surprisingly good description of molecular systems – for example, its exchange–correlation hole remains localised upon dissociation of H2, suggesting that the LDA model provides a balanced treatment of exchange and static correlation, which

is important for the calculation of spin–spin coupling constants.

Higher accuracy is achieved in the generalised gradient approximation (GGA), where density- gradient corrections are introduced in the exchange and correlation functionals. However, GGA exchange–correlation functionals are often constructed in a semi-empirical manner, by concentrating on thermochemistry and molecular equilibrium properties. In practice, GGA exchange–correlation functionals also typically provide improved spin–spin coupling constants, compared with LDA. However, we shall later see that these methods are more affected by triplet instabilities than the LDA method and are therefore less robust for the calculation of spin–spin coupling constants.

In hybrid RKS theory, orbital-dependent exact exchange is introduced in the description of the electronic system. However, it is worth noting that, if orbital-independent LDA and GGA exchange is fully replaced by orbital-dependent exact exchange, the performance of RKS theory usually deteriorates rather than improves. This behaviour occurs since GGA functionals typically benefit from error cancel- lations that are destroyed by the introduction of exact exchange [16]. In particular, use of exact exchange necessitates an accurate description of static correlation, which is neglected in, for example, the Lee–Yang–

Parr (LYP) GGA correlation functional [17,18]. This observation is particularly relevant for triplet properties, as the neglect of static correlation leads to triplet instabilities, such as those observed in the RHF theory.

In practice, exact exchange is usually introduced with a weight of about 20% relative to orbital-independent exchange. We then expect triplet-instability problems to be less severe than in the RHF theory but still noticeable. For a discussion of triplet instabilities in RKS theory, see Ref. [7].

2.3. Triplet instabilities and second-order property divergence

In the random-phase approximation (RPA) and in time-dependent DFT (TDDFT), excitation energies! are obtained from the solution of the generalised eigenvalue problem [4,5,19]

A B B A

( )

X Y ( )

¼! I 0 0 $I

( )

X Y ( )

, ð9Þ whereXandYcontain the orbital rotations associated with the excitation and where A and B are here

(5)

assumed to be real. By an orthogonal transformation, the Hessian in Equation (9) may be expressed in block diagonal form

H¼ A B B A

( )

¼ 1 ffiffiffi

p2 I I I $I

( ) AþB 0

0 A$B

( )

1ffiffiffi

p2 I I I $I

( )

, ð10Þ where the diagonal blocksA$BandAþBare the real and imaginary parts of the electronic Hessian, respectively. In spin-adapted RHF and RKS theories, the Hessian elements of spin (‘,m) are given by (using indices i and j for occupied orbitals and a and b for virtual orbitals)

½A$B(^‘,mai,bj

¼%_ij%_abð"_a$"_iÞ þ4%‘0gaibj$&ðgabijþgajibÞ þG^‘,m_ai,_bj,

ð11Þ

½AþB(^‘,mai,bj¼%ij%abð"a$"iÞ $&ðgabijþgajibÞ, ð12Þ where "_p are orbital energies, g_pqrs are two-electron integrals in Mulliken notation, & is the weight of orbital-dependent exchange and G^‘,_pqrs^m is the exchange–

correlation contribution. In the RHF theory,&¼1 and G^‘,_ai,^m_bj¼0; in pure RKS,&¼0 andG^‘,_pqrs^m 6¼0, its explicit form depending on the choice of exchange–correlation functional. In hybrid theories, some proportion &6¼0 (typically 20%) of orbital-dependent exchange is also included in DFT. For imaginary rotations, there is no difference between the singlet and triplet Hessian elements; by contrast, for real rotations, the singlet and triplet elements differ in the Coulomb contribution

%‘0g_aibj (which is absent in the triplet case) and in the exchange–correlation contribution. The imaginary Hessian is particularly simple in pure RKS theory, being a diagonal. Finally, we note that we need only consider triplet Hessian elements with spin projection m¼0, to which the other components are related by the Wigner–Eckart theorem.

Due to the special form of the eigenvalue problem in Equation (9), its complex eigenvalues occur in pairs

%!, with eigenvectors (X,Y)^T and (Y,X)^T. In the

following, we briefly review the relationship of these eigenvalues to those of the electronic Hessian in Equation (10), which are equal to the eigenvalues of A%B; for a more detailed treatment, see Ref. [7]. This relationship is particularly relevant here since time- independent second-order properties, such as the FC contribution to spin–spin coupling constants are given by expressions of the form F^yH^$¹ F. Stable RHF or

RKS solutions are characterised by a positive definite Hessian, whereas unstable solutions appear when the Hessian has one or more negative eigenvalues. We are here particularly concerned about real triplet instabilities, arising from negative eigenvalues in the triplet part of the real HessianA$B.

Performing the orthogonal transformation in Equation (10), we obtain the equivalent eigenvalue problem

A%B

ð ÞðX%YÞ ¼!ðX)YÞ, ð13Þ which upon multiplication byA)B and some simple rearrangement yields

A)B

ð ÞðA%BÞðX%YÞ ¼!²ðX%YÞ: ð14Þ Finally, if A%B and therefore (A%B)^1/2 is positive definite and symmetric, then this non-symmetric eigenvalue problem is equivalent to the following symmetric eigenvalue problem

A%B

ð Þ¹⁼²ðA)BÞðA%BÞ¹⁼²ðX_%%Y_%Þ ¼!²ðX_%%Y_%Þ, ð15Þ where we have introduced X_%¼(A%B)^1/2 X and similarly for Y_%. Thus, if one of the Hessian blocks A%B (not necessarily both) is positive definite, it follows that !²are real and hence that the excitation energies!are real or pure imaginary. WhenA%Bare both positive definite, then all excitation energies are positive. Conversely, ifA%Bare both indefinite, then the excitation energies may be complex. Finally, by setting !¼0 in Equation (9), we see that zero excitation values occur if and only if the Hessian is singular. Zero excitation values are therefore associated with divergent second-order propertiesF^yH^$¹F.

3. Computational details

To investigate the effect of triplet instabilities on the indirect nuclear spin–spin coupling constants, we have calculated these at different geometries for the carbon monoxide, nitrogen, acetylene, ethylene and ethane molecules. Since we are here concerned with the effect of triplet instabilities on the computed coupling constants, we have restricted our attention to the (usually) dominant FC triplet contribution, ignoring the small SD triplet contribution and the DSO and PSO singlet contributions. Since the FC contribution is isotropic, there is no need to consider individual tensor components. In addition to the FC contribution to the coupling constant, we have in all cases also calculated the lowest triplet excitation energy of the system, using linear response theory.

(6)

In this study, we have considered a number of popular exchange–correlation functionals, representative of three distinct approaches to the calculation of the exchange–correlation energy in molecular systems: the fully local LDA approach [20]; the gradient-corrected local GGA approach with the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional [21], the Becke–Lee–Yang–Parr (BLYP) functional [17,18] and the Keal–Tozer KT2 functional [22]; and the hybrid approach with the PBE0 functional [23], the Becke–3-parameter–Lee–Yang–Parr (B3LYP) functional [18,24], the B97-2 functional [25]

and the range-separated Coulomb-attenuated-method B3LYP (CAMB3LYP) functional [26,27]. All hybrid calculations were performed in the conventional fully variational manner, without using optimised effective potentials. For comparison, we have also applied the Hartree–Fock model but have not considered correlated wave-function models.

In all calculations, we have used Dunning’s core–valence one-electron basis set cc-pCVTZ [28], previously shown to give good results for indirect nuclear spin–spin coupling constants, the inclusion of core functions being essential for high accuracy of the FC contribution to the coupling constants. All calculations were performed with the Dalton quantum- chemistry program [29]. All reported spin–spin coupling constants are for the ¹H, ¹³C, ¹⁴N and ¹⁷O isotopomers.

4. Results and discussion 4.1. The CO and N2molecules

The multiple !-bonded CO and N2 molecules are typical examples of molecules where triplet instabilities arise upon bond stretching. Indeed, the effect of triplet instabilities on the indirect spin–spin coupling constants of these molecules has previously been studied by Auer and Gauss [6], using CC theory. These molecules are unusual in the sense that the total spin–spin coupling constant (at equilibrium) is not dominated by the FC contribution; moreover, stan- dard exchange–correlation functionals typically give poor spin–spin coupling constants for these molecules.

However, for the study of the influence of triplet instabilities on the spin–spin coupling constant, this does not matter.

In Figure 2(a), we have plotted the lowest triplet excitation energy of CO against the bond length, for the nine levels of electronic-structure theory considered in this study. With regard to the triplet instability, we discern three broad groups of methods, in agreement with the observations made in Ref. [7]. The LDA, PBE

and BLYP methods are more robust, the instability occurring at bond distances longer than 1.7 A˚, corresponding to elongations of 50% or more relative to equilibrium. The second group consists of the hybrid methods (B3LYP, PBE0, CAMB3LYP and B97-2) and the KT2 GGA method, for which the triplet instability occurs at a distance of about 1.6 A˚ (40% elongation).

Finally, the RHF theory gives a triplet instability already at 1.3 A˚, close to the experimental bond length of 1.1283 A˚.

These observations can be rationalised in terms of the proportion of exact (orbital-dependent) exchange used in the description of the electronic system. First, the RHF theory, which treats exchange exactly but neglects electron correlation, is particularly unstable.

Next, stability is significantly improved at the hybrid level of theory, which includes electron correlation but retains some proportion of orbital-dependent exact exchange (20%, 21% and 25% in the B3LYP, B97-2, and PBE0 exchange functionals, respectively). Finally, highest triplet stability is achieved by those methods that avoid orbital-dependent exact exchange alto- gether, describing exchange instead in a local or semi-local manner at the LDA and GGA levels of theory. It is worth noting that the LDA triplet excitation energy nowhere vanishes in the range of plotted bond distances, between 1.10 and 1.78 A˚. This simple exchange–correlation functional therefore appears particularly stable towards triplet instabilities, reflecting a balanced description of exchange and correlation in this model.

The only exception to the general triplet-instability pattern observed above is the KT2 GGA functional, which is as unstable as the hybrid functionals but makes no use of exact exchange. We have no expla- nation for this behaviour but speculate that it may, in some manner, be related to the fact that the KT2 functional has been developed in a different manner than the other GGA functionals in our study, by adjusting it to reproduce accurate exchange–

correlation potentials.

Turning our attention to the calculated FC contributions to the spin–spin coupling constants plotted in Figure 2(b), we note the expected divergence of the coupling constants at the onset of the triplet instability.

Thus, upon bond stretching, divergence is first observed for the RHF model, much later for the KT2 model, followed immediately by the hybrid (B97-2, PBE0, CAMB3LYP and B3LYP) models.

Finally, divergence is observed also for the local and semi-local models (BLYP, PBE0 and LDA). However, it should be emphasised that triplet stability upon bond stretching does not ensure a high accuracy in the calculated spin–spin coupling constants and excitation

(7)

energies at the equilibrium bond distance. In particular, the LDA model typically performs significantly poorer than the GGA and hybrid methods for spin–

spin coupling constants at the equilibrium geometry.

For the N2 molecule, we have plotted the lowest triplet excitation energy and the FC contribution to the spin–spin coupling constant in Figure 3. Based on our findings for CO, we have for this and the remaining molecules in this study restricted ourselves to the following representative levels of theory: the LDA theory, the GGA theory with the PBE exchange–

correlation functional, hybrid theory with the B3LYP functional and the RHF theory. The PBE functional

was selected as the best GGA functional with respect to triplet instabilities for the CO molecule, while the B3LYP functional was chosen for its wide use. As seen from Figure 3, the behaviour of N2is similar to that of the iso-electronic CO molecule. The RHF model is particularly poor, with the onset of triplet instability very close to the equilibrium bond distance. All RKS methods provide a dramatic improvement upon the RHF method, the hybrid B3LYP method diverging before the PBE GGA method. Again, the LDA model is most stable towards triplet instabilities and thus provides the best description of the indirect nuclear 0.0

1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0

R(CO) (Å)

Triplet excitation energies (eV)

LDA PBE BLYP KT2 PBE0 B97-2 B3LYP CAMB3LYP HF

–20 –10 0 10 20 30 40 50 (b) (a)

1.05 1.15 1.25 1.35 1.45 1.55 1.65 1.75 1.85

R(CO) (Å)

Spin–spin coupling constants (Hz)

LDA PBE BLYP KT2 PBE0 B97-2 B3LYP CAMB3LYP HF

Figure 2. (a) The lowest triplet excitation energy in CO as function of bond length. (b) The FC contribution to the indirect nuclear spin–spin coupling constant in CO as a function of bond length. The experimental equilibrium bond length is 1.1283 A˚.

(8)

spin–spin coupling constant in N2over a wide range of internuclear separations.

4.2. The C2H2, C2H4and C2H6molecules

In addition to bond stretching in CO and N2, we have studied three different effects in hydrocarbons. For C2H2, we have studied the symmetric bending mode (change of two HCC bond angles with the atoms in a planar cis arrangement). For C₂H₄, which provides a classical example of triplet instabilities in the RHF

theory [2], we considered twist of one –CH2 group.

Finally, for C2H6, we stretched the CC bond. In each molecule, we varied only one parameter, keeping all the other geometry parameters (including the geometry of the –CHn groups) fixed at the equilibrium or experimental values. For all molecules, we calculated the lowest triplet excitation energy and the FC contributions to the full set of spin–spin coupling constants, using the same four levels of theory as for N2.

The results for C2H2, C2H4and C2H6are shown in Figures 4, 5 and 6, respectively. For each molecule, we have plotted the lowest triplet excitation energy and 0.0

2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 (a)

(b)

0.70 0.90 1.10 1.30 1.50 1.70

LDA PBE B3LYP HF

–10 –5 0 5 10 15 20 25 30

0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50 1.60 1.70

R(NN) (Å) R(NN) (Å)

Figure 3. (a) The lowest triplet excitation energy of N2as a function of bond length. (b) The FC contribution to the indirect nuclear spin–spin coupling constant and its asymptote (dashed line) of N2as a function of bond length. The experimental bond length is 1.0953A˚. Colour scheme: LDA – dark blue, PBE – purple, B3LYP – yellow, RHF – turquoise.

(9)

one of the FC spin–spin coupling constants as a function of the geometry change, noting that all spin–

spin coupling constants are affected in the same manner. In particular, although the geometry of

$CHn groups does not change, the one- and two- bond coupling constants within these groups and the coupling constants across the carbon–carbon bond behave similarly, diverging at the same stage in the variation of the molecular geometry. This was expected as the lowest triplet excitation energy is a global molecular property whose poor description affects all coupling constants, even though its value may reflect the change of a single geometry parameter related to one bond.

For these molecules, we observe the same trends among the different levels of theory as for the N₂and CO molecules. The triplet instability appears first for the RHF theory (for C2H4, even at equilibrium) followed by the hybrid B3LYP functional, the PBE GGA functional and finally by the LDA functional.

5. Conclusions

We have studied the effect of triplet instabilities on indirect nuclear spin–spin coupling constants calculated in the RHF and RKS theories. In the RHF theory, triplet instabilities occur as a result of an

0.0 1.0 2.0 3.0 4.0 5.0 (a) 6.0

(b)

f(HCC) (degree)

Triplet excitation energies (eV)

100 150 200 250 300 350 400 450 500

0 10 20 30 40 50 60 70 80 90

f(HCC) (degree)

Spin–spin coupling constants (Hz)

Figure 4. (a) The lowest triplet excitation energy of C2H2as a function of the HCC bond angle#. (b) The FC contribution to the¹JCCcoupling constant and its asymptote (dashed line) of C2H2as a function of#. Colour scheme: LDA – dark blue, PBE – purple, B3LYP – yellow, RHF – turquoise.

(10)

imbalance in the description of the singlet ground state and the lowest triplet excited state as the molecule dissociates. Whereas the triplet state is described in a qualitatively correct manner at all geometries, the description of the ground state deteriorates upon distortion resulting in a triplet instability, with imaginary excitation energies and divergences in triplet properties, such as the FC contribution to indirect nuclear spin–spin coupling constants. In wave-function based methods, triplet instabilities are removed by allowing more configurations into the description – for example, using CI, MCSCF, SOPPA or CC theories.

In a sense, we may interpret triplet instabilities as arising from an imbalance in the treatment of exchange and correlation: in the RHF theory, triplet instabilities follow since exchange is exactly treated, whereas correlation is neglected.

In exact Kohn–Sham theory, exchange and correlation are treated exactly and there are no triplet instabilities. In approximate Kohn–Sham theory, triplet instabilities may arise following a lack of balance in the treatment of exchange and correlation. This balance is best maintained in the simple LDA model, which is less affected by triplet instabilities than all 0.0

0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 (a)

(b)

q(HCCH) (degree)

–30 –20 –10 0 10 20 30 40 50

0 10 20 30 40 50 60 70 80

0 10 20 30 40 50 60 70 8

q(HCCH) (degree)

0

Figure 5. (a) The lowest triplet excitation energy of C₂H₄as a function of the dihedral (twist) angle'. (b) The FC contribution to the²JCHcoupling constant and its asymptote (dashed line) of C2H4as a function of'. Colour scheme: LDA – dark blue, PBE – purple, B3LYP – yellow, RHF – turquoise.

(11)

other levels of theory investigated here. In the GGA theory, the improvements made in the treatment of exchange and correlation by including gradient corrections lead to important improvements in ener- getics and many molecular equilibrium properties but create an imbalance in the treatment of exchange and correlation, making this level of theory more susceptible to triplet instabilities. This imbalance becomes particularly evident in hybrid theory, where some proportion of exact exchange is introduced into the description. As a result, triplet instabilities arise more easily in the hybrid RKS theory than in the GGA

theory, sometimes affecting spin–spin coupling constants calculated at the equilibrium geometry.

Although the LDA method is the most robust method for the calculation of indirect nuclear spin–

spin coupling constants, being less affected by triplet instabilities than the other methods studied here, it is in general not the most accurate one. In previous inves- tigations, it has been found that the GGA and hybrid RKS models typically provide similar accuracies, with a slight preference for hybrid theory. However, given that the GGA models are significantly less affected by triplet instabilities than are hybrid models, the best 0.0

2.0 4.0 6.0 8.0 10.0 (a)

(b)

R(CC) (Å)

50 70 90 110 130 150 170 190 210 230 250

1.00 1.20 1.40 1.60 1.80 2.00 2.20

R(CC) (Å)

Figure 6. (a) The lowest triplet excitation energy of C₂H₆as a function of the R(CC) bond length. (b) The FC contribution to the¹J_CHcoupling constant and its asymptote (dashed line) of C₂H₆as a function of R(CC). The equilibrium CC bond length is 1.5351 A˚. Colour scheme: LDA – dark blue, PBE – purple, B3LYP – yellow, RHF – turquoise.

(12)

compromise between accuracy and robustness is pro- vided by the GGA theory. Our recommendation for the RKS calculation of indirect spin–spin coupling constants in molecular systems is the PBE exchange–

correlation functional. Given that triplet instabilities are typically less pronounced close to the optimised equilibrium geometry (at a given level of theory) than at the experimental structure, we recommend that indirect nuclear spin–spin coupling constants are calculated at optimised (rather than experimental) geometries, as observed in Ref. [8].

Acknowledgements

This work was supported by the Norwegian Research Council through the CoE Centre for Theoretical and Computational Chemistry (Grant No. 179568/V30).

We would further like to acknowledge the NOTUR com- puting facilities which have been used to conduct the calculations presented in this article.

Note

All figures can be viewed in colour online.

References

[1] N.F. Ramsey, Phys. Rev.91, 303 (1953).

[2] J. Cˇizˇek and J. Paldus, J. Chem. Phys. 47, 3976 (1967).

[3] R. Seeger and J.A. Pople, J. Chem. Phys. 66, 3045 (1977).

[4] R. McWeeny, Methods of Molecular Quantum Mechanics, 2nd ed.(Academic Press, London, 1992).

[5] P. Jørgensen and J. Simons, Second Quantization-Based Methods in Quantum Chemistry. (Academic Press, New York, 1981).

[6] A.A. Auer and J. Gauss, Chem. Phys.356, 7 (2001).

[7] R. Bauernschmitt and R. Ahlrichs, J. Chem. Phys.104, 9047 (1996).

[8] T. Helgaker, O.B. Lutnæs and M. Jaszun´ski, J. Chem.

Theory Comput.3, 86 (2007).

[9] A.A. Auer, J. Gauss and M. Pecul, Chem. Phys. Lett.

368, 172 (2003).

[10] O. Vahtras, H. A˚gren, P. Jørgensen, H.J.Aa. Jensen, S.B. Padkjær and T. Helgaker, J. Chem. Phys.96, 6120 (1992).

[11] A.A. Auer and J. Gauss, J. Chem. Phys. 115, 1619 (2001).

[12] S.A. Perera, H. Sekino and R.J. Bartlett, J. Chem. Phys.

101, 2186 (1994).

[13] J. Geertsen and J. Oddershede, Chem. Phys. 90, 301 (1984).

[14] T. Enevoldsen, J. Oddershede and S.P.A. Sauer, Theor.

Chem. Acc.100, 275 (1998).

[15] T. Helgaker, M. Jaszun´ski and M. Pecul, Progr. Nucl.

Magn. Reson. Spectrosc.53, 249 (2008).

[16] O.V. Gritsenko, P.R.T. Schipper and E.J. Baerends, J. Chem. Phys.107, 5007 (1997).

[17] A.D. Becke, Phys. Rev. A38, 3098 (1988).

[18] C. Lee, W. Yang and R.G. Parr, Phys. Rev. B37, 785 (1988).

[19] R. Bauernschmitt and R. Ahlrichs, Chem. Phys. Lett.

256, 454 (1996).

[20] S.H. Vosko, L. Wilk and M. Nusair, Can. J. Phys.58, 1200 (1980).

[21] J.P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev.

Lett.77, 3865 (1996).

[22] T.W. Keal and D.J. Tozer, J. Chem. Phys. 119, 3015 (2003).

[23] C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).

[24] A.D. Becke, J. Chem. Phys.98, 5648 (1993).

[25] P.J. Wilson, T.J. Bradley and D.J. Tozer, J. Chem. Phys.

115, 9233 (2001).

[26] T. Yanai, D.P. Tew and N.C. Handy, Chem. Phys. Lett.

393, 51 (2004).

[27] M.J.G. Peach, T. Helgaker, P. Salek, T.W. Keal, O.B. Lutnæs, D.J. Tozer and N.C. Handy, Phys.

Chem. Chem. Phys.8, 558 (2006).

[28] D.E. Woon and T.H. Dunning Jr, J. Chem. Phys. 98, 1358 (1993).

[29] DALTON, an ab initio electronic structure program, Release 1.2. (2001). http://www.kjemi.uio.no/software/

dalton/dalton.html