Nuclear Magnetic Resonance

Theory

Ref: [97]

A magnetic field $\mathbf{B}$ is applied to the system, which reacts by induced electronic ring currents. These currents produce an additional magnetic field by themselves, which is not homogeneous in space. Therefore, the actual magnetic field at the ionic positions is different for all atoms in the cell. This field determines the resonance frequency of the nuclear spin, and this resonance can be measured with a very high accuracy.

The perturbation Hamiltonian is given by

$$H^{(1)} = \frac 12 \frac em \mathbf{p}\times \mathbf{r}\cdot \mathbf{B}.$$ (23)

The difficult part of this Hamiltonian lies in the position operator which is ill defined in a periodic system. To get around this, the wavefunctions are localized and for each localized orbital, Eq. (23) is applied individually assuming the orbital being isolated in space. Around each orbital, a virtual cell is placed such that the wavefunction vanishes at the borders of that virtual cell.

The perturbation and therefore also the response are purely imaginary, so that there is no first order response density. This simplifies the equations and speeds up convergence.

NMR input

The options which control the specific NMR features are discussed below. None of them requires an argument, they can all be put in the same line.

 $&bull#bullet;$

RESTART: The run is restarted from a previous stop. The user has to take care that the required files RESTART.xxx (where xxx are NMR, p_x, p_y, p_z) exist.

 $&bull#bullet;$

CURRENT: Three density files containing current density values are written to disk. Further, the nucleus-independent chemical shift fields are written to disk. All in cube format.

 $&bull#bullet;$

NOSMOOTH: At the border of each virtual cell, the position operator is smoothened through an $\exp -r^2$. This option turns smoothing off. Can be useful if your cell is too small so that this smoothing would already occur in a region where the orbital density has not yet vanished.

 $&bull#bullet;$

NOVIRTUAL: If this keyword is specified, no individual virtual cells are used at all. All orbitals will have the same virtual cell, equal to the standard unit cell.

 $&bull#bullet;$

PSI0: With this option, the Wannier orbitals are plotted as cube files.

 $&bull#bullet;$

RHO0: With this option, the orbital densities are plotted as cube files.

#

OVERLAP: Overlapping orbitals (overlap criterion read from next line) are assigned the same virtual cells. Useful for isolated systems.

 $&bull#bullet;$

FULL A full calculation of the $\Delta j$ term in the NMR scheme is done. Relatively expensive, but quite important for accurate results in periodic systems.

NMR output

At the end of the six perturbation runs, the results are printed. This includes the magnetic susceptibility and the chemical shieldings. For the chemical shieldings, there are two values: the raw and the net ones. The raw shieldings correspond to a molecule in vacuum, without the susceptibility correction, whereas the net shieldings contains that correction for a spherical sample. It consists of an additive constant to all eigenvalues, which is printed out at the end of the net shieldings.

In more detail, the results are:

• The magnetic susceptibility tensor in both SI and cgs units. This quantity is an extensive quantity, i.e. two molecules in the simulation box will give twice the susceptibility of one.

• The raw shielding matrices of all atoms and its eigenvalues. Generally, the shielding tensor is not symmetric. To obtain unique eigenvalues, it is symmetrized ( $A_{ij} \mapsto (A_{ij}+A_{ji})/2$) before the diagonalization.

  All values are given in ppm, parts per millon. Chemical shieldings are dimensionless quantities.

• The principal values (eigenvalues) of the raw and net shielding tensors. As mentioned above, they only differ by an additive constant, the susceptibility correction for a spherical sample, which is printed out at the end of the list. The numbers are the shielding eigenvalues from most negative to most positive, the isotropic shielding (the average of the eigenvalues), and the anisotropy (the difference between the most positive and the average of the other two).

What to look at? If you search for the values which are peaked in a spectrum, you have to take the isotropic shieldings (the iso column of the output). If the system is a gas, take the raw shielding, if it is condensed matter, take the net shielding. If your system is a molecule in vacuo, but the experimentalists measure it in a solvent, add the susceptibility correction to the raw shieldings by yourself.

Why are my numbers so strange in absolute value? One more point shall be mentioned: For all nuclei except hydrogen, pseudopotentials screen the core electrons. The chemical shielding, however, is very sensitive to core and semi-core electrons. This can be corrected through a semi-empirical additive constant in many cases. This constant still needs to be added to the values given from the program. It depends on the nucleus, on the pseudopotential, and on the xc-functional.

In other cases, the calculated chemical shieldings are completely meaningless due to this deficiency. Then, you have to use a pseudopotential which leaves out the semicore states such that they are correctly taken into account. Example: carbon shieldings can be corrected very well through a constant number, silicon shieldings cannot. For Si, you have to take the $n=3$ shell completely into the valence band, requiring a cutoff larger than 300Ry.

How to compare to experiment? Usually, experimentalists measure the difference between the resonance frequency of the desired system and that of a reference system, and they call it $\delta$ (the shift) instead of $\sigma$ (the shielding). To make life more complicated, they usually define the shift of nucleus A of molecule X with respect to reference molecule ref as $\delta^{\text{\sffamily A}}_{\text{\sffamily ref}}({\text{\sffamily X}}) = \si... ... A}}({\text{\sffamily ref}})- \sigma^{\text{\sffamily A}}({\text{\sffamily X}})$. Example: To calculate $\delta^{\text{\sffamily H}}_{\text{\sffamily TMS}}({\text{\sffamily CH}}_4)$, where TMS=tetramethylsilane, the standard reference molecule for H-shifts, one would have to calculate the H-shielding of TMS and of CH$_4$ and substract them. Unfortunately, TMS is nontrivial to calculate, because it is a large molecule and the geometry is complicated (and the shieldings probably must be calculated taking into account vibrational and rotational averaging). Thus, in most cases it is better to take for instance the CH$_4$ shielding as a (computational) reference, and transform the shieldings relative to CH$_4$ to those relative to TMS through the experimental shielding of CH$_4$ with respect to TMS.

While doing so, you should not forget that the shielding is a property which is usually not yet fully converged when energies and bonding are. Therefore, the reference molecule should be calculated with the same computational parameters as the desired system (to reproduce the same convergence error). In particular, computational parameters include the type of the pseudopotential and its projectors, the xc-functional and the cutoff.

What accuracy can I expect? This is a difficult question, and there is no overall answer. First, one has to consider that on the DFT-pseudopotential level of theory, one will never reach the accuracy of the quantum chemistry community. However, for ``normal'' systems, the absolute accuracy is typically 0.5-1 ppm for hydrogen and about 5-20ppm for Carbon. The spread between extreme regions of the carbon spectrum is not reached: instead of 200ppm, one only reaches 150ppm between aliphatic and aromatic atoms, for instance. The anisotropy and the principal values of the shielding tensor can be expected to be about 10-20% too small. For hydrogen shieldings, these values are usually better, the error remains in the region of a single ppm.