This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

THE JOURNAL OF CHEMICAL PHYSICS 140, 104107 (2014)

Calculation of nuclear spin-spin coupling constants using frozen density embedding Andreas W. Götz,1,a) Jochen Autschbach,2 and Lucas Visscher3,b) 1

San Diego Supercomputer Center, University of California San Diego, 9500 Gilman Dr MC 0505, La Jolla, California 92093-0505, USA 2 Department of Chemistry, University at Buffalo, State University of New York, Buffalo, New York 14260-3000, USA 3 Amsterdam Center for Multiscale Modeling (ACMM), VU University Amsterdam, Theoretical Chemistry, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands

(Received 5 December 2013; accepted 22 January 2014; published online 11 March 2014) We present a method for a subsystem-based calculation of indirect nuclear spin-spin coupling tensors within the framework of current-spin-density-functional theory. Our approach is based on the frozen-density embedding scheme within density-functional theory and extends a previously reported subsystem-based approach for the calculation of nuclear magnetic resonance shielding tensors to magnetic fields which couple not only to orbital but also spin degrees of freedom. This leads to a formulation in which the electron density, the induced paramagnetic current, and the induced spinmagnetization density are calculated separately for the individual subsystems. This is particularly useful for the inclusion of environmental effects in the calculation of nuclear spin-spin coupling constants. Neglecting the induced paramagnetic current and spin-magnetization density in the environment due to the magnetic moments of the coupled nuclei leads to a very efficient method in which the computationally expensive response calculation has to be performed only for the subsystem of interest. We show that this approach leads to very good results for the calculation of solvent-induced shifts of nuclear spin-spin coupling constants in hydrogen-bonded systems. Also for systems with stronger interactions, frozen-density embedding performs remarkably well, given the approximate nature of currently available functionals for the non-additive kinetic energy. As an example we show results for methylmercury halides which exhibit an exceptionally large shift of the one-bond coupling constants between 199 Hg and 13 C upon coordination of dimethylsulfoxide solvent molecules. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4864053] I. INTRODUCTION

Nuclear magnetic resonance (NMR) spectroscopy is an indispensable tool in chemistry and biochemistry. It gives detailed insight into the electronic and geometric structure of molecules. The theoretical prediction of NMR parameters such as nuclear shielding tensors and indirect nuclear spinspin coupling tensors from first principles quantum chemical theory has become a routine task and is frequently used to assist in the assignment and interpretation of NMR spectra.1–4 Due to its favorable ratio of computational effort and accuracy, density-functional theory (DFT)5 is usually the method of choice for the calculation of NMR parameters. Approaches developed for DFT calculations of nuclear spin-spin couplings6–11 can be extended such that the required computation time shows a linear scaling with the system size.12 Nevertheless, the computational requirements for a quantum mechanical treatment of large systems remain high. This often prevents the investigation of systems with several hundreds or thousands of atoms such as biological macromolecules, metallo-enzymes, or molecules in solution or other complex environments. The importance of relativistic effects for heava) Electronic mail: [email protected] b) Electronic mail: [email protected]

0021-9606/2014/140(10)/104107/14/$30.00

ier nuclei13 may further increase the need for computational resources. From a conceptual point of view, however, it is neither desirable nor necessary to perform a full quantum chemical calculation for a large system if the investigated property is due to a perturbation that is very localized or probed only at a local position,14 as is the case for NMR parameters. The nuclear spin-spin coupling tensor describes the interaction of the magnetic moment of a nucleus with the induced current and spin-magnetization density of the electronic system due to the magnetic moment of a coupled nucleus. The operator corresponding to this nuclear magnetic moment decays quadratically with the distance from the nucleus and is therefore relatively short ranged. This readily calls for a subsystem based approach in which one focuses the computational effort for the property calculation on a subsystem which contains the coupled nuclei and thus avoids the quantum mechanical treatment of the full system. In the case of solvated molecules this can be achieved with continuum solvent models in which the solvent environment of the solute is described as a continuous medium characterized by its dielectric constant.15–17 Such an approximate treatment is able to describe unspecific, bulk solvent effects properly. In the presence of specific interactions like hydrogen bonding, however, the continuum description has to be combined with the explicit inclusion of a number

140, 104107-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-2

Götz, Autschbach, and Visscher

of solvent molecules.18–21 Recently, it was shown that nuclear spin-spin coupling constants can be computed for entire protein backbones with good accuracy using mixed quantum mechanics/molecular mechanics (QM/MM) with automatic fragmentation into subunits (AF-QM/MM).22 For the general description of complex environments, however, no practical method is available at present. Frozen-density embedding (FDE),23–25 which is based on a subsystem formulation of DFT,26, 27 provides a general scheme for the inclusion of environmental effects on molecular properties. In FDE, a system is partitioned into molecular subunits which are each treated by the orbital-based Kohn– Sham (KS) formalism, while the interactions between the subunits are taken into account via an embedding potential derived from orbital-free DFT. Hence, it becomes possible to focus the computational effort for a property calculation on the subsystem of interest while still keeping a quantum mechanical description of the full system. FDE is in principle an exact theory. Practical implementations, however, require not only an approximate exchange-correlation functional but also an approximate kinetic-energy functional for the construction of the embedding potential. Sufficiently accurate approximations are available which make it possible to study a large variety of systems if a judicious choice of subsystems is made.28 With this work we extend the FDE formalism to the calculation of nuclear spin-spin coupling tensors. We test the performance of the method for the calculation of solvent shifts of NMR J coupling constants using hydrogen-bonded dimers, a cluster of 17 water molecules, and dimethylsulfoxide (DMSO) complexes of dimethylmercury and methylmercury halides. In particular the latter systems represent a challenge for FDE because of the strong interaction of the coordinating solvent with mercury which is accompanied by an exceptionally large solvent shift, as shown previously by molecular dynamics (MD) calculations with explicit solvation at the quantum mechanical level.29 This work is organized as follows: In Sec. II we introduce the theory upon which the subsystem-based calculation of nuclear spin-spin coupling tensors within FDE is based. We first review current-spin-density functional theory (CSDFT) and the CSDFT-based calculation of nuclear spin-spin coupling tensors in Secs. II A and II B, respectively. In Sec. II C we extend the previously reported theory of FDE for systems in homogeneous magnetic fields30 to magnetic fields which couple not only to orbital, but also spin degrees of freedom. This formalism is then used in Sec. II D to derive workingequations for a subsystem-based calculation of nuclear spinspin coupling tensors within FDE. Computational details on the calculations performed in this work are given in Sec. III. Numerical results of our new approach along with a discussion are presented in Sec. IV and Sec. V summarizes with concluding remarks.

II. THEORY A. Current-spin-density-functional theory

The calculation of indirect nuclear spin-spin coupling tensors requires the generalization of DFT to include mag-

J. Chem. Phys. 140, 104107 (2014)

netic fields which couple to both spin and orbital degrees of freedom. For brevity, our discussion does not cover hybrid density functionals that include Hartree–Fock type exact exchange. Following Vignale and Rasolt,31 the ground state energy of a (nonrelativistic) system of interacting electrons in the presence of an external magnetic field B(r) = ∇ × A(r) can be written as a functional of the electron density ρ(r), the spin-magnetization density m(r), and the paramagnetic current density jp (r) = jorb (r) − ρ(r)A(r). The total physical current density j(r) = jorb (r) + jspin (r) is the sum of the total (gauge-invariant) orbital current density jorb (r) and the spin-current density jspin (r) = −∇ × m(r). We are using SIbased Hartree atomic units throughout this paper and make use of the Coulomb gauge for the external vector potential (∇ · A = 0). The energy functional is then given by E[ρ, m, jp ] = Ts [ρ, m, jp ] + vext (r)ρ(r) dr − 1 + 2 +

m(r) · B(r) dr +

1 2

jp (r) · A(r) dr

ρ(r)A2 (r) dr

ρ(r)ρ(r ) dr dr + Exc [ρ, m, jp ]. |r − r |

(1)

The total energy contains terms for the non-interacting kinetic energy, the interaction of the electron density with the scalar external potential, the interaction of the current and the electron density with the external vector potential, and the interaction of the spin-magnetization density with the corresponding external magnetic field, the Coulomb repulsion of the electrons, and the exchange-correlation (XC) energy, which now also depends on the paramagnetic current and the spinmagnetization density. The non-interacting kinetic energy Ts [ρ, m, jp ] is the kinetic energy of a reference system of non-interacting electrons which has the same electron density ρ, spin-magnetization density m and paramagnetic current jp as the interacting system. For such a system, the wave function is given by a single Slater determinant built from all occupied two-component single particle spinor orbitals ψ i . The electron density is given by ρ(r) =

occ

ψi∗ (r)ψi (r).

(2)

i

The spin-magnetization density is 1 ∗ m(r) = − ψ (r)σ ψi (r), 2 i i occ

(3)

where σ is the vector of Pauli matrices. The paramagnetic current density is i ∗ ψi (r)∇ψi (r) − ∇ψi∗ (r) ψi (r) . 2 i occ

