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.

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 (Σg3) O2 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-restricted4 on unrestricted5,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 Ŝ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, MS = ±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 developed7 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 (NA + NB)-electron antisymmetrizer is greatly simplified by the use of the single-exchange approximation, also called the S2 approximation as it retains terms up to second order in intermolecular overlap integrals.1 As shown in Ref. 7, the S2 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 ΨAΨB 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 MS is changed). On the other hand, the excitations involved in spin-flip SAPT do not change the value of MS 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 approximation11 relies on using the single-spin-flip approach (accessing the second-highest spin state from the high-spin one) to determine the coupling parameter JAB within the Heisenberg spin Hamiltonian model. Within this scheme, the knowledge of JAB 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 JAB 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 shown7 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 Eexchind(20) and Eexchdisp(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 ago12 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 Eexchind(20)13 and Eexchdisp(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 Eexchdisp(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 = SA + SB) can be obtained in two fully equivalent ways: using the high-spin formalism of Refs. 4–6 (which leads to the S, MS = S configuration) and using the SF-SAPT formalism (which leads to the S, MS = SASB configuration degenerate with the previous one). As the first-order exchange energy in the high-spin formalism does not need to involve the S2 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 S2 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.

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 (NA + NB)-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 MS. As such, the SAPT first-order interaction energy for a desired total spin is obtained from the expression

Eint(10)=ΨAΨB|VAPSMS|ΨAΨBΨAΨB|APSMS|ΨAΨB,
(1)

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

PSMSΨAΨB=c0ΨAΨB+c1ΨAΨB,
(2)

where c0 and c1 are the Clebsch-Gordan coefficients ⟨S(SASB)|SASASB (−SB)⟩ and ⟨S(SASB)|SA(SA − 1)SB(−SB + 1)⟩, respectively. The arrow superscripts denote a spin-flipped monomer wavefunction, defined as ΨX=(1/2SX)ŜΨX or ΨX=(1/2SX)Ŝ+ΨX. Ŝ± 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 ΨX and ΨX 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.

The combination of Eqs. (1) and (2) produces the following modified interaction energy equation:

Eint(10)=ΨAΨB|VA|ΨAΨB+c1c0ΨAΨB|VA|ΨAΨBΨAΨB|A|ΨAΨB+c1c0ΨAΨB|A|ΨAΨB.
(3)

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:

Eexch(10)=Eint(10)Eelst(10)=Eint(10)ΨAΨB|V|ΨAΨB.
(4)

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

ΨAΨB|A|ΨAΨB=NA!NB!N!S,
(5)
ΨAΨB|VA|ΨAΨB=NA!NB!N!WABS+irBirSir+jrAjrSjr+12ijrsij||rsSij,rs,
(6)

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

S=1  SAB(SAB)T  1.
(7)

In Eq. (6), Sir is a first cofactor of the determinant S, Sij,rs is a second cofactor, WAB is the nuclear repulsion between monomers A and B, Bir and Ajr 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=Sij,rs and Sij,sr=Sij,rs.

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

Dri=1SSir,
(8)

where Dri 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αDkkβ. Therefore, in the final orbital formulas for the SF-SAPT expressions, we will explicitly specify the spin block of matrix D as Drsα or Drsβ. 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

Sij,rs=S(DriDsjDsiDrj).
(9)

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 Dri matrix elements

ΨAΨB|VA|ΨAΨB=NA!NB!N!SWAB+irBirDri+jrAjrDrj+ijrs(ij|rsij|sr)DriDsj,
(10)

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

ΨAΨB|A|ΨAΨB=12SASBmnΨAΨB|A|ΨA,mmΨB,nn=NA!NB!N!12SASBmnSmm,nn
(11)

and

ΨAΨB|VA|ΨAΨB=12SASBmnΨAΨB|VA|ΨA,mmΨB,nn=NA!NB!N!12SASBmnWABSmm,nn+irBir̃Smm,nnir+jrAjr̃Smm,nnjr+12ijrsij||r̃s̃Smm,nnij,rs=NA!NB!N!12SASBmnI1+I2+I3+I4,
(12)

where Smm,nn is the determinant of the overlap matrix that results from flipping the spin of m and n in the ket, and Smm,nnir and Smm,nnij,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̃,s̃ in Eq. (12) denote the contents of the columns r, s in the spin-flipped determinant Smm,nn: thus, when r = m, r̃=m, and when r = n, r̃=n (otherwise, r̃=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 I1, …, I4 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 Sia is such that

Sia=rSraSri=SrDirSra,
(13)

where Sra 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 

Siars=Sri  s=i1SSrsSiaSriSsa  si,
(14)
Siars,tu=Srs,iu  t=iSrs,ti  u=i1SSrs,tuSiaSrs,iuStaSrs,tiSua  i{t,u}.
(15)

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

Sia,jb=1S(SiaSjbSibSja).
(16)

The first and second cofactors of doubly excited determinants can be expressed in terms of the other components as14 

Sia,jbrs=Sjbri  s=iSiarj  s=j1SSia,jbSrsSsa,jbSriSia,sbSrj  s{i,j},
(17)
Sia,jbrs,tu=Srs,ij  t=i,u=jSrs,ji  t=j,u=iSiars,tj  ti,u=jSjbrs,ti  tj,u=iSjbrs,iu  t=i,ujSiars,ju  t=j,ui1SSia,jbSrs,tuSta,jbSrs,iuSua,jbSrs,ti  Sia,tbSrs,juSia,ubSrs,tj+Sta,ubSrs,ij  i,j{t,u}.
(18)

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,

ΨAΨB|A|ΨAΨB  =NA!NB!N!12SASBmnSmm,nn  =NA!NB!N!12SASBmn1SSmmSnnSnmSmn  =NA!NB!N!12SASBSmnjDnjSjmiDmiSin,
(19)

where the cancellation of the SmmSnn 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

NA!NB!N!S1c12c0SASBmnij(DnjSjmDmiSin).
(20)

For the numerator term, we will analyze the consecutive contributions In to Eq. (12). The first of these terms, I1, 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.

I2, 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

I2=iBimSnnim+iBinSmmin+ir,r(m,n)Bir̃1S×Smm,nnSirSrm,nnSimSmm,rnSin.
(21)

We expand the first summation in this term as

iBimSnnim=iBim1SSimSnnSinSmn=iiSBimDniDmiSin,
(22)

and the analogous result for the second summation is

iBinSmmin=ijSBinDmiDnjSjm.
(23)

The third summation in I2 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̃=r for all terms in the restricted summation, this last term is equal to

ir,r(m,n)SBirDriijDmiSinDnjSjm+DmiijDriSinDnjSjm+DniijDmiSinDrjSjm.
(24)

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 I2 term is

I2=SiiBimDniDmiSin+ijBinDmiDnjSjm+irBirDriijDmiSinDnjSjmDmiijDriSinDnjSjmDniijDmiSinDrjSjm.
(25)

An analogous derivation for I3 gives

I3=SijAjmDnjDmiSin+jjAjnDmjDnjSjm+jrAjrDrjijDmiSinDnjSjmDmjijDriSinDnjSjmDnjijDmiSinDrjSjm.
(26)

The breakdown of the most complicated I4 term in Eq. (12) involves each of the cases in Eq. (18)

I4=12ij,r=m,s=nij||mnSij,mn+ij,r=n,s=mij||nmSij,nm+ij,rm,s=nij||rnSmmij,rn+ij,rn,s=mij||rmSnnij,rm+ij,r=m,snij||msSnnij,ms+ij,r=n,smij||nsSmmij,ns+ijrs,(r,s)(m,n)ij||r̃s̃1SSmm,nnSij,rsSrm,nnSij,msSsm,nnSij,rmSmm,rnSij,nsSmm,snSij,rn+Srm,snSij,mn.
(27)

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)

ij||nmSij,nm=ij||mnSij,nm=ij||mnSij,mn.
(28)

Therefore, we need to expand only one of these summations

ij||mnSij,mn=ij|mnij|nmSDmiDnjDniDmj=ij|nmSDmiDnjij|mnSDniDmj,
(29)

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,

rmij||rnSmmij,rn=smij||snSmmij,sn=smij||nsSmmij,ns
(30)

and

rnij||rmSnnij,rm=snij||smSnnij,sm=snij||msSnnij,ms.
(31)

Now, expanding the term in Eq. (30)

ij||rnSmmij,rn=ij||rnSSmmSij,rnSrmSij,mnSnmSij,rm=ij|rnSrmDmiDnj+ij|nrSrmDmiDnj+ij|rnSrmDniDmjij|nrSrmDniDmjij|rnSnmDriDmj+ij|nrSnmDriDmj+ij|rnSnmDmiDrjij|nrSnmDmiDrj=Sij|nrDmiDnjjDrjSjm+ij|rnDniDmjjDrjSjmij|rnDriDmjjDnjSjm+ij|nrDriDmjjDnjSjm+ij|rnDmiDrjjDnjSjmij|nrDmiDrjjDnjSjm
(32)

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

ij||rmSnnij,rm=Sij|mrDniDmjiDriSin+ij|rmDmiDnjiDriSinij|rmDriDnjiDmiSin+ij|mrDriDnjiDmiSin+ij|rmDniDrjiDmiSinij|mrDniDrjiDmiSin.
(33)

The restrictions in Eq. (27) of rm and rn 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 I4, we begin by noting that r̃=r,s̃=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:

Srm,nn=1S(SrmSnnSrnSnm)=Sij(DriSinDnjSjm),Ssm,nn=1S(SsmSnnSsnSnm)=Sij(DsiSinDnjSjm),Smm,rn=1S(SmmSrnSmnSrm)=Sij(DmiSinDrjSjm),Smm,sn=1S(SmmSsnSmnSsm)=Sij(DmiSinDsjSjm),Srm,sn=1S(SrmSsnSrnSsm)=Sij(DrjSjmDsiSin)ij(DriSinDsjSjm).
(34)

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

Sij|rsDriDsjij(DmiSinDnjSjm)+ij|rsDsiDrjij(DmiSinDnjSjm)+ij|srDriDsjij(DmiSinDnjSjm)  ij|srDsiDrjij(DmiSinDnjSjm)+ij|rsDmiDsjij(DriSinDnjSjm)ij|rsDsiDmjij(DriSinDnjSjm)  ij|srDmiDsjij(DriSinDnjSjm)+ij|srDsiDmjij(DriSinDnjSjm)+ij|rsDriDmjij(DsiSinDnjSjm)  ij|rsDmiDrjij(DsiSinDnjSjm)ij|srDriDmjij(DsiSinDnjSjm)+ij|srDmiDrjij(DsiSinDnjSjm)  +ij|rsDniDsjij(DmiSinDrjSjm)ij|rsDsiDnjij(DmiSinDrjSjm)ij|srDniDsjij(DmiSinDrjSjm)  +ij|srDsiDnjij(DmiSinDrjSjm)+ij|rsDriDnjij(DmiSinDsjSjm)ij|rsDniDrjij(DmiSinDsjSjm)  ij|srDriDnjij(DmiSinDsjSjm)+ij|srDniDrjij(DmiSinDsjSjm)+ij|rsDmiDnjij(DrjSjmDsiSin)  ij|rsDniDmjij(DrjSjmDsiSin)ij|srDmiDnjij(DrjSjmDsiSin)+ij|srDniDmjij(DrjSjmDsiSin)  ij|rsDmiDnjij(DriSinDsjSjm)+ij|rsDniDmjij(DriSinDsjSjm)+ij|srDmiDnjij(DriSinDsjSjm)  ij|srDniDmjij(DriSinDsjSjm).
(35)

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 I4 term can now be rewritten as

I4=Sijij|nmDmiDnjij|mnDniDmj+ijrij|nrDmiDnjj(DrjSjm)+ij|rnDniDmjj(DrjSjm)ij|rnDriDmjj(DnjSjm)+ij|nrDriDmjj(DnjSjm)+ij|rnDmiDrjj(DnjSjm)ij|nrDmiDrjj(DnjSjm)+ij|mrDniDmji(DriSin)+ij|rmDmiDnji(DriSin)ij|rmDriDnji(DmiSin)+ij|mrDriDnji(DmiSin)+ij|rmDniDrji(DmiSin)ij|mrDniDrji(DmiSin)+ijrsij|rsDriDsjij(DmiSinDnjSjm)+ij|rsDsiDrjij(DmiSinDnjSjm)+ij|rsDriDmjij(DsiSinDnjSjm)ij|rsDmiDrjij(DsiSinDnjSjm)ij|srDriDmjij(DsiSinDnjSjm)+ij|srDmiDrjij(DsiSinDnjSjm)+ij|rsDriDnjij(DmiSinDsjSjm)ij|rsDniDrjij(DmiSinDsjSjm)ij|srDriDnjij(DmiSinDsjSjm)+ij|srDniDrjij(DmiSinDsjSjm)ij|rsDniDmjij(DrjSjmDsiSin)ij|srDmiDnjij(DrjSjmDsiSin),
(36)

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 

Z(SA,SB,S)=c12c0SASB=S(S+1)+2SASBSA(SA+1)SB(SB+1)4SASB,
(37)

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¯, j¯, and r¯ represent the occupied orbitals of A, B, or either monomer, respectively, which can be combined with an α spin function, and i̲, j̲, and r̲ 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¯, j¯, r¯, i̲, j̲, and r̲ 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

ΨAΨB|VA|ΨAΨB=NA!NB!N!SWAB+i¯r¯BirDriα+i̲r̲BirDriβ+j¯r¯AjrDrjα+j̲r̲AjrDrjβ+i¯j¯r¯s¯ij|rsDriαDsjα+i¯j̲r¯s̲ij|rsDriαDsjβ+i̲j¯r̲s¯ij|rsDriβDsjα+i̲j̲r̲s̲ij|rsDriβDsjβi¯j¯r¯s¯ij|srDriαDsjαi̲j̲r̲s̲ij|srDriβDsjβ,
(38)
I1=SWABj̲DnjβSjmi¯DmiαSin,
(39)
I2=Si̲i¯BimDniβDmiαSin+i¯j̲BinDmiαDnjβSjm+i¯r¯BirDriαi¯j̲DmiαSinDnjβSjm+i̲r̲BirDriβi¯j̲DmiαSinDnjβSjmi¯r¯BirDmiαi¯j̲DriαSinDnjβSjmi̲r̲BirDniβi¯j̲DmiαSinDrjβSjm,
(40)
I3=Si¯j̲AjmDnjβDmiαSin+j¯j̲AjnDmjαDnjβSjm+j¯r¯AjrDrjαi¯j̲DmiαSinDnjβSjm+j̲r̲AjrDrjβi¯j̲DmiαSinDnjβSjmj¯r¯AjrDmjαi¯j̲DriαSinDnjβSjmj̲r̲AjrDnjβi¯j̲DmiαSinDrjβSjm,
(41)
I4=Si¯j̲ij|nmDmiαDnjβi̲j¯ij|mnDniβDmjα+i¯j̲r̲j̲ij|nrDmiαDnjβDrjβSjm+i̲j¯r̲j̲ij|rnDniβDmjαDrjβSjmi¯j¯r¯j̲ij|rnDriαDmjαDnjβSjmi̲j¯r̲j̲ij|rnDriβDmjαDnjβSjm+i¯j¯r¯j̲ij|nrDriαDmjαDnjβSjm+i¯j¯r¯j̲ij|rnDmiαDrjαDnjβSjmi¯j¯r¯j̲ij|nrDmiαDrjαDnjβSjmi¯j̲r̲j̲ij|nrDmiαDrjβDnjβSjm+i̲j¯r¯i¯ij|mrDniβDmjαDriαSin+i¯j̲r¯i¯ij|rmDmiαDnjβDriαSini¯j̲r¯i¯ij|rmDriαDnjβDmiαSini̲j̲r̲i¯ij|rmDriβDnjβDmiαSin+i̲j̲r̲i¯ij|mrDriβDnjβDmiαSin+i̲j̲r̲i¯ij|rmDniβDrjβDmiαSini̲j¯r¯i¯ij|mrDniβDrjαDmiαSini̲j̲r̲i¯ij|mrDniβDrjβDmiαSini¯j¯r¯s¯i¯j̲ij|rsDriαDsjαDmiαSinDnjβSjmi¯j̲r¯s̲i¯j̲ij|rsDriαDsjβDmiαSinDnjβSjmi̲j¯r̲s¯i¯j̲ij|rsDriβDsjαDmiαSinDnjβSjmi̲j̲r̲s̲i¯j̲ij|rsDriβDsjβDmiαSinDnjβSjm+i¯j¯r¯s¯i¯j̲ij|rsDsiαDrjαDmiαSinDnjβSjm+i̲j̲r̲s̲i¯j̲ij|rsDsiβDrjβDmiαSinDnjβSjm+i¯j¯r¯s¯i¯j̲ij|rsDriαDmjαDsiαSinDnjβSjm+i̲j¯r̲s¯i¯j̲ij|rsDriβDmjαDsiαSinDnjβSjmi¯j¯r¯s¯i¯j̲ij|rsDmiαDrjαDsiαSinDnjβSjmi¯j¯r¯s¯i¯j̲ij|srDriαDmjαDsiαSinDnjβSjm+i¯j¯r¯s¯i¯j̲ij|srDmiαDrjαDsiαSinDnjβSjm+i¯j̲r̲s¯i¯j̲ij|srDmiαDrjβDsiαSinDnjβSjm+i¯j̲r¯s̲i¯j̲ij|rsDriαDnjβDmiαSinDsjβSjm+i̲j̲r̲s̲i¯j̲ij|rsDriβDnjβDmiαSinDsjβSjmi̲j̲r̲s̲i¯j̲ij|rsDniβDrjβDmiαSinDsjβSjmi̲j̲r̲s̲i¯j̲ij|srDriβDnjβDmiαSinDsjβSjm+i̲j¯r¯s̲i¯j̲ij|srDniβDrjαDmiαSinDsjβSjm+i̲j̲r̲s̲i¯j̲ij|srDniβDrjβDmiαSinDsjβSjmi̲j¯r̲s¯i¯j̲ij|rsDniβDmjαDrjβSjmDsiαSini¯j̲r̲s¯i¯j̲ij|srDmiαDnjβDrjβSjmDsiαSin.
(42)

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

i¯r¯BirDriα=i¯r¯KLCiKBKLCrLDriα=i¯r¯KLBLKCrLDriαCiK=KLBLK(DαRI)LK=BDαRI,
(43)

where X · Y = KLXKLYKL 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, CrK are the SCF coefficients of spinorbital r, and DαRI, defined by the above equation, is the representation of Driα in the AO basis. We use capital letters in the superscript of Dα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 Bir results in a fully analogous contribution BDβRI, and we can now define

DαβRI=DαRI+DβRI,
(44)

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αRI and DβRI, yields

DαRI=DαKK+DαLK+DαMK+DαKM+DαLM+DαMM,
(45)
DβRI=DβKK+DβLK+DβNK.
(46)

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

j¯r¯AjrDrjα+j̲r̲AjrDrjβ=ADαβRJ.
(47)

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

i¯j¯r¯s¯ij|rsDriαDsjα=KLMN(DαRI)KLKM|LN(DαRJ)MN=DαRIJ[DαRJ],
(48)
i¯j¯r¯s¯ij|srDriαDsjα=KLMN(DαRI)KLKN|ML(DαRJ)MN=DαRIK[DαRJ],
(49)

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

J[X]KL=MNKM|LNXMN        K[X]KL=MNKN|MLXMN.
(50)

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

ΨAΨB|VA|ΨAΨB=NA!NB!N!SWAB+BDαβRI+ADαβRJ+DαβRIJ[DαβRJ]DαRIK[DαRJ]DβRIK[DβRJ].
(51)

Next, we move to the spin-flipped term of the numerator [Eq. (12)]. The I1 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:

I1=WABSi¯j̲KLMNDnjβCjKSKLAOCmLDmiαCiMSMNAOCnN=WABSi¯j̲KLMN(CnNDnjβCjKSKLAO)(CmLDmiαCiMSMNAO)=WABSLN(DβNJSAO)NL(DαMISAO)LN=WABS[(DβNJSAO)(DαMISAO)T].
(52)

Note that we used the fact that the orbitals are spin-restricted so that, for example, Cm↑L = Cm↓L = CmL. The term in the last line of Eq. (52) also constitutes the spin-flip part of the denominator after dividing by WAB. Using the same methods as above, the I2 and I3 terms of Eq. (12) transform as

I2=SB(DαMISAODβNI)+B(DβNJSAODαMI)+(BDαβRI)(DαMISAODβNJ)SAOB(DαRISAODβNJSAODαMI)B(DβRJSAODαMISAODβNI)
(53)

and

I3=SA(DαMISAODβNJ)+A(DβNJSAODαMJ)+(ADαβRJ)(DαMISAODβNJ)SAOA(DαRISAODβNJSAODαMJ)A(DβRJSAODαMISAODβNJ).
(54)

Finally, I4 transforms into the following form:

I4=SDαMIK[DβNJ]TDβNIK[DαMJ]T+K[DβNJ]T(DβRJSAODαMI)+K[DβNI]T(DβRJSAODαMJ)J[DαβRI](DβNJSAODαMJ)+K[DαRI]T(DβNJSAODαMJ)+K[DαRJ]T(DβNJSAODαMI)J[DαβRJ](DβNJSAODαMI)+K[DαMJ]T(DαRISAODβNI)+K[DαMI]T(DαRISAODβNJ)J[DαβRI](DαMISAODβNJ)+K[DβRI]T(DαMISAODβNJ)+K[DβRJ]T(DαMISAODβNI)J[DαβRJ](DαMISAODβNI)(DαβRIJ[DαβRJ])((DαMISAODβNJ)SAO)+(DαRIK[DαRJ]T)((DαMISAODβNJ)SAO)+(DβRIK[DβRJ]T)((DαMISAODβNJ)SAO)+J[DαβRI](DαRISAODβNJSAODαMJ)K[DαRJ]T(DαRISAODβNJSAODαMI)K[DαRI]T(DαRISAODβNJSAODαMJ)+J[DαβRJ](DαRISAODβNJSAODαMI)+J[DαβRI](DβRJSAODαMISAODβNJ)K[DβRJ]T(DβRJSAODαMISAODβNI)K[DβRI]T(DβRJSAODαMISAODβNJ)+J[DαβRJ](DβRJSAODαMISAODβNI)(DβRJSAODαMJ)K[DαRISAODβNI]T(DβRJSAODαMI)K[DαRISAODβNJ]T.
(55)

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

XJ[Y]=YJ[X]        XK[Y]=YK[X].
(56)

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

Both the MO and AO formulas for the newly developed Eexch(10) correction were implemented using PSI416 and the PSI4NUMPY framework.15 The agreement between the high-spin SF-SAPT results and the conventional ROHF-based SAPT in PSI46 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 set22–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 O2⋯O2 calculations were performed in the ANO-VTZ basis set25,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/JKFIT28 sets were used for Li and Mn atoms as well as for O2O2. All other atoms used the aug-cc-pVXZ/JKFIT29 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

ΔEexch(10)=Eexch(10)(S=SA+SB)Eexch(10)(S=|SASB|).
(57)

For the Mn⋯Mn complex, difficulties converging the ROHF iterations in PSI4 led us to use MOLPRO30 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(S2) results for this system will be presented below.

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 Δ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 S2 results. While it may come as a surprise, the S2 approximation in the commonly used form is not exact even in this case. During the derivation of the Eexch(10) correction in the S2 approximation,31 the following intermediate formula is reached:

VP+Eexch(10)+Eexch(10)P=VP,
(58)

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)P term on the lhs of Eq. (58) is normally neglected because it is at least of the order S4 [S2 in Eexch(10) and S2 in P]. 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 S2 and nonapproximated results for Li⋯H originate solely from the removal of the Eexch(10)P term.

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 S2 approximation. The deviation of the S2 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 S2 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 S2 results for the triplet also deviate slowly from the complete exchange results. At the chemical minimum of the triplet state (3.5 bohrs), the S2 approximation diverges considerably from the complete results, recovering 92% of the exchange energy for the triplet, 88% for the quintet, and 58% of ΔEexch(10).

For the Li dimer, we also investigated the significance of the Eexch(10)P term, which can easily be incorporated into the S2 approximation by solving Eq. (58) for Eexch(10). Even though the S2 approximation does not become exact for this system when the Eexch(10)P 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)P 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 S2 data. We suspect that the good performance of the Eexch(10)P 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) implementation4,6]. This provides an opportunity to gauge the relative accuracy of the 1-flip and S2 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 S2 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 S2 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 S2 values from the 1-flip ones for both the highest and lowest spin states of this system.

