Spin–phonon coupling and magnetic relaxation in single-molecule magnets

Abstract

Electron–phonon coupling is important in many physical phenomena, e.g. photosynthesis, catalysis and quantum information processing, but its impacts are difficult to grasp on the microscopic level. One area attracting wide interest is that of single-molecule magnets, which is motivated by searching for the ultimate limit in the miniaturisation of binary data storage media. The utility of a molecule to store magnetic information is quantified by the timescale of its magnetic reversal processes, also known as magnetic relaxation, which is limited by spin–phonon coupling. Several recent accomplishments of synthetic organometallic chemistry have led to the observation of molecular magnetic memory effects at temperatures above that of liquid nitrogen. These discoveries have highlighted how far chemical design strategies for maximising magnetic anisotropy have come, but have also highlighted the need to characterise the complex interplay between phonons and molecular spin states. The crucial step is to make a link between magnetic relaxation and chemical motifs, and so be able to produce design criteria to extend molecular magnetic memory. The basic physics associated with spin–phonon coupling and magnetic relaxation was outlined in the early 20th century using perturbation theory, and has more recently been recast in the form of a general open quantum systems formalism and tackled with different levels of approximations. It is the purpose of this Tutorial Review to introduce the topics of phonons, molecular spin–phonon coupling, and magnetic relaxation, and to outline the relevant theories in connection with both the traditional perturbative texts and the more modern open quantum systems methods.

---

Jon G. C. Kragskow

Jon Kragskow was awarded his MChem (Hons) in 2018 and PhD in 2022 from The University of Manchester. His PhD, supervised by Prof. Chilton, focussed on understanding the relaxation dynamics of lanthanide single molecule magnets. He is currently a postdoctoral research associate with Dr Elizaveta Suturina at The University of Bath, investigating hyperfine coupling in Nuclear Magnetic Resonance spectroscopy.

Andrea Mattioni

Andrea Mattioni studied Physics at The University of Trieste, where he received his MSc in 2016. He obtained his PhD from Ulm University, working on energy transfer dynamics in molecular light-harvesting systems in the group of Prof. Martin Plenio and Prof. Susana Huelga. In 2020 he joined the Chilton group at The University of Manchester as a postdoctoral research associate. His focus is to elucidate the impact of dissipation and strong spin–phonon coupling on the relaxation and decoherence dynamics of molecular spin systems.

Jakob K. Staab

Jakob Staab obtained his BSc degrees in Chemistry and Biomolecular Engineering at TU Darmstadt, and his MSc degree in Theoretical Chemistry at Uppsala University before starting his PhD in the Chilton group in 2020. His current work focuses on the development of methods to compute spin–phonon couplings in molecular magnetic materials as well as analysis tools for the electronic structure of d- and f-block metal complexes.

Daniel Reta

Daniel Reta is an Ikerbasque Research Fellow at The University of the Basque Country and Donostia International Physics Centre. He obtained his PhD from The University of Barcelona in 2016 working with Prof. Francesc Illas and Prof. Iberio Moreira. He then moved to The University of Manchester as a postdoctoral research associate, first with Prof. Nikolas Kaltsoyannis and subsequently with Prof. Nicholas Chilton, where he worked on the description of spin dynamics in single-molecule magnets. His current research focuses on the magnetic and reactivity properties of organic radicals, merging computational and spectroscopic approaches.

Jonathan M. Skelton

Jonathan Skelton is a UKRI Future Leaders Fellow at The University of Manchester. He completed his undergraduate and PhD degrees at The University of Cambridge in 2013, then was a postdoctoral research associate at The University of Bath with Prof. Aron Walsh before joining The University of Manchester in 2018. His research group specialises in modelling the structural dynamics and thermal physics of solids, with a particular focus on materials such as thermoelectrics for renewable-energy technologies.

Nicholas F. Chilton

Nicholas Chilton is a Professor of Computational and Theoretical Chemistry and a Royal Society University Research Fellow at The University of Manchester. He obtained his BSc Adv Hons at Monash University in 2011, under the supervision of Prof. Keith Murray, and his PhD from The University of Manchester in 2015 under the supervision of Prof. Richard Winpenny and Prof. Eric McInnes. He was a Ramsay Memorial Research Fellow from 2016–2018 and then a University of Manchester Presidential Fellow from 2018–2019. His research focusses on the spectroscopic and magnetic investigation of magnetic molecules, with a particular focus on multiconfigurational computational methods for molecular electronic structure problems, including the spin dynamics of molecular spin systems.

---

Key learning points

1. The nature of phonons in molecular crystals.

2. How spin–phonon coupling arises in molecules.

3. How the spin–phonon interaction drives magnetic relaxation.

4. How molecular spin dynamics and magnetic relaxation can be modelled with different levels of approximation.

5. How to perform first-principles calculation of the spin–phonon coupling.

1 Introduction

For the last thirty years, many inorganic, physical and computational chemists, along with experimental and theoretical physicists, have been captivated by the discovery,^1,2 study and development of single-molecule magnets (SMMs). SMMs are molecules that possess a doubly-degenerate electronic ground state, the components of which have large uniaxial magnetic moments in opposing directions that can represent binary 1 and 0.^3,4 The other low-lying excited electronic states generally have their magnetic moments oriented more towards the equatorial plane (i.e. perpendicular to the axis defined by the ground doublet); this energy difference between different orientations of the magnetic moment is the origin of magnetic anisotropy. The best-performing SMMs to-date are [DyCp^Me₅Cp^iPr₅][B(C₆F₅)₄] and [Cp^iPr₅DyI₃DyCp^iPr₅] (Cp^Me₅ = C₅(CH₃)₅ and Cp^iPr₅ = C₅(CH(CH₃)₂)₅), which both show open magnetic hysteresis up to 80 K;^5,6 the latter features strong magnetic interactions between the metal ions, which offers a new route to achieving improved performance.⁷

There have been numerous recent reviews on SMMs,^7–10 and it is not the purpose of this Tutorial Review to discuss the requisite ingredients for their design; indeed, this Review assumes a basic understanding of such matters.^4,11 Rather, here we focus on the theory and calculation of magnetisation dynamics in SMMs as arising from spin–phonon coupling. While there are many contemporary theoretical works in this area,^12–18 and some excellent reviews of the traditional methods,¹⁹ there exists a gulf between the modern parlance of these physics-based texts and the original works in the early to mid 20^th century to which the field constantly reflects.^20–23 Hence, this work intends to contextualise modern molecular spin–phonon coupling theory with a backdrop of the original works, and to provide an explanation suitable for newcomers to the field – including the basics of lattice dynamics. Herein we focus on discussion of magnetic relaxation in lanthanide (Ln) based SMMs,^8,24 which invariably possess electronic structures with the characteristic profile in Fig. 1;^7,9 that is, interelectronic repulsion dominates so that the Hund's rule ground state is a well isolated Russell-Saunders term, which is then split by spin–orbit coupling to give a well-isolated manifold of total angular momentum J, which is then split by crystal field (CF, or ligand field) interactions into linear combinations of m_J states.^7,25 The theories discussed herein are equally applicable to spin–phonon coupling in molecules with other metal ions, but the importance of the various mechanisms changes when the energy scales and electronic state multiplicities differ. Note also that while we use terms such as “spin–phonon”, “spin system”, “spin Hamiltonian”, etc., “spin” should be interpreted as a general total angular momentum.

Fig. 1 Two main mechanisms of magnetic relaxation arising from spin–phonon coupling in single-molecule magnets, Orbach (blue) and Raman (red), illustrated for the J = 15/2 multiplet of a Dy³⁺ SMM with a strong uniaxial crystal field.

SMMs are zero-dimensional superparamagnets and there is no phase transition associated with their magnetic memory, unlike for bulk ferromagnets. Thus, their magnetic dynamics are determined solely by their magnetic relaxation rates, the timescale on which an ensemble returns to thermodynamic equilibrium;^3,4 this process is often called spin-lattice relaxation and given the symbol T₁, but herein we will use the symbol τ. These rates are determined experimentally by measuring the time dependence of the magnetisation, often as a function of an external parameter such as temperature or magnetic field.^9,26 Although magnetic relaxation is an ensemble property, the underlying rates are determined by how the individual molecules exchange energy with their environment. Thus, the process of a sample returning to equilibrium is the result of a great number of possible spin–phonon interactions and relaxation pathways. There are three common mechanisms typically invoked in the interpretation of relaxation rate data (Fig. 1):¹⁹ (i) the Orbach process, in which relaxation occurs via a series of single spin–phonon absorption and emission events up and over the energy barrier; (ii) the Raman process, in which two phonons interact simultaneously with the SMM (i.e. a phonon is scattered in a collision with the molecule), where usually the transition between the ground doublet states is of most importance; and (iii) quantum tunnelling of the magnetisation (QTM) where relaxation occurs directly between the ground doublet states and does not involve absorption or emission of phonons (not shown in Fig. 1); QTM will not be discussed further herein.

This Tutorial Review will start with a discussion of the phonon degrees of freedom, followed by the basics of spin–phonon interactions, and then will discuss in detail how these ingredients can be combined to calculate magnetic relaxation rates under single-phonon (Orbach) and two-phonon (Raman) paradigms.

2 Phonons

The first point of discussion is the vibrational modes themselves, which are frequently referred to as the “bath” in the literature on open quantum systems, or the “lattice” in the context of physical chemistry. We begin with a discussion of the essential features of the vibrations of gas-phase molecules. An isolated molecule in the gas phase with n_a atoms has 3n_a degrees of freedom, made up of the three Cartesian directions x, y and z for each atom, which combine together to form translations, rotations and vibrations. The combinations of atomic displacements, and the associated vibrational frequencies (energies), are obtained within the harmonic approximation as the eigenvalues and eigenvectors of the mass-weighted Hessian matrix H, which we refer to here as the “dynamical matrix” D:where r^α_κ are the atomic degrees of freedom, the indices κ, κ′ and α, β label the atoms and Cartesian directions, respectively, and m_κ are the atomic masses.

D is a 3n_a × 3n_a matrix, and diagonalisation yields 3n_a squared frequencies ω² (eigenvalues) and associated mass-weighted displacements W in the form of 3n_a-component vectors (eigenvectors). Fig. 2 shows the energies and displacements obtained by diagonalising the dynamical matrix of CO₂, calculated with density–functional theory (DFT; see ESI,† for details). Of the 3n_a combinations of the degrees of freedom, three combinations are rigid translations, and two or three are rigid rotations depending on whether the molecule is linear or non-linear. These collective motions of the atoms do not perturb the intramolecular interactions, and therefore have ω² = 0. The remaining 3n_a − 5 or 3n_a − 6 combinations are the vibrations, which are mutually-orthogonal relative motions of the atoms that conserve the centre of mass but have an energy change associated with them (i.e. ω² > 0). For CO₂, there are three translations and two rotations with ħω = 0, and four vibrational modes comprising a doubly-degenerate bend (ħω = 632 cm⁻¹) and symmetric and antisymmetric stretches (ħω = 1315 and 2357 cm⁻¹, respectively).

Fig. 2 Energies and displacement vectors for the translations, rotations and vibrations of the CO₂ molecule, obtained by diagonalising the dynamical matrix D defined in eqn (1), from DFT calculations (see ESI,† for details). The energies are obtained by taking the square root of the eigenvalues ω² and multiplying by the reduced Planck constant ħ, and the displacements are obtained from the eigenvectors W using eqn (3) with a unit amplitude Q = 1.

Within the harmonic approximation, the potential energy of the vibrational modes is quadratic in the displacement amplitude Q according to:

The relation between Q and the Cartesian displacement of the κ^th atom is given by:where W^α_κ is the component of the eigenvector for the atom along the direction α. Solution of the Schrödinger equation for the kinetic energy and the potential energy function in eqn (2) yields allowed energy levels quantised into units of ħω:where v is the vibrational quantum number; v = 0 corresponds to the energetic ground state, and the residual of energy in this state is termed the “zero-point” energy of the vibrational mode. Fig. 3 shows the potential energy as a function of amplitude for the two stretching modes in CO₂, together with the allowed harmonic energy levels.

Fig. 3 Potential energy as a function of amplitude E(Q) for the two stretch modes in CO₂ determined from eqn (2): (a) the symmetric stretch with ħω = 1315 cm⁻¹, and (b) the antisymmetric stretch with ħω = 2357 cm⁻¹. The E(Q) for both modes are shown over a range of Q = ±0.5, and the energies are given in units of ħω. On each plot, the black lines show the harmonic energy levels E_v determined from eqn (4). Note that as the energies are given in units of ħω, the implication is that modes with small ħω will, in general, have larger displacements |Q| for a given vibrational quantum state than will modes with large ħω.

We now consider the generalisation of the harmonic approximation to periodic crystals. Vibrations in crystalline solids, termed “phonons”, take the form of travelling waves of atomic displacements that propagate through a crystal. As in molecules, the relative atomic motion has an associated frequency ω, and the wavepacket propagates through the crystal with a defined wavelength λ and velocity ν.

In principle, the phonon modes could be obtained by constructing and diagonalising a dynamical matrix in the form of eqn (1) for the whole crystal. However, since n_a approaches Avogadro's number (N_A ≈ 6.022 × 10²³), such a “brute force” approach is not feasible. Instead, the periodicity of the crystalline phase can be exploited to view a (macroscopic) crystal as an infinite repeating array of primitive unit cells with n_a atoms each. According to the Bloch theorem, a wavefunction Ψ(r) in such a periodic system can be expressed as the product of a cell-periodic function u(r) and a plane-wave exp(iq·r):²⁷

Ψ(r) = u(r) × exp(iq·r) (5)

where q is a wavevector defined in the reciprocal space (first Brillouin zone, BZ) of the crystal. We note that wavevectors are conventionally discussed in “reduced” units of fractions of the reciprocal lattice vectors, which are related to the q used herein as q = 2π(q^R₁a* + q^R₂b* + q^R₃c*), where q^R_n are the components of the reduced vector q^R, and a*, b* and c* are the reciprocal lattice vectors. In the context of a phonon, the u(r) correspond to a set of relative atomic displacements in the primitive cell, and the wavevectors q in the plane-wave terms specify the propagation direction and wavelength of the phonons, and define how the local atomic displacements are modulated across the other unit cells in the crystal (Fig. 4).

Fig. 4 Schematic illustration of how the Bloch theorem defined in eqn (5) generates four of the phonon modes of a monoatomic 2D lattice. The displacement of the atom in the primitive cell defines the cell-periodic function u(r) (blue), which is replicated across the other unit cells in the lattice with a phase pattern defined by plane-wave terms exp(iq·r) (red). In (a) and (b) the wavevector q is oriented along the real-space x direction with two different wavelengths corresponding to different modulation periods. In (c) and (d), the wavevector is aligned along the y and xy directions, respectively, with the same wavelength as in (b).

Application of the Bloch theorem yields the modified dynamical matrix D(q) given by:

where the indices l, l′ label different crystallographic unit cells. Diagonalising a given D(q) yields 3n_a squared frequencies ω_qj² and corresponding mass-weighted displacements W_qj associated with the wavevector q, where j is the phonon band index and runs from 1 to 3n_a.

Use of the Bloch theorem simplifies the problem of determining the phonons in an infinite periodic crystal to diagonalising a set of D(q) at a formally infinite number of wavevectors q. However, the ω_qj and W_qj typically do not vary significantly for small changes in q, allowing physical properties requiring integrals over the BZ to be replaced by a finite, discrete sum. Furthermore, the phonon modes at certain wavevectors can be related by crystal symmetry operations, and most properties can therefore be obtained by considering only the subset of q that lie within the “irreducible” part of the BZ; for crystals in high-symmetry spacegroups, this is a fraction of the full BZ.

To discuss the phonon spectra of solids, two quantities are commonly used: the phonon dispersion (band structure) and the density of states (DoS). The dispersion ħω_j(q) shows the evolution of the 3n_a phonon energies with the wavevector q along a specific path through the BZ, and is typically displayed with a path visiting the set of high-symmetry q-vectors comprising the zone-centre q^R = (0, 0, 0) (Γ) and the zone-boundary points, the number, coordinates and labels of which are defined by the crystallographic spacegroup. The DoS g(ħω) shows the relative number of modes as a function of phonon energy integrated over the BZ:

where Ω is the volume of the BZ, δ is the Dirac delta function (δ(x) = 0 for x ≠ 0 and ), and N is the number of wavevectors included in the summation approximating the integral. We note that the DoS is normalised so that it integrates to the number of modes in the primitive cell, i.e.. It is also common to use the eigenvector coefficients to quantify the contributions from different atoms in order to construct an atom-projected partial DoS showing how atoms contribute to modes in different energy ranges.

The dispersion and DoS for NaCl (calculated using DFT, see ESI†) highlight several features common to phonon spectra (Fig. 5). There are three “acoustic modes” that correspond to concerted motion of all of the atoms in the primitive cell in the same direction with the same amplitude; at q = Γ, this motion is fully in-phase for all unit cells in the crystal, resulting in three rigid translations of the entire crystal with ħω_j = 0 (one of the acoustic modes for NaCl is doubly degenerate along the dispersion path). [We note that, unlike for gas-phase molecules, rigid rotations of crystals are not defined.] At the zone-boundary wavevectors, here and the atomic motion in neighbouring unit cells is fully out-of-phase along particular real-space directions, and therefore the energy of the acoustic modes increases sharply away from Γ with an approximately linear dispersion. This linear dependence at small |q| leads to an approximately quadratic increase in the low-energy DoS,²⁷ which is the basis for the Debye approximation common in traditional texts on spin–phonon coupling.²² The acoustic modes have frequencies in the Hz–kHz range at q ≃ Γ and allow the transport of sound waves through the crystal, which is the origin of their name. The slope of the dispersion ∂ω_qj/∂q determines the “group velocity” v_qj of the modes, which gives the velocity and direction with which the modes travel through the crystal. For NaCl, the linear dispersion of the acoustic modes near Γ results in these modes having the largest |v_qj|, which sets the speed of sound propagation through the crystal.

Fig. 5 Phonon dispersion ħω_j(q) and density of states g(ħω) (DoS) of NaCl from DFT calculations (see ESI,† for details). The dispersion is shown along the wavevector path . There are 3n_a = 6 bands at each q, but the high symmetry of the cubic spacegroup means that a pair of acoustic modes and a pair of optic modes are degenerate along this path; the singly-degenerate bands are shown in red, and the doubly-degenerate bands in blue. The DoS is shown in black and the projections onto the Na and Cl atoms are shown in yellow and green, respectively.

The remaining 3n_a − 3 modes involve opposing motions of the atoms in the unit cell, and thus have ħω_j > 0 at all q. These modes can, in principle, lead to a change in electric polarisation (the solid-state analogue of the dipole moment) and/or electric polarisability within the primitive cell, and thus may interact with light; they are therefore termed “optic modes”. Approximately two-thirds along the path between Γ and X in the phonon dispersion of NaCl (Fig. 5) there is an avoided crossing between an acoustic and optic branch, indicating mixing between these modes; at wavevectors away from Γ, the distinction between these two classes of modes is therefore not always clear. The DoS of NaCl shows a continuous spectrum of modes up to ∼250 cm⁻¹, and the partial DoS shows that Na and Cl contribute roughly equally across the full range of modes. This is due to their similar mass (m_Na = 22.99 and m_Cl = 35.45 amu), and in systems with large mass differences between the constituent atoms the DoS often shows distinct energy regions where a particular type or group of atoms dominates the vibrations.

The phonon spectra of molecular crystals present some additional noteworthy features, which we illustrate here with the phonon dispersion and DoS of the cubic phase of crystalline NH₃ (Fig. 6; calculated using DFT, see ESI†). In the general case where there are multiple molecules in the primitive cell, the translations and rotations of each molecule (i.e. those that are present in the gas phase) combine to produce low-energy phonons involving various rigid-body motions of whole molecules that are often termed “external modes”. Three combinations of these translations preserve the intermolecular distances and give rise to the acoustic modes with ħω_j = 0 at q = Γ as in inorganic crystals. Other combinations of translations, rotations, and translations/rotations have non-zero energy, but are still relatively low-energy motions due to the generally weaker intermolecular interactions, and hence there are additionally a large number of low-energy dispersive modes that mix with the acoustic modes at q ≠ Γ; we refer to this loosely-defined set of modes as “pseudo-acoustic” modes.

Fig. 6 Phonon dispersion ħωj(q) and density of states g(ħω) (DoS) of cubic crystalline NH₃ from DFT calculations (see ESI,† for details). Panel (a) shows the full DoS and highlights four distinct groups of modes, viz. the “external” (pseudo-)acoustic modes formed from rigid-molecule translations and rotations (blue), and the “internal” modes formed from the NH₃ bending (red), H–N–H “scissoring” (orange), and N–H symmetric and antisymmetric stretch vibrations (purple). Panel (b) shows the phonon dispersion and DoS from 0–700 cm⁻¹ corresponding to the blue region in (a).

