Part 1: NMR - chemical shielding ¶

Content¶

 1 Converging chemical shielding in diamond
 2 Mapping the induced current
 3 Comparing to experimental chemical shifts
 4 Predicting experiment from theory

1 Converging chemical shielding in diamond

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• Calculate the chemical shielding
• Converge the chemical shielding with respect to plane-wave energy cutoff and k-point mesh
• Compare convergence of chemical shielding to that of energy

1.1 Task¶

Calculate the chemical shielding for 2-atomic primitive cell carbon diamond, and converge with respect to plane-wave energy cutoff and k-points.

In a magnetic field, nuclei with a non-zero spin) $\textbf{I}$ have a magnetic moment $\boldsymbol{\mu}$ that precesses (cf. the red arrow on the double inverted cone below) about an external magnetic field $\textbf{B}_{ext}$. It precesses at its Larmor frequency $\omega_L$, which is determined by the strength of the external magnetic field $\textbf{B}_{ext}$ and the nucleus' gyromagnetic ratio $\gamma$. If an oscillating magnetic field, i.e. a radio-frequency (RF) pulse, at frequency $\omega_{rf}$ is applied perpendicular (transverse frame) to the constant, external magnetic field (reference frame), then the nuclear magnetic moment changes. If RF $\omega_{rf}$ is similar to $\omega_{L}$, then magnetic resonance occurs. In other words, the nuclear magnetic moment flips from parallel to perpendicular to the magnetic field, i.e., the reference to the transverse frame. The subsequent relaxation of the magnetic moment back to the reference frame creates a signal that is measured in nuclear magnetic resonance (NMR).

Seemingly identical nuclear isotopes do not necessarily oscillate at the same frequency, despite having the same $\gamma$, due to different chemical environments. The shells of electrons surrounding the nuclei differ and shield the nuclei. For instance, the ¹H of the methylene -CH₂- group is closer in space to Br compared to the methyl groups -CH₃ in bromoethane Br-CH₂-CH₃. This shielding changes the magnetic field at the nucleus and thereby the frequency at which resonance occurs. The resultant shift in frequency relative to a standard reference is called chemical shift $\delta$ and enables the probing of the chemical structure of molecules and solid state systems using NMR.

The chemical shift $\delta$ is measured in NMR experiments and expresses the chemical shielding $\sigma$ relative to that of a reference compound $\sigma^{ref}$ (the standard reference is tetramethylsilane, TMS):

$$\tag{1} \begin{equation} \delta_{ij}(\mathbf{R}) = \sigma_{ij}^{\mathrm{ref}} - \sigma_{ij}(\mathbf{R}) \end{equation} $$

where $ij$ go over Cartesian directions, and $\mathbf{R}$ is the atomic nuclear positions.

The chemical shielding tensor may be calculated using linear response:

$$\tag{2} \begin{equation} \sigma_{ij}(\mathbf{R}) = - \frac{ \partial B^{\mathrm{ind}}_i(\mathbf{R})}{ \partial B^{\mathrm{ext}}_j} \end{equation} $$

where $B^{\mathrm{ind}}$ is the induced magnetic field.

In this exercise, you will calculate the isotropic chemical shielding for C in a diamond structure, converging for the plane-wave energy cutoff and the $\textbf{k}$-point mesh.

1.2 Input¶

The input files to run this example are prepared at $TUTORIALS/NMR/e01_shielding.

VASP looks in the current directory for four main input files, i.e., POSCAR, INCAR, KPOINTS, and POTCAR. Check them out!

POSCAR

 C diamond
   3.56700000000000
     0.0000000000  0.5000000000  0.5000000000
     0.5000000000  0.0000000000  0.5000000000
     0.5000000000  0.5000000000  0.0000000000
   C
     2
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.2500000000000000  0.2500000000000000  0.2500000000000000

INCAR.nmr

ENCUT = 400              # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-8             # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate          # Sets the "precision" mode
LASPH = .TRUE.           # Non-spherical contributions to the gradient of the density in the PAW spheres

LCHIMAG = .TRUE.         # Turns on linear response for chemical shifts
LNMR_SYM_RED = .TRUE.    # Consistent symmetry with star and k-space derivatives
NLSPLINE = .TRUE.        # Differentiable projectors in reciprocal space
KPAR = 4                 # Determines the number of k-points that are to be treated in parallel

KPOINTS.nmr

C diamond - regular N x N x N mesh centered at Gamma
0
Gamma
 4 4 4
 0 0 0

POTCAR

 Standard pseudopotential of C.

Check the tags set in the INCAR file! The first set of tags, ENCUT, ISMEAR, SIGMA, EDIFF, PREC, and LASPH are settings for a normal density-functional-theory (DFT) calculation. Note that the electronic state needs to be determined with a high accuracy (PREC=A, EDIFF≤1E-8, high ENCUT) for NMR.

