17.2 NMR Coupling Constants of Closed-Shell Systems

17.2.1 Overview of NMR Couplings

At present, the following methods are implemented for NMR spin-spin coupling constants (SSCCs) with full OpenMP support in one- or two-component calculations with escf.

HF-SCF

the coupled-perturbed Hartree–Fock (CPHF) equations in the AO basis are solved using an iterative algorithm. The RI-K and seminumerical exchange approximations are available [143]

DFT

is available with functionals up to the class of local hybrids [143]. XCFun and LibXC with range-separated functionals are supported, see also chapter 6. Meta-GGAs and local hybrid functionals require the generalized kinetic energy density [145]. The current-dependent generalization is the default and can be turned off with $curswitchdisengage.

\(GW\)-BSE

the BSE equations are solved using the precomputed \(GW\) quasiparticles [229]. Here, the RI-K approximation is mandatory for the screened Coulomb interaction of BSE.

The following Hamiltonians are available.

Schrödinger

All four terms of Ramsey’s theory [329] (Fermi-contact, spin–dipole as well as paramagnetic and diamagnetic spin–orbit contributions, abbreviated FC, SD, PSO, and DSO, respectively) can be calculated. The FC/SD cross contributions are included in the full tensor [143]. To get reliable results, a basis set with steep \(s\) functions (for example Jensen’s pcJ-\(n\) basis sets [330, 331]) is needed.

X2C

The scalar X2C Hamiltonian and its local variants are supported for all four terms [186]. Note that the FC and SD terms are coupled and consequently the FC/SD cross contribution is not explicitly evaluated. In contrast, the PSO and DSO terms can be calculated separately. The restricted kinetic balance (RKB) condition and a finite nucleus model (if $finnuc is set) for both the scalar and the vector potential are employed. In a finite nucleus model, steep \(s\) functions are not that important and the use of x2c-TZVPall [182] or Dyall’s uncontracted CVTZ [332, 333, 334] basis sets is usually sufficient. Note that the one-component approach leads to comparably large errors when the DSO term becomes pronounced, e.g., for H–H couplings. Here, complex algebra can be used to consider all X2C contributions. This is done by using the pseudo-two-component algorithm (dsop2c), which is also the default option for scalar-relativistic calculations. This does not require a two-component SCF solution [186]. Simply replace ‘dso’ in the $ncoupling group of the control file by ‘dsop2c’ (see below).

2c X2C

Spin-orbit X2C calculation. As done for the scalar X2C approach, the restricted kinetic balance (RKB) condition and a finite nucleus model for both the scalar and the vector potential can be employed. The FC, SD, and PSO terms of Ramsey’s theory are coupled and can no longer be calculated individually. Due to the X2C response, a DSO-like term is recovered (despite using RKB), but which is expensive and can therefore be turned off. Using the picture-change correction ($dsopcc) in a manner similar to Refs. [38] and [335] is a cheap alternative [336]. The use of x2c-TZVPall-2c [182] or Dyall’s uncontracted CVTZ [332, 333] basis sets is recommended.

17.2.2 NMR Coupling Constants and Spectra

The gyromagnetic ratios of the coupling nuclei can easily be changed in the atomic attributes menu of the define module. The same holds for the isotopes, see also Sec. 25.2.1.

The output will include the isotropic part of the coupling in Hz (cf. Ref. [143]),

\[\begin{equation} J_{KL}^{\text{iso}} = \frac{1}{3} \sum_{\alpha = x,y,z} (J_{KL})_{\alpha\alpha}\,, \end{equation}\](17.1)

the anisotropic contribution,

\[\begin{equation} J_{KL}^{\text{anis}} = h \frac{\gamma_K}{2\pi} \frac{\gamma_L}{2\pi} \sqrt{\frac{3}{2}\Big(\frac{1}{4}\sum_{\alpha\beta}\big((K_{KL})_{\alpha\beta} + (K_{KL})_{\beta\alpha}\big)^2 - 3(K_{KL}^{\text{iso}})^2\Big)}\,, \end{equation}\](17.2)

and the complete tensor \(J_{KL}\) according to \[\begin{equation} J_{KL} = h \frac{\gamma_K}{2\pi} \frac{\gamma_L}{2\pi} K_{KL} \end{equation}\](17.3) with the reduced coupling tensor (see Eq. 8.23 in Sec. 8, linear response) \[\begin{equation} K_{KL} = \sum_{ai} \lambda_{ai,K} R_{ai,L} = \sum_{ai} R_{ai,K} \lambda_{ai,L}\,, \end{equation}\](17.4)

17.2.3 Prerequisites for NMR Couplings

  1. mpshift needs converged MO vectors from a HF or DFT run of dscf or ridft. The 2c runs need converged spinor vectors from an ridft calculation. The flag $coulex can set to use analytical Coulomb integrals in 2c ridft and mpshift calculations. Note that the m-grids such as gridsize m3 should not be used for DFT NMR calculations. Make sure to always use the full grid for the SCF iterations, i.e. gridsize 3 or 3a.

  2. To perform an \(GW\)-BSE calculation of the NMR couplings, you have to calculate the quasiparticle energies as discussed in Sec. 14 and these have to be stored on disk.

  3. escf is parallelized by OpenMP. The corresponding environment variables have to be set previously to running escf.

