Molecular Mechanics Approach for

ARTICLE pubs.acs.org/JPCB

Density Functional Theory/Molecular Mechanics Approach for Electronic g-Tensors of Solvated Molecules Zilvinas Rinkevicius,*,†,‡ N. Arul Murugan,† Jacob Kongsted,§ Ke)stutis Aidas,† Arnfinn Hykkerud Steindal,|| and Hans Ågren† †

)

Department of Theoretical Chemistry and Biology, School of Biotechnology, Royal Institute of Technology, SE-106 91 Stockholm, Sweden ‡ Swedish e-Science Research Center (SeRC), Royal Institute of Technology, 10044 Stockholm, Sweden § Department of Physics and Chemistry, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark Centre of Theoretical and Computational Chemistry, Department of Chemistry, University of Tromsø, N-9037 Tromsø, Norway ABSTRACT: A general density functional theory/molecular mechanics approach for computation of electronic g-tensors of solvated molecules is presented. We apply the theory to the commonly studied di-tert-butyl nitroxide molecule, the simplest model compound for nitroxide spin labels, and explore the role of an aqueous environment and of various approximations for its treatment. It is found that successive improvements of the solvent shift of the g-tensor are obtained by going from the polarizable continuum model to discrete solvent models of various levels of sophistication. The study shows that an accurate parametrization of the electrostatic potential and polarizability of the solvent molecules in terms of distributed multipole expansions and anisotropic polarizabilities to a large degree relieves the need to explicitly include water molecules in the quantum region, which is the common case in density functional/continuum model approaches. It is also shown that the local dynamics of the solvent around the solute significantly influences the electronic g-tensor and should be included in benchmarking of exchange-correlation functionals for evaluation of solvent shifts of g-tensors. These findings can have important ramifications for the use of advanced hybrid density functional theory/molecular mechanics approaches for modeling spin labels in solvents, proteins, and membrane environments.

I. INTRODUCTION Modeling of electron paramagnetic resonance (EPR) spin Hamiltonian parameters, that is, electronic g-tensors, hyperfine couplings, and zero-field splitting constants, of radicals embedded in liquids or solid phase environments has become a popular subject of contemporary electronic structure theory.1
3 In such studies the effect of the solvent has been accounted for at different levels of sophistication—starting from simplistic Onsager type continuum models (see, for example, ref 4) to complex “dynamical” discrete/continuum models, like GLOB/ADMP,5 in which part of the solvent close to the radical is described as a dynamic discrete environment and the remaining part is treated as a continuum. Among a variety of available electronic structure methods, density functional theory (DFT) combined with the polarizable continuum model (PCM) appears to be the most popular choice owing to the inherent qualitative accuracy and low computational cost.2,3,6,7 The capability of the DFT/PCM methodology has been demonstrated in recent years, including, but not limited to, studies of EPR spin Hamiltonian parameters of amino acid radicals, various quinone radical anions, and nitroxides in aprotic and protic solvents.2,3,6,7 However, the applications of the r 2011 American Chemical Society

DFT/PCM approach have been less promising in more complex environments, like mixtures of solvents, liquid crystals, protein surfaces, or cellular membranes, as typically these require an explicit inclusion of a large portion of the nonhomogeneous environment around the radical, something that dramatically increases the computational cost. Due to these circumstances, only a few DFT/PCM studies of EPR parameters of radicals or other paramagnetic species in proteins have appeared.8
10 Several ways have been proposed to overcome the limitations of PCM and similar continuum models for radicals embedded in liquids or solid phase environments, like cluster approach,11,12 electron density embedding formalism,13 and hybrid quantum mechanics/molecular mechanics (QM/MM) methods.14
18 Among quantum chemistry methods that treat the interaction between the radical and the environment beyond the continuum model, the QM/MM methods have now received considerable popularity, and several works have recently appeared devoted to Received: November 14, 2010 Revised: March 8, 2011 Published: March 30, 2011 4350

dx.doi.org/10.1021/jp1108653 | J. Phys. Chem. B 2011, 115, 4350–4358

The Journal of Physical Chemistry B QM/MM modeling of electronic g-tensors and hyperfine coupling constants of radicals and other paramagnetic species in solution and protein environments.15
20 Studies, like that of the Cu(II) ion in blue copper proteins,17 have demonstrated that the QM/MM approach is of superior accuracy for electronic g-tensors compared to QM methods combined with the continuum solvent model. However, QM/MM methods suitable for computation of EPR spin Hamiltonian parameters have so far been limited to rudimentary approaches in which the interaction between the QM and MM regions is described by the QM density interacting with point charges in the MM region14,15,18 or by an approximate electrostatic interaction between fitted multipoles, determined from the QM electron density with “point charges þ higher order multipoles” in the MM region.16 Therefore, until now, to the best of our knowledge, more advanced QM/MM methods in which the interaction between the QM and MM regions in addition to the electrostatic terms also include polarization/induction as well as other important terms, like charge transfer and exchange-repulsion interactions, have not been applied in investigations of EPR spin Hamiltonian parameters. As such inclusions have turned out to be important for a variety of optical and magnetic properties in closed-shell systems,21,22 and in order to fill this gap, we here systematically investigate the applicability of the DFT/MM approach with a more advanced treatment of the coupling between the QM and MM regions. We apply this methodology to predict the electronic g-tensor of di-tert-butyl nitroxide (DTBNO) in water and explore the level of sophistication required in the description of the MM region in order to achieve an accurate treatment of the environment effects.

II. DFT/MM METHODS: THE LADDER OF SOLVENT DESCRIPTIONS Various DFT/MM methods can be distinguished at the fundamental level one from the other by the way in which the coupling between the QM and MM regions is treated. As noted above, the mainstream implementations of DFT/MM methods traditionally limit this coupling to the terms involving the electron density of the QM region interacting with point charges and higher multipoles in the MM region.14
16,18 DFT/MM methods, which go beyond this approximation, have been implemented only recently.23
27 Among such methods, the most sophisticated description of the interaction between the MM and QM regions is achieved in the so-called effective fragment potential (EFP) method in which the coupling between the QM and MM regions includes electrostatic, polarization/induction, exchange repulsion, and charge transfer interaction terms. Unfortunately, the DFT/EFP method cannot in its current implementation be used for calculation of EPR spin Hamiltonian parameters. Here we will focus on the more general DFT/MM response formalism by Olsen et al.,25 termed polarizable embedding, which is capable of computing arbitrary molecular properties expressed via linear, quadratic, and cubic response functions and their residues. Since this DFT/MM response formalism includes electrostatic and polarization/induction terms in the description of the interaction between a solute and its environment, we can establish a ladder of sophistication for the coupling between the QM and MM regions, referring to truncating the MM molecule or fragment electrostatic potential at different orders of the multipole expansion and by selecting different forms of the polarizability tensor for the MM molecule or the fragment. For the purpose of determining the

