Electron paramagnetic resonance (EPR) spectra of molecular spin centers undergoing reorientational motion are commonly simulated using the stochastic Liouville equation (SLE) with a rigid-body hindered Brownian diffusion model. Current SLE theory applies to specific spin systems such as nitroxides and to high-symmetry orientational potentials. In this work, we extend the SLE theory to arbitrary spin systems with any number of spins and any type of spin Hamiltonian interaction term, as well as to arbitrarily complex orientational potentials. We also examine the limited accuracy of the frequency-to-field conversion used to obtain field-swept EPR spectra and present a more accurate approach. The extensions allow for the simulation of EPR spectra of all types of spin labels (nitroxides, copper2+, and gadolinium3+) attached to proteins in low-symmetry environments.

Continuous-wave electron paramagnetic resonance (CW EPR) spectra reveal important information about the structural and dynamic properties of paramagnetic spin centers. In particular, the spectral line shape can be highly sensitive to the nature and time scale of rotational dynamics of the spin center. An important application in this regard is the study of the dynamics of spin labels attached to soluble or membrane proteins [see Fig. 1(a)].

FIG. 1.

(a) The motion of a nitroxide radical tethered to a protein is constrained by neighboring amino acids (T4 lysozyme T115R1, pdb 2IGC36). (b) The hindered Brownian diffusion (HBD) model consists of a rigid body (representing the spin center) stochastically rotating in an external mean-field orientational potential energy surface (representing the interaction between the spin label and the protein or lipid nanoenvironment).

FIG. 1.

(a) The motion of a nitroxide radical tethered to a protein is constrained by neighboring amino acids (T4 lysozyme T115R1, pdb 2IGC36). (b) The hindered Brownian diffusion (HBD) model consists of a rigid body (representing the spin center) stochastically rotating in an external mean-field orientational potential energy surface (representing the interaction between the spin label and the protein or lipid nanoenvironment).

Close modal

On the simplest level, the time scale of rotational dynamics of a spin center is characterized by its rotational correlation time, τc, which is related to the rotational diffusion rate constant R by τc = 1/6R. The shape of the CW EPR spectrum depends on the relation of τc to the width of the rigid-limit spectrum, Δω. If τcΔω ≫ 1, then the rotational motion has very little effect on the spin dynamics, the spectrum resembles the one in the immobile limit (τc), and the rotational motion can be neglected for spectral simulations. If, on the other hand, τcΔω ≤ 1, then the motion is fast enough to mostly average out all anisotropies, and the observed spectrum is a sum of individual lines, similar to the one in the isotropic limit (τc = 0). In this fast-motion regime, spectra can be simulated by using an isotropic Hamiltonian and treating the rotational motion as a perturbation. The intermediate regime, where approximately 1 < τcΔω < 100, is called the slow-motion regime. In this regime, the rotational motion and the spin dynamics are strongly coupled, and the spectrum is sensitive to the details of the rotational motion. For simulating spectra in the slow-motion regime, the spin dynamics and the rotational dynamics have to be treated on an equal footing. For nitroxide radicals at X-band fields (∼0.34 T), Δω/2π ≈ 200 MHz so that the slow-motion regime is around τc ≈ 1–100 ns. This is the range often observed for nitroxides attached to proteins.

Several approaches for the simulation of slow-motion CW EPR spectra have been developed. They are based on motional models that range from full deterministic atomistic molecular dynamics (MD) to simple stochastic reorientation.1–5 A simple and, therefore, widely applied model is hindered Brownian diffusion (HBD).6 As illustrated in Fig. 1(b), the HBD model represents the tumbling spin center as a single rigid body undergoing Brownian rotational diffusion with an anisotropic rotational diffusion rate R. Besides the orientation Ω of the body, no other spatial degrees of freedom are dynamic in this model. All internal degrees of freedom are considered fixed. The nano-environment (such as a protein or lipid environment), also assumed fixed, hinders the rotational diffusion. Its effect is modeled with an orientation-dependent potential energy function U(Ω) that imposes an orientation-dependent torque onto the rigid body.

The motional model is combined with spin quantum dynamics to simulate the slow-motion CW EPR spectrum. For the HBD model, the most efficient and widely used approach is a sophisticated frequency-domain method based on the stochastic Liouville equation (SLE), pioneered by Kubo and co-workers.7–18 It employs highly efficient numerical methods and has the advantage of significantly lower computational cost compared to other HBD solvers that are based on stochastic diffusive or jump trajectories.

The SLE theory was originally developed to simulate CW EPR spectra for slow-tumbling nitroxide radicals in solution or liquid crystals.10,13 The expressions derived were particular to nitroxides, and the orientational potentials employed were of high symmetry. Since then, the theory has been expanded to include two magnetic nuclei19,20 and it has been extended and applied to rigid bisnitroxides.21–23 The slowly relaxing local structure (SRLS) model was developed to include additional rotational dynamics of the cage encompassing a spin label.18,24,25 A nuclear magnetic resonance (NMR)-focused program, Spinach, can solve the SLE for large spin systems in the absence of orientational potentials.26 Simulation of slow-motion spectra for high-spin systems (S > 1/2) has also been reported. For instance, slow-motion spectra of triplets (S = 1) in the absence of magnetic nuclei have been simulated using the SLE27 or a discrete-jump approach.28–30 The SLE has been employed to calculate CW EPR linewidths for Gd3+ centers (S = 7/2) undergoing unhindered rotational dynamics by including additional internal degrees of freedom such as vibration31 or pseudo-rotation.32 Lower-symmetry orientational potentials have been utilized in work on biaxial liquid crystal phases.33 Software codes stemming from several of these works are available.

Despite these advances, there exist no comprehensive methodology and software for calculating slow-motion spectra using the SLE approach without constraints on the constitution of the spin center and/or on the complexity of the environment hindering the reorientational motion (represented by the orientational potential). Such an extended SLE theory is needed due to the increased use of spin labels to study protein dynamics and the increased use of spin labels other than nitroxides in recent years.

In this paper, we present expressions that extend SLE theory to (a) orientational potentials of arbitrary complexity and to (b) spin centers of arbitrary composition, i.e., with any number of electron and nuclear spins and spin Hamiltonian interaction terms (Zeeman, hyperfine, zero field, exchange interactions, nuclear quadrupole, etc.). Furthermore, we examine the issue of calculating the field-swept CW EPR spectrum. Traditionally, SLE solvers calculate a frequency spectrum that is then converted to a field-swept spectrum using a first-order approximation. This method is inaccurate for systems with highly anisotropic g-tensors, and we present a more accurate way to simulate field-swept spectra.

In the following, Sec. II comprehensively lays out the SLE theory. Section III briefly summarizes implementation details. Section IV illustrates the extended scope of the theory with simulations of high-spin systems, low-symmetry potentials, and multinuclear systems. A few numerical aspects are discussed. Section V contains concluding discussions. All methods presented in this paper are implemented in the open-source software package EasySpin.5,34,35

In this section, we present the key expressions of SLE theory,13,16,17,37 including our extensions. We will use a series of space-fixed and body-fixed frames, which are shown in Fig. 2(a). The laboratory frame (L) is a space-fixed frame with its z axis along the external static magnetic field and its y axis along the oscillatory magnetic-field component of the microwave radiation. The potential frame (U) is a frame that is attached to the nano-environment hindering the motion of the spin center, e.g., a protein, a membrane bilayer, a polymer matrix, or a liquid-crystal phase. We limit ourselves to situations where proteins and membranes are immobile, i.e., where the potential frame is space-fixed. Often, the frame U is also referred to as the director frame. The potential frame will be our space-fixed reference frame. The mobile spin center is associated with a series of body-fixed frames that have time-dependent orientation: There is a molecular frame (M) associated with the molecular geometry. Each of the interaction tensors (g-tensors, hyperfine tensors, etc.) has a body-fixed eigenframe. The diffusion tensor (vide infra) is body-fixed as well, and the associated eigenframe is called the diffusion frame (D). It will serve as our body-fixed reference frame. To specify the relative orientation between any two frames, we use a triplet of Euler angles, indicated by Ω = (α, β, γ) and defined in Fig. 2(b). Some of the relative orientations used in the following are indicated in Fig. 2(a).

FIG. 2.

(a) Definition of space- and body-fixed frames and their relative orientations. (b) Definition of Euler angles ΩU→D = (α, β, γ) transforming frame U (xU, yU, zU) to frame D (xD, yD, zD). The intervals of definition are [0, 2π) for α and γ, and [0, π] for β.

FIG. 2.

(a) Definition of space- and body-fixed frames and their relative orientations. (b) Definition of Euler angles ΩU→D = (α, β, γ) transforming frame U (xU, yU, zU) to frame D (xD, yD, zD). The intervals of definition are [0, 2π) for α and γ, and [0, π] for β.

Close modal

The spin state of a tumbling spin center is described by the quantum spin density operator ρ(t). Its orientational state is described classically by the orientation Ω(t). Ω indicates the orientation of a body-fixed reference frame relative to a space-fixed reference frame, as described above.

The evolution in time of ρ(t) is described by the Liouville–von Neumann equation

tρ(t)=iH(Ω(t)),ρ(t).
(1)

