Modeling spin relaxation in complex radical systems using MolSpin

Spin relaxation is an important aspect of the spin dynamics of free radicals and can have a significant impact on the outcome of their spin‐selective reactions. Examples range from the use of radicals as spin qubits in quantum information processing to the radical pair reactions in proteins that may allow migratory birds to sense the direction of the Earth's magnetic field. Accurate modeling of spin relaxation, however, is non‐trivial. Bloch–Redfield–Wangsness theory derives a quantum mechanical master equation from system‐bath interactions in the Markovian limit that provides a comprehensive framework for describing spin relaxation. Unfortunately, the construction of the master equation is system‐specific and often resource‐heavy. To address this challenge, we introduce a generalized and efficient implementation of BRW theory as a new feature of the spin dynamics toolkit MolSpin which offers an easy‐to‐use approach for studying systems of reacting radicals of varying complexity.


| INTRODUCTION
Spin effects on the reactions of radicals and other paramagnetic molecules are well established and have been studied since the 1960s.
Recent interest in this field has been fueled by the use of spincorrelated spin systems as promising qubits in quantum information processing, [1][2][3][4][5][6][7] and by the possible involvement of radical pairs as geomagnetic field sensors in migratory birds. [8][9][10][11] Due to the quantum nature of spin, the description of spin dynamics relies on fundamental quantum mechanical equations of motion, which are often challenging to solve for realistic systems. In particular, the accurate treatment of the interaction between the spins and their environment is non-trivial, but essential due to the spin relaxation it induces. This loss of spin coherence can severely limit the lifetime of particular spin states, thereby reducing the sensitivity of radical pair reactions to weak applied magnetic fields. [12][13][14] A specific example where spin relaxation may play a major role is in the mechanism by which night-migratory songbirds detect the direction of the Earth's magnetic field as an aid to navigation. This remarkable ability is thought to rely on a radical pair composed of radicals derived from flavin adenine dinucleotide (FAD) and a tryptophan residue (Trp) within a cryptochrome photoreceptor. 9 Formation of the FAD • À Trp • þ state of this protein is believed to trigger a cascade of biological reactions through spin-dependent processes which eventually constitute the magnetic compass sense. [9][10][11]15 Stochastic molecular motions of the protein cause spin interactions to fluctuate, resulting in loss of spin coherence and, in most cases, attenuation of the sensitivity to external magnetic fields. 12,16,17 Spin relaxation is often described phenomenologically, for example by using Lindblad operators. 18 However, such formalisms are not usually based on a microscopic physical picture of the actual motions and magnetic interactions of the radicals. The approach first derived by Redfield delivers a quantum master equation that describes the interaction between a spin system and its environment as a perturbation introduced through stochastic functions. 14,19 The construction of the Bloch-Redfield-Wangsness (BRW) relaxation matrix that describes the resulting spin relaxation can be a complex process. Specific implementations of BRW theory for radical pairs exist in popular spin dynamics software packages, 20,21 and as in-house codes written by research groups. 12,16,22 The latter, however, are usually specific for particular relaxation mechanisms and are not generally applicable to the broad range of situations in which spin relaxation is important.
MolSpin is an efficient modern software package that offers a variety of options for modeling spin systems ranging from spin dynamics to spin effects on the chemistry of complex reaction systems. 23 Its main features are a generalized approach leading to high flexibility, user-friendly extensibility allowing the inclusion of new modules, and intuitive input which can be further simplified by direct embedding in the web-based platform VIKING (Scandinavian Online Kit For Nanoscale Modeling). 23,24 The implementation of BRW theory into MolSpin presented here provides a toolbox for the efficient modeling of spin relaxation effects on spin-dependent chemical reactions. The algorithm requires the specification of physical parameters and implements the BRW formalism as generally as possible. We summarize the key features of Mol-Spin and our implementation of BRW theory, and discuss its current capabilities. The final part of the paper illustrates the use of BRW theory within MolSpin, highlighting the importance of accounting for spin relaxation in radical pair reactions. Figure 1 illustrates the key features of MolSpin; a complete account can be found in the MolSpin documentation. 25 The main goal of the program is to provide a generalized, easy-to-use framework that allows the investigation of arbitrary spin systems and processes involving spins ranging from qubits in quantum computing to radical pairs in biological systems. It is possible to simulate complex spin-dependent reactions involving multiple coupled spin systems. A variety of linear and bilinear spin interactions can be modeled. 23 MolSpin offers high flexibility concerning extensibility and the inclusion of new modules, for example the dynamical evolution of input parameters which permits complex time-dependent applied magnetic fields or reaction kinetics to be specified. Based on the C++ library Armadillo, 26