ARTICLE

impact of the level of detail at which molecules or fragments in the MM region are described, as well as for the sake of establishing the limits of applicability of the traditional treatment of the QM to MM coupling in which the MM region is only represented by set of point charges, we consider five different models of the water molecule. As shown in Figure 1, for the lowest level approximation of coupling between the QM and MM regions, denoted as a MM-0, we adopt the traditional QM/MM approach and describe the MM water molecule as a collection of point charges located on the oxygen and hydrogen atoms, where the point charges are taken from the TIP3P force field,28 which implicitly accounts for media polarization effects. The next natural step in improving the description of the water molecules in the MM region is the introduction of polarization/induction interaction between the MM and QM regions by assigning an isotropic polarizability of the water molecule to the oxygen atom and point charges to oxygen and hydrogen atoms according to the Ahlstr€om force field,29 denoted here as MM-1. Further improvement of the MM water description, denoted here as the MM-2 model, is achieved by describing the water polarizability by three polarizability tensors distributed over the oxygen and hydrogen atoms (see Figure 1). These water models can be considered lower level approximations, as the electrostatic potential in all of them is described at the level of only point charges. The next two models, MM-3 and MM-4, overcome this deficiency: (a) in the MM-3 model the electrostatic potential of water is described using point charges placed on the hydrogen and oxygen atoms and dipoles and quadrupoles placed on the hydrogen and oxygen atoms as well as on the midpoints of the two OH bonds; (b) the MM-4 model further improves the description of the water electrostatic potential by adding octopoles placed on the hydrogen and oxygen atoms and on the midpoints of the two OH bonds. Apart from enhancing the description of the electrostatic potential of the water molecules, the MM-3 and MM-4 models also improve the description of the water polarizability over the MM-2 model by introducing two additional expansion points (OH bond midpoints) on which a distributed polarizability tensor is placed. Thus, the MM-3 and MM-4 models describe the polarizability of the water molecule at the same level of accuracy; they slightly differ in their description of the electrostatic potential of water; i.e., MM-4 goes one order higher in multipole expansion compared to MM-3. This ladder of MM water models (see Figure 1) covers various regimes of coupling between the MM and QM regions and thus introduces a successively more sophisticated description of the electrostatic and polarization/ induction interactions between the MM and QM regions of the systems. Therefore, these MM water models constitute a general benchmarking suite for the DFT/MM response formalism, which allows to estimate the sophistication level needed in the description of a structured environment of the radical for an accurate estimation of the solvent shifts of its electronic g-tensor.

III. COMPUTATIONAL DETAILS In this work evaluation of the electronic g-tensor of DTBNO in aqueous solution at ambient temperature has been carried out using the sequential molecular dynamics/quantum mechanics (MD/QM) scheme:7,18,30,31 (a) MD simulations were performed employing the Car
Parrinello mixed QM/MM approach; (b) electronic g-tensor calculations over the set of snapshots extracted from MD trajectory were performed using the hybrid DFT/MM 4351

dx.doi.org/10.1021/jp1108653 |J. Phys. Chem. B 2011, 115, 4350–4358

The Journal of Physical Chemistry B

ARTICLE

Figure 1. Ladder of force fields for MM water molecules used in DFT/MM calculations: (a) MM-0 denotes the TIP3P force field,28 which is parametrized via three point charges (
q, q/2, q/2); (b) MM-1 denotes the Ahlstr€om force field,29 which is described by three point charges (
q, q/2, q/2) and isotropic polarizability, Riso, located on oxygen atom; (c) MM-2 denotes a modified Ahlstr€om force field29 in which Riso is substituted by three B , generated by the LoProp method63 at the B3LYP/aug-cc-pVTZ level; (d) MM-3 denotes a force field in B localized anisotropic polarizability tensors, R B B , and polarizability of water is which electrostatic potential of water described by a set of distributed point charges q, dipoles μB and quadrupoles Θ 63 described by a set of distributed anisotropic polarizability tensors (all quantities generated by the LoProp method at the B3LYP/aug-cc-pVTZ level); B B , to description of electrostatic potential of water. (e) MM-4 denotes a force which extends the MM-3 force field by addition of distributed octopoles, Ω

response method. The details of computational procedures used in each step of sequential MD/QM procedure are outlined in the following subsections. A. Hybrid Car
Parrinello Molecular Dynamics Simulation. The structure of DTBNO in aqueous solution at ambient temperature has been determined by means of Car
Parrinello mixed QM/MM molecular dynamics simulation, in which a single DTBNO molecule along with 9679 water molecules has been placed in an orthorhombic box with size of approximately, 68.4, 65.8, and 65.4 Å, and simulation has been carried out for a total production time of 21 ps employing the canonical ensemble. From the generated trajectory, 86 snapshots have been extracted at equal intervals of time to be used in DFT/MM calculations of electronic g-tensors. In these simulations, the DTBNO molecule has been treated at the density functional theory level (BLYP functional32,33 and Troullier
Martins norm conserving pseudopotentials34 for the core electrons) and water molecules have been described by the classical TIP3P force field,28 respectively. In this setup of MD simulation, we used following additional parameters: energy cutoff was set to 80 Ry for expansion of electronic wave function in a plane wave basis; the time step for the integration of the equation of motion and the fictitious electronic mass were selected to be 5 au and 600 amu, respectively; temperature control within canonical ensemble has been accomplished by Nose
Hoover thermostat.35,36 The starting equilibrated structure of “DTBNO þ 9679 waters” ensemble used in Car
Parrinello mixed QM/MM molecular dynamics simulation has been obtained from 100 ps length classical molecular dynamics simulation in isothermal
isobaric ensemble. In this simulation, the general amber force field37 along with charges generated by the CHELPG38 procedure at the B3LYP/6-31þG* level of theory has been used for DTBNO and TIP3P force field has been used for water molecules, respectively. Finally, the determination of fitted charges