Here, t indicates the time derivative, and [H, ρ] is the commutator ρH. The spin dynamics is coupled to the rotational dynamics via the orientational trajectory of the spin centers, Ω(t). H is the spin Hamiltonian operator (in angular-frequency units) summing all EPR-relevant interactions within the spin center and between the spin center and the external static magnetic field. H is implicitly time-dependent through the time dependence of Ω. During irradiation with the microwave field, H is explicitly time-dependent. However, for our purpose of simulating CW EPR spectra, we do not need to consider the interaction of the spin system with the microwave irradiation. It is possible to rewrite Eq. (1) as

tρ(t)=iH×(Ω(t))ρ(t),
(2)

where H× is called the Hamiltonian commutation superoperator and is defined by its operation on ρ, H×ρ ≡ [H, ρ]. The space spanned by all possible ρ is called the Liouville space. H× is an operator on this space.

Instead of dealing with explicit orientational trajectories Ω(t) and the associated dynamic equation, the rotational rigid-body dynamics is modeled using the orientational distribution of the spin centers in the sample at time t, P(Ω, t). The time evolution of this distribution is described by the Fokker–Planck-type differential equation6 

tP(Ω,t)=Γ(Ω)P(Ω,t),
(3)

where Γ(Ω) is the diffusion operator representing the rotational dynamics. Γ is assumed to be independent of the spin degrees of freedom, reflecting the reasonable and accurate assumption that the rotational dynamics of a spin center is independent of its spin state.

The explicit form of Γ depends on the model that is used to describe the rotational dynamics.13 Here, we focus on hindered Brownian rotational diffusion (HBD) in the presence of an orienting potential U(Ω). Brownian motion assumes the absence of inertial motion, i.e., it assumes that the diffusion process is Markovian (memoryless). U(Ω) describes an orientation-dependent effective potential energy for the spin center which is a result of the interaction of the spin center with its immediate nanoenvironment. This potential encodes that different orientations of the spin center have different energies. It leads to a systematic torque on the spin center (in the downhill direction on the potential-energy surface). The associated Γ is6 

Γ=i,j=x,y,zRijJiJj+1kBTi,j=x,y,zRijJi(JjU(Ω)).
(4)

Here, Jx, Jy, and Jz are differential angular-momentum operators around the axes x, y, and z of a body-fixed reference frame. Rij are the real-valued elements of the body-fixed diffusion tensor, which is symmetric (Rij = Rji) and is assumed to be time-independent.

During a CW EPR experiment, the orientational distribution is at thermal equilibrium at all times. The thermal-equilibrium orientational distribution Peq is stationary (tPeq = 0) and is given by

Peq(Ω)=Z1eU(Ω)/kBT
(5)

with the partition function

Z=dΩeU(Ω)/kBT,
(6)

where kB is the Boltzmann constant and T is the temperature.

The dynamical equations for the spin and the orientational degrees of freedom [Eqs. (2) and (3)] can be combined into a single equation, the stochastic Liouville equation (SLE)8,10,16

tρ(Ω,t)=iH×(Ω)ρ(Ω,t)Γ(Ω)ρ(Ω,t).
(7)

Here, ρ(Ω, t) is the total spin density operator for all spin centers with orientation Ω at time t, no matter which orientation Ω0 they had initially,

ρ(Ω,t)=dΩρ(t|Ω)P(Ω,t|Ω,0)P(Ω,0).
(8)

P(Ω′, 0) is the initial orientational distribution at time zero, P(Ω, t|Ω′, 0) is the distribution at time t given that the orientation was Ω′ initially, and ρ(t|Ω′) is the spin density matrix at time t given that the orientation was Ω′ initially. Note that ρ(Ω, t) is different from ρ(t) in Eq. (1).

In the presence of a potential, Γ of Eq. (4) is not Hermitian. It is advantageous to transform it to Hermitian form.6 This can be achieved by the transformation Γ̃=Peq1/2ΓPeq1/2. The diffusion equation then becomes tP̃=Γ̃P̃, with the scaled distribution function P̃=Peq1/2P. Defining in addition the scaled density ρ̃=Peq1/2ρ, the SLE with Γ̃ reads

tρ̃(Ω,t)=iH×(Ω)ρ̃(Ω,t)Γ̃(Ω)ρ̃(Ω,t),
(9)

where now both H× and Γ̃ are Hermitian. Since H× and Γ̃ are time-independent, the integral of this equation is

ρ̃(t)=e(iH×+Γ̃)tρ̃(0).
(10)

(We omit from now on the explicit indication of the dependence on Ω.) An explicit form of the Hermitianized diffusion superoperator is6 

Γ̃=i,j=x,y,zJi(JiU)/2kBTRijJj+(JjU)/2kBT.  
(11)

In the body-fixed diffusion frame, the diffusion tensor is diagonal, i.e., Rij = Riδij. We use this as the body-fixed reference frame without loss of generality but with significant simplifications in the expressions for Γ̃. In the diffusion frame, Γ̃ from Eq. (11) simplifies to37 

Γ̃=i=x,y,zRiJi2+i=x,y,zRi12kBTJi2U14kB2T2(JiU)2,
(12)

where we have separated the potential-independent and potential-dependent terms.

The potential function U is best expanded in a complete orthogonal basis of Wigner functions DM,KL because DM,KL have simple transformation properties under rotation—they are simultaneous eigenfunctions of J2 and Jz. The expansion is

U(ΩUD)=kBTL,M,KλM,KLDM,KL(ΩUD),
(13)

where the − sign is a matter of convention. Here, ΩU→D is the triplet of Euler angles (α, β, γ) that describes the orientation of the body-fixed diffusion frame (D) relative to the space-fixed potential frame (U) (see Fig. 2). The Wigner functions DM,KL are defined as38,39

DM,KL(ΩUD)=DM,KL(α,β,γ)=eiMαdM,KL(β)eiKγ,
(14)

where dM,KL are real-valued functions consisting of sums of products of cos(β/2) and sin(β/2). The sum in Eq. (13) runs over all possible combinations of integer rank L and projections M and K (L ≥ 0; −LML; −LKL). λM,KL are dimensionless coefficients that may be complex-valued. Since U is real-valued and dM,KL=(1)MKdM,KL, the coefficients satisfy λM,KL=(1)MKλM,KL*.

Applying the Ji operators to the Wigner expansion of U, the potential-dependent part of Γ̃ in Eq. (12) reduces to a linear combination of Wigner functions. Γ̃ becomes

Γ̃=i=x,y,zRiJi2+L,M,KX̃M,KLDM,KL(ΩUD).
(15)

The scalar expansion coefficients X̃M,KL depend on Ri and on λM,KL and are given in  Appendix C. They have the symmetry X̃M,KL=(1)MK(X̃M,KL)*. If all λM,KL are real-valued, then all X̃M,KL are real-valued as well. The largest L for which X̃M,KL is non-zero is twice the largest L for which λM,KL is non-zero. Note that the potential-dependent part of Γ̃ is a purely multiplicative operator, whereas the potential-independent part contains the differential operators Ji.

The use of the complete expansion [Eq. (13)] and the ensuing expressions for X̃M,KL constitute extensions of the existing theory, where the expansion is limited to terms with M = 0, with low even values of L and K (L = 2, 4; K = 0, 2, 4) and with real-valued coefficients λM,KL. These specializations stem from the fact that the theory was initially developed for spin probes in uniaxial liquid crystals and model membranes with nanoenvironments of cylindrical symmetry (Dh) and required only high-symmetry potentials.40 However, the nanoenvironment hindering the rotational motion of spin labels on proteins is generally of much lower symmetry (C1), and the orientational potential expansions, consequently, needs more terms. Therefore, the general expansion of Eq. (13) is necessary for being able to broadly apply the HBD model to protein-attached spin labels.

A general EPR spin Hamiltonian for a system of coupled electrons and magnetic nuclei is

H=iμBBg(i)S(i)kμNgn(k)BI(k)+iS(i)D(i)S(i)+i,kS(i)A(i,k)I(k)+i<jS(i)J(i,j)S(j)+kI(k)P(k)I(k).
(16)

The first two terms are the electron and nuclear Zeeman interactions, with the Bohr magneton μB, the external magnetic field vector B, the ith electron g-tensor g(i), the ith electron spin vector S(i), the nuclear magneton μN, the kth nuclear g factor gn(k), and the kth nuclear spin vector I(k). The third term contains all zero-field interactions for electron spins S(i) ≥ 1, with the zero-field interaction tensors D(i). The fourth term collects all hyperfine interactions between electrons and nuclei, with the hyperfine tensors A(i,k). The fifth term sums over all electron–electron interactions, with the associated coupling tensors J(i,j). The last term represents the electric quadrupole interactions for nuclei with I(k) ≥ 1, with the quadrupole tensors P(k). We assume that all vectors and tensors are represented in the same space-fixed frame. For a tumbling molecule, B is time-independent in this frame, and all the interaction tensors (g, D, A, J, P) are time-dependent.

The spin Hamiltonian of Eq. (16) is a sum over bilinear forms,

H=μaμFμbμ=μi,j=x,y,z(Fμ)ij(aμ)i(bμ)j,
(17)

where aμ and bμ are the spin vector operators or the magnetic field vector, Fμ are the time-dependent Cartesian tensors, and the axes x, y, and z refer to a space-fixed frame. The index μ runs over all interaction terms. Table I lists the correspondence between Eqs. (17) and (16).