The last system that we consider in this section is the first bimolecular complex, the O2 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 Δ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 S2 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 S2 approximation, shown in Fig. 4, stems from a cancellation of beyond-S2 effects between the two spin states of the complex.

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.16Figure 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 derivatives33 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 S2 formulation of first-order SF-SAPT, together with the (much more demanding) supermolecular complete active space self consistent field theory (CASSCF).7Figure 6 presents these results along with the new nonapproximated Δ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 S2 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 S2 and 1-flip ΔEexch(10) values agree well even at shorter distances. This observation would seem to support the idea that the removal of the S2 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 S2 results for both Δ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.

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 (S2) 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-S2) 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 S2 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 S2 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 S2 variant in other cases. Another observation made possible by the new development is that the S2 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 S2 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 N4 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 S2 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 S2 approximation. This is due to a more even effect of the S2 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.

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

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.

1.
B.
Jeziorski
,
R.
Moszyński
, and
K.
Szalewicz
,
Chem. Rev.
94
,
1887
(
1994
).
2.
M.
Bartolomei
,
E.
Carmona-Novillo
,
M. I.
Hernández
,
J.
Campos-Martínez
, and
R.
Hernández-Lamoneda
,
J. Chem. Phys.
128
,
214304
(
2008
).
3.
M.
Bartolomei
,
E.
Carmona-Novillo
,
M. I.
Hernández
,
J.
Campos-Martínez
, and
R.
Hernández-Lamoneda
,
J. Chem. Phys.
133
,
124311
(
2010
).
4.
P. S.
Żuchowski
,
R.
Podeszwa
,
R.
Moszyński
,
B.
Jeziorski
, and
K.
Szalewicz
,
J. Chem. Phys.
129
,
084101
(
2008
).
5.
M.
Hapka
,
P. S.
Żuchowski
,
M. M.
Szczȩśniak
, and
G.
Chałasiński
,
J. Chem. Phys.
137
,
164104
(
2012
).
6.
J. F.
Gonthier
and
C. D.
Sherrill
,
J. Chem. Phys.
145
,
134106
(
2016
).
7.
K.
Patkowski
,
P. S.
Żuchowski
, and
D. G. A.
Smith
,
J. Chem. Phys.
148
,
164110
(
2018
).
8.
A. I.
Krylov
,
Chem. Phys. Lett.
338
,
375
(
2001
).
10.
P. M.
Zimmerman
,
F.
Bell
,
M.
Goldey
,
A. T.
Bell
, and
M.
Head-Gordon
,
J. Chem. Phys.
137
,
164110
(
2012
).
11.
N. J.
Mayhall
and
M.
Head-Gordon
,
J. Phys. Chem. Lett.
6
,
1982
(
2015
).
12.
B.
Jeziorski
,
M.
Bulski
, and
L.
Piela
,
Int. J. Quantum Chem.
10
,
281
(
1976
).
13.
R.
Schäffer
and
G.
Jansen
,
Theor. Chem. Acc.
131
,
1235
(
2012
).
14.
R.
Schäffer
and
G.
Jansen
,
Mol. Phys.
111
,
2570
(
2013
).
15.
D. G. A.
Smith
,
L. A.
Burns
,
D. A.
Sirianni
,
D. R.
Nascimento
,
A.
Kumar
,
A. M.
James
,
J. B.
Schriber
,
T.
Zhang
,
B.
Zhang
,
A. S.
Abbott
 et al,
