- Email: [email protected]

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Solvent effects on the nitrogen NMR shielding and nuclear quadrupole coupling constants in 1-methyltriazoles Andreas Møgelhøj a, Kestutis Aidas a, Kurt V. Mikkelsen a, Jacob Kongsted b,* a b

Department of Chemistry, H. C. Ørsted Institute, University of Copenhagen, Universitetsparken 5, DK-2100 Copenhagen Ø, Denmark Department of Theoretical Chemistry, Chemical Center, University of Lund, P.O. Box 124, S-221 00 Lund, Sweden

a r t i c l e

i n f o

Article history: Received 13 May 2008 In ﬁnal form 3 June 2008 Available online 7 June 2008

a b s t r a c t In this Letter, we use a discrete polarizable solvation model for a systematic analysis of the solvent effects on the nitrogen NMR shielding and nuclear quadrupole coupling constants in a series of 1-methyltriazoles. Fairly accurate predictions are found for the solvent shifts of the nitrogen NMR shielding constants. The analysis of the relative half-height widths of the resonance signal predicted in either vacuum or aqueous solution implies that the spin relaxation time for the pyridine- and pyrrole-type nitrogen atoms possess similar magnitudes in vacuum whereas they are different in aqueous solution. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction Nuclear magnetic resonance (NMR) spectroscopy is an extremely valuable method for obtaining a detailed microscopic picture of the structure and dynamics of molecular systems [1], and has become one of the most prized tools in the ﬁelds of chemistry and biology. From the theoretical side, advanced computational protocols have been developed for calculating NMR parameters of isolated molecules [2], and it is now well known that inclusion of electron correlation and effects of nuclear motion may be mandatory to consider in order for theory to match the accuracy of experiments. For molecules in the condensed phase it is, however, only recently that ab initio modelling has approached a level with a similar predictive power. Nevertheless, many important questions concerning environmental effects on NMR parameters are still not understood. In many respects, this is related to the inherently large size of such systems making conventional ab initio theories either very computationally demanding or even prohibited. Thereby, different compromises between accuracy and computational cost have been, and still are, developed and explored. In many studies only a limited part of the environment/solvent is treated explicitly. However, such predictions made on the basis of calculations on molecular aggregates may suffer from at least two serious errors: (i) neglect of long-range electrostatic effects and (ii) lack of dynamical averaging. In addition, the geometry of the molecular aggregates is usually obtained from energy minimization procedures which is not in line with the true nature of a solvent. In our opinion, a more physically sound protocol would be to

* Corresponding author. Fax: +46 46 222 4543. E-mail address: [email protected] (J. Kongsted). 0009-2614/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2008.06.004

explore the concepts of effective Hamiltonians, i.e. to treat only a limited part of the total system at the level of ab initio theory and to invoke a more coarse-grained model for the remaining, and largest, part of the system. Along this line two distinct models have been developed. In the ﬁrst of these two approaches the solvent molecules are described as a dielectric continuum thereby neglecting any reference to the atomistic nature of the solvent. The most advanced model belonging to this class is the polarizable continuum model (PCM) [3]. Possible disadvantages associated with this model may be related to a neglect of local anisotropies around the solute, and furthermore to the deﬁnition of the cavity which contains the solute inside the dielectric medium. The former point may, however, be improved upon treating a number of solvent molecules explicitly. In the second class of effective solvent models an explicit reference to the discrete solvent is kept intact. Here, the solvent is described using molecular mechanics (MM) while the solute, and potentially some of the solvent molecules, are described using quantum mechanics (QM). This class of solvent models is generally referred to as combined quantum mechanics/molecular mechanics (QM/MM) [4–7]. In the simplest version of this approach, the solvent is treated by assigning partial point charges to the atomic sites and the electric potential due to these point charges is then introduced into the solute Hamiltonian. However, in such a procedure polarization of the solvent is neglected, i.e. only the solute is polarized. This may be reﬁned, for example, by assigning polarizable sites to the solvent giving rise to induced electric moments. Since the QM/MM model contains an atomistic reference to the solvent molecules, an explicit averaging needs to be included. This is conventionally achieved by considering a number of solute–solvent conﬁgurations based either on molecular dynamics (MD) or Monte Carlo (MC) procedures in the actual QM/MM calculations.

130

A. Møgelhøj et al. / Chemical Physics Letters 460 (2008) 129–136

In this study we use the QM/MM model in order to address solvent effects on two NMR related properties – the NMR isotropic shielding constant and the nuclear quadrupole coupling constant (NQCC). The QM/MM model we use contains an explicit term accounting for the electronic polarization of the solvent. Thus, we describe each atom in the classically treated part of the system by a point charge, an isotropic dipole–dipole polarizability and a set of Lennard–Jones parameters. However, in the actual QM/MM calculations only the point charges and polarizabilities perturb the electronic structure in question. The solute–solvent conﬁgurations considered in the QM/MM calculations are derived from classical MD simulations. In some cases we will also explore the possibility of including solvent effects using PCM. In the QM/MM calculations the ab initio level will be restricted to density functional theory (DFT). Thereby, we abbreviate our hybrid model by DFT/MM. In particular, special attention will be payed to the choice of exchange-correlation (xc) functional. The molecular systems considered in this Letter belong to the class of 1-methyltriazoles (see the chart below which also deﬁnes molecular denominations and the atomic labelling used throughout this Letter). N3

N2

C2H N1

C2H

N4

1-Me-1,2,3-triazole denoted <123>

C1H

N2 CmH3

C1H

C1H N1

N4

N5

1-Me-1,2,5-triazole denoted <125>

CmH3

