Quadrupolar NMR Relaxation from ab Initio Molecular Dynamics

Subscriber access provided by University of Florida | Smathers Libraries

Article

Quadrupolar NMR relaxation from ab-initio molecular dynamics: Improved sampling and cluster models vs. periodic calculations Adam Philips, Alex Marchenko, Lionel Alexandre Truflandier, and Jochen Autschbach J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00584 • Publication Date (Web): 18 Jul 2017 Downloaded from http://pubs.acs.org on July 20, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Quadrupolar NMR relaxation from ab-initio molecular dynamics: Improved sampling and cluster models vs. periodic calculations Adam Philips,† Alex Marchenko,† Lionel A. Truflandier,‡ and Jochen Autschbach ,∗ †

Institut des Sciences Moléculaires Université Bordeaux, CNRS UMR 5255 351 cours de la Libération 33405 Talence cedex ‡

Department of Chemistry University at Buffalo State University of New York Buffalo, NY 14260-3000, USA email: [email protected] †

Abstract: Quadrupolar NMR relaxation rates are computed for 17O, 2H nuclei of liquid water, and of 23Na+, and 35Cl– in aqueous solution, via Kohn-Sham (KS) density functional theory ab-initio molecular dynamics (aiMD) and subsequent KS electric field gradient (EFG) calculations along the trajectories. The calculated relaxation rates are within about a factor of two of experimental results, and improved over previous aiMD simulations. The relaxation rates are assessed with regard to the lengths of the simulations as well as configurational sampling. The latter is found to be the more limiting factor in obtaining a good statistical sampling and is improved by averaging over many equivalent nuclei of a system or over several independent trajectories. Further, full periodic plane-wave basis calculations of the EFGs are compared with molecular-cluster atomic-orbital basis calculations. The two methods deliver comparable results with non-hybrid functionals. With the molecular-cluster approach, a larger variety of electronic structure methods is available. For chloride, the EFG computations benefit from using a hybrid KS functional.

1 Introduction Nuclear magnetic resonance (NMR) spectroscopy is a versatile tool for studying chemical bonding and structural properties at the molecular level. In addition to the main observable parameters in the liquid state, i.e. liquid state chemical shifts and 𝐽 -coupling, relaxation processes associated with magnetic resonance phenomena can provide detailed information about the dynamics of a molecular system. Based on this fundamental connection between relaxation and dynamics, many experimental and theoretical studies have used NMR relaxation as a probe for the dynamics of molecular systems.1–15 Ions such as sodium and chloride are of particular interest, because of their ubiquity in biological systems.16, 17 Nuclei with spin quantum number 𝐼 ≥ 1 have a non-spherical distribution of nuclear charge, being either prolate or oblate, and have an associated quadrupole moment. Elements with quadrupolar nuclei make up over 70% of the periodic table18 and include isotopes such as 2H, 17O, 23Na, and 1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

35

Cl. The associated NMR relaxation mechanism is the quadrupolar interaction with the electric field gradient (EFG) at the nucleus.1, 19 The quadrupole moment interacts with an EFG brought about by a non-spherical surrounding electric charge distribution from the electrons and other nuclei. This interaction is known to dominate NMR relaxation when present.1 The theoretical expressions for the relaxation rates (or relaxation times) depend on molecular parameters such as the EFG at the nucleus of interest, the quadrupole moment, the spin quantum number, and a set of data that depend on the system’s dynamic behavior.20, 21 The dynamics is usually quantified either in terms of characteristic correlation times, or time-dependent probabilities for the EFG, or in terms of EFG autocorrelation functions (ACFs) and their integrals. Further details are provided in Section 2. The ACF-based formalism is particularly useful for theoretical and computational work, because the ACFs can be determined computationally from the time-evolution of the EFG tensor, i.e. from molecular dynamics (MD) trajectories. The EFG is a molecular probe that is very sensitive to details of structure and bonding.22, 23 This feature also implies that the EFG is sensitive to approximations in calculations and – perhaps except for weakly polarizable ions24 – it requires a first-principles quantum theoretical treatment with electron correlation, either by Kohn-Sham (KS) density functional theory or with wavefunction-based methods. These computations may use combined quantum mechanics / molecular mechanics (QM/MM) in order to render larger systems tractable. Different methods may be applied for the dynamics vs. EFG calculations. The dynamics itself can be carried out in a variety of ways, including dynamics with force-fields (MM), dynamics using QM with a classical treatment of the nuclei (often referred to as ‘ab-initio MD’, or aiMD), full quantum dynamics (which is computationally very demanding), and combinations thereof such as aiMD with corrections for nuclear quantum effects, or QM/MM MD. Most aiMD calculations are to date based on KS theory and utilize periodic boundary conditions (PBCs). Typically accessible aiMD simulation times presently range well below the nanosecond regime. Nonetheless, the fact that these calculations do not require empirical parameters renders them very attractive for simulations of physical phenomena that are not necessarily aimed at bulk properties (for which well-performing force fields are available, especially for aqueous systems). Rather, for phenomena which depend on a physically meaningful description of short-timescale dynamic processes affecting structure and bonding, aiMD is an attractive and useful approach. Quadrupolar NMR relaxation belongs to this class of phenomena. There is functionality available to calculate EFGs essentially on the fly in KS aiMD simulations with PBCs, but this functionality is typically limited to non-hybrid KS functionals. Such an approach was taken in pioneering work by Schmidt et al.25 to simulate quadrupolar relaxation rates for liquid heavy water. Our group developed an approach for quadrupolar NMR relaxation that relies on aiMD simulations and subsequent cluster model KS calculations of the EFGs along 2 ACS Paragon Plus Environment

Page 2 of 31

Page 3 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the trajectory with a molecular quantum chemistry code.12 It was shown that the relaxation rates of alkali and halide ions in aqueous solution, experimentally spanning more than 5 orders of magnitude, could be predicted from first principles. However, deviations between simulations and experiments by a factor of 3 or more were common, similar to the findings of Schmidt et al.25 Our previous study was very limited in terms of sampling: Only one aiMD simulation of between 13 to 20 ps per ion was performed, and the relaxation rates were determined from 4 to 7 ps time windows at the end of the trajectories. The limited sampling is manifest, for example, in that the calculated relaxation times 𝑇1 and 𝑇2 were not the same because of rotational anisotropy in the simulations. A disadvantage of the cluster-based approach is that the cluster size needs to be converged carefully in order to obtain meaningful results for the EFGs. A great advantage, on the other hand, is that the molecular quantum chemistry ‘toolbox’ is very versatile, potentially allowing relaxation rate simulations based on EFG calculations with hybrid KS theory or correlated wavefunctions, or with all-electron relativistic methods in the case of heavy ions such as Cs+ or I–. The limited sampling in our previous study has motivated the present work. We use multiple independent trajectories from which the EFG ACFs and relaxation parameters can be averaged to simulate the on-average spherical symmetry around a solvated ion in liquid water better, and to obtain more realistic error estimates. Further, the relaxation rates of 17O and 2H in liquid heavy water are re-examined. Car-Parinello aiMD was utilized to generate relaxation data for 23Naaq+, 35 Claq– , and for 17O and 2H nuclei in pure heavy water. The two ions were chosen because in our previous work excellent agreement was obtained with the experimental relaxation rate of 23Na+, but possibly because of error cancellation, while the calculated rates for 35Cl– overestimated the experiment by between 3 to 5 times, depending on details in the relaxation formalism used. Pure water was chosen for comparisons with Schmidt et al.25 Among the questions to be addressed in the present study are: (i) Can configurational sampling be improved and spherical isotropy be established by averaging over multiple trajectories, as far as the NMR relaxation rates are concerned? (ii) Do cluster model calculations yield comparable EFGs to those of PBC PAW calculations, based on the same aiMD trajectories? (iii) Are the dynamic events that drive nuclear quadrupolar relaxation sampled sufficiently within the time scale of the aiMD simulations? (iv) Can the poor 35Cl– data of the previous study be improved via KS calculations of the EFG tensors with hybrid functionals?

2 Theory In order to render the article self-contained, the theoretical formalism is briefly summarized here. This formalism is identical to the one detailed in our previous work,12 with only minor changes 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 31

to numerical methods and implementation. A trapezoidal (as opposed to rectangular) Riemann sum type method is used for numerical integration of ACFs, and the ACFs themselves are computed using the Wiener-Khinchin theorem; a discussion of which can be found in the Supporting Information (SI). The theory describing nuclear spin relaxation was pioneered by Bloembergen, Purcell, and Pound20 as well as Redfield.21 We use the formalism as presented by Spiess1 and Cowan19 to obtain relaxation rates from time-dependent EFG data. The longitudinal (1∕𝑇1 ) and transverse (1∕𝑇2 ) relaxation rates are given by 1 = 𝐶 𝑄 𝐺2,0 (𝜔) 𝑇1 1 = 𝐶 𝑄 𝐺2,1 (𝜔) 𝑇2

(1a) (1b)

where 𝜔 is the Larmor frequency and 𝐶 𝑄 is a constant depending on the nuclear quadrupole moment 𝑄, and the nuclear spin quantum number 𝐼: 𝐶𝑄 =

𝑒2 𝑄2 (2𝐼 + 3) 40𝐼 2 (2𝐼 − 1) ℏ2

(1c)

