Solvent Effects on the Indirect Spin–Spin Coupling Constants of Benzene: The DFT-PCM Approach

Abstract

We present an extension of the Polarizable Continuum Model (PCM) to the calculation of solvent effects on indirect spin–spin coupling constants for Hartree– Fock wave functions and Density Functional Theory. This is achieved by implementing the PCM model for singlet and triplet linear response functions. The new code is used for calculating the solvent effects on the indirect spin–spin coupling constants of benzene. For the ¹J(H¹³C) coupling constants, our calculated solvent shifts are in good agreement with experimental observations when geometry relaxation is taken into account. However, our results do not support the extrapolated gas-phase value for this coupling constant. A new experimentally derived ¹J(H ¹³C) for a vibrating benzene molecule at 300 K is proposed.

1 Introduction

The ab initio calculation of indirect spin–spin coupling constants is a challenging task. Each nuclear spin can in the non-relativistic formulation of the theory as first derived by Ramsey [1] perturb the electron density through 10 different interaction mechanisms (6 spin-dipolar (SD), 3 paramagnetic spin-dipole (PSO) and one Fermi-contact interaction), and a large number of response equations therefore need to be solved for each nucleus. The presence of the Fermi-contact interaction—which depends on the electron density close to the nuclei—puts severe demands on the quality of the basis set used. For this reason, basis set requirements for accurate calculations of indirect spin–spin coupling constants have been given much attention in the literature [2,3,4,5,6], and a few different schemes for obtaining fairly small basis sets yielding accurate spin–spin coupling constants have in recent years been devised [5,6]. The triplet nature of the Fermi-contact and the spin-dipolar operators means that ordinary spin-restricted approaches such as Hartree-Fock fails completely, often giving results that are several order of magnitudes too large as well as having incorrect signs [7,8]. The problem of triplet instabilities is often quite efficiently solved by introducing electron correlation, and several successful studies using multiconfigurational self-consistent field (MCSCF) [7,9,10,11] or coupled-cluster [12,13,14] wave functions have been presented. An alternative approach has used the second-order polarization propagator approximation (SOPPA) [2,15]. In recent years, it has been demonstrated that even in it’s spin-restricted formulation, density functional theory (DFT) is a computationally efficient route to accurate indirect spin–spin coupling constants [16,17]. Although the first DFT implementations used a finite-perturbation approach [17,18,19], often only considering the Fermi contact interaction mechanism, fully analytical implementations of all interaction mechanisms have recently been presented [16,20]. To further complicate the theoretical study of indirect spin–spin coupling constants, the vibrations of the molecular framework have been shown to give significant contributions to the indirect spin–spin coupling constants, often giving contributions that are 5-10% of the magnitude of the purely electronic contributions. Recent developments have however made it possible to obtain quite accurate estimates for the vibrational corrections to the spin–spin coupling constants [21,22], even for fairly large molecules such as benzene [23]. Even if all the factors affecting the theoretical calculation of indirect spin–spin coupling constants now can be addressed for medium-sized molecules [4], a direct comparison of the theoretical results with experimental observations is in general not possible since most experimental NMR studies take place in solution. Although there are continuous improvements in the gas-phase determination of NMR parameters, only a very small number of indirect spin–spin coupling constants have been determined for molecules in the gas phase [24]. It is generally known that spin–spin coupling constants are less influenced by solvent effects than for instance the nuclear shielding constants, but a detailed comparison of theoretical and experimental results requires that the effects of the solvent on the theoretically calculated spin–spin coupling constants are properly understood [25]. Almost all theoretical investigations of indirect spin–spin coupling constants have been limited to gas-phase investigations (see Ref. [4] and references therein for some recent examples). Only one implementation of a dielectric continuum model for studying solvent effects on spin–spin coupling constants have been presented in the literature [26]. It was initially used to study solvent effects on the spin–spin couplings of H₂Se, H₂S, and HCN [26,27], but has in recent years been used for studies on larger molecules, often together with supermolecular models in order to account for the effects of specific solute-solvent interactions [28,29]. The dielectric continuum model presented in Ref. [26] used a spherical shape for the molecular cavity. Although this model can be expected to be adequate for the study of small or nearly spherical molecules, it may not be suitable for solvent studies on larger molecules. In this paper we present a new approach for calculating solvent effects on indirect–spin coupling constants based on the polarizable continuum model (PCM) [30,31,32]. In the PCM model, the molecule is placed in a molecule-shaped cavity in the dielectric medium, the polarization effects on the solvated molecule being introduced through charges on the cavity surface. The implementation is based on our recent extensions of the PCM model to a quadratically convergent scheme for optimizing PCM-MCSCF wave functions [33] and singlet linear response functions [34] for solvated molecules. In this work we extend that implementation to also include DFT reference states and triplet linear response functions. Although we will here only discuss and employ DFT, our implementation also includes Hartree–Fock and MCSCF reference wave functions. The remainder of this paper is organized as follows: In Section 2 we present the theory for triplet linear response functions for a Kohn-Sham based DFT approach in the PCM formalism. We also briefly summarize the mechanisms that determine the indirect spin–spin coupling constants. We then apply this formalism to study the solvent dependence of the indirect spin–spin coupling constants of benzene, and the computational details of this study is summarized in Section 3. The results we have obtained is presented in Section 4, and we summarize our results in Section 5.

