Calculation of NMR shieldings using simple magnetic balance

In this tutorial we look at 4-component relativistic calculations of NMR shielding constants. For best results we recommend to use the scheme of simple magnetic balance (sMB) in conjunction with London orbitals [Olejniczak2012].

Basic theory of magnetic balance

We encourage the reader to consult the original paper [Olejniczak2012] for full details; here we just provide some key points.

In the absence of magnetic fields the coupling between the large and small components of the Dirac equation is

\[c \psi^S = \frac{1}{2m} \left[1 + \frac{E-V}{2mc^2} \right]^{-1} (\vec \sigma \cdot \vec p) \psi^L\]

The exact coupling can not be implemented per se at the basis set level due to the energy dependence of the coupling. One therefore employs the non-relativistic limit of this coupling

\[\lim_{c \to \infty} c \psi^S = \frac{1}{2m} (\vec \sigma \cdot \vec p) \psi^L,\]

valid for positive-energy solutions (and extended nuclear charge distributions). Construction of the small component basis in this manner assures the correct representation of the kinetic energy operator in the non-relativistic limit, and the coupling is therefore referred to as kinetic balance [Stanton1984]. In a 2-spinor basis kinetic balance can be implemented in such a manner that there is a 1:1 ratio between the sizes of the large and small component basis; this is referred to as restricted kinetic balance (RKB).

\[\chi_{\mu}^{2c;S} = (\vec \sigma \cdot \vec p) \chi_{\mu}^{2c;L}\]

In a scalar basis the straightforward implementation of RKB is not possible and one typically resorts to unrestricted kinetic balance (UKB)

\[\{ \chi_{\mu}^{1c;S} \} = \{ \vec p \chi_{\mu}^{1c;L} \},\]

basically generating the small component basis functions as derivative of the large component ones.

Magnetic fields are introduced through minimal substitution

\[\vec p \to \vec \pi = \vec p + e \vec A\]

which manifestly modifies the coupling between the large and the small components. The non-relativistic limit now reads

\[\lim_{c \to \infty} c \psi^S = \frac{1}{2m} (\vec \sigma \cdot \vec \pi) \psi^L\]

and is referred to as magnetic balance [Aucar1999]. Magnetic balance needs to be taken into account when constructing the first-order correction to the wave function upon inclusion of an external magnetic field and is not assured in a RKB basis [Pecul2004]. In the cited paper [Olejniczak2012] we show that magnetic balance can be obtained in a Gaussian basis combining UKB with London orbitals. However, UKB generally leads to a small component basis larger than the large component basis set from which it was generated and may increase computational cost as well as introduce linear dependencies in the basis set as compared to RKB. We therefore propose calculations of NMR shielding parameters by a simple two-step procedure described in the next section.

A sample calculation

We will consider a 4-component relativistic density functional calculation of the NMR shielding constant of hydrogen fluoride using the KT2 functional [Keal2003]. The molecular input file hf.mol reads:

DIRAC


C   2              A
        1.    1
H       0.0000               0.0000               0.0000
LARGE BASIS aug-cc-pVTZ
        9.    1
F       0.0000               0.0000               0.9168
LARGE BASIS aug-cc-pVTZ
FINISH

Step 1: SCF calculation

The first step is the generation of the Kohn–Sham wave function in a RKB basis. Although DIRAC employs scalar Gaussian basis functions restricted kinetic balance is by default imposed when transforming to orthonormal basis [Visscher2000].

The menu file scf.inp reads:

**DIRAC
.WAVE FUNCTION
**INTEGRALS
*READIN
.UNCONTRACT
**HAMILTONIAN
.DFT
 KT2
**WAVE FUNCTION
.SCF
*SCF
.CLOSED SHELL
 10
**END OF

and we make sure to recover the coefficient file DFCOEF when running the calculation:

pam --inp=scf.inp --mol=hf.mol --outcmo

Step 2: Response calculation

In the second step we start from the RKB coefficient obtained in the SCF calculation, but DIRAC will add the UKB complement and calculate the NMR shielding constant using unrestricted kinetic balance combined with London orbitals. The menu file response.inp accordingly reads:

**DIRAC
.PROPERTIES
**GENERAL
.RKBIMP
**INTEGRALS
*READIN
.UNCONTRACT
**HAMILTONIAN
.URKBAL
.DFT
 KT2
**WAVE FUNCTION
.SCF
*SCF
.CLOSED SHELL
 10
**PROPERTIES
.SHIELDING
*NMR
.LONDON
**END OF

and the pam command is:

pam --inp=response.inp --mol=hf.mol --incmo

Comparison with RKB calculations

With the above computational procedure we obtain the following NMR shielding constants:

           isotropic shielding

       ----------------------------

atom          total          dia         para

---------------------------------------------------

H            29.9502       21.8737        8.0765

F           415.8316      480.2015      -64.3699

---------------------------------------------------

By not adding the UKB complement (deleting .RKBIMP and .URKBAL), that is, using RKB, gives:

            isotropic shielding

        ----------------------------

atom          total          dia         para

--------------------------------------------------------

H            28.8028       20.7264        8.0765

F           323.2672      387.6459      -64.3787

--------------------------------------------------------

If we in addition switch off the use of London-orbitals (deleting .LONDON) we obtain:

            isotropic shielding

        ----------------------------

atom          total          dia         para

------------------------------------------------------

H            28.2920       20.2563        8.0357

F           322.9065      387.1375      -64.2310

------------------------------------------------------

The experimental value for the absolute isotropic shielding available for H in HF is 28.72 ppm [Chesnut1994].