EPR parameters from first principles

The spin Hamiltonian

The theoretical reconstruction of the EPR spectrum of a radical uses an effective operator, the spin Hamiltonian, rather than the full molecular Hamiltonian. The spin Hamiltonian acts only on the spin variables of a fictitious electron spin \(\hat{\mathbf{S}}\) and of the nuclear spins \(\hat{\mathbf{I}}^A\). A handful of parameters contain every reference to the spatial part of the many-electron wavefunction, and thus to the electronic and geometric structure [Neese2009]. The spin Hamiltonian reads

(1)\[\hat{H}_{\mathrm{SH}} = \hat{\mathbf{S}}\,\mathbf{D}\,\hat{\mathbf{S}} + \mu_\mathrm{B}\,\mathbf{B}\,\mathbf{g}\,\hat{\mathbf{S}} + \sum_A \left[ \hat{\mathbf{S}}\,\mathbf{A}^A\,\hat{\mathbf{I}}^A + \mu_\mathrm{N}\,\mathbf{B}\,\mathbf{g}_\mathrm{N}^A\,\hat{\mathbf{I}}^A + \hat{\mathbf{I}}^A\,\mathbf{Q}^A\,\hat{\mathbf{I}}^A \right] + \sum_{A<B} \hat{\mathbf{I}}^A\,\mathbf{J}^{AB}\,\hat{\mathbf{I}}^B ,\]

where the indices \(A\) and \(B\) label the magnetic nuclei, \(\mathbf{B}\) is the magnetic flux density, and \(\mu_\mathrm{B}\) and \(\mu_\mathrm{N}\) are the Bohr and nuclear magnetons. The tensors carry the chemistry of the molecule:

\(\mathbf{D}\), the zero-field splitting,
\(\mathbf{g}\), the electron Zeeman interaction,
\(\mathbf{A}^A\), the hyperfine coupling to nucleus \(A\),
\(\mathbf{Q}^A\), its quadrupole coupling, and
\(\mathbf{J}^{AB}\), the coupling between two nuclear spins.

The simplicity of this approach makes the spin Hamiltonian so useful. It acts on the product basis \(|S M_S\rangle |I_A M_{I_A}\rangle \cdots\), which for a few spins stays small enough for exact diagonalization. A set of parameters therefore yields a simulated spectrum at modest computational cost [Neese2009].

Only part of Eq. (1) matters for the radicals that EPRoducer targets, i.e. organic species with a single unpaired electron, \(S = 1/2\):

The zero-field term vanishes, because \(\mathbf{D}\) only splits the sublevels of spins with \(S > 1/2\).
The nuclear-nuclear coupling \(\mathbf{J}^{AB}\) stays far below the resolution of a continuous-wave experiment.
EPR simulations neglect the chemical shift and its anisotropy in the nuclear Zeeman term altogether [Stoll2006].

The remaining terms are the electron Zeeman interaction with the parameter \(\mathbf{g}\), the hyperfine coupling with \(\mathbf{A}^A\), and the quadrupole coupling with \(\mathbf{Q}^A\). The latter is only relevant for nuclei with \(I > 1/2\), such as ¹⁴N. The rest of this page describes how to calculate those parameters for a given molecule.

Linear response theory

Linear response theory connects the parameters of the spin Hamiltonian with the electronic structure of the radical, and it does so without the sum over excited states that ordinary perturbation theory requires [McWeeny1969]. Perturbation theory offers the most direct route to that connection, by treating the magnetic interactions as a perturbation of the Born-Oppenheimer problem, whose ground state of total spin \(S\) carries \(2S+1\) degenerate magnetic sublevels [Neese2009] [Neese2011]. Although this sum-over-states picture shows the physics with unmatched clarity, it is a poor recipe for real calculations. Its second-order terms require an infinite sum over excited many-electron states, whereas a practical calculation reaches at most a few dozen, without any general guarantee that the truncated sum converges [Neese2009]. Linear response, also known in the literature as analytic derivative, double-perturbation theory, or coupled-perturbed self-consistent field (CP-SCF) method, avoids this problem. It implicitly contains an untruncated sum over excited states, described at the same level of theory as the ground state, and therefore introduces no arbitrary truncation into the result [Neese2009].