TABLE I.

Spin Hamiltonian interaction terms, aμ · Fμ · bμ.

Interaction typeaμFμbμ
Electron Zeeman B +μBg(i) S(i) 
Nuclear Zeeman B μNgn(k) I(k) 
Zero-field splitting S(i) D(i) S(i) 
Hyperfine S(i) A(i,k) I(k) 
Electron–electron coupling S(i) J(i,j) S(j) 
Nuclear quadrupole I(k) P(k) I(k) 
Interaction typeaμFμbμ
Electron Zeeman B +μBg(i) S(i) 
Nuclear Zeeman B μNgn(k) I(k) 
Zero-field splitting S(i) D(i) S(i) 
Hyperfine S(i) A(i,k) I(k) 
Electron–electron coupling S(i) J(i,j) S(j) 
Nuclear quadrupole I(k) P(k) I(k) 

For modeling rotational dynamics, it is advantageous to reformulate H from a sum over terms with distinct physical origin to a sum over terms with distinct rotational properties. For this, we rewrite H as a sum over scalar products of irreducible spherical tensors and tensor operators,13,18,41

H=μl=02Fμ(l)T(l)(aμ,bμ)=μl=02m=ll(1)mFμ(l,m)T(l,m)(aμ,bμ).
(18)

Here, Fμ(l) are spherical tensors constructed from the matrix elements of Fμ, and T(l)(aμ, bμ) are spherical tensor operators constructed from the Cartesian components of aμ and bμ. l is the rank of the spherical tensor. Each Cartesian tensor is decomposed into three spherical tensors with ranks l = 0, 1, and 2. These three tensors have distinct rotational properties. Each spherical tensor has 2l + 1 components (m = −l, …, l), indicated by F(l,m) and T(l,m). The scalars Fμ(l,m) and the operators T(l,m)(aμ, bμ) can be constructed in a straightforward fashion from Fμ, aμ, and bμ. All the required expressions are listed in  Appendix A.

As in Eqs. (16) and (17), also in Eq. (18), all vectors and interaction matrices are represented in a space-fixed frame. Choosing the laboratory frame (L), indicating it by an additional subscript, and using (−1)mF(l,−m) = (−1)lF(l,m)* (see  Appendix A), we write the Hamiltonian as

H=μlm(1)lFμ,L(l,m)*T(l,m)(aμ,L,bμ,L).
(19)

The laboratory frame components Fμ,L(l,m) of Fμ(l) are expressed in terms of the time-independent diffusion frame components Fμ,D(l,m) using

Fμ,L(l,m)*=mFμ,D(l,m)Dm,mlΩDL*=mFμ,D(l,m)*Dm,mlΩLD,
(20)

where Dm,ml are again Wigner functions as defined in Eq. (14) and ΩL→D represents the Euler angles that parameterize the transformation from the laboratory frame L to the body-fixed diffusion frame D (see Fig. 2). We decompose this transformation into two consecutive transformations via the intermediate potential frame (U),

Dm,mlΩUD=mDm,mlΩLUDm,mlΩUD.
(21)

Combining this with the previous two equations and defining the operators

Pm,ml=1μ(1)lFμ,D(l,m)*T(l,m)(aμ,L,bμ,L)
(22)

and

Qm,ml=mDm,mlΩUDPm,ml
(23)

gives the compact expression

H=lm,mDm,mlΩUDQm,ml.
(24)

Since both L and U are space-fixed stationary frames, the only time dependence is in the stochastically varying orientation ΩU→D. The rotational time dependence is now fully isolated in the Wigner function prefactors. The time-independent operators Qm,ml, which we call rotational basis operators (RBOs), contain all the internal specifics of the spin system as well as the orientation of the potential frame relative to the lab frame. Note that this expression is similar to the potential-dependent term of the diffusion operator in Eq. (15).

Equation (24) is general in the sense that the RBOs can be constructed using the same procedure regardless of the number of spins in the spin system, the number of interaction terms, or the nature or relative size of those terms. No matter how large or complex the tumbling spin system is, all of the information about the system that is needed to calculate the EPR spectral response is collected into 35 RBOs (one rank-0, 9 rank-1, and 25 rank-2). If all the interaction matrices are symmetric, as is commonly the case, all rank-1 RBOs vanish. If the symmetry is high (for example, axial and collinear interaction tensors), then the number of non-zero RBOs reduces further. Equation (24) covers two separate situations: (i) a single potential-frame orientation ΩUD (e.g., an oriented membrane or a protein crystal) and (ii) an orientational distribution of potential-frame orientations (such as a solution of essentially immobile proteins or liposomes). The latter model with a disordered static distribution of proteins, but mobile spin labels, has been termed the MOMD (microscopic order macroscopic disorder) model.14 For calculation in such cases, the Pm,ml operators are precomputed, and the Qm,ml operators are efficiently computed from Pm,ml for each potential-frame (protein) orientation, ΩUD, without the need of recomputing any matrix elements.

The expression for the Hamiltonian presented above is quite general: it includes situations with tensor eigenframes tilted arbitrarily with respect to the diffusion frame; it includes rank-1 terms; and it permits any number of spins and types/strengths of interactions. In addition, the pre-calculation of RBOs renders simulations efficient, particularly for disordered samples.

The SLE, together with the expressions for the diffusion operator and the spin Hamiltonian, can be used to calculate the signal from any type of EPR experiment (pulse or CW). Here, we focus on CW EPR. The frequency-swept CW EPR spectrum is proportional to the Fourier–Laplace transform of the free-induction decay (FID) following a non-selective microwave pulse. (Alternatively, it can be calculated using linear-response theory.16,42) We start with the equilibrium density distribution ρeqSzPeq (corresponding to ρ̃eqSzPeq1/2), where z is the laboratory z axis. The detectable part of the state immediately after a microwave pulse with oscillating magnetic-field component along the laboratory y axis is SxPeq, where x is the laboratory x axis. The FID signal is the expectation value of Sx (or of iSx(i) if there are multiple electron spins) for this state,

Sx(t)=Sx|ρ(Ω,t)=Sx|Peq1/2ρ̃(Ω,t)=SxPeq1/2|e(iH×+Γ̃)tSxPeq1/2,
(25)

where H× is the static spin Hamiltonian superoperator in the absence of the microwave field. The notation ⟨u|v⟩ indicates integration of uv over all orientations and trace over all spin space.

The real part of the Fourier–Laplace transform of the FID gives the spectrum

I(ω,B)0Sxcos(ωt)dt=Re0Sxeiωtdt=ReSxPeq1/2|iH×(B)+Γ̃iω1SxPeq1/2.
(26)

Here, ω is the microwave angular frequency and B is the magnitude of the applied static magnetic field.

In order to evaluate the expression in Eq. (26), the quantities involved (H×, Γ̃, and SxPeq1/2) are expanded in an appropriate basis that encompasses both the quantum spin states and the classical orientational degrees of freedom. This transforms the integral expression in Eq. (26) into a linear-algebra expression (vector × matrix inverse × vector) that can be solved efficiently using numerical methods.

We use a basis where each basis function σξ is a direct product of a spin basis function σ and an orientational basis function ξ,

σξ=σξ.
(27)

The spin basis functions, σ, are given by the complete set of single-transition operators in the Zeeman (high-field) basis.43 Each single-transition operator is parameterized by a pair of projection numbers for each spin in the system

σ=m1,m1,m2,m2,=(m1m1)(m2m2),
(28)

where the spin projection quantum numbers (along the z axis of the potential frame) mi and mi for spin i run from −Si to +Si or −Ii to Ii for nuclear spins. The basis functions in Eq. (28) are orthonormal, and there are a total of i(2Si+1)2k(2Ik+1)2 of them. Instead of using mi and mi, the single-transition operators can also be indexed by pi = mimi and qi = mi + mi. This facilitates the application of spin-space truncation schemes, such as limiting the basis to single-quantum transitions (|pi| ≤ 1) and application of the high-field approximation (only pi = +1).13 

The orientational basis functions are normalized Wigner functions [see Eq. (14)] of ΩU→D,

ξ=LMK=2L+18π2DM,KL(ΩUD)
(29)

with the full range of L, M, and K. These functions are orthonormal, i.e.,

ξ1|ξ2=L1M1K1|L2M2K2=δL1,L2δM1,M2δK1,K2.
(30)

They are a convenient choice because they have simple properties under rotation and are also the coefficients for the transformation of the spherical tensors and tensor operators between coordinate frames. The complete set of Wigner functions is infinite, and in practice, the basis must be truncated. A simple truncation scheme is to impose separate upper limits on even L, odd L, |K|, and |M|. A more sophisticated pruning procedure has been developed.15,17 In general, slower rotational diffusion and larger interaction tensor anisotropies require larger orientational basis sets to produce accurate and converged simulated spectra. A crucial practical point is the confirmation that a simulated spectrum is sufficiently converged as a function of orientational basis size.

In the basis of Eq. (27), the matrix elements of the Hamiltonian superoperator in Eq. (26) are

σ1ξ1H×σ2ξ2=lσ1|Q×lΔM,ΔK|σ2ξ1|DΔM,ΔKl|ξ2
(31)