The remaining tags are used for calculating the chemical shielding. LCHIMAG calculates the chemical shielding through linear response. LNMR_SYM_RED discards symmetry operations inconsistent with how k-space derivatives in the linear response for chemical shieldings are calculated, and NLSPLINE ensures that the PAW projectors in reciprocal space are differentiable. KPAR determines how many $\textbf{k}$-points should be treated in parallel, speeding up the calculation.

The KPOINTS file specifies a 4x4x4 $\mathbf{k}$-point mesh.

1.3 Calculation¶

Open a terminal and navigate to this example's directory by entering the following:

cd $TUTORIALS/NMR/e01_*

We will first converge the chemical shielding with respect to the plane-wave energy cutoff. In order to run VASP using different energy cutoffs, the INCAR file needs to be adjusted. In the following script, ENCUT = 400 eV in INCAR is replaced with a new energy cutoff: 500 to 900 eV using a step-size of 100 eV. In this example's directory enter the following into the terminal to run the encut_shielding.sh script:

bash encut_shielding.sh

encut_shielding.sh runs the calculations and finds the shieldings for ¹³C in the OUTCAR file using the following grep command:

grep -A 5 "iso_shield" OUTCAR | tail -1 | awk '{print $6}'

You can take a look in the stdout at the different stages of the calculation. The first steps are the DFT converging the CHGCAR and WAVECAR:

 LDA part: xc-table for (Slater+PW92), standard interpolation
 found WAVECAR, reading the header
  number of k-points has changed, file:    29 present:     8
  trying to continue reading WAVECAR, but it might fail
 WAVECAR: different cutoff or change in lattice found
 POSCAR, INCAR and KPOINTS ok, starting setup
 FFT: planning ... GRIDC
 FFT: planning ... GRID_SOFT
 FFT: planning ... GRID
 reading WAVECAR
 the WAVECAR file was read successfully
 initial charge from wavefunction
 entering main loop
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.177071554064E+02   -0.17707E+02   -0.55181E+02   128   0.103E+02    0.278E+01
DAV:   2    -0.181025057223E+02   -0.39535E+00   -0.53621E+00   128   0.179E+01    0.193E+01
DAV:   3    -0.180696889204E+02    0.32817E-01   -0.24036E-01   136   0.404E+00    0.313E+00
DAV:   4    -0.180615539186E+02    0.81350E-02   -0.15850E-02   136   0.102E+00    0.901E-01
DAV:   5    -0.180614863890E+02    0.67530E-04   -0.33854E-04   168   0.175E-01    0.356E-02
DAV:   6    -0.180614871993E+02   -0.81038E-06   -0.11823E-06   128   0.120E-02    0.119E-02
DAV:   7    -0.180614877507E+02   -0.55135E-06   -0.22187E-07   168   0.537E-03    0.114E-03
DAV:   8    -0.180614877841E+02   -0.33422E-07   -0.43684E-09    96   0.781E-04    0.214E-04
DAV:   9    -0.180614877734E+02    0.10698E-07   -0.16385E-10    64   0.161E-04    0.430E-05
DAV:  10    -0.180614877719E+02    0.15249E-08   -0.72505E-12    64   0.321E-05
   1 F= -.18061488E+02 E0= -.18061488E+02  d E =-.143697E-45

Followed by the linear response along the magnetic field axes; BDIR is the applied (external) magnetic field and IDIR is the induced magnetic field along the Cartesian axes (1: x-axis, 2: y-axis, 3: z-axis); IQ are the grid points in the finite difference stencil define by ICHIBARE and runs from -ICHIBARE to +ICHIBARE. Since the default is ICHIBARE = 1, there are three IQ points:

BDIR= 1  IDIR= 2  IQ= -1
...
BDIR= 1  IDIR= 3  IQ=  0
...
BDIR= 2  IDIR= 1  IQ=  1
...
...
BDIR= 3  IDIR= 2  IQ=  1
...

Congratulations on running your first chemical shielding calculation!

The shieldings from multiple VASP calculations at different energy cutoffs have been stored in encut_shielding.dat using the encut_shielding.sh script! The chemical shieldings are found in the OUTCAR file:

 ---------------------------------------------------------------------------------------
  CSA tensor (J. Mason, Solid State Nucl. Magn. Reson. 2, 285 (1993))
  for chemical shielding  (including the isotropic core contribution)

  BDIR labels direction of applied magnetic field (i.e. B0)
  For BDIR==1, B0 along x-axis, etc.
  Induced field listed along cartesian directions for each BDIR
 ---------------------------------------------------------------------------------------
                 EXCLUDING G=0 CONTRIBUTION             INCLUDING G=0 CONTRIBUTION
            ------------------------------------   ------------------------------------
  ion  BDIR            X           Y           Z              X           Y           Z
          (   iso_shield        span        skew) (  iso_shield        span        skew)
 ---------------------------------------------------------------------------------------
    1     1        ...         ...          ...         ...           ...         ... 
          2        ...         ...          ...         ...           ...         ...
          3        ...         ...          ...         ...           ...         ... 
          (        ...         ...          ... ) (     ...           ...         ...  )
