# High Level Electronic Structure Calculation of Molecular Solid-State NMR Shielding Constants

Corentin Poidevin, Georgi L Stoychev, Christoph Riplinger, Alexander A Auer

► To cite this version:

Corentin Poidevin, Georgi L Stoychev, Christoph Riplinger, Alexander A Auer. High Level Electronic Structure Calculation of Molecular Solid-State NMR Shielding Constants. Journal of Chemical Theory and Computation, 2022, 18 (4), pp.2408-2417. 10.1021/acs.jctc.1c01095. hal-03632120

HAL Id: hal-03632120

https://hal.science/hal-03632120

Submitted on 17 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# High Level Electronic Structure Calculation of Molecular Solid-State NMR Shielding Constants

Corentin Poidevin, Georgi L. Stoychev, Christoph Riplinger, and Alexander A. Auer*

<figure><img src="image_2.png" /></figure>

Cite This: J. Chem. Theory Comput. 2022, 18, 2408–2417

<figure><img src="image_3.png" /></figure>

Read Online

## ACCESS

<figure><img src="image_4.png" /></figure>

Metrics & More

<figure><img src="image_5.png" /></figure>

Article Recommendations

<figure><img src="image_6.png" /></figure>

Supporting Information

**ABSTRACT:** In this work, we present a quantum mechanics/molecular mechanics (QM/MM) approach for the computation of solid-state nuclear magnetic resonance (SS-NMR) shielding constants (SCs) for molecular crystals. Besides applying standard-DFT functionals like GGAs (PBE), meta-GGAs (TPSS), and hybrids (B3LYP), we apply a double-hybrid (DSD-PBEP86) functional as well as MP2, using the domain-based local pair natural orbital (DLPNO) formalism, to calculate the NMR SCs of six amino acid crystals. All the electronic structure methods used exhibit good correlation of the NMR shieldings with respect to experimental chemical shifts for both ¹H and ¹³C. We also find that local electronic structure is much more important than the long-range electrostatic effects for these systems, implying that cluster

approaches using all-electron/Gaussian basis set methods might offer great potential for predictive computations of solid-state NMR parameters for organic solids.

<figure><img src="image_7.png" /></figure>

## 1. INTRODUCTION

Besides X-ray diffraction, solid-state nuclear magnetic resonance (SS-NMR) spectroscopy has become a powerful way to study molecular crystal structures. SS-NMR can be used for powder or amorphous samples and yields information about the local environment of NMR active nuclei. However, the structure cannot be resolved solely from experimental NMR data. Thus, SS-NMR results are often combined with plane wave density functional theory (DFT) electronic structure calculations with gauge-including projector augmented plane waves (GIPAW) to validate or refine structural information.¹⁻³ While periodic DFT has the advantage of simulating the crystal environment, calculations of NMR properties are in the vast majority of the cases limited to the use of (meta-)GGA functionals which have been shown to produce NMR shielding constants (SCs) of limited accuracy.⁴⁻⁷ Typically, the SCs are highly sensitive to the description of the electronic structure near the nucleus. As a consequence, in molecular calculations that aim at high accuracy, it is common to use special basis sets that include Gaussians with high exponents, and post-Hartree–Fock (HF) calculations require the inclusion of core correlation. In contrast to this, periodic boundary calculations commonly apply effective core potentials even for light atoms, and the GIPAW approach for NMR calculations combines the reconstruction of the core region using the projector augmented wave technique in combination with gauge-including orbitals to approximate a core electron density that yields satisfactory accuracy in DFT applications.

In parallel to the plane wave approach, quantum mechanics/molecular mechanics (QM/MM) models of molecular crystals have also been developed for the calculation of NMR SCs. This approach mainly consists of focusing on a central QM region, composed of an asymmetric unit surrounded by a few other asymmetric units, which is embedded in MM point charges. Here the aim is to concentrate the computational effort on obtaining the best possible local electronic environment for the central asymmetric unit by using higher level methods for the electronic structure as well as atom-centered basis sets. However, this implies that the QM region has to be of limited size, as higher-level electronic structure methods can become extremely expensive computationally. To circumvent this issue, over the last 10 years Beran et al. have developed a fragment-based method allowing large QM regions to be tackled at affordable computational costs.⁸⁻¹¹ On the other hand, cluster approaches with or without MM point charges have also been used.¹¹⁻¹⁴ In the cluster approach, the whole QM region is included in the NMR calculation while only considering the NMR SCs of the central molecule.

Received: October 31, 2021

Published: March 30, 2022

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 1.** The six studied amino acids.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 2.** QM1 corresponds to the molecule (asymmetric unit) on which the NMR properties are calculated. QM2 corresponds to the first shell of the asymmetric units around QM1 (i.e., with at least one atom within 2.5 to 3.5 Å depending on the system). The MM region corresponds to the duplication of 9 × 9 × 9 unit cells around QM1 and is composed of self-consistently optimized point charges equal to the CHELPG charges calculated for the QM1 atoms.</figcaption></figure>

As the size of the QM system in the QM/MM approach can become rather big if one wants to include all the adjacent molecules to the central asymmetric unit, the literature only reports results using (meta-)GGA or hybrid functionals for the calculation of SCs for SS-NMR. However, double-hybrid functionals such as DSD-PBEP86 or B2PLYP and post-HF methods like MP2 have been shown to lead to a great improvement of accuracy for the calculation of NMR SCs on small molecules, when compared with precise but computationally demanding methods like coupled cluster theory.<sup>15,16,5,17</sup> Recently, Dračinský et al. presented a study on the molecular crystals of six amino acids using MP2 or CCSD corrections on PBE GIPAW results.<sup>18</sup> These corrections were obtained by calculating the isolated single amino acid molecules at a higher level of theory. However, the main improvement of their results was found when they included vibrational effects via the use of path-integral molecular dynamics (PIMD).

A recent study in our group focused on using a QM/MM approach to compute the NMR SCs for ionic solids.<sup>15</sup> In this work it was found that the size of the QM region is more important than the way the embedding was done. Furthermore, MP2 and DH functionals were found to lead to a significant improvement over standard DFT for the calculation of SCs, giving more consistent results for all tested nuclei, in agreement with what has been discussed in the framework of fragmentation or composite methods to compute NMR parameters.<sup>19,20</sup> However, these improvements come at a price, as methods that include post-SCF correlation are both inherently more expensive and require larger basis sets to converge toward the complete basis set (CBS) limit, compared to other DFT approaches. Kaupp and co-workers, for example, have shown that, for NMR shielding calculations, local hybrid functionals fall between conventional functionals and DHs in terms of both accuracy and computational cost.<sup>21,22</sup>

Recently, local correlation methods for the calculation of NMR SCs at MP2 and DH functional levels have been developed and implemented.<sup>23-26</sup> This allows for the use of high-level electronic structure methods with large QM/MM clusters without drastic increase of computational effort. Thus, as a proof of concept and first benchmark, we have decided to study the same amino acid molecular crystals that Dračinský et al. have studied and use their experimental results.<sup>18</sup> However, we use nonperiodic methods with embedding as this allows us to apply methods developed for molecular systems. Namely, we use the domain-based local pair natural orbital (DLPNO) approximation as this has been shown to provide high accuracy and allows sufficiently large systems to be treated in the QM region.<sup>26-28</sup> Ultimately, we wish to investigate whether this approach yields better accuracy than GIPAW in a black-box way for practical applications. Our main focus here is on the ¹H NMR shieldings, as these are the most common experimentally but computationally challenging in SS-NMR.

## 2. STRUCTURES AND METHOD

### 2.1. Studied Systems.
As in the work of Dračinský et al., the structures of the 6 amino acids (namely, α-glycine, L-alanine, L-serine, L-aspartic acid, L-cysteine, and L-threonine) were obtained from the Cambridge Structural Database (CSD refcodes: GLYCIN29, LTHREO01, LASPRT, LCYSTN21, LSERIN01, and LALNIN12, Figure 1).<sup>29</sup> These amino acids display a large variety of H environments, i.e., C–H, S–H, N–H, and O–H, which in the latter two cases form hydrogen bonds with neighboring oxygen from carboxylates.

