We present a symmetry-adapted perturbation theory (SAPT) for the interaction of two high-spin open-shell molecules (described by their restricted open-shell Hartree-Fock determinants) resulting in low-spin states of the complex. The previously available SAPT formalisms, except for some system-specific studies for few-electron complexes, were restricted to the high-spin state of the interacting system. Thus, the new approach provides, for the first time, a SAPT-based estimate of the splittings between different spin states of the complex. We have derived and implemented the lowest-order SAPT term responsible for these splittings, that is, the first-order exchange energy. We show that within the so-called *S*^{2} approximation commonly used in SAPT (neglecting effects that vanish as fourth or higher powers of intermolecular overlap integrals), the first-order exchange energies for all multiplets are linear combinations of two matrix elements: a *diagonal exchange term* that determines the spin-averaged effect and a *spin-flip term* responsible for the splittings between the states. The numerical factors in this linear combination are determined solely by the Clebsch-Gordan coefficients: accordingly, the *S*^{2} approximation implies a Heisenberg Hamiltonian picture with a single coupling strength parameter determining all the splittings. The new approach is cast into both molecular-orbital and atomic-orbital expressions: the latter enable an efficient density-fitted implementation. We test the newly developed formalism on several open-shell complexes ranging from diatomic systems (Li⋯H, Mn⋯Mn, …) to the phenalenyl dimer.

## I. INTRODUCTION

The interactions of open-shell atoms, molecules, and ions are of great significance to many areas of molecular physics and chemistry. As a matter of fact, all chemical reactions at some stage have to experience bond breaking and rearrangement and to describe such processes one has to deal with interacting fragments of reactants or products. Many electronically excited states of molecules have open electronic shells, so in studies of the interactions of such molecules with background gas or solvent, similar problems arise. Interacting radicals are of key importance in atmosphere, where they undergo chain reactions as well as form metastable states and stable complexes.^{1,2} Collisions of radicals (such as OH, CH_{2}, CN, and many more) with hydrogen (H and H_{2}) are a subject of study for astrochemistry^{3} as they provide information about the conditions which molecules experience in the interstellar clouds. Stabilized organic radicals are relevant in organic and materials chemistry,^{4} for example, as building blocks of molecular magnets. One should also mention the high relevance of interactions of radicals in cold chemistry: open-shell molecules can be manipulated with magnetic fields to control their collisions at low temperatures, giving unique information on the interaction potential.^{5–7}

In modeling the dynamics of weakly interacting systems, accurate potential energy surfaces are crucial. In the case of open-shell systems, despite large progress in recent years, electronic structure calculations are still far from routine or robust. Although the coupled cluster (CC) method, most widely used in studying molecular interactions [in particular in its variant with single, double, and perturbative triple excitations, CCSD(T)], has been extended to the high-spin open-shell case a while ago by Knowles *et al.*^{8} and Watts *et al.*,^{9} there are many examples where the single reference coupled cluster method fails and a multireference approach is necessary. In particular, one of the most important cases for which the single-reference CC method does not work is the dissociation of a low-spin system into high-spin subsystems: for example, a breakdown of a singlet, stable molecule into two doublet molecules. However, despite many attempts and progress in theory since the introduction of the multireference coupled-cluster (MRCC) ansatz by Jeziorski and Monkhorst,^{10} the proposed variations of MRCC are far from being black-box methods (see Refs. 11 and 12 for a review). At present, to account for both static and dynamical correlation effects in molecular interactions for systems exhibiting a non-unique electronic configuration, multireference configuration interaction^{13} (MRCI) or multireference perturbation theories^{14–16} are often applied. In the case of the MRCI method, it is possible to reproduce the short-range part of the potential well when size-consistency corrections are applied. Such *a posteriori* corrections are not exact and non-unique, and the remaining size-consistency error might be on the order of the interaction energy. In the case of the perturbation methods, the same problem also arises in general; moreover, some of these methods suffer from the problem of intruder states.^{17} One should also mention that all multireference methods rely on a proper selection of the active space, which is not necessarily straightforward, and that the complexity of calculations grows exponentially with the size of the active space. This is a particularly severe problem for interacting open-shell molecules as the active space required for the complex, typically the union of the monomers’ active spaces, might be intractable.

An interesting alternative for open-shell systems is the family of the so-called spin-flip (SF) methods in which it is possible to reach arbitrary, multireference low-spin states of a given system starting from a single-reference high-spin state.^{18} The spin-flip equation of motion approach has been implemented within a variety of different theories, in particular, Hartree-Fock,^{18} density functional theory,^{19} configuration interaction doubles,^{18} and finally, using the CC framework.^{20,21} The theory has been generalized to multiple spin flips;^{22} however, the resulting calculations can be quite demanding, and it is often advantageous to describe the exchange splittings between all multiplets by a single exchange coupling constant *J*_{AB} within the Heisenberg spin Hamiltonian. In order to determine *J*_{AB}, it is sufficient to perform a single spin-flip calculation starting from the high-spin, single-determinant reference state.^{23}

Recent progress notwithstanding, the current arsenal of methods which can be used for multireference open-shell complexes is limited, and it is desirable to explore alternative approaches which provide good insight and reliable interaction potentials. The closed-shell symmetry-adapted perturbation theory (SAPT)^{24,25} has been widely recognized as a highly useful tool for calculations of interaction energies with an accuracy comparable to CCSD(T), but also as a robust analysis tool which provides insightful physical interpretation of the nature of the interaction in terms of its electrostatic, induction, dispersion, and exchange components. For general single-reference open-shell systems, Żuchowski *et al.*^{26} and Hapka *et al.*^{27} introduced the symmetry-adapted perturbation theory based on, respectively, restricted and unrestricted Slater determinants (with Hartree-Fock or Kohn-Sham orbitals). Present implementations of SAPT, however, are valid only for the single-reference high-spin case of the dimer in which the spin quantum number of the complex is equal to the sum of spin quantum numbers of monomers (*S* = *S*_{A} + *S*_{B}).

The interaction of two open-shell species with the total spin of the dimer smaller than the sum of spin quantum numbers of the monomers still poses a challenge. This is the case since low-spin configurations are typical examples of multireference systems. The knowledge of the so-called exchange splitting, which is defined as the difference between the highest and lowest possible spin state of a given dimer, is very important in interactions of small open-shell molecules (e.g., in the O_{2} $\Sigma g\u22123$ dimer^{28–30}) as well as for interactions involving stable organic radicals such as phenalenyl^{31} or unsaturated metal sites within metal organic frameworks in the O_{2} adsorption process.^{32} The new approaches of calculating low-spin potential energy surfaces could be applied to model fragments of interaction potentials in reaction dynamics, for instance, near the channels which correspond to dissociation into two radicals or near the bond-breaking geometry of the reactive complex. Spin splitting also plays an important role in cold physics since it is the term that strongly couples hyperfine states of alkali-metal atoms and drives Feshbach resonances.^{33,34} It should be stressed that in the SAPT approach exchange splitting does indeed come exclusively from exchange corrections: the SAPT terms arising from polarization theory^{24} are identical for all asymptotically degenerate states of different spin multiplicity.

This paper is the first step toward the development of SAPT for low-spin dimer states. Here, we derive the first-order exchange energy in terms of spin-restricted orbitals which preserve the correct value of the squared-spin operator $S^2$. The theory introduced in this paper is beneficial for several reasons: (i) as usually in SAPT, the interaction energy components are calculated *directly*, which contrasts with supermolecular calculations that always involve subtraction; (ii) there is no need for the active space selection for the complex: in fact, no multireference calculation is necessary (similar to the spin-flip electronic structure approaches); (iii) it is possible to calculate the energy in the monomer-centered basis set, and the overall scaling of the method [once transformed molecular-orbital (MO) integrals are available] is just *o*^{4}, where *o* is the number of occupied orbitals; and (iv) the low-spin exchange energy component, together with other SAPT corrections, can provide an insight into the nature of the interaction in a radical-radical system. The formalism presented here is valid for any spin-restricted Slater determinants. At this level, similarly to the high-spin open-shell theory,^{26} intramonomer correlation effects are not included except for the possible use of Kohn-Sham orbitals.^{35} As we will show in Sec. II, the first-order SAPT term responsible for the multiplet splittings involves matrix elements between the zeroth-order wavefunction and a function where one of the unpaired spins on one monomer has been lowered and one of the unpaired spins on the other monomer has been raised, in other words, a function which was obtained from the zeroth-order SAPT wavefunction via an intermolecular spin flip. Therefore, we will refer to our new approach as spin-flip SAPT (SF-SAPT), borrowing the nomenclature from spin-flip electronic structure theories.^{18} However, while in the latter theories the *z* component of the total spin of the system (the *M*_{S} quantum number) is changed upon the spin flip, in SF-SAPT only the $MSA$ and $MSB$ numbers for individual molecules are changed: the overall *M*_{S} value is conserved.

To date, there are only a few applications of perturbation theory to low-spin complexes. In 1984, Wormer and van der Avoird used the Heitler-London formula with appropriately spin-adapted linear combinations of Slater determinants (in the specific case of 4 active electrons) to predict the splittings between different multiplets of the O_{2}⋯O_{2} system.^{36} Extensive studies of SAPT performance up to very high order of perturbation theory exist for few-electron systems in which monomers can be represented by exact (in a given Gaussian basis) wavefunctions. Ćwiok *et al.*^{37} studied the singlet and triplet states of the H⋯H interaction, while Patkowski *et al.* reported the convergence of various SAPT variants for the lowest singlet and triplet states of the Li⋯H complex.^{38,39}

In the past few years, the exchange splitting of interaction energy has been studied for the $H2+$ system by Gniewek and Jeziorski.^{40–42} Although this is a single-electron system, its two lowest states are asymptotically degenerate and correspond to symmetric (*g*) and antisymmetric (*u*) combinations of atomic orbitals (AOs). The difference between *u* and *g* states results from the resonance tunneling of the electron between the nuclei, in a similar way as the difference between triplet and singlet states in the interaction of two doublets. Gniewek and Jeziorski used perturbation theory to study the asymptotic expansion of exchange energy in $H2+$ and discussed a new approach to exchange energy via a variational volume-integral formula.

The plan of this paper is as follows: in Sec. II, we derive the arbitrary-spin first-order exchange energy formula for two monomers described by their spin-restricted determinants and show the connection with the Heisenberg model of the scalar interaction of two spins. In Sec. III, we give the most important details regarding the implementation of the theory, and in Sec. IV we report the test calculations for several important open-shell⋯open-shell systems. In Sec. V of the paper, we summarize our results and discuss prospects for the future.

## II. THEORY

Similar to the existing closed-shell and high-spin open-shell SAPT approaches, the permutational symmetry adaptation, leading to the exchange corrections, is performed in the simplest manner, within the symmetrized Rayleigh-Schrödinger (SRS) formalism.^{24,43} Specifically, the SRS wavefunction corrections are computed from the standard Rayleigh-Schrödinger (RS) perturbation expansion without any symmetry adaptation: only later, the energy corrections are computed from a formula that involves a symmetry projector. In the spin-free formalism used, for example, in Ref. 38, this symmetry operator projects onto a particular irreducible representation of the permutation group $SNA+NB$ (*N*_{X} is the number of electrons of molecule *X*): the choice of this representation is determined by the spin multiplicity. In the more conventional, spin formalism employed here (and used in the high-spin case so far), this operator has to be replaced by the (*N*_{A} + *N*_{B})-electron antisymmetrizer and an appropriate spin projection. In the already existing high-spin open-shell SAPT theories, the monomer wavefunctions Ψ_{A} and Ψ_{B} as well as the zeroth-order dimer wavefunction Ψ_{0} = Ψ_{A}Ψ_{B} are pure spin functions so the spin projection is not needed. This is only the case if all unpaired electrons on both monomers have the same spin, that is, when *M*_{S} = ±(*S*_{A} + *S*_{B}) so that no contamination by lower-spin states is possible. Indeed, the restricted open-shell Hartree-Fock and Kohn-Sham (ROHF/ROKS) SAPT exchange corrections have been developed^{26} under the assumption that all unpaired spins point the same way.

For the low-spin case, it is known^{12,39} that the antisymmetrizer $A$ in the SRS energy expression has to be accompanied by an operator $PSMS$ that projects a dimer function onto a subspace corresponding to a particular value of the total spin *S* (with an appropriate *M*_{S}). In this way, the exchange corrections, different for different dimer spin states, are obtained by expanding the energy expression

in powers of the intermolecular interaction operator $V\u2009=\u2009\u2211i\u2208A\u2211j\u2208Bv\u0303(ij)$ and separating the polarization and exchange effects in each order. In Eq. (1), $A$ is the (*N*_{A} + *N*_{B})-electron antisymmetrizer and Φ_{AB} is the polarization (RS) expansion of the wavefunction for the complex.^{24} We will assume that Ψ_{A} and Ψ_{B} are ground-state high-spin open-shell determinants. The occupied orbitals in Ψ_{A} are $\chi i1$, $\chi i2$, …, $\chi ikA$ (inactive, doubly occupied) and $\chi a1$, $\chi a2$, …, $\chi a2SA$ (active, occupied by *α* electrons only). The occupied orbitals in Ψ_{B} are $\chi j1$, $\chi j2$, …, $\chi jkB$ (inactive, doubly occupied) and $\chi b1$, $\chi b2$, …, $\chi b2SB$ (active, occupied by *β* electrons only).

The spin projector $PSMS$ commutes with *V* and the zeroth-order Hamiltonian *H*_{0} = *H*_{A} + *H*_{B} as the latter operators are spin-independent. Moreover, the operator $\u015c2$, and thus also $PSMS$, is independent of the ordering of electrons and hence commutes with any electron permutations or their linear combinations including $A$. As a result, we can rewrite Eq. (1) in an alternative form where the spin projection acts in the bra instead of the ket,

While Eqs. (1) and (2) are fully equivalent, they have quite different interpretations. Equation (1) corresponds to a standard RS treatment where the wavefunction corrections in the expansion of Φ_{AB} are the same for all asymptotically degenerate multiplets: the spin projection enters later, at the calculation of SRS energy corrections. Equation (2) corresponds to a perturbation expansion initiated from a spin-adapted zeroth-order function $PSMS\Psi A\Psi B$, with all resulting perturbation corrections (in the expansion of $\Phi AB\u2032$) also being of pure spin. Importantly, both approaches allow the use of standard nondegenerate perturbation theory (unless degeneracy arises from the orbital part of the monomer wavefunctions, which is outside the scope of this work), albeit for quite different reasons. For Eq. (2), the projector $PSMS$ reduces the zeroth-order space of (2*S*_{A} + 1)(2*S*_{B} + 1) asymptotically degenerate product functions corresponding to a pair of interacting multiplets to a single spin-adapted combination with the requested (*S*, *M*_{S}). For Eq. (1), any of the (2*S*_{A} + 1)(2*S*_{B} + 1) product functions, including Ψ_{A}Ψ_{B}, is a valid zeroth-order state for nondegenerate RS perturbation theory as neither *V* nor *H*_{0} mix states with different $MSA$ or $MSB$. We choose Eq. (1) as the starting point for the derivation below.

We now introduce the single exchange approximation (often called the *S*^{2} approximation as it neglects terms higher than quadratic in the intermolecular overlap integrals),

where the single-exchange operator $P$ is given by

(for non-approximated alternatives in the closed-shell case; see Refs. 44–46). We now multiply both sides of Eq. (1) by the denominator, keeping only the terms of first order in *V* (thus, replacing Φ_{AB} by its zeroth-order term Ψ_{A}Ψ_{B}), and make the approximation (3). Introducing the shorthand notation ⟨*X*⟩ ≡ ⟨Ψ_{A}Ψ_{B}|*X*|Ψ_{A}Ψ_{B}⟩, we obtain (note that $A$, $AA$, and $AB$ all commute with $PSMS$)

Note that we have separated the electrostatic contribution $Eelst(10)=\u27e8V\u27e9$ from the exchange one, and we introduced the zero into the *E*^{(10)} superscript to remind the reader that no intramolecular correlation effects are included.

Before we proceed, we need to understand a bit better the action of the operator $PSMS$ on the product Ψ_{A}Ψ_{B}. The ROHF determinants Ψ_{A} and Ψ_{B} are pure spin functions corresponding to the quantum numbers $(SA,MSA=SA)$ and $(SB,MSB=\u2212SB)$, respectively. In the zeroth-order space spanned by (2*S*_{A} + 1)(2*S*_{B} + 1) product functions that correspond to total spin *S*_{A} for molecule A and total spin *S*_{B} for molecule B (and any combination of $MSA,MSB$), there exists exactly one function $\Psi S,MS$ with a total spin *S* ∈ {|*S*_{A} − *S*_{B}|, …, *S*_{A} + *S*_{B}} and its projection *M*_{S} ∈ {−*S*, …, *S*}. This function is a linear combination of only those product functions that correspond to $MSA+MSB=MS$. Therefore, the projection $PSMS\Psi A\Psi B$ picks, up to normalization, this very function $\Psi S,MS$ out of the zeroth-order space. In other words, $PSMS\Psi A\Psi B$ produces a linear combination of products of functions of A and B with the same values *S*_{A}, *S*_{B} but different $MSA,MSB$ (however, the latter two numbers add up to the same total *M*_{S}). The coefficients in this linear combination are the Clebsch-Gordan coefficients $\u27e8SMS|SAMSASBMSB\u27e9$—we will denote the coefficients ⟨*S* (*S*_{A} − *S*_{B})|*S*_{A} *S*_{A} *S*_{B} − *S*_{B}⟩, ⟨*S* (*S*_{A} − *S*_{B})|*S*_{A} (*S*_{A} − 1) *S*_{B} (−*S*_{B} + 1)⟩, and ⟨*S* (*S*_{A} − *S*_{B})|*S*_{A} (*S*_{A} − 2) *S*_{B} (−*S*_{B} + 2)⟩ by *c*_{0}, *c*_{1}, and *c*_{2}, respectively. Specifically, the spin projection can be expressed as follows:

where, for example, $\Psi A\u2193$ is a normalized wavefunction (linear combination of determinants) corresponding to the spin quantum numbers (*S*_{A}, *S*_{A} − 1). Up to a constant, this function can be obtained from Ψ_{A} by the action of the spin-lowering operator $\u015c\u2212$.

Equation (6) implies that

Therefore, Eq. (5) becomes

Now, the last term on the l.h.s. is neglected as it requires at least two electron exchanges (one in $Eexch(10)$ and one in $\u27e8PPSMS\u27e9$), and we arrive at the final formula for the SRS-like first-order exchange energy,

Let us now go back to Eq. (6) and examine the spin-flipped monomer wavefunctions such as $\Psi A\u2193$. Up to normalization, this function is equal to $\u015c\u2212$Ψ_{A}, where $\u015c\u2212=\u2211k=1NA\u015c\u2212(k)$ applies the conventional spin-lowering operator $\u015c\u2212$(*k*) = $\u015cx$(*k*) − $i\u015cy$(*k*) to each of the *N*_{A} spins. It can be shown that when Ψ_{A} is a restricted high-spin determinant like in this work, the lowering operator acts on it as follows:

where the spin orbitals present in each (normalized) determinant are explicitly listed. Note that only the active spin orbitals end up being spin-flipped. This observation is true for the ROHF determinants only: lowering of the spin for an inactive spin orbital results in a duplicate spin orbital and the determinant vanishes. If unrestricted Hartree-Fock (UHF) determinants were considered, this simplification would not take place. Analogously, the repeated application of the spin-lowering operator, $\u015c\u2212\u015c\u2212$Ψ_{A}, gives a linear combination of all possible determinants where two active spin orbitals have been flipped from *α* to *β*. The function $\u015c\u2212$Ψ_{A}, Eq. (11), is not normalized, but all 2*S*_{A} determinants entering the linear combination are clearly orthonormal. Therefore, $\Psi A\u2193=(1/2SA)\u015c\u2212\Psi A$ is normalized.

The first interesting observation from Eq. (12) is that the standard, closed-shell like contribution (the “spin-diagonal” term) $\u27e8VP\u27e9\u2212\u27e8V\u27e9\u27e8P\u27e9$ is present in the first-order exchange energy of any dimer spin state regardless of the values of the Clebsch-Gordan coefficients. It is the off-diagonal “spin-flip” terms that are responsible for the splittings (and, as we will see below, only the single spin-flip term survives). Starting with the diagonal term and using the density-matrix approach to exchange energy that is valid in monomer- as well as dimer-centered basis sets, we can use the general spin orbital expression given by Eq. (39) of Moszyński *et al.*,^{47}

In Eq. (13) and throughout the rest of the text, summation over all repeated lower and upper indices is implied. Furthermore,

$(\lambda \mu |\kappa \nu )$ is a two-electron integral in the (11|22) notation, $(vX)\nu \kappa $ is the nuclear attraction integral with all nuclei of X = A,B, $V0$ is the constant intermolecular nuclear repulsion term, and $S\mu \lambda $ is the overlap integral. In Ref. 47, the indices *α*, *α*′ (*β*, *β*′) run over all occupied spin orbitals of A (B). Equation (14) applies in the open-shell case as well (for the diagonal term) but we have to break up the summation over, e.g., the occupied spin orbitals of A into the inactive orbitals with spins *α* and *β* and active orbitals with spin *α*. When we do that and perform all spin integrations, the resulting expression for the diagonal term is

where the implicit summations now run over orbitals: *i*–(inactive A), *j*–(inactive B), *a*–(active A), and *b*–(active B). It is worth noting that an active-active exchange term $v\u0303abba$ is absent from Eq. (15) due to the fact that all active orbitals on A are paired with *α* spins and all active orbitals on B are paired with *β* spins.

For the first off-diagonal term, corresponding to a single spin flip between the interacting monomers, Eq. (12) involves two matrix elements, $\u27e8\Psi A\Psi B|VP|\Psi A\u2193\Psi B\u2191\u27e9$ and $\u27e8\Psi A\Psi B|P|\Psi A\u2193\Psi B\u2191\u27e9$. Analogous to the conventional closed-shell formalism, these elements can be expressed through the “spin-flip interaction density matrix,”

where

and, as always, $d\tau 12\u2032$ means the integration over coordinates of all electrons except 1 and 2. The spin-flip interaction density matrix can be expressed by the (also spin-flip) one- and two-electron reduced density matrices in full analogy to Eq. (99) of Ref. 48,

The spin-flip reduced density matrices on the r.h.s. of Eq. (19) are in fact no different from ordinary transition density matrices: for example, for monomer A,

Therefore, in order to obtain spin orbital expressions for these matrices, we can use Eqs. (92) and (94) of Ref. 48, except that a general *α* → *ρ* excitation is replaced by a linear combination of spin-flip excitations. Therefore, we can straightforwardly use Eq. (92) of Ref. 48 to write

The application of Eq. (94) from Ref. 48 to obtain a formula for $\Gamma A\u2193$ is a little more complicated because there is an additional summation over an occupied orbital in this equation. This summation needs to be split into three: over inactive orbitals with spin *α*, over inactive orbitals with spin *β*, and over active orbitals (which have spin *α*). The resulting expression is

Developing analogous formulas for the spin-flip reduced density matrices of monomer B, employing Eq. (19), and performing spin integration, one arrives at the following formula:

For the renormalization term, Eq. (17), we find [for example, by reanalyzing Eq. (24), replacing each $v\u0303$ by $1NANB$ times the appropriate product of overlap integrals] that

The ROHF electrostatic energy is equal to

Combining Eqs. (24)–(26), we obtain the final formula for the single spin-flip contribution to Eq. (12),

Finally, we will show that the double spin-flip term in Eq. (12), and all subsequent terms, vanishes identically. Specifically, we will prove that $\u27e8\Psi A\Psi B|VP|\Psi A\u2193\u2193\Psi B\u2191\u2191\u27e9=0$, where $\Psi A\u2193\u2193$ and $\Psi B\u2191\u2191$ are proportional to $\u015c\u2212\u015c\u2212$Ψ_{A} and $\u015c+\u015c+$Ψ_{B}, respectively. Similar to the previous term,

where the double spin-flip interaction density matrix $\rho int\u2193\u2193\u2191\u2191(12)$ is given by a formula analogous to Eq. (19). However, as a double spin flip can be treated as a special case of a double excitation, the one- and two-electron reduced spin-flip density matrices in the modified Eq. (19) are now given by Eqs. (93) and (95), respectively, of Ref. 48. Consequently,

and the two-electron density matrix is proportional to

According to Eq. (29), only the last term in the expression for $\rho int\u2193\u2193\u2191\u2191$ survives,

Such an expression vanishes upon the integration over d*τ*_{1} and d*τ*_{2}, either alone or with a spin-independent operator such as $v\u0303$(12) [as in Eq. (28)]. In particular, in the integration over d*τ*_{1}, the two spin orbitals in $\Gamma A\u2193\u2193(13|14)$ that depend on the coordinates of electron 1 always occur with opposite spins [cf. Eq. (30)], leading to a zero spin integral. This concludes the proof that the double spin-flip contribution to Eq. (12) vanishes.

To conclude, the first-order SAPT exchange correction for an interaction of two high-spin open-shell molecules with spins *S*_{A} and *S*_{B}, in an arbitrary dimer spin state |*S*_{A} − *S*_{B}| ≤ *S* ≤ *S*_{A} + *S*_{B}, is given by the formula which combines the orbital expressions for $Eexch,diag(10)$ and $Eexch,1flip(10)$ defined in Eqs. (15) and (27), respectively. By simplifying the $c1/2c0SASB$ coefficient, as outlined in the Appendix, we obtain

where

Note that *Z*(*S*_{A}, *S*_{B}, *S*_{A} + *S*_{B}) = 1, while $Z(SA,SB,|SA\u2212SB|)=\u221212max(SA,SB)$. All of the resulting matrix elements are the same for the full set of asymptotically degenerate multiplets of the dimer: the only factors in Eq. (32) that are different for different spin states are *Z*(*S*_{A}, *S*_{B}, *S*). Thus, as expected from the Heisenberg Hamiltonian model, the ratios of splittings between different multiplets are simple, system-independent numbers (expressible through the Clebsch-Gordan coefficients) as long as the single exchange approximation is applied. The observation that the single exchange approximation implies the Heisenberg picture is not new: in fact, it dates back to the work of Matsen *et al.*^{49} on spin-free quantum chemistry.