jp (r) = −

(4)

Minimization of the total energy functional (1) with respect to the basic variables leads to the KS equations for the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-3

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

determination of the KS orbitals which read 1 1 (−i∇ + As (r))2 + vs (r) + σ · Bs (r) ψj (r) = εj ψj (r). 2 2 (5) The effective potentials vs , Bs , and As are such that the ground state densities ρ, m, and jp of the interacting system are reproduced. These effective potentials are given by vs (r) =vext (r) + vH (r) + vxc (r) +

1 2 A (r) − A2s (r) , 2

Bs (r) =B(r) + Bxc (r),

(6a)

As (r) = A(r) + Axc (r),

(6b)

where vext is the external electrostatic potential and the Hartree potential is given by ρ(r ) (7) vH (r) = dr . |r − r | The XC potentials are given as functional derivatives of the XC energy with respect to the corresponding conjugate densities, vxc (r) =

δExc [ρ, m, jp ] , δρ(r)

Bxc (r) = − Axc (r) =

δExc [ρ, m, jp ] , δm(r)

δExc [ρ, m, jp ] . δjp (r)

(8a) (8b) (8c)

The one-electron equations (5) are gauge invariant and the density and current obtained from them satisfy the continuity equation ∂ρ(r) = 0. (9) ∂t For the calculation of NMR parameters, the dependence of the XC functional on the paramagnetic current jp is usually neglected. In this case, the KS equations reduce to 1 1 2 u (−i∇ + A(r)) + vs (r) + σ · Bs (r) ψj (r) = εj ψj (r). 2 2 (10) The effective potential vsu is known from standard KS theory, ∇ · j(r) +

vsu (r) = vext (r) + vH (r) + vxc (r).

(11)

Since no XC functionals are available that explicitly depend on the magnetization density m, the current-dependence arising from the electron spin is treated with standard spin-density functionals. This can be done in the collinear or non-collinear approach;32 the latter was used for spin-spin coupling constant calculations in a two-component relativistic non-hybrid and hybrid DFT framework.11, 33 B. DFT calculation of nuclear spin-spin coupling constants

To lowest order, the reduced nuclear spin-spin coupling tensor K(A, B) involving two nuclei A and B with nuclear

magnetic moments μA = γA IA and μB = γB IB is given as mixed second derivative1–4, 34 of the total energy E,

d 2 E(μA , μB )

Kst (A, B) = , (12) dμA,s dμB,t μA =0,μB =0 where γ A and γ B are the nuclear gyromagnetic ratios and IA and IB the nuclear spin angular moments. The indices s and t label individual Cartesian components. For freely tumbling molecules only the isotropic part K(A, B) = (1/3)TrK(A, B) of this matrix is observed. It is connected to the usually reported indirect coupling constants J(A, B) via ¯ (13) γA γB K(A, B). 2π The vector potential of a magnetic point-nucleus A is given as J (A, B) =

AA (r) =

1 μA × rA , c2 rA3

(14)

where we have used rA = r − RA and RA is the position of the nucleus and c is the speed of light. (For a finite-nucleus formalism see Ref. 35.) Substituting the vector potential A = AA + AB and associated magnetic field B = BA + BB due to the nuclear magnetic moments μA and μB in the energy expression of Eq. (1) and evaluating the derivative of Eq. (12) leads to the expression rA × j(μB,t ) (r) 1 dr (15) Kst (A, B) = 2 c rA3 s for the coupling tensor. It describes how the total (first-order) current induced by the magnetic moment of one nucleus interacts with the magnetic moment of the other nucleus. The induced current is given as

∂j(r)

(μB,t ) (r) = (16) j ∂μB,t μB =0 and can be decomposed into an orbital and a spin part. According to (μ

)