with ΔM = M1M2 and ΔK = K1K2. The sum runs over all l = 0, 1, 2 that also satisfy |L1L2| ≤ lL1 + L2 and l ≥ |ΔM| and l ≥ |ΔK|. Since l ≤ 2, the matrix is at most pentadiagonal in each of L, M, and K. Any zero RBO Q×lΔM,ΔK will further thin this nonzero bandedness along M and K.

The expressions for the Wigner function matrix elements needed in Eq. (31) are given in  Appendix B. The matrix representations of Q×lΔM,ΔK are constructed as follows: (1) The spin matrices Sx(i), Sy(i), Sz(i), Ix(k), Iy(k), and Iz(k) (where x, y, and z refer to the laboratory frame) in the standard Hilbert-space mi representation are constructed for each spin, and spherical tensor matrices are constructed using the expressions from Table III; (2) the Qm,ml matrices are constructed using Eqs. (22) and (23); (3) the corresponding Liouville-space matrices are constructed using the Kronecker product ⊗ according to

Q×=IQQTI
(32)

in which I is the identity matrix of the same size as the Q matrix; and (4) any spin-space truncations are applied.

The matrix elements of Γ̃ in the chosen basis are

σ1ξ1|Γ̃|σ2ξ2=δσ1,σ2ξ1|Γ0|ξ2+δσ1,σ2LX̃ΔM,ΔKL×ξ1|DΔM,ΔKL|ξ2
(33)

with the required matrix elements of Γ0 given in  Appendix C. Γ̃ is diagonal in the spin quantum numbers. The sum runs over all L that satisfy |L1L2| ≤ LL1 + L2 and L ≥ |ΔM| and L ≥ |ΔK|.

In general, the Γ̃ matrix is Hermitian. In our choice of frame where R is diagonal (the diffusion frame) and in the absence of a potential, Γ is diagonal if the diffusion tensor is isotropic or axial; otherwise, it is tridiagonal. For a non-zero potential, if the potential contains only terms with M = K = 0, or if all λ coefficients are real-valued, then Γ̃ is real-valued and symmetric; otherwise, it is complex-valued and Hermitian.

Finally, the vector elements of SxPeq1/2 in the chosen basis are

σξ|SxPeq1/2=σ|Sxξ|Peq1/2.
(34)

σ|Sx is obtained by vectorizing the mi matrix representation of Sx in column-major order, i.e., by stacking the columns of the matrix. The vector elements of Peq1/2 are given by

ξ|Peq1/2=Z1/22L+18π2  dΩDM,KL*eU/2kBT.
(35)

Analytical evaluation of the triple integral in this expression is only possible for very simple potentials.17 In general, the integral needs to be evaluated numerically. If U does not depend on α and/or γ (representing the cylindrical symmetry of the environment and/or the spin label), then the dimensionality of the integral reduces.

The norm-squared of the Peq1/2 vector equals 1 in the infinite-basis limit.38 If the basis is excessively truncated, this will be significantly less than 1. This can serve as a diagnostic of basis overtruncation.

With matrix representations of all relevant quantities in hand, we now describe how to evaluate Eq. (26). Denoting with A the matrix representation of iH×+Γ̃ in the chosen truncated basis, with b the vector representation of SxPeq1/2 in the same basis, and with I the identity matrix of the same size as A, the spectral line shape function of Eq. (26) is

I(ω)=Reb(AiωI)1b.
(36)

Evaluation of this as a function of ω for a fixed external magnetic field (i.e., fixed A) gives the frequency-swept EPR spectrum.

A variety of approaches are available to numerically evaluate the right-hand side of Eq. (36). The straightforward diagonalization method16,44 computes eigenvectors and eigenvalues of A (such that A = UΛU−1 with the diagonal matrix of eigenvalues Λ and the matrix of eigenvectors U) and transforms the expression into the eigenbasis of A, resulting in a sum over Lorentzian lines

I(ω)=RekSkΛkkiω  Sk=(bU)k(U1b)k,
(37)

which can be easily evaluated as a function of ω. In the above equation, both Sk and Λkk are in general complex-valued. The real and imaginary parts of Λkk give the linewidth and the line position, respectively. Sk determines the amplitude and phase of the line. Another method uses a linear-equation solver to calculate the solution x(ω) of (A − iωI)x(ω) = b for each ω and then obtains the spectral values via I(ω) = Re(bx(ω)). This method does not explicitly break the signal down into its line shape components. Both methods are general and work for any A without symmetry restrictions. A variety of well-established numerical algorithms (conjugate gradient, Lanczos, etc.) can be used for diagonalization and as linear solvers.

A much faster and very elegant method is applicable if A is complex symmetric. In this case, the complex symmetric Lanczos algorithm45 is used to tridiagonalize A with b as the starting vector.11,12,16,17 The spectral function is then represented as a continued-fraction expansion containing the matrix elements of the tridiagonal matrix, and it can be evaluated very efficiently left-to-right using the modified Lentz method.46 This method is fast because it converges rapidly as a function of iteration count and because the spectral function can be cheaply evaluated for the entire desired range of frequencies with only one tridiagonalization. The continued-fraction expansion is possible only if A is complex symmetric, i.e., if the matrix representations of both H× and Γ̃ are real-valued. In the σξ basis, this is only the case for a very small class of high-symmetry situations. However, under certain conditions, it is possible to convert A to complex symmetric form by transforming it to a new K-symmetrized basis with functions σLMK¯jK that are parameterized by L, M, K¯, and jK, where K¯=|K|, and jK = ±1 for K ≠ 0 and jK = (−1)L for K = 0.13,17 The transformation is achieved using A=TKATK, with the elements of the unitary transformation matrix TK given by

σLMK¯jK|σLMK=δσ,σjK2(1+δK,0)δL,LδM,M×δK¯,K+jK(1)L+KδK¯,K.
(38)

In a basis ordering where the functions with identical L, M, and K¯ are adjacent, the matrix TK is block diagonal, consisting of a sequence of 2 × 2 blocks (for K¯0) and 1 × 1 blocks (for K¯=0) along the diagonal.

Prior work derived explicit expressions for the matrix elements of H× and Γ̃ in this new basis for some special cases.13,17,37 Here, we perform this transformation numerically at the matrix level. Since all matrices are sparse (the number of non-zero elements in TK only grows linearly with basis size), this operation is very efficient.

In the K-symmetrized basis, the matrix of any Hamiltonian of the form given in Eq. (17) is real-valued, irrespective of the number, symmetry, and relative orientation of interaction tensors—as long as the tilt between the potential frame and the lab frame only involves a β angle (α = γ = 0). The situation is more complicated for Γ̃. In the LMK basis, Γ̃ is real-valued for arbitrary complicated potentials. However, in the K-symmetrized basis, it is real-valued only if all potential coefficients are real-valued and the potential contains only terms with either all M = 0 or all K = 0 (or both). Otherwise, the transformation will generate a complex-valued Hermitian Γ̃. In these cases, the transformation is not useful, and the spectral line shape function must be calculated using the diagonalization or linear-solver methods instead of the Lanczos method.

Numerical aspects of solving Eq. (36) efficiently and robustly are complicated, in light of the fact that iterative methods such as Lanczos can be numerically unstable. Previous work has carefully evaluated the merits of Lanczos vs conjugate-gradient methods.16,17

The line shape of Eq. (36) is expressed as a function of microwave frequency ω, in the presence of a constant applied magnetic field. In practice, CW EPR experiments are performed by varying the external field strength B over a range from Bmin to Bmax while irradiating the sample at fixed frequency ωmw.

For the field-swept spectrum, a separate Hamiltonian at each field point B is required, as evident from Eq. (26). Since the Hamiltonian is linear in the magnetic field, Hamiltonians for different field values can be constructed very efficiently. For this, each RBO Qm,ml is separated into a field-dependent term and a field-independent term,

Qm,ml(B)=Qm,m0,l+BQm,m1,l.
(39)

The field-independent RBOs Qm,m0,l and Qm,m1,l are pre-calculated, and then Eq. (39) is used to assemble Qm,ml for each field point, thus providing minimal overhead in Hamiltonian re-calculation. To minimize the number of field points for which the spectral function needs to be evaluated, an adaptive iterative bisection interpolation method analogous to one used for constructing Zeeman energy diagrams can be used.47 

If the complex symmetric Lanczos method is applicable (i.e., if the matrix representation of iH×+Γ̃ is complex symmetric), then the calculation of a frequency-swept spectrum is orders of magnitude more efficient than the calculation of a field-swept spectrum. In this case, the field-swept spectrum can be approximated by a frequency-swept spectrum,

I(ωmw,B)I(ω,B0).
(40)

Here, B0 = (Bmin + Bmax)/2 is the center of the desired field range, and the auxiliary frequency ω is obtained from B in one of two ways,

(i)  ω=ωmwB0B,(ii)  ω=ωmwμBgisoBB0.
(41)

giso is the isotropic g-value of the spin system. Approach (i) is exact for systems with g anisotropy and no other anisotropy (hyperfine, zero-field). For other systems, it leads to systematic errors in the resonance field positions. Approach (ii) is accurate if g is isotropic and the Zeeman interaction dominates; otherwise, there are systematic errors in the resonance field positions. At X-band, the errors in approach (ii) are very small for nitroxides, but they can be significant for copper. Approach (i) and approach (ii) give identical frequencies for B = B0 and for B = ℏωmw/μBgiso.

Finally, the field-modulated first-harmonic spectrum is obtained from the field-swept absorption spectrum via pseudo-field modulation48 or, in the limit of zero modulation amplitude, via numerical differentiation.

The theory outlined above is implemented in version 6 of the open-source MATLAB-based EPR simulation package EasySpin5,34,35 in the function chili.

All matrices are sparse and are handled in MATLAB’s compressed sparse column (CSC) format. The spectral line shape can be calculated using either the diagonalization, the linear-solver, the Lanczos method (only if applicable, given the symmetry of Γ̃), or the biconjugate gradient stabilized method as implemented in the MATLAB function bicgstab.

The integrals necessary for the vector elements in Eq. (34) are evaluated numerically using the 15-point Gauss–Kronrod quadrature as implemented in the MATLAB function integral, with absolute and relative error tolerances of 10−6. Even after using selection rules to reduce the number and dimensionality of necessary numerical integrations, the evaluation of the root-equilibrium-distribution vector [Eq. (34)] constitutes the major performance bottleneck in the overall calculation in the presence of low-symmetry potentials.

Wigner 3-j symbols for Eqs. (B4) and (C2) are evaluated using specialized expressions for small angular momenta [min(L1, L, L2) ≤ 2]38,39,49 and using a general expression otherwise.50 For large angular momenta [max(L1, L, L2) > 20], arbitrary-precision integer arithmetic (using the Java class BigInteger) is used for intermediate results. This results in a significant computational bottleneck for large angular momenta. Values of 3-j symbols calculated via these methods in EasySpin have been extensively tested for angular-momentum values up to 5000 against arbitrary-precision results obtained using Wolfram Mathematica 11.

In this section, we illustrate the scope of the extended theory by demonstrating examples of simulated CW EPR spectra for a high-spin system, for lower-symmetry orienting potentials, and for a spin system with multiple nuclei. We also discuss some numerical aspects.

As an example of a high-spin system, Fig. 3 shows a series of simulated spectra of a prototypical Gd(III) complex over a range of rotational time scales from the fast-motion regime to the quasi-rigid limit. Gd(III) is a high-spin 4f7 ion with an S = 7/2 ground state. Its magnetic properties are described by the spin Hamiltonian ℏH = μBgBS + SDS with isotropic g-value and the diagonal zero-field tensor D = diag(−D + E, −DE, 2D).

FIG. 3.

Simulated field-swept EPR spectra over a range of rotational time scales for a Gd(III) complex with isotropic g-value 2 and axial zero-field splitting D/h = 500 MHz and E/h = 0, at a microwave frequency of 9.2 GHz, with a residual Gaussian broadening of 2 mT FWHM. Black: motional spectra without orientational potential, with isotropic rotational correlation times τc = 1/6R indicated next to each spectrum. Gray: rigid-limit spectrum, calculated via spin Hamiltonian matrix diagonalization. All spectra are normalized to equal area. The orientational basis was truncated at (Lmaxeven,Lmaxodd,Mmax,Kmax)=(50,0,4,0), and the spin basis was limited to transitions with |mSmS| ≤ 1. The basis size is 9546 (222 orientational and 43 spin basis functions).

FIG. 3.

Simulated field-swept EPR spectra over a range of rotational time scales for a Gd(III) complex with isotropic g-value 2 and axial zero-field splitting D/h = 500 MHz and E/h = 0, at a microwave frequency of 9.2 GHz, with a residual Gaussian broadening of 2 mT FWHM. Black: motional spectra without orientational potential, with isotropic rotational correlation times τc = 1/6R indicated next to each spectrum. Gray: rigid-limit spectrum, calculated via spin Hamiltonian matrix diagonalization. All spectra are normalized to equal area. The orientational basis was truncated at (Lmaxeven,Lmaxodd,Mmax,Kmax)=(50,0,4,0), and the spin basis was limited to transitions with |mSmS| ≤ 1. The basis size is 9546 (222 orientational and 43 spin basis functions).

Close modal

The simulations in Fig. 3 are plotted as absorption spectra. With a rotational correlation time τc = 0.1 ns (R = 1/6τc = 1.7 rad2 ns−1), the spectrum consists of a single peak because rotational motion is so rapid that the zero-field interaction is almost averaged out. At τc = 100 ns (R = 1.7 rad2μs−1), the motional EPR spectrum is indistinguishable from the rigid-limit simulation.

Figure 4 illustrates the effect of different potential energy functions on the X-band spectrum of a nitroxide radical in the quasi-rigid limit with τc = 10 ns. Each simulation has a different orientational potential defined by the pair of coefficients λM,KL=(1)MKλM,KL = 1, with various combinations of L, M, and K. The gray spectrum is calculated in the absence of an orientational potential. The spectral line shape depends strongly on the particular term. The size of the deviation from the unhindered potential-free spectrum is dependent on (L, M, K) as well. For a given λ, this deviation is generally smaller for odd values of M and K (and L) than for even values.

FIG. 4.

The effect of various orientational potentials on the EPR spectrum of a nitroxide radical in the slow-motion regime (τc = 10 ns) at 9.5 GHz, with g-tensor principal values (2.009, 2.006, 2.002) and 14N hyperfine tensor principal values (20, 20, 100) MHz. Each simulation has a pair of nonzero potential coefficients λM,KL=(1)MKλM,KL = 1 with (L, M, K) given in the figure. The potential frame is collinear with the lab frame. The gray spectrum is calculated in the absence of an orientational potential. The orientational basis was truncated at (Lmaxeven,Lmaxodd,Mmax,Kmax)=(10,3,6,6), and the spin basis is complete. The basis size is 13 248 (368 orientational and 36 spin basis functions).

FIG. 4.

The effect of various orientational potentials on the EPR spectrum of a nitroxide radical in the slow-motion regime (τc = 10 ns) at 9.5 GHz, with g-tensor principal values (2.009, 2.006, 2.002) and 14N hyperfine tensor principal values (20, 20, 100) MHz. Each simulation has a pair of nonzero potential coefficients λM,KL=(1)MKλM,KL = 1 with (L, M, K) given in the figure. The potential frame is collinear with the lab frame. The gray spectrum is calculated in the absence of an orientational potential. The orientational basis was truncated at (Lmaxeven,Lmaxodd,Mmax,Kmax)=(10,3,6,6), and the spin basis is complete. The basis size is 13 248 (368 orientational and 36 spin basis functions).

Close modal

As an example of slow-motion simulations for spin systems with multiple magnetic nuclei, Fig. 5 shows simulations of a model copper(II) complex in which two of the ligand atoms are nitrogen-14 nuclei and the other ligand atoms are non-magnetic. The copper-63 nucleus has nuclear spin I = 3/2, and each of the nitrogen-14 nuclei has nuclear spin I = 1. An orientational potential with the coefficients λ2,22=λ2,22=1 was used.

FIG. 5.

Simulated 9.5 GHz spectra of a copper(II) complex with two nitrogen-14 ligands (exact field sweep). Parameters: g = [2.0, 2.0, 2.25], A63Cu = [50, 50, 500] MHz, A14N = [80, 80, 10] MHz. The orientational potential is λ2,22=λ2,22=1, with the potential frame aligned with the laboratory frame. The orientational basis was truncated at (Lmaxeven,Lmaxodd,Mmax,Kmax)=(10,0,4,4) and the spin basis at a maximum nuclear coherence order of 3 for 63Cu and 1 for 14N. The basis size is 153 664 (196 orientational and 784 spin basis functions).

FIG. 5.

Simulated 9.5 GHz spectra of a copper(II) complex with two nitrogen-14 ligands (exact field sweep). Parameters: g = [2.0, 2.0, 2.25], A63Cu = [50, 50, 500] MHz, A14N = [80, 80, 10] MHz. The orientational potential is λ2,22=λ2,22=1, with the potential frame aligned with the laboratory frame. The orientational basis was truncated at (Lmaxeven,Lmaxodd,Mmax,Kmax)=(10,0,4,4) and the spin basis at a maximum nuclear coherence order of 3 for 63Cu and 1 for 14N. The basis size is 153 664 (196 orientational and 784 spin basis functions).

Close modal

The simulations show how the spectrum changes with rotational correlation times from the fast-motion regime (small τc) to the quasi-rigid limit (large τc). At τc = 0.1 ns, the spectrum is close to the fast-motion regime. Most of the structure in the spectrum is due to averaged nitrogen hyperfine interactions, with some weak residual features in the low-field region from copper hyperfine splittings. As τc increases, the copper peaks become more distinct and reveal the non-isotropic orientational potential.

A major challenge with the simulation of multi-spin spectra using the SLE is the rapid factorial growth of the spin basis size with the number of spins. Three approaches at truncating this space can be applied individually or simultaneously: (a) limit the coherence orders |pk| for each nucleus k; (b) limit the total maximum nuclear coherence order |kpk|; and (c) utilize the high-field approximation where for the electron spin S only the subspace with pS = +1 is utilized. Truncation approaches based on subpartitioning based on the spin coupling topology, well established in NMR,26 are only marginally beneficial for EPR systems since the latter typically have a star topology (all nuclei coupled to a central electron spin) that cannot be partitioned into subgraphs.

