The recently proposed spin-flip symmetry-adapted perturbation theory (SF-SAPT) first-order exchange energy [Patkowski *et al.*, J. Chem. Phys. **148**, 164110 (2018)] enables the standard open-shell SAPT approach to treat arbitrary spin states of the weakly interacting complex. Here, we further extend first-order SF-SAPT beyond the single-exchange approximation to a complete treatment of the exchanges of electrons between monomers. This new form of the exchange correction replaces the single-exchange approximation with a more moderate single-spin-flip approximation. The newly developed expressions are applied to a number of small test systems to elucidate the quality of both approximations. They are also applied to the singlet-triplet splittings in pancake bonded dimers. The accuracy of the single-exchange approximation deteriorates at short intermolecular separations, especially for systems with few electrons and for the high-spin state of the complex. In contrast, the single-spin-flip approximation is exact for interactions involving a doublet molecule and remains highly accurate for any number of unpaired electrons. Because the single-exchange approximation affects the high-spin and low-spin states of pancake bonded complexes evenly, the resulting splitting values are of similar accuracy to those produced by the formally more accurate single-spin-flip approximation.

## I. INTRODUCTION

The interaction of two open-shell molecules in their high-spin states produces a whole bundle of asymptotically degenerate states corresponding to different values of the total spin quantum number *S* for the complex. The splittings between these states arise from the resonance tunneling of electrons between the two subsystems. Thus, when the intermolecular interaction is described in terms of symmetry-adapted perturbation theory (SAPT),^{1} the splittings can be attributed exclusively to the exchange terms in the interaction: the remaining electrostatic, induction, and dispersion terms are the same for the entire asymptotically degenerate bundle. On the other hand, a uniform description of all spin states is challenging for the more conventional supermolecular approach to interaction energies: while the highest-spin state might often be well described by a single-reference treatment, all the remaining, low-spin states are genuinely multireference. Therefore, in computational studies of open-shell complexes, obtaining sufficiently accurate potential energy surfaces (PESs) for low-spin states is nontrivial. For example, for the interaction of two ground-state ($\Sigma g\u22123$) O_{2} molecules, an accurate PES for the high-spin quintet state could be constructed using restricted coupled-cluster theory with singles, doubles, and perturbative triple excitations (RCCSD(T)).^{2} On the other hand, the PESs for the multireference singlet and triplet states of this complex had to be obtained by combining the RCCSD(T) quintet PES with lower-level (complete active space second order perturbation theory, CASPT2, or multireference configuration interaction with single and double excitations, MRCI) estimates of the singlet-quintet and triplet-quintet splittings.^{3}

As long as the noninteracting monomers are amenable to a single-reference description (in this work, via spin-restricted high-spin determinants Ψ_{A} and Ψ_{B}), the evaluation of SAPT corrections does not require constructing a multireference wavefunction for the complex. Instead, the perturbation series is built on top of a zeroth-order function Ψ_{0} = Ψ_{A}Ψ_{B} and the corrections are expressed in terms of single, double, … excitations out of the single Ψ_{0} reference. However, the established open-shell SAPT formulations, based on either spin-restricted^{4} on unrestricted^{5,6} Hartree-Fock (HF) or Kohn-Sham (KS) determinants, suffer from a different limitation: it is assumed that Ψ_{0} is a pure spin state (an eigenfunction of the $\u015c$^{2} operator). This is true only if the asymptotically degenerate bundle reduces to a single spin state, that is, when either one of the interacting molecules is a singlet or the complex is in the high-spin state, *M*_{S} = ±*S*. Thus, the open-shell SAPT approaches of Refs. 4–6 are only applicable to the high-spin state of the complex. When the interaction energy in a low-spin state is needed, these high-spin approaches can provide the electrostatic, induction, and dispersion corrections, but not the exchange corrections. No estimate of spin splittings can be extracted from these theories.

As the first step towards extending SAPT beyond the high-spin state of the complex, we have recently developed^{7} a new first-order exchange correction valid for an arbitrary spin state of two interacting high-spin open-shell molecules, each described by a restricted open-shell HF (ROHF) determinant. The new formalism involves explicitly projecting Ψ_{0} onto the subspace corresponding to the desired value of the spin quantum number *S*. Similarly to nearly all exchange corrections in closed-shell (and high-spin open-shell) SAPT, the evaluation of matrix elements involving the (*N*_{A} + *N*_{B})-electron antisymmetrizer is greatly simplified by the use of the *single-exchange approximation*, also called the *S*^{2} approximation as it retains terms up to second order in intermolecular overlap integrals.^{1} As shown in Ref. 7, the *S*^{2} approximation allows expressing the first-order exchange energy $Eexch(10)$ for an arbitrary spin state as a linear combination of two matrix elements: a *diagonal exchange energy* that quantifies the spin-averaged effect and a *spin-flip term* responsible for splittings between multiplets. The coefficients in this linear combination do not depend on a particular system and arise solely from the angular momentum algebra. The name “spin-flip term” reflects the fact that the matrix elements of this term are computed between Ψ_{A}Ψ_{B} and a function $\Psi A\u2193\Psi B\u2191$ 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 (that is, an intermolecular spin exchange has occurred). Thus, the new formalism, termed spin-flip SAPT (SF-SAPT), bears some similarities to the spin-flip electronic structure theories of Krylov and co-workers:^{8,9} in both approaches, a multireference low-spin state is accessed from a single reference configuration. However, the excitations in the methods of Krylov *et al.* alter the total spin of the system (the quantum number *M*_{S} is changed). On the other hand, the excitations involved in spin-flip SAPT do not change the value of *M*_{S} as a spin raise on one monomer is always accompanied by a spin lowering on the other monomer.

The spin-flip electronic structure formalism has been generalized to multiple spin flips;^{10} however, a simple and often adequate approximation^{11} relies on using the single-spin-flip approach (accessing the second-highest spin state from the high-spin one) to determine the coupling parameter *J*_{AB} within the Heisenberg spin Hamiltonian model. Within this scheme, the knowledge of *J*_{AB} is sufficient to recover the entire bundle of lower-spin states. As demonstrated in Ref. 7, in the case of spin-flip SAPT, the same single-parameter Heisenberg picture is a direct consequence of the single-exchange approximation where the value of *J*_{AB} results from the single-spin-flip matrix elements (the corresponding elements involving two or more spin flips vanish within the single-exchange approximation).