J. Chem. Theory Comput.
14
,
3504
(
2018
).
16.
R. M.
Parrish
,
L. A.
Burns
,
D. G. A.
Smith
,
A. C.
Simmonett
,
A. E.
DePrince
 III
,
E. G.
Hohenstein
,
U.
Bozkaya
,
A. Y.
Sokolov
,
R.
Di Remigio
,
R. M.
Richard
 et al,
J. Chem. Theory Comput.
13
,
3185
(
2017
).
17.
M.
Kertesz
,
Chem. Eur. J.
25
,
400
(
2019
).
18.
B.
Jeziorski
,
G.
Chałasiński
, and
K.
Szalewicz
,
Int. J. Quantum Chem.
14
,
271
(
1978
).
19.
R.
Vein
and
P.
Dale
,
Determinants and Their Applications in Mathematical Physics
(
Springer
,
New York
,
1999
).
20.
K. U.
Lao
and
J. M.
Herbert
,
J. Chem. Theory Comput.
14
,
2955
(
2018
).
21.
R. M.
Parrish
,
K. C.
Thompson
, and
T. J.
Martínez
,
J. Chem. Theory Comput.
14
,
1737
(
2018
).
22.
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
90
,
1007
(
1989
).
23.
R. A.
Kendall
,
T. H.
Dunning
, Jr.
, and
R. J.
Harrison
,
J. Chem. Phys.
96
,
6796
(
1992
).
24.
N. B.
Balabanov
and
K. A.
Peterson
,
J. Chem. Phys.
123
,
064107
(
2005
).
25.
P.-O.
Widmark
,
P.-Å.
Malmqvist
, and
B. O.
Roos
,
Theor. Chem. Acc.
77
,
291
(
1990
).
26.
M.
Bartolomei
,
M. I.
Hernández
,
J.
Campos-Martínez
,
E.
Carmona-Novillo
, and
R.
Hernández-Lamoneda
,
Phys. Chem. Chem. Phys.
10
,
5374
(
2008
).
27.
K.
Patkowski
,
T.
Korona
, and
B.
Jeziorski
,
J. Chem. Phys.
115
,
1137
(
2001
).
28.
F.
Weigend
,
J. Comput. Chem.
29
,
167
(
2007
).
29.
F.
Weigend
,
A.
Köhn
, and
C.
Hättig
,
J. Chem. Phys.
116
,
3175
(
2002
).
30.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
, and
M.
Schütz
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
242
(
2012
).
31.
R.
Moszyński
,
B.
Jeziorski
,
S.
Rybak
,
K.
Szalewicz
, and
H. L.
Williams
,
J. Chem. Phys.
100
,
5080
(
1994
).
32.
Z.
Cui
,
H.
Lischka
,
H. Z.
Beneberu
, and
M.
Kertesz
,
J. Am. Chem. Soc.
136
,
5539
(
2014
).
33.
Z.
Mou
,
T.
Kubo
, and
M.
Kertesz
,
Chem. Eur. J.
21
,
18230
(
2015
).
34.
Z.
Mou
and
M.
Kertesz
,
Angew. Chem., Int. Ed.
56
,
10188
(
2017
).
35.
P. G.
Szalay
and
R. J.
Bartlett
,
Chem. Phys. Lett.
214
,
481
(
1993
).
36.
Y.
Zhao
and
D. G.
Truhlar
,
J. Phys. Chem. A
110
,
5121
(
2006
).
37.
Z.
Mou
,
Y.-H.
Tian
, and
M.
Kertesz
,
Phys. Chem. Chem. Phys.
19
,
24761
(
2017
).

Supplementary Material