Calculating NMR Shieldings and Chemical Shifts
Background
Nuclear Magnetic Resonance (NMR) spectroscopy is a popular experimental technique as it can reveal information about the structure and bonding environment of disordered molecules where techniques like X-Ray Diffraction may not be possible. In experimental NMR, isotopes with non-zero spin reorienting themselves in a magnetic field can be observed as a spectrum of chemical shifts in relation to a reference value. The value of each shift in this spectrum is highly influenced by the electronic environment surrounding a nucleus and thus can reveal deep information about bonding environments and structure. Distinct spectra, each encoding different information, can be collected for the chemical shifts of various isotopes including 1H, 13C, and 15N. However, the interpretation of spectra is complex and often not well-defined. One major motivation of simulated spectra is improving our ability to interpret experimental results.
Theoretically, NMR shielding tensors can be simulated. This tutorial explains the basis of such computations and the "Nuclear Magnetic Resonance (NMR) Spectra" tutorial discusses transforming these results into spectra and techniques for improving the computed results. As mentioned previously, all computed chemical shifts are in relation to some reference value. The most common reference value for 1H is tetramethylsilane. This tutorial will compute shieldings in the absense of a reference to compare to benchmark CCSD(T) values, but similar steps may be followed for computing the shielding tensor of a reference compound. With a reference value, the chemical shift of an atom is defined as the computed shielding minus the reference shielding. In general, the magnitude of a NMR response is proportional to the applied magnetic field. Thus, NMR values (including the ones computed by FHI-aims) are reported in the standardized ppm (parts per million) with respect to the applied field.
Tutorial
In this tutorial, we will demonstrate how to calculate the NMR shielding tensors for pyrazole (C3H4N2). Molecules containing pyrazole rings are prevalent across medicinal, agircultural, and manufacturing applications so better understanding the structure of pyrazole has wide reaching impacts.
Figure 1: Pyrazole (Geometry from NS372 database, Schattenberg and Kaupp)
Step 1: Setup
The geometry of pyrazole was taken from the NS372 dataset, converted to the geometry.in format for FHI-aims and prepared for the shielding calculation.
There are two steps to define the exact magnetic response observables to be computed by FHI-aims.
(1) In control.in, the magnetic_response keyword defines the overall type of desired response observables. A standard control.in file for a shielding tensor computation minimally includes a choice of functional (in this tutorial, PBE), whether relativistic effects should be considered (in this tutorial, no relativity is considered - solely for the reason to enable an exact comparison to other non-relativistic literature results), and basis sets for each atom. Literature, including Schattenberg and Kaupp, discusses the performance of different density functionals for such calculations. For molecules such as pyrazole, the inclusion of scalar relativity likely has a moderately small impact on the computed values, but this impact may become more important for molecules containing heavier atoms. The keyword basis_threshold is required for calculating the magnetic shieldings, although this is generally not recommended for any other calculations because it overrides the check for a valid, well conditioned basis set. Usually this is not a problem - but please keep an eye out to ensure that your basis set is adequate.
# control.in
xc pbe
relativistic none
basis_threshold 0.0
magnetic_response s
+ [basis sets for H, C, and N]
(2) In geometry.in, the magnetic_response keyword after a particular atom line instructs FHI-aims to compute the desired response quantity (here, the shielding) for the preceding atom. Experimentally, only a single isotope NMR spectra can be collected at a single time. However, computationally, the shielding tensor of every atom may be computed simultaneously in this fashion. It is also possible to just request shieldings for a subset of all atoms (saves computer time).
# geometry.in
atom 0.000000 0.000000 0.000000 C
magnetic_response
atom 0.000000 0.000000 1.376234 C
magnetic_response
atom 1.297189 0.000000 1.754244 N
magnetic_response
atom 2.152217 0.000000 0.719064 N
magnetic_response
atom 1.364508 0.000000 -0.351110 C
magnetic_response
atom 1.664015 0.000000 2.685504 H
magnetic_response
atom -0.798075 0.000000 2.094906 H
magnetic_response
atom 1.804849 0.000000 -1.331709 H
magnetic_response
atom -0.854000 0.000000 -0.649833 H
magnetic_response
(3) Basis set choice
The short answer is: Given the benchmark results in Laasner et al. 2024 (Figure 4), reasonable precision of calculated chemical shifts associated with light elements can be expected from
-
FHI-aims "tight" / tier2 settings (around 3 ppm precision on average)
-
NAO-VCC-3Z (better than 2 ppm precision on average)
-
NAO-VCC-4Z (around 1 ppm precision on average)
The precision for 1H may in fact be better than the above average, which includes other elements as well.
(4) More extensive data on basis set convergence
Specific to FHI-aims, the work Laasner et al. 2024 discusses in depth the impact of basis sets on the computation of various magnetic response parameters. For the purposes of this tutorial, we will compute chemical shieldings using 10 different basis set choices as can be found in the following table. All paths are defined relative to the FHIaims base directory.
| Basis Set | Path to species default files |
|---|---|
| FHI-aims "light" | ./species_defaults/defaults_2020/light |
| FHI-aims "tight" | ./species_defaults/defaults_2020/tight |
| NAO-VCC-2Z | ./species_defaults/NAO-VCC-nZ/NAO-VCC-2Z |
| NAO-VCC-3Z | ./species_defaults/NAO-VCC-nZ/NAO-VCC-3Z |
| NAO-VCC-4Z | ./species_defaults/NAO-VCC-nZ/NAO-VCC-4Z |
| NAO-VCC-5Z | ./species_defaults/NAO-VCC-nZ/NAO-VCC-5Z |
| NAO-J-2 | ./species_defaults/NAO-J-n/NAO-J-2 |
| NAO-J-3 | ./species_defaults/NAO-J-n/NAO-J-3 |
| NAO-J-4 | ./species_defaults/NAO-J-n/NAO-J-4 |
| NAO-J-5 | ./species_defaults/NAO-J-n/NAO-J-5 |
Note that, in actual production calculations, it will not be necessary to probe that many basis sets.
For example, if shieldings are the only desired outcome, the NAO-VCC-4Z basis set gives shieldings with a numerical precision of around 1 ppm across a wide range of different systems. If shieldings and J-couplings are to be computed simultaneously for the atom in question, the more expensive NAO-J-4 basis set should be used instead with very similar precision for shieldings. (J-couplings should not be calculated with basis sets that were not specifically devised and tested for their computation.)
Alternatively, the NAO-VCC-3Z or NAO-J-3 basis sets will also give very reasonable results at significantly lower cost.
Step 2: Simulation and Locating Results
Now we can run the simulation! When it has completed, open the aims.out file and search for:
Results for the magnetic response
=================================
NMR shielding tensors (ppm)
---------------------------
Scrolling down you should be able to see the NMR shielding tensors for each atom (1,2,3 ...). For each atom, the diamagnetic, paramagnetic, and total shielding tensor is displayed. For the purposes of this tutorial, only the total shielding tensor, and in particular the isotropic value immediately above this tensor is relevant.
Example:
atom 9 (H):
Diamagnetic: 3.220176E+01 (isotropic)
3.701565E+01 0.000000E+00 0.000000E+00
0.000000E+00 2.648884E+01 0.000000E+00
0.000000E+00 0.000000E+00 3.310078E+01
Paramagnetic: -7.255439E+00 (isotropic)
-1.154549E+01 -6.084078E-15 -8.477800E+00
-2.396032E-13 -3.203120E+00 8.188013E-14
-8.117848E+00 -4.434096E-14 -7.017705E+00
Total: 2.494632E+01 (isotropic)
2.547016E+01 0.000000E+00 0.000000E+00
0.000000E+00 2.328572E+01 0.000000E+00
0.000000E+00 0.000000E+00 2.608308E+01
Anisotropy: 1.705140E+00
Asymmetry: 1.921638E+00
Span: 2.797360E+00
Skew: 5.617874E-01
In this case (FHI-aims "light" species defaults), the shielding value for atom 9 is 24.95 ppm. Tetramethylsilane (TMS) is the standard reference value for 1H NMR spectra. Following a similar process, the shielding value for TMS is ~31.03 ppm. Thus, the chemical shift computed for this atom is 24.95 - 31.03 = -6.08. This chemical shift is the value which could be observed in an experimental NMR spectra. In fact, in the experimental spectra for pyrazole, the observed chemical shift is 6.10ppm, validating this methodology (Chemical BooK).
In general, "light" species defaults are not always expected to yield this degree of precision - see the remarks above.
Step 3: Full Results and Comparison to Literature Values
The computed values for all 10 basis sets considered, as well as reference values from Schattenberg and Kaupp can be found in the following table. Two different reference values are given, first the "CCSD(T)" reference which is considered the gold standard for reference calculations, and second the PBE reference. Note the drastic difference in runtime between the different basis set choices. Note also that the PBE value is of course not supposed to approach the CCSD(T) value - the latter is considered to be more accurate compared to experimental results but also considerably more expensive to obtain than DFT values. All calculations were computed locally on my MacBook pro with an M1 Max processor. All values except where otherwise specified are given in ppm.
| Basis Set | Runtime (s) | Atom #1: C | Atom #2: C | Atom #3: N | Atom #4: N | Atom #5: C | Atom #6: H | Atom #7: H | Atom #8: H | Atom #9: H |
|---|---|---|---|---|---|---|---|---|---|---|
| Reference: CCSD(T) | unknown | 86.6 | 66.1 | 61.7 | -50.0 | 50.6 | 21.85 | 24.02 | 23.85 | 25.07 |
| Reference: PBE | unknown | 69.1 | 51.2 | 34.7 | -77.9 | 34.2 | 21.72 | 23.88 | 23.70 | 24.87 |
| FHI-aims "light" | 1.187 s | 70.99476 | 53.23562 | 35.68539 | -69.01069 | 33.24826 | 22.09092 | 23.98312 | 24.00648 | 24.94632 |
| FHI-aims "tight" | 12.837 s | 69.12302 | 51.40940 | 32.21111 | -75.68031 | 33.39413 | 21.71748 | 23.90323 | 23.73160 | 24.88958 |
| NAO-VCC-2Z | 8.637 s | 67.23049 | 49.36693 | 33.62777 | -79.65490 | 34.56562 | 22.22658 | 24.01489 | 23.89183 | 25.03372 |
| NAO-VCC-3Z | 19.026 s | 69.44204 | 51.83627 | 34.72863 | -76.55157 | 35.02729 | 21.79910 | 23.92784 | 23.72919 | 24.89323 |
| NAO-VCC-4Z | 47.725 s | 69.07875 | 51.16469 | 34.88937 | -77.96371 | 34.44263 | 21.73514 | 23.88442 | 23.71809 | 24.86874 |
| NAO-VCC-5Z | 110.357 s | 68.94693 | 51.10425 | 34.76173 | -77.60961 | 34.23786 | 21.72114 | 23.89316 | 23.71211 | 24.88210 |
| NAO-J-2 | 47.950 s | 67.19014 | 49.27973 | 33.49269 | -79.84363 | 34.44566 | 22.22864 | 24.00826 | 23.88575 | 25.02808 |
| NAO-J-3 | 93.680 s | 69.42041 | 51.75930 | 34.60328 | -76.74272 | 34.92740 | 21.79628 | 23.92414 | 23.72464 | 24.89047 |
| NAO-J-4 | 216.995 s | 69.03964 | 51.07729 | 34.78151 | -78.15577 | 34.33491 | 21.73727 | 23.88015 | 23.71397 | 24.86458 |
| NAO-J-5 | 515.281 s | 68.91522 | 51.02622 | 34.64721 | -77.79539 | 34.13676 | 21.72201 | 23.88978 | 23.70782 | 24.87896 |
In addition, the error with respect to the CCSD(T) reference can be computed for each computed result. The error is comparable to the error of the reference PBE values and thus can be attributed almost entirely to the difference in description between DFT with PBE and the much more computationally intensive CCSD(T). In this case, the error in the nitrogen and carbon dominates the total error whereas the hydrogen error is rather small, indicating that 1H can certainly be expected to be captured well by DFT. In general, there is a trend of increasing error with increasing atomic size when computing NMR shielding tensors.
| Basis Set | MAE (total) | MAE (H) | MAE (N) | MAE(C) |
|---|---|---|---|---|
| Reference: CCSD(T) | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| Reference: PBE | 11.5911 | 0.1550 | 27.4500 | 16.2667 |
| FHI-aims "light" | 10.1561 | 0.1395 | 22.5127 | 15.2738 |
| FHI-aims "tight" | 11.6768 | 0.1370 | 27.5846 | 16.4578 |
| NAO-VCC-2Z | 12.2582 | 0.1149 | 28.8636 | 17.3790 |
| NAO-VCC-3Z | 11.2176 | 0.1102 | 26.7615 | 15.6648 |
| NAO-VCC-4Z | 11.5524 | 0.1459 | 27.3872 | 16.2046 |
| NAO-VCC-5Z | 11.5711 | 0.1454 | 27.2739 | 16.3370 |
| NAO-J-2 | 12.3226 | 0.1170 | 29.0255 | 17.4615 |
| NAO-J-3 | 11.2763 | 0.1136 | 26.9197 | 15.7310 |
| NAO-J-4 | 11.6129 | 0.1485 | 27.5371 | 16.2827 |
| NAO-J-5 | 11.6290 | 0.1479 | 27.4241 | 16.4073 |