...
 ---------------------------------------------------------------------------------------
  IF SPAN.EQ.0, THEN SKEW IS ILL-DEFINED
 ---------------------------------------------------------------------------------------

Note that these chemical shielding include the core contributions. These are the final important terms for comparison with experiment, which you will cover in the next exercise.

Let's plot the chemical shieldings at each cutoff using py4vasp!

What does this tell us about how quickly the total energy converges compared to the chemical shielding?

Generally, we do not recommend converging with respect to the total energy. We show it here for the sake of comparison. Instead, an energy difference (e.g., formation energy) or your property of interest (e.g., chemical shielding) should taken and converged.

The energy cutoff is only one of the important parameters, alongside the $\textbf{k}$-point mesh. Here, we use a cutoff of 500 eV while testing the $\textbf{k}$-point mesh to save time. In real applications, parameters should be converged sequentially. Each technical parameter should be converged one after the other, first the energy cutoff, here a cutoff of 800 eV should be used to test the convergence of the $\textbf{k}$-point mesh. If there were further parameters, the converged energy cutoff (800 eV) and $\textbf{k}$-point mesh should be used to converge them. An example is for molecular systems, where the size of the cell (i.e., the vacuum) must be tested to achieve convergence instead of the $\textbf{k}$-point mesh.

In order to run VASP using different $\textbf{k}$-point meshes, the KPOINTS file needs to be updated. For each run, replace 4 4 4 in KPOINTS (where N N N corresponds to a NxNxN $\textbf{k}$-point mesh) with 6 6 6, 8 8 8, ..., 16 16 16, increasing by increments of 2 $\textbf{k}$-points. The script kpoints.sh does this for you. In this example's directory, enter the following into the terminal:

bash kpoints_shielding.sh

After VASP has run, the chemical shieldings can be found in the OUTCAR file. Run multiple VASP calculations with different $\mathbf{k}$-points meshes and store the shieldings in kpoints_shielding.dat! For the sake of saving time, 12x12x12, 14x14x14, and 16x16x16 have been calculated in advance and added in. You can repeat these calculations yourself by changing 4 6 8 10 to 4 6 8 10 12 14 16 and removing the last three echo statements in kpoints_shielding.sh.

Let's plot the chemical shieldings with each $\mathbf{k}$-mesh using py4vasp!

At first glance, the difference in convergence is not so apparent. It becomes clearer when the scale of the y-axes is looked at. The energy is converged to within ~10 meV with a $\textbf{k}$-point mesh of 6x6x6. To converge the chemical shielding to within 0.1 ppm, on the other hand, requires a much denser $\textbf{k}$-mesh, 14x14x14. You can see that the chemical shielding depends more on the number of $\textbf{k}$-points than on the energy cutoff. The difference between 400 and 900 eV was ~1 ppm. The difference between a 4x4x4 and a 16x16x16 $\textbf{k}$-point mesh is almost 90 ppm, nearly two orders of magnitude greater! This tells us that a correct description of the Brillouin zone with a fine $\textbf{k}$-mesh is more impactful on the chemical shielding than the energy cutoff (i.e., the completeness of the basis using plane-waves).

Converging the shielding with respect to the energy cutoff and the $\textbf{k}$-point mesh shows the importance of converging your parameter of interest and not merely the total energy. It is important to converge with respect to all necessary parameters, not just one or the other.

1.4 Questions¶

1. What is the difference between chemical shift and chemical shielding?
2. If a fine $\textbf{k}$-point is required, what does this say about the importance of describing the Brillouin zone?
3. If a high energy cutoff is required, what does this say about the completeness of the basis?

2 Mapping the induced current

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• Calculate the induced current
• Visualize the current

2.1 Task¶

Calculate the induced current for gas-phase benzene in a cubic 8x8x8 ų cell and plot the currents.

A uniform, external magnetic field $\textbf{B}_{ext}$ applied to a material induces an electric current. This induced electric current in turn induces a magnetic field. The chemical shielding tensor $\sigma$ describes the relation between this induced magnetic field $\textbf{B}_{in}^{(1)}$ and the external magnetic field $\textbf{B}_{ext}$:

$$\tag{3} \begin{equation} \textbf{B}_{in}^{(1)}(\textbf{R}) = -\sigma(\textbf{R}) \textbf{B}_{ext} \end{equation} $$

where $\textbf{R}$ is the nuclear position. The induced magnetic field $\textbf{B}_{in}^{(1)}$ can be calculated from the first-order induced current $\textbf{j}^{(1)}(\textbf{r}^{\prime})$ [Phys. Rev. B 63 (2001) 245101]:

$$\tag{4} \begin{equation} \textbf{B}_{in}^{(1)}(\textbf{R}) = \frac{1}{c} \int \textbf{j}^{(1)}(\textbf{r}^{\prime}) \times \frac{\textbf{R} - \textbf{r}^{\prime}}{|\textbf{R} - \textbf{r}^{\prime}|^3} \,d^3r^{\prime} \end{equation} $$