### 2.2. Embedded Cluster Approach and Computational Details.
Our cluster model for the calculation of NMR chemical shifts consists of (1) an asymmetric unit (QM1) of the central unit cell for which the NMR chemical shifts are calculated, (2) the first shell of asymmetric units (QM2) around QM1 (i.e., molecules with at least one atom within 2.5 Å).
<table><caption>Table 1. Correlation Parameters between Experimental ¹H Shifts and the Calculated NMR Shieldings Calculated with B3LYP and Different Levels of Embedding<sup>a,b</sup></caption><thead><tr><th></th><th>slope</th><th>intercept</th><th>R<sup>2</sup></th><th>MAE</th><th>MaxAE<sup>c</sup></th><th>SDE</th></tr></thead><tbody><tr><td>B3LYP</td><td>-1.03</td><td>31.14</td><td>0.9915</td><td>0.21</td><td>0.68</td><td>0.30</td></tr><tr><td>B3LYP (QM1/QM2)</td><td>-1.03</td><td>31.09</td><td>0.9869</td><td>0.26</td><td>0.89</td><td>0.37</td></tr><tr><td>B3LYP (QM1/MM)</td><td>-0.57</td><td>29.91</td><td>0.9313</td><td>1.06</td><td>4.85</td><td>1.46</td></tr><tr><td>B3LYP (QM1)</td><td>-0.38</td><td>29.44</td><td>0.5806</td><td>1.81</td><td>7.90</td><td>2.72</td></tr></tbody></table>

<sup>a</sup>Intercept, MAE, MaxAE, and SDE are given in ppm, while the slope and R<sup>2</sup> are unitless. <sup>b</sup>See text for embedding scheme details. <sup>c</sup>Maximum error and the corresponding hydrogen atom.

to 3.5 Å depending on the system), and (3) a point charge field corresponding to the duplication of 9 × 9 × 9 unit cells around the central cell (MM) as depicted in Figure 2. In order to include the first shell of molecules around QM1, a total of 13 molecules were used for the QM systems (QM1+QM2) in all cases with the exception of α-glycine, where 14 had to be used.

All calculations were carried out using a development version of the ORCA 5.0 program package.<sup>30,31</sup> As mentioned above, the NMR properties are only reported for the QM1 molecules in the presence of the QM2 molecules and the point charges using gauge including atomic orbitals (GIAOs). The electronic structures of QM1 and QM2 are calculated using the same method. In this study we used four classes of functionals, GGA (PBE), meta-GGA (TPSS), hybrid (B3LYP), and DH (DSD-PBEP86) as well as MP2.<sup>32–36</sup> In all cases, the RI (PBE and TPSS) or RIJCOSX approximations are used, and the atom-pairwise dispersion correction with the Becke-Johnson damping scheme (D3BJ) was used for the DFT calculations.<sup>37–39</sup> In order to reduce the computational cost, a smaller basis set is often used for the QM2 molecules (def2-TZVP or pcSseg-2 instead of pcSseg-3).<sup>40,41</sup> In all cases the DefGrid3 grid settings were used. Note that, in this work, we have applied the ad hoc gauge-invariant approach for the kinetic energy density (τ) terms in meta-GGA functionals, which is the default in ORCA 5. Recent work has demonstrated that a more rigorous treatment is provided by the Dobson ansatz for τ,<sup>42–44</sup> although the differences for some functionals, such as TPSS, are small.<sup>21</sup> We have included data calculated with TPSS and the Dobson ansatz in the Supporting Information (Table S9), and indeed, the results are very similar.

The charges of the MM region are optimized by converging them to those of the atomic (CHELPG) charges of the central asymmetric unit (QM1) through successive SCF procedures down to a threshold of 0.01.<sup>12,45</sup> For this purpose the atomic charges of the QM atoms are first computed without the surrounding point charges. These atomic charges are then mapped onto the equivalent MM atoms (the super cell consists of repeating molecules). The atomic charges are then iteratively computed in the field of the surrounding MM charges, until convergence is achieved. The property calculation is then carried out using the converged atomic charges for the MM atoms. Note that we have not included relativistic effects in the computation of the NMR chemical shifts, as we mostly focus on light atoms. However, in the presence of heavier nuclei, this might lead to deviations—in the case of sulfur, up to 0.5 ppm for ¹H SCs.<sup>21,46,47</sup>

The experimental crystal structures are used as the initial structure for the construction of our models. Then, we optimize the positions of the hydrogen atoms of the QM1 molecule at the PBE-D3BJ/def2-SVP level in the presence of the fixed QM2+MM embedding. Once optimized, the whole

QM1+QM2+MM is rebuilt using the new positions of the hydrogen atoms. This structure is then used for the calculation of the SCs.

Due to the cost of calculating such large systems at the DH or MP2 level, the DLPNO approximation is used in both cases, combined with a multilevel approach, as described in ref 48. MP2 correlation contributions from electron pairs involving at least one occupied orbital localized to the QM1 fragment were treated at the “NormalPNO” level and the rest at the “LoosePNO” level.<sup>26–28</sup> This approximation saves computational time and is expected to have a negligible effect on the calculated shieldings.

Molecular dynamic (MD) simulations were also done to estimate vibrational effects (see Section 3.5). These MD simulations were done at the GFN2-xTB level at 300 K with a time step of 1 fs.<sup>49,50</sup> Non-hydrogen atoms were held fixed. The total run length was 40 ps. After 1 ps of equilibration, 100 snapshots (evenly separated by 0.3 ps) were used for single point calculations. The successive snapshot calculations (see Section 3.5) were performed directly on the GFN2-xTB geometries, as these have been reported to yield reasonable structures also for hydrogen bonded systems.<sup>50</sup>

The correlation between experimental shifts and computed shieldings was analyzed using a linear regression fit ($\sigma_{calc} = a\delta_{exp} + b$), for which the slope (a), intercept (b), and coefficient of determination R<sup>2</sup> are reported. The predicted shifts ($\delta_{calc} = (\sigma_{calc} - b)/a$) from the fit were used to calculate the mean absolute error (MAE), maximum absolute error (MaxAE), and standard deviation of errors (SDE), see Section 1 in the SI for details. Chemical shifts were also calculated with respect to a computed reference value for the shielding in tetramethylsilane (TMS). These were also used to calculate MAE, MaxAE, SDE, and mean signed error (MSE). The structure of TMS was optimized at the B3LYP-D3BJ/def2-TZVP level using a conductor-like polarizable continuum model (CPCM) for water.<sup>51,52</sup> The experimental reference is usually an aqueous solution of sodium 3-(trimethylsilyl)-propanesulfonate (DSS); however, the ¹H chemical shift difference between TMS and DSS is less than 0.02 ppm and thus largely immaterial for the discussion here.<sup>53</sup> Note that while in SSNMR it is common to investigate linear regressions when comparing theory and experiment, choosing a suitable reference or schemes like the MSTD approach<sup>54</sup> might be advantageous for practical applications as it provides the desired cancellation of systematic errors and is not dependent on knowledge of a set of experimental values.

### 3. RESULTS AND DISCUSSION

In the following sections, we go through the various influences of the approximations made in embedded cluster electronic structure calculations on ¹H and ¹³C chemical shifts. While we will discuss the influence of basis set and electronic structure
<table><caption>Table 2. Correlation Parameters between Experimental ¹H Shifts and the Calculated NMR Shieldings Corresponding to the Values from Table 4<sup>a</sup></caption><thead><tr><th></th><th>slope</th><th>intercept</th><th>R<sup>2</sup></th><th>MAE</th><th>MaxAE<sup>b</sup></th><th>SDE</th></tr></thead><tbody><tr><td>PBE</td><td>-1.00</td><td>30.73</td><td>0.9911</td><td>0.22</td><td>0.72</td><td>SH(cys)</td><td>0.30</td></tr><tr><td>TPSS</td><td>-1.00</td><td>31.13</td><td>0.9917</td><td>0.22</td><td>0.65</td><td>SH(cys)</td><td>0.29</td></tr><tr><td>B3LYP</td><td>-1.03</td><td>31.14</td><td>0.9915</td><td>0.21</td><td>0.68</td><td>SH(cys)</td><td>0.30</td></tr><tr><td>DLPNO-DSD-PBEP86</td><td>-1.04</td><td>31.13</td><td>0.9907</td><td>0.22</td><td>0.79</td><td>SH(cys)</td><td>0.31</td></tr><tr><td>DLPNO-MP2</td><td>-1.04</td><td>31.00</td><td>0.9909</td><td>0.22</td><td>0.78</td><td>SH(cys)</td><td>0.31</td></tr><tr><td>PBE(GIPAW)<sup>c</sup></td><td>-1.10</td><td>30.75</td><td>0.9864</td><td>0.22</td><td>1.48</td><td>SH(cys)</td><td>0.38</td></tr></tbody></table>

