Quadrupolar 14N NMR Relaxation from Force-Field and Ab-Initio

Nov 21, 2018 - Adam Philips , Alex Marchenko , Lucas Colucci Ducati , and Jochen Autschbach. J. Chem. Theory Comput. , Just Accepted Manuscript...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Quadrupolar 14N NMR Relaxation from Force-Field and Ab Initio Molecular Dynamics in Different Solvents Adam Philips,† Alex Marchenko,† Lucas C. Ducati,‡ and Jochen Autschbach*,† †

Department of Chemistry, University at Buffalo, State University of New York, Buffalo, New York 14260-3000, United States Department of Fundamental Chemistry Institute of Chemistry, University of São Paulo, Av. Prof. Lineu Prestes 748, São Paulo, SP 05508-000, Brazil



Downloaded via YORK UNIV on December 7, 2018 at 13:22:13 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Quadrupolar NMR spin relaxation rates and corresponding line widths were computed for the quadrupolar nucleus 14N for neat acetonitrile as well as for 1methyl-1,3-imidazole and 1-methyl-1,3,4-triazole in different solvents. Molecular dynamics (MD) was performed with forces from the Kohn−Sham (KS) theory (ab initio MD) and force-field molecular mechanics (classical MD), followed by KS electric field gradient (EFG) calculations. For acetonitrile the agreement of the 14N line width with experiment is very good. Relative line widths for the azole nitrogens are improved over simpler approximations used previously in conjunction with single-point calculations at the multiconfigurational self-consistent field level. Overall, the NMR line widths are computed within a factor of 2 of the experimental values, giving access to reasonable estimates both of the dynamic EFG variance in the solvated systems as well as the associated correlation times that determine the relaxation rates.

1. INTRODUCTION

Δν ∝

Nuclear magnetic resonance (NMR) relaxation is of great utility to study structure and dynamics of complex chemical systems.1−15 In solution, the observed line width is proportional to the NMR relaxation rate, which can be considerable for the spin relaxation of quadrupolar nuclei. When present, the interaction of the nuclear quadrupole moment with the electric field gradient (EFG) tensors caused by the surrounding nuclear and electronic charge distribution is assumed to dominate the spin relaxation and the observed line widths.1,16,17 Among the most important quadrupolar nuclear isotopes is 14N, with 99.6% abundance and ubiquity in nearly all branches of chemistry.18,19 Liquid phase line widths for 14N can often reach on the order of 103 Hz.18 While these line widths have traditionally been in the way of 14N NMR becoming a standard experimental probe, modern high-field spectrometers are able to overcome such obstacles. The quadrupolar line widths in solution (and the quadrupolar line shapes in solid-state NMR20,21) are therefore an additional source of information about the structure and bonding around a given nucleus and the dynamics of a system. In a seminal study,22 Jaszuński et al. investigated the solvent dependence of the 14N NMR line widths for a set of azoles, i.e. small organic heterocycles with well-defined rigid structures, in various solvents, both experimentally and theoretically. Overall, a strong solvent dependence of the line width was found. The theoretical approach of Jaszuński et al. consisted of single-point calculations with an implicit solvent model to obtain the EFG tensors at the nitrogen nuclei. The NMR line widths Δν were expressed via the following equation © XXXX American Chemical Society