2 Theory

The indirect spin–spin coupling parameters describe the indirect coupling of the nuclear magnetic dipoles, mediated by the surrounding electrons. For an isolated molecule, the indirect spin–spin coupling parameters can be determined as second derivatives of the electronic energy with respect to the nuclear magnetic moments [4]. A parallel situation holds for molecules in solution described by PCM. However, in this case the indirect spin–spin coupling parameters are instead determined as second derivatives of the free-energy functional $𝒢$ of the solute-solvent system with respect the nuclear magnetic moments. The nuclear magnetic moments M_K are related to the nuclear spins I_K as

$${\mathbf{M}}_K={\gamma }_Kħ{\mathbf{I}}_K$$

(1)

where γ_K are the nuclear magnetogyric ratios and the normal and reduced nuclear indirect spin–spin coupling constant J_KL and K_KL are given as

$${\mathbf{J}}_{KL}=h\frac{{\gamma }_K}{2\pi }\frac{{\gamma }_L}{2\pi }{\mathbf{K}}_{KL}=h\frac{{\gamma }_K}{2\pi }\frac{{\gamma }_L}{2\pi }\frac{d^2𝒢}{d{\mathbf{M}}_{K}d{\mathbf{M}}_L}$$

(2)

Introducing the presence of the nuclear magnetic moments, the non-relativistic electronic Hamiltonian of the molecular solute has the form

$$\begin{array}{ll}{H}_{|\mathsf{\Psi }>}=& \frac{1}{2}{\displaystyle \sum _i[{\mathbf{p}}_i{+\mathbf{A}}_i]}\cdot [{\mathbf{p}}_i+{\mathbf{A}}_i]+{\displaystyle \sum _i{\mathbf{s}}_i\cdot \mathbf{B}\left({\mathbf{r}}_i\right)}\\ & +V_{en}+V_{ee}+V_{nn}+V\left(\mathsf{\Psi }\right)\end{array}$$

(3)

where p_i is the conjugate momentum of electron i and s_i it’s spin. V_en, V_ee, and V_nn denote the electron-nuclear attraction potential and the electron-electron and nuclear-nuclear repulsion potentials, respectively. The vector potential A(r_i) and the related induction $\mathbf{B}({\mathbf{r}}_i)=\nabla \times \mathbf{A}({\mathbf{r}}_i)$ are given by

$$\mathbf{A}({\mathbf{r}}_i)={\alpha }^2{\displaystyle \sum _K\frac{{\mathbf{M}}_K\times {\mathbf{r}}_{iK}}{r_{iK}^3}}$$

(4)

$$\begin{array}{ll}\mathbf{B}\left({\mathbf{r}}_i\right)& =\frac{8\pi {\alpha }^2}{3}{\displaystyle \sum _K\delta ({\mathbf{r}}_{iK}){\mathbf{M}}_K}\\ & +{\alpha }^2{\displaystyle \sum _K\frac{{3\mathbf{r}}_{iK}{\mathbf{r}}_{iK}^T-r_{iK}{\mathrm{I}}_3}{r_{iK}^5}{\mathbf{M}}_K}\end{array}$$

(5)

where α is the fine-structure constant, r_iK is the position of electron i relative to nucleus K and I₃ is the 3 × 3 unit matrix. The potential term V (Ψ) in Eq. 3 represents the interaction of the electrons with the reaction field produced by the solvent polarized by the charge distribution of the solute. The argument of V (Ψ) indicates that it depends on the electronic state of the solute, introducing a nonlinear character into the Hamiltonian. Introducing Eqs. 4 and 5 into the electronic Hamiltonian Eq. 3 and rearranging, we obtain

