Calculating the hyperfine coupling constant

The hyperfine coupling constant (cf. hyperfine splitting) describes the interaction between the nuclear magnetic dipole moment and the magnetic field generated by the electrons (i.e., the nuclear spin-electron spin coupling). The hyperfine coupling constant is calculated using LHYPERFINE . The hyperfine splitting often includes the interaction between the nuclear quadrupole moment and the electric field gradient (EFG), which is calculated separately using LEFG and the description in performing an EFG calculation. The hyperfine splitting can be measured using electron paramagnetic resonance (EPR), also called electron-spin resonance (ESR), and in atomic spectroscopy. The theory is covered in the NMR category page and LHYPERFINE page.

Step-by-step instructions

The hyperfine constant is calculated post-self-consistent field (post-SCF) using LHYPERFINE. A well-converged SCF calculation is therefore crucial. The hyperfine coupling constant is sensitive to several input parameters that must all be tested.

Step 1 (optional): Calculate the hyperfine constant using a previously converged calculation

Since the hyperfine constant 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 need not be provided.

Step 2a: Define the nuclear gyromagnetic ratios

The hyperfine coupling constant depends on the nuclear gyromagnetic ratios defined in NGYROMAG. Since the defaults are set to 1, the gyromagnetic ratios must be defined to obtain meaningful coupling constants. 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. . For a more complete list, see the NGYROMAG tag documentation.

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

The break condition for the self-consistency step EDIFF does not strongly influence the coupling parameter for our test systems. However, it is important to confirm this for your system before performing more expensive convergence tests.

The hyperfine coupling constant depends on the nuclear gyromagnetic ratios defined in NGYROMAG. Since the defaults are set to 1, the gyromagnetic ratios must be defined to obtain meaningful coupling constants. Each species in your POSCAR file should be defined; there is no need to define each individual ion.

Use the converged energetic break condition for all subsequent convergence tests.

Step 3: Converge the plane-wave basis

The plane-wave basis can strongly influence the coupling constant. Unconverged values should not be compared to experiment. 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 MHz, although this will not be feasible for heavier elements.

Use the converged plane-wave energy cut-off for all subsequent convergence tests.

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 to within 0.1 MHz is achieved.

Use the converged k--point mesh for all subsequent calculations.

Step 5: Compare to experiment

The purpose of these calculations is to compare to experiment. An example is given in Ref. . It is important to include core contributions , as these can account for a significant portion of the Fermi contact term. The total coupling parameter can be compared to EPR.

Step 6 (optional): Perform hybrid calculations

In the literature, HSE06 has been shown to better localize defect states, which improves comparison to experiment relative to PBE . Consider performing a hybrid calculation, if it is affordable. You will need to repeat the convergence tests described in Steps 2-4.

Mind: Significant time savings can be achieved for hybrid functionals by starting from a WAVECAR from a GGA calculation.

Recommendations and advice

The hyperfine coupling constant requires tightly converged settings. The energetic break condition EDIFF and the plane-wave energy cutoff ENCUT impact the convergence of the hyperfine coupling constant. For solid-state systems, the choice of k-point mesh KPOINTS used can also be very important. Besides these input settings, the hyperfine coupling constant is influenced by several other factors, specifically structure, POTCAR, and method.

Structure

The structure defined in POSCAR will impact the hyperfine constant in two ways. The first and most important is that cells that are too small may converge to non-magnetic systems. For example, the NV-diamond defect cannot be properly described by a 15-atom supercell (based on a 2x2x2 cell from primitive diamond). As you increase the k-point mesh, the magnetization disappears due to coupling between neighboring defect sites (i.e., increasing the k-point mesh causes the coupling to disappear (all zeros)). Be careful to use a large enough cell for your calculation, otherwise even converged settings will produce meaningless outputs. Sometimes, a non-magnetic solution is incorrectly found. If you are certain that it should be magnetic, then you can force this by using NUPDOWN to fix the number of unpaired electrons in your system during the calculation. Ensure that the energy of the magnetic state is lower than the non-magnetic state.

A second structural problem will come from the precise POSCAR that you use. Slightly different lattice parameters (10 mÅ) can change the hyperfine coupling constant by ~0.5 MHz. Make sure to use a well-optimized structure.

Hybrid functionals

The calculated hyperfine coupling parameter is strongly influenced by the chosen method. PBE tends to underestimate the coupling constant relative to experiment. HSE06 and other hybrid functionals improve this, matching well with experiment . Hybrid functional localize defect states, resulting in an improved description over GGA functionals. The INCAR tags are specified in hybrid functionals and no additional tags are required. We found that the improvement for hybrid functionals is seen for both range-separated and unscreened hybrid functionals.

PAW pseudopotential

The choice of PAW pseudopotential can be important. GW pseudopotentials (i.e. *_GW) make a difference of a few ppm to the chemical shielding.

It is expected that using GW pseudopotentials will improve the description of the hyperfine coupling parameter for heavier elements, though it has not been tested. Additionally, for heavier elements, it is expected that using additional valence electrons (i.e., *_sv, *_pv) will improve calculating the hyperfine coupling constant.

Additional tags

There are a few additional INCAR tags that should be used to ensure tight convergence, specifically, the precision should be set to PREC = Accurate, rather than Normal. We also recommend EDIFF to a minimum of 1E-6 eV to ensure well-converged WAVECAR.

LASPH = .TRUE. includes the non-spherical contributions to the gradient of the density in PAW spheres. For symmetric systems, the difference is small, a few ppm. For asymmetric systems, the difference is much larger, on the order of 10 ppm. Therefore, the non-spherical contributions must be included. However, using LASPH = .TRUE. decreases the calculated hyperfine coupling constant, weakening agreement with experiment, including in the literature where it was not previously used .

Example scripts for convergence tests

Several tests are necessary to obtain converged coupling parameters. We provide some example scripts below:

Important: Change the values of NGYROMAG to your corresponding system.

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.hyperfine

PREC = Accurate  
ENCUT = 400         
EDIFF = 1E-6  
ISMEAR = 0; SIGMA = 0.01 
LHYPERFINE = .TRUE.
NGYROMAG = 10.7084 3.077
LASPH = .TRUE.
ISPIN = 2

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

rm -f hyperfine.dat
for a in 4 5 6 7 8
do
cp INCAR.hyperfine 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.hyperfine

C
0
G
 1 1 1
 0 0 0

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

rm -f hyperfine.dat
for a in 1 2 4 6 8
do
cp KPOINTS.hyperfine 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 from at least the value of ENMAX given in the POTCAR file and then increase incrementally in steps of 100 eV:

Initial INCAR: INCAR.hyperfine

PREC = Accurate  
ENCUT = 400         
EDIFF = 1E-6  
ISMEAR = 0; SIGMA = 0.01 
LHYPERFINE = .TRUE.
NGYROMAG = 10.7084 3.077
LASPH = .TRUE.
ISPIN = 2

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

rm -f hyperfine.dat
for a in 400 500 600
do
cp INCAR.hyperfine INCAR
sed -i "s/400/$a/g" INCAR

mpirun -np 4 $PATH_TO_EXECUTABLE/vasp_std

cp OUTCAR OUTCAR.$a
done