<sup>a</sup>Intercept, MAE, MaxAE, and SDE are given in ppm, while the slope and R<sup>2</sup> are unitless. <sup>b</sup>Maximum error and the corresponding hydrogen atom.

<sup>c</sup>GIPAW values from the work of Dračinský et al.<sup>18</sup>

<table><caption>Table 3. Correlation Parameters between Experimental ¹H Shifts and the Calculated NMR Shieldings Corresponding to the Values Calculated at the DLPNO-MP2 Level Using pcSseg-3 (See also Table 4), pcSseg-2 (See also Table S2), and pcSseg-2 with NormalPNO Together with either RIJONX or RIJCOSX Approximations<sup>a</sup></caption><thead><tr><th>DLPNO-MP2 settings</th><th>slope</th><th>intercept</th><th>R<sup>2</sup></th><th>MAE</th><th>MaxAE<sup>b</sup></th><th>SDE</th></tr></thead><tbody><tr><td>Multilevel, pcSseg-3, RIJCOSX</td><td>-1.04</td><td>31.00</td><td>0.9909</td><td>0.22</td><td>0.78</td><td>SH(cys)</td><td>0.31</td></tr><tr><td>Multilevel, pcSseg-2, RIJCOSX</td><td>-1.05</td><td>30.89</td><td>0.9875</td><td>0.26</td><td>0.89</td><td>SH(cys)</td><td>0.36</td></tr><tr><td>NormalPNO, pcSseg-2, RIJONX</td><td>-1.04</td><td>31.09</td><td>0.9910</td><td>0.22</td><td>0.73</td><td>SH(cys)</td><td>0.31</td></tr><tr><td>NormalPNO, pcSseg-2, RIJCOSX</td><td>-1.04</td><td>31.10</td><td>0.9910</td><td>0.22</td><td>0.73</td><td>SH(cys)</td><td>0.31</td></tr></tbody></table>

<sup>a</sup>Intercept, MAE, MaxAE, and SDE are given in ppm, while the slope and R<sup>2</sup> are unitless. <sup>b</sup>Maximum error and the corresponding hydrogen atom.

<table><caption>Table 4. Experimental ¹H Chemical Shifts and the Corresponding Calculated NMR Shieldings in ppm</caption><thead><tr><th></th><th>hydrogen</th><th>δ<sub>exp</sub><sup>a</sup></th><th>PBE<sup>b</sup></th><th>TPSS<sup>b</sup></th><th>B3LYP<sup>b</sup></th><th>DLPNO-DSD-PBEP86<sup>c</sup></th><th>DLPNO-MP2<sup>c</sup></th></tr></thead><tbody><tr><td rowspan="3">L-alanine</td><td>H-α</td><td>3.82</td><td>26.91</td><td>27.28</td><td>27.23</td><td>27.23</td><td>27.09</td></tr><tr><td>NH<sub>3</sub></td><td>8.5</td><td>22.15</td><td>22.57</td><td>22.25</td><td>22.18</td><td>22.04</td></tr><tr><td>H-β</td><td>1.38</td><td>29.43</td><td>29.81</td><td>29.67</td><td>29.63</td><td>29.52</td></tr><tr><td rowspan="3">α-glycine</td><td>NH<sub>3</sub></td><td>8.48</td><td>22.35</td><td>22.81</td><td>22.46</td><td>22.39</td><td>22.25</td></tr><tr><td>H-α1</td><td>4.23</td><td>26.53</td><td>26.91</td><td>26.78</td><td>26.78</td><td>26.67</td></tr><tr><td>H-α2</td><td>3.06</td><td>27.51</td><td>27.90</td><td>27.79</td><td>27.79</td><td>27.69</td></tr><tr><td rowspan="6">L-serine</td><td>H-α</td><td>3.64</td><td>27.05</td><td>27.41</td><td>27.42</td><td>27.43</td><td>27.32</td></tr><tr><td>H-β1</td><td>3.75</td><td>27.03</td><td>27.45</td><td>27.38</td><td>27.34</td><td>27.21</td></tr><tr><td>H-β2</td><td>4.46</td><td>26.28</td><td>26.75</td><td>26.67</td><td>26.66</td><td>26.51</td></tr><tr><td>NH<sub>3</sub></td><td>8.37</td><td>22.37</td><td>22.84</td><td>22.50</td><td>22.43</td><td>22.30</td></tr><tr><td>OH</td><td>3.79</td><td>27.49</td><td>27.89</td><td>27.55</td><td>27.34</td><td>27.10</td></tr><tr><td>COOH</td><td>15.57</td><td>15.00</td><td>15.45</td><td>14.96</td><td>14.81</td><td>14.58</td></tr><tr><td rowspan="6">L-aspartic acid</td><td>H-β1</td><td>3.27</td><td>26.86</td><td>27.33</td><td>27.21</td><td>27.21</td><td>27.09</td></tr><tr><td>H-β2</td><td>2.54</td><td>27.73</td><td>28.13</td><td>28.00</td><td>28.01</td><td>27.81</td></tr><tr><td>NH<sub>3</sub></td><td>8.32</td><td>22.70</td><td>23.11</td><td>22.82</td><td>22.80</td><td>22.63</td></tr><tr><td>H-α</td><td>3.76</td><td>27.56</td><td>27.95</td><td>27.90</td><td>27.87</td><td>27.69</td></tr><tr><td>H-β1</td><td>3.55</td><td>27.44</td><td>27.90</td><td>27.80</td><td>27.82</td><td>27.69</td></tr><tr><td>H-β2</td><td>2.78</td><td>28.11</td><td>28.50</td><td>28.50</td><td>28.43</td><td>28.23</td></tr><tr><td rowspan="3">L-cysteine</td><td>H-α</td><td>4.28</td><td>26.74</td><td>27.13</td><td>27.10</td><td>27.14</td><td>27.00</td></tr><tr><td>NH<sub>3</sub></td><td>8.65</td><td>22.24</td><td>22.70</td><td>22.36</td><td>22.34</td><td>22.20</td></tr><tr><td>SH</td><td>1.92</td><td>28.09</td><td>28.57</td><td>28.46</td><td>28.32</td><td>28.19</td></tr><tr><td rowspan="6">L-threonine</td><td>NH<sub>3</sub></td><td>8.03</td><td>22.62</td><td>23.07</td><td>22.77</td><td>22.72</td><td>22.53</td></tr><tr><td>H-α</td><td>4.02</td><td>26.79</td><td>27.17</td><td>27.11</td><td>27.11</td><td>26.93</td></tr><tr><td>H-β</td><td>3.78</td><td>26.86</td><td>27.31</td><td>27.29</td><td>27.26</td><td>27.08</td></tr><tr><td>OH</td><td>7.95</td><td>22.39</td><td>22.84</td><td>22.44</td><td>22.31</td><td>22.17</td></tr><tr><td>H-γ</td><td>1.39</td><td>29.53</td><td>29.92</td><td>29.78</td><td>29.69</td><td>29.59</td></tr></tbody></table>

<sup>a</sup>Experimental values from ref 18. <sup>b</sup>Values calculated using the scheme described in Section 2.2 with pcSseg-3 and pcSseg-2 basis sets for QM1 and QM2, respectively. <sup>c</sup>pcSseg-3 and def2-TZVP basis sets for QM1 and QM2, respectively, were used, and NormalPNO and LoosePNO settings using the fragment scheme presented in Section 2.2 were used.

methods, the main focus lies on the question whether embedded cluster calculations allow the properties of molecular methods to be carried over to the description of NMR parameters for molecules in the solid. This includes also additional schemes like molecular higher order correlation corrections or approximate treatments of vibrational effects.