Further, 𝑒 is the elementary charge, and ℏ is the reduced Planck constant. The superscript 𝑄 denotes relaxation arising from the nuclear quadrupolar interaction. The 𝐺𝑙,𝑚 (𝜔) are formally functions of frequency and defined as linear combinations of spectral densities, 𝑔𝑙,𝑚 (𝜔), via 𝐺2,0 (𝜔) = 4𝑔2,2 (2𝜔) + 4𝑔2,−2 (−2𝜔) + 𝑔2,1 (𝜔) + 𝑔2,−1 (−𝜔)

(2a)

𝐺2,1 (𝜔) = 2𝑔2,−2 (−2𝜔) + 3𝑔2,−1 (−𝜔) + 2𝑔2,1 (𝜔) + 3𝑔2,0 (0)

(2b)

The spectral densities are defined as half-Fourier transforms of the auto-correlation functions (ACFs), 𝑓𝑙,𝑚 (𝜏), of the time-domain EFGs: 𝑔𝑙,𝑚 (𝜔) =

∫0

∞

𝑓𝑙,𝑚 (𝜏)𝑒𝑖𝜔𝜏 𝑑𝜏

(3)

𝑅∗𝑙,𝑚 (𝑡)𝑅𝑙,𝑚 (𝑡 + 𝜏)𝑑𝑡

(4)

In turn, an ACF 𝑓𝑙,𝑚 (𝜏) for the EFG is defined as 𝑓𝑙,𝑚 (𝜏) =

∫𝑡0

∞

where 𝑅𝑙,𝑚 (𝑡) is an element of the time-domain second-rank irreducible spherical EFG tensor. 4 ACS Paragon Plus Environment

Page 5 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Further, 𝑡0 is a time origin onward from which the time-dependent function, or a discrete time series in the case of discrete data, is integrated, and 𝜏 is the delay time. The ACF is often expressed as 𝑓𝑙,𝑚 (𝜏) = ⟨𝑅∗𝑙,𝑚 (𝑡)𝑅𝑙,𝑚 (𝑡 + 𝜏)⟩ (5)

where the brackets denote an ensemble average. In Equations (4) and (5), the irreducible components 𝑅2,𝑚 of the EFG tensor are defined in terms of the Cartesian EFG tensor elements 𝑉𝛼𝛽 , with 𝛼, 𝛽 ∈ {𝑥, 𝑦, 𝑧}, via √ 1 𝑅2,0 = 3 𝑉 (6a) 6 𝑧𝑧 𝑅2,±1 = ∓𝑉𝑥𝑧 − 𝑖𝑉𝑦𝑧

(6b)

1 𝑅2,±2 = (𝑉𝑥𝑥 − 𝑉𝑦𝑦 ) ± 𝑖𝑉𝑥𝑦 2

(6c)

For a random, stationary process, the ACFs 𝑓𝑙,𝑚 are real.1, 25, 26 Therefore, it is necessary and appropriate to consider only the real portion of the ACFs computed from the complex quantities 𝑅𝑙,𝑚 . Finally, the quantities 𝑉𝛼𝛽 are the Cartesian second derivatives of the electric potential 𝜙, 𝑉𝛼𝛽 =

(

𝜕2𝜙 𝜕𝛼𝜕𝛽

)

(7)

nuc

evaluated in the laboratory coordinate frame at the position of the nucleus. In MD-based computations, the relaxation rates are therefore determined from the Cartesian EFG tensor elements calculated at snapshots in time over the course of one or more MD trajectories. The 𝑉𝛼𝛽 are known to be sensitive to approximations in the electronic structure model, and for heavy elements relativistic effects can be important. The definition of the ACF for applications with discrete data is 𝑁−𝜏−1 ∑ 1 𝑓𝑙,𝑚 (𝜏) = 𝑅∗𝑙,𝑚 (𝑡)𝑅𝑙,𝑚 (𝑡 + 𝜏) 𝑁 − 𝜏 𝑡=𝑡

(8)

0

where 𝑁 is the total number of discrete data points collected. Since real-world acquirable data are discrete and finite, this is the expression used in practice. There are varying definitions and numerical methods for the calculation of discrete ACFs. For the present work, the ACFs were computed by using the Wiener-Khinchin theorem,27 as discussed in more detail in the Supporting Information (SI). The general utility of an ACF is to quantify how a signal or a set of data is correlated to itself as a function of a time delay. From it, a correlation time 𝜏𝑐 can be defined which quantifies the 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 31

time scale of the phenomena that cause the autocorrelation. The correlation time can be thought of as the time required for the system to lose its memory of the initial state, such that at 𝑡 ≫ 𝜏𝑐 the ACF vanishes. The correlation time can be quantified by the width of the decaying ACF. In the context of this work, correlation times 𝜏𝑙,𝑚 for the different EFG tensor components are calculated as ∞ 𝑓 (𝜏) ∞ ⟨𝑅∗ (𝑡)𝑅 (𝑡 + 𝜏)⟩ 𝑙,𝑚 𝑙,𝑚 𝑙,𝑚 𝑑𝜏 = 𝜏𝑙,𝑚 = 𝑑𝜏 (9) 2 ∫ ∫𝑡0 𝑓 ⟨|𝑅𝑙,𝑚 (𝑡0 )| ⟩ 𝑡0 𝑙,𝑚 (𝑡0 )

There are further approximations that simplify the relaxation theory. The ‘fast motion’ or ‘extreme narrowing’ limit is valid for non-viscous systems in which molecular motion is fast and minimally restricted.1, 19, 26 This condition is satisfied when 𝜔𝜏𝑐 ≪ 1, i.e. for systems where the correlation times are small. At this limit, the exponential of equation (3) is nearly equal to one for the time domain over which the ACF is decaying (i.e. a time domain not much larger than the correlation time). The calculation of the spectral densities then simplifies to a direct integration of the EFG ACFs: ∞ 𝑓 (𝜏)𝑑𝜏 (10) 𝑔𝑙,𝑚 = ∫𝑡0 𝑙,𝑚 Each of the 𝜏𝑙,𝑚 can be concisely expressed as 𝜏𝑙,𝑚 =

𝑔𝑙,𝑚 𝜎𝑙,𝑚

(11)

where the 𝑔𝑙,𝑚 are the spectral densities defined above. The 𝜎𝑙,𝑚 are defined as the 𝑓𝑙,𝑚 (𝑡0 ) in the denominator of equation (9) and represent the variance of the respective EFG component over the sampled configurations. Note that 𝑔𝑙,𝑚 is no longer a function of frequency so, at this extreme narrowing limit, the relaxation rates are proportional to the value of the spectral densities at zero frequency. These values are then used in the linear combinations of Equations (2). In an isotropic liquid, the analyte and its solvent shell adopt all possible orientations. This condition would presumably be met when data from many, very long, MD trajectories are averaged, but in practice the rotational isotropy can be difficult to establish. Under the condition of rotational isotropy, the ACFs of the spherical EFG tensor components 𝑓𝑙,𝑚 , and corresponding spectral densities 𝑔𝑙,𝑚 become invariant with respect to 𝑚 such that 𝑔2 = 𝑔2,±2 = 𝑔2,±1 = 𝑔2,0

(12)

𝑄 It follows then that the 𝐺𝑙,𝑚 of equations (2) are also invariant with respect to 𝑚 and can be reduced to a single 𝐺2𝑄 𝐺2𝑄 = 10𝑔2 (13)

6 ACS Paragon Plus Environment

Page 7 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The 𝑔2 in this equation can, in principle, be a choice of any of the 𝑔𝑙,𝑚 or – in the current case of finite simulations of relatively short timescale – an average of them. The ultimate consequence of this rotational isotropy condition is that the longitudinal and transverse relaxation rates become equal and an isotropic rate can be defined as 1 1 1 = = = 10𝐶 𝑄 𝑔2 𝑇iso 𝑇1 𝑇2

(14)

This isotropic relaxation rate can be expressed in terms of the correlation time 𝜏𝑐 , and the total EFG variance. The total EFG variance is defined here as the sum of the individual component variances 2 ∑ ⟨𝑉 (0)2 ⟩ = 𝜎2,𝑚 (15) 𝑚=−2

In the absence of rotational anisotropy, a single correlation time 𝜏𝑐 can be defined as 2 ∑ 1 𝑔 𝜏𝑐 = ⟨𝑉 (0)2 ⟩ 𝑚=−2 2,𝑚

(16)

Equation (14) can then be expressed such that the isotropic relaxation rate is explicitly proportional to the correlation time and total EFG variance. 1 = 2𝐶 𝑄 𝜏𝑐 ⟨𝑉 (0)2 ⟩ 𝑇iso

(17)

The factor of two arises from the factor of ten in equation (13) multiplied by 1∕5 from averaging of the five tensor components. The quadrupolar relaxation rate is proportional to the EFG autocorrelation time. This result reveals that short correlation times, which correspond to more rapid averaging of the EFG fluctuations at the nucleus due to fast molecular motion, that is, more effective averaging of fields around the nucleus, render the relaxation mechanism less effective. It should be kept in mind that with finite sampling of the configuration space, the 𝑔2,𝑚 are not necessarily numerically equivalent to each other. In this case, 1∕𝑇iso can be defined in different formally equivalent ways, leading to slightly different numerical answers. The choice presented above is used for the numerical data presented in Section 4. An alternative would be to take 1∕𝑇iso simply as the average of 1∕𝑇1 and 1∕𝑇2 .