| FEATURES OF MOLSPIN
where b ρ t ð Þ is the time-dependent density matrix of the spin system, b H t ð Þ is the Hamilton operator that includes all spin interactions, b K is the reaction operator and Â, Â ½ denotes the commutator of two operators. Here, and in the following, energies have the dimensions of angular frequency, that is, ℏ ¼ 1. The interactions within b H are implemented in MolSpin in the following generalized forms 23 : where b H SS i ð Þ is the single-spin interaction of spin i, described by the spin operator b S ! i , with an external field vector F ! coupled via a tensor A. μ B is the Bohr magneton, α is a constant factor, and g i is the g-tensor of spin i which can be specified for each spin separately. A common single-spin interaction is the Zeeman interaction of a particle with an external magnetic field. In contrast, b H DS i, j ð Þ describes the interaction between two spins i and j with a coupling tensor A. Two-spin interactions, include dipolar, exchange and hyperfine interactions. MolSpin uses the formalism in Equation (2) to guarantee versatile usability and the ability to solve specific physical problems (see Manual). 23 More specific interactions such as spin-phonon coupling, zero-field splitting and crystal-field splitting are not implemented in the current version of MolSpin. 31,32 The object-oriented framework of the program, however, allows for a straightforward expansion of the code in the future.
Once the time-dependent density matrix b ρ t ð Þ is known, the probability that the system is in a specific state j ii, p i t ð Þ, at a time t is obtained from: Here b P i ¼j iihi j is the projection operator for state j ii and Tr() is the matrix trace. The reaction term in Equation (1) in MolSpin can be defined in two ways. The so-called Haberkorn model assumes the Here k i is the rate constant of the spin-dependent reaction of a spin state of interest (e.g., singlet or triplet state in a radical pair) described by a specific projection operator and the curly braces Â, Â f g denote the anti-commutator. n is the total number of reactions included. This approach is widely used and directly derived from second-order perturbation theory. 34 One feature of the Haberkorn model is the loss of trace preservation of the density matrix. This can be avoided using a Lindblad dissipator, which leads to 35 where b Q i are projection operators of the singlet or triplet radical states onto the associated product states. The Lindblad dissipator is also commonly used to describe phenomenological spin relaxation processes induced by the interaction of a spin system with its environment. 23 Spin-dependent reactions play key roles in biological processes, for example, avian magnetoreception [9][10][11]15 or cellular damage through radicals. [36][37][38] The observable characteristic in spin-dependent reactions is the product yield ϕ i t ð Þ, [9][10][11]15 which for a specific spin state can be defined as: Several numerical methods could be used to solve Equation (6) depending on the nature of the spin system. 25 In many cases, the upper limit of the integral in Equation (6) is set to infinity.
MolSpin includes Lindblad dissipators both as reaction operators and to model spin relaxation phenomenologically. However, when the interaction between the environment and the spin system is considered in terms of a real microscopic physical model of spin interactions, the second-order stochastic theory developed by Redfield becomes F I G U R E 1 Illustration of features of MolSpin for studying arbitrary spin systems. Spin relaxation has been introduced as a new feature, as detailed in this work. MD, molecular dynamics.
applicable. 14,19,39 BRW theory yields a master equation for the iso-