**3.1. ¹H NMR Chemical Shifts—Cluster Embedding.** To study the influence of the long-range electrostatic environment on the calculated NMR shielding, the effect of the point charges (MM part) has to be assessed. First, we find that the converged MM charges do not fluctuate significantly with the level of calculation (see Table S1 in the Supporting
Information). Second, we also calculated the NMR SCs of the six systems at the B3LYP level without including the point charges, as can be seen in the second line of Table 1 ("B3LYP (QM1/QM2)"). The obtained slope is the same as with the MM embedding. On the other hand, the R², the MAE, and the MaxAE are found to be slightly worse. As an illustration we also calculated the NMR SCs of the six molecules without QM2 and MM (see the line "B3LYP (QM1)" of Table 1) for which the correlation with the experimental chemical shifts is, as expected, very poor (R² = 0.58). Finally, we performed calculations in which the QM2 region was instead treated using point charges (labeled "B3LYP (QM1/MM)"), which is a commonly used treatment in the literature. Interestingly, the R² coefficient for this fit is much better at 0.93, but the slope and intercept of the fit differ significantly from the ideal -1 and the TMS reference shielding, respectively. Thus, at the QM1/MM level of embedding, reasonably quantitative results can be obtained, but only if the systematic errors are reduced, e.g., using a linear regression model. This shows that even though the MM embedding is an important part of the model, the most important aspect is the local chemical environment and electronic structure.

### 3.2. ¹H NMR Chemical Shifts–Basis Set, COSX, and DLPNO Errors.
The basis set is less likely to be a significant source of error here, as can be seen when comparing results from Table 2 and Table S3 in the Supporting Information: when using a smaller basis set (i.e., psSseg-2 for QM1 and def2-SVP for QM2) only slightly different results are obtained. Additionally, for MP2, we also evaluated the error induced by using the RIJCOSX approximation compared to RIJONX. As can be seen from Table 3, the RIJCOSX approximation does not lead to a significant increase in the error of the NMR SCs. The multilevel scheme for the DLPNO approximations does deteriorate the results slightly, compared to the full Normal-PNO level, but not to an extent that would change any of our conclusions.

### 3.3. ¹H NMR Chemical Shifts–Level of Theory, DFT vs DH-DFT, and GIAO vs GIPAW.
The calculated ¹H NMR shieldings using the scheme described in Section 2.2 are reported in Table 4 along with the experimental ¹H chemical shifts from the work of Dračinský et al.¹⁸ The linear fit parameters and statistical analysis are presented in Table 2, and an illustration of the fitted data is shown in Figure 3 (see Figure S1 in the SI for the rest of the fits). Surprisingly, all the electronic structure methods used (GGA (PBE), Meta-GGA (TPSS), hybrid (B3LYP), double-hybrid (DSD-PBEP86), and MP2) exhibit good correlation with the experiment in terms of slope, close to the ideal -1.0, R² coefficients of determination around 0.99, MAEs of 0.22 ppm, and maximum errors of better than 0.8 ppm. Note that Dračinský et al. reports slopes around -1.1, MAEs of 0.22–0.26 ppm, and maximum errors of 1.49–1.76 ppm for PBE(GIPAW) values and values using corrections to the electronic structure from molecular calculations.

As in the study of Dračinský et al., we observe the largest deviations for S–H of L-cysteine, but they are still smaller than when using PBE GIPAW. Two positions for the S–H are found in the crystal structure of L-cysteine, either S–H···S or S–H···O.⁵⁵⁻⁵⁷ Thus, we computed both structures and took an average of the NMR shielding values. However, both structures lead to similar shieldings for all ¹H, including the S–H (±0.15 ppm), and therefore do not improve the overall correlation, which is also what was observed by Dračinský et al.

<figure><img src="image_11.png" /><figcaption>**Figure 3.** ¹H shieldings (ppm) vs experimental ¹H chemical shifts and the linear fit calculated at the DLPNO-DSD-PBEP86 level using the values from Table 4.</figcaption></figure>

Other hydrogen atoms show significant deviations from the fit, namely, the O–H from L-serine and L-threonine and the H-α, H-β₁, and H-β₂ of L-aspartic acid. Note that for several protons the deviations are mostly independent of the method used (see Figure 4). This suggests they are likely due to issues other than the electronic structure, like the reference geometries or vibrational effects. Relativistic spin–orbit (SO) effects are unlikely to be the source of these errors, although they can explain the large error for the S–H proton. Results calculated using the ΔσSO corrections from ref 18 are presented in the SI. It is apparent that none of the other nuclei (for which ΔσSO is below 0.1 ppm) benefit from the correction, while the fit parameters and other statistical quantities are only marginally improved.

Overall, the smallest deviations are obtained with the TPSS functional, while results using perturbative correlation corrections (DLPNO-DSD-PBEP86 and DLPNO-MP2) are slightly worse than those using other standard DFT functionals. This is in contrast to what has been found in previous studies for heavier nuclei, where DH-DFT often is an improvement over GGA or hybrid functional DFT,²⁶ but there are indications in the literature that error for hydrogen shifts can be quite different than for heavier nuclei, particularly when comparing to experimental, rather than ab initio, reference data.⁵⁸ It is also worth noting that, even using the DLPNO approach, MP2 and DH calculations take about 10 times longer than PBE or TPSS calculations and about 5 times longer than the B3LYP ones (see Figure S5 in the SI).

Statistical parameters for the chemical shifts calculated with respect to TMS are shown in Table S5 and errors for individual nuclei in Figure S3 in the SI. If a reference compound is chosen to evaluate and compare relative chemical shifts, rather than doing a linear fit, there are significant systematic deviations in these results, as can be expected. This can be due to the very different accuracies of the calculated shieldings for the amino acids and the reference system, for example, caused by the different electronic structure and treatment of the environment. The reference shielding at the PBE level is 31.35 ppm, which differs from the “optimal” value of 30.73 ppm, given by the intercept of the linear fit, by 0.62 ppm. This accounts for the MSE and MAE and is almost half of the MaxAE. This demonstrates that the linear fit removes most of
<figure><img src="image_12.png" /><figcaption>**Figure 4.** Errors between experimental ¹H shifts and those calculated with each method from the data in Table 4 using the linear fit parameters in Table 2.</figcaption></figure>

<table><caption>**Table 5.** Correlation Parameters between Experimental ¹H Shifts and the Calculated NMR Shieldings Corresponding to the Values from Table 4 for PBE and DLPNO-DSD-PBEP86 and PBE Values from the Cluster Scheme plus the Difference between the Calculated Values at the PBE and DSD-PBEP86 Levels of the Single Molecules (QM1)ª</caption><thead><tr><th></th><th>slope</th><th>intercept</th><th>R<sup>2</sup></th><th>MAE</th><th>MaxAE<sup>b</sup></th><th>SDE</th></tr></thead><tbody><tr><td>PBE</td><td>-1.00</td><td>30.73</td><td>0.9911</td><td>0.22</td><td>0.72</td><td>SH(cys)</td><td>0.30</td></tr><tr><td>DLPNO-DSD-PBEP86</td><td>-1.04</td><td>31.13</td><td>0.9907</td><td>0.22</td><td>0.79</td><td>SH(cys)</td><td>0.31</td></tr><tr><td>PBE (QM1)</td><td>-0.35</td><td>29.07</td><td>0.5108</td><td>2.17</td><td>8.94</td><td>OH(thr)</td><td>3.13</td></tr><tr><td>DSD-PBEP86 (QM1)</td><td>-0.39</td><td>29.41</td><td>0.6351</td><td>1.59</td><td>6.93</td><td>OH(thr)</td><td>2.43</td></tr><tr><td>PBE + Δσ<sub>DSD-PBEP86</sub></td><td>-1.05</td><td>31.06</td><td>0.9849</td><td>0.28</td><td>1.12</td><td>SH(cys)</td><td>0.40</td></tr></tbody></table>

<table><caption>**Table 6.** Correlation Parameters between Experimental ¹H Shifts and the Calculated NMR Shieldings Corresponding to the Mean of 100 Calculations at the PBE Level at Snapshot Geometries from the Hydrogen MD Simulation (Denoted “PBE-MD”)ª</caption><thead><tr><th></th><th>slope</th><th>intercept</th><th>R<sup>2</sup></th><th>MAE</th><th>MaxAE</th><th>SDE</th></tr></thead><tbody><tr><td>PBE-MD</td><td>-1.04</td><td>30.82</td><td>0.9901</td><td>0.23</td><td>0.88</td><td>OH.ser)</td><td>0.32</td></tr><tr><td>B3LYP + DIFF</td><td>-1.06</td><td>31.35</td><td>0.9929</td><td>0.20</td><td>0.68</td><td>OH.ser)</td><td>0.27</td></tr><tr><td>DLPNO-MP2 + DIFF</td><td>-1.07</td><td>31.33</td><td>0.9936</td><td>0.20</td><td>0.61</td><td>H-β1(asp)</td><td>0.26</td></tr><tr><td>DLPNO-DSD-PBEP86 + DIFF</td><td>-1.07</td><td>31.21</td><td>0.9926</td><td>0.21</td><td>0.60</td><td>H-β1(asp)</td><td>0.28</td></tr></tbody></table>