for DTBNO has been accomplished using the CHELPG38 procedure implementation in the Gaussian0339 software package. All Car
Parrinello mixed QM/MM molecular dynamic simulations have been carried out using the CPMD 3.11/GROMOS package,40 while for classical molecular dynamics simulations we resorted to the Amber 8 package.41 B. Hybrid DFT/MM Computations of Electronic g-Tensors. Calculations of the electronic g-tensor of DTBNO in this work have been performed using the hybrid DFT/MM approach in which the MM region was constructed from MD trajectory snapshot by including all water molecules within 20 Å radius from the center of the QM region and the QM region was restricted to be either DTBNO or DTBNO with two water molecules closest to the oxygen atom of the NO bond. This selection of the MM region size ensures that a converged description with respect to bulk solvent environment is achieved according to earlier studies of various molecular properties of solvated molecules using the DFT/MM approach.21,42 Five force fields of different level of sophistication, denoted as MM-X (X = 0
4), have been employed to describe the MM region in DFT/MM calculations. Detailed descriptions of these force field are given in the previous section devoted to the ladder of solvent environment description within the DFT/MM approach. In analogy to previous studies43,44 of DTBNO in various solvents, we employ the BP8632,45 and B3LYP32,33,46,47 exchange-correlation functionals in our calculations of electronic g-tensors, since these two functional are most widely used in calculations of this type and are known to give accurate values of electronic g-tensor shifts for a wide range of organic and inorganic radicals.48
50 Both functionals have been combined with the Huz-III basis set,51
54 which are frequently employed in studies of various magnetic properties. In these calculations we used a spin-restricted approach50 for evaluation of electronic g-tensor in which 4352

dx.doi.org/10.1021/jp1108653 |J. Phys. Chem. B 2011, 115, 4350–4358

The Journal of Physical Chemistry B

Figure 2. Orientation of g-tensor principle axes in di-tert-butyl nitroxide and its (SOMO) and (SOMO
1) orbitals, which defines the electronic g-tensor of di-tert-butyl nitroxide in vacuum and solution.

mass-velocity, one-electron gauge, one- and two-electron spin
orbit contributions to the g-tensor shift have been accounted for and a two-electron spin
orbit operator has been substituted by the approximate mean-field spin
orbit operator.55 The above-described methodology has been used in all calculations of the electronic g-tensor of DTBNO using the DFT/MM method implemented in a development version of the DALTON program.56

IV. SOLVENT SHIFTS OF THE ELECTRONIC G-TENSORS: DI-TERT-BUTYL NITROXIDE IN WATER A. Rationalization by Stone’s Theory. Di-tert-butyl nitroxide is a typical persistent radical from the family of nitroxides, and its EPR spin-Hamiltonian parameters, namely, the electronic g-tensor and the hyperfine coupling constant of the nitrogen atom, have been extensively studied both experimentally and theoretically.19,43,44,57,58 The behavior of its electronic g-tensor is well understood and can be rationalized by Stone’s theory59,60 according to which the electronic g-tensor of DTBNO and other nitroxides is dominated by a single contribution arising from the spin
orbit coupling between the singly occupied molecular orbital (SOMO) and the lone pair molecular orbital (SOMO
1) lying right below; see Figure 2. Thus, the electronic g-tensor of DTNBO follows a distinct pattern for the components gzz > gxx > gyy ≈ ge, where the z principal axis is oriented along the NO bond and the x principal axis is lying in the plane of the R2NO group and the y principal axis is perpendicular to the plane of the R2NO group. This picture for the electronic g-tensor of DTBNO is supported by recent DFT calculations,43,44,57 and consequently the solvent shifts of the g-tensor observed in aprotic and protic solvents can be interpreted by this rather simplistic model taking into account only the solvent influence on the SOMO and SOMO
1 orbitals. According to this picture the solvent shift is largest for the gzz component and less pronounced for the gxx component, whereas the smallest gyy component is almost not influenced at all by the solvent. The g-tensor shift is in this model caused by three mechanisms: (a) redistribution by the solvent of spin density in the SOMO orbital located predominantly on the NO bond; (b) blue shift of the (SOMO
1) f (SOMO) excitation of (n) f (π*) type, which increases the denominator of the leading perturbation theory term in Stone’s theory and therefore decreases the g-tensor shift of the gzz component; (c) delocalization of the lone pair orbital (SOMO
1) by the solvent. However, according to our previous study of the electronic g-tensor of DTBNO,43 the second mechanism, namely, the blue shift of the (SOMO
1) f (SOMO) transition, is the dominant mechanism responsible for the solvent shift, while the first and third mechanisms are of less importance. This type of rationalization has been successfully applied to interpret EPR measurements in both aprotic and

ARTICLE