While it would not be fair to expect quantitative accuracy from such a (conceptually and computationally) simple approximation as first-order perturbation theory, the SF-SAPT approach to $Eexch(10)$ has been shown^{7} to provide reasonable, qualitatively correct multiplet splittings for a number of representative complexes. The accuracy of first-order SF-SAPT is generally similar to the (much more involved) supermolecular complete active space self-consistent field (CASSCF) calculation and the method does not break down at large separations as is the case for size-inconsistent approaches such as MRCI. However, the accuracy of SF-SAPT splittings in the region of strong intermolecular overlap is limited by two issues: the lack of second- and higher-order exchange effects and the single-exchange approximation. Addressing the first issue involves deriving and implementing the SF-SAPT generalizations of the second-order exchange corrections $Eexch\u2212ind(20)$ and $Eexch\u2212disp(20)$, which is in progress in our group and will be described in a separate publication. In this work, we focus on addressing the second issue by deriving and implementing an improved expression for $Eexch(10)$ in SF-SAPT, in which the single-exchange approximation has been replaced by a much milder *single-spin-flip (1-flip) approximation*.

The full nonapproximated $Eexch(10)$ expression in standard closed-shell SAPT has been introduced a long time ago^{12} and this expression, unlike those for higher SAPT exchange corrections, is available in most SAPT implementations including the high-spin open-shell ones.^{4–6} Much more recently, Schäffer and Jansen derived and implemented nonapproximated expressions for the second-order exchange corrections $Eexch\u2212ind(20)$^{13} and $Eexch\u2212disp(20)$.^{14} We will directly adopt Schäffer and Jansen’s approach, based on the properties of singly and doubly excited determinants and their cofactors, in our development of the SF-SAPT $Eexch(10)$ correction (in fact, as we will see below, the spin-flip excitation is just another kind of double excitation such as the one in the formulas for $Eexch\u2212disp(20)$). It should be mentioned that the nonapproximated second-order closed-shell SAPT corrections of Refs. 13 and 14 have recently been implemented within the freely available PSI4NUMPY framework.^{15}

It should be stressed that while the single-exchange approximation *implies* the single-spin-flip approximation (the neglect of matrix elements involving double and higher spin flips), the two approximations are not nearly equivalent. For example, in any interaction involving an ROHF monomer in a doublet state, there is only one unpaired spin on this monomer that can be flipped, so the 1-flip approximation is exact. On the other hand, the single-exchange approximation is not exact for any spin states of the monomers, even closed-shell singlets. It should be noted that the interaction energy in the high-spin state of the complex (*S* = *S*_{A} + *S*_{B}) can be obtained in two fully equivalent ways: using the high-spin formalism of Refs. 4–6 (which leads to the *S*, *M*_{S} = *S* configuration) and using the SF-SAPT formalism (which leads to the *S*, *M*_{S} = *S*_{A} − *S*_{B} configuration degenerate with the previous one). As the first-order exchange energy in the high-spin formalism does not need to involve the *S*^{2} approximation, we will be able to verify the importance of both approximations for the high-spin state of the complex (but not for the low-spin states) by comparison with the high-spin method implemented, e.g., in the PSI4 package.^{6,16} In fact, in Ref. 7, the largest errors due to the *S*^{2} approximation in the high-spin complex occurred for the doublet-doublet Li–Li interaction, and this approximation might be far from exact for another doublet-doublet complex tested there, the pancake-bonded phenalenyl (PLY) dimer. Thus, the replacement of the single-exchange approximation by the 1-flip one (which is exact for two interacting doublets) might be expected to significantly improve the short-range splittings in some of the systems studied in Ref. 7 and possibly in pancake bonded complexes in general.^{17}

The structure of the rest of this article is as follows: In Sec. II, we develop the formalism and derive formulas for the arbitrary-spin first-order exchange correction in terms of molecular-orbital (MO) integrals. In Sec. III, we recast the MO formulas into the atomic-orbital (AO) basis to facilitate an efficient implementation that does not require integral transformation and can utilize the benefits of density fitting (DF). Sec. IV contains the results of our new methodology for the systems studied in Ref. 7 as well as several larger pancake-bonded complexes. Finally, Sec. V presents conclusions.

## II. MOLECULAR ORBITAL FORMALISM

Throughout this paper, the indices *i* and *j* denote all occupied spinorbitals of monomers A and B, respectively. Furthermore, the index *i* is split into inactive (*k*, corresponding to a doubly occupied orbital) and active (*m*) spinorbitals of A, and the index *j* is split into inactive (*l*) and active (*n*) spinorbitals of B. The indices *r* and *s* are used for arbitrary spinorbitals occupied in the zeroth-order wavefunction regardless of the monomer. We will add an arrow *↑* or *↓* to the spinorbital index whenever an explicit specification of the spin is necessary. Ψ_{A} and Ψ_{B} are the ground-state wavefunctions for the individual monomers and are assumed to be ROHF wavefunctions where the unpaired electrons in Ψ_{A} have *α* spin and the unpaired electrons of Ψ_{B} have *β* spin. The product of these two wavefunctions is considered the zeroth-order dimer wavefunction. As with the previous derivation of the SF-SAPT first-order exchange correction within the single-exchange approximation, the current correction is computed within the symmetrized Rayleigh-Schrödinger (SRS) formalism.^{1,18} In the case of low-spin states, the (*N*_{A} + *N*_{B})-electron antisymmetrizer $A$ is accompanied by the spin projector $PSMS$, which projects the dimer wavefunction onto the subspace corresponding to the spin quantum numbers *S* and *M*_{S}. As such, the SAPT first-order interaction energy for a desired total spin is obtained from the expression

where *V* is the perturbation operator that collects the interactions between the monomers. The spin projector acting on the dimer wavefunction is approximated by its expansion truncated after a single spin flip

where *c*_{0} and *c*_{1} are the Clebsch-Gordan coefficients ⟨*S*(*S*_{A} − *S*_{B})|*S*_{A}*S*_{A}*S*_{B} (−*S*_{B})⟩ and ⟨*S*(*S*_{A} − *S*_{B})|*S*_{A}(*S*_{A} − 1)*S*_{B}(−*S*_{B} + 1)⟩, respectively. The arrow superscripts denote a spin-flipped monomer wavefunction, defined as $\Psi X\u2193=(1/2SX)\u015c\u2212\Psi X$ or $\Psi X\u2191=(1/2SX)\u015c+\Psi X$. $\u015c$_{±} are spin-raising and spin-lowering operators, which act on a wavefunction by applying the one-electron spin-lowering or raising operators to all electrons in the wavefunction. The result of this operation is the sum of the wavefunctions where one of the active electrons (and only the active electrons) has had its spin flipped. The terms in parentheses in the above definitions normalize the functions $\Psi X\u2193$ and $\Psi X\u2191$ as required in Eq. (2). Due to the assumed spins of the active electrons in Ψ_{A} and Ψ_{B}, the spin lowering only makes sense for Ψ_{A} and spin raising only makes sense for Ψ_{B}. While the application of the full spin projector produces terms with multiple spin flips,^{7} the truncation in Eq. (2) to terms with no more than singly spin-flipped monomers is the essence of the single-spin-flip approximation.

