3.46 Magnetic Response: Nuclear Magnetic Resonance and Magnetizability

FHI-aims is capable of perturbatively calculating the response of a molecular system to a magnetic field. The magnetic field can be either external or that associated with the magnetic moments of atomic nuclei. The functionality is currently restricted to non-periodic geometries.

Please also consult the key reference that describes the magnetic response functionality in FHI-aims, including detailed basis set convergence studies and recommendations: Ref. [184].

Currenty functionality includes:

•

Nuclear magnetic shielding tensors and the spin-spin coupling tensors, which are the central parameters in nuclear magnetic resonance (NMR) spectroscopy;
•

The magnetizability tensor, which is the proportionality between the induced magnetic moment and the magnetic field. It can be thought of as the magnetic analogue of the polarizability;
•

Electric field gradient (experimental).

References:

•

Raul Laasner, Iuliia Mandzhieva, William P, Huhn, Johannes Colell, Victor Wen-zhe Yu, Warren S. Warren, Thomas Theis, Volker Blum, Molecular NMR shieldings, J-couplings, and magnetizabilities from numeric atom-centered orbital based density-functional calculations, Electronic Structure 6, 027002, https://dx.doi.org/10.1088/2516-1075/ad45d4 – FHI-aims implementation [184].
•

Martin Kaupp, Michael Bühl, and Vladimir G. Malkin, editors, Calculation of NMR and EPR Parameters: Theory and Applications, Wiley-VCH, Weinheim, Germany, 2004 — A comprehensive treatment of general theory and applications.
•

V. Sychrovsky et al., Nuclear magnetic resonance spin-spin coupling constants from coupled perturbed density functional theory, J. Chem. Phys. 113, 3530 (2000) — Nonrelativistic J-couplings in DFT.
•

G. Schreckenbach and T. Ziegler, Calculation of NMR Shielding Tensors Using Gauge-Including Atomic Orbitals and Modern Density Functional Theory, J. Phys. Chem. 99, 606 (1995) — Standard theory of GIAO shieldings in DFT.
•

T. Helgaker et al., Nuclear shielding constants by density functional theory with gauge including atomic orbitals, J. Chem. Phys. 113, 2983 (2000) — GIAO GGA correction (shieldings only).

The basic parameters of a magnetic response experiment are a collection of electron spins, $\bm{S}=\sum_{i}\bm{S}_{i}$ , nuclear spins, $\bm{I}=\sum_{i}\bm{I}_{i}$ , and an external magnetic field $\bm{B}$ . Considering all possible couplings between these parameters, the fully general effective Hamiltonian has the form

H_{S}=H(\bm{S},\bm{S})+H(\bm{S},\bm{I})+H(\bm{S},\bm{B})+H(\bm{I},\bm{I})+H(% \bm{I},\bm{B})+H(\bm{B},\bm{B}).

(3.232)

We have not explicitly considered coupling to the electron orbital angular momentum, $\bm{L}=\sum_{i}\bm{L}_{i}$ , whose effects are assumed to be incorporated into the electronic wavefunction. Based on which terms are dominant in a given spectroscopy, we can partition the Hamiltonian (3.232) as follows:

•

The NMR nuclear spin Hamiltonian:

$\begin{split}H^{\text{NMR}}=&H(\bm{I},\bm{B})+H(\bm{I},\bm{I})\\ =&-\sum_{A}\bm{B}(\overleftrightarrow{1}-\overleftrightarrow{\sigma_{A}})\bm{% \mu}_{A}+\sum_{A>B}h\bm{I}_{A}(\overleftrightarrow{D_{AB}}+\overleftrightarrow% {J_{AB}})\bm{I}_{B},\\ \end{split}$ (3.233)

where $\overleftrightarrow{\sigma}_{A}$ are the nuclear shielding tensors, $\overleftrightarrow{D_{AB}}$ and $\overleftrightarrow{J_{AB}}$ are the direct and indirect spin-spin coupling tensors, $\bm{I}_{A}$ is the nuclear spin of atom $A$ , and $\bm{\mu}_{A}=\hbar\gamma_{A}\bm{I}_{A}$ is the corresponding nuclear magnetic dipole moment.
•

The EPR spin Hamiltonian:

$\begin{split}H^{\text{EPR}}=&H(\bm{S},\bm{S})+H(\bm{S},\bm{I})+H(\bm{S},\bm{B}% ).\end{split}$ (3.234)
•

The effective Hamiltonian describing the second order response to an external magnetic field:

$H^{\text{Mag}}=-\frac{1}{2}\bm{B}\overleftrightarrow{\xi}\bm{B},$ (3.235)

where $\overleftrightarrow{\xi}$ is the magnetizability tensor. For a closed shell molecule, it describes the lowest order response of the energy to an external magnetic field.

The magnetic response parameters can be expressed as mixed second derivatives of the total electronic energy with respect to the external magnetic field and the magnetic moments,