1-Me-1,2,4-triazole denoted <124> N3

N1 C1H

N1

CmH3

C1H

N2

" K ij

r ¼1þ ¼1þ

#

2

d E K

dBi dmj X lm

B¼mK ¼0 2

Dlm

o hlm K oBi dmj

þ

X oDlm ohlm : K lm oBi omj

ð1Þ

Here Dlm is an element of the density matrix in AO basis and hlm is a matrix element of the effective one-electron Hamiltonian. The second term in Eq. (1) is the diamagnetic contribution (an expectation value) while the third term is the paramagnetic contribution, which involves the derivative of the density matrix with respect to the magnetic induction. The latter term is obtained by solving a set of linear response equations for three components of the magnetic induction [10–12]. The surrounding polarizable environment contributes implicitly to Eq. (1) through a perturbation of the density matrix. In addition, an explicit contribution enters through the derivative of the density matrix with respect to the magnetic induction [13]. In order to ensure origin invariance of the predicted NMR shielding tensors, the calculations are performed in the basis of gauge including atomic orbitals (GIAOs) [14–17]. Details concerning the theory and implementation of response theory in the context of magnetic perturbation can be found in Ref. [2], while a description of GIAO-DFT/MM can be found in Ref. [13]. The dominating relaxation mechanism for quadrupolar nuclei, i.e. nuclei possessing a spin grater than 1/2, is governed by the nuclear quadrupole moment. Assuming that molecular rotation can be treated isotropically and furthermore can be described adequately by a single relaxation time, sc, the measured line width is proportional to the inverse relaxation time [18], and the halfheight width of the resonance signal, Dm1 , may be derived from 2

CmH3

Dm1 2

C1H

1-Me-1,3,4-triazole denoted <134>

The analysis of solvent effects on the NMR shielding constants and NQCCs has been limited to the nitrogen atoms in these compounds, and we only consider water as a solvent. The choice of this set of molecules has been put forward due to the existence of experimental data for both solvent shifts in the NMR shielding constants as well as parameters related to the NQCCs. Furthermore, it is well known that accurate calculations of NMR shielding constants for nitrogen atoms may, due to the presence of lone pair and potentially multiple bonds, in fact represent a particularly difﬁcult task for ab initio computations (see for example the discussion in Ref. [8]). All four 1-methyltriazoles considered contain one pyrrole- and two pyridine-type nitrogen atoms. The pyridine-type nitrogen atoms are capable of performing hydrogen bonds with water, and solvent effects are therefore expected to be different for these two types of nitrogen atoms. In the calculations on h125i and h134i we assume that the two pyrrole-type nitrogen atoms are chemically equivalent. This means that we in this Letter consider 10 chemically nonequivalent nitrogen atoms. 2. Method 2.1. Ab intio calculation of NMR parameters The nuclear magnetic shielding tensor for nucleus K, rK, is expressed in terms of the second-order response of the electronic energy to an external magnetic induction, B, and a nuclear magnetic moment, mK. In the atomic orbital (AO) basis, the expression becomes [2,9]:

2 1 3 2IK þ 3 g2K eqK / ¼ V K;zz sc : 1 þ 2 T q 40 IK ð2IK 1Þ 3 h

ð2Þ

In this equation IK, is the nuclear spin, eqK is the scalar electric quadrupole moment of nucleus K, while VK,zz is a component of the electric ﬁeld gradient (EFG) tensor with the asymmetry parameter denoted as gK. The nuclear quadrupole coupling constant is deﬁned as eqK VK,zz. Thereby, the EFG tensor becomes the central quantity to obtain from ab initio electronic structure theory. The components of the electric ﬁeld gradient at nucleus K is by deﬁnition ordered as

jV K;zz j P jV K;yy j P jV K;xx j

ð3Þ

and the asymmetry parameter gK is deﬁned according to

gK ¼

V K;xx V K;yy : V K;zz

ð4Þ

The electric ﬁeld gradient at nucleus K is calculated as an expectation value of the operator

^K ¼ V

X 3r 2 1 rT riK X 3R2 1 RT RKL iK iK þ Z L KL 5 KL ; 5 r RKL iK L6¼K i

ð5Þ

and is thereby obtained more easily than the shielding constant. As discussed previously for the NMR shielding constant, solvent effects enter implicitly in the calculation of the expectation value of the electric ﬁeld gradient operator through a perturbation of the density matrix. In addition, an explicit contribution appears due to the charges and induced dipoles of the solvent molecules. However, this explicit contribution is easily obtained by considering the (converged) induced dipole moments and partial atomic charges deﬁning the solvent molecules [19]. 2.2. Computational protocol Since the solvent model used in this Letter relies on an explicit representation of the solvent molecules, an explicit averaging over

131

A. Møgelhøj et al. / Chemical Physics Letters 460 (2008) 129–136

