17.1 NMR Shieldings of Closed-Shell Systems

17.1.1 Overview of NMR Shieldings

At present, the following methods are implemented in mpshift.

HF-SCF: the coupled-perturbed Hartree–Fock (CPHF) equations in the AO basis are solved using a semi-direct iterative algorithm [40] similar to dscf. The (multipole-accelerated-) resolution-of-the-identity-fitting approximation is available for the Coulomb term (MARI-J) [45]. Due to the effective screening approach the exchange part shows a low-order scaling. Full OpenMP support.
2c HF-SCF: the generalized two-component CPHF equations are solved as outlined in [38, 44]. The RI-J and the seminumerical exchange approximation are commended for efficiency. Full OpenMP support.
DFT: using either non-hybrid functionals where no iterations are needed for the coupled-perturbed Kohn–Sham (CPKS) equations [318] or hybrid functionals where the same algorithm as at the HF-SCF level is used. 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. By default, this can be done with the external vector potential [319, 45]. However, it is also possible to use the paramagnetic current density. We strongly recommend to use the current-dependent generalization for meta-GGAs and local hybrid functionals in NMR calculations ($curswitchengage) [320, 145, 321]. Starting with V8.0, the current-dependent generalization is the default and can be turned off with $curswitchdisengage. Then, the external vector potential is used for the generalization. Full OpenMP support.
2c DFT: using semilocal or (range-separated) hybrid functionals in a relativistic two-component formalism [44]. Note that this requires to iteratively solve the CPKS equations even for non-hybrid functionals due to spin–orbit effects and the non-vanishing spin-current densities of the ground state. The XC kernel can be neglected with the keyword $nmr_ziegler as suggested in [322]. However, including the kernel is recommended for the Spin-Orbit Heavy Atom on the Light Atom (SO-HALA) effect. The XC kernel should always be included for hybrid DFT calculations, as hybrid functionals always require an iterative procedure. Note that the current-dependent formalism for meta-GGAs is not available with GIAOs. Full OpenMP support.
MP2: semi-direct method, see ref. [41] and [323]. In contrast to HF and DFT, MP2 only comes with limited OpenMP capabilities. Additionally, the RI approximation is not yet supported for MP2.

The following Hamiltonians are available in addition to the usual non-relativistic one.

ECP: In molecules with ECP-carrying atoms, chemical shieldings on all the other atoms can be computed with mpshift in the way suggested by van Wüllen [324]. ECPs can be used to treat scalar-relativistic effects on neighboring atoms [46]. Spin–orbit ECPs are not available.
X2C: A scalar-relativistic or spin-free all-electron exact two-component Hamiltonian can be used to calculate the NMR shielding tensor of heavy elements [42]. A finite nucleus model based on a Gaussian charge distribution is available for the scalar potential and the vector potential. Grids with an increased number of radial points [183] (e.g. gridsize 4a) should be used. Moreover, it is recommended to use NMR-tailored basis sets, i.e. the x2c-XVPall-s (X=S,TZ) type bases. These employ additional tight functions to accurately sample the density in the vicinity of the nuclei. A local X2C Hamiltonian is further available and recommended. The diagonal local approximation to the unitary decoupling transformation (DLU) is employed. This results in a very efficient algorithm. The error introduced by the DLU scheme is negligible.
2c X2C: Spin–orbit two-component version of the previously mentioned X2C Hamiltonian [44]. Self-consistent treatment of spin–orbit interaction. DLU scheme is also available for efficiency. The finite nucleus model and the mSNSO approximation are further recommended.

17.1.2 NMR Chemical Shifts and Spectra

NMR shifts are obtained by comparing nuclear shieldings of your test compound with a reference molecule ($\delta_{\text{subst}} = \delta_{\text{ref}} + \sigma_{\text{ref}} - \sigma_{\text{subst}}$). Therefore, you have to choose a reference molecule with a well-known shift for which you can easily calculate the absolute shielding constant. This implies a certainty about the geometry, too. Furthermore, you have to use the very same basis set for corresponding atoms to minimize the basis set influence. Please note that the output is already given in units of ppm.

17.1.3 Prerequisites for NMR Shieldings

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 NMR shieldings typically require more flexible basis sets than necessary for geometries or energies. In X2C NMR shielding calculations, the use of tailored basis sets, i.e. x2c-SVPall-s, x2c-TZVPall-s [183] and x2c-QZVPall-s [184] etc., is recommended.
For HF or DFT calculations, no NMR specifications have to be made in the control file. Note that the m-grids such as gridsize m3 should not be used for DFT NMR or EPR calculations, as the response equations are only solved for the large grid. Please make sure to always use the full grid for the SCF iterations, i.e. gridsize 3 or 3a.
To perform an MP2 calculation of the NMR shieldings, you have to prepare the input with mp2prep -c.
mpshift is parallelized by OpenMP. The corresponding environment variables have to be set previously to running mpshift.

17.1.4 How to Perform a NMR Shielding Calculation with HF/DFT

All you have to do for running mpshift is typing mpshift at the shell level.

The results of a HF or DFT calculation (the trace of the total shielding tensors, its anisotropy and the CPHF contribution for each symmetry distinct atom) are written into the control file after the keyword $nmr <rhf/uhf/ghf/dft> shielding constants.