the systematic error and gives a better estimate of the “best case” performance of each method. Thus, we mostly focus on the linear fit results in the rest of the discussion.

**3.4. ¹H NMR Chemical Shifts—Higher Level Correc- tions from Molecular Calculations.** In the study of Dračinský et al. of the same amino acids, they used the difference of shieldings between the single molecule calculated at the PBE level and those calculated at a higher level of theory (e.g., PBE0, MP2, or CCSD) as corrections to their PBE- GIPAW results. In their study, this correction alone did not lead to an improvement of the correlation with respect to the experimental shifts for ¹H NMR. In order to assess whether such procedures could be efficiently used in combination with our cluster approach, we calculated the NMR shieldings of the single molecules at the PBE and DSD-PBEP86 levels. This allows us to evaluate the correction that should be added to our QM/MM PBE results by comparison with the full embedded cluster results obtained at the DSD-PBEP86 of theory. As can be seen in Table 5, this scheme, at least with a DSD-PBEP86 correction, does not lead to an improvement of

our results, but rather to their deterioration. While the slope and R² values for the full embedded cluster approach for PBE and DSD-PBEP86 are -1.00 and 0.9911 and -1.04 and 0.9907, respectively, the corrected PBE values are -1.05 and 0.9849. We would argue here that this indicates non- transferability of correlation effects from the molecule to the solid due to the effect of the local crystal structure (intermolecular interaction) on the electronic structure. Inspection of the rows labeled “QM1” as well as Table 1 reveals that the molecular results show very large deviations, so it is reasonable that a correction obtained without embedding will yield very different results than the higher level calculation with full embedding. Note that our findings are in agreement with what Dračinský et al. report. In their work, corrections obtained from molecular MP2 and CCSD calculations actually increase the deviations from experiment.

**3.5. ¹H NMR Chemical Shifts—Assessment of Vibra- tional Corrections.** In an attempt to further improve our results, we performed molecular dynamics simulations, allowing only the hydrogen atoms of the QM1 molecule to

<footnote>ªIntercept, MAE, MaxAE, and SDE are given in ppm, while the slope and R² are unitless. ᵇMaximum error and the corresponding hydrogen atom.</footnote>

<footnote>ᵇB3LYP, DLPNO-MP2, and DLPNO-DSD-PBEP86 values correspond to the values from Table 4 plus the difference between the “PBE-MD” values and those from the static calculations at the PBE level.</footnote>
<table><caption>Table 7. Experimental ¹³C Chemical Shifts and the Corresponding Calculated NMR Shieldings in ppm</caption><thead><tr><th></th><th>carbon</th><th>δ<sub>exp</sub><sup>a</sup></th><th>PBE<sup>b</sup></th><th>TPSS<sup>b</sup></th><th>B3LYP<sup>b</sup></th><th>DLPNO-DSD-PBEP86<sup>b,c</sup></th><th>DLPNO-MP2<sup>c</sup></th></tr></thead><tbody><tr><td rowspan="3">L-alanine</td><td>C-α</td><td>50.92</td><td>122.75</td><td>128.71</td><td>122.55</td><td>134.94</td><td>137.95</td></tr><tr><td>COO</td><td>177.71</td><td>-3.28</td><td>4.56</td><td>-8.58</td><td>7.33</td><td>12.23</td></tr><tr><td>C-β</td><td>20.36</td><td>154.14</td><td>160.50</td><td>154.46</td><td>165.91</td><td>169.45</td></tr><tr><td rowspan="3">α-glycine</td><td>COO</td><td>176.25</td><td>-1.78</td><td>5.92</td><td>-7.56</td><td>8.88</td><td>14.35</td></tr><tr><td>C-α</td><td>43.58</td><td>132.05</td><td>137.97</td><td>131.45</td><td>142.82</td><td>145.77</td></tr><tr><td>C-α</td><td>55.69</td><td>118.11</td><td>124.76</td><td>118.04</td><td>130.21</td><td>132.84</td></tr><tr><td rowspan="3">L-serine</td><td>C-β</td><td>62.86</td><td>107.61</td><td>115.51</td><td>109.22</td><td>121.68</td><td>124.33</td></tr><tr><td>COO</td><td>175.05</td><td>0.13</td><td>7.86</td><td>-5.01</td><td>10.60</td><td>16.04</td></tr><tr><td>COO</td><td>175.91</td><td>-1.05</td><td>6.80</td><td>-6.28</td><td>10.43</td><td>16.14</td></tr><tr><td rowspan="3">L-aspartic acid</td><td>C-α</td><td>53.78</td><td>116.04</td><td>122.71</td><td>116.32</td><td>129.40</td><td>133.00</td></tr><tr><td>C-β</td><td>37.77</td><td>137.49</td><td>145.01</td><td>137.91</td><td>149.58</td><td>153.27</td></tr><tr><td>COOH</td><td>174.66</td><td>3.73</td><td>11.35</td><td>-0.88</td><td>15.57</td><td>19.18</td></tr><tr><td rowspan="3">L-cysteine</td><td>C-β</td><td>28.09</td><td>143.55</td><td>149.91</td><td>144.50</td><td>157.26</td><td>160.78</td></tr><tr><td>C-α</td><td>56.01</td><td>117.45</td><td>123.82</td><td>117.45</td><td>129.65</td><td>132.15</td></tr><tr><td>COO</td><td>173.37</td><td>0.39</td><td>8.07</td><td>-5.52</td><td>10.50</td><td>15.74</td></tr><tr><td rowspan="3">L-threonine</td><td>COO</td><td>172.06</td><td>-0.02</td><td>7.70</td><td>-5.85</td><td>10.90</td><td>16.65</td></tr><tr><td>C-α</td><td>61.25</td><td>111.19</td><td>118.34</td><td>111.34</td><td>123.75</td><td>126.31</td></tr><tr><td>C-β</td><td>66.93</td><td>102.23</td><td>109.40</td><td>103.91</td><td>116.63</td><td>118.81</td></tr><tr><td></td><td>C-γ</td><td>20.48</td><td>155.13</td><td>161.52</td><td>155.00</td><td>165.77</td><td>169.34</td></tr></tbody></table>

ᵃExperimental values from ref 18. ᵇValues calculated using the scheme described in Section 2.2 with pcSseg-3 and pcSseg-2 basis sets for QM1 and QM2, respectively. ᶜpcSseg-3 and def2-TZVP basis sets for QM1 and QM2, respectively, were used, and NormalPNO and LoosePNO settings using the fragment scheme presented in Section 2.2 were used.