3 Computational Details The aiMD simulations were carried out using the Quantum-Espresso (QE) software package.28 These simulations utilized ultrasoft pseudopotentials and the exchange-correlation functional of 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Perdew, Burke, and Emzerhof (PBE).29, 30 The simulation cells included 64 water molecules. The sodium and chloride simulations also included the respective analyte ion, and a counter ion of either H+ in the case of chloride or OH– in the case of sodium. The periodic cells were cubic, with dimensions chosen such that the density was ∼ 1 g/cm3 in the case of the ions in light water and slightly higher (1.11 g/cm3 ) for the heavy water box. Specifically, the box dimensions were 12.42 Å for water, 12.55 Å for sodium, and 12.54 Å for chloride. kinetic energy cutoff of 100 Ry was chosen for the plane-wave (PW) basis and a fictitious electron mass (𝜇) of 380 a.u. was used for the Car-Parinello dynamics, with a time step of 6.0 a.u. (0.145 fs). The choice for 𝜇 is consistent with the recommendation of Grossman et al.31 and Schwegler et al.32 In previous works by our group, a value of 𝜇 = 450 a.u. was employed, however, a reduced electronic fictitious mass for water simulations resulted in less structured radial distribution functions in better agreement with experiment,33, 34 as well as temperature stability over the course of the simulation because the electronic degrees of freedom are energetically more separated from those of the nuclei. Random initial cell packing and 5 ps of pre-equilibration with classical MD in NVT with the molecular modeling package Tinker35 preceded each aiMD simulation. The aiMD simulations consisted of ∼ 5 ps of equilibration in the canonical (NVT) ensemble at 350 K using a three-chain Nosé-Hoover thermostat (90, 45, and 15 THz), followed by a production phase in the microcanonical (NVE) ensemble. The first ∼ 1 ps in the NVE ensemble was also taken to be part of equilibration. In accordance with precedents set in the literature,36–40 an elevated simulation temperature of 350 K was chosen, along with semi-empirical van der Waals corrections to the energies and forces,41, 42 to correct for the over-structured, or ‘glassy’43 nature of water typically observed in aiMD simulations with non-hybrid functionals such as PBE. Electric field gradient (EFG) computations for the ions and pure water were carried out using the gauge-including projector-augmented wave (GIPAW) module of QE28, 44 with periodic boundary conditions (PBCs). The ‘GI’ functionality is needed for NMR shielding calculations, but not for EFGs, and therefore we refer to these calculations by using the ‘PAW’ acronym. EFGs for the ion-water clusters were also computed using the Amsterdam Density Functional (ADF) software package,45–47 utilizing non-periodic boundary conditions, all-electron Slater-type orbital (STO) triple-𝜁 basis sets with two sets of polarization functions (TZ2P) for the analyte ion, and a double-𝜁 polarized basis set (DZP) for oxygen and hydrogen. For consistency with our previous work,12 scalar relativistic effects were treated by the Zeroth Order Regular Approximation (ZORA),48 although for the isotopes considered herein these effects may be neglected. Calculations were carried out using both non-hybrid (PBE) and hybrid (PBE0, 25% exact exchange)49 functionals. All ion-water clusters consisted of the ion of interest and the 20 nearest neighbor solvent molecules from the instantaneous configuration taken from the MD. The convergence of the calculated EFG was studied in our previous work, and 20 explicit solvent molecules were found 8 ACS Paragon Plus Environment

Page 8 of 31

Page 9 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to be sufficient. Additional bulk solvent effects were simulated via the Conductor like Screening Model (COSMO).50 EFG tensors were calculated for evenly spaced snapshot configurations along the production phase of a simulation. Long (40 ps) trajectories utilized between 1024 and 4595 configurations (specified where appropriate) and shorter (4 to 10 ps) trajectories used between 256 and 1024 configurations. PAW EFGs were computed on the full systems at each selected time step, maintaining PBCs. These computations were accomplished in two steps: first a single-point SCF was performed in QE for the selected MD snapshot. The density was computed using a 100 Ry plane-wave kinetic energy cutoff and the same cell parameters as in the CPMD simulation. Then EFGs were computed given the previously determined electron density. The PAW norm-conserving pseudopotentials were those distributed by Ceresoli.51 These are large core pseudopotentials, which restrict the 2s and 2p electrons of sodium and chloride to the core. Given the energetic separation of these atomic orbitals, this is a good approximation, especially for chloride. As a check, PAW EFG calculations were carried out for a representative sodium simulation with a small core pseudopotential (2s and 2p electrons included in valence) and were shown to not significantly impact computed relaxation rates (see SI). The relaxation rates were calculated with an in-house re-implementation of the QE ‘DynPro’ relaxation rate module of Reference 12. The new program utilizes the IPython framework52 and is more easily extensible to other relaxation mechanisms than the original code. Trajectory analyses, cluster extractions, and visualizations were performed with an open source tool recently developed in our group (version 0.3.6).53 Nuclear quadrupole moments used for the relaxation rate calculations were taken from Reference 54. In order to distinguish between the EFG calculations with all-electron STO basis sets on molecular clusters extracted from the aiMD trajectories, versus the plane-wave basis PAW calculations with PBCs on the full simulation cells, we designate the former by ‘STO’ and the latter by ‘PAW’. Note that it is also possible to carry out STO (or other atomic orbital type basis) calculations with PBCs. However, in our previous work,12 a cluster model with a solvent shell of 20 water molecules around an ion, in conjunction with COSMO for bulk solvation, was identified as a suitably converged model for the ion relaxation calculations. We therefore consider the STO vs. PAW comparisons mainly to probe the differences in the electronic structure (in particular for the much more polarizable halide ion), and in the case of PAW also the quality of how the all-electron electronic structure is reconstructed by the projector method.



3.5

Water box (D2O)

3.0

Chloride Sodium

H2O Exp(2008) D2O Exp(2008)

2.5

2.5

2.0

2.0

gOO(r)

gOO(r)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 31

1.5

1.5 1.0

1.0 0.5

0.5 0.0

0.0 2

3

4

5

6

2

3

˚ r (A)

4

5

6

˚ r (A)

Figure 1: Left: Oxygen-oxygen radial distribution functions (RDF) for each of the three explored systems. 𝑔OO (𝑟) denotes the RDF or pair correlation as a function of pair distance in Angstrom. Right: Experimental oxygen-oxygen RDFs for heavy and light water from joint x-ray diffraction and neutron scattering experiments of Soper34 (extracted using a plot digitizer55 ).

4 Results and Discussion 4.1

Assessment of Molecular Dynamics

This subsection assesses some of the structural properties obtained from the molecular dynamics trajectories used for the current study. Sets of ‘long’ and ‘short’ trajectories were simulated to assess the problems of finite sampling and rotational anisotropy. Each set of independent trajectories consisted of one ‘long’ simulation with 40 ps of production, and ten ‘short’ simulations each with 4 ps of production (up to 10 ps for pure water). All data presented here and in the following sections were obtained from the production phases of the simulations. The radial distribution functions (RDF) reported in Figures 1–3 were generated from the long trajectory of the respective system (sodium, chloride, and pure water). The oxygen-oxygen (O-O) RDF gives a description of the structure of the simulated water and is used as a benchmark for the quality and accuracy of the aiMD. Non-hybrid functionals, such as PBE used for the present work, are known to over-structure water,32, 43 which is manifested as a higher O–O first peak, and a deeper first trough, when compared to RDFs obtained from neutron scattering and/or x-ray diffraction experiments.33, 34 This over-structuring is sometimes referred to as a ‘glassy’ character of water simulated with DFT.43 Higher simulation temperatures (350 K) and van der Waals dispersion corrections were employed in order to minimize these effects. Although the temperature is not controlled by a thermostat during the production phase (NVE) from which the pair correlation data is taken, the instantaneous temperatures remain close to the target and do not drift significantly (except in a preliminary pure 10 ACS Paragon Plus Environment

5

35

RDF Count

30

4 25 3

20

2

15

nNaO(r)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


gNaO(r)

Page 11 of 31

10 1

5 0

0 2

3

4

5

6

˚ r (A)

Figure 2: Na-O RDF from 40 ps sodium trajectory. The number of water molecules in the first solvation shell is determined from the integration of the first peak and suggests the sodium is between 5- and 6-coordinate with average coordination distance of 2.4 Å water trajectory generated with a larger value of 𝜇, in which cooling was observed in the NVE production phase). Overall, the O-O RDFs, in comparison with the literature, indicate a sufficiently accurate simulation of the intermolecular structure of liquid water at the current level of theory. Differences in the RDFs for the ion simulations originate from the disruptive effects of the ion on the hydrogen bonding network as well as from the choice to use light water in these simulations. The RDFs of figures 2–3 provide information regarding the solvation structure around sodium and chloride ions in the simulations. Here, the count function, 𝑛(𝑟), is also shown which integrates the RDF and indicates the number of solvent molecules in the solvation shell associated with a particular RDF peak. The structure and dynamics of the first aqueous solvation shell is thought to have the most influence on the EFG at the quadrupolar nucleus and therefore to be the driving mechanism of relaxation.9, 14 The mean radius of the first solvation shell for Na+ corresponds to the first peak in the RDF and is found to be at 2.42 Å. Integration of this peak up to the first local minimum gives 5.4 Na-O pairs, suggesting that the sodium transitions between 5- and 6-coordinate throughout the simulation. This coordination motif is consistent with our previous observations.12 In the case of an anion like chloride, the coordination number can be obtained either from the Cl-O or Cl-H RDF. The first solvation peak, with a maximum at a distance of 2.16 Å, then integrates to 5.3 Cl-H pairs. These paired hydrogens correspond to the water hydrogens which tend to orient directly towards the chloride anion. The second hydrogen in each of the solvating waters is directed away and does not participate in direct solvation, as shown in Figure 4 for a representative aiMD snapshot in comparison with the solvated sodium ion. The initial peak of the Cl-O RDF also integrates to approximately 5. Again, these data suggest a preference for a 5coordinate configuration with some tendency to adopt a 6-coordinate structure. This is consistent 11 ACS Paragon Plus Environment