A perturbation-based approach at reducing the spin-space dimension is applicable if the hyperfine couplings are such that a subset of nuclei have hyperfine couplings A small enough to fall into the fast-motion regime (A ≪ 1/τc). In this approach, the SLE-simulated spectrum of the spin system containing only the nuclei that fall into the slow-motion regime is convolved with the perturbation-based simulated spectrum of the remaining fast-motion-regime nuclei. This post-convolution approach factors the spin space and has been implemented for isotropic diffusion in the absence of a potential.51,52 It is not applicable if the hyperfine tensors in the system are of comparable magnitude.

The application of these truncations and approximations is crucial for making many-spin simulations feasible. Numerical experiments must be performed to carefully ascertain that any truncation does not degrade the convergence of the simulated spectrum. Practically, simulations of systems with more than about six spins are very memory- and time-consuming. Luckily, very few slow-motion CW EPR spectra resolve splittings or other features from more than a few electron or nuclear spins.

Lanczos and conjugate-gradient methods are numerically unstable under certain conditions since they do not produce a fully orthogonal basis. In EPR simulations, it is known that this instability can manifest, for instance, in the presence of anisotropic diffusion tensors and very slow motion, leading to a failure to obtain a converged spectrum. In case numerical instabilities are encountered, a slight change in the orientational basis (by adding a few rotational basis functions) often mitigates the issue. As a safe fallback, other linear solvers or diagonalization can be used as described above. However, this comes at a significant computational cost. Our implementation shows reasonable stability of the Lanczos-based approach in difficult regimes. For example, Fig. 6 illustrates simulations in the presence of a strongly anisotropic diffusion tensor, with correlation times up to 1000 ns. With the chosen basis, which is not larger than necessary to obtain smooth spectra, convergence is not an issue.

FIG. 6.

Simulation of the 350 mT EPR spectrum of a doublet spin system with rhombic g-tensor with principal values (1.8, 2.0, 2.2), an ordering potential with λ002=0.5, and a strongly anisotropic axial diffusion tensor (τz = 0.1 ns; τxy = 0.3 ns, 1 ns, 3 ns, 1 ns, 30 ns, 100 ns, 300 ns, 1000 ns). The magnetic field is aligned with the axis of the potential. The quasi-rigid-limit spectrum in gray was simulated with τz = τxy = 1000 ns. The orientational basis was truncated at (Lmaxeven,Lmaxodd,Mmax,Kmax)=(60,0,4,8) for the motional spectra and at (80, 0, 2, 60) for the gray spectrum.

FIG. 6.

Simulation of the 350 mT EPR spectrum of a doublet spin system with rhombic g-tensor with principal values (1.8, 2.0, 2.2), an ordering potential with λ002=0.5, and a strongly anisotropic axial diffusion tensor (τz = 0.1 ns; τxy = 0.3 ns, 1 ns, 3 ns, 1 ns, 30 ns, 100 ns, 300 ns, 1000 ns). The magnetic field is aligned with the axis of the potential. The quasi-rigid-limit spectrum in gray was simulated with τz = τxy = 1000 ns. The orientational basis was truncated at (Lmaxeven,Lmaxodd,Mmax,Kmax)=(60,0,4,8) for the motional spectra and at (80, 0, 2, 60) for the gray spectrum.

Close modal

Another practical numerical aspect is the time cost of simulations. It strongly depends on whether the continued-fraction expansion of the spectral function can be computed (via, e.g., the complex symmetric Lanczos method). If yes, only one tridiagonalization is necessary for simulating the entire spectrum. If not, diagonalization or linear-solver methods need to be used at every field/frequency point and result in significant runtime slowdowns. For the spectra in Fig. 3, the slowdown is >100× (0.92 s vs 123.5 s). For the simulations in Fig. 4, the slowdowns range between 30× and 7000×.

In this paper, we have laid out expressions to apply the SLE methodology for HBD models to spin systems with any number or nature of interaction terms. From the familiarly defined lab-frame spin vector operators and molecular-frame Cartesian interaction tensors, ISTOs and RBOs are constructed. From the RBOs and the choice of orientational basis, the elements of the Hamiltonian superoperator and of the rotational diffusion operator are calculated. Once the calculation of the elements of H× and Γ is complete, the spectral line shape function is evaluated.

We have also presented an expression for the diffusion operator in the presence of orienting potentials of any form, even of very low symmetry. This is useful for two scenarios: (1) for fitting an effective potential to an orientational histogram obtained from a molecular-dynamics trajectory and (2) for fitting a potential expansion to an experimental CW EPR spectrum. For the first case, many terms in the potential expansion are needed. For the second case, one has to consider the limited information content and the given signal-to-noise ratio of the CW EPR spectrum and limit modeling to potentials with relatively few terms.

Care must be taken to assure the orientational basis and the spin basis are large enough to provide sufficient convergence of the spectrum. The allowable truncation level depends on the spin Hamiltonian parameters, the rotational diffusion rates, and the orientational potential.

For systems with highly anisotropic g-tensors, care must be taken to simulate explicit field sweeps when fitting experimental spectra. Using an approximate field sweep method, it is possible that the fit parameters extracted from the best-fit simulation will be inaccurate due to the error in the peak positions and spectral width of the frequency-to-field converted simulation. Explicit field sweep simulations are more time consuming, but the construction of field-independent and field-dependent RBOs allows for reasonable simulation times in small to medium-sized spin systems.

Several important numerical challenges remain. (1) One is the efficient evaluation of the integrals for the basis representation of the root equilibrium distribution Peq1/2, as given in Eq. (34). Numerical integration for this, as done in this work, constitutes a significant computational bottleneck, particularly for large orientational basis sets where the basis contains highly oscillatory large-L basis functions. (2) Another challenge is the identification of additional symmetrizing basis transformations, similar to the K-symmetrization in Eq. (38), and the M-symmetrization/truncation, which is applicable only in the high-field regime13 for M = 0 potentials. The goal of the symmetrizations is to allow for efficient basis truncation and for the application of the complex symmetric Lanczos method. (3) The factorial growth of the spin space with the number of spins represents a challenge for large spin systems since the simulation times grow factorially, despite the spin-space truncations outlined above. (4) The Lanczos algorithm can be numerically unstable and fail to converge occasionally. Future advances addressing these challenges will improve the performance and robustness of the simulations.

It is important to reiterate the assumptions underlying the HBD model: (1) The spin dynamics does not affect the rotational dynamics. This assumption is accurate in all known cases. (2) The spin center is internally rigid, i.e., the principal values and relative orientations of the interaction tensors are time-independent in a body-fixed frame. This assumption is violated in some Gd(III) complexes, where there can be significant internal conformational dynamics.31,32 In such cases, the HBD model needs to be expanded by additional degrees of freedom, with associated potential energy functions and interaction tensor dependencies. (3) The rotational motion is Brownian, i.e., there are no inertial effects. This assumption is reasonably well satisfied for spin labels in solution and for spin labels at mobile solvent-exposed sites in proteins but is less appropriate for more restricted sites, where a multi-site Markov jump model might be more successful.4,5,53 (4) The motion can be described with a single diffusion tensor (i.e., there is a single time scale). Since the rotational motion of a spin label is a consequence of simultaneous torsional motions around several different bonds, rotational motion may occur over multiple time scales. To what degree an assumption about a single time scale is adequate remains to be explored. (5) The external potential is time-independent, at least on the time scale of the rotational diffusion. This assumption disregards the fact that different parts of the nano-environment might move at different time scales and that some of these time scales might be on the order of the spin label diffusion time scale, which would complicate the situation.

Our extensions of the SLE theory now allow assumptions (4) and (5) to be tested against atomistic molecular-dynamics simulations and against experimental data. The HBD model can then be quantitatively compared to other approaches, such as Markov state models.

This work was supported by the National Science Foundation (Grant No. CHE-1452967, S.S.), the National Institutes of Health (Grant No. GM125753, S.S.), and the Research Corporation for Science Advancement (Award No. 23447, S.S.). We thank David E. Budil (Northeastern University) and Keith A. Earle (University at Albany) for many helpful discussions as well as EasySpin users for bug reports.

Tables II and III list irreducible spherical tensor components.41 The tables give expressions for the conversion of vectors and tensors from their Cartesian form to irreducible spherical tensor form. Table III can be derived from Table II using the Cartesian tensor C with Cij = biaj (not aibj).

TABLE II.

Irreducible spherical tensor components F(l,m) in terms of the Cartesian components Fij of the interaction tensor F.

(l, m)F(l,m)
(0, 0) 13(Fxx+Fyy+Fzz) 
(1, 0) i2(FxyFyx) 
(1, ±1) 12[(FzxFxz)±i(FzyFyz)] 
(2, 0) +23[Fzz12(Fxx+Fyy)] 
(2, ±1) 12[(Fxz+Fzx)±i(Fyz+Fzy)] 
(2, ±2) +12[(FxxFyy)±i(Fxy+Fyx)] 
(l, m)F(l,m)
(0, 0) 13(Fxx+Fyy+Fzz) 
(1, 0) i2(FxyFyx) 
(1, ±1) 12[(FzxFxz)±i(FzyFyz)] 
(2, 0) +23[Fzz12(Fxx+Fyy)] 
(2, ±1) 12[(Fxz+Fzx)±i(Fyz+Fzy)] 
(2, ±2) +12[(FxxFyy)±i(Fxy+Fyx)] 
TABLE III.

