Comparing Vibrationally Averaged Nuclear Shielding Constants by

Article pubs.acs.org/JPCA

Comparing Vibrationally Averaged Nuclear Shielding Constants by Quantum Diffusion Monte Carlo and Second-Order Perturbation Theory Yee-Hong Ng and Ryan P. A. Bettens* Department of Chemistry, National University of Singapore, 3 Science Drive 3, 117543 Singapore S Supporting Information *

ABSTRACT: Using the method of modified Shepard’s interpolation to construct potential energy surfaces of the H2O, O3, and HCOOH molecules, we compute vibrationally averaged isotropic nuclear shielding constants ⟨σ⟩ of the three molecules via quantum diffusion Monte Carlo (QDMC). The QDMC results are compared to that of second−order perturbation theory (PT), to see if second-order PT is adequate for obtaining accurate values of nuclear shielding constants of molecules with large amplitude motions. ⟨σ⟩ computed by the two approaches differ for the hydrogens and carbonyl oxygen of HCOOH, suggesting that for certain molecules such as HCOOH where big displacements away from equilibrium happen (internal OH rotation), ⟨σ⟩ of experimental quality may only be obtainable with the use of more sophisticated and accurate methods, such as quantum diffusion Monte Carlo. The approach of modified Shepard’s interpolation is also extended to construct shielding constants σ surfaces of the three molecules. By using a σ surface with the equilibrium geometry as a single data point to compute isotropic nuclear shielding constants for each descendant in the QDMC ensemble representing the ground state wave function, we reproduce the results obtained through ab initio computed σ to within statistical noise. Development of such an approach could thereby alleviate the need for any future costly ab initio σ calculations.

■

INTRODUCTION Over the years, nuclear magnetic resonance (NMR) has emerged as one of the most important and powerful spectroscopic techniques for the characterization of molecules. As it is the only method presently that is able to provide detailed atomic structures in solution, much effort has been devoted to the accurate predictions of nuclear shielding constants from first principles. Unfortunately, to date, the quality of ab initio calculated σ is not yet at the stage where unambiguous assignments of experimental isotropic nuclear shielding constants of individual atoms can be performed (for an excellent review on this topic, see ref 1). Many factors influence the accuracy of ab initio calculated σ. These include electron correlation,2−5 conformers,6−8 temperature,9−13 solvent,6,14−20 relativistic,21,22 and zero-point vibrational effects.23−27 For the case of the zero-point vibrational effects, secondorder perturbation theory23 is the most widely used approach for calculating vibrational contributions to NMR properties. The method though generally fails at describing molecules with large amplitude motions. Consequently, for molecules with large displacements away from equilibrium, second-order perturbation theory may be inadequate for obtaining accurate zero−point vibrational corrections to isotropic nuclear shielding constants. Originally formulated by Anderson,28,29 quantum diffusion Monte Carlo was first applied to electronic structure problems © 2016 American Chemical Society

and subsequently used in the studies of molecular vibrations and intermolecular modes in molecular clusters.30−32 Being numerically exact, QDMC allows one to solve for the ground state properties of any quantum system to within statistical noise. The favorable scaling of computational costs with system size is another reason that it has become one of the most used methods in all of today’s quantum chemical calculations. However, for accurate results, the QDMC method requires an accurate potential energy surface. Not only has Collins’ PES construction method been used to study numerous reactions,33−42 it has also been successfully applied to obtain vibrationally averaged rotational constants of bound systems.43 Combining (1) quantum diffusion Monte Carlo, (2) quantum chemical ab initio calculations, and (3) Collin’s potential energy surface construction method, one should then in principle, be able to compute the most accurate of vibrationally averaged isotropic nuclear shielding constants. In this work, we evaluate vibrationally averaged isotropic nuclear shielding constants by means of second−order perturbation theory and quantum diffusion Monte Carlo. The results of the two approaches are compared, to see if second− order PT is adequate for obtaining accurate values of nuclear shielding constants. In the later part of the paper, we extend the Received: January 3, 2016 Revised: February 1, 2016 Published: February 2, 2016 1297

DOI: 10.1021/acs.jpca.6b00036 J. Phys. Chem. A 2016, 120, 1297−1306

Article

The Journal of Physical Chemistry A

Figure 1. Flowchart of the random walk algorithm with the discrete and descendant weighting schemes.

approach of modified Shepard’s interpolation to the construction of σ surfaces. Since the descendant weighting scheme used within QDMC to compute ⟨σ⟩ requires heavy computation of shielding constants of various conformations of the molecules, we investigate the possibility of replacing

costly ab initio isotropic nuclear shielding constants calculations with simple evaluations from such surfaces. Such a development could greatly increase the efficiency of evaluating ⟨σ⟩ via the QDMC approach, thereby allowing for practical and accurate computations of such quantities in the future. 1298


Article


Figure 2. Flowchart of the potential energy surface growing algorithm.

■

COMPUTATIONAL METHODS All quantum chemical calculations were performed using the Gaussian 09 package.44 Quantum Diffusion Monte Carlo. The random walk algorithm45 with the discrete46 and descendants47 weighting schemes used for the quantum diffusion Monte Carlo calculations have been well documented and is summarized in Figure 1. In this work, we randomly displace the 3N − 6 normal coordinates of the molecule. Defining the normal coordinate QK to be

reduced masses of the normal coordinates are then given by the Wilson G-matrix, its elements being 3N

Gbv =

J=1

∑ LJK J=1

1 ∂Q b ∂Q v mJ ∂xJ ∂xJ

(2)

Specifically, the diagonal entries have the meaning of reciprocal reduced masses along the respective normal modes of vibration while the off-diagonal elements correspond to pairwise kinetic couplings between them. With the definition of normal coordinates in eq 1,

3N

QK =

∑

mJ (xJ − xJ ,eq)

3N

(1)

Gbv =

where x denotes the 3N Cartesian coordinates of all N atoms, xeq the corresponding 3N Cartesian coordinates at equilibrium geometry, m the masses of the respective atoms, and LJK the components of the matrix made up of the normalized column eigenvectors of the mass weighted Hessian (MWH), the

∑ J=1

1 (LJb mJ )(LJv mJ ) = mJ

3N

∑ (LJbLJv) J=1

(3)

Since the eigenvectors of a symmetric matrix (MWH) are orthogonal and can be further normalized such that they are orthonormal, the column vectors of L are orthonormal. Hence 1299


Article

The Journal of Physical Chemistry A 3N

Gbv =

∑ (LJbLJv) = δbv

hw[B(j)] =

(4)

J=1

where δbv is the Kronecker delta. Thus, under the definition of Qk of eq 1, the Wilson G-matrix is the identity matrix and the reduced masses of all normal coordinates Qk equal 1. The components of the MWH are hereby given to be 1 ∂ 2U mimj ∂xi ∂xj

MWHij =

RMSP[B(j)] = γj 2 =

(5)

respectively, the algorithm of Figure 2 was then used to further grow the PESs. The hw function weighs a particular configuration against the number of existing data points in that space, as well as the number of replicas that visit that region during the course of the simulation. Conformational space in which many replicas visit indicate areas of chemical importance, and the hw function essentially chooses a point which reside in a region of chemical importance containing sparse data points to be added to the existing potential energy surface. The γ2 function on the other hand evaluates the sum of the weighted differences between the computed potential by individual data point’s truncated Taylor series and the modified Shepard’s interpolated potential. Configurations with large disagreements between them reside in areas whereby the PES is inaccurate, and the RMSP method chooses such configurations to be added to the existing PES. Alternating between the h-weight and RMSP methods then increase the number of data points in the PES, thereby improving the agreement between the modified Shepard’s interpolated potential and the actual ab initio energy for any arbitrary geometry. Convergence of the PES was checked by performing 70 QDMC calculations for every 25 data points added to the potential energy surface. We diffuse the initial 1000 replicas, starting at the equilibrium geometry, for 10000 steps, with a step−size of 1 atomic unit (a.u.). ER, being the QDMC calculated zero−point energy (ZPE) of the system, was averaged from τ = 500 onward, at 100 time−step intervals, in each run of the QDMC simulations. The final computed ZPE is then the mean of the 70 ER, with the estimated error being 3ς/√N, where N = 70, and ς is the standard deviation of the 70 energies. The method of descendants weights was used to generate geometries of a vibrational averaged population of molecules. The first set of descendants was followed after 10000 steps, for ndeswtsstep, such that ∼100 unique ancestors remain from a starting population of ∼1000. The next set of tracking begins at τ = 10000 + ndeswtsstep, with the procedure repeating until 300 sets of unique ancestors that characterizes the square of the ground state nuclear wave function |ψ0|2 have been generated. Using the gauge invariant atomic orbitals (GIAO) method,4 ab initio calculations were performed on the ∼30000 geometries, and the vibrationally averaged isotropic nuclear shielding constant ⟨σ⟩ is given by the simple average of all the ancestors’ σ. The estimated error here is similarly 3ς/√N, where N ∼ 30000, and ς is the standard deviation of the ∼30000 ancestors’ σ. Second−Order Perturbation Theory. We performed geometry optimizations followed by frequency calculations at the B3LYP/aug-cc-pVTZ level of theory. Isotropic nuclear shielding constants σ were calculated at the same level of theory by the GIAO method using the B3LYP/aug-cc-pVTZ optimized geometries.

(6)

where the summation runs over all data points Nd (importance sampled configurations). The normalized weight function wj[B(j)] takes the form vj[B(j)] Nd ∑k = 1 vk[B(j)]

(7)

while the unnormalized weight function used by Ischtwan and Collins49 is vj[B(j)] =

1 ∥B − B(j)∥2p

(8)

To guarantee convergence of the potential V(B) for an Natomic system, the value of 2p is required to be greater than the larger of 3N − 6 and the order of the Taylor series expansion. In the present work, 2p is chosen to be 18. In constructing the respective potential energy surfaces, we start off with a single point of equilibrium geometry. They are then further initialized with 2*(3N − 6) data points where Qk separately take on values of plus and minus root-mean-square (RMS) of the harmonic oscillators’ normal coordinates. These serve as walls that prevent replicas from diffusing into regions of very high V at the start of the simulation. The RMS normal coordinates are hereby given to be RMS(Bk ) =

ςk =

1 2ςk

(9)

2πωkμk ℏ

Nd

∑ wk[B(j)][Tk[B(j)] − V (B)]2 (12)

Nd

wj[B(j)] =

(11)

k=1

∑ wj[B(j)]Tj[B(j)] j=1

⎧ ⎫ vj[B(n)] ⎪ ⎪ ⎨ N ⎬ ⎪ ⎪ d n = 1/ n ≠ j ⎩ ∑k = 1 vk[B(j)] ⎭ Nd

∑

and

where U is the electronic energy of the molecule. The use of the algorithm of Figure 1 requires an efficient and accurate way for computing the potential of each replica at each time−step of the simulation. As mentioned in the introduction, we made use of modified Shepard’s interpolation48 for the construction of the potential energy surfaces of the various molecules under study. Under this formalism, the PES is represented as a sum of weighted Taylor series expansions T, truncated after second−order, about importance sampled configurations. As with the random displacements of the molecules’ internal degrees of freedom, this interpolation is also carried out in normal coordinates space. Consequently, the potential V at any configuration B is V (B) =

1 Nd − 1

(10)

where ωk and μk are respectively, the harmonic frequency and reduced mass of the kth normal mode of vibration. Defining the h-weight and RMSP functions to be 1300


Article

The Journal of Physical Chemistry A Zero-point vibrational corrections (ZPVCs) were computed using second-order perturbation theory with the molecules in their respective equilibrium geometries as the unperturbed systems. Following the works of Ruud et al.,23,50 the ZPVC to σ of a particular nucleus is given by σZPVC =

1 4

∑ K

1 ∂ 2σ 1 − ωK ∂Q K2 4

∑ K

1 ∂σ ωK 2 ∂Q K

∑ L

FKLL ωL (13)

where QK is the normal coordinate K with harmonic frequency ωK, and FKLL is the third derivative of the electronic energy E with respect to the normal coordinates QK, QL, QL. The vibrationally averaged isotropic nuclear shielding constant ⟨σ⟩ is then computed as ⟨σ ⟩ = σeq + σZPVC

(14)

The ZPVCs were evaluated at the B3LYP/aug-cc-pVTZ level of theory. All cubic force constants, first and second derivatives of the shielding constants with respect to the normal coordinates were evaluated numerically with a step-size of 2.5 atomic units (a.u) in the normal coordinates. It should be noted that, based on the above definition of the normal coordinates, a step-size of 2.5 a.u in the normal coordinates corresponds to displacements of the atoms such that 3N

∑ (xJ − xJ ,eq)2

■

J=1

≈ 0.05Å (15)

RESULTS AND DISCUSSION QDMC vs Second−Order PT ⟨σ⟩. As seen in Figures 3, 6, and 8, at the B3LYP/aug-cc-pVTZ level of theory, the potential

Figure 4. ⟨σ⟩ of H2O as a function of the number of interpolated data points. All quantum chemical calculations were performed at the B3LYP/aug-cc-pVTZ level of theory. The error bars indicate thrice the standard error of the mean of the σ of the ∼30000 ancestors. The magenta, green, and red horizontal lines indicate respectively, σeq, σexp, and the second-order PT evaluated ⟨σ⟩.

Figure 3. ZPE of H2O as a function of the number of interpolated data points. All quantum chemical calculations were performed at the B3LYP/aug-cc-pVTZ level of theory. The error bars indicate thrice the standard error of the mean of ER from 70 runs of the QDMC simulation. The green and red horizontal lines indicate respectively, the ZPE calculated under the harmonic approximation, and the current best estimate for the ZPE.51

Figure 5. ZPE (cm−1) of H2O as a function of the number of interpolated data points. All quantum chemical calculations were performed at the CCSD(T)/aug-cc-pVTZ level of theory. The error bars indicate thrice the standard error of the mean of ER from 70 runs of the QDMC simulation. The green and red horizontal lines indicate, respectively, the ZPE calculated under the harmonic approximation and the current best estimate for the ZPE.51

energy surfaces of the H2O, O3 and HCOOH molecules converge with 125, 50, and 275 data points in the interpolated PES respectively. For the case of the HCOOH molecule, missing sampled points can be found at Nd = 50, 150, 225, and 325. Addition of data points of high potential and steep negative sloping gradients had resulted in the formation of holes in the PES. The modified Shepard’s interpolation consequently returned negative values for the potentials of

many replicas. The resultant ER and error bars were wildly off the Nd = 1 point, and further data points had to be added to patch up the potential energy surface before the simulation could produce sensible results. It should be hereby mentioned that the drilling of holes is especially serious in higher 1301


Article


Figure 6. ZPE of O3 as a function of the number of interpolated data points. All quantum chemical calculations were performed at the B3LYP/aug-cc-pVTZ level of theory. The error bars indicate thrice the standard error of the mean of ER from 70 runs of the QDMC simulation. The green horizontal line indicate the ZPE calculated under the harmonic approximation while the current best estimate for the zero-point energy, 1429,52 is not shown in the plot as it widely skews the graph to the bottom.

Figure 8. ZPE of HCOOH as a function of the number of interpolated data points. All quantum chemical calculations were performed at the B3LYP/aug-cc-pVTZ level of theory. The error bars indicate thrice the standard error of the mean of ER from 70 runs of the QDMC simulation. The green horizontal line indicate the ZPE calculated under the harmonic approximation.

Table 1. Calculated and Experimental ZPE (cm−1) and ⟨σ⟩ (ppm)a Property ZPE ⟨σO⟩ ⟨σH1⟩ ⟨σH2⟩ ZPE ZPE ⟨σO(center)⟩ ⟨σO(end1)⟩ ⟨σO(end2)⟩ ZPE ⟨σC1⟩ ⟨σH2⟩ ⟨σO3⟩ ⟨σH4⟩ ⟨σO5⟩

QDMC

PT

H2O-B3LYP/aug-cc-pVTZ PES 4610 ± 12 − 312.8 ± 0.8 313.0 30.61 ± 0.07 30.56 30.60 ± 0.06 30.56 H2O-CCSD(T)/aug−cc−pVTZ PES 4636 ± 11 − O3-B3LYP/aug−cc−pVTZ PES 1583 ± 7 − −1091 ± 3 −1089 −1588 ± 6 −1584 −1588 ± 6 −1584 HCOOH-B3LYP/aug-cc-pVTZ PES 7360 ± 11 − 17.4 ± 0.4 17.5 23.16 ± 0.05 23.05 106 ± 1 106 25.26 ± 0.07 25.14 −116 ± 1 −114

Experiment 4634 323.6 30.10 30.10 4634 1429 −743 −1309 −1309 − − − − − −

a

QDMC values shown are the results of the simulations with maximum number of data points in the potential energy surfaces. Experimental values for the ZPE and σ of the HCOOH molecule unavailable from literature.

With only the equilibrium geometry, the QDMC calculated ZPEs match the analytic solution of the harmonic oscillators, as they should. However, as shown in Table 1, upon convergence of the PESs, the final values for the ZPEs of the H2O and O3 molecules still disagree with the experimental values. This disparity is attributed to the inadequacies of the B3LYP functional, or DFT in general. When the level of theory is increased to CCSD(T)/aug-cc-pVTZ (Figure 5), the QDMC obtained ZPE matches the experimental ZPE for the H2O molecule. Since the main purpose of this work is to not to obtain theoretical shielding constants that are of spectroscopic quality but to compare between QDMC and second−order PT computed ⟨σ⟩, from here on out, we make no further references to that of the experimental values. Figures 4, 7, and 10 show the QDMC and second−order PT calculated ⟨σ⟩ as a function of Nd in the interpolated PES. It can

Figure 7. ⟨σ⟩ of O3 as a function of the number of interpolated data points. All quantum chemical calculations were performed at the B3LYP/aug-cc-pVTZ level of theory. The error bars indicate thrice the standard error of the mean of the σ of the ∼30000 ancestors. The magenta, and red horizontal lines indicate respectively, σeq, and the second-order PT evaluated ⟨σ⟩. The experimental values of −74353 for Ocenter and −130953 for Oend are not shown in the plots as they widely skew the graph to the top.

dimensional surfaces, as a single data point can greatly alter the behavior of the PES in many different directions. 1302


Article


highlight the deficiency of second−order PT. For certain molecules such as HCOOH that has a large amplitude motion (internal OH rotation), vibrationally averaged isotropic nuclear shielding constants of experimental quality may only be obtainable with the use of more sophisticated and accurate methods, such as quantum diffusion Monte Carlo. σ Surfaces. In the previous section, we performed ab initio calculations on the ancestors population that characterizes |ψ0,N|2 of the respective molecules. However, as mentioned previously, each set of descendants weights has ∼100 ancestors and the averaging over 300 sets of them means that ab initio calculations have to be performed on ∼30000 unique geometries at each sampled Nd. Sampling the PES at 17 different Nd for the HCOOH molecule (11 for the H2O and O3 molecules), the total number of ab initio computations required to generate the graphs of Figures 4, 7, and 10 are then on the order of 510000! Noting the computational intensity of the process, we thereby seek an alternate method for computing the isotropic nuclear shielding constant of any molecule at any arbitrary geometry. We extend the method of modified Shepard’s interpolation to the evaluation of isotropic nuclear shielding constants. At each of the data points of the constructed potential energy surface, we perform ab initio calculations to obtain the value of σ at that geometry. Using the two and three point finite differences method, we further evaluated the first and second derivatives of σ with respect to the normal coordinates Qj. The isotropic nuclear shielding constant of any arbitrary configuration B can then be similarly written as a weighted average of the second−order Taylor series about all Nd data points:

Figure 9. Nuclei labels of formic acid. The white, gray, and red spheres represent H, C, and O atoms, respectively.

Nd

σ(B) =

∑ wj[B(j)]T jσ[B(j)] (16)

j=1

T jσ[B(j)] = σ0[B(j)] + [B − B(j)]T GB′ σ(j) +

1 [B − B(j)]T GB″(σj)[B − B(j)] 2

(17)

where σ0[B(j)] is the value of the shielding constant at B(j), σ σ G′B(j) , and G″B(j) are respectively, the gradient vector and second derivative matrix of the same property with respect to the normal coordinates. In eq 17, the superscript T denotes the transpose. We hereby use the constructed shielding constants surfaces to compute σ of the same ancestors in the sets of descendants of H2O and O3. Results of the study are as shown in Figures 11 and 12. H2O and O3. In Figures 11 and 12, the green lines result from the evaluations of the isotropic nuclear shielding constants of the ancestors by the surface with corresponding number of data points along the x−axis. Initially, it was thought that, large inaccuracies would surface in the computed shielding constants with only a single point in the σ surface. However, as seen in both figures, the blue and green points at Nd = 1 in each of the six graphs virtually overlap. This agreement (within statistical noise) between the blue and green lines continue as the number of data points in the σ surface is increased. This result then prompts the use of the σ surface with only the equilibrium geometry as a single data point to evaluate the shielding constants of all the ancestors at the various number of interpolated data points in the potential energy surface. As with previously, the resultant red lines in the two figures agree with

Figure 10. ⟨σ⟩ of HCOOH as a function of the number of interpolated data points. All quantum chemical calculations were performed at the B3LYP/aug-cc-pVTZ level of theory. The error bars indicate thrice the standard error of the mean of the σ of the ∼30000 ancestors. The magenta and red horizontal lines indicate, respectively, σeq and the second-order PT evaluated ⟨σ⟩.

be seen that, even with the equilibrium geometry as a single data point, ⟨σO⟩ of H2O, as well as ⟨σ⟩ of all nuclei of O3 and HCOOH (Figure 9) differ from the corresponding equilibrium values. General deviations of ⟨σ⟩ from σeq with various Nd in the PES highlight the importance of zero−point vibrational averaging in any theoretical calculation of isotropic nuclear shielding constants. Examining QDMC and second−order PT calculated ⟨σ⟩ revealed that the approaches yield essentially the same results for all nuclei except ⟨σH2⟩, ⟨σH4⟩, and ⟨σO5⟩ of the HCOOH molecule. Taking into account of the estimated errors, Table 1 shows that the QDMC and second−order PT calculated ⟨σ⟩ of the three nuclei differ by 0.06, 0.05, and 1 ppm, respectively. Although the differences are small, the results nevertheless 1303


Article


Figure 11. Comparisons between ab initio and σ surface evaluated isotropic nuclear shielding constants of H2O. The blue line is the result of ab initio computed shielding constants of the ancestors, while the green line is due to evaluations by the surface with corresponding number of data points in the σ surface. Finally, the red line is the consequence of evaluations by the surface with only 1 data point: the point of equilibrium geometry.

Figure 12. Comparisons between ab initio and σ surface evaluated isotropic nuclear shielding constants of O3. The blue line is the result of ab initio computed shielding constants of the ancestors, while the green line (overlaps with blue line) is due to evaluations by the surface with corresponding number of data points in the σ surface. Finally, the red line is the consequence of evaluations by the surface with only 1 data point: the point of equilibrium geometry.

the ab initio calculated results of the blue line to within statistical noise. With a sample size of ∼330000 unique configurations, the modified Shepard’s interpolated σ surface with the equilibrium geometry as a single data point can thus be said to be adequate for obtaining isotropic nuclear shielding constants that are of comparable quality to that of a proper ab initio calculation. HCOOH. We hereby construct shielding constant surfaces of the HCOOH molecule with the equilibrium geometry as a single data point in the σ landscape. The surface is then used to compute the isotropic nuclear shielding constants of the ancestors in the sets of descendants. For the case of HCOOH, the ancestors were not stored, and similar sets of ancestors had to be regenerated for the cause. Results of the study are as shown in Figure 13. It can be seen from Figure 13 that the agreement between the blue and red lines is not as perfect as in Figures 11 and 12. As mentioned, the previous ancestors of the formic acid in Figure 13 were not stored, and the geometries of the σ surface evaluated ancestors are not identical to those that made up the blue line. Nonetheless, the σ surface evaluated shielding constants still agree with the ab initio results to within statistical noise. This further supports the claim that the modified Shepard’s interpolated σ surface with the equilibrium geometry as a single data point is sufficient for obtaining

isotropic nuclear shielding constants that are of comparable quality to that of a proper ab initio calculation.

■

CONCLUSION Using the modified Shepard’s interpolated PES to calculate the potentials of various configurations of replicas, the method of descendants weighting was implemented within the QDMC approach to obtain vibrationally averaged isotropic nuclear shielding constants of the H2O, O3, and HCOOH molecules. In the cases of the hydrogens and carbonyl oxygen of HCOOH, the QDMC and PT evaluated ⟨σ⟩ disagreed. This implied the deficiency of second−order PT. For certain molecules such as HCOOH where large amplitude motions are present (internal OH rotation), vibrationally averaged isotropic nuclear shielding constants of experimental quality may only be obtainable with the use of more sophisticated and accurate methods, such as quantum diffusion Monte Carlo. The approach of modified Shepard’s interpolation had also been successfully extended to compute isotropic nuclear shielding constants of various geometries of H2O, O3, and HCOOH. Results showed that, a σ surface with the equilibrium geometry as a single data point is sufficient for computing shielding constants that are of comparable quality to that of a proper ab initio calculation. Development of such an approach 1304



■

could thereby alleviate the need for any future costly ab initio σ calculations.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b00036. Values of all QDMC obtained ⟨σ⟩ with different number of data points in the interpolated PES (PDF)

■

REFERENCES

(1) Mulder, F. A. A.; Filatov, M. NMR chemical shift data and ab initio shielding calculations: Emerging tools for protein structure determination. Chem. Soc. Rev. 2010, 39, 578−590. (2) Gauss, J. Effects of electron correlation in the calculation of nuclear magnetic resonance chemical shifts. J. Chem. Phys. 1993, 99, 3629−3643. (3) Bühl, M.; Gauss, J.; Hofmann, M.; Schleyer, P. v. R. Decisive electron correlation effects on computed 11B and 13C NMR chemical shifts. Application of the GIAO-MP2 method to boranes and carbaboranes. J. Am. Chem. Soc. 1993, 115, 12385−12390. (4) Cheeseman, J. R.; Trucks, G. W.; Keith, T. A.; Frisch, M. J. A comparison of models for calculating nuclear magnetic resonance shielding tensors. J. Chem. Phys. 1996, 104, 5497−5509. (5) Sundholm, D.; Gauss, J.; Ahlrichs, R. The electron correlation contribution to the nuclear magnetic shielding tensor of the hydrogen molecule. Chem. Phys. Lett. 1995, 243, 264−268. (6) Atieh, Z.; Allouche, A. R.; Lazariev, A.; Van Ormondt, D.; Graveron-Demilly, D.; Aubert-Frécon, M. DFT calculations of 1H chemical shifts, simulated and experimental NMR spectra for sarcosine. Chem. Phys. Lett. 2010, 492, 297−301. (7) Dumez, J.-N.; Pickard, C. J. Calculation of NMR chemical shifts in organic solids: Accounting for motional effects. J. Chem. Phys. 2009, 130, 104701. (8) Gortari, I. D.; Portella, G.; Salvatella, X.; Bajaj, V. S.; Van Der Wel, P. C. A.; Yates, J. R.; Segall, M. D.; Pickard, C. J.; Payne, M. C.; Vendruscolo, M. Time averaging of NMR chemical shifts in the MLF peptide in the solid state. J. Am. Chem. Soc. 2010, 132, 5993−6000. (9) Sundholm, D.; Gauss, J.; Schäfer, A. Rovibrationally averaged nuclear magnetic shielding tensors calculated at the coupled-cluster level. J. Chem. Phys. 1996, 105, 11051−11059. (10) Jameson, C. J.; Keith Jameson, A.; Cohen, S. M.; Parker, H.; Oppusunggu, D.; Burrell, P. M.; Wille, S. Temperature dependence of the 15N and 1H nuclear magnetic shielding in NH3. J. Chem. Phys. 1981, 74, 1608−1612. (11) Jameson, C. J.; De Dios, A. C.; Jameson, A. K. The 31P shielding in phosphine. J. Chem. Phys. 1991, 95, 9042−9053. (12) Jameson, C. J.; De Dios, A. C.; Jameson, A. K. Nuclear magnetic shielding of nitrogen in ammonia. J. Chem. Phys. 1991, 95, 1069−1079. (13) Jameson, C. J. Gas-phase NMR spectroscopy. Chem. Rev. 1991, 91, 1375−1395. (14) Han, W.-G.; Jalkanen, K. J.; Elstner, M.; Suhai, S. Theoretical study of aqueous N-acetyl-L-alanine N′-methylamide: Structures and Raman, VCD, and ROA spectra. J. Phys. Chem. B 1998, 102, 2587− 2602. (15) Poon, C.-D.; Samulski, E. T.; Weise, C. F.; Weisshaar, J. C. Do bridging water molecules dictate the structure of a model dipeptide in aqueous solution? J. Am. Chem. Soc. 2000, 122, 5642−5643. (16) Malkin, V. G.; Malkina, O. L.; Steinebrunner, G.; Huber, H. Solvent effect on the NMR chemical shieldings in water calculated by a combination of molecular dynamics and density functional theory. Chem. - Eur. J. 1996, 2, 452−457. (17) Cossi, M.; Crescenzi, O. Different models for the calculation of solvent effects on 17O nuclear magnetic shielding. J. Chem. Phys. 2003, 118, 8863−8872. (18) Mennucci, B.; Martínez, J. M.; Tomasi, J. Solvent effects on nuclear shieldings: Continuum or discrete solvation models to treat hydrogen bond and polarity effects? J. Phys. Chem. A 2001, 105, 7287− 7296. (19) Cui, Q.; Karplus, M. Molecular properties from combined QM/ MM methods. 2. chemical shifts in large molecules. J. Phys. Chem. B 2000, 104, 3721−3743. (20) Pfrommer, B. G.; Mauri, F.; Louie, S. G. NMR chemical shifts of ice and liquid water: The effects of condensation. J. Am. Chem. Soc. 2000, 122, 123−129. (21) Cromp, B.; Carrington, T., Jr.; Salahub, D. R.; Malkina, O. L.; Malkin, V. G. Effect of rotation and vibration on nuclear magnetic resonance chemical shifts: Density functional theory calculations. J. Chem. Phys. 1999, 110, 7153−7159.