where $\textbf{r}^{\prime}$ is position in space.

The induced current $\textbf{j}^{(1)}(\textbf{r}^{\prime})$ can be expressed as a sum of three parts: a contribution from the pseudized valence wavefunctions $\textbf{j}^{(1)}_{bare}(\textbf{r}^{\prime})$, a diamagnetic term $\textbf{j}^{(1)}_{\Delta d}(\textbf{r}^{\prime})$, and a paramagnetic term $\textbf{j}^{(1)}_{\Delta p}(\textbf{r}^{\prime})$. The latter two account for the deviation of the pseudized valence wavefunction from the true all-electron wavefunction close to the nucleus [Phys. Rev. B 76 (2007) 024401]:

$$\tag{5} \begin{equation} \textbf{j}^{(1)}(\textbf{r}^{\prime}) = \textbf{j}^{(1)}_{bare}(\textbf{r}^{\prime}) + \textbf{j}^{(1)}_{\Delta d}(\textbf{r}^{\prime}) + \textbf{j}^{(1)}_{\Delta p}(\textbf{r}^{\prime}) \end{equation} $$

It is possible to visualize the induced current for a material. In this exercise, you will produce plots of the current for benzene C₆H₆, which has a distinct ring-current pattern.

2.2 Input¶

The input files to run this example are prepared at $TUTORIALS/NMR/e02_induced_current.

VASP looks in the current directory for four main input files, i.e., POSCAR, INCAR, KPOINTS, and POTCAR. Check them out!

POSCAR

Benzene
1.0
        8.0000000000         0.0000000000         0.0000000000
        0.0000000000         8.0000000000         0.0000000000
        0.0000000000         0.0000000000         8.0000000000
    C    H
    6    6
Cartesian
4.000000000       4.000000000       5.397505962
4.000000000       5.210292596       4.698718710
4.000000000       5.210292596       3.301281290
4.000000000       4.000000000       2.602494038
4.000000000       2.789707404       3.301281290
4.000000000       2.789707404       4.698718710
4.000000000       4.000000000       6.489679174
4.000000000       6.156030752       5.244777584
4.000000000       6.156030752       2.755222416
4.000000000       4.000000000       1.510320826
4.000000000       1.843969248       2.755222416
4.000000000       1.843969248       5.244777584

INCAR

ENCUT = 400              # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV

EDIFF = 1E-8             # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate          # Sets the "precision" mode
LASPH = .TRUE.           # Non-spherical contributions to the gradient of the density in the PAW spheres

LCHIMAG = .TRUE.         # Turns on linear response for chemical shifts
LNMR_SYM_RED = .TRUE.    # Consistent symmetry with star and k-space derivatives
NLSPLINE = .TRUE.        # Differentiable projectors in reciprocal space
WRT_NMRCUR = 1           # Prints the currents to NMRCURBX

KPOINTS

Gamma-point k-mesh
0
Gamma
 1 1 1
 0 0 0

POTCAR

 Pseudopotential of C and H.

The INCAR files are identical to those used in Example 1, with the addition of WRT_NMRCUR that ensures that the currents are printed out as NMRCURBX on the x-axis; separate files exist for y- and z-axes.

The KPOINTS file specifies a Γ-point $\mathbf{k}$-mesh. For an isolated molecule, as large a cell is required, as well as a $\Gamma$-only $\mathbf{k}$-mesh. Note that a sufficiently large vacuum around the molecule needs to be used; this vacuum is a parameter that needs to be increased until the parameter of interest (e.g., chemical shielding) is converged.

2.3 Calculation¶

Open a terminal and navigate to this example's directory:

cd $TUTORIALS/NMR/e02_*

Run the VASP calculation using:

mpirun -np 4 vasp_std

Take a look at the benzene molecule that you have just calculate using py4vasp:

You can see that benzene is a highly symmetric molecule, consisting of a hexagon of CH units in a plane. How do you think that this will influence the currents?

The chemical shieldings can, once again, be found in the OUTCAR file. The induced current is found in the NMRCURBX, NMRCURBY, and NMRCURBZ files, depending on whether the magnetic field is applied along x, y, or z.

Let's take the currents and visualize them using the plot_nmrcur_slice.py script. plot_nmrcur_slice.py is a python script that takes four arguments:

• The currents file along the x-, y-, or z-direction: NMRCURBX, NMRCURBY, or NMRCURBZ
• How many n points to include from the current grid, i.e. 1 would be every point, 2 would be every second point, and 3 every third point.
• Where to take a yz-plane, cutting the x-axis in direct coordinates, e.g. 0.5 would be halfway through the cell.
• How much to scale the arrows showing the magnitude of the current by, where 4.5 would scale the arrow length $L$ by $(L / 10^{4.5}) \times W$, where $W$ is the width of the axes.