$$\begin{array}{ll}{H}_{|\mathsf{\Psi }>}=& \frac{1}{2}{\displaystyle \sum _i{p}_i^2+V_{en}+V_{ee}+V_{nn}+V_{\left(\mathsf{\Psi }\right)}+{\displaystyle \sum _{KL}{\mathbf{M}}_K^T{\mathbf{h}}_{KL}^{DSO}{\mathbf{M}}_L}}\\ & {\displaystyle \sum _K{\mathbf{M}}_K^T{\mathbf{h}}_K^{PSO}}+{\displaystyle \sum _K{\mathbf{M}}_K^T({\mathbf{h}}_K^{FC}+{\mathbf{h}}_K^{SD})}\end{array}$$

(6)

where the diamagnetic spin-orbit (DSO) and paramagnetic spin-orbit (PSO) operators are given by

$${\mathbf{h}}_{KL}^{DSO}=\frac{{\alpha }^4}{2}{\displaystyle \sum _i\frac{({\mathbf{r}}_{iK}^T{\mathbf{r}}_{iL}){\mathrm{I}}_3-{\mathbf{r}}_{iK}{\mathbf{r}}_{iL}^T}{r_{iK}^3{r}_{iL}^3}}$$

(7)

$${\mathbf{h}}_K^{PSO}={\alpha }^2{\displaystyle \sum _i\frac{{\mathbf{r}}_{iK}\times {\mathbf{p}}_i}{{\mathbf{r}}_{iK}^3}}$$

(8)

and the triplet Fermi-contact (FC) and spin-dipole (SD) operators by

$${\mathbf{h}}_K^{FC}=\frac{8\pi {\alpha }^2}{3}{\displaystyle \sum _i\delta ({\mathbf{r}}_{iK}){\mathbf{s}}_i}$$

(9)

$${\mathbf{h}}_K^{SD}={\alpha }^2{\displaystyle \sum _i\frac{3({\mathbf{s}}_i^T{\mathbf{r}}_{iK}){\mathbf{r}}_{iK}-r_{iK}^2{\mathbf{s}}_i}{r_{iK}^5}}$$

(10)

These operators represent the different mechanisms of the coupling between the nuclear spins mediated by the surrounding electrons. Due to the nonlinearity of the molecular Hamiltonian, the fundamental energetic quantity of a solvated molecule is the so called free-energy functional $𝒢$ defined as

$$\{𝒢\}=\frac{<\mathsf{\Psi }|H_{|\mathsf{\Psi }>}-\frac{1}{2}V(\mathsf{\Psi })|\mathsf{\Psi }>}{<\mathsf{\Psi }\left|\right|\mathsf{\Psi }>}$$

(11)

whose stationarity points correspond to the eigenfunctions of the solute Hamiltonian

$$H_{|\mathsf{\Psi }>}|\mathsf{\Psi }>=E|\mathsf{\Psi }>$$

(12)

Variational approximate solutions of the nonlinear Schrödinger equation in Eq. 12 can be obtained by searching for the stationary points of the free-energy functional with respect to the whole set of variational parameters [33]. The evaluation of the spin–spin coupling constants as second derivatives of the free-energy functional $𝒢$ requires a procedure based on the linear response theory for molecular solutes, i.e. a variational-perturbative scheme which ensures the stationarity condition of the free-energy functional at first order with respect any value of the nuclear magnetic moments. In the following we describe the procedure for the variational Density Functional Theory approach. Ignoring the details of the parameterization, we denote the Kohn-Sham free energy as $𝒢$(M_K,λ_S, λ_T) where λ_S and λ_T are two sets of variational parameters associated with singlet and triplet variations of the electronic state [16]. For the optimized free-energy state, the reduced spin–spin coupling constants may then be calculated as

$$\begin{array}{ll}{\mathbf{K}}_{KL}& =\frac{d^2𝒢}{d{\mathbf{M}}_{K}d{\mathbf{M}}_L}\\ & =\frac{{\partial }^2𝒢}{\partial {\mathbf{M}}_K\partial {\mathbf{M}}_L}+\frac{{\partial }^2𝒢}{\partial {\mathbf{M}}_K\partial {\lambda }_S}\frac{\partial {\lambda }_S}{\partial {\mathbf{M}}_L}+\frac{{\partial }^2𝒢}{\partial {\mathbf{M}}_K\partial {\lambda }_T}\frac{\partial {\lambda }_T}{\partial {\mathbf{M}}_L}\end{array}$$