Irreducible spherical tensor operator components for a product of two vector operators a and b, with a± = ax ± iay.

(l, m)T(l,m) (a, b), from CartesianT(l,m) (a, b), from polar
(0, 0) 13(axbx+ayby+azbz) 13azbz+12(a+b+ab+) 
(1, 0) +i2(axbyaybx) 122(a+bab+) 
(1, ±1) +12[(azbxaxbz)±i(azbyaybz)] 12(a±bzazb±) 
(2, 0) +23azbz12(axbx+ayby) +23azbz14(a+b+ab+) 
(2, ±1) 12[(axbz+azbx)±i(aybz+azby)] 12(a±bz+azb±) 
(2, ±2) +12[(axbxayby)±i(axby+aybx)] +12a±b± 
(l, m)T(l,m) (a, b), from CartesianT(l,m) (a, b), from polar
(0, 0) 13(axbx+ayby+azbz) 13azbz+12(a+b+ab+) 
(1, 0) +i2(axbyaybx) 122(a+bab+) 
(1, ±1) +12[(azbxaxbz)±i(azbyaybz)] 12(a±bzazb±) 
(2, 0) +23azbz12(axbx+ayby) +23azbz14(a+b+ab+) 
(2, ±1) 12[(axbz+azbx)±i(aybz+azby)] 12(a±bz+azb±) 
(2, ±2) +12[(axbxayby)±i(axby+aybx)] +12a±b± 

F(l,m) and T(l,m) satisfy the symmetry relations

F(l,m)=(1)lmF(l,m)*,  T(l,m)=(1)lmT(l,m).
(A1)

All T(l,m) are real-valued in the conventional Zeeman basis.

The frame transformation of the spherical tensors is as follows:

Fμ,B(l)=(Fμ,A(l))TDl(ΩBA)=(Dl(ΩBA))TFμ,A(l)=(Dl(ΩAB))*Fμ,A(l),
(A2)

where we have used Dl(ΩAB) = Dl(ΩBA) with ΩBA = (α, β, γ) as in Fig. 2 and ΩAB = (−γ, −β, −α).

The Wigner matrix for the composition of two transformations is given in terms of the Wigner matrices of the component transformations by

Dl(ΩCA)=Dl(ΩBA)Dl(ΩCB).
(A3)

Using the shorthand ξi=Li,Mi,Ki [see Eq. (29)], the matrix elements for the angular-momentum operators are38,39

ξ1|J2|ξ2=δLMδK1,K2L(L+1),
(B1)
ξ1|Jz|ξ2=δLMδK1,K2(K2),
(B2)
ξ1|J±|ξ2=δLMδK1,K2±1L2(L2+1)K2(K2±1),
(B3)

where J± = Jx ± iJy. x, y, and z refer to a body-fixed frame and δLM=δL1,L2δM1,M2. Note the negative sign in the equation involving Jz.

The matrix elements of DM,KL are

ξ1|DM,KL|ξ2=  (1)K1M12L1+12L2+1×L1   L   L2M1   M   M2L1   L   L2K1   K   K2,
(B4)

where the expressions in parentheses are Wigner 3-j symbols.38,39 Due to the selection rules for the 3-j symbols, the matrix elements can be non-zero only if |L1L2| ≤ LL1 + L2 and M = M1M2 and K = K1K2. They are all real-valued and possess the symmetry

ξ2|DM,KL|ξ1=ξ1|DM,KL*|ξ2=(1)KMξ1|DM,KL|ξ2.
(B5)

Additionally, the following matrix elements are useful:

ξ1|J±DM,KL|ξ2=ξ1|DM,K±1L|ξ2L(L+1)K(K±1),
(B6)
ξ1|JzDM,KL|ξ2=ξ1|DM,KL|ξ2(K).
(B7)

The matrix elements of the isotropic part of the diffusion operator in the LMK basis are

ξ1|Γ0|ξ2=δL1,L2δM1,M2δK1,K2×R(L2(L2+1)K22)+RzK22+δL1,L2δM1,M2RdδK1,K2+2cL2,K2+1+cL2,K2++δK1,K22cL2,K21cL2,K2,
(C1)

with Rd = (RxRy)/4, R = (Rx + Ry)/2, and cL,K±=L(L+1)K(K±1).

The expansion coefficients X̃M,KL in Eq. (15) are given by