Thus, in order to compute the first-order exchange energy $Eexch(10)$ within the single-spin-flip approximation, one needs to evaluate the four matrix elements present in Eq. (3) and subtract the electrostatic contribution:

The leading terms in the numerator and denominator of Eq. (3) are the previously derived components of the complete SAPT first-order interaction energy^{12} and have the following forms:

where $S$ is the determinant of the overlap matrix of occupied spinorbitals of both monomers, **S**

In Eq. (6), $Sir$ is a first cofactor of the determinant $S$, $Sij,rs$ is a second cofactor, *W*_{AB} is the nuclear repulsion between monomers A and B, *B*_{ir} and *A*_{jr} are elements of the nuclear attraction matrices for the corresponding monomer, and ⟨*ij*||*rs*⟩ are antisymmetrized two-electron integrals in the physicists’ notation. The first cofactor $Sir$ of a determinant is obtained by the deletion of row *i* and column *r* from the original determinant and multiplying the resulting determinant by (−1)^{i+r}. The second cofactor $Sij,rs$ is obtained by the deletion of two rows *i*, *j* and two columns *r*, *s* from the original determinant and multiplying the resulting determinant by (−1)^{i+j+r+s}. Furthermore, it is antisymmetric with respect to the order of deletions: $Sji,rs=\u2212Sij,rs$ and $Sij,sr=\u2212Sij,rs$.

The Cramer’s rule for the relationship between a determinant, its inverse, and its cofactors implies that

where *D*_{ri} are the elements of the inverse of the overlap matrix, **D** = **S**^{−1}. It is instructive to examine the structure of matrices **S** and **D** in more detail. Obviously, **S**, and thus also **D**, is block-diagonal with respect to spin. Now, the spin-up and spin-down blocks of **S** are not the same: the spin-up block contains overlap between orbital types *k*, *m*, and *l*, and the spin-down block contains overlap between orbital types *k*, *l*, and *n*. Thus, the spin-up and spin-down blocks of the inverse matrix **D** are completely distinct even for the common indices such as *k*: $Dkk\u2032\alpha \u2260Dkk\u2032\beta $. Therefore, in the final orbital formulas for the SF-SAPT expressions, we will explicitly specify the spin block of matrix **D** as $Drs\alpha $ or $Drs\beta $. Finally, note that the matrix **D**, as an inverse of a symmetric matrix **S**, is also symmetric.

According to Eq. (8) and the relationship between the first and second cofactors,^{19} the second cofactors are equivalent to

Note that the second cofactors are indeed antisymmetric with regard to swapping either *i* and *j* or *r* and *s*. Using Eqs. (8) and (9), Eq. (6) can be rewritten in the form making explicit use of the *D*_{ri} matrix elements

where we have used the fact that *r* and *s* span the same spinorbital space and are therefore interchangeable.

The other terms in the numerator and denominator of Eq. (3) contain the spin-flipped wavefunctions where one active electron from each monomer has had its spin flipped. These terms are a specific subset of the double excitations needed for the computation of the nonapproximated second-order SAPT exchange dispersion correction which was recently derived by Schäffer and Jansen.^{14} Adopting the approach of Ref. 14, the second terms in the numerator and denominator of Eq. (3) are rewritten as

and

where $Sm\u2191\u2192m\u2193,n\u2193\u2192n\u2191$ is the determinant of the overlap matrix that results from flipping the spin of *m* and *n* in the ket, and $Sm\u2191\u2192m\u2193,n\u2193\u2192n\u2191ir$ and $Sm\u2191\u2192m\u2193,n\u2193\u2192n\u2191ij,rs$ are the first and second cofactors of that determinant. The arrows on the indices denote spin up or spin down. The tilded indices $r\u0303,s\u0303$ in Eq. (12) denote the *contents* of the columns *r*, *s* in the spin-flipped determinant $Sm\u2191\u2192m\u2193,n\u2193\u2192n\u2191$: thus, when *r* = *m*^{↑}, $r\u0303=m\u2193$, and when *r* = *n*^{↓}, $r\u0303=n\u2191$ (otherwise, $r\u0303=r$). Note that this meaning of a tilde over an index is completely different from the notation of Refs. 13 and 14. The four consecutive terms *I*_{1}, …, *I*_{4} in this equation will be analyzed separately—see below. It is obvious at this point that $NA!NB!N!$ appears in all terms in the top and bottom of Eq. (3) and will be canceled out in the total interaction energy.

It is beneficial at this time to define the relationships between non-excited, singly excited, and doubly excited determinants and their cofactors. The relationship between a determinant $S$ and a singly excited determinant $Si\u2192a$ is such that

where *S*_{ra} are elements of the overlap matrix. The first and second cofactors of a singly excited determinant can be expressed in terms of singly excited determinants and the first and second cofactors of the non-excited determinant in the following ways:^{13}

The doubly excited determinants can be written in terms of singly excited determinants as

The first and second cofactors of doubly excited determinants can be expressed in terms of the other components as^{14}

With the relationships between various determinants and cofactors now defined, it is possible to express the remaining parts of Eqs. (11) and (12) in terms of the elements of the **D** and **S** matrices. Starting with the term from the denominator,

where the cancellation of the $Sm\u2191\u2192m\u2193Sn\u2193\u2192n\u2191$ term is due to the spin-diagonal nature of **S** and **D** [cf. Eq. (13)]. A singly excited determinant with a spin-flipping excitation requires coupling **S** and **D** matrix elements of opposite spin, resulting in a zero spin integral. This property will be used repeatedly to simplify the following equations. The new summations in the last line have been truncated from *r* since *m*^{↓} is orthogonal to all occupied spinorbitals on A and *n*^{↑} is orthogonal to all occupied spinorbitals on B (note that, e.g., the index *j* covers both the inactive spinorbitals *l*^{↑}, *l*^{↓} and active spinorbitals *n*^{↓}). The complete denominator is thus equal to

For the numerator term, we will analyze the consecutive contributions *I*_{n} to Eq. (12). The first of these terms, *I*_{1}, is simply the nuclear repulsion term multiplied by the same doubly excited determinant that appears in the denominator. It is easy to see that the terms in the numerator that contain the nuclear repulsion are, in fact, equal to the denominator multiplied by the nuclear repulsion. As such, the total interaction energy contains this term exactly once, which is to be expected.

*I*_{2}, the term containing the nuclear potential of monomer B, contains a first cofactor of a doubly excited determinant. As the summation over *r* contains both *m*^{↑} and *n*^{↓}, it has to be broken into three parts accounting for the special cases shown in Eq. (17). The resulting equation is