This data group is write-only for mpshift, but you can utilize it for graphical rendering of the calculated NMR spectra and for a quick overview of the results. A more detailed output with the complete shielding tensors can be found in the output of mpshift, so it is recommended to put the output in a file when calling the program.

For version 7.5 the solver for the CPHF equations has been reimplemented based on the Davidson algorithm used by escf and aoforce. This utilizes the keyword $shiftconv or $rpaconv (default 7, i.e. a residuum threshold of $10^{-7}$) to check for the convergence of the perturbed orbitals and density. Many new features made available with version 7.5 and later only support the Davidson solver. The 2c part of the program employs the Davidson solver of Ref. [38] and the same default for the norm of the residuum. Note that 2c calculations are only meaningful with the spin–orbit X2C Hamiltonian or its local variants. Spin-orbit ECPs are not supported.

A list of keywords for mpshift can be found in Section 25.2.31. Special features are as follows.

The mpshift program can be restarted at any stage of computing, since all intermediate results are written into the file restartcs. In case of an external program abort you have to remove the $actual step flag (by the command actual -r or using an editor). mpshift analyses this file and decides where to continue.
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 conductor-like screening model (COSMO) to account for counterions and solvation effects can be selected as done for the preceding HF or DFT calculation [45]. Note that COSMO only leads to additional terms due to the GIAOs.
NMR shielding tensors can be calculated for given nuclei only by setting $nucsel in the control file followed by the number of the nucleus of interest, e.g., $nucsel 1,3,7. Alternatively, you can set the element, e.g., $nucsel "c","h".
The default maximum number of iterations for the CPHF procedure is set to 35. It can be increased by adapting $csmaxiter in the control file.
Seminumerical exchange can be used for substantial speed-ups in HF and hybrid DFT calculations, i.e. $esenex is recommended for large-scale calculations.
The common gauge-origin variant is used with $cgo_nmr and the common gauge origin is specified with $cgo followed by the number of the atom. If atom number zero is given, the center of mass is chosen as the gauge origin. With $-1$, the Cartesian origin is selected. By default, the heaviest atom is used as the gauge origin. The GIAO version is used by default.
Nucleus-independent chemical shifts (NICS) can be calculated by setting $nics in the control file and then listing the Cartesian coordinates in atomic units just as the coordinates of the molecule itself. The NICS coordinates and isotropic as well as diagonal components of the tensor are stored in a file called nics_dft (or rhf/uhf/mp2), and the complete tensor is given in mB convention in the file nicstensor_dft and likewise for HF/MP2. For further details, see Section 25.2.31.

17.1.5 How to Perform a MP2 NMR Shielding Calculation

To perform an MP2 calculation of the NMR shieldings you have to prepare the input with mp2prep -c. mpshift will then calculate both the HF-SCF and MP2 shielding constants. The results are written into the control file after the keywords $nmr rhf shielding constants and $nmr mp2 shielding constants, respectively. The script mp2prep will create the keywords for the control file

$csmp2
$thize      .10000000E+10
$mointunit
 type=intermed unit=61 size=0 file=halfint
 type=1112     unit=63 size=0 file=moint#1
 type=1122     unit=64 size=0 file=moint#j
 type=1212     unit=65 size=0 file=moint#k
 type=1212a    unit=70 size=0 file=moint#a
 type=gamma#1  unit=71 size=0 file=gamma#1
 type=gamma#2  unit=72 size=0 file=gamma#2
 type=dtdb#1   unit=76 size=0 file=dtdb#1
 type=dtdb#2   unit=77 size=0 file=dtdb#2
$traloop 1
$statistics mpshift

and starts a statistics run of mpshift (by calling mpshift). If the resulting disk space requirement exceeds the automatically detected free disk space on your system, it will increase $traloop and run a statistics run again. This will be done as long as your free disk space is not sufficient for the calculation.

If the mp2prep script fails to run on your system, try to use the -p option or do the procedure described above by hand. Call mp2prep -h for more information about mp2prep.

17.1.6 Interfaces of Mpshift to Other Modules or Programs

mpshift can be combined with the other modules or programs to calculate various magnetic properties such as VCD, magnetizabilities, or ring currents.

Vibrational circular dichroism (VCD) spectra can be calculated using the gallier script utilzing the aoforce module [46].
The (un)perturbed density matrix can be stored on disk by setting $gimic in the control file. These matrices are required as input for the gauge-including magnetically induced currents (GIMIC) method [325, 326], see https://github.com/qmcurrents/gimic/. This method allows to study electron delocalization pathways and to estimate the degree of aromaticity [327], or to compute the magnetizability [328, 145].

17.1.7 Known Limitations of Mpshift

Molecular point groups that contain reducible e representations are not supported (C$_n$, C$_{n\text{h}}$ S$_n$, T, and T$_h$ with $n>2$).
The following features of mpshift are not available for open-shell systems: Calculations of shieldings at the MP2 level of theory, usage of the old CPHF/CPKS solver, calculation of VCD spectra. NICS calculations can only be performed for the orbital contribution of the shielding tensor. Two-component calculations are restricted to closed-shell systems.
The current-dependent generalization for meta-GGAs is not yet available for 2c NMR calculations with GIAOs. Here, a common gauge origin has to be set [43].
2c NMR calculations with local hybrid functionals are restricted to the common gauge origin ansatz [43].