Linear response treats every spin-Hamiltonian parameter as a second derivative of the total ground-state energy. Each parameter is bilinear in two perturbations, such as a component of the magnetic field, a nuclear magnetic moment or an electron magnetic moment, and therefore equals the second derivative of the energy with respect to a pair of perturbation parameters \(\lambda\) and \(\kappa\). For a self-consistent field ground state, such as the outcome of a DFT calculation, and with basis functions that do not depend on the perturbations, that derivative is

(2)\[\left.\frac{\partial^2 E}{\partial\lambda\,\partial\kappa}\right|_{0} = \sum_{\mu\nu} P_{\mu\nu} \langle \varphi_\mu | \hat{h}^{\lambda\kappa} | \varphi_\nu \rangle + \sum_{\mu\nu} \frac{\partial P_{\mu\nu}}{\partial\kappa} \langle \varphi_\mu | \hat{h}^{\lambda} | \varphi_\nu \rangle ,\]

where \(\varphi_\nu\) are basis functions, \(\hat{h}^{\lambda}\) the derivative of the one-electron operator with respect to \(\lambda\), and \(P_{\mu\nu}\) the total ground-state density matrix. The latter contains the coefficients of the expansion of the electronic density in the basis set \(\{\varphi\}\),

\[\rho(\mathbf{r}) = \sum_{\mu\nu} P_{\mu\nu}\, \varphi_\mu(\mathbf{r})\, \varphi_\nu(\mathbf{r}) .\]

By construction, \(P_{\mu\nu}\) is the sum of the spin up (\(\uparrow\)) and spin down (\(\downarrow\)) contributions [Kohn1996],

\[P_{\mu\nu} = P^\uparrow_{\mu\nu} + P^\downarrow_{\mu\nu} = \sum_i c^\uparrow_{\mu i}c^\uparrow_{\nu i} + \sum_i c^\downarrow_{\mu i}c^\downarrow_{\nu i} ,\]

where the index \(i\) runs over the occupied molecular orbitals. In the context of EPR, the spin-density matrix is also important,

\[P^{\uparrow-\downarrow}_{\mu\nu} = P^\uparrow_{\mu\nu} - P^\downarrow_{\mu\nu} .\]

The two terms of Eq. (2) replace the first- and second-order terms of the perturbation sum exactly. The first needs only the unperturbed density, whereas the second needs the response of the density to a perturbation [Neese2009].

Obtaining that response is the actual computational work. Differentiating the SCF equations and expanding the perturbed orbitals in the unperturbed ones gives a linear system, \(\mathbf{A}\,\mathbf{U}^{\lambda} = \mathbf{V}^{\lambda}\), for the occupied-virtual mixing coefficients \(\mathbf{U}^{\lambda}\), in which the perturbation enters only through \(\mathbf{V}^{\lambda}\). The structure of \(\mathbf{A}\) depends on the character of the perturbation:

Electric field-like perturbations produce real, symmetric coefficients and an \(\mathbf{A}\) that carries the Coulomb term, the exchange-correlation kernel and the Hartree-Fock exchange.
Magnetic field-like perturbations, the relevant class for EPR, produce purely imaginary, Hermitian coefficients. In this case, every contribution of a local potential cancels. Without Hartree-Fock exchange, \(\mathbf{A}\) becomes diagonal and the system solves trivially. The result is a sum-over-orbitals expression that superficially resembles the sum-over-states result but has an entirely different origin.

Hybrid functionals reinstate the coupling and leave a linear system of dimension \(N_{\mathrm{occ}} \times N_{\mathrm{virt}}\), which can reach hundreds of thousands of coefficients. Iterative solvers handle it without ever building \(\mathbf{A}\) explicitly [Neese2009].