17.2.4 How to Perform a NMR Coupling Calculation with HF/DFT

SSCCs require a converged HF or DFT calculation and the respective molecular orbitals, just like an NMR shift calculation. The SSCC calculation is triggered by the keyword $ncoupling in the control file (no $scfinstab is needed). This keyword along with several options can be set in the last menu in define. Generally, the terms of interest are inserted and a threshold for the output is given. This avoids listing coupling constants which are almost zero. The data group $ncoupling can be set using the ncoup section in define. A sample input might look like this:

$ncoupling
  fc sd pso dso nofcsdcross     or
  simple                        or
  fromfile
  reduced
  thr=0.1
$nucsel 1,3,5-8
$nucsel2 "c","h"

You can specify the contributions to be calculated by indicating the appropriate abbreviation. If both fc and sd are specified, the FC/SD cross term is calculated, unless manually switched off.

The simple method only calculates the FC and FC/SD cross terms, thus yielding the most important contributions to the isotropic and anisotropic part. Be warned that this might yield qualitatively wrong results in some cases! This method is quite fast because only response equations for the FC term are solved. However, this means that it is incompatible with $nucsel2 (see below).

If the option fromfile is set, SSCCs from the data group $coupling_reduced from a previous calculation are used. This option can be used to obtain SSCC for a different isotope without redoing the expensive part of the calculation.

If none of the keywords mentioned above is given, all terms (equivalent to fc sd pso dso) are calculated.

In a two-component calculation, fc, sd, and pso shall always be used together and nofcsdcross is not permissible. The simple method is likewise not allowed. If the DSO term is to be computed by picture-change correction, it needs to be activated here and $dsopcc needs to be set.

By default, the output is multiplied by the gyromagnetic ratios, yielding \(J\) in Hertz. With reduced, the reduced coupling constants \(K\) are printed in units of \(10^{19}\,\text{T}^2/\text{J}\). The content of $coupling_reduced is not influenced by this setting.

In the final output, couplings where both the isotropic and the anisotropic part are smaller than the threshold thr (in Hertz or \(10^{19}\,\text{T}^2/\text{J}\), default: 0.1) will not be printed.

You can specify for which atoms coupling constants shall be calculated by the data groups $nucsel and $nucsel2. As shown above, the syntax is either a list of atomic indices or of element symbols. The default is to calculate all atoms. Response equations (the most expensive step) are solved for the atoms in $nucsel and right-hand side integrals are calculated for atoms in any of the data groups. This means that a coupling tensor is obtained if at least one of the partner is in $nucsel. Notice that the very same data group $nucsel is used by the module mpshift.

It is furthermore recommended to set $rpaconv to at least 6. For large-scale calculations, please adapt the data group $maxcor in the control file. This will allow for a more efficient batching in the Davidson algorithm.

The file referenced in $coupling_reduced (default file name: coupling.reduced) is suitable for automated post-processing. It consists of the atomic indices, \(J^{\text{iso}}\), \(\gamma_K \cdot \gamma_L\), and the nine elements of the matrix \(\tilde{K} = J/\gamma_K \gamma_L\) (see eqs. 17.4 and 17.3 for an explanation of these quantities). The gyromagnetic ratios \(\gamma_K\) are given in \(10^7 \,\text{rad}\,\text{s}^{-1}\,\text{T}^{-1}\). All couplings with \(K \ge L\) (ignoring $nucsel and thr) are printed.

Starting with Version 7.6, the kinetic energy density is generalized with the paramagnetic current density by default [145]. This affects only the paramagnetic spin–orbit term. The generalization can be disabled by $curswitchdisengage.

In a relativistic two-component calculation [185], the desired 2c keywords have to be set in addition to $ncoupling. The calculation of spin–spin coupling constants is currently restricted to closed-shell systems. Thus, the keyword $kramers should be set for the ground-state calculation. The current-dependent generalization of the response of the kinetic energy density is used by default [185, 145, 229], whereas the ground-state current density is only included with $curswitchengage in the ridft calculation. fc, sd, and pso shall always be used together and nofcsdcross is not permissible. The simple method is likewise not allowed. If the DSO term is to be computed by picture-change correction, it needs to be activated here with $dsopcc.

17.2.5 How to Perform a NMR Coupling Calculation with \(GW\)-BSE

The GW-BSE formalism is available for non-relativistic, scalar-relativistic, and relativistic two-component calculations [229]. This only requires GW quasiparticle energies and the keyword $bse or $cbse in addition to the SSCC settings. Note that the RI-K approximation is mandatory for the BSE part. The data group $ncoupling is set as done for HF/DFT. Therefore, the workflow is as follows.

  1. Calculate the molecular orbitals at the HF or DFT level dscf or ridft.

  2. Perform an \(GW\) calculation to compute the quasiparticles with escf.

  3. Run the BSE or cBSE calculation for the NMR coupling constants with the RI-K approximation in escf.