Article pubs.acs.org/JCTC
Nuclear Magnetic Shielding Constants from Quantum Mechanical/ Molecular Mechanical Calculations Using Polarizable Embedding: Role of the Embedding Potential Casper Steinmann,* Jógvan Magnus Haugaard Olsen, and Jacob Kongsted Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, 5230 Odense, Denmark S Supporting Information *
ABSTRACT: We present NMR shielding constants obtained through quantum mechanical/molecular mechanical (QM/ MM) embedding calculations. Contrary to previous reports, we show that a relatively small QM region is sufficient, provided that a high-quality embedding potential is used. The calculated averaged NMR shielding constants of both acrolein and acetone solvated in water are based on a number of snapshots extracted from classical molecular dynamics simulations. We focus on the carbonyl chromophore in both molecules, which shows large solvation effects, and we study the convergence of shielding constants with respect to total system size and size of the QM region. By using a high-quality embedding potential over standard point charge potentials, we show that the QM region can be made at least 2 Å smaller without any loss of quality, which makes calculations on ensembles tractable by conventional density functional theory calculations.
1. INTRODUCTION The broad range of applications to which NMR spectroscopy can be applied to solve scientific problems were recently highlighted in a paper by Barrett et al.1 Recent years have witnessed an increased interest in using calculated NMR data in combination with experiments in order to elucidate, for example, structural properties of complex materials such as proteins.2 In this same period, methods that can be used to predict NMR properties have also greatly increased their level of sophistication. Indeed, it is today possible to use ab initio methods to calculate very accurate values for key NMR parameters for smaller molecular systems. For larger systems and usually of biological interestit has become popular to use methods relying on some sort of parametrization of the NMR parameters directly in terms of structural data; that is, the NMR parameters are predicted on the basis of simple analytical and parametrized relationships, giving as input the molecular structure. Even though such parametrized schemes can be tuned to become rather accurate, predictions based directly on ab initio methods are surely to be preferred. One class of large molecular systems, which are often met in chemistry, is that of a solvated molecule. For solute−solvent systems, which is the subject of this paper, continuum approaches have often been used in order to provide a representation of the solvent.3 However, explicit treatment of the solvent, including dynamical effects from classical molecular dynamics (MD) simulations, have gained popularity. It is well-known that it is often sufficient to use a cluster model with the solute and sometimes also its nearest-neighbor solvent molecules treated by high-level quantum mechanics (QM), while the remainder of the solvent environment is treated classically, often using simple point © 2014 American Chemical Society
charge potentials taken from standard molecular mechanics (MM) force fields. It has, however, been shown that use of ab initio-derived polarizable potentials that can adapt to the quantum mechanically treated charge distribution is preferable in order to obtain good accuracy.4,5 In the present paper we demonstrate that, by using such high-quality embedding potential, one can obtain converged nuclear magnetic shielding constants (and chemical shifts) for even quite small QM regions compared to the case where the environment is modeled by a standard point charge model. We show that an important aspect of embedding calculations is the selfconsistent account of polarization effects such as in the recently developed polarizable embedding approach by Olsen et al.4,6,7 Even though it has previously been found that large QM regions were needed to obtain converged NMR shielding constants,8 we here show that the QM size can be reduced by using a better embedding potential. Thus, we show that the polarizable embedding strategy represents a computationally efficient route if the aim is to calculate key NMR properties of only a small part of a larger molecular system. In the present paper our focus is on solute−solvent systems, but the aim is to extend the formalism presented in this work to polymeric systems such as proteins. This paper presents the underlying theory behind calculations of gauge-origin-independent NMR shielding constants within the polarizable embedding formalism. Received: October 9, 2013 Published: February 14, 2014 981
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
2. METHODOLOGY 2.1. NMR Shielding Tensors in the Polarizable Embedding Framework. The calculation of an element of the NMR shielding tensor for nucleus N, σNij , in self-consistent field theory is calculated according to9 σijN = 1 +
∑ Dμν μν
∂ 2hμν ∂Bi ∂mj , N
+
∑ μν
moment expansions and the induction term on distributed polarizabilites. In the following we will use a multi-index notation that allows more compact expressions. We use a three-dimensional multi-index defined as k = (kx, ky, kz), whose norm is given by the sum of its elements, |k| = kx + ky + kz, and factorial is the product between the factorial of the elements, k! = kx!ky!kz!. The multi-index power of a vector in Cartesian space is defined as rk = xkxykyzkz, which allows us to write a general partial derivative operator as ∇kr = (∂|k|)/(∂xkx ∂yky ∂zkz). By use of this notation, the electrostatic part of the PE operator is given as4
∂Dμν ∂hμν ∂Bi ∂mj , N
(1)
Here, Dμν is an element of the density matrix in the atomic orbital (AO) basis, hμν is an element of the one-electron Hamiltonian also in AO basis, Bi is a component of the magnetic field (induction), and mj,N is a component of the nuclear magnetic moment of nucleus N. Evaluation of the last term in eq 1, that is, the paramagnetic contribution, requires knowledge of the perturbed density matrix, which in our case is obtained by solving a set of appropriate response equations10 for each of the three components of the magnetic field. In order to remove the dependency on the gauge origin in the calculated chemical shielding tensors, gauge-including atomic orbitals (GIAOs)10−12 are employed; that is, the AO basis functions have attached a phase that carries an explicit reference to the perturbation: χμN (B)
⎤ ⎡ i = exp⎢ − (B × R μ) ·r⎥ χμN (0) ⎦ ⎣ 2
K
veŝ =
|k|= 0 s = 1
KS
ts(,kpq) = − p ∇kr
M
⎤
m=1
⎦
∑ Tsm(k)Zm⎥⎥
(5)
1 q |r − R s|
(6)
and T(k) sm is an element of the interaction tensor defined as (k) Tsm = ∇km
(2)
1 |R m − R s|
(7)
For each value of |k| in the summation in eq 5, there is an implicit summation over (|k| + 1)(|k| + 2)/2 multi-indices; for example, for |k| = 1 it is k = (1, 0, 0), k = (0, 1, 0), and k = (0, 0, 1). The summation runs up to K, which is the highest order of multipoles included in the embedding potential. We use a second quantization formulation where Ê pq is a singlet excitation operator expressed as a string of creation a†pσ and annihilation operators.14 The v̂ind operator is given as S
vind ̂ = −∑
∑ μsind,(k) ∑ ts(,kpq) Epq̂
|k|= 1 s = 1
(8)
pq
where the summation over multi-indices is restricted to |k| = 1 and μind,(k) is an element of an induced dipole at site s. Formally, s the induced dipoles are determined from a set of coupled linear equations that can be expressed in matrix form as
μind = AF
(9) 15
where A is the classical response matrix (also known as the relay matrix) and F is the combined field from the QM electrons and QM nuclei and the field from the permanent multipoles in the environment: F = Fe + Fn + Fmul
(10)
Here, μ and F are supervectors of length 3S, where S is the number of polarizable sites. Note that the number of expansion points for the multipoles and the polarizabilities need not be equal. The induced dipoles can be solved for by use of eq 9. This can be done directly, which requires explicit construction of the classical response matrix. However, due to the generally large size of the matrix (3S × 3S), eq 9 is usually solved via iterative techniques. Due to the use of GIAOs in combination with the embedding potential, specific contributions to the response equations appear. These contributions are here presented as a ind
KS
(3)
̂ where fKS vac is the gas-phase Kohn−Sham operator and v̂PE is the polarizable embedding (PE) operator given as vPE ̂ = veŝ + vind ̂
⎡ ( −1)|k| (k)⎢ ̂ + Ms ∑ ts(,kpq) Epq ⎢⎣ k! pq
Here, M(k) is the kth element of the |k|th order multipole s moment at a site s; that is, for k = (0, 0, 0) it is a charge, for k = (1, 0, 0) it is the x component of a dipole, and so on. t(k) s,pq are integrals over the |k|th order derivative of |r − Rs|−1:
Here, Rμ is the vector between the global gauge origin and the position of the field-independent AO, χNμ (0), attached to nucleus N. As shown and discussed in detail previously,13 contributions from the embedding potential enter into eq 1 in two places: in the first term (the diamagnetic contribution) through the density matrix and in the second term (the paramagnetic contribution) through the perturbed density matrix. The polarizable embedding (PE) scheme represents a model for explicit inclusion of environmental effects in quantum chemical calculations.4,6,7 Within the embedding scheme, the environment is described by use of a multicenter multipole expansion of the electric potential created by the environment. Usually the expansion centers are taken at the position of the atomic nuclei, but other choices might also be relevant. Polarization of the environment is accounted for through the use of distributed polarizabilities that give rise to induced dipoles. The PE method may be coupled to different quantum chemical formulations, and the most used of these is currently PE in combination with density functional theory (DFT), referred to as PE-DFT.4 This method is available for both energies and ground-state molecular properties as well as excited-states properties through the use of time-dependent DFT (PE-TD-DFT). In the PE-DFT scheme, an effective ̂ Kohn−Sham (KS) operator fKS eff is given as
̂ + v̂ feff̂ = fvac PE
S
∑∑
(4)
Here, v̂es and v̂ind are the one-electron operators for electrostatic interactions between environment and QM core and for induction effects (polarization) of the environment, respectively. The electrostatic term is based on distributed multipole 982
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
multipole moments and polarizabilities are distributed onto the atomic sites of the molecules that define the environment. To evaluate the performance of the embedding scheme in calculation of shielding constants determined under the influence of either PE or TIP3P potentials, two different systems based on earlier studies are investigated. For this we have used 120 evenly distributed snapshots of acetone and acrolein (Figure 1) solvated in water, obtained from a MD
generalization of the terms described in ref 13 to reflect the current implementation. In short, in the field-dependent AO basis, the elements of the v̂PE operator take the form ⟨χμ (B)|vPE ̂ |χν (B)⟩ = ⟨χμ (B)|veŝ |χν (B)⟩ + ⟨χμ (B)|vind ̂ |χν (B)⟩ K
=
S
|k|= 0 s = 1 S
−
|k |
∑ ∑ (−1) k!
Ms(k)ts(,kμν) (B)
∑ ∑ μsind,(k)ts(,kμν) (B) |k|= 1 s = 1
(11)
where the integrals now carry an explicit dependence on the magnetic field due to the use of GIAOs. The contributions needed in order to determine the first-order density matrix (eq 1) require evaluation of the first-order derivative of eq 11 with respect to a component of the magnetic field: Figure 1. Depictions of acrolein (left) and acetone (right) with atomic labels.
d⟨χμ (B)|vPE ̂ |χν (B)⟩ dBi K
=
S
|k |
∑ ∑ ( − 1) k!
|k|= 0 s = 1 S
−
∑ ∑ μsind,(k) |k|= 1 s = 1
Ms(k)
dts(,kμν) (B)
simulation described in works by Aidas et al.29,30 The 120 snapshots were used in order to perform statistical averaging. We will analyze the convergence of the calculated shielding constants with respect to size of the QM region and size of the cluster, that is, the whole system (QM + environment), as it approaches bulk. The size of the QM region and total system are defined as minimum atom-based distances from the solute to the solvent molecules. The size of the QM region is determined in such a way that if a solvent molecule has at least one atom within a distance RQM from the solute, it is included in the QM region. Otherwise it is treated classically. Similarly, the size of the molecular cluster is determined in such a way that if an atom of a solvent molecule is within Rsys it is included; otherwise it is removed. The size of the whole system is always greater than or equal to the size of the QM region (Rsys ≥ RQM). Chemical shifts were obtained by subtracting the shielding constants in the isolated molecule from those in the solvated molecule:
dBi
dts(,kμν) (B) dBi
(12)
where the integral derivatives are defined according to dts(,kμν) (B) dBi
dχμ (B)
=
dBi +
∇kr
χμ (B) ∇kr
1 χ (B ) |r − R s| ν dχν (B) 1 dBi |r − R s|
(13)
These integrals are evaluated analytically for zeroth- and firstorder derivatives of |r − Rs|−1.13 The integrals corresponding to higher-order derivatives are evaluated on the basis of an implemented finite difference scheme using the analytically calculated integrals. No problems with stability were observed in the evaluation of second-order derivative integrals. In a future implementation, we expect to perform the integral evaluation completely analytically to any order by use of the Gen1Int16 package. 2.2. Computational Details. The GIAO contributions to the NMR shielding constants calculated via PE-DFT were added to a development version of the Dalton17,18 program using the Polarizable Embedding library7 and will be available in a future release. Geometry optimizations of vacuum structures were performed with Gaussian 0919 using the B3LYP20,21 exchange−correlation (xc) functional and the augcc-pVTZ22 basis set. Shielding constants were calculated with Dalton18 using the KT323 GGA xc-functional and the pcS-1 and pcS-2 polarization-consistent basis sets by Jensen,24 optimized for calculating nuclear magnetic shielding constants. For the classical water potentials, either TIP3P25,26 or LoProp27 potentials were used. The LoProp potentials were calculated with Molcas28 at the B3LYP/aug-cc-pVDZ level of theory for all multipole moments through quadrupoles and dipole polarizabilities (this potential is denoted PE in the following). Note that we used an ANO-type recontraction of the aug-ccpVDZ basis set as required by the LoProp procedure. All
δ = σsol − σvac
(14)
The presented error bars and uncertainties are calculated as standard error of the mean (SEM) of the shielding constants for a set of snapshots as s SEMN = (15) 120 where s is the standard deviation of the calculated NMR shielding constants for a nucleus N and 120 is the number of snapshots used to calculate the shielding constants. Optimal bin widths for histograms are calculated according to the Freedman−Diaconis rule.31
3. RESULTS We shall initially be concerned with the convergence of the shielding constants as a function of size of the system (Rsys) and size of the QM region (RQM). For this we have chosen acrolein, and we will investigate in detail the convergence of carbonyl chromophore 17O and 13C shielding constants. For all 120 MD snapshots, we calculated the shielding constants for Rsys ∈ [1, 12] Å and RQM ∈ [1, 4] Å. We used the pcS-1 basis set for this part of the study because the need for large (RQM = 4 Å) QM 983
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
calculations would otherwise be computationally too demanding, given the large number of calculations that is required. As we will discuss later, the choice of pcS-1 is satisfactory for evaluation of the shielding constants from QM/MM calculations but is not enough for prediction of accurate gas-phase shielding constants. 3.1. 17O and 13C Nuclear Shielding Constants in Acrolein. Convergence of the 17O shielding constant in acrolein, calculated at the PE-KT3/pcS-1 level of theory, with respect to system size is presented in Figure 2 for several values
Figure 3. Convergence of absolute 17O NMR shielding constant in acrolein, calculated with the TIP3P potential as a function of system size and size of QM region. The results have been averaged over 120 MD snapshots. The first value of each series is a pure QM calculation where RQM = Rsys. Calculations were performed at the KT3/pcS-1 level of theory.
Figure 2. Convergence of absolute 17O NMR shielding constant in acrolein, calculated with the PE potential as a function of system size and size of QM region. Results have been averaged over 120 MD snapshots. The first value of each series is a pure QM calculation where RQM = Rsys. Calculations were performed at the KT3/pcS-1 level of theory.
of RQM. The shielding constant increases very rapidly as the size of the entire complex is increased but slows down and eventually converges at around Rsys = 8 Å. Increasing Rsys to 12 Å yields little improvement and the converged values are here −221.0, −226.3, −223.1, and −223.2 ppm for RQM = 1, 2, 3, and 4 Å, respectively. The computational cost of using Rsys = 12 Å over Rsys = 8 Å is negligible, and we have therefore used this size in the following calculations. In Figure 3 the corresponding results are shown for the TIP3P potential. In this case the converged values at Rsys = 12 Å are −256.4, −244.5, −232.1, and −228.1 ppm for RQM = 1, 2, 3, and 4 Å, respectively. We note that there is no observable difference in the convergence behavior with respect to system size for both PE and TIP3P; however, the effect of size of the QM region varies greatly between the two, with the PE results showing a much faster convergence with respect to size of the QM region. The shielding constants of 17O based on the use of the PE and TIP3P potentials as a function of size of the QM region (system size is fixed at Rsys = 12 Å) are shown in Figure 4. At RQM = 1 Å the calculated PE- and TIP3P-based average shielding constants deviate considerably from the expected convergence pattern. The reason is most likely that we are missing important quantum mechanical effects, since at RQM = 1 Å it is only acrolein that is modeled with quantum mechanics while all solvent molecules are treated classically and we therefore completely neglect exchange−repulsion with the
Figure 4. Convergence of absolute 17O shielding constants in acrolein as a function of size of the QM region. The system size is fixed at Rsys = 12 Å. Results have been averaged over 120 MD snapshots. Calculations were performed at the KT3/pcS-1 level of theory.
solvent. At RQM = 2 Å, however, we include hydrogen-bonded solvent molecules in the QM region and thereby also recover important quantum-mechanical interactions with the solvent. This is in line with the study by Fradelos and Wesolowski32 where they demonstrated that the lack of exchange−repulsion can introduce large errors in pure Coulombic embedding models especially when large or diffuse basis sets are used. The different behaviors of PE and TIP3P potentials can be explained by the fact that the PE potential offers a very accurate electrostatic potential compared to TIP3P.4,5 The same errors, due to the lack of important short-range quantum mechanical effects, are also present in the calculations of the shielding constants of 13C in the acrolein carbonyl chromophore (see Figure 5), and here we also observe more 984
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
for small QM regions are also valid for acetone. On this basis, we do not consider a detailed analysis of the RQM = 1 Å results. Note that the results below are all obtained for Rsys = 11 Å. Rather than evaluating the shielding constants for acetone, we find it more instructive to calculate the gas-to-aqueous solution chemical shifts of the carbonyl chromophore because it introduces additional requirements on the quality of the calculation, as we shall see below. 3.2.1. Chemical Shift of Carbonyl 17O in Acetone. The calculated gas-to-aqueous solution chemical shift of 17O in acetone is obtained according to eq 14 and is listed in Table 1 Table 1. Calculated Acetone 17O Chemical Shifts for Increasing QM Sizea δO (ppm) PE/pcS-1 PE/pcS-2 TIP3P/pcS-1 TIP3P/pcS-2 exptl33
Figure 5. Convergence of absolute 13C NMR shielding constants in acrolein as a function of size of the QM region. The system size is fixed at Rsys = 12 Å. Results have been averaged over 120 MD snapshots. Calculations were performed at the KT3/pcS-1 level of theory.
a
RQM = 2 Å
RQM = 3 Å
RQM = 4 Å
81.8 77.9 68.9 64.2
86.2 71.0 79.6 64.7
85.4 69.8 81.1 65.4 75.5
All results were calculated with the KT3 xc-functional.
for increasing values of RQM. For the PE potential, calculated chemical shifts for the pcS-1 basis set are 81.8 ± 2.4 ppm for RQM = 2 Å and 85.4 ± 2.1 ppm for RQM = 4 Å (see Figure 6a). This is to be compared with the TIP3P results (shown in Figure 6b) of 68.9 ± 2.2 ppm for RQM = 2 Å and 81.1 ± 2.0 ppm for RQM = 4 Å. To elucidate the effect of using a more flexible basis set, we calculated the averaged chemical shifts using the pcS-2 basis set, which is of triple-ζ quality. These results are also presented in Table 1. With the PE potential and the pcS-2 basis set, the averaged calculated chemical shifts are 77.9, 71.0, and 69.8 ppm for RQM = 2, 3, and 4 Å, respectively. The difference in calculated chemical shifts between pcS-1 and pcS-2 for 17O in acetone is 14.6 ppm at RQM = 4 Å. This difference can mainly be attributed to the calculated shielding constant in the gas phase (see Supporting Information), which is not converged at the pcS-1 level of theory but is further decreased by 19.6 ppm when the basis set is increased to pcS-2. Further increase of the basis set size yields no appreciable benefit. On the other hand, basis functions on the solvent molecules for RQM > 1 Å in the QM/MM calculations help saturate the system with basis functions, and increasing the basis set from pcS-1 to pcS-2 therefore has a smaller effect, decreasing the shielding constant by 7.5 ppm for RQM = 3 Å and 5 ppm for RQM = 4 Å (calculated only for a single snapshot; see Supporting Information). For RQM = 2 Å, the system is less saturated and a change of 11.4 ppm is obtained. The experimental value for δO is 75.5 ppm.33 In conclusion, the pcS-1 basis set provides reasonable accuracy for the embedded calculation, provided there are some solvent molecules treated by QM. However, the gas-phase result is not converged with respect to the basis set and thus the calculated chemical shift deviates considerably. Use of the pcS-2 basis set ensures convergence also for the gas-phase calculation. The statistically averaged calculated chemical shifts obtained at the KT3/pcS-1 level of theory shown in Table 1 for the PE and TIP3P potentials indicate that the obtained results converge to different values: 85.4 and 81.1 ppm, respectively, at RQM = 4 Å. This is quite unsatisfactory, as for increasing sizes of the QM region the calculated values should converge toward
pronounced errors when the TIP3P potential is used. Note that even though the shielding constants calculated with the PE potential at RQM = 1 Å are quite close to the converged values, this is due to fortuitous error cancelation. The errors due to missing exchange−repulsion might as well lead to the opposite scenario (see Supporting Information). On the basis of these results, we find that converged shielding constants may be calculated efficiently with a rather small QM region because we use a high-quality embedding potential. In fact, a QM region as small as RQM = 2 Å may be employed, provided that the most essential water molecules are treated at the QM level. It does not come without problems, however, as the effect of missing exchange−repulsion has clear implications in the results. In a recent study, Flaig et al.8 employed a linear scaling approach to investigate the convergence behavior of shielding constants of large (intractable system sizes for conventional QM methods) molecular clusters (>800 atoms). The test systems in that study were large organic systems such as an aminopyrazole peptide or the molecular clip, both solvated in explicit water. Here, they found that quite large QM regions were needed in order to obtain converged NMR shielding constants (around 7 Å for 17O and 6 Å for 13C) compared to a full QM calculation. Our results indicate that, in general, the size of the QM region may be reduced (compared to conventional QM/MM methods) when an improved embedding potential is used, thus resulting in a very efficient computational strategy. Finally, we note that the shielding constants of nonpolar hydrogens are largely unaffected by the type of embedding potential, as shown in the Supporting Information. 3.2. 17O and 13C Nuclear Shielding Constants in Acetone. Results for the carbonyl chromophore in acetone are very similar to results obtained for acrolein, in terms of convergence with respect to both system and QM region size (see Supporting Information). This is not unexpected when the similarity of the chemical environments is considered. Furthermore, the conclusions and comments made above with respect to lack of important quantum mechanical effects 985
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
from RQM = 5 to 6 Å. On this basis, we conclude that the calculated shielding constants and hence the chemical shifts obtained are very sensitive to the quality of the embedding potential and converge substantially faster with a high-quality and polarizable embedding potential compared to the TIP3P potential. 3.2.2. Chemical Shifts of Carbonyl 13C in Acetone. Calculated chemical shifts for 13C are consolidated in Table 3 Table 3. Calculated Acetone 13C Chemical Shifts for Increasing QM Sizea δC (ppm) PE/pcS-1 PE/pcS-2 TIP3P/pcS-1 TIP3P/pcS-2 exptl34 a
RQM = 2 Å
RQM = 3 Å
RQM = 4 Å
−16.0 −15.9 −13.8 −13.6
−16.6 −14.8 −15.5 −13.8
−16.8 −14.9 −16.1 −14.1 −18.9
All results were calculated with the KT3 xc-functional.
for both PE and TIP3P potentials and are shown in Figure 7a for the PE potential and Figure 7b for the TIP3P potential. When the pcS-1 basis set is employed, calculated chemical shifts for 13C obtained with the PE potential have similar convergence behavior as for 17O, in that convergence is obtained quite fast. The calculated shifts are −16.0, −16.6, and −16.8 ppm for RQM = 2, 3, and 4 Å, respectively. The corresponding calculated chemical shifts for the TIP3P potential and the pcS-1 basis set are −13.8, −15.5, and −16.1 ppm, again illustrating the slower convergence of the TIP3P potential. Upon increasing the basis set to pcS-2, as was done for oxygen, the calculated chemical shifts with PE potential are −15.9, −14.8, and −14.9 ppm for RQM = 2, 3, and 4 Å, respectively. Again, we observe around 1 ppm deviation across the tested range of distances. Correspondingly, for TIP3P we obtain −13.6, −13.8, and −14.1 ppm. Surprisingly, the deviation is only 0.5 ppm. As we have discussed in the preceding sections, this favorable result toward the TIP3P potential is due to a fortuitous error cancelation stemming from a softer potential. The experimental chemical shift δC is −18.9 ppm.34 However, the calculated gas-phase shielding constant of 13 C, contrary to what was observed with 17O, is not fully converged when the pcS-2 basis set is used. The slow convergence of the TIP3P potential observed for 17 O is repeated for 13C in the carbonyl chromophore of acetone. In Table 4, calculated values of the shielding constants for 13C are presented. The convergence pattern that was observed for oxygen is repeated here, in that the calculated shielding constants obtained with the TIP3P potential show a relatively large shift when RQM is increased from 4 to 5 Å. Finally, we note that many other sources of error must be taken into account when calculated chemical shifts are
Figure 6. Calculated chemical shift for the oxygen atom of the carbonyl chromophore in acetone by use of (a) PE potential or (b) TIP3P potential. Here, Rsys = 11 Å and two different RQM sizes are shown. The solid black line is the experimentally observed value. Calculations were performed at the KT3/pcS-1 level of theory.
the same number. Due to the computational demand of calculating averaged shielding constants from 120 snapshots for larger values of RQM, we have used snapshot 16 as a probe with respect to the convergence behavior as RQM increases beyond 4 Å. In Table 2, absolute shielding constants obtained at the KT3/pcS-1 level of theory for 17O from both PE and TIP3P potentials are shown for RQM values up to 6 Å. Here we clearly see that the PE results have converged to within 1 ppm already at RQM = 3 Å. There are no significant changes in the calculated shielding constants for RQM = 4 or 5 Å. On the contrary, the TIP3P numbers look converged at RQM = 4 Å, because the change from RQM = 3 to 4 Å is only −0.7 ppm. However, upon increasing RQM to 5 Å, there is a large change (−4.4 ppm) in the calculated shielding constants. Increasing RQM to 6 Å yields a further change of −0.7 ppm. A similar test was done on snapshot 18 (see Supporting Information) to verify that this was not a single incident. It shows that the PE number again converges to within 1 ppm of the RQM = 6 Å result at RQM = 3 Å, whereas TIP3P shows a large −3.1 ppm change when going
Table 2. Calculated Shielding Constants for 17O in Acetone for Increasing QM Sizea 17
O shielding constant (ppm)
PE/pcS-1 TIP3P/pcS-1 a
RQM = 1 Å
RQM = 2 Å
RQM = 3 Å
RQM = 4 Å
RQM = 5 Å
RQM = 6 Å
−238.3 −270.1
−243.8 −260.1
−247.0 −254.3
−246.1 −253.6
−245.7 −249.2
−246.1 −248.5
All results were calculated with the KT3 xc-functional. 986
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
4. CONCLUSIONS We have calculated the isotropic NMR shielding constants of both acrolein and acetone in solution using density functional theory and an explicit description of the solvent molecules with either a polarizable embedding (PE) potential or the standard TIP3P potential. The reported results are statistically averaged, employing 120 snapshots for each of the molecular systems. We used the KT3 exchange−correlation functional along with polarization-consistent (pc) basis sets optimized for evaluation of shielding constants to ensure fast convergence with respect to system size, basis set size, and size of QM region. The TIP3P potential consists of simple point charges, whereas the polarizable embedding potential was derived with the LoProp method using the B3LYP xc-functional and the aug-cc-pVDZ basis set. We showed that a high-quality polarizable embedding potential in QM/MM-type calculations of isotropic NMR shielding constants is a feasible solution, provided that the first solvation shell is included in the QM region. Increasing it further yields little improvement. We focused on the carbonyl chromophore of acrolein and acetone because it (and its relatives) exhibits strong solvation effects in polar solvents. While there was virtually no difference in the convergence behavior with respect to cluster size between the two types of embedding potentials, there was a large difference when convergence with respect to QM size was considered. We found that, by use of a polarizable embedding potential, the required size of the QM region was at least 2 Å smaller for both systems compared to use of a standard classical potential for water. The basis set saturation that happens when solvent molecules are included in the QM region helps to accelerate the convergence with respect to basis set size. When calculating chemical shifts, one must be careful that the gas-phase calculations are also converged with respect to the basis set. We also found that lack of important QM effects such as exchange−repulsion lead to spurious (and sometimes fortuitous) error cancelation for the smallest QM regions. This is a point we intend to address in a future study.
Figure 7. Calculated chemical shift for the carbon atom of the carbonyl chromophore in acetone by use of (a) PE potential or (b) TIP3P potential. Here, Rsys = 11 Å and two different RQM sizes are shown. The solid black line is the experimentally observed value. Calculations were performed at the KT3/pcS-1 level of theory.
compared with experimental values.24 In the case of 17O, for instance, Kongsted et al.35 previously found that the choice of the KT3 xc-functional introduces an average error of 7 ppm in the solvent shift when compared to coupled-cluster predictions. Even though this deviation is less than those found by other xcfunctionals, 7 ppm is on the same scale as the errors we observe for the pcS-1 and pcS-2 basis sets. Furthermore, vibrational effects have been found to account for roughly 5 ppm,36 and this is neglected in the current study. Finally, the choice of potential used in the MD simulations is also critical, as recently highlighted by Lindorff-Larsen et al.37 and Beauchamp et al.38 Thus, in terms of basis set, the condensed phase results are definitely converged with pcS-1 for the absolute shielding constants; however, the choice of calculating the chemical shift forces us to use pcS-2 due to convergence issues related to calculation of gas-phase shielding tensors.
■
ASSOCIATED CONTENT
S Supporting Information *
Fourteen figures, showing averaged results of the remaining two carbon atoms and four hydrogen atoms of acrolein and acetone, and 11 tables, listing calculated absolute shielding constants for the carbonyl chromophores of acrolein and acetone with large basis sets and large QM regions. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail [email protected]. Notes
The authors declare no competing financial interest.
Table 4. Calculated Shielding Constants for 13C in Acetone for Increasing QM Sizea 13
C shielding constant (ppm)
PE/pcS-1 TIP3P/pcS-1 a
RQM = 1 Å
RQM = 2 Å
RQM = 3 Å
RQM = 4 Å
RQM = 5 Å
RQM = 6 Å
−21.3 −16.8
−20.7 −17.6
−19.7 −18.5
−20.5 −19.6
−20.8 −20.3
−20.7 −20.5
All results were calculated with the KT3 xc-functional. 987
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
■
Article
Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision A.02; Gaussian, Inc., Wallingford, CT, 2009. (20) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (21) Stephens, P.; Devlin, F.; Chabalowski, C.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (22) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (23) Keal, T. W.; Tozer, D. J. J. Chem. Phys. 2004, 121, 5654. (24) Jensen, F. J. Chem. Theory Comput. 2008, 4, 719−727. (25) Jorgensen, W. L. J. Am. Chem. Soc. 1981, 103, 335−340. (26) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (27) Gagliardi, L.; Lindh, R.; Karlstrom, G. J. Chem. Phys. 2004, 121, 4494−4500. (28) Aquilante, F.; De Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P. Å.; Neogrády, P.; Pedersen, T. B.; Pitoňaḱ , M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. J. Comput. Chem. 2010, 31, 224−247. (29) Aidas, K.; Møgelhøj, A.; Nielsen, C. B.; Mikkelsen, K. V.; Ruud, K.; Christiansen, O.; Kongsted, J. J. Phys. Chem. A 2007, 111, 4199− 4210. (30) Aidas, K.; Møgelhøj, A.; Nilsson, E. J. K.; Johnson, M. S.; Mikkelsen, K. V.; Christiansen, O.; Söderhjelm, P.; Kongsted, J. J. Chem. Phys. 2008, 128, 194503. (31) Freedman, D.; Diaconis, P. Probab. Theory Relat. Fields 1981, 57, 453−476. (32) Fradelos, G.; Wesołowski, T. A. J. Phys. Chem. A 2011, 115, 10018−10026. (33) Cossi, M.; Crescenzi, O. J. Chem. Phys. 2003, 118, 8863. (34) Tiffon, B.; Dubois, J.-E. Org. Magn. Reson. 1978, 11, 295−302. (35) Kongsted, J.; Aidas, K.; Mikkelsen, K. V.; Sauer, S. P. J. Chem. Theory Comput. 2008, 4, 267−277. (36) Kongsted, J.; Ruud, K. Chem. Phys. Lett. 2008, 451, 226−232. (37) Lindorff-Larsen, K.; Maragakis, P.; Piana, S.; Eastwood, M. P.; Dror, R. O.; Shaw, D. E. PLoS One 2012, 7, No. e32131. (38) Beauchamp, K. A.; Lin, Y.-S.; Das, R.; Pande, V. S. J. Chem. Theory Comput. 2012, 8, 1409−1414.
ACKNOWLEDGMENTS We are grateful for computational resources provided by the Danish Center for Scientific Computing. J.K. thanks the Danish Councils for Independent Research (the Sapere Aude programme), the Lundbeck Foundation, and the Villum foundation for financial support.
■
REFERENCES
(1) Barrett, P. J.; Chen, J.; Cho, M.-K.; Kim, J.-H.; Lu, Z.; Mathew, S.; Peng, D.; Song, Y.; Van Horn, W. D.; Zhuang, T.; Sönnichsen, F. D.; Sanders, C. R. Biochemistry 2013, 52, 1303−1320. (2) Robustelli, P.; Stafford, K. A.; Palmer, A. G. J. Am. Chem. Soc. 2012, 134, 6365−6374. (3) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999−3094. (4) Olsen, J. M.; Aidas, K.; Kongsted, J. J. Chem. Theory Comput. 2010, 6, 3721−3734. (5) Schwabe, T.; Olsen, J. M. H.; Sneskov, K.; Kongsted, J.; Christiansen, O. J. Chem. Theory Comput. 2011, 7, 2209−2217. (6) Olsen, J. M. H.; Kongsted, J. Adv. Quantum Chem. 2011, 61, 107−143. (7) Olsen, J. M. H. Ph.D. thesis, University of Southern Denmark, Odense, Denmark, 2012; DOI 10.6084/m9.figshare.156852. (8) Flaig, D.; Beer, M.; Ochsenfeld, C. J. Chem. Theory Comput. 2012, 8, 2260−2271. (9) Helgaker, T.; Jaszunski, M.; Ruud, K. Chem. Rev. 1999, 99, 293. (10) Helgaker, T.; Jorgensen, P. J. Chem. Phys. 1991, 95, 2595−2601. (11) London, F. J. Phys. Radium 1937, 8, 397−409. (12) Wolinski, K.; Hinton, J. F.; Pulay, P. J. Am. Chem. Soc. 1990, 112, 8251−8260. (13) Kongsted, J.; Nielsen, C. B.; Mikkelsen, K. V.; Christiansen, O.; Ruud, K. J. Chem. Phys. 2007, 126, No. 034510. (14) Helgakar, T.; Jørgensen, P.; Olsen, J. Molecular ElectronicStructure Theory; Wiley: New York, 2000; pp 2−5. (15) Applequist, J.; Carl, J. R.; Fung, K.-K. J. Am. Chem. Soc. 1972, 94, 2952−2960. (16) Gao, B.; Thorvaldsen, A. J.; Ruud, K. Int. J. Quantum Chem. 2011, 111, 858−872. (17) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jansik, B.; Jensen, H. J. A.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V.; Salek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; Sylvester-Hvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; Ågren, H. WIREs Comput. Mol. Sci. 2013, DOI: 10.1002/wcms.1172. (18) Dalton, a Molecular Electronic Structure Program, development version, 2013; http://daltonprogram.org/; accessed October 14, 2013. (19) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; 988
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Nuclear Magnetic Shielding Constants from Quantum Mechanical/ Molecular Mechanical Calculations Using Polarizable Embedding: Role of the Embedding Potential Casper Steinmann,* Jógvan Magnus Haugaard Olsen, and Jacob Kongsted Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, 5230 Odense, Denmark S Supporting Information *
ABSTRACT: We present NMR shielding constants obtained through quantum mechanical/molecular mechanical (QM/ MM) embedding calculations. Contrary to previous reports, we show that a relatively small QM region is sufficient, provided that a high-quality embedding potential is used. The calculated averaged NMR shielding constants of both acrolein and acetone solvated in water are based on a number of snapshots extracted from classical molecular dynamics simulations. We focus on the carbonyl chromophore in both molecules, which shows large solvation effects, and we study the convergence of shielding constants with respect to total system size and size of the QM region. By using a high-quality embedding potential over standard point charge potentials, we show that the QM region can be made at least 2 Å smaller without any loss of quality, which makes calculations on ensembles tractable by conventional density functional theory calculations.
1. INTRODUCTION The broad range of applications to which NMR spectroscopy can be applied to solve scientific problems were recently highlighted in a paper by Barrett et al.1 Recent years have witnessed an increased interest in using calculated NMR data in combination with experiments in order to elucidate, for example, structural properties of complex materials such as proteins.2 In this same period, methods that can be used to predict NMR properties have also greatly increased their level of sophistication. Indeed, it is today possible to use ab initio methods to calculate very accurate values for key NMR parameters for smaller molecular systems. For larger systems and usually of biological interestit has become popular to use methods relying on some sort of parametrization of the NMR parameters directly in terms of structural data; that is, the NMR parameters are predicted on the basis of simple analytical and parametrized relationships, giving as input the molecular structure. Even though such parametrized schemes can be tuned to become rather accurate, predictions based directly on ab initio methods are surely to be preferred. One class of large molecular systems, which are often met in chemistry, is that of a solvated molecule. For solute−solvent systems, which is the subject of this paper, continuum approaches have often been used in order to provide a representation of the solvent.3 However, explicit treatment of the solvent, including dynamical effects from classical molecular dynamics (MD) simulations, have gained popularity. It is well-known that it is often sufficient to use a cluster model with the solute and sometimes also its nearest-neighbor solvent molecules treated by high-level quantum mechanics (QM), while the remainder of the solvent environment is treated classically, often using simple point © 2014 American Chemical Society
charge potentials taken from standard molecular mechanics (MM) force fields. It has, however, been shown that use of ab initio-derived polarizable potentials that can adapt to the quantum mechanically treated charge distribution is preferable in order to obtain good accuracy.4,5 In the present paper we demonstrate that, by using such high-quality embedding potential, one can obtain converged nuclear magnetic shielding constants (and chemical shifts) for even quite small QM regions compared to the case where the environment is modeled by a standard point charge model. We show that an important aspect of embedding calculations is the selfconsistent account of polarization effects such as in the recently developed polarizable embedding approach by Olsen et al.4,6,7 Even though it has previously been found that large QM regions were needed to obtain converged NMR shielding constants,8 we here show that the QM size can be reduced by using a better embedding potential. Thus, we show that the polarizable embedding strategy represents a computationally efficient route if the aim is to calculate key NMR properties of only a small part of a larger molecular system. In the present paper our focus is on solute−solvent systems, but the aim is to extend the formalism presented in this work to polymeric systems such as proteins. This paper presents the underlying theory behind calculations of gauge-origin-independent NMR shielding constants within the polarizable embedding formalism. Received: October 9, 2013 Published: February 14, 2014 981
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
2. METHODOLOGY 2.1. NMR Shielding Tensors in the Polarizable Embedding Framework. The calculation of an element of the NMR shielding tensor for nucleus N, σNij , in self-consistent field theory is calculated according to9 σijN = 1 +
∑ Dμν μν
∂ 2hμν ∂Bi ∂mj , N
+
∑ μν
moment expansions and the induction term on distributed polarizabilites. In the following we will use a multi-index notation that allows more compact expressions. We use a three-dimensional multi-index defined as k = (kx, ky, kz), whose norm is given by the sum of its elements, |k| = kx + ky + kz, and factorial is the product between the factorial of the elements, k! = kx!ky!kz!. The multi-index power of a vector in Cartesian space is defined as rk = xkxykyzkz, which allows us to write a general partial derivative operator as ∇kr = (∂|k|)/(∂xkx ∂yky ∂zkz). By use of this notation, the electrostatic part of the PE operator is given as4
∂Dμν ∂hμν ∂Bi ∂mj , N
(1)
Here, Dμν is an element of the density matrix in the atomic orbital (AO) basis, hμν is an element of the one-electron Hamiltonian also in AO basis, Bi is a component of the magnetic field (induction), and mj,N is a component of the nuclear magnetic moment of nucleus N. Evaluation of the last term in eq 1, that is, the paramagnetic contribution, requires knowledge of the perturbed density matrix, which in our case is obtained by solving a set of appropriate response equations10 for each of the three components of the magnetic field. In order to remove the dependency on the gauge origin in the calculated chemical shielding tensors, gauge-including atomic orbitals (GIAOs)10−12 are employed; that is, the AO basis functions have attached a phase that carries an explicit reference to the perturbation: χμN (B)
⎤ ⎡ i = exp⎢ − (B × R μ) ·r⎥ χμN (0) ⎦ ⎣ 2
K
veŝ =
|k|= 0 s = 1
KS
ts(,kpq) = − p ∇kr
M
⎤
m=1
⎦
∑ Tsm(k)Zm⎥⎥
(5)
1 q |r − R s|
(6)
and T(k) sm is an element of the interaction tensor defined as (k) Tsm = ∇km
(2)
1 |R m − R s|
(7)
For each value of |k| in the summation in eq 5, there is an implicit summation over (|k| + 1)(|k| + 2)/2 multi-indices; for example, for |k| = 1 it is k = (1, 0, 0), k = (0, 1, 0), and k = (0, 0, 1). The summation runs up to K, which is the highest order of multipoles included in the embedding potential. We use a second quantization formulation where Ê pq is a singlet excitation operator expressed as a string of creation a†pσ and annihilation operators.14 The v̂ind operator is given as S
vind ̂ = −∑
∑ μsind,(k) ∑ ts(,kpq) Epq̂
|k|= 1 s = 1
(8)
pq
where the summation over multi-indices is restricted to |k| = 1 and μind,(k) is an element of an induced dipole at site s. Formally, s the induced dipoles are determined from a set of coupled linear equations that can be expressed in matrix form as
μind = AF
(9) 15
where A is the classical response matrix (also known as the relay matrix) and F is the combined field from the QM electrons and QM nuclei and the field from the permanent multipoles in the environment: F = Fe + Fn + Fmul
(10)
Here, μ and F are supervectors of length 3S, where S is the number of polarizable sites. Note that the number of expansion points for the multipoles and the polarizabilities need not be equal. The induced dipoles can be solved for by use of eq 9. This can be done directly, which requires explicit construction of the classical response matrix. However, due to the generally large size of the matrix (3S × 3S), eq 9 is usually solved via iterative techniques. Due to the use of GIAOs in combination with the embedding potential, specific contributions to the response equations appear. These contributions are here presented as a ind
KS
(3)
̂ where fKS vac is the gas-phase Kohn−Sham operator and v̂PE is the polarizable embedding (PE) operator given as vPE ̂ = veŝ + vind ̂
⎡ ( −1)|k| (k)⎢ ̂ + Ms ∑ ts(,kpq) Epq ⎢⎣ k! pq
Here, M(k) is the kth element of the |k|th order multipole s moment at a site s; that is, for k = (0, 0, 0) it is a charge, for k = (1, 0, 0) it is the x component of a dipole, and so on. t(k) s,pq are integrals over the |k|th order derivative of |r − Rs|−1:
Here, Rμ is the vector between the global gauge origin and the position of the field-independent AO, χNμ (0), attached to nucleus N. As shown and discussed in detail previously,13 contributions from the embedding potential enter into eq 1 in two places: in the first term (the diamagnetic contribution) through the density matrix and in the second term (the paramagnetic contribution) through the perturbed density matrix. The polarizable embedding (PE) scheme represents a model for explicit inclusion of environmental effects in quantum chemical calculations.4,6,7 Within the embedding scheme, the environment is described by use of a multicenter multipole expansion of the electric potential created by the environment. Usually the expansion centers are taken at the position of the atomic nuclei, but other choices might also be relevant. Polarization of the environment is accounted for through the use of distributed polarizabilities that give rise to induced dipoles. The PE method may be coupled to different quantum chemical formulations, and the most used of these is currently PE in combination with density functional theory (DFT), referred to as PE-DFT.4 This method is available for both energies and ground-state molecular properties as well as excited-states properties through the use of time-dependent DFT (PE-TD-DFT). In the PE-DFT scheme, an effective ̂ Kohn−Sham (KS) operator fKS eff is given as
̂ + v̂ feff̂ = fvac PE
S
∑∑
(4)
Here, v̂es and v̂ind are the one-electron operators for electrostatic interactions between environment and QM core and for induction effects (polarization) of the environment, respectively. The electrostatic term is based on distributed multipole 982
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
multipole moments and polarizabilities are distributed onto the atomic sites of the molecules that define the environment. To evaluate the performance of the embedding scheme in calculation of shielding constants determined under the influence of either PE or TIP3P potentials, two different systems based on earlier studies are investigated. For this we have used 120 evenly distributed snapshots of acetone and acrolein (Figure 1) solvated in water, obtained from a MD
generalization of the terms described in ref 13 to reflect the current implementation. In short, in the field-dependent AO basis, the elements of the v̂PE operator take the form ⟨χμ (B)|vPE ̂ |χν (B)⟩ = ⟨χμ (B)|veŝ |χν (B)⟩ + ⟨χμ (B)|vind ̂ |χν (B)⟩ K
=
S
|k|= 0 s = 1 S
−
|k |
∑ ∑ (−1) k!
Ms(k)ts(,kμν) (B)
∑ ∑ μsind,(k)ts(,kμν) (B) |k|= 1 s = 1
(11)
where the integrals now carry an explicit dependence on the magnetic field due to the use of GIAOs. The contributions needed in order to determine the first-order density matrix (eq 1) require evaluation of the first-order derivative of eq 11 with respect to a component of the magnetic field: Figure 1. Depictions of acrolein (left) and acetone (right) with atomic labels.
d⟨χμ (B)|vPE ̂ |χν (B)⟩ dBi K
=
S
|k |
∑ ∑ ( − 1) k!
|k|= 0 s = 1 S
−
∑ ∑ μsind,(k) |k|= 1 s = 1
Ms(k)
dts(,kμν) (B)
simulation described in works by Aidas et al.29,30 The 120 snapshots were used in order to perform statistical averaging. We will analyze the convergence of the calculated shielding constants with respect to size of the QM region and size of the cluster, that is, the whole system (QM + environment), as it approaches bulk. The size of the QM region and total system are defined as minimum atom-based distances from the solute to the solvent molecules. The size of the QM region is determined in such a way that if a solvent molecule has at least one atom within a distance RQM from the solute, it is included in the QM region. Otherwise it is treated classically. Similarly, the size of the molecular cluster is determined in such a way that if an atom of a solvent molecule is within Rsys it is included; otherwise it is removed. The size of the whole system is always greater than or equal to the size of the QM region (Rsys ≥ RQM). Chemical shifts were obtained by subtracting the shielding constants in the isolated molecule from those in the solvated molecule:
dBi
dts(,kμν) (B) dBi
(12)
where the integral derivatives are defined according to dts(,kμν) (B) dBi
dχμ (B)
=
dBi +
∇kr
χμ (B) ∇kr
1 χ (B ) |r − R s| ν dχν (B) 1 dBi |r − R s|
(13)
These integrals are evaluated analytically for zeroth- and firstorder derivatives of |r − Rs|−1.13 The integrals corresponding to higher-order derivatives are evaluated on the basis of an implemented finite difference scheme using the analytically calculated integrals. No problems with stability were observed in the evaluation of second-order derivative integrals. In a future implementation, we expect to perform the integral evaluation completely analytically to any order by use of the Gen1Int16 package. 2.2. Computational Details. The GIAO contributions to the NMR shielding constants calculated via PE-DFT were added to a development version of the Dalton17,18 program using the Polarizable Embedding library7 and will be available in a future release. Geometry optimizations of vacuum structures were performed with Gaussian 0919 using the B3LYP20,21 exchange−correlation (xc) functional and the augcc-pVTZ22 basis set. Shielding constants were calculated with Dalton18 using the KT323 GGA xc-functional and the pcS-1 and pcS-2 polarization-consistent basis sets by Jensen,24 optimized for calculating nuclear magnetic shielding constants. For the classical water potentials, either TIP3P25,26 or LoProp27 potentials were used. The LoProp potentials were calculated with Molcas28 at the B3LYP/aug-cc-pVDZ level of theory for all multipole moments through quadrupoles and dipole polarizabilities (this potential is denoted PE in the following). Note that we used an ANO-type recontraction of the aug-ccpVDZ basis set as required by the LoProp procedure. All
δ = σsol − σvac
(14)
The presented error bars and uncertainties are calculated as standard error of the mean (SEM) of the shielding constants for a set of snapshots as s SEMN = (15) 120 where s is the standard deviation of the calculated NMR shielding constants for a nucleus N and 120 is the number of snapshots used to calculate the shielding constants. Optimal bin widths for histograms are calculated according to the Freedman−Diaconis rule.31
3. RESULTS We shall initially be concerned with the convergence of the shielding constants as a function of size of the system (Rsys) and size of the QM region (RQM). For this we have chosen acrolein, and we will investigate in detail the convergence of carbonyl chromophore 17O and 13C shielding constants. For all 120 MD snapshots, we calculated the shielding constants for Rsys ∈ [1, 12] Å and RQM ∈ [1, 4] Å. We used the pcS-1 basis set for this part of the study because the need for large (RQM = 4 Å) QM 983
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
calculations would otherwise be computationally too demanding, given the large number of calculations that is required. As we will discuss later, the choice of pcS-1 is satisfactory for evaluation of the shielding constants from QM/MM calculations but is not enough for prediction of accurate gas-phase shielding constants. 3.1. 17O and 13C Nuclear Shielding Constants in Acrolein. Convergence of the 17O shielding constant in acrolein, calculated at the PE-KT3/pcS-1 level of theory, with respect to system size is presented in Figure 2 for several values
Figure 3. Convergence of absolute 17O NMR shielding constant in acrolein, calculated with the TIP3P potential as a function of system size and size of QM region. The results have been averaged over 120 MD snapshots. The first value of each series is a pure QM calculation where RQM = Rsys. Calculations were performed at the KT3/pcS-1 level of theory.
Figure 2. Convergence of absolute 17O NMR shielding constant in acrolein, calculated with the PE potential as a function of system size and size of QM region. Results have been averaged over 120 MD snapshots. The first value of each series is a pure QM calculation where RQM = Rsys. Calculations were performed at the KT3/pcS-1 level of theory.
of RQM. The shielding constant increases very rapidly as the size of the entire complex is increased but slows down and eventually converges at around Rsys = 8 Å. Increasing Rsys to 12 Å yields little improvement and the converged values are here −221.0, −226.3, −223.1, and −223.2 ppm for RQM = 1, 2, 3, and 4 Å, respectively. The computational cost of using Rsys = 12 Å over Rsys = 8 Å is negligible, and we have therefore used this size in the following calculations. In Figure 3 the corresponding results are shown for the TIP3P potential. In this case the converged values at Rsys = 12 Å are −256.4, −244.5, −232.1, and −228.1 ppm for RQM = 1, 2, 3, and 4 Å, respectively. We note that there is no observable difference in the convergence behavior with respect to system size for both PE and TIP3P; however, the effect of size of the QM region varies greatly between the two, with the PE results showing a much faster convergence with respect to size of the QM region. The shielding constants of 17O based on the use of the PE and TIP3P potentials as a function of size of the QM region (system size is fixed at Rsys = 12 Å) are shown in Figure 4. At RQM = 1 Å the calculated PE- and TIP3P-based average shielding constants deviate considerably from the expected convergence pattern. The reason is most likely that we are missing important quantum mechanical effects, since at RQM = 1 Å it is only acrolein that is modeled with quantum mechanics while all solvent molecules are treated classically and we therefore completely neglect exchange−repulsion with the
Figure 4. Convergence of absolute 17O shielding constants in acrolein as a function of size of the QM region. The system size is fixed at Rsys = 12 Å. Results have been averaged over 120 MD snapshots. Calculations were performed at the KT3/pcS-1 level of theory.
solvent. At RQM = 2 Å, however, we include hydrogen-bonded solvent molecules in the QM region and thereby also recover important quantum-mechanical interactions with the solvent. This is in line with the study by Fradelos and Wesolowski32 where they demonstrated that the lack of exchange−repulsion can introduce large errors in pure Coulombic embedding models especially when large or diffuse basis sets are used. The different behaviors of PE and TIP3P potentials can be explained by the fact that the PE potential offers a very accurate electrostatic potential compared to TIP3P.4,5 The same errors, due to the lack of important short-range quantum mechanical effects, are also present in the calculations of the shielding constants of 13C in the acrolein carbonyl chromophore (see Figure 5), and here we also observe more 984
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
for small QM regions are also valid for acetone. On this basis, we do not consider a detailed analysis of the RQM = 1 Å results. Note that the results below are all obtained for Rsys = 11 Å. Rather than evaluating the shielding constants for acetone, we find it more instructive to calculate the gas-to-aqueous solution chemical shifts of the carbonyl chromophore because it introduces additional requirements on the quality of the calculation, as we shall see below. 3.2.1. Chemical Shift of Carbonyl 17O in Acetone. The calculated gas-to-aqueous solution chemical shift of 17O in acetone is obtained according to eq 14 and is listed in Table 1 Table 1. Calculated Acetone 17O Chemical Shifts for Increasing QM Sizea δO (ppm) PE/pcS-1 PE/pcS-2 TIP3P/pcS-1 TIP3P/pcS-2 exptl33
Figure 5. Convergence of absolute 13C NMR shielding constants in acrolein as a function of size of the QM region. The system size is fixed at Rsys = 12 Å. Results have been averaged over 120 MD snapshots. Calculations were performed at the KT3/pcS-1 level of theory.
a
RQM = 2 Å
RQM = 3 Å
RQM = 4 Å
81.8 77.9 68.9 64.2
86.2 71.0 79.6 64.7
85.4 69.8 81.1 65.4 75.5
All results were calculated with the KT3 xc-functional.
for increasing values of RQM. For the PE potential, calculated chemical shifts for the pcS-1 basis set are 81.8 ± 2.4 ppm for RQM = 2 Å and 85.4 ± 2.1 ppm for RQM = 4 Å (see Figure 6a). This is to be compared with the TIP3P results (shown in Figure 6b) of 68.9 ± 2.2 ppm for RQM = 2 Å and 81.1 ± 2.0 ppm for RQM = 4 Å. To elucidate the effect of using a more flexible basis set, we calculated the averaged chemical shifts using the pcS-2 basis set, which is of triple-ζ quality. These results are also presented in Table 1. With the PE potential and the pcS-2 basis set, the averaged calculated chemical shifts are 77.9, 71.0, and 69.8 ppm for RQM = 2, 3, and 4 Å, respectively. The difference in calculated chemical shifts between pcS-1 and pcS-2 for 17O in acetone is 14.6 ppm at RQM = 4 Å. This difference can mainly be attributed to the calculated shielding constant in the gas phase (see Supporting Information), which is not converged at the pcS-1 level of theory but is further decreased by 19.6 ppm when the basis set is increased to pcS-2. Further increase of the basis set size yields no appreciable benefit. On the other hand, basis functions on the solvent molecules for RQM > 1 Å in the QM/MM calculations help saturate the system with basis functions, and increasing the basis set from pcS-1 to pcS-2 therefore has a smaller effect, decreasing the shielding constant by 7.5 ppm for RQM = 3 Å and 5 ppm for RQM = 4 Å (calculated only for a single snapshot; see Supporting Information). For RQM = 2 Å, the system is less saturated and a change of 11.4 ppm is obtained. The experimental value for δO is 75.5 ppm.33 In conclusion, the pcS-1 basis set provides reasonable accuracy for the embedded calculation, provided there are some solvent molecules treated by QM. However, the gas-phase result is not converged with respect to the basis set and thus the calculated chemical shift deviates considerably. Use of the pcS-2 basis set ensures convergence also for the gas-phase calculation. The statistically averaged calculated chemical shifts obtained at the KT3/pcS-1 level of theory shown in Table 1 for the PE and TIP3P potentials indicate that the obtained results converge to different values: 85.4 and 81.1 ppm, respectively, at RQM = 4 Å. This is quite unsatisfactory, as for increasing sizes of the QM region the calculated values should converge toward
pronounced errors when the TIP3P potential is used. Note that even though the shielding constants calculated with the PE potential at RQM = 1 Å are quite close to the converged values, this is due to fortuitous error cancelation. The errors due to missing exchange−repulsion might as well lead to the opposite scenario (see Supporting Information). On the basis of these results, we find that converged shielding constants may be calculated efficiently with a rather small QM region because we use a high-quality embedding potential. In fact, a QM region as small as RQM = 2 Å may be employed, provided that the most essential water molecules are treated at the QM level. It does not come without problems, however, as the effect of missing exchange−repulsion has clear implications in the results. In a recent study, Flaig et al.8 employed a linear scaling approach to investigate the convergence behavior of shielding constants of large (intractable system sizes for conventional QM methods) molecular clusters (>800 atoms). The test systems in that study were large organic systems such as an aminopyrazole peptide or the molecular clip, both solvated in explicit water. Here, they found that quite large QM regions were needed in order to obtain converged NMR shielding constants (around 7 Å for 17O and 6 Å for 13C) compared to a full QM calculation. Our results indicate that, in general, the size of the QM region may be reduced (compared to conventional QM/MM methods) when an improved embedding potential is used, thus resulting in a very efficient computational strategy. Finally, we note that the shielding constants of nonpolar hydrogens are largely unaffected by the type of embedding potential, as shown in the Supporting Information. 3.2. 17O and 13C Nuclear Shielding Constants in Acetone. Results for the carbonyl chromophore in acetone are very similar to results obtained for acrolein, in terms of convergence with respect to both system and QM region size (see Supporting Information). This is not unexpected when the similarity of the chemical environments is considered. Furthermore, the conclusions and comments made above with respect to lack of important quantum mechanical effects 985
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
from RQM = 5 to 6 Å. On this basis, we conclude that the calculated shielding constants and hence the chemical shifts obtained are very sensitive to the quality of the embedding potential and converge substantially faster with a high-quality and polarizable embedding potential compared to the TIP3P potential. 3.2.2. Chemical Shifts of Carbonyl 13C in Acetone. Calculated chemical shifts for 13C are consolidated in Table 3 Table 3. Calculated Acetone 13C Chemical Shifts for Increasing QM Sizea δC (ppm) PE/pcS-1 PE/pcS-2 TIP3P/pcS-1 TIP3P/pcS-2 exptl34 a
RQM = 2 Å
RQM = 3 Å
RQM = 4 Å
−16.0 −15.9 −13.8 −13.6
−16.6 −14.8 −15.5 −13.8
−16.8 −14.9 −16.1 −14.1 −18.9
All results were calculated with the KT3 xc-functional.
for both PE and TIP3P potentials and are shown in Figure 7a for the PE potential and Figure 7b for the TIP3P potential. When the pcS-1 basis set is employed, calculated chemical shifts for 13C obtained with the PE potential have similar convergence behavior as for 17O, in that convergence is obtained quite fast. The calculated shifts are −16.0, −16.6, and −16.8 ppm for RQM = 2, 3, and 4 Å, respectively. The corresponding calculated chemical shifts for the TIP3P potential and the pcS-1 basis set are −13.8, −15.5, and −16.1 ppm, again illustrating the slower convergence of the TIP3P potential. Upon increasing the basis set to pcS-2, as was done for oxygen, the calculated chemical shifts with PE potential are −15.9, −14.8, and −14.9 ppm for RQM = 2, 3, and 4 Å, respectively. Again, we observe around 1 ppm deviation across the tested range of distances. Correspondingly, for TIP3P we obtain −13.6, −13.8, and −14.1 ppm. Surprisingly, the deviation is only 0.5 ppm. As we have discussed in the preceding sections, this favorable result toward the TIP3P potential is due to a fortuitous error cancelation stemming from a softer potential. The experimental chemical shift δC is −18.9 ppm.34 However, the calculated gas-phase shielding constant of 13 C, contrary to what was observed with 17O, is not fully converged when the pcS-2 basis set is used. The slow convergence of the TIP3P potential observed for 17 O is repeated for 13C in the carbonyl chromophore of acetone. In Table 4, calculated values of the shielding constants for 13C are presented. The convergence pattern that was observed for oxygen is repeated here, in that the calculated shielding constants obtained with the TIP3P potential show a relatively large shift when RQM is increased from 4 to 5 Å. Finally, we note that many other sources of error must be taken into account when calculated chemical shifts are
Figure 6. Calculated chemical shift for the oxygen atom of the carbonyl chromophore in acetone by use of (a) PE potential or (b) TIP3P potential. Here, Rsys = 11 Å and two different RQM sizes are shown. The solid black line is the experimentally observed value. Calculations were performed at the KT3/pcS-1 level of theory.
the same number. Due to the computational demand of calculating averaged shielding constants from 120 snapshots for larger values of RQM, we have used snapshot 16 as a probe with respect to the convergence behavior as RQM increases beyond 4 Å. In Table 2, absolute shielding constants obtained at the KT3/pcS-1 level of theory for 17O from both PE and TIP3P potentials are shown for RQM values up to 6 Å. Here we clearly see that the PE results have converged to within 1 ppm already at RQM = 3 Å. There are no significant changes in the calculated shielding constants for RQM = 4 or 5 Å. On the contrary, the TIP3P numbers look converged at RQM = 4 Å, because the change from RQM = 3 to 4 Å is only −0.7 ppm. However, upon increasing RQM to 5 Å, there is a large change (−4.4 ppm) in the calculated shielding constants. Increasing RQM to 6 Å yields a further change of −0.7 ppm. A similar test was done on snapshot 18 (see Supporting Information) to verify that this was not a single incident. It shows that the PE number again converges to within 1 ppm of the RQM = 6 Å result at RQM = 3 Å, whereas TIP3P shows a large −3.1 ppm change when going
Table 2. Calculated Shielding Constants for 17O in Acetone for Increasing QM Sizea 17
O shielding constant (ppm)
PE/pcS-1 TIP3P/pcS-1 a
RQM = 1 Å
RQM = 2 Å
RQM = 3 Å
RQM = 4 Å
RQM = 5 Å
RQM = 6 Å
−238.3 −270.1
−243.8 −260.1
−247.0 −254.3
−246.1 −253.6
−245.7 −249.2
−246.1 −248.5
All results were calculated with the KT3 xc-functional. 986
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
Article
4. CONCLUSIONS We have calculated the isotropic NMR shielding constants of both acrolein and acetone in solution using density functional theory and an explicit description of the solvent molecules with either a polarizable embedding (PE) potential or the standard TIP3P potential. The reported results are statistically averaged, employing 120 snapshots for each of the molecular systems. We used the KT3 exchange−correlation functional along with polarization-consistent (pc) basis sets optimized for evaluation of shielding constants to ensure fast convergence with respect to system size, basis set size, and size of QM region. The TIP3P potential consists of simple point charges, whereas the polarizable embedding potential was derived with the LoProp method using the B3LYP xc-functional and the aug-cc-pVDZ basis set. We showed that a high-quality polarizable embedding potential in QM/MM-type calculations of isotropic NMR shielding constants is a feasible solution, provided that the first solvation shell is included in the QM region. Increasing it further yields little improvement. We focused on the carbonyl chromophore of acrolein and acetone because it (and its relatives) exhibits strong solvation effects in polar solvents. While there was virtually no difference in the convergence behavior with respect to cluster size between the two types of embedding potentials, there was a large difference when convergence with respect to QM size was considered. We found that, by use of a polarizable embedding potential, the required size of the QM region was at least 2 Å smaller for both systems compared to use of a standard classical potential for water. The basis set saturation that happens when solvent molecules are included in the QM region helps to accelerate the convergence with respect to basis set size. When calculating chemical shifts, one must be careful that the gas-phase calculations are also converged with respect to the basis set. We also found that lack of important QM effects such as exchange−repulsion lead to spurious (and sometimes fortuitous) error cancelation for the smallest QM regions. This is a point we intend to address in a future study.
Figure 7. Calculated chemical shift for the carbon atom of the carbonyl chromophore in acetone by use of (a) PE potential or (b) TIP3P potential. Here, Rsys = 11 Å and two different RQM sizes are shown. The solid black line is the experimentally observed value. Calculations were performed at the KT3/pcS-1 level of theory.
compared with experimental values.24 In the case of 17O, for instance, Kongsted et al.35 previously found that the choice of the KT3 xc-functional introduces an average error of 7 ppm in the solvent shift when compared to coupled-cluster predictions. Even though this deviation is less than those found by other xcfunctionals, 7 ppm is on the same scale as the errors we observe for the pcS-1 and pcS-2 basis sets. Furthermore, vibrational effects have been found to account for roughly 5 ppm,36 and this is neglected in the current study. Finally, the choice of potential used in the MD simulations is also critical, as recently highlighted by Lindorff-Larsen et al.37 and Beauchamp et al.38 Thus, in terms of basis set, the condensed phase results are definitely converged with pcS-1 for the absolute shielding constants; however, the choice of calculating the chemical shift forces us to use pcS-2 due to convergence issues related to calculation of gas-phase shielding tensors.
■
ASSOCIATED CONTENT
S Supporting Information *
Fourteen figures, showing averaged results of the remaining two carbon atoms and four hydrogen atoms of acrolein and acetone, and 11 tables, listing calculated absolute shielding constants for the carbonyl chromophores of acrolein and acetone with large basis sets and large QM regions. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail [email protected]. Notes
The authors declare no competing financial interest.
Table 4. Calculated Shielding Constants for 13C in Acetone for Increasing QM Sizea 13
C shielding constant (ppm)
PE/pcS-1 TIP3P/pcS-1 a
RQM = 1 Å
RQM = 2 Å
RQM = 3 Å
RQM = 4 Å
RQM = 5 Å
RQM = 6 Å
−21.3 −16.8
−20.7 −17.6
−19.7 −18.5
−20.5 −19.6
−20.8 −20.3
−20.7 −20.5
All results were calculated with the KT3 xc-functional. 987
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
Journal of Chemical Theory and Computation
■
Article
Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision A.02; Gaussian, Inc., Wallingford, CT, 2009. (20) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (21) Stephens, P.; Devlin, F.; Chabalowski, C.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (22) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (23) Keal, T. W.; Tozer, D. J. J. Chem. Phys. 2004, 121, 5654. (24) Jensen, F. J. Chem. Theory Comput. 2008, 4, 719−727. (25) Jorgensen, W. L. J. Am. Chem. Soc. 1981, 103, 335−340. (26) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (27) Gagliardi, L.; Lindh, R.; Karlstrom, G. J. Chem. Phys. 2004, 121, 4494−4500. (28) Aquilante, F.; De Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P. Å.; Neogrády, P.; Pedersen, T. B.; Pitoňaḱ , M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. J. Comput. Chem. 2010, 31, 224−247. (29) Aidas, K.; Møgelhøj, A.; Nielsen, C. B.; Mikkelsen, K. V.; Ruud, K.; Christiansen, O.; Kongsted, J. J. Phys. Chem. A 2007, 111, 4199− 4210. (30) Aidas, K.; Møgelhøj, A.; Nilsson, E. J. K.; Johnson, M. S.; Mikkelsen, K. V.; Christiansen, O.; Söderhjelm, P.; Kongsted, J. J. Chem. Phys. 2008, 128, 194503. (31) Freedman, D.; Diaconis, P. Probab. Theory Relat. Fields 1981, 57, 453−476. (32) Fradelos, G.; Wesołowski, T. A. J. Phys. Chem. A 2011, 115, 10018−10026. (33) Cossi, M.; Crescenzi, O. J. Chem. Phys. 2003, 118, 8863. (34) Tiffon, B.; Dubois, J.-E. Org. Magn. Reson. 1978, 11, 295−302. (35) Kongsted, J.; Aidas, K.; Mikkelsen, K. V.; Sauer, S. P. J. Chem. Theory Comput. 2008, 4, 267−277. (36) Kongsted, J.; Ruud, K. Chem. Phys. Lett. 2008, 451, 226−232. (37) Lindorff-Larsen, K.; Maragakis, P.; Piana, S.; Eastwood, M. P.; Dror, R. O.; Shaw, D. E. PLoS One 2012, 7, No. e32131. (38) Beauchamp, K. A.; Lin, Y.-S.; Das, R.; Pande, V. S. J. Chem. Theory Comput. 2012, 8, 1409−1414.
ACKNOWLEDGMENTS We are grateful for computational resources provided by the Danish Center for Scientific Computing. J.K. thanks the Danish Councils for Independent Research (the Sapere Aude programme), the Lundbeck Foundation, and the Villum foundation for financial support.
■
REFERENCES
(1) Barrett, P. J.; Chen, J.; Cho, M.-K.; Kim, J.-H.; Lu, Z.; Mathew, S.; Peng, D.; Song, Y.; Van Horn, W. D.; Zhuang, T.; Sönnichsen, F. D.; Sanders, C. R. Biochemistry 2013, 52, 1303−1320. (2) Robustelli, P.; Stafford, K. A.; Palmer, A. G. J. Am. Chem. Soc. 2012, 134, 6365−6374. (3) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999−3094. (4) Olsen, J. M.; Aidas, K.; Kongsted, J. J. Chem. Theory Comput. 2010, 6, 3721−3734. (5) Schwabe, T.; Olsen, J. M. H.; Sneskov, K.; Kongsted, J.; Christiansen, O. J. Chem. Theory Comput. 2011, 7, 2209−2217. (6) Olsen, J. M. H.; Kongsted, J. Adv. Quantum Chem. 2011, 61, 107−143. (7) Olsen, J. M. H. Ph.D. thesis, University of Southern Denmark, Odense, Denmark, 2012; DOI 10.6084/m9.figshare.156852. (8) Flaig, D.; Beer, M.; Ochsenfeld, C. J. Chem. Theory Comput. 2012, 8, 2260−2271. (9) Helgaker, T.; Jaszunski, M.; Ruud, K. Chem. Rev. 1999, 99, 293. (10) Helgaker, T.; Jorgensen, P. J. Chem. Phys. 1991, 95, 2595−2601. (11) London, F. J. Phys. Radium 1937, 8, 397−409. (12) Wolinski, K.; Hinton, J. F.; Pulay, P. J. Am. Chem. Soc. 1990, 112, 8251−8260. (13) Kongsted, J.; Nielsen, C. B.; Mikkelsen, K. V.; Christiansen, O.; Ruud, K. J. Chem. Phys. 2007, 126, No. 034510. (14) Helgakar, T.; Jørgensen, P.; Olsen, J. Molecular ElectronicStructure Theory; Wiley: New York, 2000; pp 2−5. (15) Applequist, J.; Carl, J. R.; Fung, K.-K. J. Am. Chem. Soc. 1972, 94, 2952−2960. (16) Gao, B.; Thorvaldsen, A. J.; Ruud, K. Int. J. Quantum Chem. 2011, 111, 858−872. (17) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jansik, B.; Jensen, H. J. A.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V.; Salek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; Sylvester-Hvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; Ågren, H. WIREs Comput. Mol. Sci. 2013, DOI: 10.1002/wcms.1172. (18) Dalton, a Molecular Electronic Structure Program, development version, 2013; http://daltonprogram.org/; accessed October 14, 2013. (19) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; 988
dx.doi.org/10.1021/ct400880n | J. Chem. Theory Comput. 2014, 10, 981−988