η2 yzzij eQ yz2 1 3 2I + 3 ijj = + 1 j zjj V33zz τc 40 I 2(2I − 1) jjk 3 zz{k ℏ T2 {

(1)

Here, I is the nuclear spin, η is the EFG tensor asymmetry parameter, eQ is the electric quadrupole moment of the nucleus with Q being a cross section usually given in barn (10−28 m2), and V33 is the largest principal component of the EFG. Further, τc is a common correlation time associated with the relaxation process and the dynamics of the system. Under the assumption that τc is the same for all the observed nuclei and solvents, 14N line width ratios (relative line widths) are fully determined by V33 and η, the two independent parameters characterizing the traceless EFG tensor. Accordingly, Jaszuński et al. performed EFG calculations at the multiconfigurational self-consistent field (MCSCF) level with a dielectric medium as a solvent model and determined line width ratios relative to the methyl-bonded nitrogen NA in a set of compounds including the two azoles of Figure 1. This static theoretical approach gave increasing EFG magnitudes and reproduced the experimentally observed overall trend of increasing relative line widths, with an increasing dielectric constant of the solvent; but Jaszuński et al. also acknowledged the limitations of the static approach in describing solvated systems. Explicit solvent−solute interactions and their dynamics likely play a significant role in the quadrupolar relaxation. Indeed, the authors wrote22 “The approximation of identical and isotropic correlation times for molecular rotation with respect to chemically nonequivalent nitrogen atoms within the same Received: August 2, 2018 Published: November 21, 2018 A

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

representation. The longitudinal (1/T1) and transverse (1/ T2) relaxation rates are given by 1 = C Q G2,0(ω); T1

1 = C Q G2,1(ω) T2

(2a)

where ω is the Larmor frequency, and C is a constant depending on the nuclear quadrupole moment Q and the nuclear spin quantum number I: Q

Figure 1. Molecules studied in this work: (a) acetonitrile, (b) 1methyl-1,3-imidazole (short: “imidazole”), and (c) 1-methyl-1,3,4triazole (short: “triazole”).

CQ =

e 2Q 2(2I + 3) 40I 2(2I − 1)ℏ2

(2b)

Further, e is the elementary charge, and ℏ is the reduced Planck constant. The superscript Q denotes relaxation arising from the nuclear quadrupole-EFG interaction. The quantities G2,m(ω) are linear combinations of spectral densities g2,m(ω)

molecule is probably the weakest point in our analysis.”. Moreover, without calculating τc, absolute relaxation rates and NMR line widths cannot be obtained. The experimental data also showed notable outliers, for instance 1-methyl-1,3,4triazole in benzene, for which the line widths and ratios were notably greater than those in other nonpolar solvents. In addition, the line width ratios for the most polar solvents were strongly underestimated by the static theoretical approach. For these reasons, we deemed it necessary to investigate 14N quadrupolar relaxation theoretically with a dynamic approach. The present study features a molecular dynamics (MD) based approach which builds on pioneering force-field9 and ab initio MD studies23 of quadrupolar relaxation and follows a computational protocol established in two of our previous works.12,24 This approach is used to compute absolute 14N NMR line widths for neat acetonitrile as well as 1-methyl-1,3imidazole (short: “imidazole”) and 1-methyl-1,3,4-triazole (short “triazole”) shown in Figure 1 in various solvents. The molecular dynamics simulations are followed by EFG tensor calculations with explicit solvent along the trajectories to obtain the quadrupolar relaxation rates. Due to the unique trends in solvent effects discussed above, and the added sampling and error estimates gained from having two chemically equivalent nitrogens, the main focus of the current study is on the triazole. We use Kohn−Sham (KS) theorybased ab initio molecular dynamics (aiMD) to describe the solute−solvent interactions during the dynamics explicitly with quantum mechanics, as well as classical force field MD (FFMD) to compute the correlation times on the relevant time scales with acceptably small sampling errors. In accordance with previous work,24 sets of independent MD simulations for the same system are used for improved sampling in the imidazole and triazole relaxation calculations. This also allows an estimation of some of the sampling errors and proves to be crucial in approaching equivalence for the computed properties of the equivalent NB nitrogens in the triazole. Section 2 summarizes the theory and working equations for the computation of quadrupolar relaxation rates and provides computational details. Section 3 first presents a benchmark of the MD protocol for 14N NMR line widths, namely for neat acetonitrile, which is followed by a comparison of MD-based results for the azole compounds with the experimental and theoretical results of Jaszuński et al. Section 4 provides general conclusions.

G2,0(ω) = 4g2,2(2ω) + 4g2, −2(− 2ω) + g2,1(ω) + g2, −1(− ω)

(3a) G2,1(ω) = 2g2, −2(− 2ω) + 3g2, −1(− ω) + 2g2,1(ω) + 3g2,0(0)

(3b)

and, in turn, the spectral densities g2,m(ω) are the half-Fourier transforms of the time autocorrelation functions (ACFs) f 2,m(τ) = ⟨R*2,m(t)R2,m(t + τ)⟩ of the time-domain EFG rank-2 spherical tensor components R2,m(t), with the brackets denoting the ensemble average. In numerical simulations, the ensemble averages are replaced, for instance, by calculations of the ACF for discrete points along an MD trajectory. In the fast motion limit (extreme narrowing), where the ACFs decay to zero over a range τ with ωτ ≪ 1, the half-Fourier transforms reduce to a simple integration g2, m =

∫0



f2, m (τ )dτ

(4)

Under extreme narrowing and rotational isotropy, the g2,m(ω) and therefore the longitudinal and transverse relaxation rates become equal. An isotropic rate can be formulated accordingly as 1 1 1 = = = 2C Q τc⟨V (0)2 ⟩ T1 T2 Tiso

(5)

The quantity ⟨V(0) ⟩ = ∑m f l,m(0) is defined here as the total EFG variance ∑m⟨R*2,m(t)R2,m(t)⟩, and τc is a corresponding correlation time defined as 2

1 τc = ⟨V (0)2 ⟩

2



g2, m

(6)

m =−2

When there is sampling over many equivalent nuclei and, ideally, many independent trajectories, T1 and T2 may indeed become equivalent numerically, even in computationally demanding aiMD simulations.24 In this work, computed halfheight NMR line widths are computed from the transverse relaxation rate as defined in Equation 2a, under extreme narrowing, according to the definition

2. THEORY AND COMPUTATIONAL DETAILS The theoretical framework for nuclear spin relaxation was established by Bloembergen, Purcell, and Pound25 and Redfield.26 In the current work, we follow the formalism as presented by Spiess1 and by Cowan16 to obtain relaxation rates from time-dependent EFG data in the spherical tensor

Δν =

1 πT2

(7)

In the data tables of Section 3, ⟨V(0) ⟩ is reported in atomic units (au), τc in ps, and Δν in Hz. Note that one au of field gradient is approximately 9.717 × 1021 V/m2 and corresponds 2

B

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 2. 14N NMR line widths relative to NA for the triazole computed using various methods: “Jasz MCSCF” = Jaszuński et al. single-point MCSCF from ref 22, “Static KS” = nonhybrid KS single-point calculations, “saMD” = statically averaged ratios from FFMD clusters, ‘Dynamic (FFMD)’ = fully dynamic approach from classical dynamics, and ‘Expt’ = experimental results from ref 22. NB data for the triazole are averaged. Error bars for line width ratios are obtained by propagating the standard error of the absolute line widths from FFMD.

minimize coupling between electronic and nuclear degrees of freedom. A kinetic energy cutoff of 100 Ry was chosen for the plane-wave (PW) basis, and a fictitious electron mass (μ) of 380 au24,35,36 was used for the Car−Parrinello dynamics, with a time step of 6.0 au (0.145 fs). Random initial cell packing and 5 ps of pre-equilibration with FFMD in the NVT ensemble were performed with Tinker (see above). The aiMD simulations consisted of an additional 5 ps of equilibration in the NVT ensemble at 350 K using a three-chain Nosé−Hoover thermostat (90, 45, and 15 THz), followed by dynamics in the microcanonical (NVE) ensemble. The first 1 ps in the NVE ensemble was also taken to be part of equilibration. Total production times for “long” simulations were 50 ps for water, 60 ps for benzene, and 40 ps for neat acetonitrile. Sets of “short” simulations were each 20 ps for both the water and benzene systems. For ab initio water, an elevated simulation temperature of 350 K along with semiempirical van der Waals corrections37,38 to the energies and forces was employed to correct for the overstructured or “glassy” nature of water typically observed in aiMD simulations with nonhybrid functionals such as PBE.39,40 In a recent study, Willow et al. compared aiMD methods to polarizable force field models and attributed much of this overstructuring at the nonhybrid KS level to an overestimated polarizability.41 For consistency, the elevated temperature and dispersion corrections were used in the simulations with other solvents as well. EFGs for neat acetonitrile were calculated with the PAW module of QE,31,42 using periodic-boundary conditions and PAW pseudopotentials from the same pslibrary as used for the dynamics.32 EFGs for finite azole-solvent clusters were computed using the Amsterdam Density Functional (ADF)43−45 software package for snapshots at evenly spaced intervals of the simulations. These calculations utilized the PBE functional, nonperiodic boundary conditions, all-electron Slater-type orbital (STO) triple-ζ basis sets with two sets of polarization functions (TZ2P) for the atom centers in the analyte molecule, and a double-ζ basis set (DZ) for all atoms of the solvent molecules. Since only light elements are considered, and to minimize the computational cost for the large set of cluster calculations, relativistic effects were not

to a nuclear quadrupole coupling constant of 235.0 MHz per barn of nuclear quadrupole moment cross section. FFMD simulations were carried out with the Tinker molecular modeling package27 using the Optimized Potentials for Liquid Simulations all atom force field (OPLS-aa).28 This package provides parameters for acetonitrile and all solvents investigated herein. Water was modeled using the TIP3P parameter set.29 Force field parameters later developed for heterocycles by Jorgensen et al.30 were implemented to model the azoles. These force fields allow for intramolecular flexibility, including stretching, bending, and torsional motion. The periodic simulation cells are cubic, containing either neat acetonitrile, or an azole with a given number of solvent molecules. Cell dimensions were chosen such that the density of the system was equal to that of the pure solvent at ambient conditions. The specific parameters for each system were as follows: neat acetonitrile: 64 molecules, cell dimension = 18.13 Å; cyclohexane: 25 molecules, cell dimension = 18.88 Å; benzene: 32 molecules, cell dimension = 17.37 Å; acetone: 32 molecules, cell dimension = 16.54 Å; water: 64 molecules, cell dimension = 13.14 Å. In order to validate the choice of using a relatively small numbers of solvent molecules in the context of FFMD, a set of ten simulations was performed for the triazole and 500 benzene molecules. The relaxation results from this set is reported in Table S10 in the Supporting Information (SI). Within the standard error, the resulting line widths are equivalent to those generated from smaller scale simulations. For each solvent a set of randomly packed periodic boxes was optimized, followed by 50 ps of dynamics in the canonical (NVT) ensemble at 300 K. The time step was 1 fs, and the first 10 ps of each simulation were used for equilibration. The aiMD simulations were of the Car−Parrinello (CPMD) type and were carried out using the Quantum-Espresso (QE) software package.31 These simulations utilized ultrasoft pseudopotentials from pslibrary1.0.032 and the exchangecorrelation functional of Perdew, Burke, and Ernzerhof (PBE).33,34 To keep the computational effort manageable, the benzene simulations contained 16 solvent molecules, corresponding to a cell dimension of 13.58 Å. Hydrogens were replaced with deuterium in all CPMD simulations as to C

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Table 1. 14N NMR Relaxation Resultsa for Neat Acetonitrile from 100 ps FFMD and 40 ps aiMD, Using QE-PAW EFGs

treated (and shown by our static calculations to be negligible). Solute−solvent clusters were extracted from the periodic simulations using an open source tool developed in our group (version 0.3.11).46 The number of nearest neighbor solvent molecules extracted depends on the specific solvent and are as follows: cyclohexane: 12; benzene: 12; acetone: 16; water: 30. Additional bulk solvent effects were simulated via the Conductor like Screening Model (COSMO).47 Dielectric constants were the same as those used by Jaszuński et al.22 The dielectric constants and rigid sphere molecular radii respectively used for the COSMO model are as follows: cyclohexane: 2.024, 3.5 Å; benzene: 2.238, 3.28 Å; diethyl ether: 4.335, 3.46 Å; acetone: 20.70, 3.08 Å; methanol: 32.63, 2.53 Å; acetonitrile: 37.5, 2.76 Å; water: 78.54, 1.93 Å. The sampled simulation time and number of extracted clusters used for the line width calculations for each simulation are as follows: acetonitrile FFMD: 100 ps, 1000 snapshots; all azole FFMD: 40 ps, 400 snapshots; triazole short aiMD: 20 ps, 465 snapshots (benzene), 1390 snapshots (water); triazole long aiMD: 60 ps, 1368 (benzene), 50 ps, 1169 (water). In previous work24 as many as 4595 snapshots were used for simulations of comparable length; however, the relevant correlation times in the current systems are somewhat longer, and an increased time step between sampled configurations was shown to be valid (see Table S9 in the SI). For static EFG calculations of the azole compounds, reported in Tables 2 and 3 and Figure 2, geometries were optimized with PBE/TZ2P including the solvent model. Scalar relativistic effects−though negligible−were included via the use of the Zeroth Order Regular Approximation (ZORA) allelectron Hamiltonian.48 Comparisons were also made to hybrid functionals PBE049 and 3-parameter Becke Lee−Yang− Parr (B3LYP).50 The largest principal component of the EFG deviates less than 10% among these functionals, and PBE provides the best agreement with the calculated line width ratios of Jaszuński et al. Therefore, considering the large number of calculations required for the MD-based line width calculations, PBE was deemed sufficient for the EFG calculations. Tables S2−S4 in the SI provide the results obtained from nonrelativistic static calculations as well as from the other functionals.

FFMD aiMD

1 T2

τc

⟨V(0)2⟩

Δν

ΔνExptb

230.42 157.84

0.93 0.58

1.15 1.12

73(3) 50(7)

82(2) 82(2)

a

Calculated values are from this work. 64 molecules per simulation cell. EFGs were computed for 1000/1174 snapshot clusters from the FFMD/aiMD simulation. bExperimental data (“Expt”) are from ref 51. Standard error in the mean of 64 equivalent nitrogens and experimental error shown in parentheses. Simulation was performed at 300 K, and experimental results were determined at 298 K. Relaxation rates and line widths are in Hz, τc is in ps, and ⟨V(0)2⟩ is in atomic units.

provide an acceptable level of sampling.24 The agreement of the calculated line width with the experiment at 298 K51 is nearly quantitative for the FFMD simulation, with a relative deviation of only about 11%. The shorter (due to computational cost) aiMD trajectory gives not as well of an agreement with the experiment, but nonetheless both sets of calculations are superior to the factor of 2 to 3 deviations for quadrupolar relaxation rates reported for liquid heavy water calculated from aiMD simulations.23,24 With 14N being a main group element, the EFG tensors calculated from the KS theory with a nonhybrid functional are expected to be sufficiently accurate.19,52,53 The good agreement of the calculated NMR line width with the experiment also indicates that the FFMD describes the dynamics of the system well as far as the quadrupolar relaxation is concerned. The differences between the two sets of MD show up mainly in the autocorrelation times, while the EFG variances are essentially the same. It is likely that different force fields would also lead to similar variations of τc as found here in the comparison between FFMD and aiMD. Per Equation 5, in the isotropic extreme narrowing limit the relaxation rate is−up to a constant containing Q2−determined by the product of the correlation time τc and the total EFG variance ⟨V(0)2⟩. For 14N in acetonitrile, τc is under 1 ps and therefore less than the predicted correlation times of approximately 1.6 ps for quadrupolar 2H and 17O relaxation in liquid heavy water. FFMD gives τc larger by a factor of 1.5 or more than the correlation times for group 1 and 17 ions in water calculated from aiMD and FFMD.9,12,24,54 The results of single-point KS calculations of the 14N EFG for the azoles along with relative NMR line widths according to static ratios determined from Equation 1 are reported in Tables 2 and 3. The line width ratios from the KS calculations agree quite well with those reported by Jaszuński et al. based on MCSCF calculations, with a smooth trend of increasing line width relative to NA with increasing dielectric constant used for the solvent model. The calculated largest EFG principal components become less negative with increasing solvent polarity for all nitrogens, but the trend for NA is the weakest, causing the increasing Δν ratio. The asymmetry parameter, η at NA decreases with increasing solvent polarity for both the imidazole and triazole. In the triazole, η decreases for NB, while in the imidazole it is nearly zero with no consistent trend with respect to solvent polarity. An asymmetry parameter close to zero is consistent with the more symmetric bonding motif about NB in the imidazole. With a satisfactory benchmark for the MD-based methods applied to 14N in acetonitrile, and reproducibility of the line

3. RESULTS AND DISCUSSION As a proof of concept and baseline for the required sampling for 14N, a FFMD simulation of neat acetonitrile was performed. Since the extraction of solute−solvent clusters is somewhat arbitrary for a neat liquid and would unnecessarily reduce sampling, the 14N EFGs were computed for all 64 molecules in the simulation box under periodic boundary conditions, using the PAW formalism available in the QE code. We showed recently24 that the EFG tensors for solvated ions such as Na+ and Cl− calculated with a molecular cluster vs the periodic PAW approach are in good agreement. A specific approach may be chosen, depending on their relative advantages vs disadvantages, for instance sampling of equivalent nuclei within a periodic simulation cell versus the availability of a larger arsenal of computational tools for cluster-model level molecular electronic structure calculations. The acetonitrile relaxation data are summarized in Table 1. The computed NMR line width has a small standard error in the mean, due to the averaging over the 64 equiv nitrogens. Therefore, similar to our previous study of 2H and 17O relaxation in liquid heavy water, a single MD trajectory should D

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Table 2. V33, η, and Relative 14N NMR Line Widths for the Nitrogens of 1-Methyl-1,3-imidazole, from Static Nonhybrid KS Calculations with a Continuum Solvent Modela,c solvent gas-phase cyclohexane benzene diethyl ether acetone methanol acetonitrile water

label

V33

η

Δν/ ΔνNA

Δν/ b ΔνNA Jasz

Δν/ b ΔνNA Expt

NA NB NA NB NA NB NA NB NA NB NA NB NA NB NA NB

−0.562 −0.857 −0.528 −0.823 −0.524 −0.819 −0.505 −0.800 −0.485 −0.780 −0.482 −0.776 −0.482 −0.776 −0.478 −0.772

0.350 0.062 0.312 0.005 0.307 0.002 0.283 0.036 0.256 0.072 0.251 0.080 0.249 0.078 0.246 0.089

1.00 2.24 1.00 2.36 1.00 2.37 1.00 2.45 1.00 2.53 1.00 2.55 1.00 2.55 1.00 2.57

1.00 2.46 1.00 2.55 1.00 2.57 1.00 2.64 1.00 2.72 1.00 2.73

1.00 1.74 1.00 2.40 1.00 2.33 1.00 2.56 1.00 4.13

1.00 2.74

1.00 2.71

approach for reducing certain sampling errors by exploring more of the configurational space, while at the same time providing some sampling error estimates. (Systematic errors from the electronic structure calculations, and sampling errors arising from neglecting events with time scales at and beyond the simulation times, are obviously not covered.) Our most detailed data set is that for the triazole in Table 5. The best agreement with experiment is for the triazole 14N line widths in water, with deviations of about 10% or less. For the other three solvents that were investigated, the deviations reach 20−40%. Despite the fact that the NMR line widths from these simulations do not agree as well with the experimental data as those from the triazole in water and the neat acetonitrile simulations, it is worthwhile reiterating that deviations between calculated and experimental quadrupolar relaxation rates (or line widths) by factors of 2 to 3 were previously deemed quite satisfactory,12,23,24,54 i.e. when the correct order of magnitude was obtained. This reflects the difficulties of describing both the dynamic behavior and the electronic structures correctly and with accurate sampling. An estimation of the uncertainty in our present data is also provided in the form of the standard error in the mean for the line widths computed from each set of 10 simulations with a given solvent. The considerable standard errors for the NB, as well as the numerical inequivalency of the two NB nitrogens of the triazole (listed separately as NB and NB′) are representative of the challenge of obtaining sufficient sampling within the finite simulations. In all cases there is good equivalence established for the total EFG variance ⟨V(0)2⟩ between NB and NB′, i.e. the two nitrogens experience the same average magnitude and spread of the EFG during the simulations. The computed correlation time τc is more limited by the sampling. Simulation times much greater than the correlation time are necessary in order to obtain smooth ACFs and well converged τc. In the case of water as the solvent, the computed correlation times are on the order of 5 ps, which is over 10% of the sampled FFMD simulation time. In general, simulations must be long enough to sample the relevant dynamics and obtain accurate correlation times but (even in the case of FFMD) not so long as to unnecessarily increase the quantity and expense of the subsequent properties calculations. A single 1 ns FFMD simulation was run for the triazole in water (see Table S7 in the SI), and the computed correlation times in that case are indeed no greater than 5.5 ps, suggesting that there are diminishing returns for simulations of this length. One important finding that can be extracted from Tables 4 and 5 is that the assumption of a common τc for the different nitrogens underlying the static approach is not generally valid, as the correlation times vary differently with the solvent and differ markedly among the nitrogen atoms in the same molecule in some cases. The final column in Tables 4 and 5 is denoted Δν/ΔνNA saMD (the subscript stands for “statically averaged MD”) and gives the N:NA line width ratios based on Equation 1 but averaged over the full set of cluster configurations from the dynamics that were used to calculate the relaxation rates. This approach accounts for the effects of explicit solvent, and it also averages over many structures from the dynamics, but, as in the static calculations, there is the assumption of equivalent correlation times for the unique nitrogens. When compared to the static calculations with implicit solvent, the inclusion of explicit solvent has a strong effect for the triazole, increasing the ratios significantly. The NB/NA line width ratios in the limited data

a

Reported values are from this work unless otherwise specified. “Jasz” = MCSCF calculations by Jaszuński et al., ref 22. Experimental data (“Expt”) are from the same reference. Principal EFG component, V33, is in atomic units. cNA is the methyl-bound nitrogen. b

Table 3. V33, η, and Relative 14N NMR Line Widths for the Nitrogens in 1-Methyl-1,3,4-triazole from Static Nonhybrid KS Calculations with a Continuum Solvent Modela,c solvent gas-phase cyclohexane benzene diethyl ether acetone methanol acetonitrile water

label

V33

η

Δν/ ΔνNA

Δν/ b ΔνNA Jasz

Δν/ b ΔνNA Expt

NA NB NA NB NA NB NA NB NA NB NA NB NA NB NA NB

−0.546 −0.990 −0.508 −0.953 −0.503 −0.948 −0.481 −0.927 −0.458 −0.903 −0.454 −0.899 −0.454 −0.899 −0.449 −0.894

0.366 0.279 0.305 0.355 0.297 0.365 0.257 0.412 0.212 0.465 0.205 0.476 0.205 0.475 0.197 0.487

1.00 3.22 1.00 3.56 1.00 3.61 1.00 3.84 1.00 4.11 1.00 4.16 1.00 4.16 1.00 4.22

1.00 3.24 1.00 3.46 1.00 3.49 1.00 3.66 1.00 3.86 1.00 3.88

1.00 2.73 1.00 3.98 1.00 2.76 1.00 4.60 1.00 5.68

1.00 3.88

1.00 5.97

a

Reported values are from this work unless otherwise specified. Minor differences in the data between the chemically equivalent NB appeared due to not constraining the more or less arbitrary rotation of the methyl group in the structure optimizations. Averaged data shown here. See Table S1 in the SI for the individual values. b“Jasz” = MCSCF calculations by Jaszuński et al., ref 22. Experimental data (“Expt”) are from the same reference. Principal EFG component, V33, is in atomic units. cNA is the methyl-bound nitrogen.

width ratios established for the static method used by Jaszuński et al., we proceed to the MD-based relaxation calculations of the azole systems. Tables 4 and 5 collect 14N relaxation data for the imidazole and triazole, respectively, obtained from sets of ten 40 ps (production time) FFMD simulations with various solvents. Previous work in this group24 has shown averaging across multiple independent simulations to be an effective E

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Table 4. Summary of Relaxation Resultsa from Sets of Ten 40 ps FFMD Simulations and Subsequent KS EFG Calculations on Clusters of the Imidazole and Explicit Solventc solvent benzene water

label A

N NB NA NB

1 T2

τc

⟨V(0)2⟩

145.56 268.92 135.83 482.27

1.28 1.19 1.60 2.55

0.46 1.00 0.39 0.73

Δν 46 86 43 154

Δ(Δν)

(1.00) (1.87) (1.00) (3.58)

11 12 12 18

ΔνExptb 62 149 81 220

(1.00) (2.40) (1.00) (2.71)

%dev

Δν/ΔνNA saMD

−25.3 −42.6 −46.6 −30.2

1.00 2.21 1.00 2.22

a Reported values are from this work unless otherwise specified. EFGs were computed for 400 snapshot clusters from each simulation. Δ(Δν) is the standard error in the mean of the computed line widths from the set of ten simulations for a given solvent. bExperimental data (‘Expt’) are from ref 22. Relative NMR line widths are reported in parentheses for both computed and experimental values. Relaxation rates and line widths are in Hz, τc is in ps, and ⟨V(0)2⟩ is in atomic units. cNA is the methyl-bound nitrogen.

Table 5. Summary of Relaxation Resultsa from Sets of Ten 40 ps FFMD Simulations and Subsequent KS EFG Calculations on Clusters of the Triazole and Explicit Solventc solvent

label

cyclohexane

NA NB NB′ NA NB NB′ NA NB NB′ NA NB NB′

benzene

acetone

water

1 T2

τc

⟨V(0)2⟩

78.08 254.63 270.78 145.83 494.65 501.25 121.61 537.46 846.98 193.87 1376.65 1178.82

1.07 0.88 0.82 2.35 1.68 1.55 2.12 2.01 2.72 5.30 5.39 4.49

0.31 1.37 1.36 0.27 1.33 1.32 0.24 1.22 1.23 0.16 1.07 1.08

Δν 25 81 86 46 157 160 39 171 270 62 438 375

Δ(Δν)

(1.00) (3.24) (3.44) (1.00) (3.41) (3.48) (1.00) (4.38) (6.92) (1.00) (7.06) (6.05)

7 14 18 10 20 29 6 29 51 10 83 76

ΔνExptb 41 112 112 59 235 235 47 216 216 69 412 412

(1.00) (2.73) (2.73) (1.00) (3.98) (3.98) (1.00) (4.60) (4.60) (1.00) (5.97) (5.97)

%dev

Δν/ΔνNA saMD

−39.4 −27.6 −23.0 −21.3 −33.0 −32.1 −17.6 −20.8 24.8 −10.6 6.4 −8.9

1.00 5.31 5.33 1.00 4.98 4.97 1.00 4.50 4.47 1.00 7.09 7.18

a Reported values are from this work unless otherwise specified. EFGs were computed for 400 snapshot clusters from each simulation. Δ(Δν) is the standard error in the mean of the computed line widths from the set of ten simulations for a given solvent. bExperimental data (‘Expt’) are from ref 22. Relative NMR line widths are reported in parentheses for both computed and experimental values. Relaxation rates and line widths are in Hz, τc is in ps, and ⟨V(0)2⟩ is in atomic units. cNA is the methyl-bound nitrogen.

Table 6. Summary of Relaxation Dataa from Sets of Ten 40 ps FFMD Simulations and Subsequent KS EFG Calculations on the Triazole after Removal of All Solvent Moleculesc solvent (dynamics only) benzene

water

label A

N NB NB′ NA NB NB′

1 T2

τc

⟨V(0)2⟩

189.66 535.64 544.41 460.61 1470.58 1522.17

2.41 1.71 1.58 5.46 4.87 4.10

0.35 1.42 1.42 0.34 1.40 1.40

Δν 60 171 173 147 468 485

(1.00) (2.85) (2.88) (1.00) (3.18) (3.30)

Δ(Δν) 13 21 32 34 66 150

ΔνExptb 59 235 235 69 412 412

(1.00) (3.98) (3.98) (1.00) (5.97) (5.97)

%dev 1.7 −27.2 −26.4 113.0 13.6 17.7

Reported values are from this work unless otherwise specified. EFGs were computed for 400 snapshot clusters from each simulation. Δ(Δν) is the standard error in the mean of the computed line widths from the set of ten simulations for a given solvent. bExperimental data (‘Expt’) are from ref 22. Relative NMR line widths are reported in parentheses for both computed and experimental values. Relaxation rates and line widths are in Hz, τc is in ps, and ⟨V(0)2⟩ is in atomic units. cNA is the methyl-bound nitrogen. a

that it performs poorly for the nonpolar solvents but acceptably for the polar solvents. The EFG is known to be dominated by short-range interactions mainly with the first solvation shell or, in the case of a molecule, by the solute dynamics itself. In order to directly assess the relative importance of these interactions on the computed NMR line widths, EFGs and relaxation parameters were computed for the same FFMD snapshots as used for benzene and water in Table 5 but with all solvent molecules removed in the KS EFG calculations. What remains is the intramolecular influence on the EFG, but with the intramolecular dynamics influenced by the solvent interactions occurring in the FFMD. The results are reported in Table 6

set for the imidazole are underestimated by this method but not severely in comparison with experiment. For the triazole, NA ΔνNB/ΔνsaMD is severely overestimated for the nonpolar solvents. For the polar solvents, however, the line width ratios obtained via explicit treatment of the solvent in the dynamics are in better agreement with experiment than those from the simplified static approach. Figure 2 visually compares the relative line width results for the triazole from the various methods. Within the (considerable) error bars estimated for the full MD-based line width calculations, there is agreement with the experiments. The static method falls short for the polar solvents but performs acceptably for the nonpolar solvents. The “saMD” approach gives the opposite behavior, in F

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Table 7. Summary of Relaxation Resultsa from aiMD Simulations and Subsequent KS EFG Calculations on Clusters of the Triazole and Explicit Solventc solvent

1 T2

τc

⟨V(0)2⟩

219.58 678.32 526.88 182.97 1319.98 2054.77

2.21 2.50 2.26 2.38 3.90 6.33

0.40 1.42 1.41 0.30 1.15 1.17

label A

benzene

N NB NB′ NA NB NB′

water

Δν 70 216 168 58 420 654

Δνexpb

(1.00) (3.09) (2.40) (1.00) (7.24) (11.28)

59 235 235 69 412 412

(1.00) (3.98) (3.98) (1.00) (5.97) (5.97)

a

Reported values are from this work unless otherwise specified. Total simulation time for benzene was 60 ps with EFGs computed for 1368 snapshot configurations. Total simulation time for water was 50 ps with EFGs calculated for 1169 snapshot configurations. bExperimental data (‘Expt’) are from ref 22. Relative NMR line widths are reported in parentheses for both computed and experimental values. Relaxation rates and line widths are in Hz, τc is in ps, and ⟨V(0)2⟩ is in atomic units. cNA is the methyl-bound nitrogen.

Table 8. Summary of Relaxation Resultsa from Sets of Ten 20 ps aiMD Simulations and Subsequent KS EFG Calculations on Clusters of the Triazole and Explicit Solventc solvent benzene

water

label A

N NB NB′ NA NB NB′

1 T2

τc

⟨V(0)2⟩

124.26 572.29 532.46 136.05 879.18 558.72

1.39 1.87 1.65 1.73 3.25 2.43

0.41 1.42 1.41 0.31 1.16 1.16

Δν 40 182 169 43 280 178

(1.00) (4.61) (4.29) (1.00) (6.51) (4.14)

Δ(Δν) 9 35 42 9 29 30

Δνexpb 59 235 235 69 412 412

(1.00) (3.98) (3.98) (1.00) (5.97) (5.97)

a Reported values are from this work unless otherwise specified. EFGs were computed for 465 (benzene) and 1169 (water) snapshot clusters from each corresponding simulation. (Δ(Δν) is the standard error in the mean of the computed line widths from the set of ten simulations for a given solvent. bExperimental data (‘Expt’) are from ref 22. Relative NMR line widths are reported in parentheses for both computed and experimental values. Relaxation rates and line widths are in Hz, τc is in ps, and ⟨V(0)2⟩ is in atomic units. cNA is the methyl-bound nitrogen.

Figure 3. Oxygen−oxygen (left) and nitrogen(B)−hydrogen (right) pair correlation function for FFMD simulations of the triazole in water. The single line represents a 1 ns simulation, and the bands represent ±2 standard deviations around the mean across ten 40 ps simulations.

would otherwise produce τc in excess of ca. 5 ps, thus biasing the average over the “short” trajectories toward small τc. It is interesting to note, however, that the relative τc for NA and NB (with the exception of NB′ in the long simulation) is in agreement with a suggestion made by Kongsted et al. that the pyridine-type nitrogen should have a correlation time 1.6 times that of the pyrrole-type nitrogen.55 Figures 3 and 4 provide the oxygen−oxygen and NBhydrogen pair correlation functions (PCF) for all long and short FFMD and aiMD simulations of the triazole. PCFs for sets of short simulations are represented as bands with widths corresponding to the average across the simulations plus or minus 2 standard deviations. In terms of average structure, there is much more variation across simulations in aiMD than in classical dynamics. Correspondingly, in the aiMD, although

and reveal a significant solvent dependence for the NB/NA relative line widths, especially for water. The most significant contributor to this is an increase in the total EFG variance (more so for NA) as compared to the fully solvated model. Tables 7 and 8 summarize results computed from aiMD simulations of the triazole in water and in benzene. In both solvents, the total EFG variance is systematically slightly larger in the structures taken from aiMD as compared to FFMD. Overall, weaker equivalence is established between NB and NB′ for the correlation time and consequently the absolute line widths which indicates the lesser sampling. The 20 ps aiMD simulations summarized in Table 8 produce systematically smaller correlation times than the longer simulations of Table 7. Here, the 20 ps production runs−especially for water−are not long enough to converge the ACFs in simulations that G

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 4. Oxygen−oxygen (left) and nitrogen(B)−hydrogen (right) pair correlation function for aiMD simulations of the triazole in water. The single line represents a 50 ps simulation, and the bands represent ±2 standard deviations around the mean across ten 20 ps simulations.

Figure 5. Average velocity (left) and rotational (right) autocorrelation functions for solvent molecules from the longest available simulations of the triazole in water and benzene. The corresponding correlation time is noted in the legend. FFMD is excluded from the VACFs because the print frequency of the simulation was on a too slow time scale.

rotational correlation times for both solvents in FFMD coincide quite well with the corresponding EFG correlation times from the same simulations. It therefore appears that the slower rotational motion of the solvent around a solute is the major contributor to the quadrupolar relaxation for our samples. The RACFs reveal measurable differences in the dynamics from aiMD and FFMD, but the trend is not consistent for the two solvents. Even with the elevated simulation temperatures in the aiMD, the reorientational motion for water is slightly slower as compared to FFMD. The opposite is found for benzene. Additional plots of these ACFs as well as some dynamic structural analysis for the triazole itself can be found in the SI. Ultimately, for the azoles a simulation time frame of 40−60 ps appears to be at the lower end of time scales required for adequate sampling to estimate the correlation times reliably across the different solvents. Therefore, since well-parametrized force fields are available for these molecules and solvents, sets of independent simulations of this length using force field dynamics, in combination with KS calculations of the EFG tenors, are a good compromise between sampling and computational cost. However, apart from the obvious sampling errors evidenced by the nonequivalency of NB and NB′, the single long aiMD simulations do a good job of producing reasonable magnitudes for the different Δν. The line widths of the triazole NB in benzene are consistently underestimated at all investigated levels of theory. However, the correlation times for NB and NB′ are shorter than that of NA in the FFMD simulations with benzene, while the

there is significant statistical overlap, there is weaker agreement with respect to the structure about NB versus NB′. It is also interesting to note that the first peak in the N−H PCF (corresponding to direct coordination of solvent hydrogen atoms to the basic nitrogens) appears at a shorter distance and is resolved in the aiMD. In the FFMD, the water hydrogen bonding to the basic triazole nitrogens appears as a shoulder in the first main peak which corresponds to the vicinal intramolecular hydrogen atoms. Short-distance hydrogen bonding compared to experiment has been noted before for MD at the nonhybrid KS level.56 While there is clearly more variation in the structure in aiMD versus FFMD, this does not translate to a much larger standard error of the computed NMR line widths. In fact, in water, the relative standard error for the NB Δν is a bit greater with FFMD. We note that the PCFs represent descriptions of the average structure from the simulations and do not include the dynamic information necessary for computing relaxation rates/NMR line widths. To obtain a more complete picture of the dynamics in these simulations, linear velocity and rotational ACFs (VACF and RACF, respectively) were computed for the available simulations. These are vector autocorrelations, the former utilizing the center-of-mass velocity vector and the latter utilizing a vector normal to the plane of each molecule. Figure 5 provides plots of these functions as averages over the solvent molecules in the simulations. The VACFs decay on a much faster time scale than the autocorrelation of the EFG fluctuations in the azole, while the rotational correlation times are of the same order of magnitude. Namely, the H

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

about 20%, with the ratios accordingly differing only little. It indeed appears that for the azoles the main approximation in the static approach is the assumption of identical correlation times.

opposite trend is seen in the results obtained from aiMD. This suggests that aiMD may capture solute−solvent interactions differently than FFMD. One potential candidate for such an interaction in benzene is π-stacking assisted by dispersion interactions. The latter are treated as a semiempirical add-on in the aiMD simulations, while nonbonding interactions in the OPLS-aa force field were obtained via fits to Monte Carlo simulations for 34 organic liquids (but not including benzene).28 Visual inspection of the simulations did not reveal that aiMD produces more or fewer π-stacked structures than FFMD. Per Equation 5, the magnitude of the EFG fluctuations, along with the EFG tensor autocorrelation time, drives the quadrupolar relaxation process. In order to quantify the importance of the fluctuation magnitude, we consider ⟨V(0)2⟩ versus the square of the unsigned-mean, ⟨|V(0)|⟩2, of the total EFG tensor. The latter can be thought of as a dynamic analogue to the EFG used in the static approach (however, in the case of the final column in Tables 4 and 5, averaging is done across the ratios computed for the individual clusters). Differences between ⟨V(0)2⟩ and ⟨|V(0)|⟩2 may cancel to a large degree when ratios are taken for different nitrogens. Tables 9 and 10 provide these quantities for the nitrogens in

4. SUMMARY AND CONCLUSIONS Absolute NMR line widths due to nuclear quadrupole relaxation of the 14N nuclear spins in neat acetonitrile and two C−N heterocyclic (azole) compounds in several solvents were computed via classical and ab initio molecular dynamics, followed by EFG calculations on snapshot configurations at the KS DFT level. These results can be directly compared to experiment and show improvement over simpler models which ignore unique solute−solvent interactions and/or do not explicitly determine the correlation times governing the relaxation. Given that the accessible simulation times are longer for FFMD vs aiMD, averaged results over sets of FFMD simulations provided the most extensive sampling and the most accurate results, the best of which being the 14N line width of neat acetonitrile and triazole in water, with percent errors on the order of 10% with respect to the experimental NMR line widths. The NB line widths for the triazole in benzene were consistently underestimated by the calculations, and the trends in correlation times across the nitrogens changed systematically based on whether clusters were taken from FFMD or aiMD. These results suggest there may be subtle interactions in the triazole-benzene system that contribute to the relaxation. Further investigation is required in order to determine the exact nature of these interactions. While aiMD appears to sample a wider structural range, classical dynamics with well parametrized force fields−if available−provides a good compromise for the description of the subtle solute−solvent interactions, the sampling time and computational cost, and facilitates the use of multiple independent trajectories. Fundamentally, the benefit gained from the dynamics-based approach taken here as opposed to a static ratios-based approach is 3-fold. First, absolute line widths are determined that can be directly compared to experiment. Second, it allows the computation of unique correlation times for chemically unique environments. Third, the mean square fluctuations in the EFG are fully considered.

Table 9. Total EFG Variance ⟨V(0)2⟩ (Mean-Square EFG Sum) vs the Square-of-Unsigned-Mean EFG ⟨|V(0)|⟩2, Averaged for Nitrogens in the Triazole from Sets of Ten 40 ps FFMD Simulations in Various Solventsa,b solvent

label

cyclohexane

NA NB NB′ NA NB NB′ NA NB NB′ NA NB NB′

benzene

acetone

water

⟨V(0)2⟩ 0.31 1.37 1.36 0.27 1.33 1.32 0.24 1.22 1.23 0.16 1.07 1.08

(1.00) (4.39) (4.36) (1.00) (4.82) (4.82) (1.00) (5.09) (5.10) (1.00) (6.72) (6.79)

⟨|V(0)|⟩2 0.26 1.14 1.13 0.23 1.11 1.09 0.20 1.02 1.00 0.13 0.92 0.92

(1.00) (4.34) (4.29) (1.00) (4.77) (4.69) (1.00) (5.08) (4.98) (1.00) (6.87) (6.86)



a ⟨V(0)2⟩ and ⟨|V(0)|⟩2 are in atomic units. bRatios with respect to the methyl-bound nitrogen NA are in parentheses.

Table 10. Total EFG Variance ⟨V(0) ⟩ (Mean-Square EFG Sum) and the Corresponding Static Square-of-Mean EFG ⟨V(0)⟩2, Averaged for Nitrogens in the Triazole from Sets of Ten 20 ps aiMD Simulations in Benzene and Watera,b solvent

label

benzene

NA NB NB′ NA NB NB′

water

⟨V(0)2⟩ 0.41 1.42 1.41 0.31 1.16 1.16

(1.00) (3.46) (3.45) (1.00) (3.77) (3.75)

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00807. EFG and relaxation data for different nonrelativistic static data, plot of corresponding line width ratios with different functionals, relaxation data for each of the independent short trajectories (FFMD and aiMD), relaxation data from a single 1 ns FFMD simulation of triazole in water, triazole relaxation data from FFMD simulations with 500 benzene molecules, convergence of V33 with respect to the number of nearest neighbor benzenes, velocity and rotational ACFs for the triazole, comparisons of PCFs from FFMD and aiMD simulations of triazole in water, distance of nearest neighbor hydrogen to NB over the course of 50 ps aiMD simulation of triazole in water, histogram of NB−NB bond distances for sets of 40 ps FFMD simulations of the triazole in different solvents, ACFs of NB−NB bond

⟨|V(0)|⟩2 0.35 1.18 1.19 0.26 1.00 1.00

ASSOCIATED CONTENT

S Supporting Information *

2

(1.00) (3.41) (3.43) (1.00) (3.88) (3.86)

⟨V(0)2⟩ and ⟨|V(0)|⟩2 are in atomic units. bRatios with respect to the methyl-bound nitrogen NA are in parentheses. a

the triazole from the classical and ab initio simulations, respectively. As expected, the mean-square is systematically greater than the square-unsigned-mean, but not by more than I

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation



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. (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) Cowan, B. Nuclear Magnetic Resonance and Relaxation; Cambridge University Press: Cambridge, UK, 2005. (17) Kowalewski, J.; Mäler, L. Nuclear spin relaxation in liquids: Theory, experiments, and applications; Taylor & Francis: New York, 2006. (18) Mann, B. E. The common nuclei. In NMR and the Periodic Table; Harris, R. K., Mann, B. E., Eds.; Academic Press: London, 1978; pp 87−105. (19) O’Dell, L.; Schurko, R.; Harris, K.; Autschbach, J.; Ratcliffe, C. Interaction Tensors and Local Dynamics in Common Structural Motifs of Nitrogen: A Solid-state 14N NMR and DFT Study. J. Am. Chem. Soc. 2011, 133, 527−546. (20) Perras, F. A.; Widdifield, C. M.; Bryce, D. L. QUEST− QUadrupolar Exact SofTware: A fast graphical program for the exact simulation of NMR and NQR spectra for quadrupolar nuclei. Solid State Nucl. Magn. Reson. 2012, 45-46, 36−44. (21) Wasylishen, R. Characterization of NMR tensors via experiment and theory. In Calculation of NMR and EPR Parameters. Theory and Applications; Kaupp, M., Bühl, M., Malkin, V. G., Eds.; WileyVCH: Weinheim, 2004; pp 433−447, DOI: 10.1002/ 3527601678.ch27. (22) Jaszuński, M.; Mikkelsen, K. V.; Rizzo, A.; Witanowski, M. A Study of the Nitrogen NMR Spectra of Azoles and their Solvent Dependence. J. Phys. Chem. A 2000, 104, 1466−1473. (23) 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. (24) Philips, A.; Marchenko, A.; Truflandier, L. A.; Autschbach, J. Quadrupolar NMR relaxation from ab-initio molecular dynamics: Improved sampling and cluster models vs. periodic calculations. J. Chem. Theory Comput. 2017, 13, 4397−4409. (25) Bloembergen, N.; Purcell, E.; Pound, R. Relaxation effects in nuclear magnetic resonance absorption. Phys. Rev. 1948, 73, 679−715. (26) Redfield, A. G. On the Theory of Relaxation Processes. IBM J. Res. Dev. 1957, 1, 19−31. (27) Ponder, J. W. TINKER - Software Tools for Molecular Design. https://dasher.wustl.edu/tinker/ (accessed 09/18). (28) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (29) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (30) McDonald, N. A.; Jorgensen, W. L. Development of an AllAtom Force Field for Heterocycles. Properties of Liquid Pyrrole, Furan, Diazoles, and Oxazoles. J. Phys. Chem. B 1998, 102, 8049− 8059. (31) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (32) Dal Corso, A. Pseudopotentials periodic table: From H to Pu. Comput. Mater. Sci. 2014, 95, 337−350. (33) Perdew, J. P.; Burke, K.; Wang, Y. Generalized gradient approximation for the exchange-correlation hole of a many-electron

distances for sets of 40 ps FFMD simulations of the triazole in different solvents, ACFs of Vzz (Cartesian) and V33 (principal component) of the EFG at NB for sets of 40 ps FFMD simulations of the triazole in different solvents (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: jochena@buffalo.edu. ORCID

Adam Philips: 0000-0002-5742-6817 Alex Marchenko: 0000-0001-8915-0328 Lucas C. Ducati: 0000-0002-6539-4325 Jochen Autschbach: 0000-0001-9392-877X Funding

We acknowledge the Center for Computational Research (CCR) at the University at Buffalo for providing computational resources. This work has been supported by grant CHE1560881 from the National Science Foundation. L.C.D. is grateful for financial support from FAPESP (2015/08541-6 and 2017/17750-3). We thank the Center for Computational Research (CCR) in Buffalo for providing computational resources, in particular those made available via grant MRI1724891 from the National Science Foundation. Notes

The authors declare no competing financial interest.



REFERENCES

(1) Spiess, H. W. Rotation of Molecules and Nuclear Spin Relaxation. In NMR Basic Principles and Progress; Diehl, P., Fluck, R. K. E., Eds.; Springer: Berlin, 1978; Vol. 15, pp 55−214, DOI: 10.1007/978-3-642-66961-3_2. (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) Garrett, B. B.; Denison, A.; Rabideau, S. 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, 81, 1204−1221. (6) Hanni, M.; Lantto, P.; Vaara, J. Nuclear spin relaxation due to chemical shift anisotropy of gas-phase 129Xe. 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: J

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation system. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 16533− 16539. (34) Perdew, J. P.; Burke, K.; Ernzerhof, M. Reply to Comment on generalized gradient approximation made simple. Phys. Rev. Lett. 1998, 80, 891. (35) 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. (36) 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. (37) Grimme, S. Accurate description of van der Waals complexes by density functional theory including empirical corrections. J. Comput. Chem. 2004, 25, 1463−1473. (38) 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. (39) 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, 204510. (40) Lin, I.-C.; Seitsonen, A. P.; Tavernelli, I.; Rothlisberger, U. Structure and Dynamics of Liquid Water from ab Initio Molecular Dynamics−Comparison of BLYP, PBE, and revPBE Density Functionals with and without van der Waals Corrections. J. Chem. Theory Comput. 2012, 8, 3902−3910. (41) Willow, S. Y.; Zeng, X. C.; Xantheas, S. S.; Kim, K. S.; Hirata, S. Why Is MP2-Water “Cooler” and “Denser” than DFT-Water? J. Phys. Chem. Lett. 2016, 7, 680−684. (42) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (43) 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. et al. Amsterdam Density Functional, SCM, Theoretical Chemistry; Vrije Universiteit: Amsterdam, The Netherlands. https://www.scm. com (accessed 03/17). (44) 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. (45) Fonseca Guerra, C.; Snijders, J. G.; te Velde, G.; Baerends, E. J. Towards an order-N DFT method. Theor. Chem. Acc. 1998, 99, 391. (46) Duignan, T. J.; Marchenko, A. Exatomic: A unified platform for computational chemists. https://github.com/exa-analytics/exatomic (accessed Nov 30, 2018); DOI: 10.5281/zenodo.60053. (47) 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, 2, 799−805. (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) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields. J. Phys. Chem. 1994, 98, 11623−11627. (51) Tiffon, B.; Ancian, B. 14N nuclear magnetic resonance study of the rotation of acetonitrile in n-alkanes. J. Chem. Phys. 1982, 76, 1212−1216. (52) Schwerdtfeger, P.; Pernpointner, M.; Nazarewicz, W. Calculation of Nuclear Quadrupole Coupling Constants. In Calculation of NMR and EPR Parameters. Theory and Applications; Kaupp, M., Bühl, M., Malkin, V. G., Eds.; Wiley-VCH: Weinheim, 2004; pp 279−291, DOI: 10.1002/3527601678.ch17.

(53) Autschbach, J.; Zheng, S.; Schurko, R. W. Analysis of Electric Field Gradient Tensors at Quadrupolar Nuclei in Common Structural Motifs Concepts. Concepts Magn. Reson., Part A 2010, 36A, 84−126. (54) 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. (55) Møgelhøj, A.; Aidas, K.; Mikkelsen, K. V.; Kongsted, J. Solvent Effects on the Nitrogen NMR Shielding and Nuclear Quadrupole Coupling Constants in 1-Methyltriazoles. Chem. Phys. Lett. 2008, 460, 129−136. (56) 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, 084502.

K

DOI: 10.1021/acs.jctc.8b00807 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX