Ab Initio NMR Chemical Shift Calculations Using ... - ACS Publications

Oct 18, 2016 - Australian National University, Leonard Huxley Bldg 56, Mills Rd, Canberra, ACT 2601, Australia. •S Supporting Information. ABSTRACT:...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCA

Ab Initio NMR Chemical Shift Calculations Using Fragment Molecular Orbitals and Locally Dense Basis Sets Published as part of The Journal of Physical Chemistry virtual special issue “Mark S. Gordon Festschrift”. Roger Amos* and Rika Kobayashi Australian National University, Leonard Huxley Bldg 56, Mills Rd, Canberra, ACT 2601, Australia S Supporting Information *

ABSTRACT: We revisit the NMR shielding constants of a model 10-residue peptide system by investigating the use of Jensen’s NMR specialized basis sets and locally dense basis sets within the fragmentation molecular orbital scheme. It is found that this procedure can reproduce the shielding constants of a full calculation at only a fraction of the cost. Strategies for choosing fragments and complications that need to be considered within the method are discussed extensively.



where μN is the nuclear magneton, gA is the nuclear g-factor, riA is the vector from electron i to nucleus A, and riG is the vector from electron i to the gauge origin. The paramagnetic terms involves derivatives of the wave function, or the density3,4

INTRODUCTION It is now well-established that nuclear magnetic resonance (NMR) spectroscopy has become a powerful experimental tool for obtaining information about the structure and dynamics of proteins, and nucleic acids and their complexes. Computationally, the basic theory for calculating NMR shieldings was formulated by Ramsey in the early 1950s.1 However, the beginnings of modern methods were provided by Stevens, Pitzer, and Lipscomb in 1963.2 The calculations presented by Stevens et al. were simple by modern standards but they outlined the methods which would become important. The NMR shielding tensor is a second-order response property where the perturbations are the external magnetic field, B⃗ , and magnetic moment m of nucleus of interest σ=

σp = − 2 ∑

where

LG =

Hp =

The magnetic moment and the magnetic field are both vectors, so σ is actually a tensor, but only the trace of the tensor (the isotropic shielding) is considered here. It can further be written as a sum of two terms, a diamagnetic term and a paramagnetic term

⎛ gA μ N ⎞⎛ ri A × pi ⎞ ⎜ ⎟⎜ ⎟ ⎝ c 2 ⎠⎝ ri A 3 ⎠

The main contribution of Stevens et al. was to provide a way of evaluating this derivative, which is now routine on a much larger scale than their early calculations. However, there is a problem with these formulas in that the separation into diamagnetic and paramagnetic terms is not uniquethe two terms depend upon the gauge origin, i.e., the choice for the origin of the magnetic field. This is a feature of the choice of Hamiltonian used to derive the equations. Over the years, this gauge problem was tackled by many approaches, such as individualized gauge for localized orbitals (IGLO),5 local orbital/local gauge (LORG),6 continuous set of gauge

σ = σd + σ p

The diamagnetic term depends only on the unperturbed wave function and is calculated as the expectation value of the diamagnetic shielding operator σ d = Ψ0|H d|Ψ0 ⎛ gA μ N ⎞⎛ riTGri A − ri ArITG ⎞ ⎜ ⎟⎜ ⎟ ⎝ 2c 2 ⎠⎝ ri A 3 ⎠ © 2016 American Chemical Society

1 riG × pi 2

is the angular momentum operator (which depends on the gauge origin) and

∂ 2E ∂B⃗ ∂m

Hd =

Ψ0|H p|Ψn Ψn|LG|Ψ0 ∂Ψ = −2 Ψ0|H p| 0 E 0 − En ∂LG

Received: September 11, 2016 Revised: October 17, 2016 Published: October 18, 2016 8907

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915

Article

The Journal of Physical Chemistry A

Figure 1. Scheme of the 10-residue polypeptide model system.

Figure 2. 10-peptide model system, including the side chains.

transformations (CSGT),7 and gauge including atomic orbitals (GIAO).8 GIAOs incorporate a factor into the basis set that cancels the origin dependence of the Hamiltonian, is generally reckoned to be the most rigorous approach, and is our method of choice for this paper, used as implemented within density functional theory.9 Alongside method development, specialized basis sets for NMR studies have been developed, the best known being those by Jensen.4,10,11 These originated with the pcJ-n series4 for accurate calculation of nuclear spin−spin coupling constants, by analyzing the basis set requirements of the different operators that contribute, viz. diamagnetic spin−orbit operator, paramagnetic spin−orbit operator, spin−dipole operator, Fermicontact term. Similar analysis of the NMR shielding constants found the only additional basis set requirement to be for the paramagnetic spin−orbit operator, where only p-type tight functions have any significant influence. Thus, the pcS-n hierarchy10 adds a tight p-primitive and reoptimizes the contraction pattern; pcSseg-n introduced in 201511 is just a segmented contracted version for computational efficiency. Ab initio NMR functionality is now included in all major quantum chemistry programs. However, the prohibitive cost of the most accurate methods such as CCSD(T) still means that application to large systems, particularly the aforementioned biological system, remains impractical. For this reason, fragmentation methods have rapidly gained interest as a means of reducing the cost of quantum chemical problems, as cheaper methods are often unable to reach chemical accuracy. The idea behind fragmentation methods is to break a molecule into many fragments, QM calculations are performed on the fragments, and the properties for the whole molecule are reconstituted from the corresponding values for

the fragments. Note that the (1/r3) term in the equations above means that the operators depend strongly on regions of the electron density close to the nucleus, meaning that NMR shielding is a local property and thus particularly suited to evaluation by fragment methods. A variety have already been applied to the calculation of NMR shieldings, such as • fragment molecular orbitals (FMO)12 • density-matrix based algorithm, GIAO-HF13 • AF-QM/MM14,15 • adjustable density matrix assembler16 • combined fragmentation methods17 • systematic molecular fragmentation by annihiliation (SMFA)18 We take one such approach, using some of the methodology of fragment molecular orbitals, and investigate the effects of using Jensen’s NMR specialized basis sets and locally dense basis sets.19 The local nature of the chemical shielding makes locally dense basis set models particularly suitable for use with fragment methods, as it allows different fragments to use different basis sets with the less important fragments being given less weight. This has the obvious benefit of retaining accuracy while saving on cost.



METHODOLOGY We revisit the 10-residue model peptide system of Gao et al.12 for which they carried out calculations of shielding constants using their FMO-based method with the 6-31G(d) basis set. Their method uses the GAMESS package20 for the FMO calculations, in combination with the Gaussian package21 for NMR GIAO and CSGT shift calculations using the FMO fragment densities and point charges. They found values of absolute isotropic shielding constants using their FMO-based 8908

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915

Article

The Journal of Physical Chemistry A

construction and is large enough to count as the practical basis set limit. Jensen has an even larger basis, pcSseg-4, but tests on small molecules10 shows hardly any change (∼0.01 ppm or less) from pcSseg-3. A graphical illustration of the similarity between pcSseg-2 and pcSseg-3 is given in Figure 4, showing the two sets of results to be virtually identical. Equivalent graphs for 6-311+G(2d,p) and pcSseg-1 are in the Supporting Information. The RMS values in Table 1, and the full data in the Supporting Information, show that the results for N and O nuclei are more sensitive than for H and C, particularly for oxygen, which cannot be considered truly converged even with pcSseg-3 (although that may be a feature of this particular system). The C and H atoms are less sensitive to the basis set, but even for these, smaller bases, such as 6-31G(d) or pcSseg-1, have significant deviations from the basis set limit. The RMS deviations could be reduced if presented as chemical shifts rather than for absolute shieldings. For example, if the absolute shieldings for, for example, C6H6 were calculated from 631G(d) and pcSseg-3 and then the peptide results were taken with respect to these as references, there would be some apparent cancellation of errors, and the smaller basis set results would look better. To be specific, the Carbon absolute shielding in C6H6 is 68.3 ppm with the 6-31G(d) basis and 39.8 ppm with pcSseg-3, a difference of 28.5 ppm. The RMSD between 6-31G(d) and pcSseg-3 for absolute shieldings is 21.8, so shifting the values by 28.5 ppm would bring chemical shifts closer together (actually overshooting a little) and reduces the RMSD to 11.1, an apparent increase in accuracy. However, there is a hidden problem with this: such a procedure shifts all the values by the same amount and thus cannot change either the order or the relative spacing of the lines. In other words, the spectrum is actually exactly the same, just translated to a different point on the axis. Closer examination of the full set of values (Supporting Information) reveals that many of the predicted shifts change order; i.e., simply shifting the results from a small basis would ascribe the shieldings to the wrong nuclei. This particularly affects the H shieldings; the C results are less affected. This problem would also apply to approaches22 that scale the results, as these also cannot change the order. The results show that to reach the basis set limit, one must use sets at least as large as pcSseg-3. This is in accord with the papers of Jensen4,10,11 who also considered a wide range of other basis sets. Other papers23−25 have considered the Jensen basis sets in conjunction with both DFT and correlated methods and have reached the same conclusion. Note that the basis set dependency is actually larger than the difference between DFT methods. The examples given above for calculations on the whole molecule (and it should be remembered that this is just a section of a larger protein) show that (a) much larger basis sets are needed than is generally assumed to eliminate the basis set error and (b) the cost would be prohibitive if routine application to molecules of this size were required. To see this, the total time for the calculation at the pcSseg-3 level (10 338 basis functions) was 2165 h (just over 90 days!) on a 16 core cluster. About 80% of this time was in solving the coupled response equations, which can be avoided for some approaches,26 but this would still be prohibitive on a larger molecule. For this reason we have investigated fragmentation methods, in the next section.

method to be close to a conventional, full-molecule method with RMS deviations being not much greater than 1 ppm. However, the errors due to incomplete basis sets are significantly greater than this, so the first part of the present investigation is to determine the basis set limit. We have used a similar method but employed a different fragmentation scheme and investigated the use of specialist and locally dense basis sets. The model system (MetGlnIlePheValLysThrLeuThrGly), schematically shown in Figure 1, was built by taking the sequence of the first 10 residues of ubiquitin protein, chosen as a representative example of a large biological molecule. The geometry was obtained from averaged snapshots of MD simulations, and kindly provided to us by Dr Gao; coordinates are given in Table 1 of the Supporting Information. An image of the whole molecule is shown in Figure 2, primarily to illustrate the positions and complexity of the sidechains. Calculations on the Whole Molecule. The full molecule absolute shieldings have been calculated in the gas phase for the model peptide system at the B3LYP level with the Jensen pcSseg-1, pcSseg-2, pcSseg-3 basis sets, and also with the 631G(d) and 6-311+G(2d,p) bases as these are commonly recommended as being suitable for NMR calculations. The full results are given in the Supporting Information. RMS deviations for the various basis sets compared to pcSseg-3 are given in Table 1. The two terminal oxygens are omitted from the analysis as their inclusion distorts the results misleadingly. Table 1. RMS Deviation (and Maximum Error) of Shieldings (ppm) Compared to the Basis Set Limit pcSseg-1 nuclei

RMSD

max error

all C H N O

4.09 4.9 0.7 9.0 6.4

18.7 16.9 2.9 18.7 10.3

pcSseg-2 RMSD

max error

0.4 0.6 0.1 0.9 0.75

1.5 1.0 0.4 1.3 1.5

6-31G(d)

6-311+G(2d,p)

RMSD

max error

RMSD

max error

18.1 21.8 1.1 25.3 42.2

62.6 35.8 2.0 37.7 62.6

4.9 6.6 0.6 5.1 10.6

14.4 10.0 2.0 6.7 14.4

Looking at the differences between pcSseg-2 and pcSseg-3 it is seen that the results are essentially converged with respect to basis set size and we take the pcSseg-3 results as our basis set limit. A plot of the absolute shieldings, not chemical shifts, for comparison of 6-31G(d) to pcSseg-3, is given in Figure 3. Values lying on the diagonal blue line would agree exactly with pcSseg-3, which shows 6-31G(d) overestimates in general, apart from a few values, notably the oxygens which are seriously too negative. The terminal oxygen atoms in this polypeptide are extremely dependent on the basis set. This is probably due to the molecule as specified having partial zwitterionic character, which results in the charges on these oxygens varying considerably with basis set. In fact, results for all the oxygens in the molecule are more sensitive to basis set than are the other nuclei. The benefits of the specialized Jensen basis sets can be seen by comparing pcSseg-1 with 6-31G(d), both of which are essentially DZP (double-ζ plus polarization) quality, or pcSseg2 with 6-311+G(2d,p), which are TZ2P (triple-ζ with two polarization shells). The pcSseg-3 basis is QZ3P (or better) in 8909

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915

Article

The Journal of Physical Chemistry A

Figure 3. 6-31G(d) isotropic shieldings (in ppm) compared to pcSseg-3 (=basis set limit). The red dots are the 6-31G(d) values.

Calculations Using Fragmentation Approach. The initial fragmentation started with the Facio software.27 This divides the 10-peptide molecule into ten fragments each comprising one peptide unit and its associated side chain. Figure 5 shows a screenshot where the yellow and red cubes show the fragmentation points. Preliminary experimentation revealed the following: • It is not possible to get meaningful results from calculations on fragments that contain only one peptide unit. • It is not possible to divide a peptide unit; each CONH group must be treated as a single entity, with a basis set of equal quality on each atom. • To obtain meaningful results for a specific atom, or group, the minimum requirement is that its immediate neighbors must be included. These factors mean that the minimum sized fragment will/ should contain three peptide units, possibly with their side chains still attached. Figure 6 shows a typical such three-peptide fragment from the center of the molecule. Only the shielding constants from the center peptide in this fragment can reproduce the full molecule calculation. To obtain all the

shielding constants, this means that a series of fragment calculations are neededeight different calculations, in this case, as the end units are unique. NMR shielding calculations on these fragments were carried out with the Jensen pcSseg-n series of basis sets. The results quoted below are for the pcSseg-3 basis. Figure 7 shows the plot of the absolute shieldings of the model peptide system using the fragment approach described above compared to the full molecule calculation at the pcSseg-3 level. The agreement is excellent, with the exception of a few points, primarily the oxygen atoms, which, as remarked previously, are very sensitive to the details of the calculation. The full list of results is in the Supporting Information. However, the main thing to notice is that whereas the whole molecule calculation took 2165 h, the times for the fragments are much less. The times taken for the analogous fragment calculations are shown in Table 2, and the longest took 21 h. Because each fragment calculation is independent, they can all be made simultaneously, and thus the entire result for the whole molecule can be obtained in this amount of time on a supercomputer, 2 orders of magnitude faster. Note that this is 8910

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915

Article

The Journal of Physical Chemistry A

Figure 4. Shieldings (ppm) with pcSseg-2 (red dots) compared to pcSseg-3. Values on the blue line agree exactly.

Figure 5. Whole molecule with fragmentation points as calculated with Facio.

8911

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915

Article

The Journal of Physical Chemistry A

Figure 6. Example three-peptide fragment for the 10-residue polypeptide.

calculations the charges have little effect, probably because the fragments are large. Figure 9 compares accuracies of the two different fragment approaches, one with the local dense basis and one with the largest basis on each fragment. The two fragmentation schemes give almost the same results; again the outliers are the oxygen atoms. Finally, it could be criticized that the density functional used in this study, B3LYP, may be considered as getting past its useby date and indeed there are many studies where it is applied inappropriately. However, it should be noted that subsequent functional development has been in areas, e.g., range separation and dispersion, that will not significantly affect NMR property calculations. There are more recent functionals, such as the NMR-specialized KT2 and KT3 functionals,29,30 though they are not so much a functional improvement, as a reparametrization. Earlier experience31 leads us to believe that the effect of the functional is not as crucial to the accuracy of the results as having the basis set correct, which is supported by other studies.10,23,24 We should emphasize that our findings regarding fragmentation and local basis sets are generally valid for any DFT or wave function method.

less than even the pcSseg-2 calculation for the whole molecule, which took 34 h. However, though it “only” requires 18 h of calculation time, the fragment shown in Figure 6 is still quite large; it has 57 atoms and with the pcSseg-3 basis there are 2895 basis functions. It is possible to reduce the time even further using the idea of a locally dense basis. Following our earlier study28 on systematic investigation of locally dense basis sets, where it was found that to achieve chemical accuracy, by which we mean within 0.1 ppm for H, 1 ppm for C, and 10 ppm for the other nuclei, required at least pcS-1 for all parts of the molecule coupled with pcS-3 on the group of interest, the three-peptide fragments were further divided, as schematically shown in Figure 8. The basis set on the central unit is kept as pc0Sseg-3 but reduced to pcSseg-1 on the neighboring units Table 3 shows how the calculation time for the individual fragment jobs using pcSseg-3 on the central peptide unit and pcSseg-1 on adjacent peptide units, compared to those in Table 2. The timings show that using a locally dense basis reduces the cost by almost another order of magnitude for most fragments. The end fragments take longest as these require two of the components to retain the full pcSseg-3 basis. It is also possible to include long-range effects by embedding the fragments in charges representing the “missing” part of the molecule. The charges used were from an FMO calculation, as this procedure iterates the monomer fragments until they are self-consistent with each other, and thus the derived charges represent those in the full molecule. This is preferable to using charges obtained from isolated fragments. In the present



CONCLUSION We have investigated the use of specialized basis sets for the calculation of NMR shielding constants within the FMO fragmentation scheme on a model peptide system acting as a prototype for the ubiquitin molecule. The factors that determine the accuracy of the calculated NMR shielding constants include the basis set, the methodology, e.g., DFT, MP2, etc., vibrational corrections, solvent effects, and in a large molecule, the probability that the geometry can only be defined 8912

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915

Article

The Journal of Physical Chemistry A

Figure 7. Shieldings (ppm) for pcSseg-3 fragments versus pcSseg-3 for the whole molecule.

Table 2. Calculation Time for Individual Fragments at the B3LYP/pcSseg-3 Level To Be Compared with 2165 h for the Whole Molecule Calculation

Table 3. Reduction in Calculation Time for the Individual Fragment Jobs in Table 2 Using pcSseg-3 on the Central Peptide Unit and pcSseg-1 on Adjacent Peptide Units

fragment

no. of atoms

no. of basis functions

calculation time (h)

fragment

no. of basis functions

calculation time (h)

1−2−3 2−3−4 3−4−5 4−5−6 5−6−7 6−7−8 7−8−9 8−9−10

54 58 57 60 54 57 49 44

2683 3006 2895 3033 2679 2817 2475 2310

16.1 19.2 18.1 13.5 15.5 16.5 13.1 21.5

1−2−3 2−3−4 3−4−5 4−5−6 5−6−7 6−7−8 7−8−9 8−9−10

1907 1374 2213 1261 1405 1154 1257 1657

6.0 2.2 3.3 1.4 2.1 1.1 1.7 3.7

Figure 8. Scheme of the fragment-based NMR calculation model for a 10-residue polypeptide using locally dense basis sets. 8913

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915

Article

The Journal of Physical Chemistry A

Figure 9. Comparison of shielding results (ppm) from locally dense fragments against fragments with the full basis.

as a thermal average. Of these factors only the basis set dependence is investigated here. It has been found that to reach chemical accuracy requires basis sets of quadruple-ζ quality, which are very expensive. However, the cost of the calculation can be dramatically reduced (by 2 orders of magnitude or more), without significant loss of accuracy, by use of fragmentation methods, in combination with judicious placement of locally dense basis sets. Strategies regarding effective fragmentation, such as cleavage points and fragment size, are being clarified, and with this refined methodology we are in a position to tackle real systems of real interest.





6-311+G(2d,p) versus pcSseg-3; Figure S2, shieldings (in ppm) for pcSseg-1 versus pcSseg-3; Figure S3, shieldings (in ppm) from fragment approach with pcSseg-3 fragments compared to fragments with pcSseg-3/ pcSseg-1 locally dense basis (PDF)

AUTHOR INFORMATION

Corresponding Author

*R. Amos. E-mail: [email protected]. Phone: +612 6125 3436. Notes

ASSOCIATED CONTENT

The authors declare no competing financial interest.

S Supporting Information *



The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b09190. Table S1, geometry of the 10-peptide model system (in Ångstrom); Table S2, full molecule absolute shieldings (not chemical shifts), in ppm, for 6-31G(d),6-311+G(2d,p), pcSseg-1, pcSseg-2, pcSseg-3 and for the fragment approach; Figure S1, shieldings (in ppm) for

ACKNOWLEDGMENTS

This research was undertaken on the NCI National Facility in Canberra, Australia, which is supported by the Australian Commonwealth Government. The authors thank Dr. Qi Gao for kindly providing the geometry for the model system. 8914

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915

Article

The Journal of Physical Chemistry A



(23) Flaig, D. M.; Maurer, M.; Hanni, M.; Braunger, K.; Kick, L.; Thubauville, M.; Ochsenfeld, C. Benchmarking hydrogen and carbon NMR chemical shifts at HF, DFT, and MP2 levels. J. Chem. Theory Comput. 2014, 10, 572−578. (24) Magyarfalvi, G.; Pulay, P. Assessment of density functional methods for nuclear magnetic resonance shielding calculations. J. Chem. Phys. 2003, 119, 1350−1357. (25) Jensen, S. R.; Flå, T.; Jonsson, D.; Monstad, R. S.; Ruud, K.; Frediani, L. Phys. Chem. Chem. Phys. 2016, 18, 21145−21161. (26) Helgaker, T.; Wilson, P. J.; Amos, R. D.; Handy, N. C. Nuclear shielding constants by density functional theory with gauge including atomic orbitals. J. Chem. Phys. 2000, 113, 2983−2989. (27) Suenaga, M. Facio: New Computational Chemistry Environment for PC GAMESS. J. Comput. Chem., Jpn. 2005, 4, 25−32. (28) Reid, D. M.; Kobayashi, R.; Collins, M. A. Systematic study of locally dense basis sets for NMR shielding constants. J. Chem. Theory Comput. 2014, 10, 146−152. (29) Keal, T. W.; Tozer, D. J. The exchange-correlation potential in Kohn−Sham nuclear magnetic resonance shielding calculations. J. Chem. Phys. 2003, 119, 3015−3024. (30) Keal, T. W.; Tozer, D. J. A semiempirical generalized gradient approximation exchange-correlation functional. J. Chem. Phys. 2004, 121, 5654−5660. (31) Reid, D. M.; Collins, M. A. Approximating CCSD (T) Nuclear Magnetic Shielding Calculations Using Composite Methods. J. Chem. Theory Comput. 2015, 11, 5177−5181.

REFERENCES

(1) Ramsey, N. F. Electron Coupled Interactions between Nuclear Spins in Molecules. Phys. Rev. 1953, 91, 303−307. (2) Stevens, R. M.; Pitzer, R. M.; Lipscomb, W. N. Perturbed HartreeFock Calculations. I. Magnetic Susceptibility and Shielding in the LiH Molecule. J. Chem. Phys. 1963, 38, 550−560. (3) Helgaker, T.; Jaszunski, M.; Ruud, K. Ab Initio Methods for the Calculation of NMR Shielding and Indirect Spin-Spin Coupling Constants. Chem. Rev. 1999, 99, 293−352. (4) Jensen, F. The Basis Set Convergence of Spin-Spin Coupling Constants Calculated by Density Functional Methods. J. Chem. Theory Comput. 2006, 2, 1360−1369. (5) Kutzelnigg, W. Theory of magnetic susceptibilities and NMR chemical shifts in terms of localized quantities. Isr. J. Chem. 1980, 19, 193−200. (6) Hansen, A. E.; Bouman, T. D. Localized orbital/local origin method for calculation and analysis of NMR shieldings. Applications to 1 3C shielding tensors. J. Chem. Phys. 1985, 82, 5035−5047. (7) Keith, T. A.; Bader, R. F. W. Calculation of magnetic response properties using a continuous set of gauge transformations. Chem. Phys. Lett. 1993, 210, 223−231. (8) Ditchfield, R. Molecular orbital theory of magnetic shielding and magnetic susceptibility. J. Chem. Phys. 1972, 56, 5688−5691. (9) 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. (10) Jensen, F. Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods. J. Chem. Theory Comput. 2008, 4, 719−727. (11) Jensen, F. Segmented Contracted Basis Sets Optimized for Nuclear Magnetic Shielding. J. Chem. Theory Comput. 2015, 11, 132− 138. (12) Gao, Q.; Yokojima, S.; Kohno, T.; Ishida, T.; Fedorov, D. G.; Kitaura, K.; Fujihira, M.; Nakamura, S. Ab initio NMR chemical shift calculations on proteins using fragment molecular orbitals with electrostatic environment. Chem. Phys. Lett. 2007, 445, 331−339. (13) Flaig, D.; Beer, M.; Ochsenfeld, C. Convergence of electronic structure with the size of the QM region: example of QM/MM NMR shieldings. J. Chem. Theory Comput. 2012, 8, 2260−2271. (14) He, X.; Wang, B.; Merz, K. M., Jr. Protein NMR chemical shift calculations based on the automated fragmentation QM/MM approach. J. Phys. Chem. B 2009, 113, 10380−10388. (15) Zhu, T.; He, X.; Zhang, J. Z. H. Fragment density functional theory calculation of NMR chemical shifts for proteins with implicit solvation. Phys. Chem. Chem. Phys. 2012, 14, 7837−7845. (16) Exner, T. E.; Frank, A.; Onila, I.; Möller, H. M. Toward the quantum chemical calculation of NMR chemical shifts of proteins. 3. Conformational sampling and explicit solvents model. J. Chem. Theory Comput. 2012, 8, 4818−4827. (17) Tan, H.-J.; Bettens, R. P. A. Ab initio NMR chemical-shift calculations based on the combined fragmentation method. Phys. Chem. Chem. Phys. 2013, 15, 7541−7547. (18) Reid, D. M.; Collins, M. A. Calculating nuclear magnetic resonance shieldings using systematic molecular fragmentation by annihilation. Phys. Chem. Chem. Phys. 2015, 17, 5314−5320. (19) Chesnut, D. B.; Moore, K. D. Locally dense basis sets for chemical shift calculations. J. Comput. Chem. 1989, 10, 648−659. (20) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347−1363. (21) Gaussian 09, Revision E.01, 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, Inc.: Wallingford, CT, 2009. (22) Kupka, T.; Stachów, M.; Chełmecka, E.; Pasterny, K.; Stobińska, M.; Stobiński, L.; Kaminský, J. Efficient modeling of NMR parameters in carbon nanosystems. J. Chem. Theory Comput. 2013, 9, 4275−4286. 8915

DOI: 10.1021/acs.jpca.6b09190 J. Phys. Chem. A 2016, 120, 8907−8915