(13)

The derivatives of the variational parameters λ_S and λ_T with respect M_L are obtained by solving the response equations [11,35]

$$\frac{{\partial }^2𝒢}{\partial {\lambda }_S\partial {\lambda }_S}\frac{\partial {\lambda }_S}{\partial {\mathbf{M}}_L}=-\frac{{\partial }^2𝒢}{\partial {\lambda }_S\partial {\mathbf{M}}_L}$$

(14)

$$\frac{{\partial }^2𝒢}{\partial {\lambda }_T\partial {\lambda }_T}\frac{\partial {\lambda }_S}{\partial {\mathbf{M}}_L}=-\frac{{\partial }^2𝒢}{\partial {\lambda }_T\partial {\mathbf{M}}_L}$$

(15)

where the second derivatives on the left-hand sides are, respectively, the singlet and triplet electronic Hessian. The solution of the singlet linear response equations Eq. 14 gives the perturbing effect of the imaginary singlet paramagnetic spin-orbit (PSO) operator, while the solution of the triplet linear response equations Eq. 15 represents the perturbing effect of the real triplet Fermi-contact (FC) and of the spin-dipole (SD) operator. The perturbing effect of the second-order real singlet diamagnetic spin-orbit (DSO) operators enter the reduced coupling constant via the first term of Eq.2 and do not require solution of any response equation [36]. In many cases, the Fermi-contact FC contribution dominates the spin–spin couplings, but it is impossible to predict with any certainty when the non-FC contributions may be neglected [4], and all terms should therefore in general be included in the calculation.

3 Computational details

The calculation of accurate indirect spin–spin coupling constants is a challenging task, putting severe demands on the quality of the one-electron basis and the treatment of electron correlation effects [4]. Since we will here mainly be concerned with the changes in the spin–spin coupling constants of benzene induced by different solvents, we do not need to have a basis set or a correlation description that gives highly accurate results, although it is important that the method and basis sets used is able to give a qualitatively correct description of the electronic structure of the molecule in the region close to the nuclei. DFT has in recent years been demonstrated to be a very computationally efficient way of obtaining accurate indirect spin–spin coupling constants [16,17], and this is the approach we will use here. Based on a recent implementation of indirect spin–spin coupling constants at the DFT level [16], we have included the solvent contributions using the PCM model to the singlet and linear response functions as described in Section 2. A study by Helgaker et al. [5] demonstrated that basis sets based on Huzinaga’s compilation [37] with polarization functions by van Wüllen [38], decontracted in the s space and with additional tight s functions, provide a good compromise between the accuracy of the results and the size of the basis set. We will in this work employ the Huz-IIsu2 basis set [5] where two tight s functions have been added to the atomic Huz-II basis set in a geometric progression. There are in general two distinct contributions to an observed solvent effect: (1) The direct change in the electronic wave function due to the polarization of the molecular wave function by the solvent, and (2) An indirect effect caused by changes in the molecular geometry as the molecule is placed in the solvent. To account for the latter effect, we have optimized the geometry of the molecules in the different solvents using the gradient formalism recently presented [39]. All calculations have used D_2h symmetry following the algorithm described in Ref. [40]. It is important to emphasize that our algorithm for generating the molecular cavity will not generate a cavity possessing the full D_6h symmetry, but rather a cavity of D_2h symmetry, the highest point group used in the Dalton program [41]. For this reason we obtain a slight symmetry-breaking of the solvent-optimized geometry. Although the symmetry-breaking effects on the molecular geometry is very small (vide infra), they do give rise to noticeable differences in the calculated indirect spin–spin coupling constants, and we return to this point in the next section. The solvent cavity was built from six interlocking spheres centered on the carbon atoms of the benzene molecule, each with a radius of 1.90 Å. All calculations have been done with a local version of the Dalton program [41].

4 Results