X̃M,KL=12RdλM,K+2LcL,K+1cL,K+2+RdλM,K2LcL,K1+cL,K2++RλM,KLL(L+1)K2+RzλM,KLK214(2L+1)(1)KML1,M1,K1L2,M2,K2λM1,K1L1λM2,K2L2×L1   L   L2M1   M   M2RdcL1,K1+cL2,K2+L1   L   L2K1+1   K   K2+1+RdcL1,K1cL2,K2L1   L   L2K11   K   K21+RcL1,K1+cL2,K2L1   L   L2K1+1   K   K21+RzK1K2L1   L   L2K1   K   K2.
(C2)
1.
H.-J.
Steinhoff
and
W. L.
Hubbell
, “
Calculation of paramagnetic resonance spectra from Brownian dynamics trajectories: Application to nitroxide side chains in proteins
,”
Biophys. J.
71
,
2201
2212
(
1996
).
2.
P.
Håkansson
,
P. O.
Westlund
,
E.
Lindahl
, and
O.
Edholm
, “
A direct simulation of EPR slow-motion spectra of spin labelled phospholipids in liquid crystalline bilayers based on a molecular dynamics simulation of the lipid dynamics
,”
Phys. Chem. Chem. Phys.
3
,
5311
5319
(
2001
).
3.
D.
Sezer
,
J. H.
Freed
, and
B.
Roux
, “
Simulating electron spin resonance spectra of nitroxide spin labels from molecular dynamics and stochastic trajectories
,”
J. Chem. Phys.
128
,
165106
(
2008
).
4.
D.
Sezer
,
J. H.
Freed
, and
B.
Roux
, “
Using Markov models to simulate electron spin resonance spectra from molecular dynamics trajectories
,”
J. Phys. Chem. B
112
,
11014
11027
(
2008
).
5.
P. D.
Martin
,
B.
Svensson
,
D. D.
Thomas
, and
S.
Stoll
, “
Trajectory-based simulation of EPR spectra: Models of rotational motion for spin labels on proteins
,”
J. Phys. Chem. B
123
,
10131
10141
(
2019
).
6.
L. D.
Favro
, “
Rotational Brownian motion
,” in
Fluctuation Phenomena in Solids
, edited by
R. E.
Burgess
(
Academic Press
,
1965
), pp.
79
101
.
7.
R.
Kubo
, “
Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems
,”
J. Phys. Soc. Jpn.
12
,
570
586
(
1957
).
8.
J. H.
Freed
,
G. V.
Bruno
, and
C.
Polnaszek
, “
Electron spin resonance line shapes and saturation in the slow motional regime
,”
J. Phys. Chem.
75
,
3385
3399
(
1971
).
9.
C. F.
Polnaszek
,
G. V.
Bruno
, and
J. H.
Freed
, “
ESR line shapes in the slow-motional regime: Anisotropic liquids
,”
J. Chem. Phys.
58
,
3185
3199
(
1973
).
10.
J. H.
Freed
, “
Theory of slow tumbling ESR spectra for nitroxides
,” in
Spin Labeling: Theory and Applications
, edited by
L. J.
Berliner
(
Academic Press
,
New York
,
1976
), pp.
53
133
.
11.
G.
Moro
and
J. H.
Freed
, “
Efficient computation of magnetic resonance spectra and related correlation functions from stochastic Liouville equation
,”
J. Phys. Chem.
84
,
2837
2840
(
1980
).
12.
G.
Moro
and
J. H.
Freed
, “
Calculation of ESR spectra and related Fokker–Planck forms by the use of the Lanczos algorithm
,”
J. Chem. Phys.
74
,
3757
3773
(
1981
).
13.
E.
Meirovitch
,
D.
Igner
,
E.
Igner
,
G.
Moro
, and
J. H.
Freed
, “
Electron-spin relaxation and ordering in smectic and supercooled nematic liquid crystals
,”
J. Chem. Phys.
77
,
3915
3938
(
1982
).
14.
E.
Meirovitch
,
A.
Nayeem
, and
J. H.
Freed
, “
Analysis of protein-lipid interactions based on model simulations of electron spin resonance spectra
,”
J. Phys. Chem.
88
,
3454
3465
(
1984
).
15.
K. V.
Vasavada
,
D. J.
Schneider
, and
J. H.
Freed
, “
Calculation of ESR spectra and related Fokker–Planck forms by the use of the Lanczos algorithm. II. Criteria for truncation of basis sets and recursive steps utilizing conjugate gradients
,”
J. Chem. Phys.
86
,
647
661
(
1987
).
16.
D. J.
Schneider
and
J. H.
Freed
, “
Spin relaxation and motional dynamics
,”
Adv. Chem. Phys.
73
,
387
527
(
1989
).
17.
D. J.
Schneider
and
J. H.
Freed
, “
Calculating slow motional magnetic resonance spectra: A user’s guide
,”
Biol. Magn. Reson.
8
,
1
76
(
1989
).
18.
A.
Polimeno
and
J. H.
Freed
, “
Slow motional ESR in complex fluids: The slowly relaxing local structure model of solvent cage effects
,”
J. Phys. Chem.
99
,
10995
11006
(
1995
).
19.
S. K.
Misra
, “
Simulation of slow-motion CW EPR spectrum using stochastic Liouville equation for an electron spin coupled to two nuclei with arbitrary spins: Matrix elements of the Liouville superoperator
,”
J. Magn. Reson.
189
,
59
77
(
2007
).
20.
A.
Collauto
,
M.
Zerbetto
,
M.
Brustolon
,
A.
Polimeno
,
A.
Caneschi
, and
D.
Gatteschi
, “
Interpretation of cw-ESR spectra of p-methyl-thio-phenyl-nitronyl nitroxide in a nematic liquid crystalline phase
,”
Phys. Chem. Chem. Phys.
14
,
3200
3207
(
2012
).
21.
A.
Polimeno
,
M.
Zerbetto
,
L.
Franco
,
M.
Maggini
, and
C.
Corvaja
, “
Stochastic modeling of CW-ESR spectroscopy of [60]fulleropyrrolidine bisadducts with nitroxide probes
,”
J. Am. Chem. Soc.
128
,
4734
4741
(
2006
).
22.
M.
Zerbetto
,
S.
Carlotto
,
A.
Polimeno
,
C.
Corvaja
,
L.
Franco
,
C.
Toniolo
,
F.
Formaggio
,
V.
Barone
, and
P.
Cimino
, “
Ab initio modeling of CW-ESR spectra of the double spin labeled peptide Fmoc-(Aib-Aib-TOAC)2-Aib-OMe in acetonitrile
,”
J. Phys. Chem. B
111
,
2668
2674
(
2007
).
23.
M.
Gerolin
,
M.
Zerbetto
,
A.
Moretto
,
F.
Formaggio
,
C.
Toniolo
,
M.
van Son
,
M. H.
Shabestari
,
M.
Huber
,
P.
Calligari
, and
A.
Polimeno
, “
Integrated computational approach to the electron paramagnetic resonance characterization of rigid 310-helical peptides with TOAC nitroxide spin labels
,”
J. Phys. Chem. B
121
,
4379
4387
(
2017
).
24.
Z.
Liang
,
J. H.
Freed
,
R. S.
Keyes
, and
A. M.
Bobst
, “
An electron spin resonance study of DNA dynamics using the slowly relaxing local structure model
,”
J. Phys. Chem. B
104
,
5372
5381
(
2000
).
25.
Z.
Liang
,
Y.
Lou
,
J. H.
Freed
,
L.
Columbus
, and
W. L.
Hubbell
, “
A multifrequency electron spin resonance study of T4 lysozyme dynamics using the slowly relaxing local structure model
,”
J. Phys. Chem. B
108
,
17649
17659
(
2004
).
26.
H. J.
Hogben
,
M.
Krzystyniak
,
G. T. P.
Charnock
,
P. J.
Hore
, and
I.
Kuprov
, “
Spinach—A software library for simulation of spin dynamics in large systems
,”
J. Magn. Reson.
208
,
179
194
(
2011
).
27.
J. H.
Freed
,
G. V.
Bruno
, and
C.
Polnaszek
, “
ESR line shapes for triplets undergoing slow rotational reorientation
,”
J. Chem. Phys.
55
,
5270
5281
(
1971
).
28.
J. R.
Norris
and
S. I.
Weissman
, “
Studies of rotational diffusion through the electron-electron dipolar interaction
,”
J. Phys. Chem.
73
,
3119
3124
(
1969
).
29.
F. B.
Bramwell
, “
ESR studies of phosphorescent corannulene; evidence for pseudorotation
,”
J. Chem. Phys.
52
,
5656
5661
(
1970
).
30.
A.
Blank
and
H.
Levanon
, “
Triplet line shape simulation in continuous wave electron paramagnetic resonance experiments
,”
Concepts Magn. Reson., Part A
25A
,
18
39
(
2005
).
31.
A.
Borel
,
R. B.
Clarkson
, and
R. L.
Belford
, “
Stochastic Liouville equation treatment of the electron paramagnetic resonance line shape of an S-state ion in solution
,”
J. Chem. Phys.
126
,
054510
(
2007
).
32.
D.
Kruk
,
J.
Kowalewski
,
D. S.
Tipikin
,
J. H.
Freed
,
M.
Mośinski
,
A.
Mielczarek
, and
M.
Port
, “
Joint analysis of ESR lineshapes and 1H NMRD profiles of DOTA-Gd derivatives by means of the slow motion theory
,”
J. Chem. Phys.
134
,
024508
(
2011
).
33.
E.
Berggren
and
C.
Zannoni
, “
Rotational diffusion of biaxial probes in biaxial liquid crystal phases
,”
Mol. Phys.
85
,
299
333
(
1995
).
34.
S.
Stoll
and
A.
Schweiger
, “
EasySpin: A comprehensive software package for spectral simulation and analysis in EPR
,”
J. Magn. Reson.
178
,
42
55
(
2006
).
35.
S.
Stoll
and
A.
Schweiger
, “
EasySpin: Simulating cw ESR spectra
,”
Biol. Magn. Reson.
27
,
299
321
(
2007
).
36.
Z.
Guo
,
D.
Cascio
,
K.
Hideg
,
T.
Kálái
, and
W. L.
Hubbell
, “
Structural determinants of nitroxide motion in spin-labeled proteins: Tertiary contact and solvent-inaccessible sites in helix G of T4 lysozyme
,”
Protein Sci.
16
,
1069
1086
(
2007
).
37.
K. A.
Earle
,
D. E.
Budil
, and
J. H.
Freed
, “
250-GHz EPR of nitroxides in the slow-motional regime: Models of rotational diffusion
,”
J. Phys. Chem.
97
,
13289
13297
(
1993
).
38.
D. A.
Varshalovich
,
A. N.
Moskalev
, and
V. K.
Khersonskii
,
Quantum Theory of Angular Momentum
(
World Scientific
,
1988
).
39.
R. N.
Zare
,
Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics
(
Wiley
,
1991
).
40.
C. F.
Polnaszek
and
J. H.
Freed
, “
Electron spin resonance studies of anisotropic ordering, spin relaxation, and slow tumbling in liquid crystalline solvents
,”
J. Phys. Chem.
79
,
2283
2306
(
1975
).
41.
M.
Mehring
,
Principles of High Resolution NMR in Solids
(
Springer
,
1983
).
42.
A.
Abragam
,
Principles of Nuclear Magnetism
(
Clarendon Press
,
1983
).
43.
R. R.
Ernst
,
G.
Bodenhausen
, and
A.
Wokaun
,
Principles of Nuclear Magnetic Resonance in One or Two Dimensions
(
Clarendon Press
,
1990
).
44.
G.
Binsch
, “
Unified theory of exchange effects on nuclear magnetic resonance line shapes
,”
J. Am. Chem. Soc.
91
,
1304
1309
(
1969
).
45.
Z.
Bai
,
J.
Demmel
,
J.
Dongarra
,
A.
Ruhe
, and
H.
van der Vorst
, “
Templates for the solution of algebraic eigenvalue problems
,” in
A Practical Guide
(
Society of Industrial and Applied Mathematics
,
2000
).
46.
W. H.
Press
,
S. A.
Teukolsky
,
W. T.
Vetterling
, and
B. P.
Flannery
,
Numerical Recipes: The Art of Scientific Programming
, 3rd ed. (
Cambridge University Press
,
2007
).
47.
S.
Stoll
and
A.
Schweiger
, “
An adaptive method for computing resonance fields for continuous-wave EPR spectra
,”
Chem. Phys. Lett.
380
,
464
470
(
2003
).
48.
J. S.
Hyde
,
M.
Pasenkiewicz-Gierula
,
A.
Jesmanowicz
, and
W. E.
Antholine
, “
Pseudo field modulation in EPR spectroscopy
,”
Appl. Magn. Reson.
1
,
483
496
(
1990
).
49.
A. R.
Edmonds
,
Angular Momentum in Quantum Mechanics
(
Princeton University Press
,
1957
).
50.
S.-T.
Lai
and
Y.-N.
Chiu
, “
Exact computation of the 3-j and 6-j symbols
,”
Comput. Phys. Commun.
61
,
350
360
(
1990
).
51.
M.
Pasenkiewicz-Gierula
,
W. K.
Subczynski
, and
W. E.
Antholine
, “
Rotational motion of square planar copper complexes in solution and phospholipid bilayer membranes
,”
J. Phys. Chem. B
101
,
5596
5606
(
1997
).
52.
G.
Della Lunga
,
M.
Pezzato
,
M. C.
Baratto
,
R.
Pogni
, and
R.
Basosi
, “
A new program based on stochastic Liouville equation for the analysis of superhyperfine interaction in CW-ESR spectroscopy
,”
J. Magn. Reson.
164
,
71
77
(
2003
).
53.
D.
Sezer
,
J. H.
Freed
, and
B.
Roux
, “
Multifrequency electron spin resonance spectra of a spin-labeled protein calculated from molecular dynamics simulations
,”
J. Am. Chem. Soc.
131
,
2597
2605
(
2009
).