For practical applications to large systems, it is advantageous to recast the newly derived $Eexch(10)$ expressions from the molecular-orbital (MO) basis into the atomic orbital (AO) one, similar to the AO first-order exchange expressions without the single-exchange approximation for closed-shell SAPT^{50} and UHF-based high-spin open-shell SAPT.^{27} The AO approach is advantageous as it avoids the AO-MO transformation of many different types of two-electron integrals present in Eqs. (15) and (27). Moreover, as we will see below, the AO expressions make heavy use of the generalized Coulomb and exchange operators, whose evaluation in the psi4 code^{51} is highly optimized both with and without density fitting (DF).

In the following, we will denote AO indices by capital letters assuming, for simplicity, that the same AO basis set has been used to expand molecular orbitals of A and B (according to the so-called dimer-centered basis set formalism^{52}). The SCF coefficients of the molecular spin orbital *λ* will be denoted by $C\lambda K$—obviously, in the ROHF formalism, the SCF coefficients are the same for spin *α* and spin *β*, for example, $Ci\alpha K=Ci\beta K$. We will use boldface letters for matrices and denote by $A\u22c5B=\u2211KLAKLBKL$ the scalar product of two matrices. The inactive and active parts of density matrices for monomer A are $PKLiA=Ci\alpha KCi\alpha L=Ci\beta KCi\beta L$, $PKLaA=Ca\alpha KCa\alpha L$, and similarly for monomer B (note that the active electrons on A all have spin *α*, the active electrons on B all have spin *β*). Therefore, the total density matrix is **P**^{A} = 2**P**^{iA} + **P**^{aA}. Now, the generalized Coulomb and exchange matrices are defined as

and reduce to standard Coulomb and exchange matrices, or their inactive/active contributions, in special cases, for example, **J**^{iA} = **J**[**P**^{iA}], **K**^{aA} = **K**[**P**^{aA}], …. The electrostatic potential matrix for monomer A is

where **v**_{A} is the matrix of the nuclear attraction operator. By adding the exchange matrices, we obtain the following spin-dependent Fock matrices **h**:

With this notation, the diagonal and single spin-flip terms of the ROHF-based $Eexch(10)$ formula, Eqs. (15) and (27), respectively, can be rewritten as follows:

where **S**^{AO} is the overlap matrix in the AO basis.

## III. COMPUTATIONAL DETAILS

The MO-based formulas for the arbitrary-spin first-order exchange energy, Eqs. (15) and (27), have been implemented in two versions to verify the numerical correctness of the codes: into the ROHF-based SAPT code of Ref. 26 forming a part of the SAPT2012 package^{53} and into the development version of the psi4 package^{51} using the straightforward psi4numpy framework^{54} in which each term corresponds to a single line of Python code. The exchange energy for the high-spin state of the dimer could also be computed with the code of Ref. 26 and we have verified that our new code gives identical results in the high-spin case. In fact, for two interacting doublets, the formulas for $Eexch,diag(10)$ and $Eexch,1flip(10)$ [Eqs. (15) and (27), respectively] add up exactly to the formula for $Eexch(10)$ in the high-spin (triplet) state in accordance with Eq. (33) for *S*_{A} = *S*_{B} = 1/2 and *S* = 1. The AO-based first-order exchange expressions, Eqs. (41) and (42), have been implemented 2into the development version of psi4, making use of its efficient evaluation of generalized Coulomb and exchange matrices, both with and without the DF approximation. Whenever the latter approximation was applied, the standard aug-cc-pV*X*Z/JKFIT set^{55} was used as the auxiliary basis accompanying the orbital set aug-cc-pV*X*Z.^{56} As shown in Sec. IV E, the DF approximation to $Eexch(10)$ performs very well, and the density fitted AO-based code can be applied to much larger systems than the ones studied here. Unless stated otherwise, we used the aug-cc-pVTZ basis set.

## IV. NUMERICAL RESULTS

In this section, we study the applications of theory developed in Sec. II to several test systems. For the first one, the Li⋯H complex, it is possible to compare the spin-flip theory with exact SAPT calculations of Patkowski *et al.*^{38} Then, we focus on other atom-atom complexes including the lithium dimer, Li⋯N, and N⋯N systems. We also examine the Mn⋯Mn system, very extensively studied in the past, for which the spin-exchange splitting is surprisingly low. Finally, we also show a test of exchange splittings for molecular cases: the O_{2}($\Sigma g\u22123$) dimer and the phenalenyl dimer.

For the atomic Li⋯Li, Li⋯N, and N⋯N systems, we can compare the energy differences between low-spin and high-spin dimers with results obtained from multireference *ab initio* methods: complete active space self-consistent field (CASSCF) and internally contracted multireference configuration interaction (MRCI) with the Davidson correction implemented in molpro 2012.^{57} The exchange splitting Δ*E* between the interaction energies for the highest (*E*_{int}(*S*_{A} + *S*_{B})) and lowest (*E*_{int}(|*S*_{A} − *S*_{B}|)) multiplicities from *ab initio* calculations can be compared with the spin-flip SAPT term, Eq. (27), by multiplying the latter by a following factor:

where from now on we use Δ*E*^{(10)} to denote the difference between the highest- and lowest-multiplicity interaction energy of a complex at the SAPT level introduced in this work.

### A. Li⋯H

The first system we investigate is the Li⋯H complex, where the interaction between two doublets gives rise to a triplet state and a singlet state. For both of them, full configuration interaction (FCI) interaction energies as well as SAPT corrections to high order (in several variants including SRS) have been obtained before.^{38} Therefore, to facilitate our comparisons with the data of Ref. 38, we use the same (4s4p1d/5s2p) basis set in our first-order exchange calculations. The exchange energies of both spin states as functions of the interatomic separation *R* are presented in Fig. 1. The excellent agreement between the results of Ref. 38 (where the electron correlation within the lithium atom was described at the FCI level) and our new values (where the Li atom was described by the ROHF theory) not only validates our implementation but also indicates that the intramonomer correlation effect on $Eexch(1)$ is negligible for this system. One should note that, as one of the monomers has only one electron, the single exchange approximation is exact in this case.