We expand the first summation in this term as

and the analogous result for the second summation is

The third summation in *I*_{2} contains the same doubly excited determinant that was previously described, as well as two others. Following the same logic as above to eliminate the singly excited determinants that vanish by spin integration and expand the remaining determinants, and noting that $r\u0303=r$ for all terms in the restricted summation, this last term is equal to

The restriction in the summation in Eq. (24) to *r* not equal to either *m*^{↑} or *n*^{↓} can be lifted since the result of the additional *r* = *m*^{↑} and *r* = *n*^{↓} terms is equal to zero in these cases. For example, when *r* = *m*^{↑}, the first two terms in parentheses cancel each other and the third term is zero due to spin. A similar result can be obtained for *r* = *n*^{↓}.

With these results in hand, the complete *I*_{2} term is

An analogous derivation for *I*_{3} gives

Before tackling these summations, we can take advantage of antisymmetry relations of integrals and cofactors to reduce the number of terms that need to be expanded in the first line of Eq. (27)

Therefore, we need to expand only one of these summations

where the two terms that are omitted in the second line are zero due to spin.

When dealing with the summations in the second and third line of Eq. (27), it is again possible to equate some terms to each other using their antisymmetry. Additionally, *r* and *s* can be swapped arbitrarily since they both span the complete occupied space. Therefore,

and

Now, expanding the term in Eq. (30)

with the terms that are removed vanishing due to spin integration. The analogous treatment of the other term containing single excitations yields

The restrictions in Eq. (27) of *r* ≠ *m*^{↑} and *r* ≠ *n*^{↓} for the summations involving Eqs. (32) and (33), respectively, can be lifted since it can now be seen that these conditions reduce the given term to zero.

For the last summation of *I*_{4}, we begin by noting that $r\u0303=r,s\u0303=s$ under the restrictions of this summation. Next, we expand the doubly excited determinants. The first of them has been handled above and the others expand as follows:

Now we make the necessary replacements and expand the integrals and second cofactors for the last summation in Eq. (27)

Inspection of the terms in Eq. (35) shows that four of them vanish upon spin integration. Moreover, upon the summation over *r*, *s*, the remaining 24 terms can be collected into 12 equal pairs when we again take advantage of the fact that *r* and *s* can be swapped. Again, the restrictions (*r*, *s*) ≠ (*m*^{↑}, *n*^{↓}) on the summation can be lifted because the sum of the terms reduces to zero when either *r* or *s* is equal to *m*^{↑} or *n*^{↓}. The total *I*_{4} term can now be rewritten as

where the common $S$ has been factored out and the $12$ canceled by the factor of 2 from the term pairings.

With all necessary terms in hand, there remain a few minor points to consider. First, it can now be seen that $S$ appears in all terms in the numerator and denominator of the total energy equation and will cancel out. Second, we introduce the following function:^{7}

which replaces both the ratio of the Clebsch-Gordan coefficients and normalization factors from the spin-flipped wavefunctions. Through a proper choice of the factor *Z*, the desired spin state is obtained.^{7} Finally, the equations up to this point remain in the spinorbital form and spin integration has only been accounted for as a means of eliminating terms that vanish. We will now make the spin integration explicit and specify the exact spin blocks of the **D** matrix that give nonzero contributions in the resulting orbital expressions. This step is necessary in preparation for the atomic-orbital equivalents of these expressions (which will be derived in Sec. III), as the latter would otherwise lose any information about the spin combinations that result in nonvanishing contributions. In the resulting expressions below, all sums run over the respective orbitals and not spinorbitals. However, the spin integration results in different ranges of summation for orbitals occupied by spin-up and spin-down electrons. These ranges will be specified by an overbar for spin-up indices and an underbar for the spin-down ones: specifically, $i\xaf$, $j\xaf$, and $r\xaf$ represent the occupied orbitals of A, B, or either monomer, respectively, which can be combined with an *α* spin function, and $i\u0332$, $j\u0332$, and $r\u0332$ denote the corresponding orbital types that can be combined with a *β* spin function. In terms of the inactive and active orbitals on both monomers, the summations over $i\xaf$, $j\xaf$, $r\xaf$, $i\u0332$, $j\u0332$, and $r\u0332$ break up into summations over (*k*, *m*), (*l*), (*k*, *l*, *m*), (*k*), (*l*, *n*), and (*k*, *l*, *n*), respectively. In this notation, the orbital equivalents of Eqs. (10), (19), (25), (26), and (36) are

## III. ATOMIC ORBITAL FORMALISM

We now recast the equations derived in Sec. II from molecular orbitals to atomic orbitals, which provides significant computational benefits as integral transformation is avoided and efficient generalized Coulomb and exchange matrix codes can be utilized.^{7,16,20,21} Let us start from the expression for the spin diagonal component of the numerator [Eq. (38)]. In the AO form, the leading nuclear repulsion term stays the same, so we move to the first attractive potential term

where **X** · **Y** = *∑*_{KL}**X**_{KL}**Y**_{KL} is the matrix dot product, the capital letter indices span the AO set, **B** is the nuclear potential matrix of monomer B in the AO basis, *C*_{rK} are the SCF coefficients of spinorbital *r*, and $D\alpha RI$, defined by the above equation, is the representation of $Dri\alpha $ in the AO basis. We use capital letters in the superscript of $D\alpha RI$ to remind the reader that *R*, *I* are no longer the indices of the matrix but merely specify the spinorbital spaces over which the contributions have been summed. The second term involving *B*_{ir} results in a fully analogous contribution $B\u22c5D\beta RI$, and we can now define

as the AO representations of the **D** matrices, unlike the MO ones, are all the same size and can be added. Breaking the indices into the inactive and active ones, corresponding to a further restriction of the summations in the definition of $D\alpha RI$ and $D\beta RI$, yields

Note that the different **D** matrices are in general not symmetric. The derivation of the remaining nuclear attraction terms follows an analogous path, giving

The two-electron terms give way to similar transformations of the **D** matrices. Some representative spin-block contributions are

where **J**[**X**] and **K**[**X**] are the generalized Coulomb and exchange matrices, defined as follows:

The other spin blocks give rise to analogous contributions and the complete spin diagonal term can now be written in the AO form as

Next, we move to the spin-flipped term of the numerator [Eq. (12)]. The *I*_{1} term in Eq. (12) is the nuclear repulsion part for the spin-flip term and is explicitly defined in Eq. (39). This term is transformed into the AO basis as follows:

Note that we used the fact that the orbitals are spin-restricted so that, for example, *C*_{m↑L} = *C*_{m↓L} = *C*_{mL}. The term in the last line of Eq. (52) also constitutes the spin-flip part of the denominator after dividing by *W*_{AB}. Using the same methods as above, the *I*_{2} and *I*_{3} terms of Eq. (12) transform as