3.0

2.5

RDF Count

2.5

RDF Count

30

60

25

2.0

50

20

1.5

40

1.0

30 1.0

10

0.5

20 0.5

5

0.0 4

5

10

0.0

0 3

nClH(r)

15

gClH(r)

1.5

nClO(r)

2.0

gClO(r)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 31

6

0 2

˚ r (A)

3

4

5

6

˚ r (A)

Figure 3: Cl-O (left) and Cl-H (right) RDFs from 40 ps chloride trajectory. Both indicate a preferred coordination number of 5 water molecules. First solvation peak of the Cl-H RDF corresponds to 5.3 Cl-H pairs at average distance of 2.16 Å.

Figure 4: Configuration snapshot from the 23Na+ trajectory (left) and 35Cl– trajectory (right). Six nearest neighbor water molecules are opaque while all others are rendered translucent. with previous simulations.12, 56 4.2

Relaxation Results

The sets of relaxation data from the computations are presented in Tables 1–4. Plots of ACFs are shown in Figures 5, 8, and 9. A summary of relaxation data from our formally best calculations, i.e. with EFGs obtained with the hybrid functional for 23Na+ and 35Cl–, are compared with available experimental relaxation rates in Table 5. For pure water, the cluster model is less straightforward 12 ACS Paragon Plus Environment

Page 13 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and we used the PAW calculations of the EFG tensors with the PBE non-hybrid functional, in order to take advantage of the improved sampling that comes from averaging over all molecules in the simulation box. While far from a complete statistical sampling, the set of ten shorter trajectories allows for an assessment of the sampling errors associated with the calculations. There are several sources of sampling error in the calculation of the relaxation rates and as a consequence the results reported in Tables 1–4 may not be indicative of the true accuracy and precision of these results. One source of statistical error is the number of configurations chosen along the trajectory to compute the EFGs. It is important that the grid of time steps to generate the relaxation data is chosen such that the time between selected frames is very much shorter than the correlation time. This condition is taken to be sufficiently achieved here, even with the most sparse grid spacing used for a calculation at around 0.04 ps. Furthermore, using sparser grid sampling in calculations for the sodium and chloride systems does not significantly affect the relaxation results, at least down to 1024 configurations across a 40 ps trajectory and 256 configurations across a 4 ps trajectory (see SI Tables S5 – S7). This suggests that the selected grid spacing is not a limiting factor in the precision of the results. More fundamental sources of statistical error originate from the finite nature of the aiMD simulations. A consequence is that a single simulation on a time scale on the order of 10 to 100 ps may not sample a representative configurational space of the real system. This configurational sampling is important in validating the spherical isotropy condition inherent in the relaxation calculations (see section 2). Averaging over a set of independent trajectories allows for an assessment of this source of error but the sampling remains incomplete, of course. It is important to note that on top of the various issues related to the sampling, the simulations and EFG calculations also afford various types of systematic errors related to the approximations made in the calculations, chiefly among them the chosen electronic structure methods. For the simulations of pure (heavy) water, Figure 5 shows the ACFs of the EFG components for a single long trajectory (left) and averages over ten short trajectories (right). Both sets of ACFs decay quickly within a few picoseconds. As noted in previous studies on quadrupolar relaxation,9, 10, 24 a key feature of the MD derived EFGs is a multi-step decay process, as opposed to a simple exponential decay which has sometimes been assumed.20, 57, 58 This type of decay is seen in the ACFs in Figure 5, as well as in the ACFs for the other nuclei studied in this work, consisting of an initial rapid falloff followed by a slower decay to the baseline. The initial steep decay, suggested by Schmidt et al.25 to be due to libration modes, makes up much of the loss of correlation for chloride and sodium (less so for 17O and 2H in water), but the secondary decay accounts for most of the spectral density and therefore relaxation.14 In the region of the EFG ACF which transitions between these two decay regimes, there is a small local maximum or ‘notch’ which, to the best of our knowledge, is a feature not addressed in previous studies on water25 or 13 ACS Paragon Plus Environment


1.0

1.0

f2,0 f2,1 f2,2

0.8

f2,0 f2,1 f2,2

0.8

0.6

0.6

ACF

ACF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 31

0.4

0.4

0.2

0.2

0.0

0.0 0

5

10

15

20

25

0

2

Time (ps)

4

6

8

10

Time (ps)

Figure 5: Autocorrelation functions (ACFs) of the second rank irreducible tensor components of the electric field gradient (EFG) averaged over the 64 oxygen nuclei of the water box. EFGs computed using PAW method. Left: 4595 configurations from single, 40 ps trajectory. The horizontal axis is cut short in the figure. Right: Average ACF over ten 10 ps trajectories, each computed from 1151 configurations. the monovalent ions.9, 14 This feature becomes more apparent when averaging over a set of trajectories and is particularly prominent for 35Cl–. It is possible that this feature is an artifact of the Table 1: 17O relaxation data𝑎 for the set of pure heavy water trajectories at 350 K ⟨ ⟩ 1 (s−1 ) 𝑇1 (s−1 ) 𝑇1 (s−1 ) 𝜏𝑐 (ps) 𝑉 (0)2 𝑇 1

2

iso

01 78.03 75.60 76.41 0.95 3.51 02 140.00 149.65 146.44 1.91 3.36 03 123.10 132.35 129.27 1.68 3.36 04 104.69 102.33 103.12 1.30 3.47 05 133.65 127.00 129.22 1.65 3.43 06 156.00 124.31 134.87 1.76 3.35 07 102.46 92.66 95.93 1.20 3.49 08 133.28 130.78 131.61 1.71 3.36 09 140.12 120.86 127.28 1.68 3.32 10 140.54 128.43 132.47 1.79 3.24 Average 125.19 118.40 120.66 1.56 3.39 St.Dev 22.20 20.72 20.48 0.29 0.08 St.Err 7.02 6.55 6.48 0.09 0.03 long(40ps) 89.82 89.48 89.59 1.11 3.54 𝑎 EFGs computed using PAW method. All quantities are averages over 64 water oxygens in each trajectory. The numbers 01–10 in the first column refer to the different short trajectories. 4595 configurations of the 40 ps trajectory are sampled. Each of the 10 short trajectories are 10 ps and 1151 configurations are sampled. ⟨𝑉 (0)2 ⟩ in atomic units. 14 ACS Paragon Plus Environment

Page 15 of 31

1.0

1.0

f2,0 f2,1 f2,2

0.8

f2,0 f2,1 f2,2

0.8

0.6

0.6

ACF

ACF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.4

0.4

0.2

0.2

0.0

0.0 0

5

10

15

20

25

0

2

Time (ps)

4

6

8

10

Time (ps)

Figure 6: ACFs of the second rank irreducible tensor components of the EFG averaged over the 128 deuterium nuclei of the water box. EFGs computed using PAW method. Left: 4595 configurations from single, 40 ps trajectory. The horizontal axis is cut short in the figure. Right: Average ACF over ten 10 ps trajectories, each computed from 1151 configurations. finite sampling. However, the same feature can be seen in the force-field MD based EFG-ACFs of Carof et al. for aqueous chloride.14 The relative smoothness of the ACFs for 17O and 2H (Figures 5 and 6) compared to those of Table 2: 2H relaxation data𝑎 for set of trajectories of heavy water at 350 K ⟨ ⟩ 1 1 1 −1 −1 −1 2 (s ) (s ) (s ) 𝜏 (ps) 𝑉 (0) 𝑐 𝑇 𝑇 𝑇 1

2

iso