\begin{split}\overleftrightarrow{\sigma_{A}}=&\left.\frac{\partial^{2}E(\bm{B}% ,\bm{\mu})}{\partial\bm{B}\partial\bm{\mu}_{A}}\right|_{\bm{B}=0;\bm{\mu}=0},% \\ \overleftrightarrow{J_{AB}}=&\frac{\gamma_{A}\gamma_{B}}{2\pi}\left.\frac{% \partial^{2}E(\bm{B},\bm{\mu})}{\partial\bm{\mu}_{A}\partial\bm{\mu}_{B}}% \right|_{\bm{B}=0;\bm{\mu}=0},\\ \overleftrightarrow{\xi}_{AB}=&-\left.\frac{\partial^{2}E(\bm{B},\bm{\mu})}{% \partial^{2}\bm{B}}\right|_{\bm{B}=0;\bm{\mu}=0},\end{split}

(3.236)

where $\bm{\mu}$ denotes the collection of all magnetic moments. In the first of Eqs. (3.236), the total energy should not contain the classical nuclear Zeeman term. The derivatives are, to a very good approximation, taken at zero values of $\bm{B}$ and $\bm{\mu}$ , which allows us to consider only the linear response of the system.

In time-independent second-order perturbation theory, the general expression for a property $\overleftrightarrow{X}$ is

\overleftrightarrow{X}=\sum_{i}^{\text{occ}}\left[\braket{\psi_{i}}{H^{MN}}{% \psi_{i}}+2\mathrm{Re}\braket{\psi_{i}^{M}}{H^{N}}{\psi_{i}}\right],

(3.237)

where $\ket{\psi_{i}}$ are the unperturbed wavefunctions, $\ket{\psi^{M}}$ are the first-order wavefunctions with respect to the perturbation $M$ , and $H^{M}$ and $H^{MN}$ are the first and second derivatives of the Hamiltonian. The derivatives are calculated at zero values of the perturbation ( $\bm{B}$ and/or $\bm{\mu}$ ). The total number of terms depends on the type of perturbation. For instance, in case of J-couplings, the paramagnetic part is a sum of three separate terms, $H^{N}=H^{\text{FC}\bm{\mu}_{A}}+H^{\text{PSO}\bm{\mu}_{A}}+H^{\text{SD}\bm{\mu% }_{A}}$ , which are listed below.

Density functional perturbation theory

The first-order wavefunctions are calculated according to density functional perturbation theory (DFPT) (also called coupled-perturbed Kohn-Sham or coupled-perturbed Hartree-Fock). The DFPT steps are the following:

•

Update the self-consistent first-order Hamiltonian:

$H_{\sigma}^{(1)}(\bm{r})=H_{\text{NSC},\sigma}^{(1)}(\bm{r})+\left[\frac{% \partial v_{xc,\sigma}[n]}{\partial n_{\sigma}}n_{\sigma}^{(1)}(\bm{r})+\frac{% \partial v_{xc,\sigma}[n]}{\partial n_{\sigma^{\prime}}}n_{\sigma^{\prime}}^{(% 1)}(\bm{r})\right]_{n=n(\bm{r})}+H_{\text{Ha}}^{(1)}n_{\sigma}(\bm{r}),$ (3.238a)

where $H_{\text{NSC}}^{(1)}$ is the non-self-consistent part of the first-order Hamiltonian, $v_{xc,\sigma}$ is the xc-potential, $H_{\text{Ha}}^{(1)}$ is the first-order Hartree potential, and $n_{\sigma}^{(1)}$ is the first-order density. First-order Hartree potential is computed only for an open-shell system and vanishes otherwise. For a closed shell system, only $n_{\sigma}^{(1)}$ in Eq. (3.238a) is updated, the rest need not be calculated more than once. In fact, the xc-kernel, which depends only on the unperturbed density, is the same for any perturbation and is calculated only once during the entire calculation. In Eq. (3.238a), only the LDA contribution is shown, but GGA and hybrids are also available. Note that for purely imaginary operators the perturbed density vanishes, which means that in LDA and GGA there is no contribution from the xc-kernel and it is sufficient to perform only a single DFPT step (this is the sum-over-states approach). With hybrids or Hartree-Fock, this does not hold because while the first-order density vanishes, the first-order density matrix does not.
•

Compute the first-order wavefunctions:

$\psi_{n}^{(1)}(\bm{r})=\sum_{i}^{\text{occ}}C_{in}^{(1)}\phi_{i}(\bm{r})=\sum_% {i}^{\text{occ}}\phi_{i}(\bm{r})\sum_{k}C_{ik}U_{kn},$ (3.238b)

where $\phi_{i}(\bm{r})$ are the basis functions, $\bm{C}$ are the unperturbed wavefunction coefficients, $\bm{C}^{(1)}$ are the first-order coefficients, and

$U_{mn}=-\frac{1}{2}\sum_{ij}C_{im}S_{ij}^{(1)}C_{jn}$ (3.238c)