and

Finally, *I*_{4} transforms into the following form:

Note that the final formula (55) can be written in many equivalent ways due to the “Hermitian-like” symmetry of the Coulomb and exchange matrices

The form adopted in Eq. (55) minimizes the number of Coulomb and exchange matrix evaluations necessary.

## IV. RESULTS

Both the MO and AO formulas for the newly developed $Eexch(10)$ correction were implemented using PSI4^{16} and the PSI4NUMPY framework.^{15} The agreement between the high-spin SF-SAPT results and the conventional ROHF-based SAPT in PSI4^{6} was verified for the test systems where the single-spin-flip approximation is exact, e.g., the systems containing Li. For other systems, the 1-flip values are compared with the exact high-spin $Eexch(10)$ to determine the adequacy of the proposed approximation.

The aug-cc-pVTZ basis set^{22–24} was used for the Li⋯Li, Li⋯N, N⋯N, and Mn⋯Mn complexes. The results for the pancake bonded systems were obtained in the aug-cc-pVDZ basis set.^{22,23} To allow for comparison with the previous results, the O_{2}⋯O_{2} calculations were performed in the ANO-VTZ basis set^{25,26} and the Li⋯H calculations were performed in the basis set of Ref. 27. When the density fitting (DF) approximation was utilized, the def2-QZVPP/JKFIT^{28} sets were used for Li and Mn atoms as well as for $O2\cdots O2$. All other atoms used the aug-cc-pV*X*Z/JKFIT^{29} sets with *X* being the same as for the orbital basis set. Throughout the discussion below, we will refer to the exchange splitting between the highest and lowest spin states of the complexes, which is defined as

For the Mn⋯Mn complex, difficulties converging the ROHF iterations in PSI4 led us to use MOLPRO^{30} for the monomer calculations. The ROHF orbital energies and vectors were read from the MOLPRO output and passed into PSI4NUMPY. It should also be noted that in the initial SF-SAPT work (Ref. 7), the Mn⋯Mn results were computed from incorrectly converged monomer ROHF wavefunctions. Therefore, the corrected SF-SAPT(*S*^{2}) results for this system will be presented below.

### A. Diatomics and the O_{2} dimer

The first of the smaller test complexes that we consider is Li⋯H, where only one pair of electrons can be exchanged. As already mentioned, the single-spin-flip approximation is exact whenever at least one monomer in the complex is a doublet. The values of $Eexch(10)$ calculated in both approximations for the singlet and triplet states of the Li⋯H complex are provided in Table I, along with the $\Delta Eexch(10)$ value and the full configuration interaction (FCI) based SAPT results of Ref. 27. While the 1-flip results agree remarkably well with the FCI-based values, confirming that the effects of intramolecular correlation are minuscule for this system as observed earlier,^{7} there is a noticeable discrepancy between these values and the *S*^{2} results. While it may come as a surprise, the *S*^{2} approximation in the commonly used form is not exact even in this case. During the derivation of the $Eexch(10)$ correction in the *S*^{2} approximation,^{31} the following intermediate formula is reached:

where $P$ is the single-exchange operator. A related formula is obtained during the derivation of the $Eexch(10)(S2)$ correction in SF-SAPT [Eq. (9) of Ref. 7]. At that point, the $Eexch(10)\u27e8P\u27e9$ term on the lhs of Eq. (58) is normally neglected because it is at least of the order *S*^{4} [*S*^{2} in $Eexch(10)$ and *S*^{2} in $\u27e8P\u27e9$]. However, this term is not zero even for systems such as Li⋯H where only a single electron exchange is possible. We have verified that the differences between the *S*^{2} and nonapproximated results for Li⋯H originate solely from the removal of the $Eexch(10)\u27e8P\u27e9$ term.

. | $Eexch(10)(S=0)$ . | $Eexch(10)(S=1)$ . | $\Delta Eexch(10)$ . | ||||||
---|---|---|---|---|---|---|---|---|---|

R
. | FCI . | S^{2}
. | 1-flip . | FCI . | S^{2}
. | 1-flip . | FCI . | S^{2}
. | 1-flip . |

6.00 | −1434.15 | −1484.93 | −1434.71 | 1562.92 | 1509.18 | 1564.05 | 2997.07 | 2994.11 | 2998.76 |

7.00 | −602.08 | −611.89 | −602.71 | 624.05 | 615.34 | 624.87 | 1226.13 | 1227.23 | 1227.58 |

8.00 | −230.99 | −232.81 | −231.39 | 234.31 | 233.31 | 234.76 | 465.30 | 466.12 | 466.15 |

9.00 | −83.09 | −83.48 | −83.28 | 83.54 | 83.55 | 83.75 | 166.63 | 167.03 | 167.03 |

10.00 | −28.39 | −28.49 | −28.47 | 28.44 | 28.50 | 28.53 | 56.83 | 56.99 | 56.99 |

11.00 | −9.24 | −9.27 | −9.27 | 9.25 | 9.28 | 9.28 | 18.49 | 18.55 | 18.55 |

11.50 | −5.18 | −5.19 | −5.19 | 5.18 | 5.19 | 5.19 | 10.35 | 10.39 | 10.39 |

12.00 | −2.86 | −2.87 | −2.87 | 2.86 | 2.87 | 2.87 | 5.72 | 5.74 | 5.74 |

12.50 | −1.56 | −1.56 | −1.56 | 1.55 | 1.56 | 1.56 | 3.11 | 3.12 | 3.12 |

13.00 | −0.83 | −0.84 | −0.84 | 0.83 | 0.84 | 0.84 | 1.66 | 1.67 | 1.67 |

14.00 | −0.23 | −0.23 | −0.23 | 0.23 | 0.23 | 0.23 | 0.45 | 0.45 | 0.45 |

. | $Eexch(10)(S=0)$ . | $Eexch(10)(S=1)$ . | $\Delta Eexch(10)$ . | ||||||
---|---|---|---|---|---|---|---|---|---|

R
. | FCI . | S^{2}
. | 1-flip . | FCI . | S^{2}
. | 1-flip . | FCI . | S^{2}
. | 1-flip . |

6.00 | −1434.15 | −1484.93 | −1434.71 | 1562.92 | 1509.18 | 1564.05 | 2997.07 | 2994.11 | 2998.76 |

7.00 | −602.08 | −611.89 | −602.71 | 624.05 | 615.34 | 624.87 | 1226.13 | 1227.23 | 1227.58 |

8.00 | −230.99 | −232.81 | −231.39 | 234.31 | 233.31 | 234.76 | 465.30 | 466.12 | 466.15 |