The g-tensor

With the response density available, and considering \(\lambda=B_k\) and \(\kappa=S_l\) in Eq. (2), every spin-Hamiltonian parameter becomes a closed expression over the basis functions. The \(g\)-tensor, in particular, separates into four contributions [Neese2009] [Neese2011],

(3)\[g_{kl} = g_\mathrm{e}\,\delta_{kl} + \Delta g^{\mathrm{RMC}}\,\delta_{kl} + \Delta g^{\mathrm{GC}}_{kl} + \Delta g^{\mathrm{OZ/SOC}}_{kl} ,\]

where \(k, l \in \{x, y, z\}\) label Cartesian components and, as in all expressions of this page, atomic units apply. The first term is the free-electron value \(g_\mathrm{e} = 2.002319\). It usually dominates, but as a natural constant it carries no information about the geometric or electronic structure of the molecule. That information resides entirely in the three remaining terms.

Relativistic mass correction (RMC). This term is isotropic,

(4)\[\Delta g^{\mathrm{RMC}} = -\frac{\alpha^2}{S} \sum_{\mu\nu} P^{\uparrow-\downarrow}_{\mu\nu}\, \langle \varphi_\mu | \hat{T} | \varphi_\nu \rangle ,\]

where \(\alpha \approx 1/137\) is the fine-structure constant and \(\hat{T}\) the kinetic energy operator. It weighs the kinetic energy of the unpaired spin density, and it corrects for the relativistic increase of the electron mass at high velocity.

Gauge correction (GC). This is a first-order correction to the spin-orbit coupling in the presence of the external field,

(5)\[\Delta g^{\mathrm{GC}}_{kl} = \frac{1}{2S} \sum_{\mu\nu} P^{\uparrow-\downarrow}_{\mu\nu} \left\langle \varphi_\mu \left| \sum_A \xi(\hat r_A) \left[ \delta_{kl}\,\hat{\mathbf{r}}_A \cdot \hat{\mathbf{r}} - \hat r_{A,k}\, \hat r_l \right] \right| \varphi_\nu \right\rangle ,\]

where the sum runs over all nuclei \(A\), \(\hat{\mathbf{r}}\) is the position operator and \(\hat{\mathbf{r}}_A = \hat{\mathbf{r}} - \mathbf{R}_A\) the position operator relative to nucleus \(A\) at \(\mathbf{R}_A\). The scalar operator \(\xi(\hat r_A)\) is the radial factor of the spin-orbit coupling near nucleus \(A\), which in the one-electron form of the operator reads \(\xi(\hat r_A) = \alpha^2 Z_A / (2 \hat r_A^3)\), with \(Z_A\) the nuclear charge [Neese2009]. Because the RMC and GC terms are both expectation values over the unperturbed spin density, they need no response calculation. Both usually stay small.

Orbital Zeeman / spin-orbit coupling cross term (OZ/SOC). The dominant contribution to \(\Delta\mathbf{g}\) comes from the cross term between the spin-orbit coupling and the orbital Zeeman interaction,

(6)\[\Delta g^{\mathrm{OZ/SOC}}_{kl} = \frac{1}{2S} \sum_{\mu\nu} \frac{\partial P^{\uparrow-\downarrow}_{\mu\nu}}{\partial B_k} \langle \varphi_\mu | \hat{z}^{\mathrm{SOMF}}_{l} | \varphi_\nu \rangle .\]

Here, \(\partial P^{\uparrow-\downarrow}_{\mu\nu}/\partial B_k\) is the response of the spin-density matrix to the component \(k\) of the external magnetic field. The operator \(\hat{z}^{\mathrm{SOMF}}_l\) is the component \(l\) of the spatial part of the spin-orbit mean-field (SOMF) operator. This effective one-electron operator combines the spin-orbit coupling of each electron in the field of the nuclei with a mean-field average of the two-electron spin-same-orbit and spin-other-orbit terms, which screen the nuclear contribution. It recovers about 99% of the full two-electron spin-orbit coupling [Neese2009] [Neese2005]. Physically, the external field induces orbital angular momentum in the electron cloud, and the spin-orbit coupling transfers it to the spin of the unpaired electron. The term is therefore a genuine response property that grows with the spin-orbit coupling of the atoms that carry spin density.