| BASICS OF BRW THEORY
A comprehensive derivation of BRW theory can be found in textbooks and publications. 19,39,40 A brief description of the essential formalism is presented here, emphasizing the numerical implementation in MolSpin. BRW theory describes the interaction between a spin system S and the environment or bath B. Mathematically speaking, the total Hamiltonian b H of the spin system and the bath can be factorized as 39 : Here, b H S is the Assuming that the ensemble-averaged expectation value of b H I is zero, that the bath is unaffected by the system dynamics, and that the time-dependent perturbation has a significantly lower amplitude than the static Hamiltonian b H S , the evolution of the system can be formulated in the interaction picture, up to second-order in b H I , as 39 : where b ρ S t ð Þ is the density operator of the spin system, b ρ B is the density operator of the bath, Tr B is the trace over the bath states. In the following all derivations are in the interaction picture unless otherwise stated. Equation (8) rests on the assumption (often called the Born approximation) that the complete density operator b ρ t ð Þ (system and bath) can be factorized as 39 : Equation (8) is non-Markovian insofar as the change in the density matrix at a time t depends on the previous time instances t 0 < t. A further simplification of Equation (8) postulates a Markovian behavior of the spin system (Markov approximation). This is realized, by first assuming that b ρ S t ð Þ changes slowly with t, so that Equation (8) can be rewritten 39 : Although this equation is local in time, it is still not Markovian due to the time evolution of the density matrix b ρ S t ð Þ, which depends on the explicit choice of the initial state b ρ S 0 ð Þ at t ¼ 0.
Equation (10) can be formulated in Markovian form if the substitution t 0 ! t À t 0 is made and the upper limit of the integral is extended In general, the interaction operator b H I can be decomposed into a sum of direct products of operators b S α describing the spin system and in which b S α and b B α can be chosen arbitrarily as long as b Equation (11) can now be rewritten by formulating the commutators where Equation (13) can now be formulated in matrix form in the eigenbasis of b H S : where ω mn ¼ ε m À ε n and ε m is the eigenvalue of b H S corresponding to the eigenstate j mi. Equation (13) can now be rewritten in the Schrödinger picture as 39 (the subscript S is omitted here): The first term on the right hand arises from the static Hamiltonian b H S while the integral accounts for the spin relaxation induced by the system-bath interactions. Equation (16) can be rewritten in terms of the spectral densities J αβ ω mn ð Þ 39 : which represent the strength of the system-bath coupling at frequency ω mn . Equation (16) thus assumes the form 16 : where R abcd is given by: Including the reaction operator and transforming Equation (18) into Liouville space, one obtains the following equation of motion for the spin density operator:

| COMPUTATIONAL FEATURES OF BRW THEORY IN MOLSPIN
Depending on the problem of interest different spin operators S α are most suitable as a basis for the construction of R abcd in Equation (19).
For instance, when relaxation results from rotational motion, a spherical tensor basis b T for the interaction of interest is recommended, while for more general problems, the usual Cartesian spin operators b S x , b S y and b S z may be used. 16,40,41 MolSpin provides both basis sets. The spherical tensor formalism relies on the decomposition of linear and bilinear interactions into spatial tensors a k q and spherical tensors b T k q of rank k with q ¼ Àk, ::,0,…, Here usually zero due to A being symmetric. 40 In MolSpin, the spatial and spherical tensors are as follows:  31 In the secular approximation, specific matrix elements of R abcd are set to zero when the difference ω ab À ω cd is greater than a set cutoff value. 31 Even though the use of a secular approximation can lead to a speedup, the selection of the matrix elements to be neglected is non-trivial and requires care in order to produce meaningful results.
Often the correlation functions in Equation (17) are simple exponential decays of the form: Here τ c is the correlation time and g 0 ð Þ is the amplitude of the correlation. By integrating Equation (26) according to Equation (17), the spectral density J ω ð Þ then reads: The imaginary term on the right-hand side of the equation is the dynamic frequency shift, which is often neglected when studying relaxation phenomena due to its small contribution. MolSpin has the option to neglect this term.
In many realistic situations, such as the motion of atoms in proteins, the spectral density cannot accurately be described by a single exponential correlation function. 12 A more realistic description is achieved by constructing the spectral density from a linear combination of correlation functions 12 : where g q 0 ð Þ and τ q c are the amplitudes and correlation times, respectively. 12