protic solvents. Therefore, the well-established mechanism of solvent shifts of the g-tensor component and the availability of previous DFT/PCM and DFT/COSMO results make DTBNO an ideal candidate for testing the performance of DFT/MM methods. B. DFT/MM Results for Di-tert-butyl Nitroxide. In aqueous solution DTBNO and other similar nitroxides typically form hydrogen bonds between the NO oxygen atom and one or two water molecules.13,57,61,62 Based on this, we designed two models for the QM region: one in which only the DTBNO molecule is included in the QM region and the other in which the DTBNO radical and two water molecules hydrogen bonded to NO in DTBNO are included in the QM region. Both models have been employed in previous studies of DTNBO as well as other nitroxides using the DFT/PCM and DFT/COSMO methods.43,44,57 Our results for DTBNO obtained using the DFT/PCM and DFT/MM-X (X = 0
4) methods are tabulated in Table 1, where values of the electronic g-tensor shift components are averaged over 86 snapshots extracted from a trajectory obtained using hybrid Carr
Parrinello molecular dynamics simulation. First we consider the results obtained using the simpler QM region model in which only the DTBNO molecule is included in the QM region and the aqueous environment is modeled using the polarizable continuum model or the discrete MM-X, X = 0
4, models (see Figure 1). In agreement with previous findings, the DFT/PCM method fails to predict the electronic g-tensor of DTBNO in aqueous solution; this method gives an isotropic g-tensor shift, Δgiso, that is 3853 ppm for the B3LYP functional and 3741 ppm for BP86 functional, respectively. These results are in rather good agreement with our previous static DFT/PCM results (with Δgiso being 3744 and 3526 ppm for B3LYP and BP86 functionals) but severely overestimates the experimental electronic g-tensor shift, which is equal to 3241 ( 10 ppm. The poor performance of DFT/PCM or DFT/COSMO methods is due to the electrostatic nature of the polarizable or conductorlike continuum models, which prevents these models from describing the hydrogen bonding between solute and protic solvent molecules without explicit inclusion of solvent molecules into the QM region. The major part of solvent shift of Δgiso originates from the largest, namely Δgzz, component of the electronic g-tensor shift, and consequently the mechanisms responsible for the solvent shift of Δgzz will determine the overall behavior of the electronic g-tensor of DTBNO in aqueous solution. In particular, the solvent blue shift of the (n) f (π*) excitation between the (SOMO
1) and (SOMO) orbitals is important in this context as it almost fully amounts to the solvent shift, according to our previous findings.43 The widely employed DFT/MM approach, in which the MM region is described by a collection of point charges (the MM-0 model in Figure 1), and in which the interaction between the QM and MM regions is described only by means of Coulomb electrostatic interaction, predicts Δgiso to be 3653 ppm (B3LYP) or 3592 ppm (BP86), which gives better agreement with experimental Δgiso by 197 ppm or by 149 ppm compared to DFT/PCM results. Therefore, the introduction of a structured environment in electronic gtensor calculations when only the solute is treated at the QM level provides an improved description of the solvent shift of Δgiso of DTBNO in aqueous solution, and this finding is expected to be valid also for other nitroxides solvated in protic solvents. A similar superior performance of the QM/MM methods with the MM region described by point charges, over 4353

dx.doi.org/10.1021/jp1108653 |J. Phys. Chem. B 2011, 115, 4350–4358

The Journal of Physical Chemistry B

ARTICLE

Table 1. Electonic g-Tensor of Di-tert-butyl Nitroxide in Water Computed with DFT/PCM and DFT/MM Methods Using Two QM Regions, Only Di-tert-butyl Nitroxide and Di-tert-butyl Nitroxide with Two Water Molecules Hydrogen Bonded to NO Bonda,b,c DTBNO method

Δgxx

Δgyy

Δgzz

Δgiso

method

Δgxx

Δgyy

Δgzz

Δgiso

B3LYP

4217


60

8133

4096

BP86

3924


38

7724

3870

B3LYP/PCM

4079


65

7544

3853

BP86/PCM

3839


38

7422

3741

B3LYP/MM-0

3954


67

7082

3656

BP86/MM-0

3721


38

7094

3592

B3LYP/MM-1

3932


67

7068

3644

BP86/MM-1

3709


37

7085

3585

B3LYP/MM-2 B3LYP/MM-3

3854 3817


69
68

6823 6765

3536 3505

BP86/MM-2 BP86/MM-3

3658 3615


38
36

6936 6881

3519 3487

B3LYP/MM-4

3819


68

6729

3493

BP86/MM-4

3636


40

6819

3471

DTBNO þ 2H2O Δgxx

Δgyy

Δgzz

B3LYP

4006


84

7146

3690

B3LYP/PCM

3887


87

6804

3538

B3LYP/MM-0

3881


88

6764

B3LYP/MM-1

3878


88

B3LYP/MM-2 B3LYP/MM-3

3848 3834

B3LYP/MM-4

method

exp.d

Δgiso

Δgxx

Δgyy

Δgzz

BP86

3786


37

6970

3573

BP86/PCM

3699


49

6780

3477

3519

BP86/MM-0

3685


49

6732

3456

6757

3516

BP86/MM-1

3680


49

6725

3452


89
89

6676 6648

3478 3464

BP86/MM-2 BP86/MM-3

3654 3641


51
52

6675 6651

3427 3413

3833


89

6632

3459

BP86/MM-4

3639


52

6648










3241 ( 10

exp.d










method

Δgiso

3412 3241 ( 10

a Electronic g-tensor shift principal values, Δgii (i = x, y, z), are averaged over the 86 snapshots taken from Carr
Parrinello molecular dynamics simulations. b Electronic g-tensor calculations have been carried out using the Huz-III basis set. c In the DFT/MM calculations the MM region includes all water molecules within a 20 Å radius sphere around di-tert-butyl nitroxide. d Experimental data are taken from the work of Kawamura et al.58 EPR measurements are carried out in water at 20 ( 3 °C.

the QM/(PCM or COSMO) methods, has been observed in ab initio studies of electronic g-tensors of blue copper proteins.17 However, we would like to point out here that despite the improved accuracy of DFT/MM-0 over the DFT/PCM method, it still overestimates the experimental measured electronic g-tensor shift by around 11
13% depending on the exchangecorrelation functional used in the calculations. After establishing the performance of the simplest DFT/MM-0 method, we turn to the DFT/MM-1 method which uses more sophisticated descriptions of the MM region (see Figure 1); i.e., water molecules are represented not only by three point charges but also by an isotropic polarizabilty. This method predicts Δgiso to be 3644 ppm (B3LYP) and 3585 ppm (BP86) and therefore provides only a marginal improvement over MM-0 results that indicate that similar accuracy of the interaction between QM and MM regions is achieved by using explicitly isotropic polarizabilities in the force field parametrization (MM-1 model) or by using point charges adjusted to a solid state environment; i.e., polarization effects are accounted for implicitly in the force field parametrization (MM-0 model). Taking this into account, the next natural step of the MM description is to introduce “more structure” by substituting the isotropic polarizability in the force field with a distributed set of anisotropic polarizability tensors (see MM-2 model in Figure 1), which leads to the DFT/MM-2 method. Going from DFT/MM-1 to DFT/ MM-2, the isotropic g-tensor shifts become more accurate by around 60
100 ppm, where a pronounced improvement is obtained for the B3LYP functional due to the nonvanishing polarization/induction contribution to the linear response

equations from the perturbed electron density, which is equal to zero for pure DFT functionals, like BP86, in the case of imaginary perturbations. Thus, the inclusion of the polarization/ induction interaction between the QM and MM regions improves the theoretical to experimental agreement, and at the DFT/MM-2 level of theory we are able to reduce the overestimation of the experimental Δgiso value to around 9%. The last two flavors of the DFT/MM method, with more accurate descriptions of the electrostatic potential and polarizability of the water molecules in the MM region compared to the MM-2 model (for details, see Figure 1), provide only limited improvement of the computed isotropic g-tensor shift of DTBNO, with values of 3505 ppm (B3LYP) and 3487 ppm (BP86) for DFT/ MM-3 method, and 3493 ppm (B3LYP) and 3471 ppm (BP86) for DFT/MM-4 method, respectively. Therefore, going from the MM-2 to the MM-3 model, which includes an addition of two new midbond expansion points for the polarizability and a multipole charge expansion up to quadrupoles instead of point charges, only a slight improvement of the solvent shift is obtained. The discrepancy between theoretical and experimental is then reduced roughtly by 30 ppm. Furthermore, the MM-3 model which truncates the multipole expansion at octopoles (see MM-4 model in Figure 1) has negligible effect on the solvent shift; Δgiso changed only by 12
16 ppm depending on exchangecorrelation functional used in calculations (see DFT/MM-3 and DFT/MM-4 results in Table 1). Taking this into account, we conclude that a converged representation of the electrostatic potential of the water molecules in the MM region is already achieved by using a distributed multipole expansion terminated 4354

dx.doi.org/10.1021/jp1108653 |J. Phys. Chem. B 2011, 115, 4350–4358

The Journal of Physical Chemistry B at the quadrupoles. This finding is furthermore in good agreement with observations made in investigations of other molecular properties in which an accurate representation of the electrostatic potential of the solvent molecules is typically achieved using a distributed multipole expansion truncated at the quadrupole level. Summarizing this part, with only DTBNO included into the QM region, Δgiso is predicted with increased accuracy going from DFT/PCM to DFT/MM in the following order: DFT/PCM < DFT/MM-0 ≈ DFT/MM-1 < DFT/MM-2 < DFT/MM-3 ≈ DFT/MM-4, where DFT/MM-3 and DFT/MM-4 overestimate the experimental Δgiso by around 230
260 ppm depending on the exchange-correlation functional. Thus already with a simplest possible model for the QM region including only the solute, the DFT/MM method combined with an accurate MM model has a potential to predict electronic g-tensor of radicals in protic solvents with an accuracy comparable to ordinary DFT methods for molecules in gas phase or noble gas matrices. A more sophisticated model for the QM region includes, in addition to the DTBNO molecule, also the two water molecules that bind by hydrogen bonds to the NO oxygen of DTBNO. This type of model for the QM region is conventionally used in connection with DFT/PCM or DFT/COSMO studies of molecules solvated in protic solvents or incorporated in proteins. The idea is to include solvent molecules required to describe the dominant hydrogen bonds between solvent and solute. The DFT/PCM method in combination with this model for the QM region has been successfully applied to investigate electronic gtensors and hyperfine coupling constants of organic radicals in various protic solvent by Barone et al.2,3,61,62 and other researchers.43,44 The main reason for the good performance of this approach is that it leaves only the bulk part of the solvent, beyond the solvent molecules strongly hydrogen bonded to solute, to be accounted for by the continuum model, where the main contribution comes from long-range electrostatic effects which can be accurately described by the PCM or COSMO models. For example, in our case the addition of two water molecules to the DTBNO radical in the QM region already in vacuum calculations allows to obtain a Δgiso equal to 3690 ppm (B3LYP) and 3573 ppm (BP86), which is in better agreement with experimental data than DFT/PCM results for the simpler QM region with single DTBNO radical in it. Thus, the hydrogen bonding is responsible for the larger part of the solvent shift of the DTBNO electronic gtensor as it exerts influence on the g-tensor via all three mechanisms discussed in the beginning of this section. Now wrapping the QM region into a dielectric continuum using the PCM model, we further improve the accuracy and Δgiso becomes equal to 3538 and 3533 ppm for B3LYP and BP86 functionals, respectively. Thus, within the limits of the polarizable continuum model, the longrange electrostatic effects amounts to around 23
27% of the total solvent shift computed using the DFT/PCM method. This result is in good agreement with our previous findings43 for DTBNO in methanol and water solutions. Contrary to the DFT/MM-0 results obtained using the simplest QM region model which includes only DTBNO, the DFT/MM-0 method combined with the more extended QM region model only slightly outperforms the DFT/PCM method, where more pronounced gain in accuracy going from DFT/PCM to DFT/MM-0 is obtained for the B3LYP functional. Thus, the better performance of the DFT/MM-0 method over DFT/PCM using the simpler QM region is caused mainly by the lack of a structured environment in the PCM model, which is an inherent

ARTICLE

limitation of this model. Beyond DFT/MM-0, adding to the Coulomb electrostatic interaction term between the QM and MM regions also the polarization/induction interaction term (see the DFT/MM-1 model in Figure 1), the isotropic g-tensor shift remains almost the same for both functional (see DFT/ MM-0 and DFT/MM-1 results in Table 1). This result shows that the polarization effects of the QM region by the MM region is described almost at the same level by the TIP3P force field implicitly (DFT/MM-0) and by the Ahlstr€om force field explicitly (DFT/MM-1). Further improvement of the Δgiso is obtained using the distributed anisotropic representation of the MM water polarizability, i.e., the MM-2 model, in DFT/MM calculations, and within this approach the Δgiso becomes 3478 and 3427 ppm for B3LYP and BP86 functionals, respectively. Based on these results, we can decompose the solvent shift of Δgiso of DTBNO in aqueous solution into two parts (assuming the additivity of contributions): (a) explicit hydrogen bonding contribution varies from around 300 to around 400 ppm depending on exchange-correlation functional; (b) sum of electrostatic Coulomb interaction and polarization/induction contributions equal to 212 ppm (B3LYP) or 146 ppm (BP86). Here, we find that the Coulomb electrostatic and polarization/ induction contributions to the solvent shift from the MM region; using both functionals is significant and amounts up to 34% and 49% of the solvent shift for B3LYP and BP86 functionals. Furthermore, we find that DFT/MM-2 with the simpler QM region predicts the Δgiso to be 3536 ppm (B3LYP) or 3519 ppm (BP86) in rather good agreement with the results for DFT/ MM-2 with the extended QM region, which are 3478 and 3427 ppm for the B3LYP and BP86 functionals, respectively. The improved description of the electrostatic potential and the polarizability of water in the MM region introduced by the MM-3 and MM-4 models over the MM-2 model leads to negligible changes in Δgiso and its solvent shift indicating a smooth convergence of the Δgiso going from MM-2 to MM-4. For example, the DFT/MM-3 and DFT/MM-4 methods combined with the B3LYP functional predict the Δgiso to be 3464 and 3459 ppm, respectively. Thus, DFT/MM-4 combined with the QM region, which includes DTBNO and two hydrogen bonded water molecules, overestimates the measured Δgiso of DTBNO in aqueous solution by 218 ppm (B3LYP) or 171 ppm (BP86), where slightly better performance of the BP86 functional is observed. These results are only slightly more accurate than the one obtained with DFT/MM-4 using the simpler QM solute model, which includes only DTBNO, and therefore we conclude that the difference between the two QM models diminishes going from simpler water models like MM-1 to more advanced models like MM-4. We conclude that MM-3 and MM-4 allow to achieve almost the same accuracy as an explicit treatment of hydrogen bonded waters in the evaluation of the electronic g-tensor that is within the accuracy of description of hydrogen binding by the B3LYP functional. This finding can have significant ramifications for practical calculations of electronic g-tensors of organic radicals in more complex protic solvents, their mixtures, and proteins, as it will allow to restrict the QM region to the radical of interest and treat the remaining part of the system using an (advanced) MM model in which solvent molecules or the protein residues are represented via distributed electrostatic potentials (truncated at quadrupoles) and distributed anisotropic polarizability tensors. We are now in the position to carry out research along this direction to verify these findings for various nitroxide spin labels in proteins and 4355

dx.doi.org/10.1021/jp1108653 |J. Phys. Chem. B 2011, 115, 4350–4358

The Journal of Physical Chemistry B

ARTICLE

Table 2. Electonic g-Tensor of Di-tert-butyl Nitroxide Solvated in Water: Static and Dynamic DFT Calculation Results modela

typeb

Δgxx

Δgyy

Δgzz

B3LYP/PCMc

DTBNO þ 2H2O

dynamic

3887


87

6804

B3LYP/MM-4c

DTBNO

dynamic

3819


68

6729

3493

B3LYP/MM-4c

DTBNO þ 2H2O

dynamic

3833


89

6632

3459

B3LYP/PCMd,e

DTBNO þ 2H2O

static

3624


78

6487

3344

B3LYP/COSMO f,g

DTBNO þ 2H2O

static

3609


198

6276

3229

B3LYP/D-COSMO-RSf,g

DTBNO þ 2H2O

static










3240

BP86/PCMc

DTBNO þ 2H2O

dynamic

3894


88

6793

3533

BP86/MM-4c BP86/MM-4c

DTBNO DTBNO þ 2H2O

dynamic dynamic

3636 3639


40
52

6819 6648

3471 3412

BP86/PCMd,e

DTBNO þ 2H2O

static

3458


58

6176

3192

BP86/COSMO f,g

DTBNO þ 2H2O

static

3394


206

5968

3052

BP86/D-COSMO-RSf,g

DTBNO þ 2H2O

static










3027

method

Δgiso 3538

3241 ( 10

Exp.h a

QM region model used in calculations: single di-tert-butyl nitroxide denoted as DTBNO and di-tert-butyl nitroxide with two waters hydrogen bonded to NO bond denoted as DTBNO þ 2H2O. bType of calculations: “dynamic” stands for electronic g-tensor calculations carried out over 86 snapshots taken from Carr
Parrinello molecular dynamics simulation, and “static” stands for single point electronic g-tensor calculation over the optimized QM region geometry. cElectronic g-tensor shifts computed in this work. dElectronic g-tensor shifts taken from Rinkevicius et al.43 eDi-tert-butyl nitroxide geometry optimized in vacuum at the B3LYP/6-31G(d,p) level has been used in g-tensor calculations. fElectronic g-tensor shifts taken from Sinnecker et al.44 g Di-tert-butyl nitroxide geometry optimized in solvent at the B3LYP/6-31G(d,p) level has been used in g-tensor calculations. hExperimental data are taken from the work of Kawamura et al.58 EPR measurements are carried out in water at 20 ( 3 °C.

membranes, where the direct inclusion of all hydrogen bonded residues to the spin label in the QM region rapidly becomes computationally expensive. Summarizing, the quality of the DFT/MM results obtained with the QM region, which includes DTBNO and two waters hydrogen bonded to it, give the order of the models as Δgiso: DFT/PCM < DFT/MM-0 ≈ DFT/MM-1 < DFT/MM-2 ≈ DFT/MM-3 ≈ DFT/MM-4, with respect to increasing ability to predict the isotropic g-tensor shift of DTBNO. Here, compared to the previously given order of DFT/MM and DFT/PCM methods for the QM solute only model, two important changes are seen: first, DFT/MM-0 is similar to DFT/PCM indicating that the PCM approach is well suited for a description of the protic solvent provided that the most important solvent molecules interacting with the solute are included in the QM region; second, advanced models of the solvent molecules converge more rapidly with the extended QM region model, since the solvent description beyond the first solvation sphere puts smaller demands on the accuracy of the electrostatic potential and polarizability of the individual solvent molecules. Taking these findings into account, we conclude that the sophistication of the description of the MM region can be varied depending on the complexity of the QM region, where the simpler QM region models require a more sophisticated MM description. We hope that these findings will be helpful for practical setup of calculations of DFT/MM of EPR spin Hamiltonian parameters as well as other molecular properties. C. Role of Local Solvent Dynamics. The analysis given above effectively includes the local solvent dynamics via averaging over a set of snapshots obtained from hybrid CPMD simulations. In order to present the complete picture of the electronic g-tensor in aqueous solution, we will compare those results with results from static DFT/PCM and DFT/COSMO calculations by Rinkevicius et al.43 and by Sinnecker et al.,44 respectively; see Table 2. We note that the previous “static” results are obtained using a solvent optimized geometry of the selected QM region, which corresponds to a hypothetical structure at 0 K temperature, and that the “dynamic” calculations account for local solvent

dynamics at 298 K temperature. We estimate the effects of the solvent dynamics on the results by comparing with the DFT/ PCM results. Assuming the extended QM region model, one can estimate that the change of the static to the dynamic picture increases Δgiso by around 200 ppm independently on the exchange-correlation functional used in our calculations. Therefore, a comparison between static and dynamic results is not entirely appropriate and can be made only on qualitative grounds. Nevertheless, a quick inspection of Table 2 indicates that the static DFT/COSMO and DFT/PCM calculations with the QM region that includes two water molecules predict on average the electronic g-tensor in better agreement with experiment than the DFT/MM-4 results for both B3LYP and BP86 functionals. In view of this finding, one should ask if it is worth the effort to use a combined approach which involves averaging of the g-tensor component over the molecular dynamics trajectory carried out using hybrid DFT/MM or DFT/PCM methods. The answer to this question is nontrivial—on one hand the static DFT/PCM calculations can seemingly reproduce the experimental measurements with better accuracy than more sophisticated dynamic DFT/PCM or DFT/MM calculations at the same time as they show strong dependence on the exchange-correlation functional used; on the other hand, DFT/MM or DFT/ PCM calculations in which the local dynamics of the solvent at finite temperature is accounted for are physically more appealing and show less pronounced dependence of the results on the exchange-correlation functional. Furthermore, in recent studies of electronic g-tensors of nitroxides, Pavone and co-workers demonstrated that dynamical effects and electronic g-tensors of prototypical nitroxides exhibit a strong dependence on the NO bond length,61 which can be hardly accounted for in a static DFT/PCM or “DFT/COSMO” calculation. Our results confirm the observations made by Pavone et al.61 that Δgiso shows almost no dependence on the value of the out-of-plane angle Θ but exhibits a rather strong dependence on the NO bond length. Taking this into account, the ability of an averaged structure used in static DFT/PCM or DFT/COSMO calculations becomes 4356

dx.doi.org/10.1021/jp1108653 |J. Phys. Chem. B 2011, 115, 4350–4358

The Journal of Physical Chemistry B doubtful, and probably the better agreement with experimental data achieved in these calculations is due to a fortuitous cancellation of errors. In this context, one finds that benchmarking studies of various exchange-correlation functionals show that the BP86 functional is more accurate than the B3LYP functional in the DFT calculations of electronic g-tensors of organic radicals48
50 and the same trend is observed in this work, but the opposite has been found for DTBNO in various solvents in our previous DFT/PCM study,43 where the BP86 functional underestimated the experimental Δgiso compared to B3LYP functional. Therefore, it is obvious that further investigations of local solvent dynamics effects in various protic solvents must be carried out in order to establish more rigorously the performance of various exchange-correlation functionals in computations of electronic g-tensors of organic radicals as well as to estimate the importance of dynamic effects for this EPR spin-Hamiltonian parameter.

V. CONCLUSIONS The hybrid density functional theory/molecular mechanics approach, which allows for an advanced treatment of the coupling between quantum mechanical (QM) and molecular mechanical (MM) regions of a system under study, has been applied for the first time to investigate solvent shifts of electronic g-tensors in aqueous solution and is used to benchmark the accuracy of various models to construct the MM region of the total system. The DFT/MM method is capable of consistently accounting for electrostatic Coulomb and polarization/induction interactions between the QM and MM parts for linear and nonlinear response properties, where the contributions due to the MM region are included in both the Kohn
Sham and response equations. The MM region can be treated at different levels of sophistication: the electrostatic potential of the solvent molecule can be reconstructed using distributed multipole expansions in which an arbitrary set of expansion points can be selected and multipoles up to octopoles can be included; the polarizability of the solvent molecule can be represented by a simple isotropic polarizability or a set of anisotropic polarizability tensors located on various expansion points in the molecules. Results for ditert-butyl nitroxide (DTBNO) in aqueous solution show that the polarization/induction interaction, which is typically not included in conventional implementations of QM/MM methods, provides a significant contribution, around one-third of the solvent shift of gtensor, due to the MM region, and must therefore be included in practical DFT/MM calculations along with a “realistic” representation of the polarizability of the MM molecules. Even more important, the results indicate that with accurate models of the solvent molecules, which describe the electrostatic potential of the solvent molecule via a distributed multipole expansion truncated at quadrupoles and its polarizability via a distributed set of anisotropic polarizability tensors, allows to effectively reproduce the structured environment of the solute in protic solvents and reduce the QM region to the solute itself, thus avoiding explicit inclusion of solvent molecules in the QM region. This result can be of profound importance for large scale DFT/MM simulations of spin labels and other radicals in proteins and cellular membranes, as it opens the way for the use of a restricted QM region in calculations of EPR spin Hamiltonian parameters as well as other molecular properties. That in a turn makes calculations of large nonhomogeneous systems possible provided an acceptable quality of the MM region description is achieved. Here, we would like to point out that our DFT/ MM method, which allows to include electrostatic Coulomb and

ARTICLE

polarization/induction interactions between the QM and MM regions, clearly demonstrates the importance of various interactions between these two regions of the system beyond the simple electrostatic picture. This will focus further developments into a direction for an accurate account of long and short-range van der Waals interactions and charge transfer between solvent and solute, which is needed in order to expand the applicability of DFT/MM to systems like ionic liquids, as well as to increase its accuracy for protic solvents and protein environments. Systematic studies of DFT/MM for other EPR spin Hamiltonian parameters—hyperfine coupling constants and zero-field splitting parameters—of organic radicals in various environments will follow after this work.

’ AUTHOR INFORMATION Corresponding Author

*E-mail [email protected].

’ ACKNOWLEDGMENT This work was partially funded by the EU Commission (contract INFSO-RI-261523) under the ScalaLife collaboration and was supported by a grant from the Swedish Infrastructure Committee (SNIC) for the project “Multiphysics Modeling of Molecular Materials”, SNIC 022/09-25. J.K. thanks The Danish Councils for Independent Research and the Lundbeck Foundation for financial support. We thank Xing Chen for technical assistance in preparation of figures. ’ REFERENCES (1) Neese, F. Curr. Opin. Chem. Biol. 2003, 7, 125–135. (2) Improta, R.; Barone, V. Chem. Rev. 2004, 104, 1231–1253. (3) Barone, V.; Cimino, P.; Pedone, A. Magn. Reson. Chem. 2010, 48, S11–S22. (4) Fernandez, B.; Christiansen, O.; Bludsky, O.; Jørgensen, P.; Mikkelsen, K. V. J. Chem. Phys. 1996, 104, 629–635. (5) Brancato, G.; Rega, N.; Barone, V. J. Am. Chem. Soc. 2007, 129, 15380–15390. (6) Ciofini, I.; Reviakine, R.; Arbuznikov, A.; Kaupp, M. Theor. Chem. Acc. 2004, 111, 132–140. (7) Asher, J. R.; Kaupp, M. Theor. Chem. Acc. 2008, 119, 477–487. (8) Engstr€om, M.; Vaara, J.; Schimmelpfennig, B.; Ågren, H. J. Phys. Chem. B 2002, 106, 12354–12360. (9) Svistunenko, D. A.; Jones, G. A. Phys. Chem. Chem. Phys. 2009, 11, 6600–6613. (10) Bernini, C.; Sinicropi, A.; Basosi, R.; Pogni, R. Appl. Magn. Reson. 2010, 37, 279–288. (11) Pauwels, E.; Van Speybroeck, V.; Waroquier, M. Spectrochim. Acta A 2006, 63, 795–801. (12) Pauwels, E.; Verstraelen, T.; De Cooman, H.; Van Speybroeck, V.; Waroquier, M. J. Phys. Chem. B 2008, 112, 7618–7630. (13) Neugebauer, J.; Louwerse, M. J.; Belanzoni, P.; Wesolowski, T. A.; Baerends, E. J. J. Chem. Phys. 2005, 123, 114101. (14) Moon, S.; Patchkovskii, S.; Salahub, D. R. J. Mol. Struct. Theochem 2003, 632, 287–295. (15) Sinnecker, S.; Neese, F. J. Comput. Chem. 2006, 27, 1463–1475. (16) Houriez, C.; Ferre, N.; Masella, M.; Siri, D. J. Chem. Phys. 2008, 128, 244504. (17) Vancoillie, S.; Pierloot, K. J. Phys. Chem. A 2008, 112, 4011– 4019. (18) Pauwels, E.; Declerck, R.; Verstraelen, T.; De Sterck, B.; Kay, C. W. M.; Van Speybroeck, V.; Waroquier, M. J. Phys. Chem. B 2010, 114, 16655–16665. (19) Houriez, C.; Ferre, N.; Siri, D.; Masella, M. J. Phys. Chem. B 2009, 113, 15047–15056. 4357

dx.doi.org/10.1021/jp1108653 |J. Phys. Chem. B 2011, 115, 4350–4358

The Journal of Physical Chemistry B (20) Houriez, C.; Ferre, N.; Masella, M.; Siri, D. J. Mol. Struct. Theochem 2009, 898, 49–55. (21) Aidas, K.; Møgelhøj, A.; Kjær, H.; Nielsen, C. B.; Mikkelsen, K. V.; Ruud, K.; Christiansen, O.; Kongsted, J. J. Phys. Chem. A 2007, 111, 4199–4210. (22) Aidas, K.; Møgelhøj, A.; Nilsson, E. J. K.; Johnson, M. S.; Mikkelsen, K. V.; Christiansen, O.; S€oderhjelm, P.; Kongsted, J. J. Chem. Phys. 2008, 128, 194503. (23) Adamovic, I.; Freitag, M. A.; Gordon, M. S. J. Chem. Phys. 2003, 118, 6725–6732. (24) Nielsen, C. B.; Christiansen, O.; Mikkelsen, K. V.; Kongsted, J. J. Chem. Phys. 2007, 126, 154112. (25) Olsen, J. M.; Aidas, K.; Kongsted, J. J. Chem. Theory Comput. 2010, 6, 3721–3734. (26) Yoo, S.; Zahariev, F.; Sok, S.; Gordon, M. S. J. Chem. Phys. 2008, 129, 144112. (27) Jensen, L.; van Duijnen, P. T.; Snijders, J. G. J. Chem. Phys. 2003, 119, 12998–13006. (28) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (29) Ahlstr€om, P.; Wallqvist, A.; Engstr€om, S.; J€onsson, B. Mol. Phys. 1989, 68, 563–581. (30) Worthington, S. E.; Roitberg, A. E.; Krauss, M. J. Phys. Chem. B 2001, 105, 7087–7095. (31) Ramalho, T. C.; da Cunha, E. F. F.; de Alencastro, R. B. J. Phys.: Condens. Matter 2004, 16, 6159–6170. (32) Becke, A. D. Phys. Rev. A 1988, 38, 3098–3100. (33) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37 785–789. (34) Troullier, N.; Martins, J. L. Phys. Rev. B 1991, 43, 1993– 2006. (35) Nose, S. J. Chem. Phys. 1984, 81, 511–519. (36) Hoover, W. G. Phys. Rev. A 1985, 31, 1695–1697. (37) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157–1174. (38) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361–373. (39) Frisch, M. J.; et al. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (40) Hutter, J.; et al. Computer code CPMD, version 3.11, copyright IBM Corp. and MPI-FKF, Stuttgart, Germany, 1990
2002. (41) Case, D. A.; et al. AMBER 8, University of California, San Francisco, CA, 2004. (42) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Åstrand, P.-O.; Christiansen, O. J. Chem. Phys. 2004, 121, 8435–8445. (43) Rinkevicius, Z.; Telyatnyk, L.; Vahtras, O.; Ruud, K. J. Chem. Phys. 2004, 121, 5051–5060. (44) Sinnecker, S.; Rajendran, A.; Klamt, A.; Diedenhofen, M.; Neese, F. J. Phys. Chem. A 2006, 110, 2235–2245. (45) Perdew, J. P. Phys. Rev. B 1986, 33, 8822–8824. (46) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (47) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200–1211. (48) Malkina, O. L.; Vaara, J.; Schimmelpfennig, B.; Munzarova, M.; Malkin, V. G.; Kaupp, M. J. Am. Chem. Soc. 2000, 122, 9206–9218. (49) Neese, F. J. Chem. Phys. 2001, 115, 11080–11096. (50) Rinkevicius, Z.; Telyatnyk, L.; Sazek, P.; Vahtras, O.; Ågren, H. J. Chem. Phys. 2003, 119, 10489–10496. (51) Ruden, T. A.; Lutnæs, O. B.; Helgaker, T.; Ruud, K. J. Chem. Phys. 2003, 118, 9572–9581. (52) Huzinaga, S. J. Chem. Phys. 1965, 42, 1293–1302. (53) Kutzelnigg, W.; Fleischer, U.; Schindler, M. NMR Basic Principles and Progress; Springer Verlag: Berlin, Germany, 1990; Vol 23, pp 165
262. (54) Van W€ullen, C. Ph.D. Thesis, Ruhr-Universit€at: Bochum, Germany, 1992. (55) Hess, B. A.; Marian, C. M.; Wahlgren, U.; Gropen, O. Chem. Phys. Lett. 1996, 251, 365–371.

ARTICLE

(56) DALTON, a molecular electronic structure program, release 2.0, 2005; see www.daltonprogram.org. (57) Pavone, M.; Cimino, P.; De Angelis, F.; Barone, V. J. Am. Chem. Soc. 2006, 128, 4338–4347. (58) Kawamura, T.; Matsunami, S.; Yonezawa, T. Bull. Chem. Soc. Jpn. 1967, 40, 1111–1115. (59) Stone, A. J. Mol. Phys. 1963, 6, 509–515. (60) Stone, A. J. Mol. Phys. 1964, 7, 311–316. (61) Pavone, M.; Cimino, P.; Crescenzi, O.; Sillanp€a€a, A.; Barone, V. J. Phys. Chem. B 2007, 111, 8928–8939. (62) Stendardo, E.; Pedone, A.; Cimino, P.; Cristina Menziani, M.; Crescenzi, O.; Barone, V. Phys. Chem. Chem. Phys. 2010, 12, 11697– 11709. (63) Gagliardi, L.; Lindh, R.; Karlstr€om, G. J. Chem. Phys. 2004, 121, 4494–4500.

4358

dx.doi.org/10.1021/jp1108653 |J. Phys. Chem. B 2011, 115, 4350–4358