B,t ) (r) + ρ(r) jorbB,t (r) = j(μ p

∂AB (r) , ∂μB,t

(17)

the first-order induced orbital current can be further decomposed into a paramagnetic and a diamagnetic part. The latter is responsible for the diamagnetic spin-orbit (DSO) part of the coupling tensor which is given as (μA,s ,μB,t ) ρ(r) dr, (18) KstDSO (A, B) = hDSO with the DSO operator (μ

A,s hDSO

,μB,t )

=

1 (rA · rB )δst − rB,s rA,t . c4 rA3 rB3

(19)

The DSO term can directly be computed from the unperturbed ground state electron density. The paramagnetic spin-orbit (PSO) part of the coupling tensor is given as (μ ) 1 rA × jp B,t (r) PSO dr (20) Kst (A, B) = 2 c rA3 s

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-4

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

and requires the determination of the first-order induced para(μ ) magnetic current jp B,t = ∂jp /∂μB,t . The sum of the so-called Fermi-contact (FC) and spin-dipole (SD) contributions is given as (μ ) rA × jspinB,t (r) 1 FC+SD (A, B) = 2 dr. (21) Kst c rA3 s

It requires the determination of the induced spin-current den(μ ) sity jspinB,t = ∂jspin /∂μB,t which in turn requires knowledge of the induced spin-magnetization density m(μB,t ) = ∂m/∂μB,t . The induced orbital current and spin-magnetization density can be determined from the response of the unperturbed KS orbitals to the magnetic field generated by the nuclear magnetic moment μB . The KS equations in the presence of a magnetic field and in the approximation of uncoupled DFT are given by Eq. (10). The first-order perturbed KS orbitals (μ ) ψi B,t = ∂ψi /∂μB,t are thus obtained from the first-order KS equations (μ )

(0) (μ ) F − εi(0) ψi B,t + F (μB,t ) − εi B,t ψi(0) = 0, (22) where the superscript (0) has been used to refer to the KS operator, the KS orbitals and KS orbital eigenvalues of the unperturbed system while the first-order perturbed KS operator is given as

∂F (μB )

F (μB,t ) = . (23) ∂μB,t μB =0 To first order in the magnetic moment μB , the perturbing operator is given as a sum of three contributions which can be treated separately, (μ

)

(μ

)

(μ

B,t + hFCB,t + hSDB,t h(μB,t ) = hPSO

with (μ

)

B,t =− hPSO

(μ

i c2

rB ×∇ rB3

(μ

,

(25a)

t

4π δ(rB )σt , 3c2 (σ · rB )rB,t 1 σt . = 2 3 − 2c rB5 rB3

)

)

(24)

hFCB,t = hSDB,t

)

(25b) (25c)

The first-order perturbed orbitals are usually determined by expanding them in terms of the virtual canonical KS orbitals of the unperturbed system as (μB,t )

ψi

(r) =

virt

(μ

)

uia B,t ψa(0) (r)

(26)

a

and it can easily be seen from Eq. (22) that the coefficients (μ ) uia B,t are determined from (0) (μ ) (0) ψ F B,t ψi (μB,t ) . (27) uia = a (0) εi − εa(0) The PSO perturbation operator is purely imaginary, which implies that the first-order perturbed orbitals will also be purely imaginary and that the first-order change in the electron density and spin-magnetization density vanishes. Thus,

in the approximation of uncoupled DFT, the first-order per(μ ,t) turbed KS operator FPSOB is just the perturbing operator (μB,t ) (μB,t ) can directly be obtained hPSO and the coefficients uia,PSO from Eq. (27). Choosing the unperturbed KS orbitals as real so that the first-order perturbed orbitals are purely imaginary, the first-order induced paramagnetic current can be expressed in terms of KS orbitals as B,t ) (r) j(μ p

occ (0) (μB,t ) (μB,t ) ψi (r)∇ψi,PSO = −i (r) − ψi,PSO (r)∇ψi(0) (r) . i

(28) The FC and SD perturbation operators, on the other hand, induce a first-order change in the spin-magnetization density, but not in the total electron density. As a consequence the per(μ ) (μ ) turbed KS operator FFCB,t or FSDB,t contains, in addition to (μB,t ) (μ ) the perturbing operator hFC or hSDB,t , a contribution from the first-order induced XC potential δBxc (r) (μB,t ) (μB,t ) m (r ) dr . (29) Bxc (r) = δm(r ) The first-order induced magnetization density, expressed in terms of (real) KS orbitals, becomes m(μB,t ) (r) = −

occ

(μ

)

B,t ψi(0) σ ψi,FC+SD

(30)

i

and it is clear that the expansion coefficients from Eq. (27) have to be determined in an iterative fashion. It is noted that the formalism summarized thus far is based on a nonrelativistic KS Hamiltonian. For relativistic approaches, see Ref. 13. With the zeroth-order regular approximation (ZORA) relativistic method36–38 used for some of the calculations in this work, the formalism appears very similar to the nonrelativistic framework but employs suitably modified operators and complex two-component KS orbitals.10, 11, 33 In particular in a spin-free (scalar) variant, the ZORA-based perturbation equations closely resemble those presented here, but with non-singular replacements of the operators in Eq. (25). Functional-specific details in the following discussion of FDE are straightforwardly transferred to the relativistic domain as long as nonrelativistic spin-density functionals are employed (as it is commonplace). C. Frozen-density embedding for systems in external magnetic fields

The original formulation of FDE23, 24 within DFT is based upon a division of the total system into two subsystems with electron densities ρ1 (r) and ρ2 (r) such that the total electron density ρ(r) is obtained as the sum of these subsystem electron densities. The subsystem electron densities are determined separately from a set of one-electron equations in which the effect of the density in the other subsystem is represented in terms of an effective embedding potential. This formalism has been extended30 to include external, homogeneous magnetic fields in the framework of current-DFT, which requires in addition a partitioning of the paramagnetic current density jp (r). In the case of a general external magnetic field, which couples not only to orbital, but also to spin

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-5

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

degrees of freedom, also the total spin-magnetization density m(r) needs to be partitioned. We therefore have ρ(r) = ρ1 (r) + ρ2 (r),

(31a)

jp (r) = jp,1 (r) + jp,2 (r),

(31b)

m(r) = m1 (r) + m2 (r).

E[ρ1 , m1 , jp,1 , ρ2 , m2 , jp,2 ]

+ Tsnad [ρ1 , m1 , jp,1 , ρ2 , m2 , jp,2 ] + − +

1 2

1 + 2

vext (r)ρ(r) dr

m(r) · B(r) dr +

jp (r) · A(r) dr

ρ(r)A2 (r) dr

ρ(r)ρ(r ) dr dr + Exc [ρ, m, jp ]. |r − r |

(32)

The non-additive kinetic-energy functional is defined as Tsnad [ρ1 , m1 , jp,1 , ρ2 , m2 , jp,2 ] = Ts [ρ, m, jp ] − Ts [ρ1 , m1 , jp,1 ] − Ts [ρ2 , m2 , jp,2 ]. (33) The non-interacting kinetic energies Ts [ρ1 , m1 , jp,1 ] and Ts [ρ2 , m2 , jp,2 ] of the two subsystems are the kinetic energies of reference systems of non-interacting electrons which have electron densities, spin-magnetization densities and paramagnetic currents such that Eq. (31) for the total, interacting system holds. The subsystem densities and currents can therefore be represented in terms of KS orbitals of the individual subsystems and the corresponding non-interacting kinetic energies can be evaluated exactly from these orbitals. The KS orbitals of the total system, however, are not known. Therefore, the total kinetic energy cannot be evaluated in the standard way. In practical applications of the FDE scheme, one has to resort to an approximate, orbital-independent kineticenergy functional to evaluate Tsnad . Minimization of the total energy functional (32) with respect to the basic variables of one subsystem leads, for given basic variables of the other subsystem (dubbed as frozen subsystem or embedding environment), to one-particle equations for the determination of the KS orbitals which read 2 1 1 −i∇ + As,1 (r) + vs,1 (r) + σ · Bs,1 (r) ψi,1 (r) 2 2 = εi,1 ψi,1 (r).

nad (r) + vs,1 (r) = vsu (r) + vT,1

1 2 A (r) − A2s,1 (r) , 2

As,1 (r) = As (r) + Anad T,1 (r),

(35a)

Bs,1 (r) = Bs (r) + Bnad T,1 (r),

(35b)

(31c)

With these definitions, the total current is given as the sum of the currents in the two subsystems, j(r) = j1 (r) + j2 (r). Furthermore, if the continuity equation is satisfied for each of the two subsystems individually, it is also satisfied for the total system. The FDE-CSDFT total energy can now be expressed as a functional of the electron densities, paramagnetic current densities and spin-magnetization densities of the subsystems,

= Ts [ρ1 , m1 , jp,1 ] + Ts [ρ2 , m2 , jp,2 ]

In the spirit of earlier works on FDE, these equations will be called Kohn–Sham equations with constrained electron density (KSCED). The effective potentials in the KSCED equations are given by

(34)

and differ from the standard effective potentials of CSDFT by contributions due to the non-additive kinetic energy nad vT,1 (r) =

δTsnad [ρ1 , m1 , jp,1 , ρ2 , m2 , jp,2 ] , δρ1 (r)

(36a)

Anad T,1 (r) =

δTsnad [ρ1 , m1 , jp,1 , ρ2 , m2 , jp,2 ] , δjp,1 (r)

(36b)

Bnad T,1 (r) =

δTsnad [ρ1 , m1 , jp,1 , ρ2 , m2 , jp,2 ] . δm1 (r)

(36c)

In general, it will not be possible to find a frozen environment density such that ρ − ρ 2 is positive. As a consequence the solution of Eq. (34) will not yield the same electron density, spin-magnetization density, and paramagnetic current density as the solution of the standard CSDFT equations (5), even if the exact non-additive kinetic-energy functional would be used. In this case, the subsystem densities and currents can be iteratively refined in so-called “freezeand-thaw”39 (FT) cycles by exchanging the role of the frozen and non-frozen subsystems, that is, by solving two coupled sets of KSCED equations (34) for subsystems 1 and 2. The currently available approximations for the nonadditive kinetic-energy functional do not depend on the paramagnetic current. It has already been argued30 that for the calculation of NMR shielding parameters in weakly interacting systems, where these approximations are applicable, the neglect of the current dependence is likely to be smaller than the intrinsic error of the approximate functionals. Neglecting also the current dependence of the XC functional one finally arrives at the uncoupled KSCED equations 1 1 nad (−i∇+A(r))2 +vsu (r)+vT,1 (r)+ σ · Bs (r) + Bnad (r) T,1 2 2 · ψi,1 (r) = εi,1 ψi,1 (r).

(37)

Just as for the XC energy, the available functionals for the non-additive kinetic energy do not explicitly depend on the magnetization density m and spin-density functionals have to be employed in the collinear or non-collinear approach. As a consequence of the neglect of the current dependence of Exc and Tsnad , the KSCED equations of the two subsystems are coupled only via the electron densities and spinmagnetization densities, but not via the paramagnetic current densities. Calculations on the non-frozen subsystem can thus be carried out without knowledge of the induced paramagnetic current jp,2 in the frozen subsystem.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-6

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

D. Calculation of nuclear spin-spin coupling tensors with frozen-density embedding

The FDE formalism described in Sec. II C can be applied to the calculation of nuclear spin-spin coupling tensors by decomposing the first-order current density induced by the magnetic field of a nucleus into contributions of two subsystems, (μB,t )

j(μB,t ) (r) = j1

(μB,t )

(r) + j2

(r).

(38)

Kst (A, B) = Kst,1 (A, B) + Kst,2 (A, B),

1 Kst,n (A, B) = 2 c

(μ

)

rA × jn B,t (r) rA3

nad,(μB,t )

BT,n

=

(43a)

(r)

δBnad T,n [ρ1 , m1 , ρ2 , m2 ](r) (μB,t ) mn (r ) dr , δmn (r )

(n = 1, 2). (43b)

Then also the coupling tensor, Eq. (15), can be expressed in terms of subsystem contributions as

with

additive kinetic-energy potentials δBxc,n [ρ, m](r) (μB,t ) (μB,t ) mn (r ) dr , (r) = Bxc,n δmn (r )

(39)

dr,

(n = 1, 2).

s

(40) These contributions can again be split up into DSO, PSO, FC, and SD terms. For their determination, one proceeds as follows. First, the unperturbed ground state densities of the two subsystems are obtained from the KSCED equations, possibly applying the FT procedure. The DSO contributions of the two subsystems to the coupling tensor can then directly be evaluated from (μA,s ,μB,t ) DSO (A, B) = hDSO ρn (r) dr, (n = 1, 2) (41) Kst,n since they depend only on the unperturbed subsystem electron densities. In order to calculate the PSO contribution to the coupling tensor, the induced paramagnetic current in each of the two subsystems needs to be determined. In the approximation of uncoupled DFT, the first-order perturbed KS operator (μB ,t) for a subsystem n is just given by the PSO operator of FPSO,n Eq. (25a). It is independent of the densities and currents in the subsystems. Thus, there is no coupling between the subsystems and, using Eq. (27), the first-order perturbed KS orbitals (μB,t ) (μ ) and paramagnetic currents jp,nB,t are obtained indeψi,PSO,n pendently for each subsystem from its unperturbed ground state KS orbitals and orbital eigenvalues. The PSO contributions of the two subsystems to the coupling tensor can then be evaluated from (μ ) 1 rA × jp,nB,t (r) PSO dr, (n = 1, 2). Kst,n (A, B) = 2 c rA3 s (42) The FC and SD contributions to the coupling tensor require the determination of the induced spin-magnetization (μ ) densities mn B,t in the subsystems, which is a bit more involved than the determination of the induced paramagnetic current densities in the uncoupled approximation. The KSCED equations (37) of the two subsystems are coupled via the spin-magnetization densities in the subsystems. Thus, in addition to the FC and SD operators of Eqs. (25b) and (25c), (μB ,t) (μB ,t) and FSD,n conthe first-order perturbed KS operators FFC,n tain contributions from the first-order induced XC and non-

Since these induced potentials depend on the induced spin(μ ) magnetization densities mn B,t of both subsystems, the determination of the latter needs to proceed via an iterative procedure. Once the first-order spin-magnetization densities induced in the subsystems have been determined, the FC and SD contributions to the coupling tensor are obtained from (μB,t ) rA × jspin,n (r) 1 FC+SD dr, (n = 1, 2), Kst,n (A, B) = 2 c rA3 s (44) (μB,t ) (μ ) where jspin,n = −∇ × mn B,t are the induced spin-currents in the subsystems. The current density which is induced in the electronic system by the magnetic moment of a nucleus is probed by the nuclear spin of the coupled nucleus. As can be seen, for example, from Eqs. (15) and (40), only the induced current in the vicinity of the coupled nucleus contributes to the coupling tensor. Just as is the case for NMR shielding tensors,30 the subsystem-based formalism for the calculation of nuclear spin-spin coupling tensors presented above readily allows to exploit this near-sightedness. With an appropriate partitioning in which both coupled nuclei are located in the same subsystem, the major contribution to the coupling tensor is due to the induced current of the subsystem that contains the coupled nuclei. It then becomes possible to neglect the induced current in the other subsystem and its electron and spin density can be kept frozen for the calculation of the coupling tensor. At this point it also becomes possible to employ additional approximations for the construction of the frozen density,25 for instance, by using a sum-of-fragments density, or by applying the FT procedure only to parts of the system that are close to the subsystem of interest. III. COMPUTATIONAL DETAILS

All calculations have been performed with a development version of the Amsterdam Density Functional (ADF)25, 40–42 program package. Nuclear spin-spin coupling tensors have been calculated with the CPL program10, 11, 33 which is part of this package and which we have modified for FDE as described in Sec. III A. To set up and execute all calculations as well as to retrieve the data from the output files, we have used PyADF43 which is a scripting framework for quantum chemistry implemented in the Python programming language. The ZORA scalar relativistic method36–38 has been employed for the mercury compounds. The (non-relativistic) gradient-corrected PBE44, 45 XC functional was applied to determine the unperturbed molecular orbitals while the local

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-7

Götz, Autschbach, and Visscher

spin density approximation (LSDA) in the parametrization V of Vosko, Wilk and Nusair (VWN)46 was used for the first-order perturbed XC potential due to the perturbation operators. Additional results using the PBE functional for the first-order perturbed XC potential are collected in the supplemental material.47 The FDE calculations have been done with the Thomas–Fermi functional48, 49 and with the gradientcorrected PW91k50 functional for the non-additive kinetic energy. In particular, the latter has been shown to yield reliable results for different electronic properties in weakly interacting systems including van der Waals51–54 and hydrogenbonded30, 54–58 ones. In consistency with the use of the LSDA for the first-order perturbed XC potential, the LSDA was employed for the first-order perturbed non-additive kinetic energy potential (see Sec. III A). For the mercury compounds, the ZORA analogs of the FC, DSO, and PSO operators have been included in the response calculations. The ZORA analog of the SD operator has not been taken into account because its contribution to the J coupling constant has been shown to be negligible for these compounds.59 In non-relativistic calculations all relevant operators have been included. The default point charge model was used for the nuclear charge distribution and the calculations were carried out with the TZP basis of the ADF basis set library40 which is a tripleζ valence/double-ζ core all-electron Slater basis augmented with one set of polarization functions (6p for Hg, 2p for H, nd for C, Cl, Br, I). The basis of hydrogen and first- and second-row elements was augmented with tight s-functions60 while the basis of mercury was augmented with tight sand p-functions,10 which is essential to obtain accurate results for the FC contribution to the nuclear spin-spin couplings. The FDE calculations have been performed both with a monomolecular expansion basis, denoted as FDE(m), in which only basis functions located on the non-frozen subsystem are used for the expansion of its KS orbitals, and a supermolecular (global) expansion basis, denoted as FDE(s), in which the basis functions of both subsystems are used for the expansion of the KS orbitals of the non-frozen subsystem. From the point of view of computational efficiency FDE(m) is clearly superior to FDE(s). However, while FDE(m) allows for polarization of the subsystem electron densities, no charge transfer between the subsystems can be described with FDE(m) and more accurate results should be expected from FDE(s). In FDE(m) calculations, the numerical integration grid was restricted to the active subsystem.57 A. Implementation notes

The calculation of nuclear spin-spin coupling tensors using FDE does not require major modifications to existing NMR couplings programs. For the DSO and PSO terms, Eqs. (41) and (42), it merely has to be ensured that the electron density and KS orbitals obtained from a FDE calculation can be used in the program. Some work is required for the FC and SD terms, Eq. (44). First, the frozen electron density of the embedding environment needs to be known as well. The first-order perturbed XC potential of Eq. (43a), which is required to determine the first-order induced spin-magnetization density, de-

J. Chem. Phys. 140, 104107 (2014)

pends on the total electron density of the system. Next, the non-additive kinetic energy contribution to the first-order induced potential, Eq. (43b), needs to be evaluated. In order to do so, conventional spin-density functionals can be used with the collinear or non-collinear approach for the definition of the spin densities,32 just as is done for the XC potential.11 For this work we have adopted the LSDA, that is we make use of TTF [ρσ ], (45) Ts [ρ↑ , ρ↓ ] ≈ 22/3 σ

where σ = ↑, ↓ labels the spin and TTF [ρ] = CTF ρ 5/3 (r) dr

(46)

is the Thomas–Fermi functional48, 49 with Thomas–Fermi constant CTF = (3/10)(3π 2 )2/3 . The kernel of the non-additive kinetic energy which determines the induced potential then takes the form δ 2 Tsnad [ρ↑,1 , ρ↓,1 , ρ↑,2 , ρ↓,2 ] δρσ,n (r)δρσ ,n (r ) ≈ 22/3 δσ σ δ(r − r )

δ 2 Ts [ρ]

δ 2 Ts [ρ]

· − , δρ 2 ρ=ρσ (r) δρ 2 ρ=ρσ,n (r)

(n = 1, 2), (47)

where δ 2 Ts [ρ] 10 = CTF ρ −1/3 (r). 2 δρ 9

(48)

Finally, to calculate the contribution to the couplings due to induced currents in the embedding environment, the nuclear positions of the frozen subsystem need to be known. For obvious reasons, one will always choose a partitioning such that the coupled nuclei are in the same subsystem. Then, however, as argued at the end of Sec. II D, the response calculation has to be performed only for the subsystem which contains these nuclei. B. Solvent induced shift of J coupling constants

The change of an intramolecular J coupling constant in a molecule P upon complexation with a molecule Q is obtained as the difference between the coupling constant in the complex PQ and the coupling constant in the monomer P with the deformed complex geometry, Q J KS = JPKS,P (P Q) − JPKS,P Q (P ). Q

(49)

Here, the super- and subscript denote the basis set and geometry, respectively. The system for which the calculation is performed is given in parenthesis. The above given expression contains a basis set superposition error (BSSE), KS,P Q JPBSSE (P ) − JPKS,P Q (P ) = JP Q Q (P )

(50)

and the BSSE corrected solvent shift of the J coupling constant is given as Q KS JBSSE = JPKS,P (P Q) − JPBSSE Q (P ) Q Q Q = JPKS,P (P Q) − JPKS,P (P ). Q Q

(51)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-8

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

(a)

(b)

(c)

FIG. 1. Structures of hydrogen bonded dimers used in the J coupling calculations: (a) hydrogen fluoride dimer, (b) ammonia–water dimer, (c) water dimer. 1 J(F,H), 1 J(O,H), and 1 J(N,H) couplings were calculated for the hydrogen bond donor subsystems (for the hydrogen atom that is involved in the hydrogen bond) and for the hydrogen bond acceptor subsystems (in the case of ammonia for the hydrogen atom in the symmetry plane). For the FDE calculations the dimers have been partitioned into donor (left molecule) and acceptor (right molecule) subsystems.

In the case of FDE calculations with a supermolecular expansion basis (denoted as FDE(s)) the solvent shift is given in analogy to Eq. (49) as Q (P Q) − JPKS,P J FDE(s) = JPFDE,P Q (P ) Q

(52)

and the BSSE corrected solvent shift as FDE(s) Q Q JBSSE = JPFDE,P (P Q) − JPKS,P (P ). Q Q

(53)

If a monomolecular expansion basis is employed in the FDE calculations (denoted as FDE(m)), there is no BSSE because the basis set employed is the same in the embedding calculation and the reference calculation. The complexation induced shift in the J couplings obtained with FDE(m) is then given as (P Q) − JPKS,P J FDE(m) = JPFDE,P Q Q (P ).

(54)

The shift of the J coupling constants defined in this way does not include effects due to the structural changes upon complexation (dimerization). IV. RESULTS A. Hydrogen-bonded dimers

As a first test we have investigated the ability of FDE to reproduce the shift of J couplings upon hydrogen bond formation in dimers of small molecules. For this purpose we have computed the 1 J(O,H), 1 J(N,H), and 1 J(F,H) couplings in the

hydrogen bond donor and hydrogen acceptor subsystems of the ammonia–water dimer, the water dimer, and the hydrogen fluoride dimer, respectively. The structures of the dimers have been taken from the HB6/04 dataset61 and are depicted in Figure 1. The results are collected in Tables I and II.47 The shift in the J coupling constants upon hydrogen bond formation is largest in the hydrogen fluoride dimer while it is of comparable size in the ammonia–water dimer and the water dimer. The results obtained with FDE agree remarkably well with the KS reference results. There is no strong dependence of the values on the non-additive kinetic energy functional, although in most cases the PW91k functional leads to slightly better results than the Thomas–Fermi functional. One freeze-thaw cycle is sufficient to converge the J couplings. For the ammonia–water dimer and the water-dimer FDE(m,0) and FDE(m,3) yield similar results. Already a simple embedding is able to capture the majority of the solvent shift and polarization with freeze-thaw cycles does not have a pronounced effect on the results. Similarly, the use of a supermolecular expansion basis instead of a monomolecular expansion basis does not have a significant effect on those results. However, FDE(s,3) leads to overpolarization of the J couplings in the hydrogen bond donor subsystems in all cases as compared to FDE(m,3). The situation is different for the hydrogen fluoride dimer. Due to significant charge transfer upon dimer formation, the results improve significantly both upon polarization

TABLE I. Shifts of the 1 J couplings in the hydrogen bond donor and hydrogen bond acceptor subsystems upon formation of small hydrogen-bonded dimers (see Figure 1). The coupling constant for the isolated hydrogen bond donor subsystem is 341.9 Hz, −59.2 Hz, and −60.4 Hz and the coupling constant for the isolated hydrogen bond acceptor subsystem is 347.9 Hz, 37.2 Hz, and −62.0 Hz for HF–HF, NH3 –H2 O, and H2 O–H2 O, respectively (at the dimer geometry). The expansion basis (monomolecular or supermolecular) and the number of FT cycles is given in parentheses (J couplings in Hz). HF–HF

NH3 –H2 O

H2 O–H2 O

donor 1 J(F,H)

acceptor 1 J(F,H)

donor 1 J(O,H)

acceptor 1 J(N,H)

donor 1 J(O,H)

acceptor 1 J(O,H)

KS

17.9

22.2

− 5.3

0.9

− 4.1

− 2.7

Thomas–Fermi FDE(m,0) FDE(m,1) FDE(m,3) FDE(s,3)

6.6 10.6 10.6 19.9

24.3 26.8 26.9 25.2

− 5.6 − 6.1 − 6.2 − 6.8

0.4 0.6 0.6 0.7

− 3.7 − 4.2 − 4.2 − 4.8

− 2.1 − 2.7 − 2.7 − 2.7

PW91k FDE(m,0) FDE(m,1) FDE(m,3) FDE(s,3)

5.7 10.0 10.1 19.8

25.9 27.9 27.9 26.2

− 5.0 − 5.7 − 5.7 − 6.4

0.6 0.8 0.8 0.8

− 3.2 − 3.7 − 3.7 − 4.4

− 2.6 − 3.0 − 3.0 − 3.0

Method

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-9

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

TABLE II. Solvent shift of the 1 J(H,F) coupling constant in the hydrogen bond donor system of the hydrogen fluoride dimer. The expansion basis (monomolecular or supermolecular) and the number of FT cycles is given in parentheses (J couplings in Hz). 1 J(F,H)

FC

SD

KS

42.0

− 1.8

Thomas–Fermi FDE(m,0) FDE(m,3) FDE(s,3)

30.2 37.3 42.3

PW91k FDE(m,0) FDE(m,3) FDE(s,3)

22.7 30.6 36.3

DSO

PSO

Total

0.7

− 23.0

17.9

− 0.8 − 1.0 0.0

− 0.2 − 0.3 − 0.3

− 22.7 − 25.3 − 22.2

6.6 10.6 19.9

0.3 − 0.1 1.0

− 0.2 − 0.2 − 0.2

− 17.2 − 20.2 − 17.2

5.7 10.0 19.8

with freeze-thaw cycles and upon use of a supermolecular expansion basis. This is, in particular, true for the J coupling constant in the hydrogen bond donor subsystem while freezethaw leads to slightly worse results in the hydrogen bond acceptor subsystem. Due to the large discrepancy in the shift 1 J(F,H) for the hydrogen bond donor subsystem of the hydrogen–fluoride dimer with FDE(m,0) and also FDE(m,3) it is instructive to break down the individual contributions from the different J coupling operators. Table II reveals that the shift 1 J(F,H) in the hydrogen fluoride dimer is mainly due to two large contributions with opposite sign, a large positive contribution from the FC term and a negative contribution of the same order of magnitude from the PSO term. For FDE(m,0) and FDE(m,3) the error in the total shift is quite large (around 70% and 40%, respectively). However, the relative error in the individual contributions is much smaller which is an encouraging result. FDE(s,3) with the Thomas–Fermi functional is very close to the reference KS results for all terms. The good result in the total shift with the PW91k functional, however, turns out to be partly due to error cancellation between the FC and the PSO term.

B. Water cluster

We now turn to the calculation of the solvent shift of the J coupling constants of a water molecule upon formation of a cluster of 17 water molecules. Such a cluster could, for example, have been extracted from a snapshot of a molecular dynamics or Monte Carlo simulation of water aimed at computing the NMR parameters of water at thermodynamic equilibrium in solution, as has been done for NMR shieldings.62 Here we study only a single cluster of 17 water molecules with coordinates taken from the work of Cybulski et al.63 as depicted in Figure 2. It can be expected that the majority of the solvent shift is due to the directly hydrogen bonded water molecules, in this case the two hydrogen bond acceptor molecules (red spheres in Figure 2) and the two hydrogen bond donor molecules (blue spheres in Figure 2). We thus have performed supermolecular KS-DFT calculations for clusters of three water molecules (including the hydrogen bond acceptors), five

FIG. 2. The shift of the J coupling constants of the central water molecule upon formation of a cluster of 17 water molecules was calculated. Cluster sizes of 3 (including the closest hydrogen bond acceptors, red spheres), 5 (additionally including the closest hydrogen bond donors, blue spheres), and 17 were investigated with KS-DFT. FDE calculations included only the central water molecule or the cluster of 3 and 5 water molecules in the active region and freeze-thaw was performed either for some of the closest water molecules or all water molecules.

water molecules (additionally including the hydrogen bond donors), and 17 water molecules (the entire cluster). FDE calculations were performed with different partitionings into active region, freeze-thaw region, and completely frozen region. We denote these partitionings with a triple of numbers [i/j/k] in square brackets to indicate the number of water molecules in the active region (i), the freeze-thaw region (j), and the frozen region (k). For instance, FDE[1/0/16](m,0) means that only the central water molecule is in the active region while the surrounding 16 water molecules are kept frozen at their gas phase densities. In a FDE[3/2/12](m,1) calculation, the active region consists of the central water molecule and the two closest hydrogen bond acceptors (red), while one freezethaw cycle is performed with the two closest hydrogen bond donors (blue) and the remaining 12 water molecules are kept frozen at their gas phase densities. The results obtained with the different approaches are listed in Table III. We discuss the results for the 1 J(O,H) coupling constant since the results for the other intramolecular coupling constants follow the same trend. The solvent shift for the cluster of 17 water molecules is 1 J(O,H) = −14.4 Hz. A large fraction of this shift (−9.1 Hz) is captured by the small cluster model of three water molecules that includes the two hydrogen bond acceptors. However, even the larger cluster that contains all directly hydrogen bonded water molecules leads to a solvent shift of only −11.9 Hz and thus remains far from the reference value. Conventional KS-DFT based calculations for a smaller cluster therefore would not be a reliable approximation. The computationally most efficient FDE approximation, FDE[1/0/16](m,0) leads to a solvent shift of 1 J(O,H) = −15.0 Hz that is suspiciously close to the reference value of −14.4 Hz. It includes only the central water molecule in the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-10

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

TABLE III. Shift of the J coupling constants of the central water molecule upon formation of a cluster of 17 water molecules (see Figure 2). The coupling constants for the isolated water molecule (at the cluster geometry) are 1 J(O,H) = −48.5 Hz, 1 J(O,H ) = −48.9 Hz, and 2 J(H,H ) = −4.1 Hz. The triple of numbers [i/j/k] in square brackets indicate the number of water molecules in the active region (i), the freeze-thaw region (j), and the frozen region (k). The expansion basis and the number of FT cycles are given in parentheses. FDE calculations were performed with the PW91k functional for the non-additive kinetic energy (J couplings in Hz). Method KS cluster of 3 KS cluster of 5 KS cluster of 17 FDE[1/0/16](m,0) FDE[1/4/12](m,1) FDE[1/4/12](m,3) FDE[1/16/0](m,1) FDE[1/16/0](m,3) FDE[3/0/14](m,0) FDE[3/2/12](m,1) FDE[3/2/12](m,3) FDE[3/14/0](m,1) FDE[3/14/0](m,3) FDE[5/0/12](m,0) FDE[5/12/0](m,1) FDE[5/12/0](m,3)

1 J(O,H)

1 J(O,H )

2 J(H,H )

− 9.1 − 11.9 − 14.4 − 15.0 − 18.4 − 18.7 − 19.6 − 20.3 − 13.6 − 14.0 − 14.1 − 14.9 − 14.8 − 13.6 − 14.2 − 14.3

− 6.9 − 10.9 − 14.4 − 14.2 − 17.8 − 18.0 − 19.2 − 19.9 − 11.7 − 13.7 − 13.7 − 14.4 − 14.8 − 13.1 − 14.1 − 14.2

− 0.2 − 0.3 − 0.5 − 1.3 − 1.6 − 1.6 − 1.6 − 1.6 − 0.5 − 0.7 − 0.7 − 0.8 − 0.8 − 0.4 − 0.5 − 0.5

active region and keeps the electron density of all surrounding water molecules frozen at their individual gas phase values. This good result clearly is due to fortuitous error cancellation and similar observations were made in the case of NMR shieldings.62 Relaxing the electron density of the surrounding subsystems leads to successively worse results with a value of −20.3 Hz for FDE[1/16/0](m,3), in which the electron density of all water molecules has been iteratively optimized with three freeze-thaw cycles. Clearly, with the PW91k functional for the non-additive kinetic energy, such a partitioning does not lead to satisfactory results. The FDE results can be significantly improved by extending the size of the active subsystem through inclusion of the hydrogen bond acceptor water molecules. We obtain a result of −13.6 Hz for simple embedding with FDE[3/0/14](m,0), that is an error of 0.8 Hz, while all results that relax the environment with freeze-thaw are within ±0.5 Hz of the reference value. For comparison, a supermolecular KS-DFT calculation with the cluster model consisting of 5 water molecules is in error by 2.5 Hz. Including all hydrogen bonded water molecules into the active subsystem, finally, brings FDE with freeze-thaw into quantitative agreement with the reference supermolecular KS-DFT results. A single freeze-thaw iteration is sufficient to converge the results. C. Solvated mercury complexes

Environmental effects on nuclear spin-spin coupling constants in complexes with open coordination sites such as dimethylmercury and methylmercury halides, square planar Pt complexes, and similar systems can be substantial and are often amplified by relativistic effects.29, 59, 64, 65 Solvent

molecules may coordinate to the heavy atom and in order to achieve agreement of computed J couplings with experimental data obtained from solution an explicit treatment of the solvent coordination may become necessary.59 An analysis of the solvent coordination effect in terms of individual orbital contributions has shown that it is sometimes due to charge donation from solvent lone pairs to the heavy metal which is not accounted for if the solvent is only implicitly treated as a polarizable continuum.59 In other cases there is no obvious evidence of solvent-to-solute electron donation, but more subtle interactions involving metal valence orbitals and solvent molecules can nonetheless cause dramatic solvent effects on metal–ligand coupling constants.29 Here we want to assess the performance of FDE to capture solvent effects of the latter type. We focus on the 1 J(199 Hg,13 C) couplings in Hg(CH3 )X (X = CH3 , Cl, Br, I) and study the influence of metal coordination by DMSO. The computation of J coupling constants involving heavy nuclei like Hg requires a relativistic theoretical framework. For isotropic J coupling constants involving Hg atoms, it was shown that scalar ZORA calculations with standard nonhybrid functionals are able to achieve good agreement with experiment for a variety of systems.3, 59, 66 Optimized geometries for the complexes were taken from reference 59 and are depicted in Figure 3. It is clear that the conformational space of the coordinated solvent molecules comprises many close-lying local minima which are in thermodynamic equilibrium in solution. In addition, the coordination number may vary on the time scale of NMR measurements. A complete description of solvent effects thus would have to capture these dynamic effects by proper sampling methods as has been done, for example, for the determination of solvent shifts of NMR chemical shifts of acetonitrile from FDE calculations.62 Indeed, Hg–C spin-spin coupling has been modeled by configurational averaging of ab initio MD configurations with explicit solvent molecules,29 which delivers excellent agreement with experiment. However, such an exhaustive modeling of solvent effects is not the purpose of the present work, as our main focus is assessing the performance of FDE as compared to full supermolecular KS

(a)

(b)

(c)

(d)

FIG. 3. Structures of solvated Hg complexes used in the J coupling calculations: (a) Hg(CH3 )2 + 3 DMSO, (b) Hg(CH3 )Cl + 3 DMSO, (c) Hg(CH3 )Br + 3 DMSO, (d) Hg(CH3 )I + 3 DMSO. For the FDE calculations the complexes have been partitioned into four fragments: the Hg(CH3 )X complex and each of the three DMSO solvent molecules.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-11

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

TABLE IV. Shift of the 1 J(199 Hg,13 C) coupling constant in the compounds Hg(CH3 )X (X = CH3 , Cl, Br, I) upon coordination of three DMSO molecules. The coupling constant for the isolated compounds (at the geometries of the DMSO complexes) is 402.3 Hz, 779.7 Hz, 747.4 Hz, and 691.2 Hz, respectively. The expansion basis (monomolecular or supermolecular) and the number of FT cycles is given in parentheses (J couplings in Hz). 1 J(199 Hg,13 C)

Hg(CH3 )2

Hg(CH3 )Cl

Hg(CH3 )Br

Hg(CH3 )I

KS

233.3

636.6

660.1

698.5

Thomas–Fermi FDE(m,0) FDE(m,1) FDE(m,3) FDE(s,3)

157.4 148.3 148.7 131.1

504.1 524.4 525.1 491.5

518.3 541.0 541.9 507.8

522.2 561.8 563.7 538.8

PW91k FDE(m,0) FDE(m,1) FDE(m,3) FDE(s,3)

114.1 112.1 112.6 90.6

436.6 467.9 468.8 440.6

453.0 486.9 487.9 460.1

465.8 515.2 517.3 504.5

calculations. Moreover, the optimized solvent-solute clusters from Ref. 59 were shown to produce a very good approximation of the solvent shift of the Hg–C J-coupling. For the FDE calculations the solvated Hg complexes of Figure 3 have been partitioned into four fragments: the Hg(CH3 )X complex and each of the three solvent molecules. Table IV summarizes the results for the shifts of the 1 199 J( Hg,13 C) coupling constants of the Hg(CH3 )X compounds upon DMSO complexation,47 as obtained from Eq. (49) for supermolecular KS calculations and from Eqs. (54) and (52) for FDE(m) and FDE(s), respectively. We first note that the J coupling constants of the isolated compounds as computed from KS-DFT are quite large, ranging from 400 Hz to 780 Hz, and that the solvent shifts are almost as large as the coupling constants in the isolated complexes, which has been noted frequently in prior work on metal-ligand J-coupling.29, 59, 64, 65 This indicates a strong perturbation of the electron density of the mercury complexes due to the coordinating DMSO molecules. A simple embedding with the Thomas–Fermi nonadditive kinetic energy functional and without polarization of the solvent subsystems, that is FDE(m,0), is able to recover between 67% and 79% of the solvent shift. Depending on the compound, this amounts to a deviation between 76 Hz (dimethylmercury) and 176 Hz (methylmercuryiodide) from the supermolecular KS calculations. This is a surprisingly good result since DMSO is a rather strongly coordinating solvent. Polarization of the electron density via freeze-thaw improves upon the FDE(m,0) values for the methylmercury halides, while a small deterioration is found for dimethylmercury. One FT cycle is sufficient to converge the results and FDE(m,3) yields between 64% (dimethylmercury) and 82% (methylmercurychloride) of the reference values. The agreement with the reference calculations deteriorates if a supermolecular expansion for the basis set is employed. Depending on the compound, FDE(s,3) recovers between 56% (dimethylmercury) and 77% (methylmercuryhalides) of the solvent shift.

FIG. 4. Computational efficiency of FDE as compared to KS-DFT for the ground state DFT calculation (SCF) and the calculation of the 1 J(199 Hg,13 C) couplings (CPL) in Hg(CH3 )Cl(DMSO)3 . SCF times for FDE include the time required for the SCF of all frozen subsystems (fragments). FDE(s,3) required a total of 3634 s (3515 s SCF and 119 s CPL) and is off the chart. Wall clock times were obtained on an 8-core Intel Xeon E5430 CPU (2.66 GHz).

The solvent shifts obtained with the PW91k non-additive kinetic energy functional are consistently underestimated as compared to the values obtained with the Thomas–Fermi functional. Otherwise, the same trends as discussed for the Thomas–Fermi functional are observed. Thus, the results which presumably should be the best, that is FDE(s,3) with the PW91k functional, recover between 39% (dimethylmercury) and 72% (methylmercuryiodide) of the solvent shift. For the investigated mercury complexes we conclude that FDE(m,1) with the Thomas–Fermi functional can be used to obtain a reasonable estimate of the magnitude and trends in solvent shifts. Improvements hinge on the development of better functionals for the non-additive kinetic energy. D. Computational performance

The aim of subsystem DFT methods is to reduce the computational effort that is required to compute molecular properties. In order to quantify the computational cost savings we have collected wall clock times for KS-DFT and FDE based calculations of J coupling constants. Results for the 1 J(199 Hg,13 C) couplings in Hg(CH3 )Cl(DMSO)3 are collected in Figure 4, while Figure 5 collects timings for the water cluster. On the computer which we have used for the timings, the self-consistent field (SCF) procedure to obtain the ground

FIG. 5. Computational efficiency of FDE as compared to KS-DFT for the ground state DFT calculation (SCF) and the calculation of the 1 J(O,H) couplings (CPL) in the cluster of 17 water molecules. SCF times for FDE include the time required for the SCF of all frozen subsystems (fragments). Wall clock times were obtained on an 8-core Intel Xeon E5430 CPU (2.66 GHz).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-12

Götz, Autschbach, and Visscher

state electron density of Hg(CH3 )Cl(DMSO)3 with supermolecular KS-DFT required 317 s. The subsequent calculation of the J couplings with the CPL program required 121 s. FDE calculations can significantly outperform the supermolecular KS-DFT calculations. FDE(m,0) required only 86 s for the SCF and 7 s for the property calculation, achieving a total speed-up of 4.7 over KS-DFT. Note that all FDE timings reported in Figure 4 include the time required for the SCF calculations of the frozen subsystems, that is the three DMSO solvent molecules. Using freeze-thaw to mutually polarize the active and frozen subsystems increases the SCF time significantly and FDE(m,3) calculations are almost as expensive as KS-DFT calculations. However, FDE(m,1) yields numerical results that are virtually identical to FDE(m,3) and is a factor of 2 faster than supermolecular KS-DFT for the combined SCF and property calculation. FDE calculations with a supermolecular expansion basis merely serve as a reference and for obvious reasons are considerably more expensive since they use the same basis set and numerical quadrature grid as supermolecular KS-DFT. The most obvious approach to improve the computational efficiency of FDE(s), if required, would be to only include basis functions from the environment subsystems that are in the vicinity of the active subsystem. These timings illustrate the potential that FDE has for a computationally efficient calculation of NMR J coupling constants. For our small Hg(CH3 )Cl(DMSO)3 test system we obtained a speed-up of 2 to 5 with FDE(m). This performance advantage of FDE increases for larger systems, as illustrated by the results for the water cluster in Figure 5. Supermolecular KS-DFT is significantly more expensive due to its nonlinear scaling computational cost with the system size, in particular for the NMR property calculation. The cost for the FDE property calculation, however, is independent of the total system size. Using FDE, we obtain speed-ups between 3 and 7 for the water cluster with setups that yield NMR J couplings that are very close to the reference supermolecular KS results or virtually exact. It is clear that the computational advantage of FDE will be more prominent for larger water clusters as has been shown, for example, for the calculation of NMR shielding tensors.62 V. CONCLUSIONS

We have derived the theoretical framework for a subsystem-based calculation of nuclear spin-spin coupling tensors in the framework of current-spin-DFT. In this subsystem formulation, FDE-CSDFT, we partition the electron density, the paramagnetic current, and the spin-magnetization density into individual contributions from each subsystem. The influence of the other subsystems (the environment) is thereby included via an effective embedding potential which depends on both the environment density and (induced) current and spin density. By neglecting the current dependence of the exchange–correlation and kinetic energy functionals, the DSO and PSO contributions to the nuclear spin-spin coupling tensor can be obtained directly from the subsystems’ ground state electron densities and KS orbitals, respectively. The calculation of the induced spin current, and consequently the FC and SD contributions to the nuclear spin-spin coupling

J. Chem. Phys. 140, 104107 (2014)

tensor, has to proceed in an iterative fashion across all subsystems, however, because the embedding potential depends on the spin-magnetization density of the environment. For the practical calculation of nuclear spin-spin coupling tensors in the presented FDE-CSDFT framework it is reasonable to exploit the fact that the perturbation due to the magnetic moment of a nucleus is very localized. With a proper choice of subsystems, the induced paramagnetic current and spin-magnetization density outside the active subsystem (the one which contains the coupled nuclei) can be made negligible so that it can be safely neglected. This leads to a very efficient computational scheme in which the calculation of nuclear spin-spin coupling tensors is restricted to the subsystem of interest while the influence of the environment is included via an effective embedding potential. This FDECSDFT based calculation of nuclear spin-spin coupling tensors is implemented in the CPL program which is part of the ADF software package. Our test calculations for hydrogen bound dimers demonstrate that the FDE scheme is able to reproduce the shift of 1 J(O,H) and 1 J(N,H) couplings upon hydrogen bond formation in close agreement with supermolecular KS-DFT calculations. For the water dimer and the ammonia-water dimer, already a simple embedding without freeze-thaw cycles and with a monomolecular expansion basis leads to good results. One freeze-thaw cycle is sufficient to converge the shift in the J couplings. Due to a more pronounced charge transfer between the hydrogen bond donor and hydrogen bond acceptor subsystems in the hydrogen fluoride dimer, a supermolecular basis set expansion provides a better agreement with the reference KS-DFT results for 1 J(F,H) couplings. Use of such basis sets is, however, impractical and may also lead to deterioration of results in case of significant covalent or coordination bonding between the subsystems.54, 67 For a larger cluster consisting of 17 water molecules we have shown that near-quantitative accuracy can be obtained with FDE for the 1 J(O,H) couplings in the central water molecule if the directly hydrogen-bonded water molecules are included in the active subsystem. A single freeze-thaw cycle is also in this case sufficient to converge the results. Restricting the active subsystems to the central water molecule that contains the coupled nuclei does not give satisfactory results leading to an error of around 25%. Surprisingly good results were obtained with FDE for the large shift of 1 J(199 Hg,13 C) couplings in Hg(CH3 )X complexes (X = CH3 , Cl, Br, I) upon coordination of three DMSO solvent molecules, where the complexes and each of the three DMSO molecules were treated as individual subsystems. The interaction of the coordinating solvent with mercury is rather strong and accompanied by exceptionally large solvent shifts between 233 Hz and 699 Hz. FDE recovers around 75% of this shift, somewhat depending on the compound and the details of the FDE calculation. The good performance of FDE supports the conclusions from Ref. 29 that changes in Hg–C bonding leading to the strong increase of J-couplings are simply due to the presence of nearby solvent molecules. Solvent-to-metal charge donation, which cannot be described by the frozen–density ansatz, only plays a secondary role in this case.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-13

Götz, Autschbach, and Visscher

We have also shown that the FDE(m) calculation of J couplings does significantly reduce the computational cost as compared to a regular, supermolecular calculation. For subsystems that are close to the coupled nuclei or that have otherwise a strong effect on the electron density and induced currents in proximity of the coupled nuclei, one freezethaw cycle is thereby required to converge the results. For Hg(CH3 )Cl(DMSO)3 , we obtained a speed-up of 2 to 5 with and without freeze-thaw, respectively. For the water cluster we achieve speed-ups of 3 to 7 for our most accurate approaches that include freeze-thaw. Our discussion of the timings includes not only the time required for the property calculation but also the SCF time required for the determination of the ground state density and orbitals. The speed-up for the property calculation itself is much larger. For bigger systems additional approximations can be invoked for subsystems that are sufficiently far from the subsystem of interest to enable further gains in computational efficiency. One could use, for example, an approximate electron density from an electronic structure method that is computationally cheaper than DFT. The success of FDE and the subsystem based calculation of nuclear spin-spin coupling tensors within the FDE-CSDFT framework hinges on appropriate partitioning of the system and the availability of sufficiently accurate functionals for the non-additive kinetic energy. Good results can be obtained with the currently available approximations such as Thomas– Fermi and PW91k for weakly interacting subsystems. Exact embedding schemes are emerging68–70 that can also be used for strongly interacting subsystems. This will make it possible to use the subsystem based DFT calculation of NMR J couplings presented in this work for a variety of interesting applications, ranging from a realistic modeling of solvent effects including a thorough sampling of solvent structures to the calculation of J couplings in proteins and other complex macromolecules.

ACKNOWLEDGMENTS

A.W.G. and L.V. would like to thank The Netherlands Organization for Scientific Research (NWO) for financial support in form of a VICI grant for L.V. and gratefully acknowledge computer time provided by the Dutch National Computing Facilities (NCF). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation, Grant No. OCI-1053575. Computer time was provided by the San Diego Supercomputer Center through NSF XSEDE Award No. TG-CHE130010 to A.W.G. J.A. acknowledges support from the National Science Foundation (NSF), Grant No. CHE-1265833. 1 T.

Helgaker, M. Jaszunski, and K. Ruud, Chem. Rev. 99, 293 (1999). Autschbach and T. Ziegler, Coord. Chem. Rev. 238–239, 83 (2003). 3 J. Autschbach, “The calculation of NMR parameters in transition metal complexes,” in Principles and Applications of DFT in Inorganic Chemistry, Structure and Bonding, Vol. 112, edited by N. Kaltsoyannis and J. E. McGrady (Springer, Heidelberg, 2004), pp. 1–48. 4 Calculation of NMR and EPR Parameters. Theory and Applications, edited by M. Kaupp, M. Bühl, and V. G. Malkin (Wiley-VCH, Weinheim, 2004). 2 J.

J. Chem. Phys. 140, 104107 (2014) 5 R.

G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford University Press, Oxford, 1989). 6 V. G. Malkin, O. L. Malkina, and D. R. Salahub, Chem. Phys. Lett. 221, 91 (1994). 7 R. M. Dickson and T. Ziegler, J. Phys. Chem. 100, 5286 (1996). 8 V. Sychrovský, J. Gräfenstein, and D. Cremer, J. Chem. Phys. 113, 3530 (2000). 9 T. Helgaker, M. Watson, and N. C. Handy, J. Chem. Phys. 113, 9402 (2000). 10 J. Autschbach and T. Ziegler, J. Chem. Phys. 113, 936 (2000). 11 J. Autschbach and T. Ziegler, J. Chem. Phys. 113, 9410 (2000). 12 M. A. Watson, P. Salek, P. Macak, M. Jaszunski, and T. Helgaker, Chem. Eur. J. 10, 4627 (2004). 13 J. Autschbach and S. Zheng, Annu. Rep. NMR Spectrosc. 67, 1 (2009). 14 A. S. P. Gomes and C. R. Jacob, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 108, 222 (2012). 15 P.-O. Åstrand, K. V. Mikkelsen, P. Jørgensen, and K. Ruud, J. Chem. Phys. 108, 2528 (1998). 16 K. Ruud, L. Frediani, R. Cammi, and B. Mennucci, Int. J. Mol. Sci. 4, 119 (2003). 17 J. Tomasi, B. Menucci, and R. Cammi, Chem. Rev. 105, 2999 (2005). 18 K. V. Mikkelsen, K. Ruud, and T. Helgaker, J. Comput. Chem. 20, 1281 (1999). 19 D. Zaccari, V. Barone, J. E. Peralta, R. H. Contreras, O. E. Taurian, E. Diez, and A. Esteban, Int. J. Mol. Sci. 4, 93 (2003). 20 V. Sychrovsky, B. Schneider, P. Hobza, L. Zidek, and V. Sklenae, Phys. Chem. Chem. Phys. 5, 734 (2003). 21 M. Pecul and K. Ruud, Magn. Reson. Chem. 42, S128 (2004). 22 B. Wang, X. He, and K. M. Merz, J. Chem. Theory Comput. 9, 4653 (2013). 23 T. A. Wesolowski and A. Warshel, J. Phys. Chem. 97, 8050 (1993). 24 T. A. Wesolowksi, “One-electron equations for embedded electron density: Challenge for theory and practical payoffs in multi-level modelling of complex polyatomic systems,” in Computational Chemistry: Reviews of Current Trends, edited by J. Leszczynski (World Scientific, 2006), Vol. 10, Chap. 1, pp. 1–82. 25 C. R. Jacob, J. Neugebauer, and L. Visscher, J. Comput. Chem. 29, 1011 (2008). 26 P. Cortona, Phys. Rev. B 44, 8454 (1991). 27 G. Senatore and K. R. Subbaswamy, Phys. Rev. B 34, 5754 (1986). 28 C. R. Jacob and J. Neugebauer, “Subsystem density-functional theory,” Wiley Interdiscip. Rev.: Comput. Mol. Sci. (published online). 29 S. Zheng and J. Autschbach, Chem. Eur. J. 17, 161 (2011). 30 C. R. Jacob and L. Visscher, J. Chem. Phys. 125, 194104 (2006). 31 G. Vignale and M. Rasolt, Phys. Rev. B 37, 10685 (1988). 32 C. van Wüllen, J. Comput. Chem. 23, 779 (2002). 33 J. Autschbach, J. Chem. Phys. 129, 094105 (2008); Erratum 130, 209901 (2009). 34 T. Helgaker, M. Jaszu´ nski, and M. Pecul, Prog. Nucl. Magn. Reson. Spectrosc. 53, 249 (2008). 35 J. Autschbach, ChemPhysChem 10, 2274 (2009). 36 E. van Lenthe, E. J. Baerends, and J. G. Snijders, J. Chem. Phys. 99, 4597 (1993). 37 E. van Lenthe, E. J. Baerends, and J. G. Snijders, J. Chem. Phys. 101, 9783 (1994). 38 E. van Lenthe, R. van Leeuwen, E. J. Baerends, and J. G. Snijders, Int. J. Quantum Chem. 57, 281 (1996). 39 T. A. Wesolowski and J. Weber, Chem. Phys. Lett. 248, 71 (1996). 40 See http://www.scm.com for ADF2013.01, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands (accessed 09/21/2013). 41 C. F. Guerra, J. Snijders, G. te Velde, and E. J. Baerends, Theoret. Chem. Acc. 99, 391 (1998). 42 G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. F. Guerra, S. J. A. van Gisbergen, J. G. Snijders, and T. Ziegler, J. Comput. Chem. 22, 931 (2001). 43 C. R. Jacob, S. M. Beyhan, R. E. Bulo, A. S. P. Gomes, A. W. Götz, K. Kiewisch, J. Sikkema, and L. Visscher, J. Comput. Chem. 32, 2328 (2011). 44 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 45 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997). 46 S. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). 47 See supplementary material at http://dx.doi.org/10.1063/1.4864053 for solvent shifts of J coupling constants in hydrogen-bonded dimers and mercury complexes using the PBE functional for the first-order perturbed XC potential. 48 L. H. Thomas, Proc. Cambridge Philos. Soc. 23, 542 (1927).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10