a number of solute–solvent conﬁgurations is needed. In this respect, MD simulations have been performed from which a number of statistically uncorrelated solute–solvent conﬁgurations has been extracted. The NMR parameters are calculated using these conﬁgurations and the ﬁnal predictions of the spectroscopic parameters are obtained by statistical averaging. In both the MD simulations and the QM/MM calculations, the explicitly polarizable potential is used. The QM and MM subsystems are thus allowed to mutually polarize each other. This gives rise to an additional iterative procedure for obtaining self-consistency for both the electronic density and the induced dipole moments [12,20]. The MD simulations were performed using the MOLSIM program [21]. For the 1-methyltriazoles (possibly complexed with a few water molecules), geometries were obtained by optimization using B3LYP/aug-cc-pVTZ method [22,23] in combination with PCM [24–27] as implemented in the GAUSSIAN 03 suite of programs [28]. We use default settings of PCM as implemented in GAUSSIAN 03. Partial atomic point charges were calculated using the CHelpG procedure [29] at the level of B3LYP/aug-cc-pVTZ. Atomic polarizabilities were calculated at the B3LYP/aug-cc-pVTZ level using the LoProp approach [30] implemented in the MOLCAS program [31]. The Lennard–Jones parameters from Ref. [32] have been used. For water, a polarizable force ﬁeld by Ahlström et al. [33] has been used. The parameters for the 1-methyltriazole and water force ﬁelds are compiled in Table 1. From Table 1 we observe that the isotropic atomic polarizabilities for nitrogen fall in the interval of 1.0–1.2 Å3, while the isotropic carbon polarizabilities vary between 1.0 and 1.4 Å3. Even upon comparison of the same type of atoms within the different molecules, a fairly large variation in the value of the polarizability is found. This observation suggests that it might be important to calculate the polarizability for the actual atom and not rely on a common atom type polarizability. All four MD simulations where performed in a cubic box of side lengths 24.92 Å containing 1 rigid solute and 511 rigid water molecules. The box length has been chosen so as to reproduce the experimental density of water. A spherical cut-off radius equal to half of the box length was used. The simulations where performed at a temperature of 298.15 K utilizing periodic boundary conditions and a time step of 2 fs. The initial equilibration ran for 400 ps followed by a production run of 1.2 ns. Solute-solvent conﬁgurations where dumped every 10 ps. In the MD simulations we consider exclusively the 14N isotope for all nitrogen atoms. The DFT/MM calculations have been performed using the Dalton quantum chemistry program package [34]. The MidasCpp program [35] was used to generate inputs to Dalton as well as to perform the statistical analysis of the results. The DFT/MM scheme is applied for the molecular solute–solvent conﬁgurations extracted from the classical MD simulation. The solute, or potentially the solute and a number of closest solvent molecules, is considered using DFT, whereas the rest of the solvent molecules are described

by the MM potential. In the DFT/MM calculations we used the same force ﬁeld for the water molecules as was employed in the MD simulation. A spherical cut-off radius, Rcut, has been applied to all molecular conﬁgurations. In line with our previous investigations [36,37], we here use Rcut = 12 Å, which amounts to include around 240 water molecules in the DFT/MM calculations. 3. Results and discussion 3.1. Method analysis The large-scale calculations presented in this Letter will be based on DFT. In Table 2 we show the results of a basis set analysis of the NMR shielding constants and NQCCs for vacuum optimized geometries of h125i and h134i. The results in Table 2 were obtained using the Dunning basis sets aug-cc-pVXZ [23] (denoted axz) and aug-cc-pCVXZ [38] (denoted acxz), where X = D,T,Q,5. In addition, we have considered the Pople basis sets 6-311++G**, 6311++G(2d,2p) and 6-311++G(3df,3pd) [39] abbreviated by P1, P2 and P3, respectively. All calculations refer to DFT employing the B3LYP xc-functional. As seen from Table 2, a systematic decrease in the value of the isotropic shielding constants is observed in the sequences a(c)dz to a(c)5z. Since the changes observed when increasing the basis set from acqz to ac5z are fairly small, the results are most likely converged using the ac5z basis set. However, a deﬁnite conclusion concerning this issue can only be addressed by exploring the results obtained by a basis set of ac6z quality. For these compounds, the ac5z basis set amounts to 1486 contracted basis functions, and taking into consideration that in the liquid a number of conﬁgurations (>100) needs to be considered, the use of the ac5z basis set is, from a computational point of view, prohibited. The variations in the predicted NMR shielding constants using the Pople basis sets are rather modest, and provide predictions which are of atz quality, clearly representing non-converged results. Since our interest is more in the solvation shift than in the absolute values of the NMR shielding constant (which for high-level predictions also would require detailed analysis of the effects due to nuclear motions), we consider instead the solvation shift, dN. The solvent shifts have also been included in Table 2. These shifts are deﬁned as dN = rN (1-methyltriazole + 2w/ PCM)-rN (1-methyltriazole), i.e. we geometry-optimize a speciﬁc 1-methyltriazole in complex with two water molecules (all embedded in a dielectric continuum), evaluate the NMR shielding constant (the dielectric continuum was not present in the NMR calculations) and subtract from this the corresponding result obtained for the isolated 1-methyltriazole molecule. The latter has been geometry-optimized in isolated form. This means that the shifts reported in Table 2 include the effects of geometric relaxation. For these solvent shifts we observe a steady increase in magnitude using the Dunning basis set series. However, the changes are much

Table 1 The atomic parameters used in MD simulations h1XYi

Atom

N1

NX

NY

Cm

C1

C2

Hm1

Hm2

H1

H2

Ow

Hw

h123i

a (Å3)

1.012 0.5312 0.991 0.3517 0.995 0.6537 1.134 0.0615

1.045 0.3058 1.049 0.6028 1.078 0.5094 1.083 0.3659

1.101 0.3464 1.162 0.6371 1.078 0.5094 1.083 0.3659

1.005 0.2488 1.002 0.3798 1.036 0.2942 1.410 0.3756

1.261 0.4100 1.190 0.1886 1.243 0.0762 1.199 0.2516

1.276 0.2867 1.171 0.5447 1.243 0.0762 1.199 0.2516

0.362 0.1031 0.362 0.1598 0.347 0.1089 0.311 0.1460

0.335 0.0858 0.337 0.1072 0.376 0.1141 0.428 0.1491

0.378 0.1719 0.377 0.0922 0.403 0.0875 0.379 0.0508