9.00 | −83.09 | −83.48 | −83.28 | 83.54 | 83.55 | 83.75 | 166.63 | 167.03 | 167.03 |

10.00 | −28.39 | −28.49 | −28.47 | 28.44 | 28.50 | 28.53 | 56.83 | 56.99 | 56.99 |

11.00 | −9.24 | −9.27 | −9.27 | 9.25 | 9.28 | 9.28 | 18.49 | 18.55 | 18.55 |

11.50 | −5.18 | −5.19 | −5.19 | 5.18 | 5.19 | 5.19 | 10.35 | 10.39 | 10.39 |

12.00 | −2.86 | −2.87 | −2.87 | 2.86 | 2.87 | 2.87 | 5.72 | 5.74 | 5.74 |

12.50 | −1.56 | −1.56 | −1.56 | 1.55 | 1.56 | 1.56 | 3.11 | 3.12 | 3.12 |

13.00 | −0.83 | −0.84 | −0.84 | 0.83 | 0.84 | 0.84 | 1.66 | 1.67 | 1.67 |

14.00 | −0.23 | −0.23 | −0.23 | 0.23 | 0.23 | 0.23 | 0.45 | 0.45 | 0.45 |

For the Li dimer and the Li⋯N complex, the single-spin-flip $Eexch(10)$ result is exact due to the doublet nature of a lithium atom, meaning there is only one spin to be flipped. Figure 1 shows how severely the $Eexch(10)$ values for the lithium dimer are affected by the *S*^{2} approximation. The deviation of the *S*^{2} results from exact $Eexch(10)$ is already visible at the van der Waals minimum of the triplet state (7.9 bohrs), and is on its way to catastrophic failure at the chemical minimum of the singlet state (5.0 bohrs). At around 4.2 bohrs, the splitting predicted by $Eexch(10)(S2)$ changes sign, leading to an unphysical energetic ordering of the two states. It can also be seen that the *S*^{2} value for the singlet deviates from its complete $Eexch(10)$ counterpart at a slower rate than for the triplet. Similar results are seen for the Li⋯N complex in Fig. 2. The two methods give equivalent results near the 10.2 bohr van der Waals minimum for the quintet state of this system (excluded from the figure) and the *S*^{2} results for the triplet also deviate slowly from the complete exchange results. At the chemical minimum of the triplet state (3.5 bohrs), the *S*^{2} approximation diverges considerably from the complete results, recovering 92% of the exchange energy for the triplet, 88% for the quintet, and 58% of $\Delta Eexch(10)$.

For the Li dimer, we also investigated the significance of the $Eexch(10)\u27e8P\u27e9$ term, which can easily be incorporated into the *S*^{2} approximation by solving Eq. (58) for $Eexch(10)$. Even though the *S*^{2} approximation does not become exact for this system when the $Eexch(10)\u27e8P\u27e9$ term is taken into account (one is still missing the effects of double and triple electron exchanges), the inclusion of this term resulted in a much improved agreement with the complete exchange results. Therefore, the inclusion of $Eexch(10)\u27e8P\u27e9$ while neglecting all other multiple exchanges was tested on a number of other small complexes. Unfortunately, such a treatment was found to drastically overcorrect for the effects missing in $Eexch(10)(S2)$ so that the resulting values were in worse agreement with the complete exchange results than the *S*^{2} data. We suspect that the good performance of the $Eexch(10)\u27e8P\u27e9$ inclusion for Li⋯Li is merely a consequence of each monomer having only one valence electron.

The N⋯N complex is the first of our test systems where the exact and single-spin-flip $Eexch(10)$ values are distinct. Both of these results can be calculated for the high-spin state [the exact one using the standard high-spin SAPT(ROHF) implementation^{4,6}]. This provides an opportunity to gauge the relative accuracy of the 1-flip and *S*^{2} approximations, as shown in Fig. 3. The 1-flip treatment is a much milder approximation that slightly overestimates the exact $Eexch(10)$ value, as opposed to the single-exchange approximation which underestimates it. Figure 3 also shows the low-spin exchange energies, where again the singlet state is less affected by the *S*^{2} approximation than the highest spin state. In this way, the difference in the splittings predicted by the two approximations is primarily driven by the high-spin state error. The range presented in Fig. 3 falls between the septet van der Waals minimum of 7.2 bohrs and the singlet chemical minimum of 2.1 bohrs.

Table II provides the SF-SAPT $Eexch(10)$ values for the manganese dimer. For the undecaplet exchange energies, the single-spin-flip results maintain perfect agreement with the exact high-spin $Eexch(10)$ throughout nearly the entire range presented here, while the *S*^{2} results begin to deviate at a relatively long-range distance of 9 bohrs. That said, the energy splittings from the two approximations are in close agreement due to the more even deviation of the *S*^{2} values from the 1-flip ones for both the highest and lowest spin states of this system.

. | $Eexch(10)(S=0)$ . | $Eexch(10)(S=5)$ . | $\Delta Eexch(10)$ . | ||||
---|---|---|---|---|---|---|---|

R
. | S^{2}
. | 1-flip . | S^{2}
. | 1-flip . | Exact . | S^{2}
. | 1-flip . |

4.50 | 115.32 | 157.35 | 117.51 | 159.52 | 159.51 | 2.19 | 2.17 |

5.00 | 79.10 | 99.44 | 79.89 | 100.23 | 100.23 | 0.79 | 0.79 |

5.50 | 52.96 | 62.39 | 53.25 | 62.69 | 62.69 | 0.29 | 0.30 |

6.00 | 34.62 | 38.84 | 34.73 | 38.95 | 38.95 | 0.11 | 0.11 |

6.50 | 22.16 | 23.98 | 22.20 | 24.02 | 24.02 | 0.04 | 0.04 |

7.00 | 13.92 | 14.68 | 13.94 | 14.70 | 14.70 | 0.02 | 0.02 |

7.50 | 8.60 | 8.91 | 8.61 | 8.92 | 8.92 | 0.01 | 0.01 |

8.00 | 5.25 | 5.37 | 5.25 | 5.37 | 5.37 | 0.00 | 0.00 |

9.00 | 1.88 | 1.90 | 1.88 | 1.90 | 1.90 | 0.00 | 0.00 |

10.00 | 0.65 | 0.65 | 0.65 | 0.65 | 0.65 | 0.00 | 0.00 |

11.00 | 0.22 | 0.22 | 0.22 | 0.22 | 0.22 | 0.00 | 0.00 |

12.00 | 0.07 | 0.07 | 0.07 | 0.07 | 0.07 | 0.00 | 0.00 |