01 1.44 1.35 1.38 1.08 0.30 02 2.53 2.48 2.50 2.02 0.29 03 2.18 2.29 2.25 1.82 0.29 04 1.77 1.70 1.72 1.36 0.30 05 2.19 2.21 2.20 1.75 0.29 06 2.38 2.24 2.28 1.85 0.29 07 1.69 1.67 1.68 1.32 0.30 08 2.35 2.21 2.25 1.82 0.29 09 2.14 2.11 2.12 1.73 0.29 10 2.20 2.25 2.24 1.85 0.28 Average 2.09 2.05 2.06 1.66 0.29 St.Dev 0.33 0.34 0.33 0.28 0.01 St.Err 0.10 0.11 0.10 0.09 0.00 long(40ps) 1.45 1.35 1.39 1.08 0.30 𝑎 EFGs computed using PAW method. All quantities are averages over 128 deuterium nuclei in each trajectory. The numbers 01–10 in the first column refer to the different short trajectories. 4595 configurations of the 40 ps trajectory are sampled. Each of the 10 short trajectories are 10 ps and 1151 configurations are sampled. ⟨𝑉 (0)2 ⟩ in atomic units. 15

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the ions (Figures 8 and 9) arises from the averaging of the unique 64 oxygens and 128 deuterium nuclei, respectively, per trajectory. For the water simulations, there is close agreement between the ACFs of EFG components with varying the 𝑚 index (cf. equation 4). This lack of a dependence on 𝑚 suggests that the rotational isotropy condition is satisfied to a good degree. Based on a visual inspection of the ACFs for the long vs. short trajectories, it is not immediately evident that averaging over a set of ten trajectories further enforces rotational isotropy. However, the data from an analogous set of trajectories using light water shows more pronounced deviations in the different ACF components for the long trajectory, whereas virtually no differences are seen for the 10-trajectory average (see SI). This means that features of a particular trajectory may lead to rotational anisotropy, even if the EFGs are averaged over many nuclei, due to the finite sampling of the configurational space. The discrete nature of the simulations also introduces statistical error in the computations. The discrete ACF for a finite time series is an approximation of the true ACF. At longer time delays, the discrete, finite ACF acquires statistical errors, due to a smaller sample of data points, and is necessarily forced to zero at the maximum time delay. Such finite sampling effects are not expected to introduce a significant amount of error to the relaxation rates as long as the correlation time is much shorter than the length of the simulation. For the sodium and chloride trajectories, the computed correlation time is well under 1 ps, and the ACFs decay to zero very quickly (see Figures 8 and 9). The water trajectories have correlation times on the order of 1 ps. It can be seen in Figures 5 and 6 that the time it takes for the ACFs to decay to zero is much longer than for the ions, prompting us to longer simulation times for the short water trajectories. The 17O relaxation data for water are collected in Table 1. For both the single long and the ten short trajectory averages, 1∕𝑇1 , 1∕𝑇2 , and 1∕𝑇iso agree within the standard error, as one would expect from the ACF plots in Figure 5. The individual 17O relaxation rates, per trajectory, show large variations, as quantified by the sizable standard deviations of ∼20 Hz relative to the rates of ∼100 Hz. Likewise, in the corresponding 2H data in Table 2, the standard deviations are about 1/6 of the rates. Here, due to the even larger number of nuclei (128) available for the averaging, 1∕𝑇1 , 1∕𝑇2 , and 1∕𝑇iso are in excellent agreement in the 10-trajectory average. Both for 17O and 2 H, the variations of the relaxation rates among the ten short trajectories are almost entirely due to the variations in the correlation times, as the ⟨𝑉 (0)2 ⟩ remain fairly consistent. This is also the case when the short and long trajectories are compared, meaning that in all of the simulations the dynamic fluctuations of the EFG tensor components are within a comparable range. For the ion relaxation rate calculations, we employed both the molecular cluster STO basis set (short: ‘STO’) and the infinite-periodic PAW approaches (short: ‘PAW’). Figure 7 assesses the agreement between EFGs calculated for 23Na+ and 35Cl– with the different methods. Each point on a given scatter plot corresponds to the 𝑉𝑧𝑧 component of an ion’s Cartesian EFG tensor for one of the selected configurations along the trajectory. When the STO and PAW methods are 16 ACS Paragon Plus Environment

Page 16 of 31

Page 17 of 31

a

0.2

0.1

Vzz ADF PBE0

Vzz ADF PBE

0.2

0.0

0.0

−0.1

−0.2

−0.2

0.6

−0.1

0.0

0.1

Vzz QE-PAW

b

0.1

−0.1

−0.2

0.2

−0.2

c

0.6

−0.1

0.0

Vzz ADF PBE

0.1

0.2

d

0.4

Vzz ADF PBE0

0.4

Vzz ADF PBE

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.2 0.0 −0.2

0.2 0.0 −0.2

−0.4

−0.4

−0.6

−0.6 −0.6

−0.4

−0.2

0.0

0.2

Vzz QE-PAW

0.4

0.6

−0.6

−0.4

−0.2

0.0

0.2

Vzz ADF PBE

0.4

0.6

Figure 7: 𝑉𝑧𝑧 EFG component (a.u.) at 23Na+ and 35Cl– nuclei from the respective 40 ps trajectory. Scattered points each represent a different configuration along the trajectory. Pearson correlation coefficient reported as 𝑟 and 𝑦 = 𝑥 is plotted for comparison. (a) 23Na+ PAW versus STO, with PBE functional, 𝑟 = 0.9344. (b) 23Na+ PBE functional versus PBE0 hybrid functional, STO basis, 𝑟 = 0.9998. (c) 35Cl– PAW versus STO, with PBE functional, 𝑟 = 0.9796. (d) 35Cl– PBE functional versus PBE0 hybrid functional, STO basis, 𝑟 = 0.9996. compared using the same non-hybrid density functional, the scatter in the plots is testament to the technically very different ways by which the EFG is calculated in the two programs. For instance, scalar relativistic effects are accounted for explicitly in the STO calculations, but implicitly via the pseudopotential approximation for the PAW approach. However, the nuclear charges are relatively small for the ions and approximations in treating relativistic effects are not a major source of error. The plane-wave basis used in the periodic calculations is near-complete for the valence region. The all-electron STO basis is far from complete, but it should give a more balanced description of the valence vs. outer core polarization. The PAW calculations treat the long-range contributions to the EFG correctly, but the required frozen-core approximations may lead to deviations from all-electron results. This is apparent when comparing the larger relative deviations observed for 23Na+ with respect to 35Cl–. For both ions, in the plane-wave calculations the 𝐾 and



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

𝐿 core shells are represented by a pseudopotential. The 2s2p orbitals of the chlorine atom are much lower in energy than the 3s3p orbitals (on the order of 160 eV). For sodium, the energy gap is smaller (about 40 eV). Therefore, it was hypothesized that use of a small core pseudopotential for Na with 2s2p in the valence space may reduce the discrepancies between the two types of calculations. Calculations with the small core pseudopotential for sodium were carried out and shown to produce comparable scatter with respect to the all-electron STO data as the large-core potential and very similar relaxation rates (See SI). However, the EFG tensor asymmetry may be more strongly affected by the type of pseudopotential. We finally note that for the 23Na+ trajectory, the EFG fluctuations are smaller than for 35Cl–, as this ion is rather weakly polarizable, and therefore long-range contributions play a relatively larger role. Nonetheless, the range of EFGs sampled by the two methods is comparable, and even for sodium the correlation coefficient of 0.93 supports the idea that the EFG is dominated by shortrange interactions of the ion with surrounding water molecules. The scatter plot does not suggest a systematic bias toward higher or lower EFGs with either of the methods. The chloride ion is much easier polarizable than sodium, which leads to a larger range of EFG values sampled along the trajectory. There is also better agreement between the STO and PAW calculations (𝑟 = 0.98). There is near perfect agreement between the EFG calculations with the PBE and PBE0 functionals for 23Na+ (𝑟 = 1.00), meaning that both functionals treat the EFG equally well (or equally poorly). For 35Cl–, the correlation between PBE0 vs. PBE is still quite good (𝑟 = 0.98), but the plot reveals that the data lie almost perfectly on a straight line with a slope of less than 1. Since the hybrid functional should be more accurate, this finding suggests that the non-hybrid functional over-estimates the EFGs at the most positive and most negative extremes. By comparison of the PBE0 with the PBE data in Table 3 (the full set of PBE data is provided in Table S2), this overestimation is clearly reflected in larger ⟨𝑉 (0)2 ⟩ values for the PBE functional. The autocorrelation times, on the other hand, are very similar. Consequently, the smaller range of chloride EFG values sampled with the PBE0 functional results in smaller relaxation rates, by about 18% for 1∕𝑇iso in the 10-trajectory average STO calculations, which improves the agreement with experiment (Table 5). This sensitivity to the functional is expected to be even more important for larger and easily polarizable ions such as bromide and iodide. The ACFs for 35Cl– and 23Na+ (Figures 8 and 9) are significantly noisier than those from the water simulations, due to the limited sampling of a single solvated ion per simulation. Compared to the single, long, trajectory, the 10-trajectory averages provide some improvements, but it is likely that hundreds of trajectories would need to be averaged in order to satisfy the rotational isotropy limit to a similar degree as in the water simulations. As observed for the pure water system, a two-part ACF decay is apparent for both ions. Sodium exhibits a significantly more prominent initial fast decay than chloride. That is, this fast decay regime represents most of the 18 ACS Paragon Plus Environment

Page 18 of 31

Page 19 of 31

1.0

1.0

f2,0 f2,1 f2,2

0.8

0.6

ACF

ACF

f2,0 f2,1 f2,2

0.8

0.6 0.4

0.4

0.2

0.2 0.0

0.0 0

5

10

15

0

20

1

2

3

4

Time (ps)

Time (ps)

Figure 8: ACFs of the second rank irreducible tensor components of the EFG at the 35Cl– nucleus. EFGs computed using STO method with PBE0 hybrid functional. Left: 4096 configurations from single, 40 ps trajectory (horizontal axis cut short in figure). Right: Average ACF over ten 4 ps trajectories, each computed from 1024 configurations 1.0

1.0

f2,0 f2,1 f2,2

0.8

f2,0 f2,1 f2,2

0.8

0.6

0.6

ACF

ACF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.4 0.2

0.4 0.2

0.0

0.0 0

5

10

15

20

0

1

2

3

4

Time (ps)

Time (ps)