As a non-linear molecule, NH₃ has three translations, three rotations and six vibrations in the gas phase. There are four NH₃ molecules in the primitive cell of the cubic crystal structure (n_a = 4 × 4 = 16) and hence 48 phonon modes at each wavevector q. The rigid-molecule translations and rotations combine to form a dense group of 24 phonon bands spanning a range of 0–700 cm⁻¹ (Fig. 6). In this case the translations and rotations are separated by a small gap of ∼100 cm⁻¹ between the 12 translational (pseudo-)acoustic modes from around 0–200 cm⁻¹, and the 12 rotational pseudo-acoustic modes from ∼300–700 cm⁻¹. Both sets of phonons are pure external modes, and the wide dispersion is due to the H-bonding network formed between the molecules in the unit cell. The low-energy DoS remains approximately quadratic as in the Debye approximation, as the region where ħω_j → 0 is dominated by the acoustic modes, but the DoS ceases to behave quadratically above the “Debye limit” due to the pseudo-acoustic modes.²⁸

The mid- and high-energy modes in molecular solids tend to be made up of combinations of intramolecular vibrations (“internal modes”). Particularly at higher energies these tend to be localised vibrations (e.g. bond stretches) that often do not significantly disrupt the intermolecular interactions, resulting in a weak dependence on q and hence flat (“dispersionless”) bands and sharper features in the DoS. For cubic NH₃, there are three distinct groups of internal modes (Fig. 6a) corresponding to combinations of the bending modes (1100–1150 cm⁻¹, four modes), scissoring modes (1600–1700 cm⁻¹, eight modes), and the symmetric and antisymmetric stretches (3250–3450 cm⁻¹, 12 modes). The two types of N–H stretches do not mix, and the higher-energy group of stretching modes is notably split into two distinct features centred around 3290 and 3400 cm⁻¹, with four and eight modes, respectively. The features corresponding to the internal modes are much narrower than those from the external modes, indicating a narrower energy dispersion.

To conclude this section, we cover three more aspects of phonons that are relevant to the material in the following sections. Firstly, while the phonon energies are quantised into units of ħω_qj as in the gas phase, for most purposes it is the average harmonic energy level of each mode that is of interest, which is determined by the Bose–Einstein occupation number n̄_qj:

Secondly, the analogous expression for the atomic displacements as a function of the mode amplitude in eqn (3) for phonons in solids is:

where the factor of is required to maintain the relation in eqn (2). This can be used to compute the coupling between phonon modes and other material properties, in particular the CF parameters (CFPs), as discussed in the following sections.

Finally, in the harmonic approximation the phonon modes are strictly independent oscillators and therefore have infinite lifetimes τ_qj. However in real crystals, the phonon modes interact with one-another via collisions and decay, which give rise to finite lifetimes. As a consequence of the uncertainty principle, the phonons therefore have finite linewidths Γ_qj = ħ/τ_qj (full-width-at-half-maximum in energy units). [n.b. The symbol used for the phonon linewidths Γ_qj should not be confused with the centre of the BZ q = Γ.] In the case of lifetime broadening such as this, the phonon lineshape is Lorentzian (eqn (10)). However, as the phonon DoS (eqn (7)) must go to zero at zero energy (i.e. g(0) = 0), the anti-Lorentzian lineshape should be used instead (eqn (11)). In practice, however, there is little difference if Lorentzian or even Gaussian lineshapes are used.

Phonon–phonon scattering is typically the main contributor to the lifetimes τ_qj in semiconductors and insulators, although in some systems coupling between the electrons and phonons can also be important.²⁹ The lowest-order scattering processes that conserve energy and (crystal) momentum are three-phonon processes; these can be calculated and used in conjunction with the harmonic phonon frequencies and eigenvectors to compute the τ_qj from first principles,³⁰ albeit at a considerable computational cost, especially for large systems such as molecular crystals. Hence, phonon linewidths are often treated as parameters³¹ or by using statistical approximations.³² However, we have recently shown that there is little impact on the spin dynamics when either fixed, approximated, or calculated linewidths are used, provided the BZ is integrated with a sufficiently dense q-point grid.³³

4 Electronic structure of lanthanide complexes

Now we turn to the spin states of the molecule, where herein we focus on the low-lying electronic states of Ln complexes. For most cases, the trivalent Ln³⁺ oxidation state is most relevant, with a well-isolated ground configuration 4fⁿ5s²5p⁶; the 6s and 5d may become relevant depending on the oxidation state and the ligands surrounding the metal.³⁴ Trivalent Ln³⁺ ions are unique in the periodic table in that their valence electrons reside in the 4f orbitals which are strongly contracted and shielded beneath a filled set of 5s and 5p orbitals. Given the inert nature of the 4f orbitals, the electronic structures of Ln³⁺ ions in molecules are well-described by the free-ion Hunds rule ground term ^2S+1L_J, where S is the total spin angular momentum, L is the total orbital angular momentum, and J = L ± S is the total angular momentum after spin–orbit coupling (the ± sign refers to either >7 or <7 electrons in the 4f shell). While there are exceptions (e.g. Sm³⁺ and Eu³⁺), for most Ln³⁺ ions the population of excited J multiplets is essentially zero at room temperature and below, which is the temperature range of interest for SMMs. Coordination of a set of ligands to a Ln³⁺ ion gives rise to a CF splitting of the 2J + 1 m_J states of the ground J multiplet, and is described by the CF Hamiltonian (eqn (12)).³⁵

Here, B^q_k are the Stevens CF parameters of rank k and order q, the Ô^q_k are Stevens operator equivalents in the |J,m_J〉 basis which are polynomials of the angular momentum operators Ĵ_x, Ĵ_y, Ĵ_z, and the θ_k are operator equivalent factors which relate the multi-electron CF Hamiltonian in the |J,m_J〉 basis to the single-electron CF Hamiltonian in the Slater determinant basis.^36,37

5 Spin–phonon coupling

Under an initial approximation that the electronic spin system is non-interacting with the phonons, we can define the uncoupled equilibrium Hamiltonian (eqn (13)). Here, the electronic part is the CF Hamiltonian and the phonon part is the harmonic oscillator Hamiltonian; note the distinction between the phonon wavevector q and the order of the CF operator q.

The operators â_qj and â^†_qj denote phonon annihilation and creation operators, which act as and and represent reducing or increasing the vibrational quantum number of phonon mode qj. These allow us to define the dimensionless mass- and frequency-weighted normal coordinates as²⁷

We note that the Q̂_qj operator in eqn (14) is non-Hermitian, unless q coincides with the Γ point or at special points on the boundary of the BZ (see ESI†). This occurs because normal mode displacements are in general complex quantities characterised by amplitude and phase. The Hamiltonian in eqn (13) can be written in the direct (or tensor) product basis and is solved by diagonalisation to give eigenstates where ψ denotes an electronic eigenstate and are non-negative integers specifying the occupation number of each phonon mode.

The fundamental nature of the spin–phonon interaction considers how vibrations modulate electronic states. To this end, we describe spin–phonon coupling for the case of Ln SMMs using a modified form of the CF Hamiltonian, where the parameters B^q_k are dependent on the displacement along the vibrational mode coordinates Q_qj (eqn (15)). In the case of a transition metal, one may alternatively define the coupling Hamiltonian to arise from vibrational modulation of hyperfine coupling, the g-matrix, or zero-field splitting, for example;^32,38 however for systems, it is important to account for changes in the nature of the basis states that is not accommodated in a simple g-matrix model.³⁹ Though in any case, we note that a model Hamiltonian need not be assumed and that the derivatives of the Hamiltonian matrix elements can be determined directly and employed in a similar vein.^40,41 In eqn (15) we describe the modulation of the CF parameters using a Taylor series centred at the equilibrium structure Q_qj = 0, where B^q_k(0) are the same equilibrium CF parameters appearing in eqn (13) and thus are not relevant for spin–phonon coupling. Thus, we can define a spin–phonon coupling Hamiltonian Ĥ_SP up to second-order as eqn (16), where V̂⁽¹⁾_qj and describe first- and second-order spin–phonon coupling, respectively.

We note that, as a Hamiltonian operator, Ĥ_SP in eqn (16) must be Hermitian. However, due to the non-Hermitian definition of Q̂_qj in eqn (14), this implies that the derivatives in eqn (16) are taken using complex-valued displacements along the normal mode coordinates.¹⁴ Alternatively, one can choose to re-write eqn (16) such that each term in the sum over qj are all individually Hermitian; we do this by introducing a new set of creation/annihilation operators b̂_qj as linear combinations of â_qj and â_−qj, and defining Hermitian dimensionless normal displacement operators as X̂_qj = b̂_qj + b̂^†_qj (see ESI†). For the rest of this review, we only consider the Hermitian representation of the normal displacements in terms of b̂_qj and b̂^†_qj, denoting their corresponding number states as |v_qj〉.

Eqn (16) suggests that the magnitude of the spin–phonon coupling is dependent on the number of q points included in the sum within the BZ, as the sums over q collect a growing number of terms. This apparent paradox is resolved by considering eqn (9) which includes a factor thus decreasing the nuclear displacement amplitudes, and hence the coupling of individual modes, when the sum over q increases, resulting in size-consistent and convergent behaviour under these definitions.

The spin–phonon coupling Hamiltonian is evaluated in the product basis, and for the first order term V̂⁽¹⁾_qj the non-zero matrix elements are of the form:

where mode qj has either gained or lost a quantum of vibrational energy. Evaluation of the electronic part is straightforward as it is simply a matrix element of CF operators. The vibrational matrix element can be evaluated with reference to the definition of the position operator in the basis of harmonic eigenstates, giving:Considering the second-order coupling with modes qj and q′j′, the non-zero matrix elements are given by:andwhereWe will see in Section 6 that crucial quantities for determining the spin dynamics turn out to be the square modulus of the matrix elements of Ĥ_SP (eqn (17) and (20)), which consist of an electronic and a vibrational component. When considering the vibrational terms |〈v_qj ± 1|X̂_qj|v_qj〉|², we will be concerned with cases where the phonon modes are in thermal equilibrium; i.e. the v_qj harmonic vibrational levels are statistically occupied with a probability proportional to the Boltzmann factor (where Z is the partition function ). Thus, averaging over thermally populated vibrational levels, giveswhere n̄_qj is the Bose–Einstein occupation number defined in eqn (8). While calculating the vibrational and electronic component of the matrix elements of Ĥ_SP is straightforward, obtaining the spin–phonon coupling parameters and is not straightforward; see Section 7. With the total Hamiltonian of a coupled spin and phonon system up to second-order in the spin–phonon coupling defined as Ĥ_T = Ĥ_eq + Ĥ_SP, one may examine how the spin–phonon coupling affects molecular properties. For instance, we could calculate vibronic spectra directly in the coupled spin-vibrational basis,^13,42 or we could calculate how the phonon bath leads to magnetic relaxation of the electronic spin system;^14,31 the latter will be our focus here.