if $m n$ is an occupied-occupied or virtual-virtual pair, and

$U_{mn}=\sum_{ij}\frac{C_{im}H_{ij}^{(1)}C_{jn}-C_{im}S_{ij}^{(1)}C_{jn}% \epsilon_{n}}{\epsilon_{n}-\epsilon_{m}}$ (3.238d)

otherwise. $S^{(1)}$ is the first-order overlap matrix.
•

Construct a new density matrix from the first-order changes in wavefunctions:

$n_{ij\sigma}^{(1)}=\sum_{m\sigma}^{\text{occ}}\braket{\phi_{i}}{\psi_{m\sigma}% }\braket{\psi_{m\sigma}^{(1)}}{\phi_{j}}+\braket{\phi_{i}}{\psi_{m\sigma}^{(1)% }}\braket{\psi_{m\sigma}}{\phi_{j}}=\sum_{m\sigma}^{\text{occ}}c_{im}c_{jm}^{(% 1)*}+c_{im}^{(1)}c_{jm}.$ (3.238e)

The unperturbed coefficients are always real.
•

Check convergence and, if necessary, perform density matrix mixing (default is to do Pulay mixing):

$n_{ij\sigma}^{(1)}\rightarrow n_{ij\sigma}^{X,\text{next}}.$ (3.238f)
•

Construct a new density from the density matrix:

$n_{\sigma}^{(1)}(\bm{r})=\sum_{ij}^{\text{occ}}\phi_{i}(\bm{r})n_{ij\sigma}^{(% 1)}\phi_{j}(\bm{r}).$ (3.238g)

Nonrelativistic integrals

In the presence of an external magnetic field and NMR active nuclei (those with a nonzero magnetic moment) in the system, the magnetic vector potential takes the form

\bm{A}=\frac{1}{2}\bm{B}\times\bm{r}+\alpha^{2}\sum_{A}\frac{\bm{\mu}_{A}% \times\bm{r}_{A}}{r_{A}^{3}},

(3.239)

where $\bm{r}_{A}=\bm{r}-\bm{R}_{A}$ is the position relative to atom $A$ . Inserting this into the nonrelativistic Hamiltonian and taking the first and second derivatives with respect to the magnetic field and/or the nuclear magnetic moments, we arrive at the following terms:

•

Orbital Zeeman:

$H^{\bm{B}}=-\frac{i}{2}(\bm{r}\times\bm{\nabla}),$ (3.240a)
•

Spin Zeeman:

$H^{\text{SZ}\bm{B}}=\bm{S},$ (3.240b)
•

Paramagnetic spin-orbit (PSO):

$H^{\text{PSO}\bm{\mu}_{A}}=-i\alpha^{2}\frac{\bm{r}_{A}\times\bm{\nabla}}{r_{A% }^{3}},$ (3.240c)
•

Fermi contact (FC):

$H^{\text{FC}\bm{\mu}_{A}}=\frac{8\pi\alpha^{2}}{3}\delta(\bm{r}_{A})\bm{S},$ (3.240d)
•

Spin-dipole (SD):

$H^{\text{SD}\bm{\mu}_{A}}=\alpha^{2}\left[\frac{3(\bm{S}\bm{r}_{A})\bm{r}_{A}}% {r_{A}^{5}}-\frac{\bm{S}}{r_{A}^{3}}\right],$ (3.240e)
•

Diamagnetic shielding (DS):

$H^{\bm{B}\bm{\mu}_{A}}=\frac{\alpha^{2}}{2}\frac{\bm{r}\bm{r}_{A}-\bm{r}_{A}% \bm{r}^{T}}{r_{A}^{3}},$ (3.240f)
•

Diamagnetic spin-orbit (DSO):

$H^{\bm{\mu}_{A}\bm{\mu}_{B}}=\frac{\alpha^{4}}{2}\frac{\bm{r}_{A}\bm{r}_{B}-% \bm{r}_{B}\bm{r}_{A}^{T}}{r_{A}^{3}r_{B}^{3}},$ (3.240g)
•

Diamagnetic magnetizability:

$H^{\bm{B}\bm{B}}=\frac{1}{4}(r^{2}-\bm{r}\bm{r}^{T}).$ (3.240h)

In case of shieldings and magnetizability, the GIAO formalism is used to overcome the slow convergence related to the nonuniqueness of the gauge origin. Each basis function $\phi_{n}$ is equipped with a phase factor of the form

\phi_{n}(\bm{r})\rightarrow e^{-i\bm{A}_{n}\bm{r}}\phi_{n}(\bm{r})=e^{-\frac{i% }{2}(\bm{R}_{n}\times\bm{r})\bm{B}}\phi_{n}(\bm{r}),

(3.241)

where $\bm{A}_{n}=\frac{1}{2}\bm{B}\times(\bm{r}-\bm{R}_{n})$ is the vector potential with the gauge origin shifted to the origin of basis function $n$ , $\bm{R}_{n}$ . This modifies some of the integrals (3.240) and leads to a few new ones:

•

GIAO orbital Zeeman:

$H_{mn}^{\text{GIAO}\bm{B}}=\frac{i}{2}\Braket{\phi_{m}}{(\bm{R}_{mn}\times\bm{% r})H_{0}-\bm{r}_{n}\times\bm{\nabla}}{\phi_{n}}+\braket{\phi_{m}}{V_{\text{GGA% },\alpha}^{\bm{B}}}{\phi_{n}}.$ (3.242a)
•

GIAO diamagnetic shielding:

$H_{mn}^{\text{GIAO}\bm{B}\bm{\mu}_{A}}=\frac{\alpha^{2}}{2}\Braket{\phi_{m}}{% \frac{(\bm{R}_{mn}\times\bm{r})(\bm{r}_{A}\times\bm{\nabla})^{T}}{r_{A}^{3}}}{% \phi_{n}}.$ (3.242b)

•

GIAO diamagnetic magnetizability:

\begin{split}H_{mn}^{\text{GIAO}\bm{B}\bm{B}}=&\frac{1}{4}\Braket{\phi_{m}}{2(% \bm{R}_{mn}\times\bm{r})(\bm{r}_{n}\times\bm{\nabla})^{T}}{\phi_{n}}+\frac{1}{% 4}\Braket{\phi_{m}}{r_{n}^{2}-\bm{r}_{n}\bm{r}_{n}^{T}}{\phi_{n}}\\ &-\frac{1}{4}\Braket{\phi_{m}}{(\bm{R}_{mn}\times\bm{r})(\bm{R}_{mn}\times\bm{% r})^{T}H_{\text{LDA}}}{\phi_{n}}+\Braket{\phi_{m}}{V_{\text{GGA}}^{\bm{B}\bm{B% }}}{\phi_{n}}.\end{split}

(3.242c)

•

First-order overlap matrix:

$S_{mn}^{\text{GIAO}\bm{B}}=\frac{i}{2}\braket{\phi_{m}}{\bm{R}_{mn}\times\bm{r% }}{\phi_{n}}.$ (3.242d)
•

Second-order overlap matrix:

$S_{mn}^{\text{GIAO}\bm{B}\bm{B}}=-\frac{1}{4}\braket{\phi_{m}}{(\bm{R}_{mn}% \times\bm{r})(\bm{R}_{mn}\times\bm{r})^{T}}{\phi_{n}}.$ (3.242e)

In Eqs. (3.242), $H_{0}$ is the unperturbed DFT Hamiltonian, $V_{\text{GGA}}^{\bm{B}}$ and $V_{\text{GGA}}^{\bm{B}\bm{B}}$ are the first and second $\bm{B}$ derivatives of the xc-potential. In LDA, the latter terms vanish.

Table 3.2 illustrates which terms are required for calculating which properties.

Table 3.2: The number of different types of integrals required for J-couplings, shieldings, and magnetizabilities is 4. Without GIAOs, the number is 3 for shieldings and 2 for magnetizabilities.

	$\overleftrightarrow{J_{AB}}$				$\overleftrightarrow{\sigma_{A}}$ (no GIAO)		$\overleftrightarrow{\xi}$ (no GIAO)
$\ket{\psi^{M}}$	$H^{\text{FC}\bm{\mu}_{A}}$	$H^{\text{PSO}\bm{\mu}_{A}}$	$H^{\text{SD}\bm{\mu}_{A}}$		$H^{\bm{B}}$		$H^{\bm{B}}$
$H^{M}$	$H^{\text{FC}\bm{\mu}_{A}}$	$H^{\text{PSO}\bm{\mu}_{A}}$	$H^{\text{SD}\bm{\mu}_{A}}$		$H^{\text{PSO}\bm{\mu}_{A}}$		$H^{\bm{B}}$
$H^{MN}$				$H^{\bm{\mu}_{A}\bm{\mu}_{B}}$		$H^{\bm{B}\bm{\mu}_{B}}$		$H^{\bm{B}\bm{B}}$

	$\overleftrightarrow{\sigma_{A}}$ (GIAO)		$\overleftrightarrow{\xi}$ (GIAO)
$\ket{\psi^{M}}$	$H^{\text{GIAO}\bm{B}},S^{\text{GIAO}\bm{B}}$		$H^{\text{GIAO}\bm{B}},S^{\text{GIAO}\bm{B}}$
$H^{M}$	$H^{\text{PSO}\bm{\mu}_{A}}$		$H^{\text{GIAO}\bm{B}}$
$H^{MN}$		$H^{\text{GIAO}\bm{B}\bm{\mu}_{A}}$		$H^{\text{GIAO}\bm{B}\bm{B}},S^{\text{GIAO}\bm{B}\bm{B}}$

Relativistic integrals

Following J. Autschbach and T. Ziegler, J. Chem. Phys. 113, 936 (2000), the diamagnetic term is almost the same as in the non-relativistic case except for the ZORA factor $\mathcal{K}_{\text{at}(j)}$ :