<table><caption>Table 8. Correlation Parameters between Experimental ¹³C Shifts and the Calculated NMR Shieldings Corresponding to the Values from Table 7 and Various Approximate Schemes (See Text)ᵃ</caption><thead><tr><th></th><th>slope</th><th>intercept</th><th>R<sup>2</sup></th><th>MAE</th><th>MaxAE<sup>b</sup></th><th>SDE</th></tr></thead><tbody><tr><td>PBE</td><td>-0.99</td><td>172.66</td><td>0.9989</td><td>1.66</td><td>4.22</td><td>2.18</td></tr><tr><td>TPSS</td><td>-0.98</td><td>178.96</td><td>0.9990</td><td>1.57</td><td>3.94</td><td>2.10</td></tr><tr><td>B3LYP</td><td>-1.03</td><td>174.90</td><td>0.9993</td><td>1.27</td><td>4.41</td><td>1.74</td></tr><tr><td>DLPNO-DSD-PBEP86</td><td>-1.00</td><td>185.54</td><td>0.9994</td><td>1.14</td><td>4.76</td><td>1.64</td></tr><tr><td>DLPNO-MP2</td><td>-0.99</td><td>187.96</td><td>0.9993</td><td>1.35</td><td>3.40</td><td>1.70</td></tr><tr><td>B3LYP (QM1/QM2)</td><td>-1.03</td><td>174.54</td><td>0.9990</td><td>1.56</td><td>4.99</td><td>2.05</td></tr><tr><td>B3LYP (QM1/MM)</td><td>-1.01</td><td>173.96</td><td>0.9987</td><td>1.88</td><td>4.61</td><td>2.35</td></tr><tr><td>B3LYP (QM1)</td><td>-0.95</td><td>171.02</td><td>0.9977</td><td>2.29</td><td>7.09</td><td>3.14</td></tr><tr><td>PBE (QM1)</td><td>-0.92</td><td>170.76</td><td>0.9956</td><td>3.03</td><td>10.29</td><td>4.31</td></tr><tr><td>DSD-PBEP86 (QM1)</td><td>-0.92</td><td>181.88</td><td>0.9979</td><td>2.31</td><td>6.41</td><td>2.96</td></tr><tr><td>PBE + Δσ<sub>DSD-PBEP86</sub></td><td>-0.99</td><td>183.78</td><td>0.9984</td><td>2.07</td><td>5.09</td><td>2.58</td></tr><tr><td>PBE(GIPAW)<sup>c</sup></td><td>-1.02</td><td>173.05</td><td>0.9995</td><td>1.09</td><td>3.21</td><td>1.49</td></tr></tbody></table>

ᵃIntercept, MAE, MaxAE, and SDE are given in ppm, while the slope and R² are unitless. ᵇMaximum error and the corresponding carbon atom.

ᶜGIPAW values from the work of Dračinský et al.¹⁸

move. We then took the average of NMR shielding calculations from 100 snapshots for each system at the PBE level of theory. This approach has been chosen as it does not require extensive full dynamics simulations or the expensive computation of anharmonic effects but still allows the impact of the vibrational correction to the calculated SCs to be assessed (see PBE-MD in Table 6). Note that zero-point and anharmonic vibrational effects on NMR shieldings should typically be captured using a quantum treatment like perturbation theory.¹⁶,⁵⁹,⁶⁰ However, several studies suggest that MD+snapshot calculations allow for fairly robust estimates of dynamic effects with high efficiency.⁶¹⁻⁶⁴

The difference between the NMR shieldings calculated for the rigid system (Table 1) and the average of the 100 snapshots is then used as a correction (DIFF in Table 6). This correction, calculated at the PBE level, is added to results of the previously calculated NMR shieldings using other levels of theory as can be seen in Table 6. One can see a small

improvement of the MAE and the MaxAE in all cases but for the PBE results, which might hint as some level of error compensation. However, in all cases the slope deviates more from the ideal -1 compared to the static results.

### 3.6. ¹³C NMR Chemical Shifts—Similarities and Differences for Heavier Nuclei.
In the following, we will discuss the calculated carbon NMR shieldings. The results, using the scheme described in Section 2.2, are reported in Table 7 along with the experimental carbon chemical shifts from the work of Dračinský et al.¹⁸ The correlation between experimental shifts and computed shieldings is again obtained with a linear fit, and the fit parameters and statistical analysis are presented in Table 8. A good correlation between the experimental shifts and the calculated shieldings is found for all tested computational methods. In all cases, the slope is found to be very close to the ideal -1. The MAEs are between 1.1 and 1.7 ppm, and the MaxAEs are between 4 and 5 ppm, which is acceptable for molecular solids. Here we find that both the
<figure><img src="image_13.png" /><figcaption>**Figure 5.** Errors between experimental ¹³C shifts and those calculated with each method from the data in Table 7 using the linear fit parameters in Table 8.</figcaption></figure>

MAE and the SDE values improve as the level of theory increases. While the number of values is relatively small, and one has to be careful not to overanalyze the trends observed here, it appears that, especially for the heavier nuclei, post-HF correlation effects improve the results, which is in line with what is known from molecular systems. Our values are slightly worse than those obtained by Dračinský et al. using PBE and GIPAW. This could be due to the fact that, contrary to their study, the position of the carbon atoms was not relaxed in our cluster calculations as relaxing the whole geometry of the crystal is less straightforward in cluster-based calculations. A look at the error of the chemical shifts calculated for individual nuclei in Figure 5 also shows indications for this. As for the ¹H case, large errors for some nuclei are observed, especially in L-aspartic acid, which are mostly independent of the method used. In the future, we might consider combining structures calculated with a periodic approach and compute the NMR SCs with our cluster approach. Note that just as for the hydrogen SCs, molecular DSD-PBEP86 corrections on embedded PBE results do not lead to an improvement of the agreement between computed values and the experiment. Note that, for heavier elements, work by Dračinský et al. reports that such a correction can be beneficial.⁶⁵ However, here the authors corrected GIPAW GGA results with molecular hybrid functional results, and one could also speculate that the observed improvement is rather due to correcting the GIPAW error than the functional error. In our correction scheme, however, we correct embedded Gaussian basis set models using higher level nonembedded Gaussian basis set results, so that mostly the effect of transferability of correlation effects is observed.

## 4. CONCLUSIONS

For ¹H NMR chemical shifts of molecular crystals, we find that an embedded cluster approach yields robust and accurate results compared to experimental values. While a sufficiently large cluster with explicit neighboring molecules is essential, MM embedding beyond that only slightly improves the results. In agreement with previous studies with experimental reference data, we find that functionals like PBE and TPSS yield fairly accurate results, while DH-DFT only yields an improvement for ¹³C NMR chemical shifts. Standard basis sets like Jensen's pcSseg-2 and -3 yield sufficiently converged results, and an

estimate of vibrational corrections based on a simple MD sampling approach shows that these effects are small even for light nuclei like ¹H.

An interesting finding for ¹H NMR chemical shifts is that results from the cluster approach are actually superior to the GIPAW results discussed by Dračinský et al., while this is not entirely true for the ¹³C NMR chemical shifts. This raises the question whether this is due to effects of the partially optimized structures or if the reconstruction of core orbitals in the GIPAW scheme introduces larger errors for the light nuclei. Here, more work comparing appropriately embedded GIAO and periodic GIPAW results for identical systems is needed.

Using the cluster approach in combination with local correlation approximations allows us to test transferability corrections like the one proposed in the work by Dračinský et al., where molecular calculations at the Coupled Cluster level have been used to derive post-DFT corrections to the DFT-GIPAW results. Comparing pure PBE embedded cluster results with results obtained by combining embedded cluster PBE with molecular DSD-PBEP86 corrections shows that the obtained accuracy is notably worse than the actual embedded cluster DSD-PBEP86 values for both ¹H and ¹³C. This hints at a nonadditivity for details of the electronic structure in the molecule and its interaction with the local environment. Hence, from our results it seems advisable to choose the best possible level of theory for an embedded cluster model rather than investing resources in highly accurate molecular calculations to derive corrections.

When new local correlation methods like DLPNO-based Coupled Cluster can be combined with the embedded cluster GIAO scheme presented in this work, this might provide a path to highly accurate calculations with predictive power even for systems as complex as molecular crystals.

### ■ ASSOCIATED CONTENT

**SI** Supporting Information

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.1c01095.

Data analysis, CHELPG charges as a function of the method, ¹H isotropic NMR SCs calculated with looser parameters and correlation fits, linear fit calculated using settings from Table 4, ¹H chemical shifts with spin-orbit
corrections, ¹H chemical shifts using the calculated reference shielding in TMS, ¹³C isotropic NMR SCs calculated with looser parameters and correlation fits, ¹³C chemical shifts using the calculated reference shielding in TMS, shieldings calculated with TPSS and the Dobson treatment of the kinetic energy density, and timings for the calculation needed for the NMR CS calculations using the settings from Table 4 in the main text (PDF)

## ■ AUTHOR INFORMATION

### Corresponding Author

Alexander A. Auer – Max-Planck-Institut für
Kohlenforschung, 45470 Mülheim an der Ruhr, Germany;
[orcid.org/0000-0001-6012-3027];
Email: alexander.auer@kofo.mpg.de

### Authors

Corentin Poidevin – Institut des Sciences Chimiques de
Rennes, 357000 Rennes, France

Georgi L. Stoychev – Max-Planck-Institut für
Kohlenforschung, 45470 Mülheim an der Ruhr, Germany