104107-14

Götz, Autschbach, and Visscher

J. Chem. Phys. 140, 104107 (2014)

49 E.

61 Y.

50 A.

62 R.

Fermi, Rend. Accad. Lincei 6, 602 (1927). Lembarki and H. Chermette, Phys. Rev. A 50, 5328 (1994). 51 T. A. Wesolowski, Y. Ellinger, and J. Weber, J. Chem. Phys. 108, 6078 (1998). 52 T. A. Wesolowski and F. Tran, J. Chem. Phys. 118, 2072 (2003). 53 C. R. Jacob, T. A. Wesolowski, and L. Visscher, J. Chem. Phys. 123, 174104 (2005). 54 A. W. Götz, S. M. Beyhan, and L. Visscher, J. Chem. Theory Comput. 5, 3161 (2009). 55 T. A. Wesolowski, J. Chem. Phys. 106, 8516 (1997). 56 J. Neugebauer, M. J. Louwerse, P. Belanzoni, T. A. Wesolowski, and E. J. Baerends, J. Chem. Phys. 123, 114101 (2005). 57 J. Neugebauer, C. R. Jacob, T. A. Wesolowski, and E. J. Baerends, J. Phys. Chem. A 109, 7805 (2005). 58 M. Dulak, J. W. Kaminski, and T. A. Wesolowski, J. Chem. Theory Comput. 3, 735 (2007). 59 J. Autschbach and T. Ziegler, J. Am. Chem. Soc. 123, 3341 (2001). 60 M. A. Watson, N. C. Handy, A. J. Cohen, and T. Helgaker, J. Chem. Phys. 120, 7252 (2004).