. | $Eexch(10)(S=0)$ . | $Eexch(10)(S=5)$ . | $\Delta Eexch(10)$ . | ||||
---|---|---|---|---|---|---|---|

R
. | S^{2}
. | 1-flip . | S^{2}
. | 1-flip . | Exact . | S^{2}
. | 1-flip . |

4.50 | 115.32 | 157.35 | 117.51 | 159.52 | 159.51 | 2.19 | 2.17 |

5.00 | 79.10 | 99.44 | 79.89 | 100.23 | 100.23 | 0.79 | 0.79 |

5.50 | 52.96 | 62.39 | 53.25 | 62.69 | 62.69 | 0.29 | 0.30 |

6.00 | 34.62 | 38.84 | 34.73 | 38.95 | 38.95 | 0.11 | 0.11 |

6.50 | 22.16 | 23.98 | 22.20 | 24.02 | 24.02 | 0.04 | 0.04 |

7.00 | 13.92 | 14.68 | 13.94 | 14.70 | 14.70 | 0.02 | 0.02 |

7.50 | 8.60 | 8.91 | 8.61 | 8.92 | 8.92 | 0.01 | 0.01 |

8.00 | 5.25 | 5.37 | 5.25 | 5.37 | 5.37 | 0.00 | 0.00 |

9.00 | 1.88 | 1.90 | 1.88 | 1.90 | 1.90 | 0.00 | 0.00 |

10.00 | 0.65 | 0.65 | 0.65 | 0.65 | 0.65 | 0.00 | 0.00 |

11.00 | 0.22 | 0.22 | 0.22 | 0.22 | 0.22 | 0.00 | 0.00 |

12.00 | 0.07 | 0.07 | 0.07 | 0.07 | 0.07 | 0.00 | 0.00 |

The last system that we consider in this section is the first bimolecular complex, the O_{2} dimer. This system is interesting due to the large effect that geometry has on the splitting between the singlet and quintet states. Figure 4 shows the $\Delta Eexch(10)$ values for the four representative geometries at several center-of-mass distances. The spin splitting is much larger in the L (linear) configuration and much smaller in the X configuration. For the H, T, and X structures, the 1-flip and *S*^{2} splittings are virtually indistinguishable, while the two approximations slightly deviate from each other at short range for the L configuration. We report in the supplementary material that the 1-flip results show great agreement with the exact high-spin $Eexch(10)$ values for all four geometries. The very good recovery of the singlet-quintet splitting by the *S*^{2} approximation, shown in Fig. 4, stems from a cancellation of beyond-*S*^{2} effects between the two spin states of the complex.

### B. Pancake bonded systems

An interesting possible application of SF-SAPT, already initiated in Ref. 7, is the spin-state splittings of pancake bonded dimers.^{17} These systems are composed of radicals with highly delocalized singly occupied orbitals which interact in a fashion that is intermediate between covalent and noncovalent bonding. The pancake bonded singlet minimum is separated by a relatively small gap from the van der Waals-bonded triplet state. SF-SAPT is one of the simplest approaches to investigate these splittings, and the single-spin-flip approximation is exact for this doublet-doublet interaction. We have selected a number of pancake-bonded complexes to test our method, illustrating the favorable performance of our AO implementation with density-fitted generalized Coulomb and exchange matrices [Eq. (50)] from PSI4.^{16} Figure 5 shows the monomers selected for our investigation. These include the phenalenyl (PLY) radical,^{32} one of the prototypical examples of pancake bonding systems, as well as four of its derivatives^{33} and the larger trioxotriangulene (TOT) radical.^{34} The orientations of these systems in the homodimer prefer a maximum overlap of the delocalized singly occupied orbitals, i.e., the monomers stack directly on top of one another with their atoms lining up. With the central carbons of the monomers aligned, these dimers can have a staggered or eclipsed conformation as illustrated in Fig. 5. Both conformations of the TOT dimer are considered, but only the staggered conformation is used for PLY and its derivatives.

The spin state splitting in the PLY dimer was analyzed previously within the *S*^{2} formulation of first-order SF-SAPT, together with the (much more demanding) supermolecular complete active space self consistent field theory (CASSCF).^{7} Figure 6 presents these results along with the new nonapproximated $\Delta Eexch(10)$ ones, the multireference averaged quadratic coupled cluster (MR-AQCC)^{35} benchmark from Ref. 32, and the density functional theory (DFT) results obtained with the M05-2X functional.^{36,37} Similar to the manganese dimer case, the splittings produced by the two versions of SF-SAPT do not differ greatly within this range. Some deviation can be observed at the shortest ranges in Fig. 6, but it is not drastic compared to the differences between various methods presented. The effect of the *S*^{2} approximation appears to be relatively more consistent between the two spin states than observed for smaller systems.

The results for the homodimers of the PLY derivatives, shown in Fig. 7, further indicate that the *S*^{2} and 1-flip $\Delta Eexch(10)$ values agree well even at shorter distances. This observation would seem to support the idea that the removal of the *S*^{2} approximation is less important to the improvement of the SF-SAPT splittings than the inclusion of higher-order exchange terms (which is in progress in our group). For comparison, we look at the M05-2X splitting values from Ref. 33. This functional was chosen as the most suitable for pancake bonded systems based on previous benchmarking.^{37} The tri(tert-butyl)phenalenyl (TBPLY) dimer has the smallest splitting value at the minimum, where the SF-SAPT values for both approximations agree with each other and represent about two-thirds of the DFT result. For the other derivatives, the SF-SAPT results are about half the magnitude of the DFT results at the minimum. Also, for the triaminophenalenyl (TAPLY), trifluorophenalenyl (TFPLY), and trimethylphenalenyl (TMPLY) dimers, SF-SAPT shows a different ordering of the splitting values. The DFT results show the order of the splittings as TAPLY > TFPLY > TMPLY, while both SF-SAPT methods provide the order TMPLY > TAPLY > TFPLY for the M05-2X-optimized structures.

Finally, we look at the trioxotriangulene dimer in its staggered and eclipsed conformations in Table III. The SF-SAPT methods agree with the DFT results on the ordering of the splittings of the two conformations, showing a larger splitting for the staggered one. Again, the SF-SAPT values are between half and two-thirds of the M05-2X values.^{34} The *S*^{2} results for both $\Delta Eexch(10)$ and the individual $Eexch(10)$ values are quite accurate in the presented range for this system compared to the nonapproximate 1-flip data.

. | $Eexch(10)(S=0)$ . | $Eexch(10)(S=1)$ . | $\Delta Eexch(10)$ . | |||
---|---|---|---|---|---|---|

R/R_{0}
. | S^{2}
. | 1-flip . | S^{2}
. | 1-flip . | S^{2}
. | 1-flip . |

Eclipsed | ||||||