Extract the current and plot it using plot_nmrcur_slice_c6h6.py, which is the same as plot_nmrcur_slice.py but shows the positions of the atoms in the benzene molecule:

python3 plot_nmrcur_slice_c6h6.py NMRCURBX 2 0.5 0.9

Make sure to adjust the 0.9 to find what length of arrows is easiest to visualize the currents for you. If the arrows are visible, generally only small adjustments are required.

You can then find the plotted currents in fig_nmrcurbx-slice.png.

You can also plot the NMR currents with py4vasp using the following:

You can see that the current moves around the plane of the molecule in a ring. Rings of electrons lie above and below the plane of the molecule, known as aromaticity. In the presence of a magnetic field, they flow around the ring, creating a current. This ring current generates a strong magnetic field aligned with the external field outside the ring and against within. The protons (H atoms) experience very strong deshielding as a result. This is a characteristic feature of the NMR spectra of aromatic compounds.

2.4 Questions¶

1. Which tag enables writing the induced current and where are they written to?
2. Are the H atoms of benzene shielded or deshielded?

3 Comparing to experimental chemical shifts

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• Calculate the chemical shielding for an O-containing series
• Plot the experimental chemical shifts against the calculated chemical shieldings
• Visualize the relation between experimental chemical shifts and calculated chemical shieldings

3.1 Task¶

Calculate chemical shieldings for BaSnO₃, MgO, SrO, and SrTiO₃ and then perform a linear fit to experimental chemical shifts.

The chemical shielding is not an observable. In NMR, the frequency of precession $\omega_L$ is measured. $\omega_L$ depends on the strength of the magnetic field $\textbf{B}$ and the gyromagnetic ratio $\gamma$:

$$\tag{6} \begin{equation} \omega_L = -\gamma B_{z}. \end{equation} $$

Each NMR instrument generates a strong magnetic field varying between 1-20 T in modern NMR spectrometers. The stronger the field, the greater the difference in energy between spin states, increasing the signal-to-noise ratio by increasing the population bias of the more stable spin state.

However, the stronger field also changes the measured sample's frequency $\omega_{sample}$. To allow for comparison between spectrometers, the frequency of a reference compound $\omega_{ref}$ is taken [Pure Appl. Chem. 80 (2008) 59], and the relative chemical shift $\delta$ is calculated:

$$\tag{7} \begin{equation} \delta = \frac{\omega_{ref} - \omega_{sample}}{\omega_{sample}} \times 10^6. \end{equation} $$

Computationally, the chemical shielding $\sigma$, rather than the chemical shift $\delta$, is computed using linear response. The induced currents (cf. Example 2) are used to derive, via the Biot-Savart law, the corresponding induced magnetic fields ($\textbf{B}^{(1)}_{bare}(\textbf{r}^{\prime})$, $\textbf{B}^{(1)}_{\Delta d}(\textbf{r}^{\prime})$, and $\textbf{B}^{(1)}_{\Delta p}(\textbf{r}^{\prime})$ for the bare, diamagnetic, and paramagnetic terms, respectively). The corresponding chemical shieldings are then found ($\sigma_{bare}$, $\sigma_{\Delta d}$, and $\sigma_{\Delta p}$ for the bare, diamagnetic, and paramagnetic terms, respectively) and summed together, along with the chemically invariant core contributions $\sigma_{core}$, to find the total chemical shielding $\sigma_{tot}$:

$$\tag{8} \begin{equation} \sigma_{tot} = \sigma_{bare} + \sigma_{\Delta d} + \sigma_{\Delta p} + \sigma_{core}. \end{equation} $$

In this exercise, you will calculate the chemical shielding for a series of O-containing compounds and plot against the experimental chemical shifts to find the relationship between the two.

3.2 Input¶

The input files to run this example are prepared at $TUTORIALS/NMR/e03_experiment.

VASP looks in the current directory for four main input files, i.e., POSCAR, INCAR, KPOINTS, and POTCAR. Check them out!

BaSnO3/POSCAR

BaSnO3
   4.11442843484200
   1.00000000 0.00000000 0.00000000
   0.00000000 1.00000000 0.00000000
   0.00000000 0.00000000 1.00000000
   Ba Sn O
   1 1 3
Direct
0.50000000 0.50000000 0.50000000
0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.50000000
0.00000000 0.50000000 0.00000000
0.50000000 0.00000000 0.00000000

MgO/POSCAR

MgO
   4.21699828357500
     0.5000000000000000    0.5000000000000000    0.0000000000000000
     0.0000000000000000    0.5000000000000000    0.5000000000000000
     0.5000000000000000    0.0000000000000000    0.5000000000000000
   Mg O
   1 1
Cartesian
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.5000000000000000  0.5000000000000000  0.5000000000000000

SrO/POSCAR

SrO
   5.10399789167400
     0.5000000000000000    0.5000000000000000    0.0000000000000000
     0.0000000000000000    0.5000000000000000    0.5000000000000000
     0.5000000000000000    0.0000000000000000    0.5000000000000000
   Sr O
   1 1
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.5000000000000000  0.5000000000000000  0.5000000000000000