Christoph Riplinger – FAcTs GmbH, 50677 Köln, Germany

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.jctc.1c01095

### Funding

Open access funded by Max Planck Society.

### Notes

The authors declare no competing financial interest.

## ■ ACKNOWLEDGMENTS

We would like to acknowledge the Max-Planck-Institut für
Kohlenforschung and Frank Neese for generous support.

## ■ REFERENCES

(1) Baias, M.; Widdifield, C. M.; Dumez, J. N.; Thompson, H. P. G.; Cooper, T. G.; Salager, E.; Bassil, S.; Stein, R. S.; Lesage, A.; Day, G. M.; Emsley, L. Powder Crystallography of Pharmaceutical Materials by Combined Crystal Structure Prediction and Solid-State 1H NMR Spectroscopy. *Phys. Chem. Chem. Phys.* **2013**, *15*, 8069–8080.

(2) Pickard, C. J.; Mauri, F. All-Electron Magnetic Response with Pseudopotentials: NMR Chemical Shifts. *Phys. Rev. B - Condens. Matter Mater. Phys.* **2001**, *63*, 245101–245103.

(3) Yates, J. R.; Pickard, C. J.; Mauri, F. Calculation of NMR Chemical Shifts for Extended Systems Using Ultrasoft Pseudopotentials. *Phys. Rev. B - Condens. Matter Mater. Phys.* **2007**, *76*, 1–11.

(4) Teale, A. M.; Lutnæs, O. B.; Helgaker, T.; Tozer, D. J.; Gauss, J. Benchmarking Density-Functional Theory Calculations of NMR Shielding Constants and Spin–Rotation Constants Using Accurate Coupled-Cluster Calculations. *J. Chem. Phys.* **2013**, *138*, 024111.

(5) Flaig, D.; 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.

(6) Stoychev, G. L.; Auer, A. A.; Izsák, R.; Neese, F. Self-Consistent Field Calculation of Nuclear Magnetic Resonance Chemical Shielding Constants Using Gauge-Including Atomic Orbitals and Approximate Two-Electron Integrals. *J. Chem. Theory Comput.* **2018**, *14*, 619–637.

(7) Stoychev, G. L.; Auer, A. A.; Neese, F. Efficient and Accurate Prediction of Nuclear Magnetic Resonance Shielding Tensors with Double-Hybrid Density Functional Theory. *J. Chem. Theory Comput.* **2018**, *14*, 4756–4771.

(8) Hartman, J. D.; Beran, G. J. O. Fragment-Based Electronic Structure Approach for Computing Nuclear Magnetic Resonance

Chemical Shifts in Molecular Crystals. *J. Chem. Theory Comput.* **2014**, *10*, 4862–4872.

(9) Hartman, J. D.; Kudla, R. A.; Day, G. M.; Mueller, L. J.; Beran, G. J. O. Benchmark Fragment-Based 1H, 13C, 15N and 17O Chemical Shift Predictions in Molecular Crystals. *Phys. Chem. Chem. Phys.* **2016**, *18*, 21686–21709.

(10) Hartman, J. D.; Day, G. M.; Beran, G. J. O. Enhanced NMR Discrimination of Pharmaceutically Relevant Molecular Crystal Forms through Fragment-Based Ab Initio Chemical Shift Predictions. *Cryst. Growth Des.* **2016**, *16*, 6479–6493.

(11) Hartman, J. D.; Balaji, A.; Beran, G. J. O. Improved Electrostatic Embedding for Fragment-Based Chemical Shift Calculations in Molecular Crystals. *J. Chem. Theory Comput.* **2017**, *13*, 6043–6051.

(12) Bjornsson, R.; Bühl, M. Modeling Molecular Crystals by QM/MM: Self-Consistent Electrostatic Embedding for Geometry Optimizations and Molecular Property Calculations in the Solid. *J. Chem. Theory Comput.* **2012**, *8*, 498–508.

(13) Holmes, S. T.; Iuliucci, R. J.; Mueller, K. T.; Dybowski, C. Density Functional Investigation of Intermolecular Effects on 13C NMR Chemical-Shielding Tensors Modeled with Molecular Clusters. *J. Chem. Phys.* **2014**, *141*, 164121.

(14) Holmes, S. T.; Iuliucci, R. J.; Mueller, K. T.; Dybowski, C. Critical Analysis of Cluster Models and Exchange-Correlation Functionals for Calculating Magnetic Shielding in Molecular Solids. *J. Chem. Theory Comput.* **2015**, *11*, 5229–5241.

(15) Dittmer, A.; Stoychev, G. L.; Maganas, D.; Auer, A. A.; Neese, F. Computation of NMR Shielding Constants for Solids Using an Embedded Cluster Approach with DFT, Double-Hybrid DFT, and MP2. *J. Chem. Theory Comput.* **2020**, *16*, 6950–6967.

(16) Auer, A. A.; Gauss, J.; Stanton, J. F. Quantitative Prediction of Gas-Phase 13C Nuclear Magnetic Shielding Constants. *J. Chem. Phys.* **2003**, *118*, 10407.

(17) Orrego-Hernández, J.; Hölzel, H.; Quant, M.; Wang, Z.; Moth-Poulsen, K. Scalable Synthesis of Norbornadienes via in Situ Cracking of Dicyclopentadiene Using Continuous Flow Chemistry. *Eur. J. Org. Chem.* **2021**, 2021, 5337–5342.

(18) Dračinský, M.; Vícha, J.; Bártová, K.; Hodgkinson, P. Towards Accurate Predictions of Proton NMR Parameters in Molecular Solids. *ChemPhysChem* **2020**, *21*, 2075.

(19) Reid, D. M.; Collins, M. A. Approximating CCSD(T) Nuclear Magnetic Shielding Calculations Using Composite Methods. *J. Chem. Theory Comput.* **2015**, *11*, 5177–5181.

(20) 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.

(21) Schattenberg, C. J.; Kaupp, M. Extended Benchmark Set of Main-Group Nuclear Shielding Constants and NMR Chemical Shifts and Its Use to Evaluate Modern DFT Methods. *Cite This J. Chem. Theory Comput* **2021**, *17*, 7602–7621.

(22) Schattenberg, C. J.; Reiter, K.; Weigend, F.; Kaupp, M. An Efficient Coupled-Perturbed Kohn-Sham Implementation of NMR Chemical Shift Computations with Local Hybrid Functionals and Gauge-Including Atomic Orbitals. *J. Chem. Theory Comput.* **2020**, *16*, 931–943.

(23) Gauss, J.; Werner, H.-J. NMR Chemical Shift Calculations within Local Correlation Methods: The GIAO-LMP2 Approach. *Phys. Chem. Chem. Phys.* **2000**, *2*, 2083–2090.

(24) Loibl, S.; Schütz, M. NMR Shielding Tensors for Density Fitted Local Second-Order Møller-Plesset Perturbation Theory Using Gauge Including Atomic Orbitals. *J. Chem. Phys.* **2012**, *137*, 084107.

(25) Maurer, M.; Ochsenfeld, C. A Linear- and Sublinear-Scaling Method for Calculating NMR Shieldings in Atomic Orbital-Based Second-Order Møller-Plesset Perturbation Theory. *J. Chem. Phys.* **2013**, *138*, 174104.

(26) Stoychev, G. L.; Auer, A. A.; Gauss, J.; Neese, F. DLPNO-MP2 Second Derivatives for the Computation of Polarizabilities and NMR Shieldings. *J. Chem. Phys.* **2021**, *154*, 164110.
(27) Pinski, P.; Riplinger, C.; Valeev, E. F.; Neese, F. Sparse Maps— A Systematic Infrastructure for Reduced-Scaling Electronic Structure Methods. I. An Efficient and Simple Linear Scaling Local MP2 Method That Uses an Intermediate Basis of Pair Natural Orbitals. *J. Chem. Phys.* **2015**, *143*, 034108.

(28) Pinski, P.; Neese, F. Analytical Gradient for the Domain-Based Local Pair Natural Orbital Second Order Møller-Plesset Perturbation Theory Method (DLPNO-MP2). *J. Chem. Phys.* **2019**, *150*, 164102.

(29) Allen, F. H. The Cambridge Structural Database: A Quarter of a Million Crystal Structures and Rising. *Acta Crystallogr. Sect. B Struct. Sci.* **2002**, *S8*, 380–388.