Gauge origin

Every reported \(g\) value seems to depend on the choice of the gauge origin, but there are methods to remove this dependence. The position operator \(\hat{\mathbf{r}}\) in Eq. (5) and the angular momentum \(\hat{\mathbf{l}}\), implicit in the response density of Eq. (6), refer to the origin of the coordinate system. The computed tensor therefore depends on where that origin sits. This dependence disappears only in the basis-set limit, which no practical calculation reaches.

Gauge-including atomic orbitals (GIAOs) remove the artefact rigorously, at extra computational cost [Ditchfield1974].
Placing the origin at the centre of electronic charge (CoEC) keeps results comparable while introducing an error well below those of the functional, the basis set or the molecular model, because the residual gauge error in EPR is small [Neese2009].

Two double-zeta basis sets are common for these calculations. The cc-pVDZ set belongs to Dunning’s family of correlation-consistent basis sets [Dunning1989]. EPR-II adds a single set of polarization functions and an enhanced s part, optimized for hyperfine couplings with DFT and particularly B3LYP [Barone2002].

In EPRoducer

GIAOs remain computationally affordable for the molecules that EPRoducer targets, so gval_workflow computes the \(g\) tensor with GIAOs and B3LYP/EPR-II by default (flag -x to change the functional and basis set). The principal values are read from the ORCA output with eproducer.postpro.extract_gvals().

The hyperfine coupling

The hyperfine tensor of nucleus \(A\) consists of three parts, the Fermi contact term (c), the spin-dipolar term (d) and the spin-orbit correction (SO), such that \(\mathbf{A}^A = \mathbf{A}^{A;\mathrm{c}} + \mathbf{A}^{A;\mathrm{d}} + \mathbf{A}^{A;\mathrm{SO}}\) [Neese2009] [Neese2011]. All three share the prefactor \(P_A = g_\mathrm{e}\,g_\mathrm{N}^A\,\mu_\mathrm{B}\,\mu_\mathrm{N}\).

Fermi contact term. It is isotropic,

(7)\[A^{A;\mathrm{c}}_{kl} = \delta_{kl}\,\frac{8\pi}{3}\,\frac{P_A}{2S}\, \rho^{\uparrow-\downarrow}(\mathbf{R}_A) ,\]

and depends only on the spin density at the position of the nucleus,

\[\rho^{\uparrow-\downarrow}(\mathbf{R}_A) = \sum_{\mu\nu} P^{\uparrow-\downarrow}_{\mu\nu}\, \varphi_\mu(\mathbf{R}_A)\,\varphi_\nu(\mathbf{R}_A) .\]

Since only s-type functions remain finite at a nucleus, this term probes the s character of the spin density, including the spin polarization of the core shells. That sensitivity to a single point explains why the isotropic coupling, and not the anisotropic part, drives the basis-set requirements (see Practical recommendations).

Spin-dipolar term. It describes the through-space dipolar interaction between the electron spin and the magnetic moment of the nucleus,

(8)\[A^{A;\mathrm{d}}_{kl} = \frac{P_A}{2S} \sum_{\mu\nu} P^{\uparrow-\downarrow}_{\mu\nu} \left\langle \varphi_\mu \left| r_A^{-5} \left( r_A^2\,\delta_{kl} - 3\,r_{A,k}\,r_{A,l} \right) \right| \varphi_\nu \right\rangle .\]

Since the operator in this equation has zero trace, this term contributes nothing to the isotropic coupling and alone forms the anisotropic part of the tensor. Its integrand decays as \(r_A^{-3}\). The term therefore samples the spin density in the vicinity of nucleus \(A\) and reports on its shape. The same integrals, contracted with the total electron density instead of the spin density, give the electric field gradient at the nucleus [Neese2009].