\begin{split}\bm{H}^{\text{ZORA},\bm{\mu}_{A}\bm{\mu}_{B}}=&\frac{\partial^{2}% }{\partial\bm{\mu}_{A}\partial\bm{\mu}_{B}}\bm{A}\mathcal{K}_{\text{at}}\bm{A}% =2\alpha^{4}\mathcal{K}_{\text{at}}\frac{\bm{r}_{A}\bm{r}_{B}-\bm{r}_{B}\bm{r}% _{A}^{T}}{r_{A}^{3}r_{B}^{3}},\\ \mathcal{K}_{\text{at}(j)}=&\frac{c^{2}}{2c^{2}-v_{\text{at}(j)}^{\text{free}}% }.\end{split}

(3.243)

Note that because of $\mathcal{K}_{\text{at}(j)}$ Eq. (3.243) needs to be explicitly symmetrized which is different from the non-relativistic case. The paramagnetic spin-orbit (PSO) term is similar to its non-relativistic equivalent,

\bm{H}^{\text{ZORA,PSO}\bm{\mu}_{A}}=\frac{\partial}{\partial\bm{\mu}_{A}}% \left(\bm{p}\mathcal{K}_{\text{at}}\bm{A}+\bm{A}\mathcal{K}_{\text{at}}\bm{p}% \right)=-i\alpha^{2}\left(\mathcal{K}_{\text{at}}\frac{\bm{r}_{A}\times\bm{% \nabla}}{r_{A}^{3}}+\frac{\bm{r}_{A}\times\bm{\nabla}}{r_{A}^{3}}\mathcal{K}_{% \text{at}}\right),

(3.244)

but since the ZORA factor is orbital-dependent, it needs to be evaluated as two separate integrals,

\Braket{\phi_{m}}{\bm{H}^{\text{ZORA,PSO}\bm{\mu}_{A}}}{\phi_{n}}\\ =-i\alpha^{2}\bm{r}_{A}\times\left(\int\mathcal{K}_{\text{at(n)}}\frac{\phi_{m% }(\bm{r})\bm{\nabla}\phi_{n}(\bm{r})}{r_{A}^{3}}\mathrm{d}\,\bm{r}-\int% \mathcal{K}_{\text{at}(n)}\frac{\left[\bm{\nabla}\phi_{m}(\bm{r})\right]\phi_{% n}(\bm{r})}{r_{A}^{3}}\mathrm{d}\,\bm{r}\right),

(3.245)

where both terms need to be additionally symmetrized. The ZORA spin-orbit (ZSO) term is

H_{j}^{\text{ZORA,ZSO}\bm{\mu}_{A}}=\left[\frac{\partial}{\partial\bm{\mu}_{A}% }i\bm{\sigma}\left(\bm{p}\mathcal{K}_{\text{at}}\times\bm{A}+\bm{A}\mathcal{K}% _{\text{at}}\times\bm{p}\right)\right]_{j}\\ =\alpha^{2}\left[\sigma_{j}\bm{\nabla}\left(\mathcal{K}\frac{\bm{r}_{A}}{r_{A}% ^{3}}\right)-\bm{\sigma}\nabla_{j}\left(\mathcal{K}\frac{\bm{r}_{A}}{r_{A}^{3}% }\right)\right]

(3.246)

in the $j$ direction. The corresponding matrix elements are

\begin{split}\Braket{\phi_{m}}{H_{j}^{\text{ZORA,ZSO}\bm{\mu}_{A}}}{\phi_{n}}=% &-\alpha^{2}\sigma_{j}\int\mathcal{K}_{\text{at}(n)}\frac{\bm{r}_{A}}{r_{A}^{3% }}\left[(\bm{\nabla}\phi_{m}(\bm{r}))\phi_{n}(\bm{r})+\phi_{m}(\bm{r})\bm{% \nabla}\phi_{n}(\bm{r})\right]\mathrm{d}\,\bm{r}\\ &+\alpha^{2}\bm{\sigma}\int\mathcal{K}_{\text{at}(n)}\frac{\bm{r}_{A}}{r_{A}^{% 3}}\left[(\nabla_{j}\phi_{m}(\bm{r}))\phi_{n}(\bm{r})+\phi_{m}(\bm{r})\nabla_{% j}\phi_{n}(\bm{r})\right]\mathrm{d}\,\bm{r}.\end{split}

(3.247)

If Eq. (3.247) were evaluated directly, the spin-orbitals would need to possess two independent complex components. However, we can decompose the ZSO operator as follows:

\begin{split}H_{i}^{\text{ZORA,ZSO}\bm{\mu}_{A}}=&\alpha^{2}\sum_{j}\sigma_{j}% H_{ij}^{\text{ZORA,ZSO}\bm{\mu}_{A}},\\ H_{ij}^{\text{ZORA,ZSO}\bm{\mu}_{A}}=&\delta_{ij}\bm{\nabla}\left(\mathcal{K}% \frac{\bm{r}_{A}}{r_{A}^{3}}\right)-\nabla_{i}\left(\mathcal{K}\frac{r_{Aj}}{r% _{A}^{3}}\right)\end{split}

