Calculating the chemical shieldings
It is possible to calculate the chemical shielding, the magnetic susceptibility, and the electric field gradient. The theory is covered in the NMR category page and corresponding pages.
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 [1][2]. The chemical shielding is directly related to the chemical shift δ recorded in nuclear magnetic resonance (NMR), 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. [3].
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 [4][5].
| 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. |
Magnetic susceptibility
The magnetic susceptibility is the degree of magnetization of a material in response to an applied magnetic field. It is a bulk property, in contrast to the chemical shielding, which is for each nucleus. Like the chemical shielding, the magnetic susceptibility is calculated by linear response using LCHIMAG [1][2], so they will both be shown in the same OUTCAR file. The magnetic susceptibility is measured using the Guoy balance (or method); alternatively, an Evans or Faraday balance can be used.
Step-by-step instructions
The magnetic susceptibility is calculated post-self-consistent field (post-SCF) using LCHIMAG. A well-converged SCF calculation is therefore crucial. The magnetic susceptibility can be sensitive to several input parameters that must all be independently tested.
Step 1 (optional): Calculate the magnetic susceptibility using a previously converged calculation
Since the magnetic susceptibility 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 magnetic susceptibility. A setting of EDIFF = 1E-8 eV is generally recommended. Convergence is taken to be within 0.01 (dimensionless units).
Step 3: Converge the plane-wave basis
A larger than standard plane-wave energy cutoff is required to fully converge the magnetic susceptibility. 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.01 (dimensionless units). The magnetic susceptibility is less dependent on the energy cutoff than the chemical shielding is.
Step 4: Converge the k point mesh
Similar to the basis, the k point mesh can strongly influence the magnetic susceptibility. The k point mesh should be increased incrementally, i.e., 1x1x1, 2x2x2, 3x3x3, until convergence within 0.01 (dimensionless units) is achieved. It is slightly more dependent on k point mesh than the chemical shieldings are.
Step 5: Compare to experiment
The purpose of these calculations is to compare to the experiment. The computed magnetic susceptibilities can be directly compared to the measured magnetic susceptibility, in contrast to the chemical shielding [6].
Recommendations and advice
Calculating the magnetic susceptibilities 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. Since the same tag is used, much of the advice for chemical shieldings is applicable for the magnetic susceptibility.
PAW pseudopotentials
The standard PAW pseudopotentials POTCAR used are sufficient for calculating the magnetic susceptibility. Small differences on the order of 0.1 (dimensionless units) are seen when using slightly different types of POTCAR, e.g., GW (*_GW).
Additional tags
To ensure tight precision, the precision should be set to PREC = Accurate, rather than Normal. There is one additional tag, ICHIBARE that can be used, though the default is usually sufficient and increases the computational load significantly.
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 [7]. By including the quadrupole moment of the isotopes, the quadrupole coupling constants Cq 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 , up to 50 % [7]; 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. [8].
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 Cq are.
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.
Mind: Be aware of some specifics relevant to the implementation used:
|
Structure
The electric field gradient can be extremely dependent on structure, to the extent that using the 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 TiO2 rutile, a shift in position from 0.305 in internal coordinates to 0.3025 made a difference of 50 % to for the Ti [7]. 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.
PAW pseudopotentials
The standard PAW pseudopotentials POTCAR used are sufficient for calculating the electric field gradient. Using GW pseudopotentials can significantly improve results. Semi-core electrons can be important, so POSCAR files with *_pv or *_sv can improve the results, as will the explicit inclusion of augmentation channels with -projectors.
Additional tags
To ensure tight precision, the precision should be set to PREC = Accurate, rather than Normal.
The LASPH should be set to .TRUE., turning on the non-spherical contributions to the gradient of the density inside the PAW spheres.
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
References
- ↑ a b C. J. Pickard and F. Mauri, All-electron magnetic response with pseudopotentials: NMR chemical shifts, Phys. Rev. B 63, 245101 (2001).
- ↑ a b J. R. Yates, C. J. Pickard, and F. Mauri, Calculation of NMR chemical shifts for extended systems using ultrasoft pseudopotentials, Phys. Rev. B 76, 024401 (2007).
- ↑ G. A. de Wijs, R. Laskowski, P. Blaha, R. W. A. Havenith, G. Kresse, and M. Marsman, NMR shieldings from density functional perturbation theory: GIPAW versus all-electron calculations, J. Chem. Phys. 146, 064115 (2017).
- ↑ F. Vasconcelos, G.A. de Wijs, R. W. A. Havenith, M. Marsman, and G. Kresse, Finite-field implementation of NMR chemical shieldings for molecules: Direct and converse gauge-including projector-augmented-wave methods, J. Chem. Phys. 139, 014109 (2013).
- ↑ G.A. de Wijs, G. Kresse, R. W. A. Havenith, and M. Marsman, Comparing GIPAW with numerically exact chemical shieldings: The role of two-center contributions to the induced current, J. Chem. Phys. 155, 234101 (2021).
- ↑ F. Mauri, S. G. Louie, Magnetic Susceptibility of Insulators from First Principles, Phys. Rev. Lett. 76, 4246 (1996).
- ↑ a b c H. M. Petrilli, P. E. Blöchl, P. Blaha, and K. Schwarz, Electric-field-gradient calculations using the projector augmented wave method, Phys. Rev. B 57, 14690 (1998).
- ↑ P. Pyykkö, Year-2017 nuclear quadrupole moments, Mol. Phys. 116, 1328-1338 (2018).