Spin-orbit correction. The Fermi contact and the spin-dipolar terms are both expectation values over the unperturbed spin density, whereas the spin-orbit correction is a response property,

(9)\[A^{A;\mathrm{SO}}_{kl} = -\frac{P_A}{S} \sum_{\mu\nu} \frac{\partial P^{\uparrow-\downarrow}_{\mu\nu}}{\partial \hat{I}^{(A)}_k} \langle \varphi_\mu | \hat{z}^{\mathrm{SOMF}}_l | \varphi_\nu \rangle ,\]

in which \(\partial P^{\uparrow-\downarrow}_{\mu\nu}/\partial \hat{I}^{(A)}_k\) is the response of the spin-density matrix to component \(k\) of the nuclear spin of \(A\). The perturbing operator is now the nuclear-orbit interaction, the coupling between the magnetic moment of the nucleus and the orbital angular momentum of the electrons around it, which decays as \(r_A^{-3}\) [Neese2009]. The term is thus the hyperfine analogue of \(\Delta g^{\mathrm{OZ/SOC}}\), with the nuclear moment in place of the external field. It requires its own set of coupled-perturbed equations for every nucleus of interest, a cost that restricts it in practice to a few selected heavier nuclei. For light nuclei the correction is usually negligible [Neese2009], which simplifies the treatment of organic radicals, since they contain mainly H, C, N and O.

The quadrupole coupling

The quadrupole coupling follows from the electric field gradient (EFG) that the electrons create at a nucleus with \(I \geq 1\) (nitrogen, for example), and it is the only spin-Hamiltonian parameter that depends on the total electron density rather than on the spin density [Neese2009]. The EFG tensor at nucleus \(A\) uses the same integrals as the spin-dipolar term of Eq. (8), but contracts them with the total density matrix,

(10)\[V^{(A)}_{kl} = \sum_{\mu\nu} P_{\mu\nu} \left\langle \varphi_\mu \left| r_A^{-5} \left( r_A^2\,\delta_{kl} - 3\,r_{A,k}\,r_{A,l} \right) \right| \varphi_\nu \right\rangle .\]

The numerically largest eigenvalue of this tensor, \(V_{\mathrm{max}}\), defines the field gradient \(q = V_{\mathrm{max}}\) and, together with the quadrupole moment \(Q\) of the nucleus, the quadrupole splitting parameter \(e^2qQ\). In the eigensystem of the electric field gradient, the quadrupole term of Eq. (1) reads

(11)\[\begin{split}\hat{H}_\mathrm{Q} = \hat{\mathbf{I}}^A\,\mathbf{Q}^A\,\hat{\mathbf{I}}^A = \frac{e^2qQ}{4I(2I-1)}\,\hat{\mathbf{I}}^A \begin{pmatrix} -(1-\eta) & 0 & 0 \\ 0 & -(1+\eta) & 0 \\ 0 & 0 & 2 \end{pmatrix} \hat{\mathbf{I}}^A ,\end{split}\]

where the asymmetry parameter

(12)\[\eta = \frac{\left| V_{\mathrm{mid}} - V_{\mathrm{min}} \right|}{\left| V_{\mathrm{max}} \right|}\]

measures how far the field gradient departs from axial symmetry, with \(V_{\mathrm{mid}}\) and \(V_{\mathrm{min}}\) the two remaining eigenvalues [Neese2009]. The factor \(I(2I-1)\) in Eq. (11) reflects that only nuclei with \(I \geq 1\) carry a quadrupole moment. In organic molecules, ¹⁴N (\(I = 1\)) is therefore the only common quadrupolar nucleus, since ¹H has \(I = 1/2\) and the most abundant isotopes of carbon and oxygen, ¹²C and ¹⁶O, carry no nuclear spin.

Practical recommendations

Predicting the parameters of Eq. (1) is a well-defined but not automatic procedure. The following recommendations come from the advice of Frank Neese [Neese2009] [Neese2011] and from our own experience with EPR calculations. The defaults of EPRoducer follow them.