0.406 0.0292 0.403 0.0157 0.403 0.0875 0.379 0.0508

1.440 0.6690 1.440 0.6690 1.440 0.6690 1.440 0.6690

0.000 0.3345 0.000 0.3345 0.000 0.3345 0.000 0.3345

0.7113 3.250

0.7113 3.250

0.7113 3.250

0.4410 3.400

0.3598 3.400

0.3598 3.400

0.0657 2.650

0.0657 2.650

0.0628 2.600

0.0628 2.600

0.6506 3.166

h124i h125i h134i

(kcal/mol) r (Å)

q (au) a (Å3) q (au) a (Å3) q (au) a (Å3) q (au)

NX and NY refer to the nitrogen atoms in the molecule h1xyi. Notice that Hm1 denotes two equivalent hydrogen atoms.

0.000 0.000

132

A. Møgelhøj et al. / Chemical Physics Letters 460 (2008) 129–136

Table 2 B3LYP basis set analysis for h125i and h134i Molecule

Atom

No. of basis funcs.

adz 183

atz

aqz

391

710

a5z

acdz

1162

207

actz

acqz

469

884

23.1 111.2 75.0 121.3

27.3 116.6 71.8 126.8

ac5z 1486

P1

P2

P3

167

212

17.0 105.3 78.7 114.2

18.2 104.1 77.5 114.1

324>

rN h125i h134i

h125i

N1 N2 N1 N2

N1

g

3.407 0.49 5.011 0.35 2.730 0.32 5.143 0.21

N1 N2 N1 N2

4.10 10.10 3.26 17.76

4.14 11.04 3.35 19.14

4.26 11.33 3.56 19.46

N2 N1

g h134i

h125i h134i

24.9 113.5 73.7 123.7

3.345 0.49 4.884 0.36 2.684 0.31 5.006 0.22

g h134i

17.3 103.3 79.2 113.1

3.335 0.46 4.452 0.43 2.684 0.30 4.537 0.28

g h125i

3.4 74.8 95.3 85.6

N2

28.4 118.0 71.0 128.3

5.2 86.9 88.5 97.7 NQCC/g

3.460 0.49 5.033 0.36 2.770 0.32 5.165 0.22

– – – –

3.432 0.47 4.620 0.42 2.749 0.31 4.717 0.27 dN

3.395 0.48 4.864 0.37 2.729 0.31 4.990 0.23

3.402 0.48 4.890 0.37 2.731 0.31 5.019 0.23

4.19 10.52 3.35 18.35

4.22 11.29 3.43 19.44

4.30 11.44 3.49 19.64

28.9 118.7 70.6 129.0

3.433 0.48 4.916 0.37 2.755 0.31 5.046 0.23

– – – –

19.0 105.3 77.8 114.9

3.578 0.47 4.900 0.42 2.897 0.30 5.022 0.28

3.521 0.48 4.905 0.40 2.836 0.32 5.041 0.26

3.404 0.48 4.908 0.37 2.759 0.30 5.037 0.24

4.82 12.01 2.67 20.09

4.32 11.21 3.43 19.21

4.38 11.04 3.47 18.89

The isotropic shielding constants, rN, are given in ppm and the 14N NQCCs are given in MHz. Included is also the asymmetry parameter, g. P1, P2 and P3 denote the Pople basis sets 6-311++G**, 6-311++G(2d,2p) and 6-311++G(3df,3pd), respectively. Also shown is the shift (in ppm) in the NMR shielding constant, dN, upon complexation with two water molecules.

smaller than the ones found for the absolute NMR shielding constants and the results are seen to be fairly converged using the atz basis set. The solvation shifts predicted using the Pople basis sets are in very good agreement with the ones obtained using the Dunning basis sets. Especially, the 6-311++G(2d,2p) basis set is found to provide an optimized performance/cost ratio. Therefore, we choose to use this basis set for the DFT/MM calculations of the NMR shielding constants. Note that this basis set performs almost as well as the very large acqz basis set. However, the number of basis functions is reduced by a factor of around 7. In contrast to what was observed for the NMR shielding constants, the NQCCs show a more clear convergence pattern with respect to the basis set. In fact, fairly small changes are observed upon increasing the basis set from acqz to ac5z, and the results based on the acqz basis set are fairly converged. It is noteworthy that inclusion of core basis functions contributes signiﬁcantly, which, based on the expression for the operator in Eq. (5), is to be expected. The Pople basis set 6-311++G(3df,3pd) reproduces the results for Vzz using the largest Dunning basis set to within a deviation of 0.8%. Due to the signiﬁcantly reduced number of basis functions in the 6-311++G(3df,3pd) basis set, we choose to use the 6-311++G(3df,3pd) basis set for all NQCC calculations. The size of the compounds studied in this Letters puts constraints on the level of ab initio theory that can be used in order to address the accuracy of the xc-functional used to predict the ef-

fect of solvation. In fact, it has recently been shown that different xc-functionals may predict very different solvation shifts in NMR shielding constants [40], and it is therefore of crucial importance to validate the accuracy of the different xc-functionals against high-level ab initio methods. The nitrogen atoms in the 1-methyltriazoles are either pyridine- or pyrrole-like. However, any molecule containing both types of nitrogens is too large in order to perform a reliable method analysis. Therefore, in order to address the accuracy of the different xc-functionals, and thereby the ability of DFT to predict accurately the speciﬁc properties in question, we follow an indirect approach. Based on microsolvated structures of NH3 and HCN, we calculate the shift due to complexation in the electric ﬁeld gradients and in the nitrogen NMR shielding constants using the xc-functionals B3LYP [22], KT3 [41] and PBE0 [42,43]. These results are then benchmarked against high-level coupled cluster predictions using CCSD and CCSD(T). In addition, we also explore the performance of the wave function methods HF and MP2. Since the nitrogens in HCN and NH3 are triple and single bounded, respectively, we expect that if a particular xc-functional is capable of producing reliable results for both HCN and NH3 compared to CCSD(T), it will also be able to perform well for the nitrogen atoms in the 1-methyltriazoles. Concerning the NQCCs, we here consider the solvent shifts in the EFG tensor components. Since a component of this tensor depends on the geometry and orientation of the molecule, we here keep the geometry of