| Radical pair mechanism in avian magnetoreception
The light-induced formation of a radical pair comprising the cofactor FAD and a Trp residue in a cryptochrome photoreceptor protein is believed to play a key role in sensing the direction of the Earth's magnetic field in migratory birds. 9,11,15,34,43 The singlet and triplet states of the radical pair undergo spin-selective reactions which are thought to trigger a cascade of spin-independent biochemical reactions that constitute the signal transduction pathway. Only the singlet state can additionally recombine to the ground state of the protein, a reaction that is spin forbidden for the triplet state.
The quantum yield of the singlet recombination reaction ϕ S is generally used to assess the effects of the Earth's magnetic field on the electron spins in the radical pair in cryptochrome 9,10 which also interact with nuclear spins via hyperfine coupling. 12 The thermal motions of the protein perturb these interactions leading to loss of spin coherence. 12,16 Figure 2 introduces a simple model radical pair, containing a single nitrogen nucleus ( 14 N, spin I ¼ 1), with a system-environment coupling that leads to spin relaxation via stochastic fluctuations in the hyperfine coupling tensor.
The singlet probability of the radical pair was calculated as a function of time for different rotational correlation times. The hyperfine interaction parameters for the nitrogen nucleus were those of the FAD • À radical as studied by Hiscock et al. 43 (see SM S2) and are modulated by isotropic rotational diffusion. Note, that the model is used for illustrative purposes and is not suitable for an accurate description of the radical pair within cryptochrome.
The magnetic field B ! had strength 50 μT (typical of the geomagnetic field) and made an angle of 45 with the z-axis (see Figure 2B and Section S2). The spherical tensor basis was used here to describe the modulation of the hyperfine interaction of the nitrogen with one of the electrons of the radical pair. The spin system's initial state was chosen to be purely singlet. (via the anisotropy of the hyperfine couplings) which could be used to sense the direction of the field, that is, to provide a radical-pair-based sensor. 12 Figure 3 illustrates the magnetic field effects arising in a rather simplified version of the radical pair comprising FAD • À and Trp • þ radicals in cryptochrome. This case study aims to show the difference between a radical pair with a purely static Hamiltonian and one in which the spins relax. The magnetic field in this example is rotated around the y-axis in the zx-plane, and ϕ S is calculated according to Equation (6) with and without spin relaxation. The upper limit of Equation (6) was set to infinity. Four nuclei (see Figure 3B) were included in the FAD • À radical, with hyperfine parameters taken from Hiscock et al. 43 (see SM S3). None of the nuclear spins in Trp • þ was included. The correlation time for isotropic rotational diffusion was arbitarily chosen to be 5 ps, where only the N 5 hyperfine interaction was modulated for illustration purposes. The N 5 hyperfine interaction was found previously to have a significant impact on the quantum yield. 12 For illustrative purposes, it was assumed that the singlet and the triplet states react to form products with equal rate constants (10 4 s À1 ). 43 Without spin relaxation, the quantum yield ϕ S decreases significantly as θ is changed from zero to 90 as a result of the combined effects of the nitrogen hyperfine interaction and the external B ! -field. 43 Introducing spin relaxation for just the hyperfine interaction of N 5 leads to a notable decrease in the singlet reaction yield ϕ S for every magnetic field orientation, as seen in Figure 3A. The modulation of this interaction promotes thermalization toward equally populated states, leading to the decrease in ϕ S . The tumbling in this scenario is fast enough to average out the anisotropy of N 5 so that the other three nuclei become responsible for the θ-dependence of the red line in Figure 3A.  44,45 ; the procedure is discussed in detail by Kattnig et al. 12 The amplitudes g 0 ð Þ q and the corresponding correlation times τ q c can be directly used in MolSpin (SM, S4), which in turn computes the spectral density automatically. Figure 4C illustrates the magnetic field effect obtained as in Figure 3, but using the spectral density calculated from Figure 4A. The parameters are found in SM, section S4. As in the previous example ( Figure 3) motion-induced spin relaxation leads to a decrease of the overall singlet quantum yield. Again, the anisotropy of the N 5 hyperfine interaction is mostly averaged out through the rapid motion of F I G U R E 3 (A) Magnetic field effects as a function of the direction of the magnetic field vector B ! around the y-axis characterized by the angle θ. A radical pair comprising FAD • À and Trp • þ radicals, with and without spin relaxation, was considered. Four nuclei (N 5 , H b1 , N 10 , H 6 ) were included to model the hyperfine interactions. The reaction rate constant k i was set to 10 4 s À1 for all states, while the yield was calculated according to Equation (6) with the upper limit of the integral set to infinity. Only the hyperfine interaction with N 5 (red) was subject to spin relaxation using spatial and spherical tensors and a correlation time of 5 ps. (B) Schematic illustration of the spin system and the included nuclei.
F I G U R E 2 (A) Time-evolution of the singlet state population of a pair of electron spins one of which is coupled to a spin-1 nitrogen nucleus. The hyperfine interaction gives rise to spin relaxation assuming isotropic rotational diffusion with different correlation times (5, 50, and 500 ps). The time-dependence for the static case, with no relaxation, is included for comparison. (B) Illustration of the radical pair (light blue) in relation to the magnetic field B ! and the nitrogen hyperfine interaction, modulation of which leads to spin relaxation. the FAD • À which promotes evolution of the radical pair toward equally populated states. This example shows that it is possible to include complex motion-induced fluctuations of spin interactions using MolSpin with the procedure developed by Kattnig et al. 12