Molecular model. The choice of the molecular model matters more than any later technical decision. A large, realistic model gives the most accurate parameters, yet its low symmetry and complexity can leave the computed numbers as opaque as fitted ones. A deliberately small model, studied systematically, instead reveals the structure-spectrum relationships, and it often improves the subsequent calculation on the large system.

Geometries. Optimized geometries serve better than experimental ones unless a strong reason argues otherwise. They are reproducible, they do not inherit the resolution limits of a diffraction experiment, and they place hydrogens far more reliably. Even good protein crystal structures carry bond-length uncertainties of about 0.1 Å, which quantum chemical parameter calculations do not tolerate. When a truncated model must respect constraints from its surroundings, freezing a few coordinates helps, whereas a QM/MM treatment of the full protein handles both the geometric constraints and the long-range electrostatics. Where dynamics matters, averaging over snapshots of a molecular dynamics trajectory yields the mean parameters and, through their statistics, an estimate of the inhomogeneous line width. In the same line, for a complete study of a molecule to compare with experiments, we recommend including as many metastable configurations as possible.

Functionals. B3LYP is the best-calibrated functional for predicting EPR properties, with PBE0 as a comparable alternative. BP86 is the tested non-hybrid option. Since it lacks Hartree-Fock exchange, it runs about an order of magnitude faster, which pays off during geometry optimization. Tuning the fraction of Hartree-Fock exchange to match an experiment does not justify the computer time, because the optimal fraction varies from system to system and from property to property. Such a fit consequently adds no physical insight.

Basis sets. Double-\(\zeta\) basis sets with polarization functions (SVP, 6-31G*) are the minimum for structures, and polarized triple-\(\zeta\) sets (TZVP or cc-pVTZ) come close to the basis-set limit. Among the EPR parameters, only the isotropic hyperfine coupling is genuinely basis-set critical, because Eq. (7) samples the spin density at the nucleus and therefore needs extra flexibility in the core region. EPR-II and EPR-III cover the first-row elements, IGLO-II and IGLO-III most of the first and second rows, and CP(PPP) the first-row transition metals. Beyond the first transition row, relativistic effects demand an all-electron treatment of hyperfine and quadrupole couplings.

Expected accuracy. Computed parameters and measured spectra rarely match one to one, because both carry errors. The calculations rest on the approximations of this page. One of the largest errors falls on the largest contribution to the \(g\) values, \(\Delta g^{\mathrm{OZ/SOC}}\). Without Hartree-Fock exchange, this term reduces to a sum over pairs of occupied orbitals \(i\) and virtual orbitals \(a\) (see Linear response theory), in which every term is \(\propto 1/(\varepsilon_a - \varepsilon_i)\), with \(\varepsilon\) the orbital energies. Hybrid functionals such as B3LYP couple the pairs but keep this dependence. A small gap between the frontier orbitals therefore makes the term large, and with it any error of the functional in the orbital energies. The experiments add their own uncertainties. Deviations of 10–20% between computed and measured parameters are acceptable when fitting spectra [Neese2011]. For the \(g\) tensor, this tolerance refers to the shift \(\Delta\mathbf{g}\), since \(g\) itself never strays far from \(g_\mathrm{e}\).

Principal values. We report the three principal values of the \(g\) tensor rather than its nine components, which has two consequences. First, since the principal values carry no natural order, we sort them by convention as \(g_1 > g_2 > g_3\), also written \(g_x > g_y > g_z\). Second, unlike the tensor components, they do not depend on the orientation of the molecule, although they still carry the gauge-origin dependence discussed in Gauge origin.

Orientations. The simulated spectra assume that the molecules take every orientation with equal probability. If, for example, more molecules than average point one of their principal axes along the magnetic field, the feature of the spectrum at the corresponding \(g\) value gains intensity.

The next page explains how to compute spectra from these parameters and how the spectra respond when the parameters vary.