Figure 9: ACFs of the second rank irreducible tensor components of the EFG at the 23Na+ nucleus. EFGs computed using STO method with PBE0 hybrid functional. Left: 4096 configurations from single, 40 ps trajectory. Right: Average over ten 4 ps trajectories, each computed from 1024 configurations. EFG de-correlation for 23Na+, and the relaxation rate, whereas for 35Cl– the slower de-correlation regime plays a more important role. As mentioned already, the chloride ACF exhibits a small ‘notch’ between the two decay regimes. Further investigation of this and other features is left for future work, where we will compare chloride with its heavier halide counterparts. Relaxation data for 35Cl– and 23Na+ are provided in Tables 3 and 4 respectively. Due to the more limited sampling, as compared to the water simulations, the computed longitudinal and transverse relaxation rates for both systems differ significantly. This discrepancy is taken to represent an incomplete realization of the spherical isotropy condition, and is alleviated when sampling is 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 31

improved by averaging across the ten independent trajectories. In the case of sodium, differences in the relaxation rates across different trajectories are entirely determined by variations in the correlation times, with the total EFG variance virtually independent of the particular trajectory. However, for chloride, the total EFG variance is found to vary measurably across trajectories. This finding can be attributed to the larger polarizability of the chloride ion, which can be sensitive to the solvation environment sampled along the short trajectories. Averaging across independent trajectories allows for an estimation of the statistical error based on sampling. The standard error in the average of 1∕𝑇iso over ten short trajectories amounts to 28% for chloride, 22% for sodium, and 5% for both 17O and 2H. Again, this greater precision for oxygen and deuterium is likely attributed to sampling over many nuclei in these aiMD simulations. Table 3: 35Cl– relaxation data𝑎 for set of chloride trajectories at 350 K. ⟨ ⟩ 𝑉 (0)2 STO-PBE0 𝑇1 (s−1 ) 𝑇1 (s−1 ) 𝑇1 (s−1 ) 𝜏𝑐 (ps) 1

01 02 03 04 05 06 07 08 09 10 Average St.Dev St.Err long(40ps)

23.21 11.10 67.37 44.81 102.45 29.46 19.64 22.49 19.16 143.10 48.28 41.23 13.04 127.67

2

iso

15.21 9.11 114.22 21.38 112.28 76.12 40.84 29.59 10.27 190.82 61.98 57.33 18.13 53.93

17.88 9.77 98.61 29.19 109.01 60.57 33.77 27.22 13.23 174.92 57.42 51.21 16.19 78.51

0.11 0.06 0.41 0.17 0.56 0.31 0.17 0.14 0.08 0.56 0.26 0.18 0.06 0.43

0.16 0.18 0.25 0.17 0.20 0.20 0.21 0.20 0.17 0.32 0.21 0.05 0.01 0.19

STO-PBE Average 59.15 75.08 69.77 0.27 0.24 St.Dev 50.94 68.75 61.85 0.19 0.05 St.Err 16.11 21.74 19.56 0.06 0.02 long(40ps) 161.37 63.74 96.28 0.45 0.22 𝑎 EFGS are computed using the STO method with both the PBE0 and PBE functionals. For PBE0, 1024 configurations of the 40 ps trajectory are sampled while each of the 10 short trajectories are 4 ps and 256 configurations are sampled. For PBE, 4096 configurations of the 40 ps trajectory are sampled while each of the 10 short trajectories are 4 ps and 1024 configurations are sampled. ⟨𝑉 (0)2 ⟩ in atomic units. Here, only PBE0 results for individual short trajectories are reported and are denoted by the numbers 01–10. Complete results from PBE are in the SI. 20 ACS Paragon Plus Environment

Page 21 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Except for 23Na+, the long trajectories yield a relaxation rate that falls within the range of those from the ten short ones. For sodium, the largest 1∕𝑇iso value from the short 23Na+ trajectories is one standard deviation below the 1∕𝑇iso obtained from the long one. I.e., either the long trajectory represents an atypical situation, or the short simulations have a bias toward small relaxation rates. By dividing the long sodium aiMD simulation into 4 ps time windows, it can be seen that EFGs taken from the first 4, and in particular, from the last 4 ps result in a relaxation rate much greater than what is seen in the rest of the trajectory or from any one of the independent short simulations (SI Table S8). This may suggest that relatively rare dynamic events on this time scale play a role in driving the quadrupolar relaxation of 23Na. For chloride, the long trajectory data fall within the range of those obtained from the short ones. Table 4: 23Na+ relaxation data𝑎 for set of sodium trajectories at 350 K. ⟨ ⟩ 𝑉 (0)2 STO-PBE0 𝑇1 (s−1 ) 𝑇1 (s−1 ) 𝑇1 (s−1 ) 𝜏𝑐 (ps) 01 02 03 04 05 06 07 08 09 10 Average St.Dev St.Err long(40ps)

1

2

iso

1.65 9.79 15.50 8.63 9.91 3.05 14.37 2.14 2.55 1.63 6.92 5.12 1.62 6.76

3.20 2.68 4.16 4.25 18.17 3.57 5.82 1.61 2.69 3.40 4.95 4.53 1.43 25.72

2.68 5.05 7.94 5.71 15.42 3.40 8.67 1.79 2.64 2.81 5.61 3.95 1.25 19.40

0.06 0.10 0.20 0.13 0.31 0.08 0.25 0.04 0.06 0.06 0.13 0.09 0.03 0.42

0.03 0.03 0.03 0.03 0.03 0.03 0.02 0.03 0.03 0.03 0.03 0.00 0.00 0.03

STO-PBE Average 7.00 4.98 5.65 0.13 0.03 St.Dev 5.21 4.39 3.91 0.09 0.00 St.Err 1.65 1.39 1.24 0.03 0.00 long(40ps) 6.78 25.88 19.51 0.42 0.03 𝑎 EFGs are computed using the STO method with both the PBE0 and PBE functionals. For PBE0, 4096 configurations of the 40 ps trajectory are sampled while each of the 10 short trajectories are 4 ps and 512 configurations are sampled. For PBE, 4096 configurations of the 40 ps trajectory are sampled while each of the 10 short trajectories are 4 ps and 1024 configurations are sampled. ⟨𝑉 (0)2 ⟩ in atomic units. Here, only PBE0 results for individual short trajectories are reported and are denoted by the numbers 01–10. Complete results from PBE are in the SI. 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In comparison with experiments (reported in summary Table 5), the short trajectory average and the long trajectory relaxation rate for 17O fall below the experiment (141 Hz), with the former being quite close (121 Hz). For 2H the long trajectory result is very close to the experiment (1.4 Hz) whereas the short trajectory average is too high by about a factor of 1.5. However, the water results are within a factor of 2 from experiment, which is very acceptable, given the difficulties of calculating these relaxation rates from first principles, and improved over the aiMD simulations of Schmidt et al.25 The computed rate for 35Cl– from the ten trajectories is approximately a factor of 2 too high compared to experiment (29 Hz), but much improved over our previous results.12 In the case of 23Na+, the results bracket the experimental relaxation rate of 16 Hz, with the short trajectory average nearly a factor of 3 lower and the long trajectory result just over. Overall, these results reflect an improvement toward experiment when compared to previous computations of quadrupolar relaxation rates from first principles.12, 25

5 Summary and Conclusions Quadrupolar relaxation rates of 17O, 2H, 23Na+, and 35Cl– in water were computed from aiMD simulations and subsequent first principles calculations of EFG tensors and their ACFs. Results were overall improved over previous studies by averaging over sets of independent trajectories. While not necessarily always closer to experiment, when compared to a single longer simulation, this averaging allows for some quantification of the statistical errors. In addition, sampling from multiple trajectories more closely represents the condition of rotational isotropy in the liquid system, as exemplified by the better agreement between computed longitudinal and transverse relaxation Table 5: Summary of relaxation data𝑎 for the different trajectories ⟨ ⟩ 1 −1 𝜏𝑐 (ps) 𝑉 (0)2 (s ) 𝑇1 , exp (s−1 ) 𝑇 iso

35

–

Cl (40ps) 0.43 – Cl (10 × 4ps) 0.26(6) 23 Na+(40ps) 0.42 23 Na+(10 × 4ps) 0.13(3) 17 O(40ps) 1.11 17 O(10 × 10ps) 1.56(9) 2 H(40ps) 1.08 2 H(10 × 10ps) 1.66(9) 35

𝑎

0.19 0.21(1) 0.03 0.028(1) 3.54 3.39(3) 0.30 0.29()

79 57(16) 19 6(1) 90 121(6) 1.4 2.1(1)

1

29.25 16.25 140.83 1.44 -

Results reported for 23Na+ and 35Cl– are from the STO method with PBE0. Results for 17O and 2 H are from the PAW method with PBE. ⟨𝑉 (0)2 ⟩ in atomic units. References indicated for experimental values. 22 ACS Paragon Plus Environment

Page 22 of 31

Page 23 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