6 Spin dynamics and magnetic relaxation

In order to discuss magnetic relaxation, we must consider how the spin and phonon systems evolve in time and influence one-another; from here we will refer to phonons as “the bath” interchangeably. The goal of our interrogation is to assess the probability that a spin–phonon interaction leads to a spin flip in the electronic system, i.e. magnetic relaxation. In general, the dynamics of the phonon bath is assumed to be much faster than the spin dynamics of the SMM (though this assumption does not have to be made, see below), and the overall phonon–molecule interactions leading to a spin flip include: (i) phonons are absorbed by the SMM, (ii) phonons are emitted by the SMM, and (iii) phonons are scattered inelastically by the SMM. This section is separated into two subsections: in the first, we discuss open quantum systems methods which contain the most accurate approaches for dealing with the dynamics of quantum systems, starting from a minimal set of approximations and relaxing the accuracy of our approach. In the second, we discuss methods based on perturbation theory, which are the traditional approaches taken in this field, relying on a set of approximations from the outset.

6.1 Open quantum systems methods

As we are concerned with calculating the spin-flip probabilities, we must work with probabilistic measures of the quantum states; this necessitates use of the density matrix, [small rho, Greek, circumflex]. The density matrix encodes the probabilities of finding the total system in a particular state (diagonal matrix elements) and also allows us to describe quantum coherences between states (off-diagonal matrix elements); note that here the states of the total system include both spin and phonon degrees of freedom. At the most exact level of theory, we could construct a density matrix representation of the spin-vibrational system described by the total Hamiltonian Ĥ_T starting in a well-defined state, e.g. [small rho, Greek, circumflex](0) = |Ψ_spin–vib〉 〈Ψ_spin–vib|, and consider its unitary evolution according to the Liouville–von Neumann equation [small rho, Greek, circumflex](t) = e^−iĤ_Tt/ħ[small rho, Greek, circumflex](0)e^iĤ_Tt/ħ (the analogue of the time-dependent Schrödinger equation).⁴³ Performing such a calculation would include all interactions between the spin and phonons and give us the knowledge of the system at all times. However, this is prohibitively expensive in practice owing to the very large number of phonon modes in any real system, and hence there have arisen many ways to find approximate solutions to this problem: this is a field of research called open quantum systems.^44–47 The first simplification in solving this problem is to identify that we are not very interested in the state of the bath in the end, rather, we only care about the probability of the SMM undergoing a spin flip. Hence, we can perform a partial trace of [small rho, Greek, circumflex] over the phonon modes: this is a mathematical operation to remove the phonon degrees of freedom from the total density matrix by averaging over them (often referred to as “tracing out the bath”) to give the reduced density matrix of the spin system:Note that the expectation values of [small rho, Greek, circumflex] in the above equation are taken with respect to the phonon degrees of freedom only, leaving us with a reduced density matrix for the spin system. This is made evident in the second line of eqn (26), where the matrix elements of [small rho, Greek, circumflex]_S in the electronic eigenbasis are given explicitly in terms of the matrix element of [small rho, Greek, circumflex]. This procedure could be done either before or after interrogation of spin system dynamics. For instance, chain-mapping methods such as time-evolving density using orthogonal polynomial algorithm (TEDOPA)^48,49 employ a linearisation of the bath Hamiltonian to approximate the dynamics of the whole system with density matrix renormalisation group (DMRG) methods,⁵⁰ after which a trace over the bath degrees of freedom gives the desired reduced density matrix. Alternatively, another class of methods aims at finding approximations for the time evolution of the reduced density matrix of the spin system alone (eqn (26)) after performing the partial trace over the bath degrees of freedom analytically. In this class we find methods such as the quasi-adiabatic path integral (QUAPI),^51,52 hierarchy equations of motion (HEOM),^53,54 or the time-evolving matrix product operator (TEMPO).⁵⁵ While all of the above are examples of methods that can approach high-degrees of accuracy that may be required when the spin and phonon systems are strongly coupled, in the domain of molecular magnetism it is usually safe to assume that the spin and phonon systems are weakly coupled. Moreover, magnetic relaxation is typically much slower than the underlying phonon dynamics, thus allowing for a clear separation of system and bath timescales; we have recently confirmed this separation of timescales by calculating the phonon lifetimes for a molecular crystal of a Dy(III) SMM.³³ These additional constraints simplify the problem considerably, allowing one to treat the phonons as a weakly coupled bath with no memory (on the timescale of the spin system). Under these assumptions, commonly referred to as the Born–Markov approximation, the dynamics of the reduced density matrix is described by the Redfield equation:^44,56where the square brackets denote the commutator ([Â,B̂] = ÂB̂−B̂Â), and and are the Redfield superoperators arising from the first- and second-order spin–phonon coupling in eqn (16). Denoting the electronic part of the first- and second-order spin–phonon coupling operators as:the Redfield superoperators can be written as:where:The time-dependent coefficients in eqn (32) and (33) are defined aswhere is the total vibrational Hamiltonian in eqn (13). Eqn (34) and (35) represent the bath correlation functions calculated for a thermal distribution of phonons (eqn (8)). The thermal averages over phonons appearing in eqn (34) and (35) can be calculated as 〈…〉_eq = Tr_vib{…[small rho, Greek, circumflex]_eq}, where [small rho, Greek, circumflex]_eq = e^{−Ĥ_vib/k_BT}/Tr_vib{e^{−Ĥ_vib/k_BT}} is the thermal density matrix for phonons. The CF Hamiltonian appearing in eqn (27) includes a correction arising from the second-order spin–phonon coupling, accounting for the average energy shift induced by the quadratic terms.

The Redfield superoperators and in eqn (27) describe one- and two-phonon processes corresponding to the linear and quadratic terms of the expansion in eqn (16) (i.e. Orbach and Raman-II processes; see Section 6.2). In fact, Redfield theory incorporates the effect of the spin–phonon coupling in eqn (16) only to lowest order; higher-order interactions, i.e. two-phonon processes induced by the linear spin–phonon coupling (Raman-I mechanism) require either extending Redfield theory or resorting to other strategies that neglect the influence of electronic coherences (see Section 6.2 and ESI†).¹⁴ In principle, the effect of arbitrarily high-order spin–phonon interactions on the dynamics of the reduced spin density matrix [small rho, Greek, circumflex]_S can be included systematically in a quantum master equation of a similar form to eqn (27)via the time-convolutionless projection operator method,⁴⁴ which amounts to expanding the reduced system dynamics in powers of the system-bath coupling, giving further corrections to the equation of motion for [small rho, Greek, circumflex]_S (eqn (27)). In practice, however, going beyond second-order is not very common for systems with more than a few electronic states or with highly structured vibrational environments, since the calculation of higher-order contributions to the dynamics will likely require computational resources comparable to the ones required for numerically exact methods (such as the ones mentioned above), which are more generally applicable, since they do not rely on a perturbative expansion in the system-bath coupling. However, the Redfield equation (and its secular version) remain extremely valuable tools to investigate the dissipative dynamics of quantum systems, especially when used in conjunction with other prescriptions to extend its validity beyond the Born–Markov approximation. A common strategy of this type is to redefine the boundary between system and bath by partitioning the total Hamiltonian Ĥ_T = Ĥ_eq + Ĥ_SP in a different way, in general such as where the revised system Hamiltonian now contains the “important” spin-bath dynamics and the residual system-bath coupling is ammenable to description with the Born–Markov approximation. As an example, the reaction coordinate mapping method⁵⁷ aims to capture non-Markovian effects on the electronic system dynamics by identifying a single collective coordinate of the bath, given by a suitable linear combination of individual coordinates X̂_qj, which is treated exactly, while the rest of the bath, which only couples to the collective mode, is treated with Redfield theory. Different criteria for selecting the most relevant bath pseudo-modes may be chosen, leading to different methods.^58,59 Another strategy is to describe the coupled spin–phonon system in terms of polarons, displacing the phonon coordinates conditionally on the state of the spin system.⁴⁵ The residual spin–phonon coupling in the polaron frame describes small harmonic displacements around the redefined equilibrium configurations, which can then be treated within Redfield formalism.^60,61

We have discussed how the Redfield equation (eqn (27)) describes the dynamics of the reduced spin density matrix under the Born–Markov approximation, capturing the coupled time evolution of both electronic populations 〈ψ_i|[small rho, Greek, circumflex]_S|ψ_i〉 and coherences 〈ψ_i|[small rho, Greek, circumflex]_S|ψ_k〉 (where i ≠ k). A further approximation is often performed at this point, which amounts to averaging the oscillations induced by the time evolution of the electronic spin–phonon coupling operators in eqn (32) and (33). This procedure, known as the secular approximation,⁴⁴ decouples the evolution of populations from coherences, and is justified when magnetic relaxation is much slower than the typical timescale of the unitary spin system dynamics, determined by energy differences between different eigenvalues of the CF Hamiltonian. Since SMMs exhibit either exactly or nearly doubly-degenerate electronic states (Fig. 1), this approximation is not always applicable, as the interaction between populations and coherences within a degenerate doublet cannot be decoupled, and other strategies must be used instead.¹⁴ If we decide to neglect electronic coherences entirely (which would not be appropriate if we cared about the decoherence of the quantum spin system induced by the phonons (i.e. the T₂ timescale),⁶² but is usually a good approximation if we are only concerned with magnetic relaxation (i.e. the T₁ timescale)), we can calculate the rates for a transition between two eigenstates of the CF Hamiltonian, ψ_i → ψ_f (given symbol γ_fi), by taking the expectation value of the Redfield superoperators (eqn (30) and (31)) on the final state ψ_f, and singling out the contributions from the terms 〈ψ_i|[small rho, Greek, circumflex]_S|ψ_i〉. For instance, applying this procedure to gives:

which describes single-phonon emission (arising from the first term in the square braces) or absorption (second term in the square braces) by the spin system. In the above, the phonon occupation numbers arise from expanding the mode coordinates X̂_qj in terms of the phonon creation and annihilation operators (as performed in eqn (14) and considering the thermal averages (eqn (24) and (25)) where n̄_qj is defined in eqn (8):

〈b̂_qjb̂^†_qj〉_eq = n̄_qj + 1 (37)

〈b̂^†_qjb̂_qj〉_eq = n̄_qj. (38)

In eqn (36), the phonon modes are considered in the harmonic limit as having infinitely-long lifetimes and hence the lineshape is given by the Dirac delta function; we could augment this expression by considering a finite linewidth, for example described by eqn (11).

Eqn (36) should be familiar to readers who have read traditional texts on magnetic relaxation or spin–phonon coupling: this is almost identical to the expression for single-phonon absorption or emission obtained by applying first-order time-dependent perturbation theory, viz. Fermi's Golden Rule, to the problem of magnetic relaxation; this is the route we take in the next section.

6.2 Perturbation theory methods

When we want to calculate the rates of magnetic relaxation in SMMs, we are not concerned with electronic coherences; if, on the other hand one wants to calculate decoherence rates, read Section 6.1. In this case, we can adopt a simpler theoretical framework and describe transitions between different states of the spin system using a classical rate equation. Furthermore, unlike for the single-phonon rate equation obtained from the first-order Redfield superoperator (eqn (36)), it is tedious to obtain rate-like equations for electronic populations considering the higher-order superoperators from the time-convolutionless expansion, and using perturbation theory (viz. Fermi's golden rule) is far more straightforward for this task.¹⁴ These two approaches lead to the same rate expressions, provided that some reasonable assumptions on the spin–phonon coupling strength and on the timescales of spin and phonon dynamics are made (see ESI† for a detailed derivation of one- and two-phonon rates). Using first-order time-dependent perturbation theory, the transition rate γ_fi between any two pairs of states |ψ_i〉 and |ψ_f〉 caused by a periodic perturbation V̂ can be calculated as:⁴³where ρ(E_f) is the density of final states such that ρ(E_f)dE_f is the number of final states in the energy range E_f → E_f + dE_f. We first consider cases where the absorption or emission of a single phonon by the spin system drives transitions between electronic eigenstates |ψ_i〉 and |ψ_f〉 of the equilibrium electronic Hamiltonian (eqn (12)). For the case of magnetic relaxation, we only care about the states of the electronic spin system and not the bath states, and so we can perform the analogous procedure to “tracing out the bath” by first separating 〈ψ_f|V̂|ψ_i〉 into an electronic and vibrational part (viz.eqn (17) and (20)), then summing the contribution of all modes qj independently and assuming a thermal state of their phonon occupations (eqn (8)), which gives eqn (40) and (41), where the coupling operators are defined in eqn (28), E_i and E_f are the equilibrium electronic eigenstate energies. The only difference between this expression and eqn (36) is that we allow the phonon modes to have a finite linewidth, described by the lineshape ρ_qj (eqn (11)); note that the integration variable ħω is an argument of the phonon occupation terms, the delta function, and the lineshape), and have explicitly separated the phonon absorption (γ_fi⁻) and emission (γ_fi⁺) terms.Since phonons have a finite linewidth, we allow their occupation numbers to vary continuously with the energy, i.e. n̄(ħω) = 1/[exp(ħω/k_BT) − 1]. Note that for these single-phonon transitions, the integrals can be performed analytically to give only non-zero contributions where ħω = ±(E_f − E_i), and we remind the reader that normalisation for the number of q-points included in the summation is accounted for in the definition of the spin–phonon coupling coefficients in V̂^(1e)_qjvia the displacements in eqn (9). The approximation here is that the electronic spin–phonon coupling matrix element 〈ψ_f|V̂^(1e)_qj|ψ_i〉 does not change over the phonon lineshape: this is a fair approximation for the case of narrow lines, but n.b. the zero-point displacement of a phonon with a different energy will be different from the central value, and this is not accounted for in this approximation.

Considering now two-phonon interactions, there are three possibilities: (i) two phonons are absorbed by the spin system (γ_fi⁻⁻); (ii) two phonons are emitted by the spin system (γ_fi⁺⁺); or (iii) a phonon is scattered inelastically by the spin system (γ_fi⁻⁺ and γ_fi⁺⁻). Using first-order time-dependant perturbation theory as before, a similar derivation to the first-order rates can be carried out using the second-order spin–phonon coupling perturbation : this gives transition rates corresponding to the Raman-II mechanism (eqn (42)–(45)).^14,21 Alternatively, one can use the first-order spin–phonon coupling Hamiltonian and take the time-dependent perturbation theory to second-order, thus defining the Raman-I mechanism (eqn (46)–(49)).^14,20,23

Note that the terms appearing in the double sums over qj,qj′ in eqn (42)–(49) are symmetric with respect to exchange of the mode indices qj ↔ q′j′. The terms involving simultaneous phonon emission and absorption (labelled by superscripts +− and −+) simply transform into each other when swapping mode indices. Therefore, the double summation can be replaced by twice the restricted sum over pairs of modes (qj ≥ q′j′) and a double-counting correction factor for the case in which the mode indices are the same inserted between the summation and double integral, for computational efficiency; we note that all terms should have the factor of and not some with as given in ref. 14.

With these transition rates calculated (n.b. total transition rates γ_fi are the sum of all contributions, eqn (40)–(49)), simulation of the spin dynamics is then a case of solving a classical rate equation for the electronic eigenstate populations:⁶³

For this method to be valid, the γ_fi must obey the principle of detailed balance such that in equilibrium the forwards and backwards probability currents γ_fip_i and γ_ifp_f are equal due to microscopic reversibility:

Solving the coupled differential eqn (50) can be achieved by constructing the matrix Γ with elements (note that Γ here is the rate matrix and not the origin of the BZ, nor is it the phonon linewidth) and diagonalising‡ it (i.e. finding the matrix ϕ, that transforms Γ into a diagonal form Λ = ϕ⁻¹Γϕ;).⁶³ The eigenvalues of Γ (diagonal entries of Λ) are a set of 2J + 1 characteristic rates −τ_k⁻¹, each corresponding to a so-called “normal mode” of relaxation (n.b. these are not the normal modes of vibration, nor are they different relaxation mechanisms). Of these normal modes, one rate corresponds to the system in equilibrium and is identically zero. When considering single-phonon processes alone for SMMs, often the Orbach mechanism dominates, and one of the normal modes of relaxation corresponds to the rate of this over-barrier process, which is many orders of magnitude slower than the remaining modes consisting of fast fluctuations of the populations of the states on either side of the barrier.³ Hence, often this rate is taken to be the same as the magnetic relaxation rate. In principle the spin dynamics of the system evolve under the influence of all normal modes of relaxation, so the magnetic relaxation rate should be determined by first simulating a magnetisation decay process and then fitting the trace with an exponential decay to give τ:

To do so, one must first calculate the population of each state |ψ_i〉 as a function of time t:with λ = ϕ⁻¹p₀, where p₀ is a vector specifying the initial populations of the states; to simulate magnetic relaxation, one of the degenerate ground states would have unit population and all other states zero population. In this approach, the initial population evolves under the influence of each normal mode of relaxation (eigenvector of Γ). Then, the time-dependence of the total magnetisation is given as the sum over the magnetisation of each state, weighted by its population:

7 Practical considerations for ab initio calculation of spin–phonon coupling

Following the above theoretical descriptions, a practical ab initio calculation of spin–phonon coupling and spin dynamics requires four steps: (i) optimisation of the structure and calculation of the phonon modes; (ii) calculation of the equilibrium electronic structure; (iii) calculation of the spin–phonon coupling terms; and (iv) simulation of the spin dynamics.

The first step is usually performed using (semi-)local DFT methods, which are currently the only tractable methods for optimisation and lattice-dynamics calculations on periodic molecular crystals. Periodic DFT calculations commonly describe the valence electronic wavefunctions with a plane-wave basis set and use pseudopotentials to describe the nuclei and core electrons (e.g. in VASP,^64,65 CASTEP⁶⁶ or CP2K⁶⁷), but all-electron calculations and atom-centred Gaussian-type atomic orbital basis sets are also an option (e.g. in CRYSTAL⁶⁸). Optimisation of molecular crystal structures should generally be performed using a method including dispersion corrections (e.g. DFT-D⁶⁹) in order to obtain qualitatively correct intermolecular potential energies.

Even when starting the optimisation of a molecular crystal from an experimental structure obtained with X-ray diffraction, the optimisation process is not necessarily straightforward owing to the presence of crystallographic disorder and shallow potential energy surfaces that can result in the optimisation algorithm getting stuck in a local maximum. The hallmark of a geometry optimisation landing in a local maximum on the PES is the presence of imaginary phonon modes (ω² < 0) at the BZ centre (Γ point). These “dynamical instabilities” imply that a distortion of the atoms in the unit cell along the imaginary modes will lower the total energy, and hence that the unit cell contents are not fully optimised. Just like for molecular structure optimisations in the gas-phase, the imaginary modes can be removed by distorting the geometry along the eigenvector(s) corresponding to the imaginary phonon mode(s) and re-starting the optimisation at the minimum energy point along the distortion path. Doing so may result in a lowering of the crystallographic spacegroup symmetry. Typically one would distort along the largest magnitude imaginary mode, although this does not necessarily lead to the deepest energy minimum. After optimisation, the Γ-point phonon spectrum should be recalculated and checked for imaginary modes, and this process repeated until they are all removed.

When all phonon modes have ω² ≥ 0 at the Γ point, then the unit cell can be considered optimised. Next, the phonon modes and their dispersion in reciprocal space should be calculated. For this, one chooses a path in reciprocal space visiting the unique high-symmetry q-points in the BZ, which are determined by the crystallographic spacegroup. If all phonon branches have non-imaginary energies at all the high-symmetry wavevectors in the BZ, then the phonon calculation can be considered adequate. If there are imaginary modes at non-zero wavevectors, then these should be investigated. Off-Γ imaginary modes can arise in two cases: (i) interpolation artefacts for wavevectors that are not commensurate with the supercell expansion used in a finite-displacement calculation (or are not included in the “coarse” q-point mesh in a perturbation-theory calculation); or (ii) “true” imaginary modes that indicate a dynamical instability for which an expanded unit-cell is required to contain the energy-lowering distortion. To check which case is causing the imaginary mode, one must perform the dispersion calculation with an expanded supercell (or a larger “coarse” mesh) that includes the wavevector(s) with imaginary mode(s). If the mode(s) becomes real, then it is case (i) and the phonons must simply be calculated in the larger supercell. If the mode remains imaginary, then it is a case of (ii) and the unit cell must be enlarged to contain the wavevector in question and the structure re-optimised. Checks for imaginary modes must then be repeated on the enlarged unit cell, and further calculations performed as necessary to remove them. In general, using a larger supercell expansion when calculating the phonon modes, which is equivalent to finer q-point sampling, will improve the accuracy with which the frequencies of off-Γ modes are evaluated, but this is limited in practice by computational cost.

In the case of low-symmetry molecular crystals, it may not always be computationally tractable to remove all imaginary modes. As an example, we recently calculated the phonon spectra of [Dy(Cp^ttt)₂][B(C₆F₅)₄] (Cp^ttt = 1,2,4-^tBu₃-C₅H₂; ^tBu = C(CH₃)₃) at various pressures, and found that a 2 × 2 × 1 supercell with 1104 atoms was required to remove imaginary modes at most pressure points. Despite this, imaginary modes with ħω ∼ 10i cm⁻¹ persisted at two pressures, and phonon calculations with the 2 × 2 × 2 supercell (2208 atoms) required to remove them were simply not computationally feasible.⁷⁰ A detailed review of the physical significance of imaginary modes, possible causes, and treatments can be found in ref. 71.

This workflow for determining phonon spectra applies equally to all molecular crystals, however there are extra considerations for metal ions. For 4f elements, the near-degenerate, core-like 4f electrons, which present a problem for the single-determinant ansatz of DFT (see below) can be resolved by subsuming the 4f electrons into the pseudopotential core using “f-in-core” potentials – this is an acceptable approximation as the potential energy surface is largely unaffected by the 4f electrons. For d-block or 5f elements, the substantial bonding character of the d and 5f electrons makes this approach inappropriate, and hence these calculations are more challenging. Here, the self-interaction error in (semi-)local DFT functionals tends to overly delocalise the d- or f-states, which can be empirically corrected with an additional on-site term known as the Hubbard U parameter using the DFT+U method.^72,73 However, even this approach does not guarantee success, and often simply converging the electronic wavefunction is a challenge and can make it difficult both to optimise the structure and, subsequently, to calculate the accurate forces required for phonon calculations. In these cases, sometimes substitution of the open-shell metal ion for an isovalent, closed-shell analogue with similar chemistry (e.g. replacing Co(II) with Zn(II) in octahedral environments to remove complications from the ground orbital triplet) is a suitable approximation.

The second step, calculation of the equilibrium electronic structure, requires explicit treatment of ground and excited electronic states, as well as inclusion of spin–orbit coupling, in order to accurately represent the magnetic spin states in molecules. Compared to closed-shell molecules, for which DFT is by far the dominant approach, such “single-determinant” electronic structure methods applied to open-shell molecules fail to correctly describe the ground state, where a single electron configuration fails as a qualitatively correct description of the electronic wavefunction.⁷⁴ Instead, a multi-determinantal method is required to include the ground and excited electron configurations relevant for the magnetic properties. Often in the case of monometallic metal complexes, the magnetic states are well-localised and the d- or f-orbitals are the most important (this is obviously not true in the presence of strong magnetic coupling between multiple spin centres, e.g. for [Cp^iPr₅DyI₃DyCp^iPr₅]⁶). Here, the complete active space self-consistent field (CASSCF)⁷⁵ method provides a suitable approximation, where all electron configurations are included within a given “active” orbital space; this method is available in many codes (e.g. in OpenMolcas,^41,76 Orca⁷⁷ and Q-Chem⁷⁸). In the case of Ln³⁺ complexes, the active space typically includes the seven 4f orbitals occupied by n electrons. The number of states and different spin multiplicities included in the calculation is user-determined, and for magnetic properties it is advised to be guided by the low-lying Russell–Saunders atomic terms ^2S+1L.^22,35 Spin–orbit coupling can then be added in a second step by mixing the CASSCF states of different spin multiplicities,⁷⁹ leading to states representing the ^2S+1L_J multiplets split into 2J + 1 components by the CF potential. Owing to the relatively small CF splitting of the 4f orbitals, the low-lying electronic states of Ln complexes can be well-described by a model crystal field Hamiltonian (eqn (13)), where the CFPs can be projected directly from the ab initio calculation.⁸⁰ In general, any model Hamiltonian can be projected from such calculations for d- or 5f-containing molecules.⁴¹

While the CASSCF method provides an effective treatment of the so-called “static” electron correlation (arising from the multi-determinantal character of the electronic states), it does not do a good job of describing the “dynamic” electron correlation missing from the mean-field description of electron–electron repulsion (DFT has been so successful because it can approximate this second type of electron correlation very efficiently). This can be overcome by corrections on top of a CASSCF reference wave function either by using variational methods such as multi-reference configuration interaction (MRCI),⁸¹ or through the application of many-body perturbation theory, including methods such as second-order complete active space perturbation theory (CASPT2)^82,83 and second-order n-electron valence state perturbation theory (NEVPT2),⁸⁴ which often provide increased accuracy for magnetic properties.⁸⁰

Correlated wave function methods such as CASSCF usually do not allow the explicit inclusion of extended environments (i.e. beyond the molecule of interest) due to their undesirable scaling behaviour with system size. As an alternative, environments can be included through implicit continuum models or at an atomistic level using hybrid approaches.⁸⁵ We have recently demonstrated that an electrostatic potential derived from atomic point charges of environmental molecules⁸⁶ can provide an appropriate description of the environment. In the case of crystalline environments, the electrostatic potential imposed by finite-size cluster models converges slowly towards the true Madelung potential of an infinite crystal, and has a non-trivial dependence on the shape of the point-charge cluster model.⁸⁷ We have recently shown that the Madelung potential can be closely-approximated in a consistent manner by employing a spherical cluster of unit cells embedded in a conductor-like reaction field to screen the non-physical surface charges that arise from finite-size unit-cell expansions.³³

The third step, computation of the spin–phonon coupling parameters, requires knowledge of the electronic response to nuclear distortions which is encoded in derivatives of the CFPs as presented in Section 5. The most obvious method to obtain these parameters is by computing the electronic structure and CFPs at distorted geometries along the normal mode coordinates and fitting the parameters to a polynomial (Fig. 7); such methods have been used by numerous authors.^{16,31,38,88,89} Clearly, however, for large and/or low symmetry molecules, many expensive ab initio calculations must be performed, and this is not yet even considering the second-order terms; to our knowledge only Sourav and Lunghi have attempted the latter, where they have used machine-learning methods to greatly simplify the calculation.^62,90

Fig. 7 Change in each B^q_k from the equilibrium values for [Yb(trensal)] as a function of mass- and frequency-weighted dimensionless displacement Q_qj along a vibrational mode.⁴² Note that this mode is a totally-symmetric vibrational mode so that the C₃ symmetry of [Yb(trensal)] is not broken and hence only certain B^q_k are affected; hence, some are overlapping with zero magnitude.

Recently, we have shown that the spin–phonon coupling constants can alternatively be obtained analytically using a linear vibronic coupling (LVC) approximation^91–94 based on a single ab initio calculation at the equilibrium geometry.⁴⁰ The LVC method parametrises a first-order expansion of the CASSCF Hamiltonian matrix elements V in the atomic degrees of freedom r^α_κl (for both the SMM and the environment):^95,96

V = W⁽⁰⁾ + W⁽¹⁾ + … (55)

W⁽⁰⁾ = diag(ε₀,ε₁,ε₂,…) (56)

where W^(k) collects terms of order k in the atomic coordinates, which are parameterised in case of k = 1 by electronic energies {ε_n}, gradients {g^α_κl(n)} and non-adiabatic couplings (NACs) {λ^α_κl(n,m)} with respect to the CASSCF states n, m, which have recently become accessible for larger systems through the implementation of density fitting gradients in OpenMolcas.^92,97,98 Although terms of k > 1 can be included in a straight-forward manner to give a higher-order extension of the method, the higher-order derivative couplings become increasingly expensive to compute, though we note that a scheme to parametrise a quadratic model from first-order derivatives has been reported.⁹⁹ Furthermore, in cases where the inclusion of dynamical electron correlation is indispensable, the LVC model can be parameterised based on CASPT2 calculations.^100,101

The final step is calculation of the spin dynamics and magnetic relaxation rates themselves. Here, one must either make their own implementation of the techniques described in Section 6 or rely on tools developed by others. We have developed a suite of Python tools (molcas_suite, angmom_suite and spin–phonon_suite, all available on the PyPI repository§) to facilitate steps (i)–(iii) in conjunction with the VASP, phonopy and OpenMolcas codes, and the program Tau (available on GitLab¶) to calculate magnetic relaxation rates; we have shown that this methodology gives quantitative accuracy with respect to experimental data.³³ Lunghi and co-workers have developed their own methods in the MolForge package, available on GitHub.||

8 Conclusions and outlook

Herein we have discussed the theory and practicalities of the calculation of phonon spectra, spin–phonon coupling and spin dynamics for molecular crystals, with a particular focus on magnetic relaxation in lanthanide single-molecule magnets. This area of research is at the cutting-edge of developments in electronic structure theory, quantum dynamics, machine learning and synthetic chemistry. The depth of knowledge required to take such calculations from start to finish is immense, and many works often make numerous assumptions on the background knowledge of the reader; we hope this Tutorial Review has lifted the veil on at least one aspect for the interested reader. We expect that as packages for these calculations become more advanced and automated, that exploration of chemical space with computational methods could start to make significant inroads into molecular design, both for improving timescales of classical memory storage in single-molecule magnets and coherence times for molecular qubits. Further, we expect that such methods will become crucial for understanding spin–phonon coupling and spin dynamics at the microscopic level for a broader class of nanomaterials in quantum science, such as spintronics and solid-state defect qubits, as well as in biological contexts such as enzyme catalysis and energy transfer.

Author contributions

NFC proposed and supervised the work. All authors equally contributed to writing and editing the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

NFC thanks the ERC for a Starting Grant (STG-851504) and The Royal Society for a University Research Fellowship (URF191320). JMS thanks UKRI for a Future Leaders Fellowship (MR/T043121/1). JGCK thanks The University of Manchester for a Presidential Doctoral Scholarship.

Notes and references

1. A. Caneschi, D. Gatteschi and R. Sessoli, J. Am. Chem. Soc., 1991, 113, 5873–5874.
2. R. Sessoli, D. Gatteschi, A. Caneschi and M. A. Novak, Nature, 1993, 365, 141.
3. D. Gatteschi, R. Sessoli and J. Villain, Molecular Nanomagnets, Oxford University Press, Oxford, 2007, pp. 1–408.
4. N. F. Chilton, Ann. Rev. Mater. Res., 2022, 52, 79–101.
5. F.-S. Guo, B. M. Day, Y.-C. Chen, M.-L. Tong, A. Mansikkamäki and R. A. Layfield, Science, 2018, 362, 1400–1403.
6. C. A. Gould, K. R. McClain, D. Reta, J. G. C. Kragskow, D. A. Marchiori, E. Lachman, E.-S. Choi, J. G. Analytis, R. D. Britt, N. F. Chilton, B. G. Harvey and J. R. Long, Science, 2022, 375, 198–202.
7. Y.-C. Chen and M.-L. Tong, Chem. Sci., 2022, 13, 8716–8726.
8. R. A. Layfield and M. Murugesu, Lanthanides and Actinides in Molecular Magnetism, John Wiley & Sons, 2015.
9. J. L. Liu, Y. C. Chen and M. L. Tong, Chem. Soc. Rev., 2018, 47, 2431–2453.
10. K. L. Harriman, D. Errulat and M. Murugesu, Trends Chem., 2019, 1, 425–439.
11. C. A. P. Goodwin, F. Ortu and D. Reta, Int. J. Quantum Chem., 2020, 120, e26248.
12. M. Atanasov, D. Aravena, E. Suturina, E. Bill, D. Maganas and F. Neese, Coord. Chem. Rev., 2015, 289–290, 177–214.
13. M. Atanasov and F. Neese, J. Phys.: Conf. Ser., 2018, 1148, 012006.
14. A. Lunghi, Sci. Adv., 2022, 8, eabn7880.
15. L. T. A. Ho and L. F. Chibotaru, Phys. Rev. B, 2018, 97, 024427.
16. A. Lunghi and S. Sanvito, Sci. Adv., 2019, 5, eaax7163.
17. N. J. Higdon, A. T. Barth, P. T. Kozlowski and R. G. Hadt, J. Chem. Phys., 2020, 152, 204306.
18. L. Gu and R. Wu, Phys. Rev. Lett., 2020, 125, 117203.
19. L. Escalera-Moreno, J. J. Baldoví, A. Gaita-Ariño and E. Coronado, Chem. Sci., 2018, 9, 3265–3275.
20. R. Orbach, Proc. R. Soc. London, Ser. A, 1961, 264, 458–484.
21. K. W. H. Stevens, Rep. Prog. Phys., 1967, 30, 189.
22. A. Abragam and B. Bleaney, Electron Paramagnetic Resonance of Transition Ions, Oxford University Press, 1970.
23. K. N. Shrivastava, Phys. Status Solidi B, 1983, 117, 437–458.
24. N. Ishikawa, M. Sugita, T. Ishikawa, S.-Y. Koshihara and Y. Kaizu, J. Am. Chem. Soc., 2003, 125, 8694–8695.
25. D. Parker, E. A. Suturina, I. Kuprov and N. F. Chilton, Acc. Chem. Res., 2020, 53, 1520–1534.
26. W. J. A. Blackmore, G. K. Gransbury, P. Evans, J. G. C. Kragskow, D. P. Mills and N. F. Chilton, Phys. Chem. Chem. Phys., 2023 10.1039/D3CP01278F.
27. M. T. Dove, Introduction to lattice dynamics, Cambridge University Press, Cambridge, New York, 1993.
28. A. Chiesa, F. Cugini, R. Hussain, E. Macaluso, G. Allodi, E. Garlatti, M. Giansiracusa, C. A. P. Goodwin, F. Ortu, D. Reta, J. M. Skelton, T. Guidi, P. Santini, M. Solzi, R. De Renzi, D. P. Mills, N. F. Chilton and S. Carretta, Phys. Rev. B, 2020, 101, 174402.
29. L. Cheng, Q.-B. Yan and M. Hu, Phys. Chem. Chem. Phys., 2017, 19, 21714–21721.
30. A. Togo, L. Chaput and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 094306.
31. D. Reta, J. G. Kragskow and N. F. Chilton, J. Am. Chem. Soc., 2021, 143, 5943–5950.
32. A. Lunghi, F. Totti, R. Sessoli and S. Sanvito, Nat. Commun., 2017, 8, 14620.
33. R. Nabi, J. K. Staab, A. Mattioni, J. G. C. Kragskow, D. Reta, J. M. Skelton and N. F. Chilton, Accurate and efficient spin-phonon coupling and spin dynamics calculations for molecular solids, ChemRxiv, 2023, preprint,  DOI:10.26434/chemrxiv-2023-6q06z-v2.
34. D. H. Woen and W. J. Evans, in Handbook on the Physics and Chemistry of Rare Earths, ed. J.-C. G. Bünzli and V. K. Pecharsky, Elsevier, 2016, vol. 50 of Including Actinides, pp. 337–394.
35. B. N. Figgis and M. A. Hitchman, Ligand Field Theory and Its Applications, Wiley VCH, New York, 2000.
36. K. W. H. Stevens, Proc. R. Soc. London, Ser. A, 1952, 65, 209–215.
37. I. Ryabov, J. Magn. Reson., 1999, 140, 141–145.
38. L. Escalera-Moreno, N. Suaud, A. Gaita-Ariño and E. Coronado, J. Phys. Chem. Lett., 2017, 8, 1695–1700.
39. N. P. Kazmierczak and R. G. Hadt, J. Am. Chem. Soc., 2022, 144, 20804–20814.
40. J. K. Staab and N. F. Chilton, J. Chem. Theory Comput., 2022, 18, 6588–6599.
41. G. Li Manni, I. F. Galván, A. Alavi, F. Aleotti, F. Aquilante, J. Autschbach, D. Avagliano, A. Baiardi, J. J. Bao, S. Battaglia, L. Birnoschi, A. Blanco-González, S. I. Bokarev, R. Broer, R. Cacciari, P. B. Calio, R. K. Carlson, R. Carvalho Couto, L. Cerdán, L. F. Chibotaru, N. F. Chilton, J. R. Church, I. Conti, S. Coriani, J. Cuéllar-Zuquin, R. E. Daoud, N. Dattani, P. Decleva, C. de Graaf, M. G. Delcey, L. De Vico, W. Dobrautz, S. S. Dong, R. Feng, N. Ferré, M. Filatov(Gulak), L. Gagliardi, M. Garavelli, L. González, Y. Guan, M. Guo, M. R. Hennefarth, M. R. Hermes, C. E. Hoyer, M. Huix-Rotllant, V. K. Jaiswal, A. Kaiser, D. S. Kaliakin, M. Khamesian, D. S. King, V. Kochetov, M. Krośnicki, A. A. Kumaar, E. D. Larsson, S. Lehtola, M.-B. Lepetit, H. Lischka, P. López Ríos, M. Lundberg, D. Ma, S. Mai, P. Marquetand, I. C. D. Merritt, F. Montorsi, M. Mörchen, A. Nenov, V. H. A. Nguyen, Y. Nishimoto, M. S. Oakley, M. Olivucci, M. Oppel, D. Padula, R. Pandharkar, Q. M. Phung, F. Plasser, G. Raggi, E. Rebolini, M. Reiher, I. Rivalta, D. Roca-Sanjuán, T. Romig, A. A. Safari, A. Sánchez-Mansilla, A. M. Sand, I. Schapiro, T. R. Scott, J. Segarra-Martí, F. Segatta, D.-C. Sergentu, P. Sharma, R. Shepard, Y. Shu, J. K. Staab, T. P. Straatsma, L. K. Sørensen, B. N. C. Tenorio, D. G. Truhlar, L. Ungur, M. Vacher, V. Veryazov, T. A. Voß, O. Weser, D. Wu, X. Yang, D. Yarkony, C. Zhou, J. P. Zobel and R. Lindh, J. Chem. Theory Comput., 2023 DOI:10.1021/acs.jctc.3c00182.
42. J. G. C. Kragskow, J. Marbey, C. D. Buch, J. Nehrkorn, M. Ozerov, S. Piligkos, S. Hill and N. F. Chilton, Nat. Commun., 2022, 13, 825.
43. J. J. Sakurai, Modern Quantum Mechanics Revised Edition, Addison-Wesley, Reading, Massachusetts, 1st edn, 1994.
44. H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, 2002.
45. U. Weiss, Quantum Dissipative Systems, WORLD SCIENTIFIC, 4th edn, 2012.
46. Á. Rivas and S. Huelga, Open Quantum Systems: An Introduction, Springer Berlin Heidelberg, 2011.
47. T. C. Berkelbach and M. Thoss, J. Chem. Phys., 2020, 152, 020401.
48. J. Prior, A. W. Chin, S. F. Huelga and M. B. Plenio, Phys. Rev. Lett., 2010, 105, 050404.
49. A. W. Chin, Á. Rivas, S. F. Huelga and M. B. Plenio, J. Math. Phys., 2010, 51, 092109.
50. U. Schollwöck, Ann. Phys., 2011, 326, 96–192.
51. N. Makri, Chem. Phys. Lett., 1992, 193, 435–445.
52. M. Topaler and N. Makri, J. Chem. Phys., 1994, 101, 7500–7519.
53. Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn., 1989, 58, 101–114.
54. Y. Tanimura, Phys. Rev. A: At., Mol., Opt. Phys., 1990, 41, 6676–6687.
55. A. Strathearn, P. Kirton, D. Kilda, J. Keeling and B. W. Lovett, Nat. Commun., 2018, 9, 3322.
56. A. Redfield, in Advances in Magnetic Resonance, ed. J. S. Waugh, Academic Press, 1965, vol. 1 of Advances in Magnetic and Optical Resonance, pp. 1–32.
57. A. Nazir and G. Schaller, in The Reaction Coordinate Mapping in Quantum Thermodynamics, ed. F. Binder, L. A. Correa, C. Gogolin, J. Anders and G. Adesso, Springer International Publishing, Cham, 2018, pp. 551–577.
58. J. Roden, W. T. Strunz, K. B. Whaley and A. Eisfeld, J. Chem. Phys., 2012, 137, 204110.
59. F. Mascherpa, A. Smirne, A. D. Somoza, P. Fernández-Acebal, S. Donadi, D. Tamascelli, S. F. Huelga and M. B. Plenio, Phys. Rev. A, 2020, 101, 052108.
60. R. Silbey and R. A. Harris, J. Chem. Phys., 1984, 80, 2615–2617.
61. D. P. S. McCutcheon and A. Nazir, New J. Phys., 2010, 12, 113042.
62. S. Mondal and A. Lunghi, Spin-phonon decoherence in solid-state paramagnetic defects from first principles, 2023.
63. M. H. Alexander, G. E. Hall and P. J. Dagdigian, J. Chem. Educ., 2011, 88, 1538–1543.
64. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
65. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
66. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson and M. C. Payne, Zeitschrift für Kristallographie, 2009, 220, 567–570.
67. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöβ, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103.
68. A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith-Daga, A. Cossard, C. Ribaldone, E. Ascrizzi, N. L. Marana, J.-P. Flament and B. Kirtman, J. Chem. Theory Comput., 2022 DOI:10.1021/acs.jctc.2c00958.
69. S. Grimme, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 211–228.
70. V. S. Parmar, A. M. Thiel, R. Nabi, G. K. Gransbury, M. S. Norre, P. Evans, S. C. Corner, J. M. Skelton, N. F. Chilton, D. P. Mills and J. Overgaard, Chem. Commun., 2023, 59, 2656–2659.
71. I. Pallikara, P. Kayastha, J. M. Skelton and L. D. Whalley, Electron. Struct., 2022, 4, 033002.
72. V. I. Anisimov, J. Zaanen and O. K. Andersen, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 44, 943–954.
73. M. Yu, S. Yang, C. Wu and N. Marom, npj Comput. Mater., 2020, 6, 1–6.
74. A. Ghosh, J. Biol. Inorg. Chem., 2006, 11, 671–673.
75. B. O. Roos, P. R. Taylor and P. E. Siegbahn, Chem. Phys., 1980, 48, 157–173.
76. I. F. Galván, M. Vacher, A. Alavi, C. Angeli, F. Aquilante, J. Autschbach, J. J. Bao, S. I. Bokarev, N. A. Bogdanov, R. K. Carlson, L. F. Chibotaru, J. Creutzberg, N. Dattani, M. G. Delcey, S. S. Dong, A. Dreuw, L. Freitag, L. M. Frutos, L. Gagliardi, F. Gendron, A. Giussani, L. González, G. Grell, M. Guo, C. E. Hoyer, M. Johansson, S. Keller, S. Knecht, G. Kovačević, E. Källman, G. Li Manni, M. Lundberg, Y. Ma, S. Mai, J. P. Malhado, P. A. Malmqvist, P. Marquetand, S. A. Mewes, J. Norell, M. Olivucci, M. Oppel, Q. M. Phung, K. Pierloot, F. Plasser, M. Reiher, A. M. Sand, I. Schapiro, P. Sharma, C. J. Stein, L. K. Sørensen, D. G. Truhlar, M. Ugandi, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, O. Weser, T. A. Wesołowski, P. O. Widmark, S. Wouters, A. Zech, J. P. Zobel and R. Lindh, J. Chem. Theory Comput., 2019, 15, 5925–5964.
77. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
78. E. Epifanovsky, A. T. B. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons, A. L. Dempwolff, Z. Gan, D. Hait, P. R. Horn, L. D. Jacobson, I. Kaliman, J. Kussmann, A. W. Lange, K. U. Lao, D. S. Levine, J. Liu, S. C. McKenzie, A. F. Morrison, K. D. Nanda, F. Plasser, D. R. Rehn, M. L. Vidal, Z.-Q. You, Y. Zhu, B. Alam, B. J. Albrecht, A. Aldossary, E. Alguire, J. H. Andersen, V. Athavale, D. Barton, K. Begam, A. Behn, N. Bellonzi, Y. A. Bernard, E. J. Berquist, H. G. A. Burton, A. Carreras, K. Carter-Fenk, R. Chakraborty, A. D. Chien, K. D. Closser, V. Cofer-Shabica, S. Dasgupta, M. de Wergifosse, J. Deng, M. Diedenhofen, H. Do, S. Ehlert, P.-T. Fang, S. Fatehi, Q. Feng, T. Friedhoff, J. Gayvert, Q. Ge, G. Gidofalvi, M. Goldey, J. Gomes, C. E. González-Espinoza, S. Gulania, A. O. Gunina, M. W. D. Hanson-Heine, P. H. P. Harbach, A. Hauser, M. F. Herbst, M. Hernández Vera, M. Hodecker, Z. C. Holden, S. Houck, X. Huang, K. Hui, B. C. Huynh, M. Ivanov, D. Jász, H. Ji, H. Jiang, B. Kaduk, S. Kähler, K. Khistyaev, J. Kim, G. Kis, P. Klunzinger, Z. Koczor-Benda, J. H. Koh, D. Kosenkov, L. Koulias, T. Kowalczyk, C. M. Krauter, K. Kue, A. Kunitsa, T. Kus, I. Ladjánszki, A. Landau, K. V. Lawler, D. Lefrancois, S. Lehtola, R. R. Li, Y.-P. Li, J. Liang, M. Liebenthal, H.-H. Lin, Y.-S. Lin, F. Liu, K.-Y. Liu, M. Loipersberger, A. Luenser, A. Manjanath, P. Manohar, E. Mansoor, S. F. Manzer, S.-P. Mao, A. V. Marenich, T. Markovich, S. Mason, S. A. Maurer, P. F. McLaughlin, M. F. S. J. Menger, J.-M. Mewes, S. A. Mewes, P. Morgante, J. W. Mullinax, K. J. Oosterbaan, G. Paran, A. C. Paul, S. K. Paul, F. Pavošević, Z. Pei, S. Prager, E. I. Proynov, D. Rák, E. Ramos-Cordoba, B. Rana, A. E. Rask, A. Rettig, R. M. Richard, F. Rob, E. Rossomme, T. Scheele, M. Scheurer, M. Schneider, N. Sergueev, S. M. Sharada, W. Skomorowski, D. W. Small, C. J. Stein, Y.-C. Su, E. J. Sundstrom, Z. Tao, J. Thirman, G. J. Tornai, T. Tsuchimochi, N. M. Tubman, S. P. Veccham, O. Vydrov, J. Wenzel, J. Witte, A. Yamada, K. Yao, S. Yeganeh, S. R. Yost, A. Zech, I. Y. Zhang, X. Zhang, Y. Zhang, D. Zuev, A. Aspuru-Guzik, A. T. Bell, N. A. Besley, K. B. Bravaya, B. R. Brooks, D. Casanova, J.-D. Chai, S. Coriani, C. J. Cramer, G. Cserey, A. E. DePrince, R. A. DiStasio, A. Dreuw, B. D. Dunietz, T. R. Furlani, W. A. Goddard, S. Hammes-Schiffer, T. Head-Gordon, W. J. Hehre, C.-P. Hsu, T.-C. Jagau, Y. Jung, A. Klamt, J. Kong, D. S. Lambrecht, W. Liang, N. J. Mayhall, C. W. McCurdy, J. B. Neaton, C. Ochsenfeld, J. A. Parkhill, R. Peverati, V. A. Rassolov, Y. Shao, L. V. Slipchenko, T. Stauch, R. P. Steele, J. E. Subotnik, A. J. W. Thom, A. Tkatchenko, D. G. Truhlar, T. Van Voorhis, T. A. Wesolowski, K. B. Whaley, H. L. Woodcock, P. M. Zimmerman, S. Faraji, P. M. W. Gill, M. Head-Gordon, J. M. Herbert and A. I. Krylov, J. Chem. Phys., 2021, 155, 084801.
79. P. A. Malmqvist, B. O. Roos and B. Schimmelpfennig, Chem. Phys. Lett., 2002, 357, 230–240.
80. L. Ungur and L. F. Chibotaru, Chem. – Eur. J., 2017, 23, 3708–3718.
81. J. P. Malrieu, R. Caballol, C. J. Calzado, C. de Graaf and N. Guihéry, Chem. Rev., 2014, 114, 429–492.
82. K. Andersson, P.-Å. Malmqvist, B. O. Roos, A. J. Sadlej and K. Wolinski, J. Phys. Chem., 1990, 94, 5483–5488.
83. K. Andersson, P.-Å. Malmqvist and B. O. Roos, J. Chem. Phys., 1992, 96, 1218–1226.
84. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger and J.-P. Malrieu, J. Chem. Phys., 2001, 114, 10252–10264.
85. in Handbook of Computational Chemistry, ed. J. Leszczynski, Springer, Netherlands, 2012.
86. C. M. Breneman and K. B. Wiberg, J. Comput. Chem., 1990, 11, 361–373.
87. E. R. Smith and J. S. Rowlinson, Proc. R. Soc. London, Ser. A, 1981, 375, 475–505.
88. C. A. P. Goodwin, F. Ortu, D. Reta, N. F. Chilton and D. P. Mills, Nature, 2017, 548, 439–442.
89. P. Evans, D. Reta, G. F. Whitehead, N. F. Chilton and D. P. Mills, J. Am. Chem. Soc., 2020, 141, 19935–19940.
90. A. Lunghi and S. Sanvito, J. Phys. Chem. Lett., 2020, 11, 6273–6278.
91. F. Bernardi, M. Olivucci and M. A. Robb, Chem. Soc. Rev., 1996, 25, 321.
92. I. F. Galván, M. G. Delcey, T. B. Pedersen, F. Aquilante and R. Lindh, J. Chem. Theory Comput., 2016, 12, 3636–3653.
93. F. Plasser, S. Gómez, M. F. S. J. Menger, S. Mai and L. González, Phys. Chem. Chem. Phys., 2019, 21, 57–69.
94. J. P. Zobel, M. Heindl, F. Plasser, S. Mai and L. González, Acc. Chem. Res., 2021, 54, 3760–3771.
95. H. Köppel, W. Domcke and L. S. Cederbaum, Advances in Chemical Physics, John Wiley & Sons, Inc., 1984, pp. 59–246.
96. G. A. Worth and L. S. Cederbaum, Ann. Rev. Phys. Chem., 2004, 55, 127–158.
97. B. H. Lengsfield III and D. R. Yarkony, Advances in Chemical Physics, John Wiley & Sons, Inc., 1992, pp. 1–71.
98. H. Lischka, M. Dallos, P. G. Szalay, D. R. Yarkony and R. Shepard, J. Chem. Phys., 2004, 120, 7322–7329.
99. M. S. Schuurman and D. R. Yarkony, J. Chem. Phys., 2007, 127, 094104.
100. Y. Nishimoto, J. Chem. Phys., 2021, 154, 194103.
101. Y. Nishimoto, S. Battaglia and R. Lindh, J. Chem. Theory Comput., 2022, 18, 4269–4281.

---

Footnotes

† Electronic supplementary information (ESI) available: Details of DFT calculations for CO₂, NaCl and crystalline NH₃, derivation of magnetic relaxation rates. See DOI: https://doi.org/10.1039/d2cs00705c

‡ Note that care must be taken with numerical diagonalisation of Γ where eigenvalues span many orders of magnitude and limitations due to finite numerical precision become important. For example, magnetic relaxation rates are usually on the order of 10⁵ s⁻¹ to 10⁻⁵ s⁻¹, while the “fast fluctuations” can be ∼10¹³ s⁻¹, which sets the numerical value of “zero” as ∼10⁻³ s⁻¹ in double precision arithmetic, which clearly poses a problem for calculation of rates at low temperature. However, using quadruple precision arithmetic, numerical “zero” would be ∼10⁻²¹ s⁻¹ in the presence of eigenvalues of ∼10¹³ s⁻¹.

§ https://pypi.org/project/molcas-suite/, https://pypi.org/project/angmom-suite/, https://pypi.org/project/spin-phonon-suite/

¶ https://gitlab.com/chilton-group/tau

|| https://github.com/LunghiGroup/MolForge

---