SrTiO3/POSCAR

SrTiO3
   3.94204796873800
   1.00000000 0.00000000 0.00000000
   0.00000000 1.00000000 0.00000000
   0.00000000 0.00000000 1.00000000
   Sr Ti O
   1 1 3
Direct
0.00000000 0.00000000 0.00000000
0.50000000 0.50000000 0.50000000
0.50000000 0.50000000 0.00000000
0.00000000 0.50000000 0.50000000
0.50000000 0.00000000 0.50000000

BaSnO3/INCAR MgO/INCAR SrO/INCAR SrTiO3/INCAR

ENCUT = 500              # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV

EDIFF = 1E-8             # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate          # Sets the "precision" mode
LASPH = .TRUE.           # Non-spherical contributions to the gradient of the density in the PAW spheres

LCHIMAG = .TRUE.         # Turns on linear response for chemical shifts
LNMR_SYM_RED = .TRUE.    # Consistent symmetry with star and k-space derivatives
NLSPLINE = .TRUE.        # Differentiable projectors in reciprocal space

BaSnO3/KPOINTS MgO/KPOINTS SrO/KPOINTS SrTiO3/KPOINTS

C
0
G
 6 6 6
 0 0 0

BaSnO3/POTCAR

 Pseudopotential of Ba, Sn, and O.

MgO/POTCAR

 Pseudopotential of Mg and O.

SrO/POTCAR

 Pseudopotential of Sr and O.

SrTiO3/POTCAR

 Pseudopotential of Sr, Ti, and O.

experiment.dat

MgO 47
BaSnO3 143
SrO 390
SrTiO3 465

The INCAR files are identical to those used in Example 1.

The KPOINTS file specifies a Γ-centered $\mathbf{k}$-mesh. An 6x6x6 $\mathbf{k}$-mesh is used.

3.3 Calculation¶

Open a terminal and navigate to this example's directory:

cd $TUTORIALS/NMR/e03_*

Take a look inside e03_* with ls and you will see four directories: BaSnO3, MgO, SrO, and SrTiO3. Each of these is named with a different O-containing compound and contains the corresponding INCAR, KPOINTS, POSCAR, and POTCAR files.

From the central e03_* directory, enter MgO and calculate the chemical shielding by entering the following into the terminal:

cd $TUTORIALS/NMR/e03_*/MgO
mpirun -np 4 vasp_std

Congratulations on calculating the chemical shielding!

Running all of these calculations is very time-consuming, so we recommend only running MgO.

Optional: If you have enough time and want to try running all the calculations, you can using the instructions below:

cd $TUTORIALS/NMR/e03_*
bash O_shielding.sh