We have summarized our results for the benzene molecule in gas phase in Table 1. Our results have been calculated at the optimized geometry of the molecule using the B3LYP functional and the Huz-IIsu2 basis set. The results are compared to recently derived “experimental” equilibrium values for the indirect spin–spin coupling constants [23]. These values are meant to be more directly comparable to theoretical results, as they have been derived by subtracting a theoretically estimated zero-point vibrational contribution to the indirect spin–spin coupling constants from the experimentally observed couplings. As such, they are meant to represent the experimental spin–spin coupling constant of a non-vibrating benzene molecule, which is what we obtain in our calculation. In Table 1 we have also collected the results of other recent theoretical investigations. We note that our results in general are in quite good agreement with the vibrationally corrected experimental results, although a difference of 15 Hz (about 10% ) exists in the case of ¹J(H¹³C). Only minor differences between our results and the DFT/B3LYP results of Ref. [23] obtained with the HuzIII-su3 basis set are observed. Our results corroborate the observation [23,42] that the MCSCF results significantly overestimate most of the coupling constants. In a recent paper, an attempt was made to extract a vibrationally averaged ¹J(H¹³C) coupling constant for benzene in the gas phase by extrapolation of accurately measured coupling constants in four relatively non-polar solvents to a dielectric medium with a dielectric constant of 1 (which corresponds to the dielectric constant of vacuum) [43]. In addition to depend strongly on the solvent being well represented as a structureless dielectric continuum, such an approach also relies heavily on the assumption that the spin–spin coupling constant varies smoothly with the dielectric constant. Although this may be the case for the direct solvent effects, it is much less likely that this will be the case for the changes induced by the changes in molecular geometry. We have reinvestigated this extrapolation scheme by PCM calculations using the same solvents as in Ref. [43]. Our results for ¹J(H¹³C) are collected in Table 2. In order to make the most consistent comparison with the available experimental data of Ref. [44], we refer all solvent shifts to the results obtained in fully deuterated cyclohexane. We note from Table 2 that the agreement between our results and experiment is fair. There are two important observations that can be made from the results in Table 2: (1) The indirect contributions to the solvent shifts arising from the change in the molecular geometry are significant, being about 20-30% of the total calculated solvent shift. (2) Our calculations slightly underestimate the solvent shift in CS₂ relative to cyclohexane, whereas we overestimate the effects for pyridine and acetone. However, this is consistent with our model, in which we consider an infinitely diluted benzene molecule in a pure solvent, whereas a more reasonable representation would be a 5% solution of benzene in solvent [43]. This could be modeled using an “effective” dielectric constant [43,45], which would lead to an increase of the dielectric constant of the cyclohexane solution and reduce the effective dielectric constants of the other solutions, thus probably leading to a better agreement with the experimental observations. However, we will not explore this approach here. A final observation that can be made from the data in Table 2 is that there is a large discrepancy between theory and experiment in the predicted gas-to-cyclohexane solvent shift, the experimental shift being almost twice as large as the theoretical prediction. If we try to apply the linear extrapolation done in Ref. [43] to our theoretical data, we get a predicted gas-to-cyclohexane solvent shift of 0.869 Hz and 1.112 Hz for the coupling constants calculated at the vacuum geometry and the solvent-relaxed geometry, respectively. Both of these numbers are almost twice as large as the calculated solvent shift, clearly illustrating the inadequacy of the linear regression for obtaining the gas-phase value of the spin–spin coupling constant. Assuming that the experimentally extrapolated gas-phase value has used a cyclohexane-to-gas shift that is too large by approximately a factor of 2, we therefore propose a revised gas-phase value of ¹J(H¹³C) for a vibrating benzene molecule at 300 K of 157.46 Hz, to be compared with the original proposal of 156.99 Hz. Subtracting from this the recently calculated zero-point vibrational corrections of 4.8 Hz [23], we obtain a derived experimental gas-phase value for the equilibrium geometry spin–spin coupling constant of 152.7 Hz. This is likely to be an upper limit, as it can be expected that the zero-point vibrational correction underestimates the true vibrational corrections that would be obtained at the experimental temperature of 300 K. Let us also address the issue of symmetry breaking of our molecular cavity. As can be seen from Table 2, we report two slightly different numbers for the solvent-relaxed geometries. This symmetry-breaking arises because we in our construction of the cavity only enforce symmetries of D_2h and subgroups thereof [40]. We will therefore not be able to enforce the full D_6h symmetry of the benzene molecule in the tessellation of the cavity. As such we will get two classes of carbon and hydrogen atoms, one containing two symmetry-related carbon/hydrogen atoms, and one containing four symmetry-related carbon/hydrogen atoms. Although the differences in the bond lengths and angles are for all practical purposes negligible (CC bond lengths differing by less than 0.05 pm, CH bond lengths by less than 0.04 pm, and bond angles by about 0.04°), the sensitivity of the spin–spin coupling constants to geometry changes makes these small differences quite visible in our calculated coupling constants. Obviously, we do not experience any symmetry-breaking in the gas-phase optimized geometry, but since the solvent shifts are given relative to the value in cyclohexane, we obtain two different values for the gas-phase shifts since we have two different ¹J(H¹³C) coupling constants in cyclohexane. Similar differences can be expected for the other coupling constants, but for these couplings we have only solved response equations for one symmetry-class of atoms [46], and we therefore do not observe these differences. It is generally known that solvent effects are of minor importance for couplings other than one-bond couplings [25], and this is largely confirmed for the solvent effects calculated for the spin–spin coupling constants ⁿJ(H¹³C), ⁿJ(¹³C¹³C), and ⁿJ(HH), collected in Table 3, Table 4 and Table 5, respectively. Indeed, sizeable solvent effects are only observed for ¹J(H¹³C) and ¹J(¹³C¹³C), where the calculated solvent shifts going from gas-phase to an acetone solution are about 1.5 Hz and -1.3 Hz, respectively. Interestingly, the ¹J(¹³C¹³C) couplings display only a negligible geometry effect, the solvent-induced changes for this coupling arising only from polarization of the electronic structure of the molecule by the surrounding dielectric medium. The two-bond carbon-hydrogen coupling constant ²J(H¹³C) also has a non-negligible solvent effect of about 0.2 Hz going from gas-phase to acetone, with the geometry changes accounting for about 20% of the total shift. Indeed, in relative terms this is the largest solvent effect of all the calculated coupling constants, as the solvent effect is more than 10% of the gas-phase value. In comparison, the solvent shift for ¹J(¹³C¹³C) is only about 2.5% of the gas-phase value. The effects of the dielectric medium and the solvent-induced geometry changes are negligible for all other coupling constants, in no instances exceeding 0.1 Hz even for the most polar solvent. The two largest coupling constants, ¹J(¹³C¹³C) and ¹J(¹³CH), are both dominated by the Fermi contact interaction. This predominance of the Fermi-contact interaction is also reflected in the contributions to the solvent shifts observed for these coupling constants. Even when going from gas-phase to the solvent-optimized geometry of benzene in acetone, the contribution from the other coupling mechanisms are -0.068 and -0.040 Hz for ¹J(¹³CH) and ¹J(¹³C¹³C), respectively. Even for the one-bond carbon-carbon coupling constant, where the PSO contribution is about 10% of the Fermi-contact contribution, but with opposite sign, the solvent effect on the PSO contribution is a negligible 0.007 Hz. However, more studies on coupling constants with a smaller contribution from the Fermi-contact mechanism is needed in order to establish whether the dielectric medium effects in general only affect the Fermi contact contribution. In comparison with experiment, our results for the ⁿJ(¹³CH) (n = 2, 3, 4) and ⁿJ(HH) (n=3,4,5) couplings are not very satisfactory. For ²J(¹³CH) the agreement is quite good when considering the solvent shifts of pyridine and acetone relative to CS₂. However, whereas we obtain a rather sizeable solvent shift of 0.103 Hz from cyclohexane to carbon disulfide, no noticeable difference in the coupling constants is observed in experiment in these two solvents. It is interesting to note that the experimental solvent shifts for ³J(¹³CH), ⁴J(¹³CH), ⁴J(HH) relative to cyclohexane pass through a maximum for pyridine. Only in the case of ⁴J(¹³CH) are we able to reproduce this trend, whereas a similar maximum also occurs for ³J(HH) in our calculations. Our results indicate that this maximum in the solvent shifts arise from the geometry relaxation, compensating and opposing the effects of the direct polarization of the electronic structure of the solvated molecules. However, our results does not appear reliable enough to enable us to draw a definite conclusion on this point.