(3.248)

Each term $\sigma_{j}H_{ij}^{\text{ZORA,ZSO}\bm{\mu}_{A}}$ is calculated separately by taking the spin quantization axis to be in the $i$ direction in each perturbation calculation. This makes the ZSO operator diagonal in spin space and allows us to use spin-free real orbitals to represent the perturbations. After integrating by parts, each of those terms is expressed as

\begin{split}\delta_{ij}&\braket{\phi_{m}}{\bm{\nabla}\left(\mathcal{K}\frac{% \bm{r}_{A}}{r_{A}^{3}}\right)}{\phi_{n}}+\braket{\phi_{m}}{-\nabla_{i}\left(% \mathcal{K}\frac{r_{Aj}}{r_{A}^{3}}\right)}{\phi_{n}}=\\ &-\int\mathcal{K}_{\text{at}(n)}\frac{\bm{r}_{A}}{r_{A}^{3}}\left[(\bm{\nabla}% \phi_{m}(\bm{r}))\phi_{n}(\bm{r})+\phi_{m}(\bm{r})\bm{\nabla}\phi_{n}(\bm{r})% \right]\mathrm{d}\,\bm{r}\\ &+\int\mathcal{K}_{\text{at}(n)}\frac{r_{Aj}}{r_{A}^{3}}\left[(\nabla_{i}\phi_% {m}(\bm{r}))\phi_{n}(\bm{r})+\phi_{m}(\bm{r})\nabla_{i}\phi_{n}(\bm{r})\right]% \mathrm{d}\,\bm{r},\end{split}

(3.249)

where both terms need to be explicitly symmetrized because of $\mathcal{K}_{\text{at}(n)}$ . The observable (the J-coupling) is then computed as

\begin{split}&\overleftrightarrow{J_{AB}^{\text{ZSO}}}=\frac{\gamma_{A}\gamma_% {B}}{2\pi}\sum_{i}^{\text{occ}}2\mathrm{Re}\braket{\psi_{i}^{\text{ZSO}\bm{\mu% }_{A}}}{\bm{H}^{\text{ZSO}\bm{\mu}_{A}}}{\psi_{i}}\\ &=\frac{\gamma_{A}\gamma_{B}}{2\pi}\sum_{i}^{\text{occ}}2\mathrm{Re}\big{(}% \braket{\psi_{xx}^{\text{ZSO}\bm{\mu}_{A}}}{H_{xx}^{\text{ZSO}\bm{\mu}_{A}}}{% \psi_{i}}+\braket{\psi_{yy}^{\text{ZSO}\bm{\mu}_{A}}}{H_{yy}^{\text{ZSO}\bm{% \mu}_{A}}}{\psi_{i}}+\braket{\psi_{zz}^{\text{ZSO}\bm{\mu}_{A}}}{H_{zz}^{\text% {ZSO}\bm{\mu}_{A}}}{\psi_{i}}\\ &+2\braket{\psi_{xy}^{\text{ZSO}\bm{\mu}_{A}}}{H_{xy}^{\text{ZSO}\bm{\mu}_{A}}% }{\psi_{i}}+2\braket{\psi_{xz}^{\text{ZSO}\bm{\mu}_{A}}}{H_{xz}^{\text{ZSO}\bm% {\mu}_{A}}}{\psi_{i}}+2\braket{\psi_{yz}^{\text{ZSO}\bm{\mu}_{A}}}{H_{yz}^{% \text{ZSO}\bm{\mu}_{A}}}{\psi_{i}}\big{)}.\end{split}

(3.250)

Equation (3.250) along with the PSO and DSO terms forms the complete ZORA J-coupling.

Following R. Bouten et al., J. Phys. Chem. A 104, 5600 (2000), the ZORA shielding integrals are

\begin{split}h_{mn}^{\text{ZORA,GIAO}\bm{B}\bm{\mu}_{A}}=&\frac{\alpha^{2}}{2}% \Braket{\phi_{m}}{\mathcal{K}_{\text{at}(n)}\frac{(\bm{R}_{mn}\times\bm{r})(% \bm{r}_{A}\times\bm{\nabla})^{T}+\bm{r}_{A}\bm{r}_{n}-\bm{r}_{A}\bm{r}_{n}^{T}% }{r_{A}^{3}}}{\phi_{n}},\\ h_{mn}^{\text{ZORA,GIAO}\bm{B}}=&\frac{i}{2}\Braket{\phi_{m}}{(\bm{R}_{mn}% \times\bm{r})H_{\text{LDA}}}{\phi_{n}}+\Braket{\phi_{m}}{V_{\text{GGA}}^{\bm{B% }}}{\phi_{n}}\\ &-\frac{i}{2}\Braket{\phi_{m}}{\mathcal{K}_{\text{at}(n)}\bm{r}_{n}\times\bm{% \nabla}}{\phi_{n}},\\ h_{mn}^{\text{ZORA,GIAO}\bm{\mu}_{A}}=&-i\alpha^{2}\Braket{\phi_{m}}{\mathcal{% K}_{\text{at}(n)}\frac{\bm{r}_{A}\times\bm{\nabla}}{r_{A}^{3}}}{\phi_{n}}.\end% {split}