rm -f O_shielding.dat
for a in *O*/
do
cd $a
mpirun -np 4 vasp_std
if [[ $a ==  *"O3"* ]]; then
shielding=$(grep -A 15 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
else
shielding=$(grep -A 10 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
fi
cd ..
echo ${a///} $shielding >> O_shielding.dat
done

The chemical shieldings can be found in the OUTCAR file. Using what you learned in Example 1, take a look inside and see if you can find the isotropic shielding for the O atom in MgO.

 ---------------------------------------------------------------------------------------
  CSA tensor (J. Mason, Solid State Nucl. Magn. Reson. 2, 285 (1993))
  for chemical shielding  (including the isotropic core contribution)

  BDIR labels direction of applied magnetic field (i.e. B0)
  For BDIR==1, B0 along x-axis, etc.
  Induced field listed along cartesian directions for each BDIR
 ---------------------------------------------------------------------------------------
                 EXCLUDING G=0 CONTRIBUTION             INCLUDING G=0 CONTRIBUTION
            ------------------------------------   ------------------------------------
  ion  BDIR            X           Y           Z              X           Y           Z
          (   iso_shield        span        skew) (  iso_shield        span        skew)
 ---------------------------------------------------------------------------------------

    2     1     187.2125      0.0000      0.0000       198.9458      0.0000      0.0000
          2       0.0000    187.2125     -0.0000         0.0000    198.9458     -0.0000
          3      -0.0000      0.0000    187.2125        -0.0000      0.0000    198.9458
          (     187.2125      0.0000      0.0000) (    198.9458      0.0000      0.0000)

You can extract the shielding with the following command:

cd $TUTORIALS/NMR/e03_*/MgO
shielding=$(grep -A 10 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
echo "MgO" $shielding >> ../O_shielding.dat

Once you have extracted the chemical shieldings for MgO add those for BaSnO3, SrO, and SrTiO3 by entering the following into the terminal:

cd $TUTORIALS/NMR/e03_*
echo "BaSnO3 74.0956" >> O_shielding.dat
echo "SrO -219.9072" >> O_shielding.dat
echo "SrTiO3 -281.9145" >> O_shielding.dat

$$\tag{9} \begin{equation} \delta_{\mathrm{iso}}[N] = \delta_{\mathrm{ref}}[N] - \sigma_{\mathrm{iso}}[N]. \end{equation} $$

In an experiment, the reference has been chosen, e.g., tetramethylsilane (TMS) for ¹H, due to it having a distinct NMR peak or useful chemical properties.

In a calculation, there is no such standard reference. Any individual calculation used as a reference has the potential to distort any derived chemical shift due to limitations in the computational method or error. Instead, several experimental chemical shifts can be plotted against calculated chemical shieldings to find the relationship between the two, averaging out the effects of any individual calculation [J. Chem. Phys. 146 (2017) 064115].

$$\tag{10} \begin{equation} \delta_{\mathrm{iso}}^{\mathrm{exp}}[N] = \delta_{\mathrm{ref}}[N] + m \sigma_{\mathrm{iso}}^{\mathrm{calc}}[N] \end{equation}, $$

where $\delta_{\mathrm{iso}}^{\mathrm{exp}}$ is the experimental chemical shift, $\delta_{\mathrm{ref}}$ is that of the reference chemical shift (i.e. the y-intercept), $m$ is the gradient, $\sigma_{\mathrm{iso}}^{\mathrm{calc}}$ is the calculating chemical shielding, and $N$ is the element of interest, e.g. O. Try plotting it using the scripts below:

Note that the Pearson coefficient $r$ allows you to assess the quality of your fit. If its absolute value is below ~0.95, check that the order of your compounds is the same in O_shielding.dat and experiment.dat to ensure that you are plotting the correct points against one another.

Now that you have a plot, you can try to take the chemical shielding for a compound and predict what the experimental chemical shift should be. You will do this in the next exercise and compare to the literature reference.

3.4 Questions¶

1. What is the relation between chemical shielding and chemical shift?
2. What is the purpose of a reference in the analysis of experimental data?
3. Can you compare the chemical shielding of a single compound to experimental values?

4 Predicting experiment from theory

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• Calculate the chemical shielding for two O-containing compounds
• Predict their chemical shifts and compare to experimental chemical shifts
• Compare your series against a larger literature series

4.1 Task¶

Calculate the chemical shieldings for primitive cell BaZrO₃ and CaO. Compare that to the predicted chemical shift of the linear fit based on the BaSnO₃, MgO, SrO, and SrTiO₃ series.

As has been shown in the previous exercise based on BaSnO₃, MgO, SrO, and SrTiO₃, there is a linear relationship between experimental chemical shift and the calculated chemical shieldings. In Example 3, O chemical shieldings were fitted against shifts. BaSnO₃, MgO, SrO, and SrTiO₃ correspond to a subset of the Laskowski series [J. Chem. Phys. 146 (2017) 064115]. In the Laskowski paper, they also compared the chemical shieldings for fluorine to experiment and validated the plane-augmented-wave (PAW) approach in comparison to an all-electron (AE) basis, obtaining fits for the following linear equation in each case:

$$\tag{10} \begin{equation} \delta_{\mathrm{iso}}^{\mathrm{exp}}[N] = \delta_{\mathrm{ref}}[N] + m \sigma_{\mathrm{iso}}^{\mathrm{calc}}[N]. \end{equation} $$

From this fitted equation, it is possible to predict the experimental chemical shift from a calculated chemical shielding.

In this exercise, you will calculate the chemical shieldings for two O-containing compounds, BaZrO₃ and CaO. You will then predict their chemical shifts from this linear fit based on BaSnO₃, MgO, SrO, and SrTiO₃ and compare them to their experimental values.

4.2 Input¶

The input files to run this example are prepared at $TUTORIALS/NMR/e04_comparison.

VASP looks in the current directory for four main input files, i.e., POSCAR, INCAR, KPOINTS, and POTCAR. Check them out!

BaZrO3/POSCAR

BaZrO3
   4.19728855703300
   1.00000000 0.00000000 0.00000000
   0.00000000 1.00000000 0.00000000
   0.00000000 0.00000000 1.00000000
   Ba Zr O
   1 1 3
Direct
0.50000000 0.50000000 0.50000000
0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.50000000
0.00000000 0.50000000 0.00000000
0.50000000 0.00000000 0.00000000

CaO/POSCAR

CaO
   4.80709784405400
   0.50000000 0.50000000 0.00000000
   0.00000000 0.50000000 0.50000000
   0.50000000 0.00000000 0.50000000
   Ca O
   1 1
Direct
0.00000000 0.00000000 0.00000000
0.50000000 0.50000000 0.50000000

INCAR

ENCUT = 500              # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV

EDIFF = 1E-8             # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate          # Sets the "precision" mode
LASPH = .TRUE.           # Non-spherical contributions to the gradient of the density in the PAW spheres

LCHIMAG = .TRUE.         # Turns on linear response for chemical shifts
LNMR_SYM_RED = .TRUE.    # Consistent symmetry with star and k-space derivatives
NLSPLINE = .TRUE.        # Differentiable projectors in reciprocal space

KPOINTS

C
0
G
 6 6 6 
 0 0 0

BaZrO3/POTCAR

 Pseudopotential of Ba, Zr, and O.

CaO/POTCAR

 Pseudopotential of Ca and O.

experiment.dat

BaSnO3 143
MgO 47
SrO 390
SrTiO3 465
BaZrO3 376
CaO 294

The INCAR files are identical to those used in Example 1.

The KPOINTS file specifies a Γ-centered $\mathbf{k}$-mesh. An 6x6x6 $\mathbf{k}$-mesh is used.

4.3 Calculation¶

Open a terminal and navigate to this example's directory:

cd $TUTORIALS/NMR/e04_*

Take a look in the directory with ls and you will see two directories: BaZrO3 and CaO. In these two directories, there are the corresponding INCAR, KPOINTS, POSCAR, and POTCAR files that will calculate the chemical shieldings.

From the central e04_* directory, enter CaO and calculate the chemical shielding by entering the following into the terminal:

cd $TUTORIALS/NMR/e04_*/CaO
mpirun -np 4 vasp_std

BaZrO3 can be very time-consuming, so we recommend only running CaO.

Optional: If you have enough time and want to try running both calculations, you can use the instructions below:

cd $TUTORIALS/NMR/e04_*
bash compare.sh

rm -f O_shielding.dat
cat ../e03_*/O_shielding.dat >> O_shielding.dat
for a in *O*/
do
cd $a
mpirun -np 4 vasp_std
if [[ $a ==  *"O3"* ]]; then
shielding=$(grep -A 15 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
else
shielding=$(grep -A 10 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
fi
cd ..
echo ${a///} $shielding >> O_shielding.dat
done

Before extracting the shieldings for CaO, make sure to put the calculated chemical shielding from the other compounds using the following:

cd $TUTORIALS/NMR/e04_*
echo "BaSnO3 74.0956" >> O_shielding.dat
echo "MgO 194.8040" >> O_shielding.dat
echo "SrO -219.9072" >> O_shielding.dat
echo "SrTiO3 -281.9145" >> O_shielding.dat
echo "BaZrO3 -161.6974" >> O_shielding.dat

Extract the chemical shielding for the O in CaO by entering the following into the terminal:

cd $TUTORIALS/NMR/e04_*/CaO
shielding=$(grep -A 10 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
echo "CaO" $shielding >> ../O_shielding.dat

In Example 3, you already obtained a linear fit connecting experiment and calculation. You have just calculated the shieldings for two new O-containing compounds: BaZrO₃ and CaO. Now, take the calculated chemical shifts and predict what you expect their experimental chemical shieldings to be by entering them into your linear fit:

The difference in chemical shifts between predicted and experimental chemical shifts is larger than for the compounds that the fit was based on.

What do you think might have caused this?

Let's plot them to visualize this better:

For BaZrO₃, this deviation is due to unfrozen 5s and 5p semi-core states, i.e. the shallow core states are allowed to be polarized in the crystal field. CaO is further from the line due to a known effect, the empty Ca 3d states are lying too close to the valence band maximum in DFT. The closeness of the 3d states results in the O shift deviation [J. Chem. Phys. 146 (2017) 064115].

The y-intercept $\sigma_{ref}$, gradient $m$, and Pearson coefficient $r$ can be obtained from your fit and allow you to assess its quality. Run the following script to determine these for your fit and compare to the literature values:

The y-intercept $σ_{ref}$, gradient $m$, and Pearson coefficient $r$ are presented in a table above and allow you to compare the quality of your fit.

The Laskowski O-series contains 10 compounds. The "standard" set uses a larger energy cutoff of 700 eV and harder "_h" POTCARs (smaller core radii) than in our calculations. The "optimal" set used an energy cutoff of 900 eV and specially improved POTCARs. The results from both of these series are, perhaps unsurprisingly, an improvement on our quick calculations. Your 4-series contains a "well-behaved" subset of the Laskowski series, while the 6-series introduces some additional complexity.

It is worth noting that, ideally, $m$ would be -1, with calculation and experiment differing by only the reference chemical shift (i.e. $σ_{ref}$). This deviation of $m$ from -1 is due to limitations of DFT, e.g. to properly capture exchange-correlation effects.

When calculating chemical shieldings, not only the convergence parameters must be considered but also whether suitable PAW pseudopotentials have been selected in the first place. In most cases, the ordinary pseudopotentials are sufficient. In a few cases, like for BaZrO₃, harder pseudopotentials, with smaller core radii are required.

4.4 Questions¶

1. What impact does adding BaZrO₃ and CaO do to your series?
2. Does the choice of PAW pseudopotential affect the chemical shielding?
3. Why does the gradient $m$ deviate from the ideal -1?

Well done! You have finished Part 1 of the NMR tutorial!¶

Go to Top $\uparrow$