Table 1. Calculated gas-phase values for benzene at the optimized geometry. Comparison is made with available literature data. All numbers reported in Hz.

Coupling	DFT/B3LYP	Exp. Jeqa	MCSCF	DFT/B3LYPb
	This work	Ref. [23]	Ref. [42]	Ref. [23]
1J(H13C)	164.636	153.8	176.7	166.3
2J(H13C)	1.798	1.4	-7.4	2.0
3J(H13C)	7.877	7.0	11.7	8.0
4J(H13C)	-1.329	-1.0	-1.3	-1.2
3J(HH)	8.571	7.0		8.7
4J(HH)	0.962	1.2		1.3
5J(HH)	0.661	0.6		0.8
1J(13C13C)	59.244	56.1	70.9	60.0
2J(13C13C)	-1.768	-1.7	-5.0	-1.8
3J(13C13C)	10.822	9.4	19.1	11.2

^aDerived experimental equilibrium value. Obtained by taking the experimental observations from Ref. [42] and subtracting the zero-point vibrational corrections as calculated in Ref. [23].

^bCalculated at the optimized geometry using the Huz-IIIsu3 basis set and the B3LYP functional.

Table 1. Calculated gas-phase values for benzene at the optimized geometry. Comparison is made with available literature data. All numbers reported in Hz.