0.80 | 225.43 | 228.41 | 238.31 | 241.22 | 12.88 | 12.81 |

0.90 | 79.48 | 79.74 | 84.67 | 84.89 | 5.19 | 5.15 |

1.00 | 27.11 | 27.11 | 29.05 | 29.05 | 1.94 | 1.93 |

1.10 | 9.00 | 9.00 | 9.70 | 9.69 | 0.70 | 0.69 |

1.20 | 2.92 | 2.92 | 3.16 | 3.16 | 0.24 | 0.24 |

1.40 | 0.28 | 0.28 | 0.31 | 0.31 | 0.03 | 0.03 |

1.60 | 0.02 | 0.02 | 0.02 | 0.02 | 0.00 | 0.00 |

Staggered | ||||||

0.80 | 417.36 | 426.73 | 443.09 | 453.31 | 25.74 | 26.58 |

0.90 | 167.05 | 168.34 | 179.42 | 180.70 | 12.37 | 12.35 |

1.00 | 64.48 | 64.60 | 69.86 | 69.94 | 5.38 | 5.34 |

1.10 | 24.17 | 24.17 | 26.36 | 26.35 | 2.19 | 2.18 |

1.20 | 8.85 | 8.84 | 9.70 | 9.69 | 0.85 | 0.85 |

1.40 | 1.10 | 1.10 | 1.22 | 1.22 | 0.12 | 0.12 |

1.60 | 0.12 | 0.12 | 0.13 | 0.13 | 0.01 | 0.01 |

. | $Eexch(10)(S=0)$ . | $Eexch(10)(S=1)$ . | $\Delta Eexch(10)$ . | |||
---|---|---|---|---|---|---|

R/R_{0}
. | S^{2}
. | 1-flip . | S^{2}
. | 1-flip . | S^{2}
. | 1-flip . |

Eclipsed | ||||||

0.80 | 225.43 | 228.41 | 238.31 | 241.22 | 12.88 | 12.81 |

0.90 | 79.48 | 79.74 | 84.67 | 84.89 | 5.19 | 5.15 |

1.00 | 27.11 | 27.11 | 29.05 | 29.05 | 1.94 | 1.93 |

1.10 | 9.00 | 9.00 | 9.70 | 9.69 | 0.70 | 0.69 |

1.20 | 2.92 | 2.92 | 3.16 | 3.16 | 0.24 | 0.24 |

1.40 | 0.28 | 0.28 | 0.31 | 0.31 | 0.03 | 0.03 |

1.60 | 0.02 | 0.02 | 0.02 | 0.02 | 0.00 | 0.00 |

Staggered | ||||||

0.80 | 417.36 | 426.73 | 443.09 | 453.31 | 25.74 | 26.58 |

0.90 | 167.05 | 168.34 | 179.42 | 180.70 | 12.37 | 12.35 |

1.00 | 64.48 | 64.60 | 69.86 | 69.94 | 5.38 | 5.34 |

1.10 | 24.17 | 24.17 | 26.36 | 26.35 | 2.19 | 2.18 |

1.20 | 8.85 | 8.84 | 9.70 | 9.69 | 0.85 | 0.85 |

1.40 | 1.10 | 1.10 | 1.22 | 1.22 | 0.12 | 0.12 |

1.60 | 0.12 | 0.12 | 0.13 | 0.13 | 0.01 | 0.01 |

## V. SUMMARY

We derived a new formula for the first-order spin-flip SAPT exchange energy,^{7} valid for an arbitrary spin state of a weakly interacting complex, replacing the conventional single-exchange (*S*^{2}) approximation by the much milder single-spin-flip approximation. The new SF-SAPT correction, with the noninteracting monomers described by their ROHF determinants, involves terms similar to those appearing in the expressions for the complete (non-*S*^{2}) first-order exchange and second-order exchange-dispersion energies.^{12,14} In this way, the spin flips are treated as a subset of the double excitations found in second-order dispersion and exchange-dispersion corrections. The resulting equations were implemented in both their molecular orbital and atomic orbital forms, where the latter allows for the application of this method to much larger systems.

The newly enhanced first-order SF-SAPT approach was applied to the same selection of diatomic and small molecular test systems as in Ref. 7. The *S*^{2} approximation is not exact even for the Li⋯H system where only a single electron exchange is possible, as this approximation neglects a term that is a product of two single exchanges. This approximation is also particularly poor for the lithium dimer, leading to an unphysical crossing of the singlet and triplet curves at short range. In contrast, the single-spin-flip approximation was demonstrated to be much milder than the *S*^{2} one. The 1-flip treatment is formally exact for any dimer with at least one doublet monomer and deviates from the exact $Eexch(10)$ much slower than the *S*^{2} variant in other cases. Another observation made possible by the new development is that the *S*^{2} approximation is generally better for the low-spin states of small complexes than for the high-spin state. As the size of the system increases, at least for the complexes considered here the effect of the *S*^{2} approximation on the high and low-spin states becomes comparable.

The new $Eexch(10)$ formulation was further applied to the determination of the singlet-triplet splittings for the phenalenyl radical dimer and other pancake bonded systems. These calculations are feasible, thanks to the recasting of the MO formulas into their AO form, which allows us to take advantage of PSI4’s efficient tools for producing density-fitted generalized Coulomb and exchange matrices.^{16} The resulting implementation is only somewhat more expensive than the $Eexch(10)(S2)$ one of Ref. 7, exhibiting the same *N*^{4} scaling with the basis set size with a somewhat larger prefactor. In fact, the computation time is dominated by the construction of the Coulomb and exchange matrices, and the number of such matrices is 8 for the *S*^{2} correction and 11 for the 1-flip one. Despite the formal exactness of the 1-flip treatment, it does not provide especially different first-order splittings compared to the *S*^{2} approximation. This is due to a more even effect of the *S*^{2} approximation on the singlet and triplet states of these systems. The SF-SAPT splittings were compared to a selection of literature values, in particular, the M05-2X results from Refs. 33, 34, and 37. In comparison to the DFT results, the first-order SF-SAPT treatment underestimates the splittings and shows a different ordering of the PLY derivatives. While these splittings are generally difficult to calculate and the literature results are far from benchmark quality, we cannot expect a high quantitative accuracy from a simple first-order perturbation theory that also neglects intramolecular electron correlation. The work to extend the new SF-SAPT formalism to arbitrary-spin second-order exchange induction and exchange dispersion energies is in progress in our group.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the numerical data presented in all figures and the Cartesian coordinates for the complexes considered in this work.

## ACKNOWLEDGMENTS

This work was supported by the U.S. National Science Foundation CAREER Award No. CHE-1351978. We thank Bogumił Jeziorski and Piotr Żuchowski for stimulating discussions.