(30) Neese, F. The ORCA Program System. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2012**, *2*, 73–78.

(31) Neese, F.; Wennmohs, F.; Becker, U.; Riplinger, C. The ORCA Quantum Chemistry Program Package. *J. Chem. Phys.* **2020**, *152*, 224108.

(32) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. *Phys. Rev. Lett.* **1996**, *77*, 3865–3868.

(33) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. *J. Chem. Phys.* **1993**, *98*, 5648–5652.

(34) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Constantin, L. A.; Sun, J. Workhorse Semilocal Density Functional for Condensed Matter Physics and Quantum Chemistry. *Phys. Rev. Lett.* **2009**, *103*, 10–13.

(35) 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.

(36) Kozuch, S.; Martin, J. M. L. Spin-Component-Scaled Double Hybrids: An Extensive Search for the Best Fifth-Rung Functionals Blending DFT and Perturbation Theory. *J. Comput. Chem.* **2013**, *34*, 2327–2344.

(37) Kendall, R. A.; Früchtl, H. A. The Impact of the Resolution of the Identity Approximate Integral Method on Modern Ab Initio Algorithm Development. *Theor. Chem. Acc.* **1997**, *97*, 158–163.

(38) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, Approximate and Parallel Hartree–Fock and Hybrid DFT Calculations. A ‘Chain-of-Spheres’ Algorithm for the Hartree–Fock Exchange. *Chem. Phys.* **2009**, *356*, 98–109.

(39) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. *J. Comput. Chem.* **2011**, *32*, 1456–1465.

(40) Weigend, F.; Ahlrichs, R. *Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy Electronic Supplementary Information (ESI) Available: [DETAILS]* **2005**, *7*, 3297–3305.

(41) Jensen, F. Segmented Contracted Basis Sets Optimized for Nuclear Magnetic Shielding. *J. Chem. Theory Comput.* **2015**, *11*, 132–138.

(42) Dobson, J. F. Alternative Expressions for the Fermi Hole Curvature. *J. Chem. Phys.* **1993**, *98*, 8870–8872.

(43) Reimann, S.; Ekström, U.; Stopkowicz, S.; Teale, A. M.; Borgoo, A.; Helgaker, T. The Importance of Current Contributions to Shielding Constants in Density-Functional Theory. *Phys. Chem. Chem. Phys.* **2015**, *17*, 18834–18842.

(44) Schattenberg, C. J.; Kaupp, M. Effect of the Current Dependence of Tau-Dependent Exchange-Correlation Functionals on Nuclear Shielding Calculations. *J. Chem. Theory Comput.* **2021**, *17*, 1469–1479.

(45) Dittmer, A.; Izsák, R.; Neese, F.; Maganas, D. Accurate Band Gap Predictions of Semiconductors in the Framework of the Similarity Transformed Equation of Motion Coupled Cluster Theory. *Inorg. Chem.* **2019**, *S8*, 9303–9315.

(46) Kaupp, M.; Malkina, O. L.; Malkin, V. G.; Pyykkö, P. How Do Spin-Orbit-Induced Heavy-Atom Effects on NMR Chemical Shifts Function? Validation of a Simple Analogy to Spin-Spin Coupling by Density Functional Theory (DFT) Calculations on Some Iodo Compounds. *Chem. Eur. J.* **1998**, *4*, 118–126.

(47) Vícha, J.; Novotný, J.; Komorovsky, S.; Straka, M.; Kaupp, M.; Marek, R. Relativistic Heavy-Neighbor-Atom Effects on NMR Shifts: Concepts and Trends Across the Periodic Table. *Chem. Rev.* **2020**, *120*, 7065–7103.

(48) Sparta, M.; Retegan, M.; Pinski, P.; Riplinger, C.; Becker, U.; Neese, F. Multilevel Approaches within the Local Pair Natural Orbital Framework. *J. Chem. Theory Comput.* **2017**, *13*, 3198–3207.

(49) Bannwarth, C.; Caldeweyher, E.; Ehlert, S.; Hansen, A.; Pracht, P.; Seibert, J.; Spicher, S.; Grimme, S. Extended Tight-Binding Quantum Chemistry Methods. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2021**, *11*, e1493.

(50) Bannwarth, C.; Ehlert, S.; Grimme, S. GFN2-XTB—An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions. *J. Chem. Theory Comput.* **2019**, *15*, 1652–1671.

(51) Barone, V.; Cossi, M. Conductor Solvent Model. *J. Phys. Chem. A* **1998**, *102*, 1995–2001.

(52) Garcia-Ratés, M.; Neese, F. Effect of the Solute Cavity on the Solvation Energy and Its Derivatives within the Framework of the Gaussian Charge Scheme. *J. Comput. Chem.* **2020**, *41*, 922–939.

(53) *High Resolution NMR; Elsevier: 2000; DOI: 10.1016/b978-0-12-084662-7.x5044-3.*

(54) Jameson, A. K.; Jameson, C. J. Gas-Phase 13C Chemical Shifts in the Zero-Pressure Limit: Refinements to the Absolute Shielding Scale for 13C. *Chem. Phys. Lett.* **1987**, *134*, 461–466.

(55) Kerr, K. A.; Ashmore, J. P.; Koetzle, T. F. A Neutron Diffraction Study of L-Cysteine. *Acta Crystallogr. Sect. B Struct. Crystallogr. Cryst. Chem.* **1975**, *31*, 2022–2026.

(56) Kolesov, B. A.; Minkov, V. S.; Boldyreva, E. V.; Drebushchak, T. N. Phase Transitions in the Crystals of L- and DL-Cysteine on Cooling: Intermolecular Hydrogen Bonds Distortions and the Side-Chain Motions of Thiol-Groups. I. L-Cysteine. *J. Phys. Chem. B* **2008**, *112*, 12827–12839.

(57) Moggach, S. A.; Clark, S. J.; Parsons, S. L-Cysteine-1 at 30 K. *Acta Crystallogr. Sect. E Struct. Reports Online* **2005**, *61*, o2739–o2742.

(58) de Oliveira, M. T.; Alves, J. M. A.; Braga, A. A. C.; Wilson, D. J. D.; Barboza, C. A. Do Double-Hybrid Exchange-Correlation Functionals Provide Accurate Chemical Shifts? A Benchmark Assessment for Proton NMR. *J. Chem. Theory Comput.* **2021**, *17*, 6876–6885.

(59) Ishii, Y.; Terao, T.; Hayashi, S. Theory and Simulation of Vibrational Effects on Structural Measurements by Solid-State Nuclear Magnetic Resonance. *J. Chem. Phys.* **1997**, *107*, 2760.

(60) Ruud, K.; Astrand, P. O.; Taylor, P. R. Zero-Point Vibrational Effects on Proton Shieldings: Functional-Group Contributions from Ab Initio Calculations. *J. Am. Chem. Soc.* **2001**, *123*, 4826–4833.

(61) Benzi, C.; Crescenzi, O.; Pavone, M.; Barone, V. Reliable NMR Chemical Shifts for Molecules in Solution by Methods Rooted in Density Functional Theory. *Magn. Reson. Chem.* **2004**, *42*, S57–S67.

(62) Dračinský, M.; Hodgkinson, P. A Molecular Dynamics Study of the Effects of Fast Molecular Motions on Solid-State NMR Parameters. *CrystEngComm* **2013**, *15*, 8705–8712.

(63) Grigoleit, S.; Bühl, M. Thermal Effects and Vibrational Corrections to Transition Metal NMR Chemical Shifts. *Chem. - A Eur. J.* **2004**, *10*, 5541–5552.

(64) Auer, A. A.; Tran, V. A.; Sharma, B.; Stoychev, G. L.; Marx, D.; Neese, F. A Case Study of Density Functional Theory and Domain-Based Local Pair Natural Orbital Coupled Cluster for Vibrational Effects on EPR Hyperfine Coupling Constants: Vibrational Perturbation Theory versus Ab Initio Molecular Dynamics. *Mol. Phys.* **2020**, *118*, e1797916.

(65) Dračinský, M.; Unzueta, P.; Beran, G. J. O. Improving the Accuracy of Solid-State Nuclear Magnetic Resonance Chemical Shift Prediction with a Simple Molecular Correction. *Phys. Chem. Chem. Phys.* **2019**, *21*, 14992–15000.