Coupling	DFT/B3LYP	Exp. Jeqa	MCSCF	DFT/B3LYPb
	This work	Ref. [23]	Ref. [42]	Ref. [23]
1J(H13C)	164.636	153.8	176.7	166.3
2J(H13C)	1.798	1.4	-7.4	2.0
3J(H13C)	7.877	7.0	11.7	8.0
4J(H13C)	-1.329	-1.0	-1.3	-1.2
3J(HH)	8.571	7.0		8.7
4J(HH)	0.962	1.2		1.3
5J(HH)	0.661	0.6		0.8
1J(13C13C)	59.244	56.1	70.9	60.0
2J(13C13C)	-1.768	-1.7	-5.0	-1.8
3J(13C13C)	10.822	9.4	19.1	11.2

^aDerived experimental equilibrium value. Obtained by taking the experimental observations from Ref. [42] and subtracting the zero-point vibrational corrections as calculated in Ref. [23].

^bCalculated at the optimized geometry using the Huz-IIIsu3 basis set and the B3LYP functional.

Table 2. Calculated and experimental solvent shifts of ¹J(H¹³C) for benzene in four different solvents. Both relaxed and unrelaxed geometries reported. All numbers reported in Hz.

Solvent	Unrelaxed geometry	Relaxed geometry	Experiment (Ref. [44])
Gas	-0.426	-0.583/-0.534	-0.927
C6D12	0.000	0.000	0.000
CS2	0.135	0.182/0.165	0.215
Pyridine-d5	0.661	0.896/0.836	0.725
Acetone-d6	0.734	1.043/0.938	0.818

Table 2. Calculated and experimental solvent shifts of ¹J(H¹³C) for benzene in four different solvents. Both relaxed and unrelaxed geometries reported. All numbers reported in Hz.

Solvent	Unrelaxed geometry	Relaxed geometry	Experiment (Ref. [44])
Gas	-0.426	-0.583/-0.534	-0.927
C6D12	0.000	0.000	0.000
CS2	0.135	0.182/0.165	0.215
Pyridine-d5	0.661	0.896/0.836	0.725
Acetone-d6	0.734	1.043/0.938	0.818

Table 3. Calculated solvent shifts of ⁿJ(H¹³C) for benzene in four different solvents. Both relaxed and unrelaxed geometries reported. Solvent shifts reported relative to the cyclohexane value. Experimental data (in italics) from Ref. [44]. All numbers reported in Hz.

Solvent	2J(H13C)		3J(H13C)		4J(H13C)
	Unrelaxed	Relaxed	Unrelaxed	Relaxed	Unrelaxed	Relaxed
Gas-phase	-0.072	-0.089	0.019	0.020	-0.010	-0.004
CS2	0.084	0.103	-0.022	-0.023	0.011	0.005
	0.001		-0.016		0.011
Pyridine-d5	0.156	0.180	-0.042	-0.043	0.021	0.011
	0.077		-0.039		0.015
Acetone-d6	0.166	0.204	-0.044	-0.048	0.022	0.010
	0.095		-0.031		0.007

Table 3. Calculated solvent shifts of ⁿJ(H¹³C) for benzene in four different solvents. Both relaxed and unrelaxed geometries reported. Solvent shifts reported relative to the cyclohexane value. Experimental data (in italics) from Ref. [44]. All numbers reported in Hz.

Solvent	2J(H13C)		3J(H13C)		4J(H13C)
	Unrelaxed	Relaxed	Unrelaxed	Relaxed	Unrelaxed	Relaxed
Gas-phase	-0.072	-0.089	0.019	0.020	-0.010	-0.004
CS2	0.084	0.103	-0.022	-0.023	0.011	0.005
	0.001		-0.016		0.011
Pyridine-d5	0.156	0.180	-0.042	-0.043	0.021	0.011
	0.077		-0.039		0.015
Acetone-d6	0.166	0.204	-0.044	-0.048	0.022	0.010
	0.095		-0.031		0.007