Fig. 1. The orientation of the (a) HCN H2O and (b) NH3 H2O complexes.

133

A. Møgelhøj et al. / Chemical Physics Letters 460 (2008) 129–136

3.2. The combined (GIAO)-DFT/MM results

Table 3 Basis set analysis for NH3 and HCN using CCSD Molecule

adz

atz

aqz d

NH3 HCN NH3 HCN

Vxx Vzz Vxx Vzz

4.0 9.8

5.7 10.6

0.038 0.022 0.043 0.021

0.043 0.025 0.050 0.025

acdz

actz

acqz

4.4 10.0

5.8 10.9

– –

0.038 0.022 0.043 0.022

0.042 0.025 0.049 0.024

0.042 0.025 0.050 0.025

N

5.8 11.1 DV 0.042 0.025 0.051 0.025

The shift in the isotropic shielding constants, dN, is given in ppm and the shift in the EFG tensor components, DV, is given in au.

the isolated HCN or NH3 molecule as found in the complex. Since only one nitrogen atom is present in either HCN or NH3, we include only one water molecule in the complex. The complexes have been geometry-optimized at the CCSD(T)/aug-cc-pVTZ level using Dalton [34]. The geometries and the orientations of the complexes are shown in Fig. 1. Table 3 contains a basis set analysis of the shift due to complexation in either the nitrogen NMR shielding constant or the EFG tensor. This basis set analysis is carried out using CCSD. Since basis set requirements are more severe for methods like CCSD than DFT (correlated wave function methods require a proper description of the virtual orbitals in order to describe correlation effects accurately), the DFT results are also converged upon using the same basis set as required by (converged) CCSD. From Table 3 we ﬁnd, that the atz basis set is sufﬁcient for both properties. Therefore, this basis set is chosen for the xc-functional analysis. In Fig. 2 we have illustrated the shifts due to complexation with respect to the CCSD(T) predictions for the nitrogen NMR shielding constant and the Vxx component of the EFG tensor. For the shifts in the NMR shielding constant, the KT3 xc-functional is seen to give slightly improved results compared to B3LYP or PBE0. However, for the Vxx component of the EFG tensor, all DFT based methods are found to slightly overestimate the shift and none of the three xcfunctionals considered can be said to perform better. However, the errors introduced in the DFT based calculations are relatively small as compared to the CCSD(T) predictions. Therefore, we choose the KT3 xc-functional to use in the DFT/MM calculations of both the NMR shieldings and NQCCs. We note that we have not included effects of basis set superposition errors (BSSE) in these calculations. However, in the case of DFT(B3LYP) we found this to be rather small (around 0.5 ppm for the NMR shielding constant and negligible for the EFG tensor components) and have thereby neglected this in the presented results.

a

Having determined a proper basis set and xc-functional, we now turn to the large-scale DFT/MM calculations. An important question to address in this respect is the number of solute–solvent conﬁgurations needed to consider in the DFT/MM calculations in order to provide statistically converged results. Fig. 3 contains an illustration of the convergence proﬁle for a nitrogen (N2) NMR shielding constant and a component of the EFG tensor in h125i. As seen from Fig. 3, both the shielding constant and EFG tensor component converge very quickly, and 120 conﬁgurations are sufﬁcient. The rest of the DFT/MM calculations are thereby based on this number of conﬁgurations. In the simplest DFT/MM scheme, the solute is only perturbed electrostatically by the MM solvent molecules. However, in many situations special attention has to be payed to the hydrogen-bonding water molecules [36,37]. A simple solution to this problem is to include explicitly the water molecules with a direct coordination to the pyridine-type nitrogen atoms of the solute in the QM region of the model. Inspection of the radial distribution functions reveals that one pyridine-type nitrogen on average coordinates approximately one water molecule. This number is to be considered as an upper limit to the number of hydrogen-bonding water molecules. In the DFT/MM calculations presented below, we therefore include in the region treated using DFT both the solute and two water molecules, while the rest of the solvent water molecules are described classically. The effect of including further water molecules into the QM region is very small. The inclusion of four water molecules – the two closest to each pyridine-type nitrogen atom – leads to the change in the solvent shift of the NMR shielding constant of h125i by less than 1 ppm, whereas the results for the NQCCs are, within the statistical errors, unaffected. The results presented below thereby refer to the case where two water molecules have been included into the region treated using DFT. For a few selected conﬁgurations we have tested the effect of BSSE. As in the case of the smaller clusters discussed previously in this Letter, BSSE are found to be very small (in fact around the statistical uncertainty), and we have thereby neglected the effect of BSSE in the ﬁnal DFT/MM calculations. The DFT/MM based results for the solvent shift in the isotropic NMR shielding constants of the 1-methyltriazoles are given in Table 4. In this Table we have also included experimental shifts. Note, however, that these experimental shifts are based on measurements in cyclohexane and water solutions. For all the 1-methyltriazoles we observe a qualitatively different effect of solvation on the pyrrole- and pyridine-type nitrogens. This reﬂects a differential increase or decrease, respectively, in the electronic charge