(3.251)

Dipolar couplings

The direct spin-spin coupling tensor, $\overleftrightarrow{D_{AB}}$ , also called the dipolar coupling tensor, is expressed as

\overleftrightarrow{D_{AB}}=\frac{\hbar\alpha^{2}}{2\pi}\gamma_{A}\gamma_{B}% \left(\frac{1}{R_{AB}^{3}}-\frac{3\bm{R}_{AB}\bm{R}_{AB}^{T}}{R_{AB}^{5}}% \right),

(3.252)

where $\bm{R}_{AB}$ is the vector connecting atoms $A$ and $B$ . Because it does not depend on the electronic structure, its computation takes no time and it is always automatically included whenever a J-coupling calculation is requested

Current limitations

•

Periodic systems not supported;
•

Only the nonrelativistic formalism has been implemented (scalarrelativistic under development);
•

Hybrid functionals not supported for shieldings and magnetizability;
•

Meta functionals not supported.

For developers

Most of the magnetic response (MR) related source code resides in the directory
MagneticResponse. The parent subroutine for doing MR calculations is MR_main. It sets up the environment depending on user input, calls MR_core which performs the actual calculations, and prints various information including the results of the calculation. MR_core is called separately for every type term that needs to be computed. If the first-order wavefunctions are required, a DFPT cycle is performed first. If a diamagnetic property is calculated, the DFPT part is skipped.

The module integration provides a subroutine that is used for all integrations. The only exception is the nonrelativistic FC term, which needs a completely separate subroutine because of the delta function. The bodies of most functions that are integrated over real space are found in integrands.f90. The result of each integration is written into a 2D block cyclic matrix (the matrix_BC_out argument of the integrate). Outside the integration subroutine, all linear algebra involving large arrays is done with Scalapack routines.

Example input

The following is a minimal example for running a magnetic response calculation.

control.in:

  # H2O
  xc                 pw-lda
  magnetic_response  # Default is to calculate all magnetic
                     # response quantities
  # Basis sets
  ...

geometry.in:

  atom  0.00  0.00  0.00  O
  magnetic_response
  atom -0.96  0.00  0.00  H
  magnetic_response
  atom  0.32 -0.90  0.00  H
  magnetic_response

Notes on performance

For best performance, it is recommended to always include the following flags in control.in:

  load_balancing          .true.
  use_local_index         .true.
  collect_eigenvectors    .false.

In the timings section, the following terms are shown:

•

“Fermi contact”, “paramagnetic spin-orbit”, “spin-dipole”, “paramagnetic shielding”, “paramagnetic magnetizability” — these are the integration times of the non-self-consistent parts of the Hamiltonian [ $H_{\text{NSC}}^{(1)}$ in Eq. (3.238a)].
•

“diamagnetic spin-orbit”, “diamagnetic shielding”, “diamagnetic magnetizability” — integration of a diamagnetic magnetic property.
•

“First-order xc” — integration of the first-order response of the xc-potential.
•

“First-order density” — constructing the first-order electron density from the first-order density matrix (this is an integral over real space).
•

“First-order Hartree” — integration of the first-order response of the Hartree potential.
•

“First-order Ha update” — constructing the first-order Hartree potential from the first-order density. It mainly times calls to update_hartree_potential_p1 and sum_up_whole_potential_p1.
•

“Mat-mat multiplications” — total wall time spent on pdgemm and pzgemm.
•

“packed <-> block cyclic” — total wall time spent on subroutines that do conversion from one matrix type to another. These are packed_to_block_cyclic and block_cyclic_to_packed.
•

“Total” — total wall time spent on the magnetic response calculation. This is higher than the sum of the above, because it contains additional overhead from stuff that is in between the subroutines that are timed.
•

“Total minus individual terms” — Difference between the total wall time and the sum of individual terms. If the total time is much higher than the sum of the individual terms then there is an unexpected bottleneck somewhere in the code that does not scale. Note that the difference is between the total time of the slowest task and a sum where each term individually corresponds to that of the slowest task. Thus, do not be alarmed if this number becomes slightly negative (unlikely, but possible depending on how the load is balanced).

Maximum and minimum cpu times are shown where possible. In case the subroutine to be timed contains and explicit or implicit MPI barrier (e.g., pdgemm), individual cpu times are not meaningful and only the maximum wall time is shown.

Tags for general section of control.in:

Tag: magnetic_response(control.in)