Figure 13. Comparisons between ab initio and σ surface evaluated isotropic nuclear shielding constants of HCOOH. The blue line is the result of ab initio computed σ of the ancestors while the red line is the consequence of evaluations by the surface with only 1 data point: the point of equilibrium geometry.

■

Article

AUTHOR INFORMATION

Corresponding Author

*(R.P.A.B.) E-mail: [email protected]. Telephone: +65 6516 2846. Fax: +65 6779 1691. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thanks the National University of Singapore’s support from the Academic Research Fund, Grant Number R143-000-549-112. The authors also thank the Centre for Computational Science and Engineering for the use of their computers. 1305


Article


(42) Bettens, R. P. A.; Collins, M. A. Potential Energy Surfaces and Dynamics for the Reactions between C(3P) and H+3 (1A′1). J. Chem. Phys. 1998, 108, 2424−2433. Cited by ref 39. (43) Bettens, R. P. A. Bound State Potential Energy Surface Construction: Ab Initio Zero-Point Energies and Vibrationally Averaged Rotational Constants. J. Am. Chem. Soc. 2003, 125, 584−587. (44) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09, Revision A.1; Gaussian Inc.: Wallingford CT, 2009. (45) Gregory, J. K.; Clary, D. C. Calculations of the Tunneling Splittings in Water Dimer and Trimer Using Diffusion Monte Carlo. J. Chem. Phys. 1995, 102, 7817−7829. (46) McCoy, A. Diffusion Monte Carlo approaches for investigating the structure and vibrational spectra of fluxional systems. Int. Rev. Phys. Chem. 2006, 25, 77−107. (47) Suhm, M. A.; Watts, R. O. Quantum Monte Carlo Studies of Vibrational States in Molecules and Clusters. Phys. Rep. 1991, 204, 293−329. (48) Jordan, M. J. T.; Thompson, K. C.; Collins, M. A. Convergence of Molecular Potential Energy Surfaces by Interpolation: Application to the OH+H2 → H2O+H Reaction. J. Chem. Phys. 1995, 102, 5647− 5657. (49) Ischtwan, J.; Collins, M. A. Molecular Potential Energy Surfaces by Interpolation. J. Chem. Phys. 1994, 100, 8080−8088. (50) Åstrand, P.-O.; Ruud, K.; Taylor, P. R. Calculation of the Vibrational Wave Function of Polyatomic Molecules. J. Chem. Phys. 2000, 112, 2655−2667. (51) Barone, V. Vibrational Zero-Point Energies and Thermodynamic Functions Beyond the Harmonic Approximation. J. Chem. Phys. 2004, 120, 3059−3065. (52) Herzberg, G. Electronic Spectra and Electronic Structure of Polyatomic Molecules; D. Van Nostrand Company, Inc.: New York, 1966. (53) Solomon, I. J.; Keith, J. N.; Kacmarek, A. J.; Raney, J. K. Additional Studies Concerning the Existence of “O3F2. J. Am. Chem. Soc. 1968, 90, 5408−5411.