At the range of distances presented in Fig. 1, the diagonal exchange energy term, Eq. (41), is much smaller than the spin-flip term, Eq. (42). The $Eexch,diag(10)$ contribution stems from the difference between the “apparent Coulomb energy,” the arithmetic mean of the singlet and triplet energies, and the actual mathematical Coulomb energy, defined in SAPT as a weighted average of the energies of all asymptotically degenerate states contributing to the zeroth-order wavefunction, including the Pauli-forbidden ones.^{38,58} It can be shown that, in the first order, the mathematical Coulomb energy is identical to the electrostatic correction $Eelst(10)$ [Eq. (26)] up to terms that vanish as the fourth or higher powers of intermolecular overlap integrals.^{59} The apparent Coulomb energy does not have this property and the difference between the two (that is, $Eexch,diag(10)$) does not vanish in the *S*^{2} approximation. The only exceptions are the interactions involving one- and two-electron systems, such as H⋯H,^{37} He (^{3}*S*)⋯H,^{60} and He (^{3}*S*)⋯He (^{3}*S*),^{61} when the mathematical Coulomb energy coincides with the mean energy of physical states. In these cases, $Eexch,diag(10)$ [Eq. (15)] is identically zero as there are no inactive orbitals (the range of summation over *i* and *j* is empty). For the Li⋯H complex, as already observed in Ref. 38, the two Coulomb energies are particularly close and the $Eexch,diag(10)$ term is much smaller than $Eexch,1flip(10)$. As a result, the singlet and triplet exchange energies are almost exactly the negatives of each other. Assuming that this is precisely the case, one could define the exchange energy at the FCI level as one half of the difference between the FCI interaction energies of triplet and singlet. The FCI exchange energies defined in this way are also displayed in Fig. 1. Interestingly, $Eexch(10)$ represents about two thirds of the total exchange energy, showing that the first-order description of exchange is qualitatively correct. This result should be contrasted with the first-order electrostatic correction $Eelst(10)$ that is nearly negligible for this dispersion-bound complex. As a result, the sum $Eelst(10)+Eexch,diag(10)$ recovers only a small fraction (8% at the triplet van der Waals minimum distance of 11.5*a*_{0}) of the FCI apparent Coulomb energy, while $Eexch,1flip(10)$ recovers 71% of the FCI exchange energy.

### B. Li⋯Li, Li⋯N, and N⋯N

The total spin for the Li⋯Li complex can be either singlet or triplet, for Li⋯N—triplet or quintet, and for two interacting nitrogen atoms it can take multiplicities from singlet to septet. For these systems, unlike for the Li⋯H interaction, the benchmark SAPT corrections with exact monomer wavefunctions are not available. In Fig. 2, we show the comparison of exchange splittings from SF-SAPT, CASSCF, and MRCI (Davidson corrected). As expected, the first-order result fails to recover the exchange splitting around the chemical minima, which are at 5*a*_{0} for Li⋯Li, 3.5*a*_{0} for Li⋯N, and 2.1*a*_{0} for the nitrogen dimer. Nonetheless, in all these cases the SF-SAPT exchange splitting is very close to the CASSCF one for the separations corresponding to the van der Waals minima of the highest spin states: 7.9*a*_{0} for the lithium dimer,^{62} 10.2*a*_{0} for Li⋯N [this value was obtained in the present work using the spin-restricted CCSD(T) method with the aug-cc-pVQZ basis set], and 7.2*a*_{0} for the nitrogen dimer.^{63} Clearly, the exchange splitting from first-order SAPT exhibits correct asymptotic behavior. It is also worth noting that for Li⋯Li, unlike Li⋯H, the diagonal and spin-flip exchange terms are of the same order as the apparent and mathematical Coulomb energies are no longer close. Moreover, the splittings, as well as $Eelst(10)$ and $Eexch,diag(10)$, show very little basis set dependence: our Li⋯Li tests at *R* = 7.9*a*_{0} show that each of these three first-order contributions agrees between the (dimer-centered) cc-pVDZ and aug-cc-pVQZ basis sets to below 1%.

In the case of the lithium dimer, we could also compare the exchange splittings to the experimentally derived values given by Côté *et al.*^{64} which confirm the correct behavior of the MRCI method. Since the 1*s* shells were frozen in MRCI calculations, we experienced no problems with a size-inconsistent behavior of this method. Quite large differences between splittings obtained from MRCI and CASSCF/SAPT can possibly be attributed to strong influence of intramonomer dynamical correlation on the interaction effects. For the Li⋯N and N⋯N systems, the differences between MRCI, CASSCF, and first-order SAPT are smaller. On the other hand, it is worth noting that the MRCI method fails to reproduce the exponential decay of the exchange splitting for large interatomic separations for Li⋯N and N⋯N, due to its lack of size-consistency.

It might be quite useful to investigate the quality of the *S*^{2} approximation by inspecting the ratio of $Eexch(10)(S2)$ and $Eexch(10)$ for the high-spin state (the nonapproximated values were computed using the high-spin code of Ref. 26). From Fig. 3 it is clear that while for the lithium dimer the single exchange approximation is quite drastic, even in the van der Waals minimum region, it works very well for the nitrogen dimer and the Li⋯N system. Such behavior can be attributed to a very small ionization potential *I*_{p} of the Li atom (4.2 eV) compared to the N atom (14.5 eV) which directly affects the decay rates of the wavefunctions [at long range, this decay rate is proportional to $exp(\u22122Ipr)$^{65}].

### C. O_{2}⋯O_{2}

Oxygen dimer is one of the best studied van der Waals complexes with relevance to the chemistry of atmosphere.^{28} Hence, this was one of the very first systems for which the exchange interaction was studied. Wormer and van der Avoird calculated the exchange splitting within the Heisenberg model and the *J*_{AB} parameter was obtained with variation-perturbation theory.^{36} Later, the exchange splitting was studied by several *ab initio* multireference methods. In particular, global CASPT2 surfaces for all multiplicities were obtained by Bartolomei and co-workers.^{29} To facilitate comparisons to Ref. 29, our SAPT calculations for this system use the same atomic natural orbital-valence triple zeta (ANO-VTZ) basis set and the same bond length in the oxygen molecules (2.28*a*_{0}).

In Fig. 4 we compare the splitting between the highest and lowest spin states of the O_{2}⋯O_{2} complex for four basic angular geometries: H-shape, linear, T-shape, and X-shape with previous studies of Wormer and van der Avoird^{36} and Bartolomei and co-workers,^{29} the latter including CASSCF, MRCI, and CASPT2 calculations. The overall agreement of the spin-flip SAPT exchange energy with these references is very satisfactory. The exchange splittings for the studied geometries are in a very good agreement with the CASSCF of Bartolomei *et al.* and Wormer and van der Avoird perturbation calculations for the H-, T-, and L-shape complexes. For the X-shape geometry, the first-order exchange splittings perform similarly to MRCI and are significantly smaller compared to Heitler-London calculations of Ref. 36 and to CASSCF. In the CASPT2 calculations, the exchange splittings exhibit a node for the X-shape orientation, which is also the case for the SAPT exchange at about 7.5*a*_{0} but it is not clearly seen in Fig. 4. This node can, however, be seen on the logarithmic plot of the diagonal and spin-flip parts of the exchange splitting, Fig. 5, as a dip in the absolute value of the splitting. Note that the diagonal part of the exchange energy is very large in the oxygen dimer for all configurations including X-shape, almost two orders of magnitude larger than the spin-flip part. Such a big domination of the diagonal part is responsible for the fact that the oxygen molecules strongly repel each other at short range and form only weakly bound complexes for all spin states of the dimer. Finally, let us remark that the single exchange approximation works very well for the oxygen dimer, as shown by the comparison of full and *S*^{2}-approximated $Eexch(10)$ values for the high-spin (quintet) complex displayed in Fig. 6. This behavior is consistent with the observation made by Wormer and van der Avoird^{36} that the single-parameter Heisenberg model recovers the quintet-triplet-singlet splittings accurately.