b

4

0.015

NH3 HCN

3

NH3 HCN

0.01 Δ Vxx (au)

δ(ppm)

2 1 0

0.005

0

-1 -0.005 -2 -3

-0.01 HF

B3LYP PBE0

KT3

MP2 CCSD

HF B3LYP PBE0 KT3

MP2 CCSD

Fig. 2. The shifts in the (a) nitrogen NMR shielding constant and (b) nitrogen Vxx component of the electric ﬁeld gradient tensor due to the water molecule in the HCN H2O and NH3 H2O complexes with respect to the corresponding CCSD(T) predictions. The basis set used is aug-cc-pVTZ.

134

A. Møgelhøj et al. / Chemical Physics Letters 460 (2008) 129–136

a -54

b

-0.33

-0.34

-56

Vzz (au)

N shielding constant (ppm)

-0.335 -55

-57

-0.345 -0.35 -0.355 -0.36

-58

-0.365 -59

-0.37 0

20

40

60

80

100

120

0

20

No. of configurations

40

60

80

100

120

No. of configurations

Fig. 3. The convergence proﬁles of the (a) N2 NMR shielding constant and (b) the N2 Vzz component of the electric ﬁeld gradient tensor of h125i with respect to the number of molecular conﬁgurations included in the statistical averaging.

Table 4 The solvent shifts in the nitrogen NMR isotropic shielding constants (in ppm) for the 1-methyltriazoles based on the DFT/MM calculations using the KT3 xc-functional and 6-311++G(2d,2p) basis Molecule

Model

Atom

dN teo

dN exp

h123i

QM/MM

h124i

QM/MM

h125i

QM/MM

N1 N2 N3 N1 N2 N4 N1 N2 N1 N2 N1 N3

10.7 ± 0.3 13.6 ± 0.6 17.5 ± 0.6 7.8 ± 0.3 10.2 ± 0.7 14.8 ± 0.6 2.2 ± 0.3 7.5 ± 0.7 1.5 ± 0.1 6.5 ± 0.4 15.2 ± 0.3 25.7 ± 0.7

7.79 19.98 22.05 4.79 13.11 15.81 0.19 10.88 0.19 10.88 11.67 30.71

PCM h134i

QM/MM

The two water molecules closest to the pyridine-type nitrogen atoms have been included in the QM region. The experimental results are obtained from Ref. [49] and refer to the difference between the measurements in water and cyclohexane solution.

distribution around the nitrogen atoms upon solvation. The experimental results for the solvent shifts for the different nitrogen atoms within a given compound are rather scattered. For instance, for h123i a relative solvent shift evaluated by comparing the solvent shifts for the individual nitrogen atoms within a given compound of 30 ppm is found. Our calculations reproduce this spread very well. In the case of h123i, we predict a relative solvent shift of 28 ppm, in very good agreement with the experimental data. Such a good agreement is found for all the 1-methyltriazoles and represents a signiﬁcant improvement as compared to a previous study using a dielectric continuum model based on spherical cavity for the description of solvent effects [44]. Comparing, however, the experimental and calculated solvent shifts for each nuclei directly, we ﬁnd that the magnitude of the predicted solvent shift is systematically overestimated for the pyrrole-type nitrogen atoms, while it is underestimated for the pyridine-type nitrogens. However, since the relative effect of solvation is well reproduced, this systematic error is most likely due to missing effects in our model, most notably solvent effects on the zero-point vibrational contribution to the nitrogen NMR shielding constants. In fact, it has recently been found that this effect might be substantial [45]. In addition, the experimental solvent shifts are evaluated with respect to cyclohexane solution while the predicted ones are calculated with respect to vacuum. Note that for the nitrogen atoms within each compound there is a fairly good correlation between the atomic charge

located on the nitrogen atoms used in the MD simulations and the solvation shift in NMR shieldings, i.e. the more negative charges exhibit larger (positive) shifts. This is due to formation of stronger hydrogen bonds when the partial atomic charge becomes more negative. Due to steric effects, this correlation is not seen across the individual molecules. For one of the compounds, h125i, we have in addition to the DFT/MM model also used the PCM model in order to study solvent effects on the nitrogen NMR shielding constants. For the PCM calculations we use the same solute–solvent structures as in the DFT/ MM predictions, but replace all MM molecules by the dielectric continuum. The radii used for the formation of the cavity are the same as used in the previously discussed PCM calculations. These calculations have been performed using the PCM model implemented in Dalton [46–48]. Note that it is only the bulk solvent that is described using the PCM, while the hydrogen-bonding water molecules are considered explicitly using DFT. From Table 4 we obTable 5 The 14N NQCCs (in MHz) for the 1-methyltriazoles based on the DFT/MM calculations using the KT3 xc-functional and 6-311++G(3df,3pd) basis Molecule

Model

Atom

Vacuum

Exp.

Solution

Exp.

h123i

QM/MM

N1 N2 N3 DmN1 2 =DmN1 1

2.790 4.246 4.705 2.74

– – – 2.47

2.426 ± 0.002 3.926 ± 0.002 4.312 ± 0.003 3.10

– – – 4.98

2.76

2.33

2.92

4.42

2

2

DmN1 3 =DmN1 1 2

2

N2

N3

2

2

Dm1 =Dm1 h124i

QM/MM

N1 N2 N4 DmN1 2 =DmN1 1 2

2

2

N2

N3

2

2

Dm1 =Dm1 QM/MM

PCM

h134i

QM/MM

N1 N2 DmN1 2 =DmN1 1 2 2 N1 N2 DmN1 2 =DmN1 1 2 2 N1 N3 DmN1 3 =DmN1 1 2