(22) Minaev, B.; Vaara, J.; Ruud, K.; Vahtras, O.; Ågren, H. Internuclear distance dependence of the spin-orbit coupling contributions to proton NMR chemical shifts. Chem. Phys. Lett. 1998, 295, 455−461. (23) Ruud, K.; Åstrand, P.-O.; Taylor, P. R. An efficient approach for calculating vibrational wave functions and zero-point vibrational corrections to molecular properties of polyatomic molecules. J. Chem. Phys. 2000, 112, 2668−2683. (24) Åstrand, P.-O.; Ruud, K.; Sundholm, D. A modified variationperturbation approach to zero-point vibrational motion. Theor. Chem. Acc. 2000, 103, 365−373. (25) Böhm, M. C.; Schulte, J.; Ramírez, R. On the influence of nuclear fluctuations on calculated NMR shieldings of benzene and ethylene: A Feynman path integral-ab initio investigation. Int. J. Quantum Chem. 2002, 86, 280−296. (26) Vaara, J.; Lounila, J.; Ruud, K.; Helgaker, T. Rovibrational effects, temperature dependence, and isotope effects on the nuclear shielding tensors of water: A new 17O absolute shielding scale. J. Chem. Phys. 1998, 109, 8388−8397. (27) Wigglesworth, R. D.; Raynes, W. T.; Sauer, S. P. A.; Oddershede, J. Calculated nuclear shielding surfaces in the water molecule; prediction and analysis of σ(O), σ(H) and σ(D) in water isotopomers. Mol. Phys. 1999, 96, 1595−1607. (28) Anderson, J. B. A Random-Walk Simulation of the Schrödinger Equation: H+3 . J. Chem. Phys. 1975, 63, 1499−1503. (29) Anderson, J. B. Quantum Chemistry by Random Walk. H 2P, H+3 D3h 1A1′ , H2 3Σ+u , H4 1Σ+g , Be 1S. J. Chem. Phys. 1976, 65, 4121− 4127. (30) Dykstra, C. E.; Van Voorhis, T. A. Quantum Monte Carlo Vibrational Dynamics in a Property-Based Interaction Potential Scheme for Weakly Bound Clusters. J. Comput. Chem. 1997, 18, 702−711. (31) Buch, V. Treatment of Rigid Bodies by Diffusion Monte Carlo: Application to the Para-H2···H2O and Ortho-H2···H2O Clusters. J. Chem. Phys. 1992, 97, 726−729. (32) Gregory, J. K.; Clary, D. C. Structure of Water Clusters. The Contribution of Many-Body Forces, Monomer Relaxation, and Vibrational Zero-Point Energy. J. Phys. Chem. 1996, 100, 18014− 18022. (33) Thompson, K. C.; Jordan, M. J. T.; Collins, M. A. Polyatomic Molecular Potential Energy Surfaces by Interpolation in Local Internal Coordinates. J. Chem. Phys. 1998, 108, 8302−8316. (34) Bettens, R. P. A.; Collins, M. A. Learning to Interpolate Molecular Potential Energy Surfaces with Confidence: A Bayesian Approach. J. Chem. Phys. 1999, 111, 816−826. (35) Zhang, D. H.; Collins, M. A.; Lee, S.-Y. First-Principles Theory for the H + H2O, D2O Reactions. Science 2000, 290, 961−963. (36) Bettens, R. P. A.; Collins, M. A.; Jordan, M. J. T.; Zhang, D. H. Ab Initio Potential Energy surface for the Reactions between H2O and H. J. Chem. Phys. 2000, 112, 10162−10172. (37) Fuller, R. O.; Bettens, R. P. A.; Collins, M. A. Interpolated Potential-Energy Surface and Reaction Dynamics for BH+ + H2. J. Chem. Phys. 2001, 114, 10711−10716. (38) Collins, M. A.; Petrie, S.; Chalk, A. J.; Radom, L. ProtonTransport Catalysis and Proton-Abstraction Reactions: An Ab Initio Dynamical Study of X + HOC+ and XH+ + CO (X = Ne, Ar, and Kr). J. Chem. Phys. 2000, 112, 6625−6634. (39) Bettens, R. P. A.; Hansen, T. A.; Collins, M. A. Interpolated Potential Energy Surface and Reaction Dynamics for O(3P) + H+3 (1A′1) and OH+(3Σ−) + H2(1Σ+g ). J. Chem. Phys. 1999, 111, 6322−6332. Cited by ref 32. (40) Collins, M. A.; Bettens, R. P. A. Potential Energy Surface for the Reactions BeH2 + H ⇌ BeH + H2. Phys. Chem. Chem. Phys. 1999, 1, 939−945. (41) Bettens, R. P. A.; Collins, M. A. Interpolated Potential Energy Surface and Dynamics for the Reactions between N(4S) and H+3 (1A′1). J. Chem. Phys. 1998, 109, 9728−9736. 1306


Comparing Vibrationally Averaged Nuclear Shielding Constants by

Recommend Documents