Calculating the chemical shieldings

There are several different options available to calculate NMR properties. It is possible to calculate the chemical shielding, the two-center contributions, and the electric field gradient. The theory is already covered in the NMR category page and corresponding pages, so it will not be reiterated here.

For all calculations, tighter convergence settings than typical are required, e.g. for a structure relaxation. No additional files are required besides the four standard POSCAR, POTCAR, INCAR, and KPOINTS, unless specifically mentioned. It is important to have a well-converged structure, as all of these calculations described below can be very sensitive to it. For each of the following calculations, the NMR property is calculated post-SCF.

Chemical shielding

The chemical shielding tensor σ is the relation between the induced and external magnetic fields and describes how much the electrons shield the nuclei from an external field. The absolute chemical shielding is calculated by linear response using LCHIMAG . The chemical shielding is directly related to the chemical shift δ (cf. NMR category page and LCHIMAG page for details) and, indirectly, to the resonance frequency.

Step-by-step instructions

The chemical shielding is calculated post-self-consistent field (post-SCF) using LCHIMAG. A well-converged SCF calculation is therefore crucial. The chemical shielding is very sensitive to several input parameters that must all be independently tested.

Step 1 (optional): Calculate the chemical shielding using a previously converged calculation

Since the chemical shielding is calculated post-SCF, you can use a previously converged WAVECAR with ISTART = 1 and NELM = 1. The corresponding density, CHGCAR is calculated from the WAVECAR file before the first elementary step so it need not be included.

Step 2 (optional): Determine a suitable energetic break value

The break condition for the self-consistency step EDIFF strongly influences the chemical shielding. A setting of EDIFF = 1E-8 eV is generally recommended. Convergence is taken to be within 0.1 ppm.

Step 3: Converge the plane-wave basis

A large plane-wave energy cutoff is required to fully converge the chemical shieldings. Perform multiple calculations while increasing the basis set size, as defined in ENCUT, incrementally (e.g., by 100 eV intervals). Convergence should be aimed to be within 0.1 ppm, although this will not be feasible for heavier elements.

Step 4: Converge the k point mesh

Similar to the basis, the k point mesh can strongly influence the chemical shielding. The k point mesh should be increased incrementally, i.e., 1x1x1, 2x2x2, 3x3x3, until convergence within 0.1 ppm is achieved. It is only necessary to converge the k point mesh for crystals, gas-phase molecules should use the Γ-point only.

Step 5: Compare to experiment

The purpose of these calculations is to compare to experiment. However, the calculated absolute chemical shieldings are not directly comparable to the measured chemical shift due to the lack of a reference. To avoid bias from any single calculation, a series of calculated and their corresponding experimental values are used. The experimental chemical shifts are plotted against the calculated chemical shieldings as is found in Fig. 3 of Ref. .

Recommendations and advice

Calculating the chemical shielding requires tightly converged settings. As described in the step-wise introduction above, converging with respect to EDIFF, ENCUT, and the k point mesh is very important. There are a few additional settings that should be considered.

PAW pseudopotentials

The standard PAW pseudopotentials POTCAR used are sufficient for calculating the chemical shielding. The GIPAW is applied using the projector functions and partial waves that are stored in the regular POTCAR files. The completeness of these projector functions and partial waves determines the quality of the results. Using slightly different types of POTCAR, e.g., GW (*_GW) or with additional valence (*_sv, *_pv), can change the calculated shielding by a few ppm for the first and second row sp-bonded elements (except for H).

The PAW reconstruction with all-electron partial waves is crucial for calculating the field on the nucleus. It is therefore important to use a consistent exchange-correlation functional and so LEXCH in the POTCAR should not be overwritten with an explicit GGA tag in the INCAR if possible.

Insufficient memory

For calculating the chemical shieldings, speed had been favored over saving memory, resulting in insufficient memory occasionally. Since the linear response calculation is parallel over k points, this can be used to economize on memory by performing a regular SCF calculation at high accuracy on the full k point mesh and saving the CHGCAR file. Using ICHARG = 11 start a chemical shielding calculation for each individual k point in the first Brillouin zone (IBZ) separately, starting from CHGCAR. The shieldings can then be calculated as a k point weighted average of the symmetrized shieldings of the individual k points.

Additional tags

To ensure tight precision, the precision should be set to PREC = Accurate, rather than Normal.

Several additional INCAR tags should be considered. Specifically, LASPH should be set to .TRUE., turning on the non-spherical contributions to the gradient of the density inside the PAW spheres. Occasionally, e.g. for systems containing H or first-row elements, and short bonds, the two-center contributions to the augmentation currents in the PAW spheres are important. In this case, LLRAUG = .TRUE. should be used .