| Spin relaxation in formic acid formation
Another situation where spin relaxation may play an important role is in the conversion of CO 2 into organic molecules, which offers a way to reduce the concentration of atmospheric CO 2 while at the same time producing useful chemicals. 46,47 Our final example involves the formation of formic acid through electrocatalytic reduction of CO 2 using a tin electrode in the presence of a magnetic field of up to 900 mT. 13 It was hypothesized that the magnetic field-induced increased yield of formic acid arises from the interconversion of the electronic singlet and triplet states of a radical pair comprising CO • À 2 and H • radicals. Figure 5 shows a schematic representation of the mechanism behind the proposed reaction.
The recombination of the radical partners is only allowed for singlet radical pairs due to the spin selection rules, leading to the F I G U R E 4 (A) Autocorrelation function of the a zz (N 5 ) hyperfine component modulated by changes in the dihedral angle Ω and threedimensional librational motions of FAD • À . The blue line shows the extracted autocorrelation function using a 436 ns MD trajectory from Grüning et al. 44 The orange line shows a multi-exponential fit using 50 exponential functions according to Equations (26)-(28) with τ q c logarithmically distributed between 0.5 and 91.5 ns. Following the approach in Reference [12], the auto-correlation function was fitted up to its first zerocrossing. (B) Illustrative representation of the included nuclei, the dihedral angle Ω and librational motion around the y-axis. The magnetic field B ! (50 μT) is rotated in the xz-plane through an angle θ. The radical pair starts in a singlet state and reaction rate constants for singlet and triplet states were 10 6 s À1 , Reference [44]. Four nuclei (N 5 , H b1 , N 10 , H 6 ) were included to model the hyperfine interactions. All parameters for hyperfine interactions and autocorrelation functions were taken from Grüning et al. 44 (C) Magnetic field effect without relaxation (black) and with spin relaxation (red) through the fluctuating hyperfine component a zz of N 5 . Singlet quantum yields were calculated according to Equation (6) with the upper limit of the integral set to infinity. The corresponding input file can be found in the SM S4.
F I G U R E 5 Proposed spin-dependent mechanism of the reaction of H • and CO • À 2 to form formic acid due to Pan et al. 13 The rate constants k D and k s refer to the diffusive separation of the two radicals and their reaction to produce formate, respectively. production of the diamagnetic HCO À 2 molecule with a reaction rate constant k s . Competing with this process, the radicals may separate by translational diffusion (to undergo non-productive follow-up reactions) characterized by the diffusion rate constant k D .
An earlier study 48 found that the g-tensor anisotropy spin relaxation mechanism (rotational modulation of the anisotropic electron Zeeman interaction) was essential to reproduce the experimental observations. The relaxation was described by a relaxation supero- where S ! Aj denotes the spin operator of the electron in the CO • À 2 radical, b E ! is the identity operator, g jj are the principal g-tensor components (g xx ¼ 2:0032, g yy ¼ 1:9975, g zz ¼ 2:0014), g iso ¼ 2:0007 (the isotropic g-value of CO • À 2 ), and ω 0 is the electron Larmor frequency defined as ω 0 ¼ g iso μ B B ! ℏ . The spectral densities J ω ð Þ have the form in Equation (27) neglecting the dynamic frequency shift. This relaxation superoperator allowed the authors of Reference [48] to explain the experimental observations of Reference [13] for 30-900 mT magnetic fields by calculating the singlet quantum yield ϕ SQY at different magnetic field strengths with respect to the radical pair 48 : where the singlet yield ϕ S is defined in Equation (6) assuming the rate constants k s ¼ 10 7 s À1 and k D ¼ 10 4 s À1 . To simulate the g-tensor anisotropy relaxation mechanism in MolSpin the spherical tensor operator basis of rank 2 was used to construct the relaxation matrix starting from the interaction Hamiltonian 20,41 : b where D 2 m,n are the Wigner functions, and Ax ¼ 3a zz À a xx þ a yy þ a zz ð Þ and Rh ¼ a xx À a yy are directly associated with the spatial tensors a 2 0 and a 2 AE2 in Equations (22)-(25). The correlation function can be constructed as: Due to the properties of Wigner functions, their product  48 which agree with the experiments of Pan et al. 13 The dotted lines show results without spin relaxation and the solid lines show reaction yields including relaxation caused by rotational modulation of the Zeeman interaction using MolSpin.
F I G U R E 7 Computational speedup of the BRW algorithm obtained using multiple CPU threads. As a test system a radical pair comprising 1-5 spin-1 2 nuclei was considered, in a 500 mT magnetic field. All interactions were subject to spin relaxation. R ele denotes the number of elements of the R-matrix. The scaling of t #1 with respect to the number of spins included can be found in Table 1. The present results obtained with the use of BRW theory are identical to the earlier ones. 48 Without relaxation, the singlet yields ϕ SQY are drastically different and cannot satisfactorily reproduce the experimental observations. 13,48 In all simulations where the H and C isotopes were considered, relaxation effects appear to be more significant when compared with the 1 H, 12 C system.