1.06

1.06

1.13

3.137 4.471 4.033 2.28

– – – 2.04

3.007 ± 0.002 4.100 ± 0.003 3.532 ± 0.003 2.13

– – – 3.94

1.55

1.61

1.25

2.50

2

DmN1 3 =DmN1 1

h125i

0.99

1.47

1.27

1.70

1.58

3.472 4.604 1.77 3.472 4.604 1.77 2.789 4.676 2.79

– – 1.74 – – 1.74 – – 2.73

3.401 ± 0.002 4.417 ± 0.003 1.68 3.423 ± 0.002 4.428 ± 0.003 1.68 2.315 ± 0.002 4.210 ± 0.003 3.57

– – 2.84 – – 2.84 – – 5.97

2

The two water molecules closest to the pyridine-type nitrogen atoms have been included in the QM region. Relative half-height widths have also been included. The experimental results are obtained from Ref. [44]. The experimental results for vacuum refer to cyclohexane solution.

A. Møgelhøj et al. / Chemical Physics Letters 460 (2008) 129–136

serve that the PCM tends to underestimate slightly the solvent shifts as compared to DFT/MM. However, the discrepancies are small. Thereby, the effect of the bulk solvation may be accurately described using the PCM, which is in line with the ﬁndings in Refs. [36,37]. The results for the solvent effects on the NQCCs are given in Table 5. Here we observe for all 1-methyltriazoles a decrease in the magnitude of the NQCCs upon solvation. This is in agreement with a previous study [44] although the model employed in Ref. [44] in fact neglected the explicit solvent contribution discussed in Section 2.1. In contrast to the solvent shift for the NMR shielding constants, no general picture emerges when comparing the shifts for the pyrrole- or pyridine-type nitrogen atoms. This indicates that solvent effects on the NQCCs are not as dependent on hydrogen bonds as the corresponding NMR shielding constants. Since the protocol used in this work is unable to provide the actual correlation time, only relative half-height widths may be compared to experiment. For the isolated molecules we obtain a relatively good agreement between the calculated and experimentally measured relative half-height widths. In fact, the agreement is signiﬁcantly better than found in a previous study [44] based on a RASSCF wave function description. The good agreement between measured and predicted relative half-height widths indicates that the correlation time for the pyrrole- or pyridine-type nitrogen atoms in vacuum can be assumed to be fairly similar. Moving to the aqueous solution, however, the agreement between the measured and predicted relative half-height widths becomes less satisfactory. In fact, in some cases the measured and predicted relative half-height widths (with respect to the pyrrole-type nitrogen atom) behave qualitatively different upon solvation. Most likely, this is due to the assumption of equal correlation times for the pyrrole- or pyridine-type nitrogen atoms. As discussed above, these two atomtypes behave very differently when exposed to a polar protic solvent, and it is therefore expected that the correlation times will be different in an aqueous solution. We note that the agreement between the experimental and predicted relative half-height widths may generally be improved by scaling the correlation times according to spyridine-type 1:6 spyrrole-type , although this scaling c c does not consistently yield better results. For the h123i and h124 i compounds the aspect related to different correlation times may be further explored by considering the relative half-height widths for only the pyridine-type nitrogen atoms, e.g. the N2 width relatively to the N3 width. As seen from Table 5, the agreement between the predicted and experimental data now becomes much better and in fact fairly satisfactory. We thereby conclude that the apparent disagreement between theory and experiment in this case is related to different correlation times for the pyrrole- and pyridine-type nitrogen atoms and that valuable interpretation of the data obtained in aqueous solution can only be made by considering relative widths between the same atom types. We ﬁnally notice the very similar performance between the DFT/MM and PCM models. This is expected since, as discussed above, the NQCCs are not as sensitive to effects of hydrogen-bonding – and thereby directional electrostatic interactions – as the NMR shielding constants.

4. Summary In this Letter we have explored solvent effects on the nitrogen NMR shielding and nuclear quadrupole coupling constants in 1methyltriazoles. We have presented a detailed investigation of the effects of both basis sets and xc-functionals. The majority of the calculations have been performed using a discrete polarizable solvation model based on solute–solvent structures extracted from MD simulations. Fairly good results are found for the solvent shifts

135