rates for the ten short trajectory averages versus a single longer simulation. These relaxation rates must be equivalent in the fast motion isotropic limit. EFGs computed using PAW vs. STO methods agree with correlation coefficients greater that 0.9 when using the 20 nearest neighbor solvent clusters suggested in our group’s previous work.12 Due to its larger polarizability, the EFG at the 35Cl– nucleus exhibits greater variance than the other nuclei studied, and furthermore some functional dependence is observed for its EFGs. It would seem that hybrid functionals such as PBE0 are needed to sufficiently describe the nuclear EFGs of this and other even more polarizable ions such as bromide, iodide, and neutral heavy atoms such as xenon. Based on the accuracy of our computed relaxation rates and the observation that the EFG ACFs are well converged, our data suggest that the relevant dynamic events that drive the quadrupolar relaxation of these nuclei are well sampled on the picosecond timescale. While this appears to be the case even for the 4 ps ion trajectories, increased correlation times are observed in short time windows of the longer sodium trajectory, which may suggest that rare and subtle dynamic events around the analyte can have an important influence on the quadrupolar relaxation rate. Pure water exhibits longer correlation times and simulations of at least 10 ps are required to converge the ACFs. Outside of the pure water simulations, a major limitation is the configurational sampling that is accessible when the EFG at only one nuclear site is available from the simulations. Averaging over many nuclei in a single trajectory or over several (preferably many) independent trajectories appears to be an effective strategy for improving this configurational sampling. Based on our findings, the five questions posed in the Introduction can be answered as follows: (i) Sampling of multiple independent trajectories appears to be advantageous for the purpose of calculating quadrupolar relaxation rates, even if the trajectories are relatively short. Sampling of rare events requires a larger number of trajectories as employed here. (ii) Cluster models of the solvated ions perform well, in particular for more polarizable ions for which the greater variety of molecular electronic structure methods may be advantageous. (iii) Most of the dynamics driving the relaxation rates of the systems studied herein appears to take place over a few ps. We cannot rule out that longer time-scale events contribute much to the observed rates, since the agreement with experiment is not perfect. However, there are numerous approximations in the calculations, and incomplete sampling remains a concern. (iv) The 35Cl– relaxation is much improved over previous aiMD results by averaging multiple trajectories and using a hybrid functional for the EFG calculations. A hybrid functional would also be desirable for the dynamics itself,59 but at much increased computational cost. With the exception of the short trajectory average for 23Na+, the computed relaxation rates are within a factor of two of the experimental data, representing an important improvement over our previous studies12 and those by Schmidt et al.25 Relaxation rates for selected ions computed in other works from classical force-field dynamics and using the Sternheimer approximation for 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the EFGs,9, 24 have reported similar agreement with experiments. Even though classical MD can access much larger time scales, it relies on uniquely tailored force fields to reproduce the structural features of the analyte-solvent interactions and is unlikely to provide accurate EFGs for heavy ions that are very polarizable. In addition the validity of the Sternheimer approximation for computing EFGs of aqueous ions has been brought into question.60 Ab-initio methods offer the potential for transferable and accurate simulations and properties calculations for systems in which the details of the electronic structure play a significant role. Recent steps have been taken by groups such as DiStasio et al.36 in vastly improving the quality of methods such as aiMD simulations for reproducing the experimentally observed structure and dynamics of liquid water. As they become more computationally accessible, these kinds of ab-initio methods may be used in the context of predicting NMR relaxation rates driven by various mechanisms such as dipolar interactions and chemical shift anisotropy (CSA) as well as many other dynamically driven molecular processes.

Acknowledgments We acknowledge the Center for Computational Research (CCR) at the University at Buffalo for providing computational resources. This work has been supported by grant CHE-1560881 from the National Science Foundation.

Supporting Information Available Explanation of Wiener-Khinchin theorem, full relaxation data for the 23Na+ and 35Cl– ions computed using the PAW method as well as the STO method with PBE functional, 23Na+ PAW EFGs from small core versus large core pseudopotentials and resultant relaxation rates, dependence of relaxation data on chosen number of configurations/grid spacing, ACFs and relaxation data for light water trajectory (using 1H), 23Na+, 17O, and 2H relaxation data computed from time windows of varying length over 40 ps heavy water trajectory, plots of integrals and log of EFG ACFs for heavy water 10 ps trajectories, statistical analysis of isotropic relaxation rates computed for individual 17O nuclei over the trajectories. This material is available free of charge via the Internet at http://pubs.acs.org.

References (1) Spiess, H. W., ‘Rotation of Molecules and Nuclear Spin Relaxation’, in P. Diehl, R. K. E. Fluck (editor), ‘NMR Basic Principles and Progress’, Vol. 15, Springer, Berlin, 1978, 55–214. 24 ACS Paragon Plus Environment

Page 24 of 31

Page 25 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(2) Hansen, M. R.; Graf, R.; Spiess, H. W., ‘Interplay of Structure and Dynamics in Functional Macromolecular and Supramolecular Systems As Revealed by Magnetic Resonance Spectroscopy’, Chem. Rev. 2016, 116, 1272–1308. (3) Denison, A.; Rabideau, S.; Garrett, B. B., ‘Oxygen-17 Relaxation in Water’, J. Phys. Chem. 1967, 71, 2606–2611. (4) Powles, J.; Rhodes, M.; Strange, J., ‘Deuteron spin lattice relaxation in benzene, bromobenzene, water and ammonia’, Mol. Phys. 1966, 11, 515–523. (5) Weingärtner, H.; Hertz, H., ‘Magnetic Relaxation by Quadrupolar Interaction of Ionic Nuclei in Non-Aqueous Electrolyte Solutions. Part I. Limiting Values for Infinite Dilution’, Ber. Bunsen-Ges. Phys. Chem. 1977, 1204–1221. (6) Hanni, M.; Lantto, P.; Vaara, J., ‘Nuclear spin relaxation due to chemical shift anisotropy of gas-phase 129 Xe’, Phys. Chem. Chem. Phys. 2011, 13, 13704–13708. (7) Engström, S.; Jönsson, B.; Bengt, J., ‘A molecular approach to quadrupole relaxation. Monte Carlo simulations of diluite Li+, Na+, and Cl– aqueous solutions’, J. Magn. Reson. 1982, 50, 1–20. (8) Engström, S.; Jönsson, B., ‘Molecular dynamic simulation of quadrupole relaxation of atomic ions in aqueous solution’, J. Chem. Phys. 1984, 80, 5481–5486. (9) Roberts, J. E.; Schnitker, J., ‘Ionic quadrupolar relaxation in aqueous solution: dynamics of the hydration sphere’, J. Phys. Chem. 1993, 97, 5410–5417. (10) Chen, S. W.; Rossky, P. J., ‘Influence of Solvent and Counterion on 23Na+ Spin Relaxation in Aqueous Solution’, J. Phys. Chem. 1993, 97, 10803–10812. (11) Laage, D.; Hynes, J. T., ‘On the Molecular Mechanism of Water Reorientation’, J. Phys. Chem. B 2008, 112, 14230–14242. (12) Badu, S.; Truflandier, L. A.; Autschbach, J., ‘Quadrupolar NMR Spin Relaxation Calculated Using Ab-initio Molecular Dynamics: Group 1 and Group 17 Ions in Aqueous Solution’, J. Chem. Theory Comput. 2013, 9, 4074–4086. (13) Stirnemann, G.; Wernersson, E.; Jungwirth, P.; Laage, D., ‘Mechanisms of Acceleration and Retardation of Water Dynamics by Ions’, J. Am. Chem. Soc. 2013, 135, 11824–11831.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14) Carof, A.; Salanne, M.; Charpentier, T.; Rotenberg, B., ‘On the microscopic fluctuations driving the NMR relaxation of quadrupolar ions in water’, J. Chem. Phys. 2015, 143, 194504. (15) Lindgren, M.; Laaksonen, A.; Westlund, P.-O., ‘A theoretical spin relaxation and molecular dynamics simulation study of the Gd(H2O)93+ complex’, Phys. Chem. Chem. Phys. 2009, 11, 10368–10376. (16) Nordenskiold, L.; Chang, D. K.; Anderson, C. F.; Jr., M. T. R., ‘Sodium-23 NMR relaxation study of the effects of conformation and base composition on the interactions of counterions with double-helical DNA’, Biochemistry 1984, 23, 4309–4317. (17) Hald, M.; Jacobsen, J., ‘23Na Relaxation in DNA solutions. Influence of Intercalation on corrleation Times and Quadrupolar Coupling Constants’, J. Chem. Phys. 1992, 159, 257– 267. (18) Man, P., Quadrupolar interactions, Encyclopedia of Nuclear Magnetic Resonance, John Wiley & Sons, Chichester, U.K., 1996. (19) Cowan, B., Nuclear Magnetic Resonance and Relaxation, Cambridge University Press, Cambridge, UK, 2005. (20) Bloembergen, N.; Purcell, E.; Pound, R., ‘Relaxation effects in nuclear magnetic resonance absorption’, Phys. Rev. 1948, 73, 679–715. (21) Redfield, A. G., ‘On the Theory of Relaxation Processes’, IBM J. Res. Develop. 1957, 1, 19–31. (22) Schwerdtfeger, P.; Pernpointner, M.; Nazarewicz, W., ‘Calculation of Nuclear Quadrupole Coupling Constants’, in Kaupp, M.; Bühl, M.; Malkin, V. G. (editors), ‘Calculation of NMR and EPR Parameters. Theory and Applications’, Wiley-VCH, Weinheim, 2004, 279–291. (23) Autschbach, J.; Zheng, S.; Schurko, R. W., ‘Analysis of Electric Field Gradient Tensors at Quadrupolar Nuclei in Common Structural Motifs’, Concepts Magn. Reson. A 2010, 36A, 84–126. (24) Carof, A.; Salanne, M.; Charpentier, T.; Rotenberg, B., ‘Accurate Quadrupolar NMR Relaxation Rates of Aqueous Cations from Classical Molecular Dynamics’, J. Phys. Chem. B 2014, 118, 13252–13257.


Page 26 of 31