Table 4. Calculated solvent shifts of ¹J(¹³C¹³C) for benzene in four different solvents. Both relaxed and unrelaxed geometries reported. Solvent shifts reported relative to the gas-phase values. All numbers reported in Hz.

Solvent	1J(13C13C)		2J(13C13C)		3J(13C13C)
	Unrelaxed	Relaxed	Unrelaxed	Relaxed	Unrelaxed	Relaxed
C6D12	-0.550	-0.565	0.000	-0.015	-0.029	-0.037
CS2	-0.638	-0.653	-0.001	-0.017	-0.034	-0.042
Pyridine-d5	-1.223	-1.266	0.002	-0.025	-0.067	-0.082
Acetone-d6	-1.304	-1.346	0.002	-0.031	-0.071	-0.090

Table 4. Calculated solvent shifts of ¹J(¹³C¹³C) for benzene in four different solvents. Both relaxed and unrelaxed geometries reported. Solvent shifts reported relative to the gas-phase values. All numbers reported in Hz.

Solvent	1J(13C13C)		2J(13C13C)		3J(13C13C)
	Unrelaxed	Relaxed	Unrelaxed	Relaxed	Unrelaxed	Relaxed
C6D12	-0.550	-0.565	0.000	-0.015	-0.029	-0.037
CS2	-0.638	-0.653	-0.001	-0.017	-0.034	-0.042
Pyridine-d5	-1.223	-1.266	0.002	-0.025	-0.067	-0.082
Acetone-d6	-1.304	-1.346	0.002	-0.031	-0.071	-0.090

Table 5. Calculated and experimental solvent shifts of ⁿJ(HH) for benzene in four different solvents. Both relaxed and unrelaxed geometries reported. Solvent shifts reported relative to the cyclohexane value. Experimental data (in italics) from Ref. [44]. All numbers reported in Hz.

Solvent	3J(HH)		4J(HH)		5J(HH)
	Unrelaxed	Relaxed	Unrelaxed	Relaxed	Unrelaxed	Relaxed
Gas-phase	-0.022	-0.002	-0.004	-0.006	-0.002	-0.006
CS2	0.026	0.003	0.004	0.007	0.003	0.006
	-0.022		0.004		-0.009
Pyridine-d5	0.049	0.008	0.008	0.014	0.005	0.011
	-0.021		0.013		-0.003
Acetone-d6	0.052	0.002	0.008	0.015	0.005	0.013
	-0.009		0.006		0.005

Table 5. Calculated and experimental solvent shifts of ⁿJ(HH) for benzene in four different solvents. Both relaxed and unrelaxed geometries reported. Solvent shifts reported relative to the cyclohexane value. Experimental data (in italics) from Ref. [44]. All numbers reported in Hz.

Solvent	3J(HH)		4J(HH)		5J(HH)
	Unrelaxed	Relaxed	Unrelaxed	Relaxed	Unrelaxed	Relaxed
Gas-phase	-0.022	-0.002	-0.004	-0.006	-0.002	-0.006
CS2	0.026	0.003	0.004	0.007	0.003	0.006
	-0.022		0.004		-0.009
Pyridine-d5	0.049	0.008	0.008	0.014	0.005	0.011
	-0.021		0.013		-0.003
Acetone-d6	0.052	0.002	0.008	0.015	0.005	0.013
	-0.009		0.006		0.005

5 Summary

In this paper we have presented the theory and implementation of the polarizable continuum model for singlet and triplet linear response functions in the Dalton quantum chemistry program [41] using density functional theory. The implementation uses D_2h symmetry in all parts of the calculation, including the generation of the PCM cavity. This fully analytical implementation of the PCM model for the calculation of solvent effects on indirect spin–spin coupling constants have been applied to the study of the solvent effects on the coupling constants in benzene. For the largest coupling constants—¹J(¹³CH) and ¹J(¹³C¹³C)—our results for the solvent shifts are in quite good agreement with experiment when the molecular geometry is allowed to relax in the solvent. For the other coupling constants, the solvent effects are too small for other predicted solvent shifts to be considered accurate, and for these coupling constants, agreement with experiment is only fair. Our results do not support the use of extrapolation schemes to estimate coupling constants of gas-phase molecules from the solvent shifts observed in solution. Our calculations have demonstrated that such extrapolation schemes significantly overestimate the cyclohexane-to-gas shift of ¹J(¹³CH), and we propose instead a new “experimental” gas-phase value of the ¹J(¹³CH) spin–spin coupling constant of a vibrating molecule at 300 K to be 157.46 Hz.