Important: The treatment of the orbital magnetism is non-relativistic. This is suitable for light nuclei. The standard POTCARs are scalar-relativistic and account partially for relativistic effects. The accuracy can be improved using LBONE, which restores the small B-component of the wave function inside the PAW spheres. Spin-orbit coupling is not implemented for chemical shift calculations.

Electric field gradient

Nuclei with a spin > ± ½ are called quadrupolar nuclei. They have a non-spherical shape and therefore a non-zero electric field gradient (EFG) at the nucleus. The EFG is calculated using LEFG . By including the quadrupole moment of the isotopes, the quadrupole coupling constants C_q can be calculated (multiple definitions exist in the literature, ensure that you are correctly comparing). These are measured using nuclear quadrupole resonance (NQR) spectroscopy, a type of zero- to ultralow-field (ZULF) NMR.

Step-by-step instructions

The electric field gradient is calculated post-self-consistent field (post-SCF) using LEFG. A well-converged SCF calculation is therefore crucial. The electric field gradient is very sensitive to several input parameters that must all be independently tested. In particular, small differences in the structure can make big differences to $V_{zz}$, up to 50 % ; see the Advice section for more details. Make sure to have a well-optimized structure before you begin the convergence tests.

Step 1 (optional): Calculate the chemical shielding using a previously converged calculation

Since the chemical shielding is calculated post-SCF, you can use a previously converged WAVECAR with ISTART = 1 and NELM = 1. The corresponding density, CHGCAR is calculated from the WAVECAR file before the first elementary step so it need not be included.

Step 2a: Define the nuclear quadrupolar moments

The calculated electric field gradients are not observable in experiment. Instead, the quadrupolar coupling constant can be calculated so long as the nuclear quadrupolar moments are defined in QUAD_EFG. Each species in your POSCAR file should be defined; there is no need to define each individual ion. A short table of values can be found in Ref. .

Step 2b (optional): Determine a suitable energetic break value

The break condition for the self-consistency step EDIFF strongly influences the chemical shielding. A setting of EDIFF = 1E-8 eV is generally recommended. Convergence is taken to be within 0.1 ppm.

Step 3: Converge the plane-wave energy cutoff

A large plane-wave energy cutoff is required to fully converge the electric field gradient. Perform multiple calculations while increasing the basis set size, as defined in ENCUT, incrementally (e.g., by 100 eV intervals). Convergence should be aimed to be within 3 significant figures, although this will not be feasible for heavier elements.

Step 4: Converge the k point mesh

Similar to the basis, the k point mesh can strongly influence the coupling constant. The k point mesh should be increased incrementally, i.e., 1x1x1, 2x2x2, 3x3x3, until convergence within 3 significant figures is achieved. It is only necessary to converge the k point mesh for crystals, gas-phase molecules should use the Γ-point only.

Step 5: Compare to experiment The purpose of these calculations is to compare directly to experiment. The EFG that has been calculated is not directly measurable but the quadrupolar coupling constants C_q are.

Input

A typical INCAR file is given below:

 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

 LEFG = .TRUE.            # Electric field gradient calculations
 QUAD_EFG = 0. -696. 20.44 0. 2.860  # Nuclear quadrupolar moments for Pb I N O D

Important: Make sure to replace the QUAD_EFG in the INCAR with the values for the isotopes in your system.

Output

The EFG is listed atom-wise after the SCF cycle has been completed. First, the full 3x3 tensor is printed:

  Electric field gradients (V/A^2)
 ---------------------------------------------------------------------
  ion       V_xx      V_yy      V_zz      V_xy      V_xz      V_yz
 ---------------------------------------------------------------------
    1        -         -         -         -         -         -       

The tensor is then diagonalized and reprinted:

  Electric field gradients after diagonalization (V/A^2)
  (convention: |V_zz| > |V_xx| > |V_yy|)
 ----------------------------------------------------------------------
  ion       V_xx      V_yy      V_zz     asymmetry (V_yy - V_xx)/ V_zz
 ----------------------------------------------------------------------
    1       -         -         -             -         

The corresponding eigenvectors are printed atom-wise. Finally, the quadrupolar parameters are presented, which, unlike the EFG, may be measured by experiment.

            NMR quadrupolar parameters

  Cq : quadrupolar parameter    Cq=e*Q*V_zz/h
  eta: asymmetry parameters     (V_yy - V_xx)/ V_zz
  Q  : nuclear electric quadrupole moment in mb (millibarn)
 ----------------------------------------------------------------------
  ion       Cq(MHz)       eta       Q (mb)
 ----------------------------------------------------------------------
    1        -             -         -                      