### D. Mn⋯Mn

The manganese dimer has attracted many studies of its exceptional exchange splitting. Its outermost electronic shell 4*s* is doubly filled and its zero orbital angular momentum (*S* state) originates from a cancellation of half-filled *d*-shell momenta of electrons. For this reason, the open shells in the Mn atom are to a large extent screened by the 4*s*^{2} shell and the resulting spin exchange splitting is very small, on the order of 10 cm^{−1} in the minimum of the potential energy curve, which is 2 orders of magnitude less than the binding energy. The first *ab initio* study of exchange interaction in the Mn_{2} system was initiated by Nesbet^{66} who used the Heisenberg model and estimated the *J*_{AB} parameter to be small (up to 62 K at an interatomic separation of 4.5*a*_{0}). More advanced methods were introduced to this system after significant progress in the multireference methods was made.^{67–72} This system is, however, very challenging for multireference methods: there are 5 electrons in the submerged *d*-shell plus two electrons on the outermost 4*s* shell per atom, which makes the active space needed for dimer calculations quite large. Moreover, the perturbation theory approaches like CASPT2 exhibit problems for Mn⋯Mn related to a presence of intruder states.^{71} Clearly, a development of new methods for systems similar to the manganese dimer is warranted.

The ROHF method produces correct orbitals for the Mn atom (with a correct degeneracy of *d* orbitals) and can be straightforwardly applied to first-order SF-SAPT. In Fig. 7, the spin-flip SAPT exchange term is very small, nearly two orders of magnitude smaller than the diagonal exchange energy. Similar to the oxygen dimer, this causes a very small exchange splitting and very small differences between the potential well depths and equilibrium distances for different multiplets. It is also worthwhile to inspect the quality of the *S*^{2} approximation by comparing the high-spin exchange energies. Near the van der Waals minimum of undecaplet Mn⋯Mn, the single exchange approximation reproduces about 90% of the full exchange. In order to assess the performance of first-order SF-SAPT, we have computed the *J*_{AB} parameter and compared it with previous literature data in panel (c) of Fig. 7. In the literature, the *J*_{AB} parameter is usually given at the minimum separation for the high-spin complex (in the case of the work of Buchachenko *et al.*,^{72} the full curve was provided). Since the spin-flip splitting is obtained here at the ROHF level, it works remarkably well. In particular, our results agree very well with the experimental result of Cheeseman *et al.*^{73} (although we realize this might be somewhat fortuitous). Except for the value by Negodaev *et al.*,^{69} all presented values of *J*_{AB} are smaller than the result derived from SF-SAPT.

### E. Phenalenyl dimer

The dimer of the doublet phenalenyl radical is an example of “pancake bonding,” a strong interaction between *π*-stacked radicals that has gathered significant interest in recent years.^{74} The ground state of the phenalenyl dimer is a multireference singlet that exhibits pancake bonding with a binding energy of 11.5 kcal/mol, while its asymptotically degenerate triplet state exhibits only a van der Waals minimum with a depth of 3.6 kcal/mol.^{31} At the interplanar separation of 3.104 Å corresponding to the pancake bonded minimum depicted in Fig. 8, the reference singlet-triplet splitting, computed using the high-level multireference averaged quadratic coupled cluster (MR-AQCC) approach^{75} in Ref. 31 [with a (2,2) active space and the 6-31G(d) basis set], amounts to 17.2 kcal/mol. It should be noted that an accurate description of pancake bonding remains a challenge to many quantum-chemical approaches, most notably those based on density functional theory.^{76}

In this section, we will examine how well the singlet-triplet splitting in the phenalenyl dimer is recovered by the first-order SF-SAPT approach. In this way, we expect to find out whether (1) the simple first-order treatment of the splitting remains valid for such a strong intermolecular attraction, (2) our AO implementation in psi4 is efficient enough to enable $Eexch(10)$ calculations for large complexes, and (3) the DF approximation is just as accurate as for conventional high-spin SAPT. For this purpose, we select the pancake-bonded minimum geometry of the singlet state as established in Ref. 31, vary the intermolecular separation *R*, and perform MO- and AO-based SF-SAPT calculations in a number of basis sets. We did verify that, in the absence of the DF approximation, the SF-SAPT exchange energies from the MO and AO formalisms are identical.

The singlet-triplet splittings obtained in the aug-cc-pVDZ basis set are presented in Fig. 8. In addition to the first-order SF-SAPT calculations, we perform CASSCF(2,2)/aug-cc-pVDZ computations across the same potential energy curve. At the pancake-bonded minimum, the SAPT/aug-cc-pVDZ corrections $Eelst(10)$, $Eexch,diag(10)$, and $Eexch,1flip(10)$ amount to −19.9, 50.6, and 4.3 kcal/mol, respectively. A basis set increase to aug-cc-pVTZ changes these values to −19.7, 50.2, and 4.3 kcal/mol, respectively, confirming that the *E*^{(10)} corrections converge quickly with the basis set size. Thus, the singlet-triplet splitting from first-order SF-SAPT amounts to 8.6 kcal/mol or about half of the benchmark value. On the other hand, the CASSCF method overestimates the benchmark splitting, giving a value of 23.4 kcal/mol. Figure 8 shows that at larger intermonomer separations the SF-SAPT and CASSCF splittings get closer to each other and both quantities exhibit the same long-range decay. Thus, the underestimated splitting from SF-SAPT results likely from the single-exchange approximation applied in our calculations breaking down for the very short intermolecular distances that are characteristic of pancake bonding.