Zhao and D. G. Truhlar, J. Chem. Theory Comput. 1, 415 (2005). E. Bulo, C. R. Jacob, and L. Visscher, J. Phys. Chem. A 112, 2640 (2008). 63 H. Cybulski, M. Pecul, and J. Sadlej, Chem. Phys. 326, 431 (2006). 64 K. Sutter, L. A. Truflandier, and J. Autschbach, ChemPhysChem 12, 1448 (2011). 65 J. Autschbach and M. Sterzel, J. Am. Chem. Soc. 129, 11093 (2007). 66 J. Autschbach, and T. Ziegler, “Relativistic calculation of spin-spin coupling constants,” in Calculation of NMR and EPR Parameters. Theory and Applications, edited by M. Kaupp, M. Bühl, and V. G. Malkin (Wiley-VCH, Weinheim, 2004). 67 S. Fux, K. Kiewisch, C. R. Jacob, J. Neugebauer, and M. Reiher, Chem. Phys. Lett. 461, 353 (2008). 68 J. D. Goodpaster, N. Ananth, F. R. Manby, and T. F. Miller III, J. Chem. Phys. 133, 084103 (2010). 69 J. D. Goodpaster, T. A. Barnes, and T. F. Miller III, J. Chem. Phys. 134, 164108 (2011). 70 F. R. Manby, M. Stella, J. D. Goodpaster, and T. F. Miller III, J. Chem. Theory Comput. 8, 2564 (2012).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 77.99.133.148 On: Fri, 02 May 2014 07:37:10