Page 27 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(25) Schmidt, J.; Hutter, J.; Spiess, H.-W.; Sebastiani, D., ‘Beyond Isotropic Tumbling Models: Nuclear Spin Relaxation in Liquids from First Principles’, ChemPhysChem 2008, 9, 2313– 2316. (26) Abragam, A., ‘The Principles of Nuclear Magnetism’, Am. J. Phys. 1961, 29, 860. (27) Wiener, N., Time Series, M.I.T. Press, Cambridge, Ma, 1964. (28) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Corso, A. D.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; MartinSamos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M., ‘QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials’, J. Phys. Cond. Mat. 2009, 21, 395502. (29) Perdew, J. P.; Burke, K.; Wang, Y., ‘Generalized gradient approximation for the exchangecorrelation hole of a many-electron system’, Phys. Rev. B 1996, 54, 16533–16539. (30) Perdew, J. P.; Burke, K.; Ernzerhof, M., ‘Reply to Comment on generalized gradient approximation made simple’, Phys. Rev. Lett. 1998, 80, 891. (31) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G., ‘Towards an assessment of the accuracy of density functional theory for first principles simulations of water.’, J. Chem. Phys. 2004, 120, 300–11. (32) Schwegler, E.; Grossman, J. C.; Gygi, F.; Galli, G., ‘Towards an assessment of the accuracy of density functional theory for first principles simulations of water. II’, J. Chem. Phys. 2004, 121, 5400–5409. (33) Soper, A., ‘Joint structure refinement of x-ray and neutron diffraction data on disordered materials: application to liquid water’, J. Phys.: Condens. Matter 2007, 19, 335206. (34) Soper, A. K.; Benmore, C. J., ‘Quantum Differences between Heavy and Light Water’, Phys. Rev. Lett. 2008, 101, 065502. (35) TINKER - Software Tools for Molecular Design by J. W. Ponder. URL https:// dasher.wustl.edu/tinker/ Accessed 06/17. (36) Distasio, R. A.; Santra, B.; Li, Z.; Wu, X.; Car, R., ‘The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water’, J. Chem. Phys. 2014, 141, 84502. 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(37) Jonchiere, R.; Seitsonen, A. P.; Ferlat, G.; Saitta, A. M.; Vuilleumier, R., ‘Van der Waals effects in ab initio water at ambient and supercritical conditions’, J. Chem. Phys. 2011, 135, 154503. (38) Zhang, C.; Wu, J.; Galli, G.; Gygi, F., ‘Structural and Vibrational Properties of Liquid Water from van der Waals Density Functionals’, J. Chem. Theory Comput. 2011, 7, 3054–3061. (39) Schmidt, J.; Vandevondele, J.; William Kuo, I. F.; Sebastiani, D.; Ilja Siepmann, J.; Hutter, J.; Mundy, C. J., ‘Isobaric-isothermal molecular dynamics simulations utilizing density functional theory: An assessment of the structure and density of water at near-ambient conditions’, J. Phys. Chem. B 2009, 113, 11959–11964. (40) Tuckerman, M. E.; Parrinello, M., ‘Integrating the Car-Parrinello equations. I. Basic integration techniques’, J. Chem. Phys. 1994, 101, 1302. (41) Grimme, S., ‘Accurate description of van der Waals complexes by density functional theory including empirical corrections’, J. Comput. Chem. 2004, 25, 1463–1473. (42) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H., ‘A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu’, J. Chem. Phys. 2010, 132, 154104. (43) Sit, P. H.-L.; Marzari, N., ‘Static and dynamical properties of heavy water at ambient conditions from first-principles molecular dynamics’, J. Chem. Phys. 2005, 122. (44) Blöchl, P. E., ‘Projector augmented-wave method’, Phys. Rev. B 1994, 50, 17953–17979. (45) Baerends, E. J.; Ziegler, T.; Autschbach, J.; Bashford, D.; Bérces, A.; Bickelhaupt, F. M.; Bo, C.; Boerrigter, P. M.; Cavallo, L.; Chong, D. P.; Deng, L.; Dickson, R. M.; Ellis, D. E.; van Faassen, M.; Fan, L.; Fischer, T. H.; Fonseca Guerra, C.; Ghysels, A.; Giammona, A.; van Gisbergen, S. J. A.; Götz, A. W.; Groeneveld, J. A.; Gritsenko, O. V.; Grüning, M.; Gusarov, S.; Harris, F. E.; van den Hoek, P.; Jacob, C. R.; Jacobsen, H.; Jensen, L.; Kaminski, J. W.; van Kessel, G.; Kootstra, F.; Kovalenko, A.; Krykunov, M. V.; van Lenthe, E.; McCormack, D. A.; Michalak, A.; Mitoraj, M.; Neugebauer, J.; Nicu, V. P.; Noodleman, L.; Osinga, V. P.; Patchkovskii, S.; Philipsen, P. H. T.; Post, D.; Pye, C. C.; Ravenek, W.; Rodríguez, J. I.; Ros, P.; Schipper, P. R. T.; Schreckenbach, G.; Seldenthuis, J. S.; Seth, M.; Snijders, J. G.; Solà, M.; Swart, M.; Swerhone, D.; te Velde, G.; Vernooijs, P.; Versluis, L.; Visscher, L.; Visser, O.; Wang, F.; Wesolowski, T. A.; van Wezenbeek, E. M.;


Page 28 of 31

Page 29 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Wiesenekker, G.; Wolff, S. K.; Woo, T. K.; Yakovlev, A. L., ‘Amsterdam Density Functional, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands.’, URL http://www.scm.com. Accessed 03/17. (46) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; van Gisbergen, S. J. A.; Fonseca Guerra, C.; Snijders, J. G.; Ziegler, T., ‘Chemistry with ADF’, J. Comput. Chem. 2001, 22, 931–967. (47) Fonseca Guerra, C.; Snijders, J. G.; te Velde, G.; Baerends, E. J., ‘Towards an order-N DFT method’, Theor. Chem. Acc. 1998, 99, 391. (48) van Lenthe, E.; Baerends, E. J.; Snijders, J. G., ‘Relativistic regular two-component Hamiltonians’, J. Chem. Phys. 1993, 99, 4597–4610. (49) Adamo, C.; Barone, V., ‘Toward reliable density functional methods without adjustable parameters: The PBE0 model’, J. Chem. Phys. 1999, 110, 6158–6170. (50) Klamt, A.; Schüürmann, G., ‘COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient’, J. Chem. Soc. Perkin Trans. 2 1993, 799–805. (51) Ceresoli, D. URL https://sites.google.com/site/dceresoli/ pseudopotentials Accessed 06/17. (52) Pérez, F.; Granger, B. E., ‘IPython: a System for Interactive Scientific Computing’, Computing in Science and Engineering 2007, 9, 21–29. (53) Duignan, T. J.; Marchenko, A., ‘Exatomic: A unified platform for computational chemists’, URL https://github.com/exa-analytics/exatomic. Accessed 06/17, doi:10.5281/zenodo.60053. (54) Pyykkö, P., ‘Year-2008 nuclear quadrupole moments’, Mol. Phys. 2008, 106, 1965–1974. (55) Rohatgi, A., ‘WebPlotDigitizer’, 2016, URL WebPlotDigitizer. Accessed 06/17, v3.10.

http://arohatgi.info/

(56) Ge, L.; Bernasconi, L.; Hunt, P., ‘Linking electronic and molecular structure: insight into aqueous chloride solvation’, Phys. Chem. Chem. Phys. 2013, 15, 13169–13183. (57) Sacco, A., ‘Structure and dynamics of electrolyte solutions. A NMR relaxation approach’, Chem. Soc. Rev. 1994, 23, 129. 29 ACS Paragon Plus Environment


(58) Perng, B.-C.; Ladanyi, B. M., ‘A dielectric theory of spin-lattice relaxation for nuclei with electric quadrupole moments’, J. Chem. Phys. 1998, 109, 676–684. (59) Zhang, C.; Pham, T. A.; Gygi, F.; Galli, G., ‘Communication: Electronic structure of the solvated chloride anion from first principles molecular dynamics’, J. Chem. Phys. 2013, 138, 181102. (60) Aidas, K.; Ågren, H.; Kongsted, J.; Laaksonen, A.; Mocci, F., ‘A quantum mechanics/molecular dynamics study of electric field gradient fluctuations in the liquid phase. The case of Na+ in aqueous solution’, Phys. Chem. Chem. Phys. 2013, 15, 1621–31.

for Table of Contents use only Quadrupolar NMR relaxation from ab-initio molecular dynamics: Improved sampling and cluster models vs. periodic calculations Adam Philips,† Alex Marchenko,† Lionel A. Truflandier,‡ and Jochen Autschbach ,∗ †

Institut des Sciences Moléculaires Université Bordeaux, CNRS UMR 5255 351 cours de la Libération 33405 Talence cedex ‡

Department of Chemistry University at Buffalo State University of New York Buffalo, NY 14260-3000, USA email: [email protected] †

1.0 0.8 0.6

ACF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 31

0.4 0.2 0.0 0

1

2

3

4

Time(ps)

The electric field gradient (EFG) at an analyte nucleus (35Cl– above) is computed for many configurations along an ab-initio molecular dynamics trajectory. Quadrupolar NMR relaxation rates are computed from the autocorrelation functions of these EFG tensor components. 30 ACS Paragon Plus Environment

1.0

Page 31 of 31

ACF

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43


0.8 0.6 0.4 0.2 0.0 0

1

2

ACS Paragon Plus Environment

Time(ps)

3

4

Quadrupolar NMR Relaxation from ab Initio Molecular Dynamics

Recommend Documents