The largest error from the DF approximation, as computed in a smaller, cc-pVDZ basis set, appears for $Eexch,diag(10)$ and amounts to 0.065 kcal/mol (using the cc-pVDZ/JKFIT auxiliary basis); thus, the DF approximation is fully adequate for the first-order SF-SAPT approach. The cost of AO-based SF-SAPT expressions is dominated by the evaluation of generalized Coulomb and exchange (JK) matrices, of which a total of 8 are required. As such, SF-SAPT takes typically about as long as 4 ROHF iterations, each of them requiring two JK builds. For the phenalenyl dimer on a six-core i7-5930K Intel central processing unit (CPU), the density-fitted evaluation of Eqs. (41) and (42) takes 26 s for the cc-pVDZ basis set (674 functions) and 187 s for the cc-pVTZ basis set (1556 functions), making this method tractable for very large molecules. As the algorithm uses the native JK builders of the psi4 code,^{51} any improvements made to these core psi4 objects will also be extended to the SF-SAPT calculations.

## V. SUMMARY

We derived and implemented the molecular-orbital and atomic-orbital formulas for the first-order SAPT exchange energy for two high-spin open-shell species (described by their ROHF determinants) combined to give an arbitrary spin state of the complex (the previously existing open-shell SAPT approaches^{26,27} were restricted to the high-spin state of the complex except for a few system-specific studies). Within the single exchange approximation, the resulting exchange energies for different asymptotically degenerate multiplets are linear combinations of two common terms: the diagonal exchange (common to all spin states) and the spin-flip term (responsible for the multiplet splittings). Thus, the single exchange approximation (which makes double- and higher-spin-flip terms vanish identically) is equivalent to the Heisenberg Hamiltonian model where all splittings within an asymptotically degenerate set of multiplets are expressible by a single parameter *J*_{AB}. This equivalence was proven long ago by Matsen *et al.*;^{36,49} however, this work is the first one to give an explicit, general SAPT expression for the splitting parameter.

We investigated the behavior of the diagonal and spin-flip components of the first-order SAPT exchange correction on a number of interatomic and intermolecular complexes. For the Li⋯H system, we compared the exchange energies with existing FCI-based SAPT calculations^{38} and found a very good agreement. In particular, ROHF-based first-order SF-SAPT reproduces 71% of the FCI exchange splitting in the van der Waals minimum. For several other diatomic complexes, Li⋯Li, Li⋯N, and N⋯N, we compared the difference between $Eexch(10)$ for the highest and lowest spin state to the splittings obtained with the CASSCF and MRCI methods. The $Eexch(10)$ results agree very well with CASSCF for a wide range of distances, from an infinite separation to roughly half the distance between the high- and low-spin minima. The power of the perturbative approach is particularly impressive in the asymptotic region: the $Eexch(10)$ splitting very accurately reproduces the CASSCF values and, unlike MRCI, ensures proper exponential decay.

For the oxygen dimer in four characteristic configurations, we compared the exchange splittings with literature data. Again, the first-order SAPT predictions agree very well with the MRCI and CASSCF results of Bartolomei *et al.*^{29} A particularly interesting test of our new theory was the manganese dimer, which poses a great challenge for supermolecular calculations with a variety of multireference methods. The Heisenberg exchange coupling constant *J*_{AB} derived from $Eexch(10)$ compares qualitatively well with literature data and predicts small spin splitting due to screening from the outermost doubly occupied (4*s*^{2}) shell. However, it should be noted that for such a challenging system the overall agreement between existing supermolecular data is far from satisfactory. Interestingly, first-order spin-flip SAPT predicts the value of *J*_{AB} within the experimental bounds. Finally, a highly efficient AO-based implementation (with density fitting) of $Eexch(10)$ allowed us to apply SF-SAPT to a larger system, the phenalenyl dimer. In that case, we obtained a qualitative agreement with reference CASSCF calculations. The agreement improves at larger separations which might be attributed to the importance of terms beyond the *S*^{2} approximation. The new first-order exchange corrections have been implemented into the development version of the psi4 code^{51} as well as the sapt2012 package.^{53} The AO-based, density fitted SF-SAPT calculation is particularly efficient, leading to very fast (and entirely single-reference) qualitatively correct estimates of the strength of multiplet splittings.

Our new development is merely the first step toward extending SAPT to arbitrary spin states of the interacting complex. However, even at the present level of theory, SF-SAPT could be highly useful as a complementary method for transition metal complexes or potential energy surfaces near dissociation. The method introduced here can be also used with Kohn-Sham orbitals, provided they were obtained from a density functional which yields asymptotically correct exchange-correlation potentials.^{77} The ideas presented here can be generalized to the second-order exchange-induction and exchange-dispersion corrections as well as to a multireference, CASSCF-based description of the noninteracting wavefunctions. Moreover, while the single exchange approximation implies that double- and higher-spin-flip terms vanish, it is a different approximation than merely neglecting multiple spin flips. In fact, going beyond the *S*^{2} approach (following the ideas of Refs. 44–46) might be useful even for two interacting doublets at short range (such as in the phenalenyl dimer) where there is only a single active spin on each monomer to be flipped. The work in all of these directions is in progress in our groups.

## ACKNOWLEDGMENTS

We are indebted to Bogumił Jeziorski for our stimulating discussions on the topic of this work. We thank Massimiliano Bartolomei for sharing the oxygen dimer data and Alexiei Buchachenko for his comments regarding the Mn⋯Mn dimer. K.P. and D.G.A.S. were supported by the U.S. National Science Foundation CAREER Award No. CHE-1351978. P.S.Z. is grateful to the Polish National Science Center (NCN), Grant No. 015/19/B/ST4/02707.

### APPENDIX: PROOF OF EQ. (33)

In this appendix, we calculate the coefficient $Z(SA,SB,S)=c12c0SASB$ that appears when inserting Eq. (27) into Eq. (12), proving Eq. (33). Recall that

In the following, we will extend our previous notation of spin raising and lowering operators $\u015c+$ and $\u015c\u2212$, stating explicitly the molecule (A or B) that the operator is acting on.

We start by noting that the action of the spin-flip operator $\u015cA\u2212\u015cB+$ on the initial product state |*S*_{A} *S*_{A} *S*_{B} − *S*_{B}⟩ produces

Since $\u015c2=\u015cA2+\u015cB2+2(\u015cA)z(\u015cB)z+\u015cA\u2212\u015cB++\u015cA+\u015cB\u2212$, one might rewrite the action of the operator $\u015cA\u2212\u015cB+$ as

The last term in the above equation yields zero since the reverse spin-flip operator $\u015cA+\u015cB\u2212$ cannot raise any more spins on A or lower any more spins on B. Now, by projecting Eqs. (A3) and (A4) onto ⟨*S* (*S*_{A} − *S*_{B})| and using the fact that ⟨*S* (*S*_{A} − *S*_{B})|$\u015c2$ = *S*(*S* + 1)⟨*S* (*S*_{A} − *S*_{B})|, one finds that

which, after a simple rearrangement, gives the formula (33) for *Z*(*S*_{A}, *S*_{B}, *S*).