CASTEP > Theory in CASTEP > Calculating NMR shielding tensors using pseudopotentials

# Calculating NMR shielding tensors using pseudopotentials

A uniform external magnetic field, B, applied to a sample induces an electric current. In an insulating non-magnetic material, only the orbital motion of the electrons contributes to the this current. In addition, for the field strengths typically used in NMR experiments, the induced electric current is proportional to the external field, B. This first-order induced current, j(1)(r), produces a non-uniform magnetic field:

Eq. CASTEP 42

The shielding tensor, σ(r), connects the induced magnetic field to the applied magnetic field:

Eq. CASTEP 43

and the isotropic shielding is given by:

Eq. CASTEP 44

## All-electron magnetic response with pseudopotentials

Information about σ(r) at nuclear positions can be obtained from NMR experiments. In CASTEP, the first-order induced current, j(1)(r), is calculated and then used to evaluate the induced magnetic field, B(1) via Eq. 42. Calculations are performed with the gauge-including projector augmented-wave method (GIPAW) developed by Pickard and Mauri (2001). The projector augmented-wave formalism originally introduced by Blöchl (1994) makes it possible to obtain expectation values of all-electron operators in terms of pseudo-wavefunctions coming from a pseudopotential calculation and the GIPAW formalism reconciles the requirement of translational invariance of the crystal in a uniform magnetic field with the localized nature of Blöchl's projectors. The theory was initially applied with Troullier-Martins norm-conserving pseudopotentials (Pickard and Mauri, 2001) and has recently been extended to the case of ultrasoft pseudopotentials (Yates, 2003).

The first-order perturbation wavefunction, Ψ(1):

Eq. CASTEP 45

is proportional to the magnetic perturbation, H(1), which, in pseudopotential calculations, contains a contribution due to the non-locality of the pseudopotential, in addition to the usual term:

Eq. CASTEP 46

## Calculation of perturbed wavefunctions

In order to obtain the first-order perturbation wavefunction Ψ(1), the application of the Green function in Eq. 46 is performed in CASTEP by solving a system of linear equations:

Eq. CASTEP 47

where Q is the projection operator onto the unoccupied Hilbert subspace. The system is solved by a conjugate gradient minimization scheme, whereby the tedious summation over unoccupied states is avoided. The method is, thus, comparable in computational cost to the calculation of the ground-state electronic structure.

## NMR properties of molecules and crystals

The calculation of the induced current involves evaluation of expectation values of the position operator, whose definition for extended systems is not straightforward. Two approaches have been developed (Pickard and Mauri, 2001) and implemented in CASTEP. For translationally invariant systems with well-defined Blöch vectors, k, the application of the position operator can be reduced to calculating the gradient in k-space. Alternatively, for finite systems, which are represented by supercell structures (a molecule in a box) in CASTEP, one can apply a real-space representation of the position operator. The artificial periodicity imposed on the system causes the position operator to acquire a sawtooth shape. Thus, in the vicinity of the box boundaries, the position operator is distorted. Because of this, a sufficiently large supercell must be constructed so that the electronic density does not reach the regions where the position operator behaves unphysically. This condition is easily satisfied by the usual physical requirement that the molecules in neighboring boxes do not interact.

## NMR shielding tensor

Finally, the induced magnetic field, B(1), is calculated from the induced current and its value at the nuclear position determines the NMR shielding tensor. The calculated tensor is rather sensitive to numerical errors in both the unperturbed wavefunction, Ψ(0), and its first-order perturbation, Ψ(1), so the convergence requirements for NMR calculations are more stringent than, for example, routine band structure calculations. To obtain reliable results, the convergence of the wavefunction with respect to the number of basis functions (plane waves) must be ensured and a well-converged conjugate gradient procedure put in place to determine the coefficients of the plane waves. Apart from that, in calculating solid-state NMR properties (a small unit cell as opposed to a molecule in a box) the convergence of the k-space sum over the Brillouin zone must be achieved. The NMR observables are not particularly sensitive to uncertainties in the ground-state density distribution in the system, which determines the unperturbed Hamiltonian. Here, the usual convergence criteria for self-consistent DFT calculations are sufficient.

## Sensitivity of NMR to calculation parameters

NMR properties calculations are very sensitive to the number of k-points in the Brillouin zone. It is hard to construct a general rule for the number of k-points required, so it is recommended that you carry out convergence tests. An example convergence test for gallium phosphide is given below.

Table 1. Convergence test for GaP

Iso Aniso Cq
P 2 × 2 × 2 -1022.76 -6986.64 2.487e+02
4 × 4 × 4 217.24 -1300.93 4.045e+01
8 × 8 × 8 315.17 -165.40 1.796e+00
12 × 12 × 12 311.75 -23.34 1.010e-01
16 × 16 × 16 311.01 -3.36 -8.199e-03
Ga 2 × 2 × 2 943.83 -3323.73 1.913e+01
4 × 4 × 4 1488.73 -661.68 2.900e+00
8 × 8 × 8 1512.39 -68.91 1.240e-01
12 × 12 × 12 1509.21 -8.47 1.840e-02
16 × 16 × 16 1508.57 -1.13 1.256e-02

It follows from Table 1 that a grid of at least 12 × 12 × 12 is required to achieve convergence within 1-2 ppm, and finer grids might be needed to achieve higher accuracy.

NMR properties are very sensitive to atomic positions, which makes NMR such a useful experimental tool for structure analysis. This implies that it is highly recommended to perform a geometry optimization run prior to the NMR calculation. If the structure of a crystal or molecule under investigation is obtained from an experimental study, it is advisable to optimize positions of hydrogen atoms first of all. If the forces on heavy atoms as reported by CASTEP are high (in excess of 1 eV/Å) then complete structure optimization is recommended.