Usage: magnetic_response options
Purpose: Primary keyword for doing magnetic response calculations. This is the minimum that is required in control.in and by default leads to the full calculation of J-couplings, shieldings, and the magnetizability.
This keyword can be followed by a number of options:

•

J_coupling or J or j — Calculate J-couplings only.
•

shielding or s — Calculate the shieldings only.
•

magnet or m — Calculate the magnetizability only.
•

fc — Calculate the FC contribution to J-couplings only.
•

po — Calculate the PSO contribution to J-couplings only.
•

sd — Calculate the SD contribution to J-couplings only.
•

do — Calculate the DSO contribution to J-couplings only.
•

shielding_p or s_p — Calculate the paramagnetic contribution to the shieldings only.
•

shield_d or s_d — Calculate the diamagnetic contribution to the shieldings only.
•

magnet_p or m_p — Calculate the paramagnetic contribution to the magnetizability only.
•

magnet_d or m_d — Calculate the diamagnetic contribution to the magnetizability only.
•

no_giao — Calculate the shieldings and magnetizability using the standard formalism without GIAOs (default: false, i.e., with GIAOs).
•

full — Calculate the full tensors (default: only diagonal elements).

Comment: any combination of options works. For example,
magnetic_response fc s po full
calculates the FC and PSO terms of J-couplings and the shielding tensors including off-diagonal elements. The default is equivalent to
magnetic_response J_coupling shielding magnet

Tag: dfpt_accuracy_n1(control.in)

Usage: dfpt_accuracy_n1 value
Purpose: Convergence criterion for the DFPT self-consistency cycle, based on the RMS change in the first-order density matrix. Specifically, the unmixed output density matrix is checked against the input density matrix corresponding to the same iteration. The RMS value is calculated over all directions and spins that are being processed.
value is a real positive number (in electrons). Default: 1d-9.

Tag: dfpt_iter_limit(control.in)

Usage: dfpt_iter_limit number
Purpose: Maximum number of DFPT cycles.
number is an integer number. Default: 40.

Tag: dfpt_linear_mix_param(control.in)

Usage: dfpt_linear_mix_param value
Purpose: Parameter for linear mixing of first-order electron density. Used in Pulay and simple linear mixing.
value is a real number between 0 and 1. Default: 1.0

Tag: dfpt_pulay_steps(control.in)

Usage: dfpt_pulay_steps number
Purpose: Number of steps kept in memory for Pulay mixing. Value of 1 corresponds to simple linear mixing.
number is a positive integer number. Default: 8

Tag: mr_gauge_origin(control.in)

Usage: mr_gauge_origin x y z
Purpose: Gauge origin for the calculation of the shieldings or the magnetizability. This has no effect with GIAOs.
x y z are real numbers (in Å) that specify the position of the gauge origin. Default: the center of mass.
Comment: The center of mass is based on the natural abundance of elements. If this is not desirable, the user can set the gauge origin manually.

Tag: output_sxml(control.in)

Usage: output_sxml [name]
Purpose: If present, the spin system and the results of calculations are printed into an sxml file [Biternas et al., J. Mag. Res. 240, 124 (2014)], which can serve as input to other software, e.g., Spinach [spindynamics.org].
The name of the output file is name if present. Otherwise, the default name is aims.xsml.

Tag: mr_experimental(control.in)

Usage: mr_experimental
Purpose: Overrides any safety checks that would cause the run to stop. This allows the user to test features that are still in development and even avoid simple sanity checks such calculating the shieldings without specifying any atoms in geometry.in.

Tags for geometry.in:

Tag: magnetic_response(geometry.in)

Usage: magnetic_response
Purpose: Includes the current atom in the magnetic response calculations. If only the magnetizability is required, this keyword need not be used in geometry.in. Otherwise, the calculation of the shieldings or J-couplings is aborted if no atoms are flagged for MR calculations in geometry.in

Tag: magnetic_moment(geometry.in)

Usage: magnetic_moment value
Purpose: Overrides the default magnetic moment for the given atom. The default values (in units of the nuclear magneton $\mu_{N}$ ) can be found in MagneticResponse/MR_nuclear_data.f90. In case of J-couplings, the isotopes used are also printed in the output.
value, a real number, is the magnetic moment in atomic units (-5.157d-4 for O-17).

Tag: nuclear_spin(geometry.in)

Usage: nuclear_spin value
Purpose: Overrides the default nuclear spin for the given atom. The default values can be found in MagneticResponse/MR_nuclear_data.f90 and are also printed in the output for J-coupling calculations.
value, a real number, is the nuclear spin (2.5 for O-17).

Tag: isotope(geometry.in)

Usage: isotope number
Purpose: Overrides the default isotope mass number for the given atom. For more flexibility, the magnetic moment and spin can be specified separately with the above keywords. The default isotopes numbers can be found in MagneticResponse/MR_nuclear_data.f90.
number, a positive integer, is the mass of the given atom (17 for O-17).