of the nitrogen NMR shielding constants thereby validating the solvent model. For the nuclear quadrupole coupling constants a less clear picture emerges. However, based on a detailed analysis, we argue that the discrepancies between measured and predicted quantities are mainly due to the assumption of similar correlation times for the pyrrole- and pyridine-type nitrogen atoms. In fact, a fairly good agreement between theory and experiment is recovered by considering only relative widths between the same type of nitrogen atoms. Finally, we have for a selected compound considered the possibility to include solvent effects using a dielectric continuum model. We ﬁnd that the continuum model, when combined with a statistical description of the solvent and applied to the solute including the nearest solvent molecules explicitly, describes the effect of bulk solvation at a level comparable to a discrete environment. Acknowledgements The authors thank Prof. Michal Jaszunski for valuable comments on this manuscript. We thank Prof. Kenneth Ruud for providing us with the PCM implementation in Dalton. The authors thank the Danish Center for Scientiﬁc Computing (DCSC) for computational resources. A.M. thanks the Novo Scholarship program for ﬁnancial support. K.A. and K.V.M acknowledge support from Forskningsrådet for Natur og Univers through Copenhagen Center for Atmospheric Research. J.K. thanks the Villum Kann Rasmussen foundation for ﬁnancial support. References [1] D.M. Grant, R.K. Harris (Eds.), Encyclopedia of Nuclear Magnetic Resonance, Wiley, New York, 1996. [2] T. Helgaker, M. Jaszunski, K. Ruud, Chem. Rev. 99 (1999) 293. [3] J. Tomasi, B. Mennucci, R. Cammi, Chem. Rev. 105 (2005) 2999. [4] A. Warshel, M. Levitt, J. Mol. Biol. 103 (1976) 227. [5] U.C. Singh, P.A. Kollman, J. Comput. Chem. 7 (1986) 718. [6] M.J. Field, P.A. Bash, M. Karplus, J. Comput. Chem. 11 (1990) 700. [7] M.A. Thompson, J. Phys. Chem. 100 (1996) 14492. [8] E. Fadda, M.E. Casida, D.R. Salahub, J. Phys. Chem. A 107 (2003) 9924. [9] J. Gauss, J. Chem. Phys. 99 (1993) 3629. [10] J. Olsen, P. Jørgensen, J. Chem. Phys. 82 (1985) 3235. [11] K. Ruud, T. Helgaker, R. Kobayashi, P. Jørgensen, K.L. Bak, H.J.Aa. Jensen, J. Chem. Phys. 100 (1994) 8178. [12] C.B. Nielsen, O. Christiansen, K.V. Mikkelsen, J. Kongsted, J. Chem. Phys. 126 (2007) 154112. [13] J. Kongsted, C.B. Nielsen, K.V. Mikkelsen, O. Christiansen, K. Ruud, J. Chem. Phys 126 (2007) 034510. [14] F. London, J. Phys. Radium 8 (1937) 397. [15] H.F. Hameka, Rev. Mod. Phys. 34 (1962) 87. [16] R. Ditchﬁeld, Mol. Phys. 27 (1974) 789. [17] K. Wolinski, J.F. Hinton, P. Pulay, J. Am. Chem. Soc. 112 (1990) 8251. [18] A. Abragam, The Principles of Nuclear Magnetic Resonance, Oxford University Press, Oxford, 1961. [19] A.J. Stone, Theory of intermolecular forces, Oxford University Press Inc., New York, 1996. [20] A. Osted, J. Kongsted, K.V. Mikkelsen, O. Christiansen, Mol. Phys. 101 (2003) 2055. [21] P. Linse, MOLSIM is an integrated MD/MC/BD simulation program belonging to the MOLSIM package, Version 3.3.0, December 05, 2001. [22] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [23] R.A. Kendall, T.H. Dunning, R.J. Harrison, J. Chem. Phys. 96 (1992) 6796. [24] S. Miertuš, E. Scrocco, J. Tomasi, Chem. Phys. 55 (1981) 117. [25] E. Cancés, B. Mennucci, J. Math. Chem. 23 (1998) 309. [26] E. Cancés, B. Mennucci, J. Tomasi, J. Chem. Phys. 107 (1997) 3031. [27] B. Mennucci, E. Cancés, J. Tomasi, J. Phys. Chem. B 101 (1997) 10506. [28] M.J. Frisch et al., GAUSSIAN 03, Revision B.05, Gaussian, Inc., Wallingford, CT, 2004. [29] C.M. Breneman, K.B. Wiberg, J. Comp. Chem. 11 (1990) 361. [30] L. Gagliardi, R. Lindh, G. Karlström, J. Chem. Phys. 121 (2004) 4494. [31] G. Karlström et al., Comput. Mater. Science 28 (2003) 222. [32] J. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, J. Comput. Chem. 25 (2004) 1157. [33] P. Ahlström, A. Wallqvist, S. Engström, B. Jönsson, Mol. Phys. 68 (1989) 563. [34] DALTON, a molecular electronic structure program, release 2.0, 2005.

136

A. Møgelhøj et al. / Chemical Physics Letters 460 (2008) 129–136

[35] O. Christiansen, M.B. Hansen, J. Kongsted, P. Seidler, D. Toffoli, MidasCpp, Molecular Interactions, Dynamics and Simulation Chemistry program package in C++, 2006. [36] K. Aidas, et al., J. Phys. Chem. A 111 (2007) 4199. [37] J. Kongsted, B. Mennucci, J. Phys. Chem. A 111 (2007) 9890. [38] D.E. Woon, T.H. Dunning, J. Chem. Phys. 103 (1995) 4572. [39] M.J. Frisch, J.A. Pople, J.S. Binkley, J. Chem. Phys. 80 (1984) 3265. [40] J. Kongsted, K. Aidas, K.V. Mikkelsen, S.P.A. Sauer, J. Chem. Theory Comput. 4 (2008) 267. [41] T.W. Keal, D.J. Tozer, J. Chem. Phys. 121 (2004) 5654. [42] M. Ernzerhof, G.E. Scuseria, J. Chem. Phys. 110 (1999) 5029.

[43] C. Adamo, V. Barone, J. Chem. Phys. 110 (1999) 6158. [44] M. Jaszunski, K.V. Mikkelsen, A. Rizzo, M. Witanowski, J. Phys. Chem. A 104 (2000) 1466. [45] J. Kongsted, K. Ruud, Chem. Phys. Lett. 451 (2008) 226. [46] R. Cammi, L. Frediani, B. Mennucci, K. Ruud, J. Chem. Phys. 119 (2003) 5818. [47] M. Pecul, D. Marchesan, K. Ruud, S. Coriani, J. Chem. Phys. 122 (2005) 024106. [48] L. Ferrighi, D. Marchesan, K. Ruud, L. Frediani, S. Coriani, J. Chem. Phys. 123 (2005) 204104. [49] M. Witanowski, W. Sicinska, Z. Biedrzycka, Z. Grabowski, G.A. Webb, J. Magn. Res. A 112 (1995) 66.