Recommendations and advice

Calculating the electric field gradient requires tightly converged settings. As described in the step-wise introduction above, converging with respect to EDIFF, ENCUT, and the k point mesh is very important. There are a few additional settings that should be considered.

Structure

The electric field gradient is extremely dependent on structure, to the extent that experimental structure can improve results. A small difference in the positions of atoms can make a huge difference to the EFG. For the O in TiO₂ rutile, a change from 0.305 in internal coordinates to 0.3025 made a difference of 50 % to $V_{zz}$ for the Ti . This is an atypical case but highlights the importance of using a well-optimized structure, ideally the experimental structure if available. This extreme sensitivity to the structure is indicative of why the quadrupolar coupling constant is so useful for determining information about a system's chemical environment.

The same tight settings for chemical shielding are required, alongside a stronger dependence on the structure and the chosen POTCAR used:

• A larger ENCUT value than usual, generally much higher than the value given by ENMAX in the POTCAR file, e.g. 800 eV for C.
• A small EDIFF is typically required to provide converged chemical shifts, e.g. 1E-8 eV.
• Tighter precision, e.g. PREC = Accurate.
• The structure is extremely important so using the experimental structure can improve results. EFG is very sensitive to structure so differences of 2.5 mÅ can make differences of 40 % to $V_{zz}$ .
• The use of PAW potentials has a strong influence, GW POTCAR files often improve values.
• Semi-core electrons can be important (check the POSCAR files with *_pv or *_sv) as well as explicit inclusion of augmentation channels with $d$-projectors.

In addition, test whether non-spherical contributions are important for your system:

• Non-spherical contributions to the gradient of the density inside PAW spheres, i.e. LASPH = .TRUE.

Be aware of some specifics relevant to the implementation used:

• Several definitions of $C_q$ are used in the NMR community, ensure that you are comparing between the same definitions in calculation and experiment.
• For heavy nuclei inaccuracies are to be expected because of an incomplete treatment of relativistic effects.

For each system, make sure to test that the chemical shieldings calculated are converged with respect to ENCUT, EDIFF, and KPOINTS mesh. Convergence is considered to be typically within 3 significant figures for the EFG in the z-direction $V_{zz}$ (though this may be infeasible for heavier elements).

Example scripts for convergence tests

Several tests are necessary to obtain various NMR parameters. Make sure to change the example INCAR files to include the tags for your desired calculation. We provide some example scripts below:

Energetic break criterion tests

For converging the energetic break criterion for a single ionic step (EDIFF), start with the 1E-4 and then increase by orders of magnitude:

Energetic break criterion: INCAR.nmr

 PREC = Accurate  
 ENCUT = 400         
 EDIFF = 1E-4  
 ISMEAR = 0; SIGMA = 0.01 

Script to loop through EDIFF from 1E-4 eV to 1E-8 eV:

for a in 4 5 6 7 8
do
cp INCAR.nmr INCAR
sed -i "s/1E-4/1E-$a/g" INCAR

mpirun -np 4 $PATH_TO_EXECUTABLE/vasp_std

cp OUTCAR OUTCAR.$a
done

k-points tests

For converging k points, start with the Γ-point and increase the k-point mesh incrementally:

Initial Γ-only mesh: KPOINTS.nmr

C
0
G
 1 1 1
 0 0 0

Script to go through k-point meshes from Γ-only to 8x8x8:

for a in 1 2 4 6 8
do
cp KPOINTS.nmr KPOINTS
sed -i "s/1 1 1/$a $a $a/g" KPOINTS

mpirun -np 4 $PATH_TO_EXECUTABLE/vasp_std

cp OUTCAR OUTCAR.$a
done

Energy cutoff tests

For converging the energy cutoff, start with the value of ENMAX given in the POTCAR file and then increase incrementally in steps of 100 eV:

Initial INCAR: INCAR.nmr

 PREC = Accurate  
 ENCUT = 400         
 EDIFF = 1E-6  
 ISMEAR = 0; SIGMA = 0.01 

Script to loop through ENCUT from 400 eV to 600 eV:

for a in 400 500 600 700 800
do
cp INCAR.nmr INCAR
sed -i "s/400/$a/g" INCAR

mpirun -np 4 $PATH_TO_EXECUTABLE/vasp_std

cp OUTCAR OUTCAR.$a
done