| PERFORMANCE
One of the major disadvantages of BRW theory is the exponential growth of the R-matrix, Equation (19), with the size of the Hilbert space. The demand for memory rises rapidly and the time required to calculate all matrix elements becomes a challenge when Equation (20) is solved serially. To overcome the computational bottlenecks, the algorithm implemented in MolSpin was parallelized using OpenMP, 49 on top of the already parallelized Armadillo library, 26 which permits an additional reduction in the computational time required for the construction of the relaxation matrix. Due to the numerical structure of Equation (19) it is natural to parallelize each operator pair segment S α S β , adding each term separately to R abcd . Figure 7A illustrates the reduction in the computational time required to solve a quantum yield calculation involving a spherical tensor Redfield matrix with different number of matrix elements R ele as a function of the number of computational threads. To compare the performance of MolSpin when constructing R of different sizes, the computational speedup was calculated as the ratio of the time required for a serial calculation (t #1 ) to the time required for calculations using more than one thread (t #n ).
The calculations were performed on a workstation with 504 GB RAM DDR4, Intel Xenon Gold 6226 CPU @ 2.90 GHz. The benchmark spin system in this scenario was a radical pair containing up to five H nuclei (spin- 1 2 ). Especially for the larger systems with up to $268 million matrix elements, the implemented parallelization gives a significant performance boost for up to 18 threads. The largest R-matrix with $268 million elements only requires 3057 s = 0.85 h instead of 26,738 s = 7.43 h (see Table 1). Thus, it is possible to efficiently calculate large R-matrices using MolSpin.   Note: As a test system, a radical pair interacting with 1-5 spin-1 2 H nuclei was considered, exposed to a magnetic field of 500 mT pointing in the z-direction. All interactions have been subject to spin relaxation. Calculations were performed on a wokstation with 504 GB RAM DDR4, Intel(r) Xenon(R) Gold 6226R CPU @ 2.90 GHz.