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

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

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

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

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

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

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

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

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

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

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

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

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

in powers of the intermolecular interaction operator V=iAjBṽ(ij) and separating the polarization and exchange effects in each order. In Eq. (1), A is the (NA + NB)-electron antisymmetrizer and ΦAB is the polarization (RS) expansion of the wavefunction for the complex.24 We will assume that ΨA and ΨB are ground-state high-spin open-shell determinants. The occupied orbitals in ΨA are χi1, χi2, …, χikA (inactive, doubly occupied) and χa1, χa2, …, χa2SA (active, occupied by α electrons only). The occupied orbitals in ΨB are χj1, χj2, …, χjkB (inactive, doubly occupied) and χb1, χb2, …, χb2SB (active, occupied by β electrons only).

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

Eint=PSMSΨAΨB|VA|ΦABPSMSΨAΨB|A|ΦAB.
(2)

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

We now introduce the single exchange approximation (often called the S2 approximation as it neglects terms higher than quadratic in the intermolecular overlap integrals),

A(1+P)AAAB,
(3)

where the single-exchange operator P is given by

P=iAjBPij
(4)

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

V+Eexch(10)PSMS+PPSMS=VPSMS+VPPSMS.
(5)

Note that we have separated the electrostatic contribution Eelst(10)=V from the exchange one, and we introduced the zero into the E(10) superscript to remind the reader that no intramolecular correlation effects are included.

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

PSMSΨAΨB=c0ΨAΨB+c1ΨAΨB+c2ΨAΨB+,
(6)

where, for example, ΨA is a normalized wavefunction (linear combination of determinants) corresponding to the spin quantum numbers (SA, SA − 1). Up to a constant, this function can be obtained from ΨA by the action of the spin-lowering operator Ŝ.

Equation (6) implies that

PSMS=c0,
(7)
VPSMS=c0V.
(8)

Therefore, Eq. (5) becomes

VPPSMS+c0Eexch(10)+Eexch(10)PPSMS=VPPSMS.
(9)

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

Eexch(10)=c01VPPSMSVPPSMS.
(10)

Let us now go back to Eq. (6) and examine the spin-flipped monomer wavefunctions such as ΨA. Up to normalization, this function is equal to ŜΨA, where Ŝ=k=1NAŜ(k) applies the conventional spin-lowering operator Ŝ(k) = Ŝx(k) − iŜy(k) to each of the NA spins. It can be shown that when ΨA is a restricted high-spin determinant like in this work, the lowering operator acts on it as follows:

ŜΨA=n=12SAχi1αχi1βχikAαχikAβχa1α×χan1αχanβχan+1αχa2SAα,
(11)

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

Inserting Eq. (6) into Eq. (10), we obtain

Eexch(10)=VP+c1/c0ΨAΨB|VP|ΨAΨB+c2/c0ΨAΨB|VP|ΨAΨB++VP+c1/c0ΨAΨB|P|ΨAΨB+c2/c0ΨAΨB|P|ΨAΨB+.
(12)

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

VPVP=ṽαββα+Sαβ(ṽαβααṽαβαα)+Sβα(ṽαβββṽαβββ)SαβSβαṽαβαβSαβSβαṽαβαβ+SαβSβαṽαβαβ.
(13)

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

ṽμνλκ=(λμ|κν)+(vA)νκSμλ/NA+(vB)μλSνκ/NB+V0SμλSνκ/NANB,
(14)

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

Eexch,diag(10)=VPVP=ṽajjaṽibbi2ṽijji+ṽajaaSajṽajaaSaj+ṽajjjSja2ṽajjjSja+ṽijaiSaj2ṽijiaSajṽabjbSja+ṽibbbSbiṽibbbSbi+ṽibiiSib2ṽibiiSib+ṽijjbSbi2ṽijbjSbiṽabaiSib+ṽajiaSij+ṽibbjSji2ṽajaiSij2ṽibjbSji+2ṽijiiSij+2ṽijjjSji4ṽijiiSij4ṽijjjSji+ṽababSajSja+ṽababSibSbi+2ṽibabSajSji+2ṽajabSijSbi+2ṽajajSajSja+2ṽajajSijSji+ṽajajSajSjaṽajajSajSja+2ṽibibSibSbi+2ṽibibSijSji+ṽibibSibSbiṽibibSibSbi+4ṽijajSajSji2ṽijajSajSji+4ṽijibSijSbi2ṽijibSijSbi+2ṽijijSibSbi+4ṽijijSijSji+2ṽijijSajSja+4ṽijijSijSji2ṽijijSijSji,
(15)

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

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

ΨAΨB|VP|ΨAΨB=ṽ(12)ρint(12)dτ12,
(16)
ΨAΨB|P|ΨAΨB=1NANBρint(12)dτ12,
(17)

where

ρint(12)=NANBΨAΨBPΨAΨBdτ12
(18)

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

ρint(12)=ρA(1|2)ρB(2|1)ρA(1|4)ΓB(24|21)dτ4ΓA(13|12)ρB(2|3)dτ3ΓA(13|14)ΓB(24|23)dτ3dτ4.
(19)

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

ρA(1|1)=NAΨA*(1,2,,NA)ΨA(1,2,,NA)dτ1,
(20)
ΓA(12|12)=NA(NA1)ΨA*(1,2,3,,NA)×ΨA(1,2,3,,NA)dτ12.
(21)

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

ρA(1|1)=12SAn=12SAχanα(1)χanβ(1).
(22)

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

ΓA(12|12)=12SAn=12SAm=1kAχanα(1)χimα(2)χanβ(1)χimα(2)χimα(1)χanα(2)χanβ(1)χimα(2)χanα(1)χimα(2)χimα(1)χanβ(2)+χimα(1)χanα(2)χimα(1)χanβ(2)+m=1kAχanα(1)χimβ(2)χanβ(1)χimβ(2)χimβ(1)χanα(2)χanβ(1)χimβ(2)χanα(1)χimβ(2)χimβ(1)χanβ(2)+χimβ(1)χanα(2)χimβ(1)χanβ(2)+m=12SAχanα(1)χamα(2)χanβ(1)χamα(2)χamα(1)χanα(2)χanβ(1)χamα(2)χanα(1)χamα(2)χamα(1)χanβ(2)+χamα(1)χanα(2)χamα(1)χanβ(2).
(23)

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

ΨAΨB|VP|ΨAΨB=12SASBṽabba+ṽabbjSja+ṽajjbSba2ṽajbjSba+ṽabbbSbaṽabbbSba+ṽabiaSib+ṽibaiSab2ṽiaibSba+ṽabaaSabṽaaabSba2ṽabijSibSja2ṽajibSijSba+4ṽajijSibSba+4ṽibijSabSja4ṽijijSabSba+2ṽibibSabSba2ṽabibSibSba2ṽibibSabSba+2ṽibabSabSbi+2ṽajajSabSba2ṽabajSabSja2ṽajajSabSba+2ṽabajSabSja+ṽababSabSbaṽababSabSbaṽababSabSba+ṽababSabSba.
(24)

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

ΨAΨB|P|ΨAΨB=12SASBSabSba.
(25)

The ROHF electrostatic energy is equal to

Eelst(10)=V=4ṽijij+2ṽajaj+2ṽibib+ṽabab.
(26)

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

12SASBEexch,1flip(10)=ΨAΨB|VP|ΨAΨBVΨAΨB|P|ΨAΨB=12SASBṽabba+ṽabbjSja+ṽajjbSba2ṽajbjSba+ṽabbbSbaṽabbbSba+ṽabiaSib+ṽibaiSab2ṽiaibSba+ṽabaaSabṽaaabSba2ṽabijSibSja2ṽajibSijSba+4ṽajijSibSba+4ṽibijSabSja+2ṽibibSabSba2ṽabibSibSba+2ṽibabSabSbi+2ṽajajSabSba2ṽabajSabSja+2ṽabajSabSja+ṽababSabSbaṽababSabSba+ṽababSabSba.
(27)

Finally, we will show that the double spin-flip term in Eq. (12), and all subsequent terms, vanishes identically. Specifically, we will prove that ΨAΨB|VP|ΨAΨB=0, where ΨA and ΨB are proportional to ŜŜΨA and Ŝ+Ŝ+ΨB, respectively. Similar to the previous term,

ΨAΨB|VP|ΨAΨB=ṽ(12)ρint(12)dτ12,
(28)

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

ρA(1|1)=0,
(29)

and the two-electron density matrix is proportional to

ΓA(12|12)m,n=12SAχanα(1)χamα(2)χanβ(1)χamβ(2)χamα(1)χanα(2)χanβ(1)χamβ(2)χanα(1)χamα(2)χamβ(1)χanβ(2)+χamα(1)χanα(2)χamβ(1)χanβ(2).
(30)

According to Eq. (29), only the last term in the expression for ρint survives,

ρint(12)=ΓA(13|14)ΓB(24|23)dτ3dτ4.
(31)

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

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

Eexch(10)=VPVP+c1/c0ΨAΨB|VP|ΨAΨBVΨAΨB|P|ΨAΨB=Eexch,diag(10)+Z(SA,SB,S)Eexch,1flip(10),
(32)

where

Z(SA,SB,S)=S(S+1)+2SASBSA(SA+1)SB(SB+1)4SASB.
(33)

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

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

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

J[X]KL=(KL|MN)XMN,
(34)
K[X]KL=(KM|NL)XMN
(35)

and reduce to standard Coulomb and exchange matrices, or their inactive/active contributions, in special cases, for example, JiA = J[PiA], KaA = K[PaA], …. The electrostatic potential matrix for monomer A is

ωA=vA+2JiA+JaA,
(36)

where vA is the matrix of the nuclear attraction operator. By adding the exchange matrices, we obtain the following spin-dependent Fock matrices h:

hAα=vA+2JiA+JaAKiAKaA,
(37)
hAβ=vA+2JiA+JaAKiA,
(38)
hBα=vB+2JiB+JaBKiB,
(39)
hBβ=vB+2JiB+JaBKiBKaB.
(40)

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

Eexch,diag(10)=PiB(2KiA+KaA)PaBKiA(PiASAOPiB)(hAα+hAβ+hBα+hBβ)(PaASAOPiB)(hAα+hBα)(PiASAOPaB)(hAβ+hBβ)+2(PiBSAOPiASAOPaB)ωA+2(PiBSAOPiASAOPiB)ωA+(PaBSAOPiASAOPaB)ωA+(PiBSAOPaASAOPiB)ωA+2(PiASAOPiBSAOPiA)ωB+2(PiASAOPiBSAOPaA)ωB+(PiASAOPaBSAOPiA)ωB+(PaASAOPiBSAOPaA)ωB2(PiASAOPiB)K[PiASAOPiB]2(PaASAOPiB)K[PiASAOPiB]2(PiASAOPaB)K[PiASAOPiB](PaASAOPiB)K[PaASAOPiB](PiASAOPaB)K[PiASAOPaB],
(41)
Eexch,1flip(10)=PaBKaA(PaASAOPaB)(hAα+hBβ)+(PaASAOPiB)KaB+(PiASAOPaB)KaA+2(PiBSAOPaASAOPaB)ωA+(PaBSAOPaASAOPaB)ωA+2(PiASAOPaBSAOPaA)ωB+(PaASAOPaBSAOPaA)ωB2(PaASAOPaB)K[PiASAOPiB]2(PaASAOPaB)K[PiASAOPaB]2(PaASAOPiB)K[PaASAOPaB]2(PaASAOPiB)K[PiASAOPaB](PaASAOPaB)K[PaASAOPaB],
(42)

where SAO is the overlap matrix in the AO basis.

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

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

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

ΔE(10)=(1+12max(SA,SB))Eexch,1flip(10),
(43)

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

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

FIG. 1.

First-order exchange energy for the singlet and triplet states of the Li⋯H complex, computed within the ROHF-based approach of this work (points) and the FCI-based SAPT of Ref. 38 (lines). The singlet-triplet splitting from the FCI calculations is shown for comparison.

FIG. 1.

First-order exchange energy for the singlet and triplet states of the Li⋯H complex, computed within the ROHF-based approach of this work (points) and the FCI-based SAPT of Ref. 38 (lines). The singlet-triplet splitting from the FCI calculations is shown for comparison.

Close modal

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

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

FIG. 2.

Comparison of the ΔE(10) energy splittings for the Li⋯Li, Li⋯N, and N⋯N systems with the supermolecular CASSCF and (Davidson corrected) MRCI results. For the Li dimer, experimentally derived values for the splitting given by Côté et al.64 were also plotted.

FIG. 2.

Comparison of the ΔE(10) energy splittings for the Li⋯Li, Li⋯N, and N⋯N systems with the supermolecular CASSCF and (Davidson corrected) MRCI results. For the Li dimer, experimentally derived values for the splitting given by Côté et al.64 were also plotted.

Close modal

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

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

FIG. 3.

Quality of the single exchange approximation for the high-spin first-order exchange energy in the Li⋯Li, Li⋯N, and N⋯N systems measured as the ratio Eexch(10)(S2)/Eexch(10). Color-coded vertical dashed lines correspond to the positions of minima for the high-spin complexes.

FIG. 3.

Quality of the single exchange approximation for the high-spin first-order exchange energy in the Li⋯Li, Li⋯N, and N⋯N systems measured as the ratio Eexch(10)(S2)/Eexch(10). Color-coded vertical dashed lines correspond to the positions of minima for the high-spin complexes.

Close modal

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

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

FIG. 4.

Comparison of the first-order exchange splittings (defined here as the quintet energy minus the singlet energy) with literature data for two interacting ground-state oxygen molecules (3Σg) for four basic geometric configurations. Since the exchange splittings for the X-shape configurations are much smaller compared to other orientations and they sometimes change sign, they are plotted on a non-logarithmic scale. The data marked “Wormer et al.” are from Ref. 36 and the CASSCF/MRCI/CASPT2 results are from Ref. 29.

FIG. 4.

Comparison of the first-order exchange splittings (defined here as the quintet energy minus the singlet energy) with literature data for two interacting ground-state oxygen molecules (3Σg) for four basic geometric configurations. Since the exchange splittings for the X-shape configurations are much smaller compared to other orientations and they sometimes change sign, they are plotted on a non-logarithmic scale. The data marked “Wormer et al.” are from Ref. 36 and the CASSCF/MRCI/CASPT2 results are from Ref. 29.

Close modal
FIG. 5.

Diagonal and spin-flip parts of the first-order exchange energy for the dimer of ground-state oxygen molecules (3Σg) for four basic geometric configurations. For the X-shape configurations, the absolute value of the spin-flip term is shown, hence the dip around 7a0 corresponding to a change of sign from positive (for small R) to negative.

FIG. 5.

Diagonal and spin-flip parts of the first-order exchange energy for the dimer of ground-state oxygen molecules (3Σg) for four basic geometric configurations. For the X-shape configurations, the absolute value of the spin-flip term is shown, hence the dip around 7a0 corresponding to a change of sign from positive (for small R) to negative.

Close modal
FIG. 6.

The quality of the S2 approximation for the high-spin quintet state of the molecular oxygen dimer. The oxygen molecules are in their ground 3Σg state and their mutual orientation corresponds to four basic geometric configurations.

FIG. 6.

The quality of the S2 approximation for the high-spin quintet state of the molecular oxygen dimer. The oxygen molecules are in their ground 3Σg state and their mutual orientation corresponds to four basic geometric configurations.

Close modal

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

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

FIG. 7.

First-order exchange energy for the Mn⋯Mn complex: (a) diagonal and spin-flip contributions to the exchange energies compared with the high-spin exchange energy; the extremely small contribution of spin-flip exchange is evident; (b) the quality of the S2 approximation for the high-spin (undecaplet) Mn⋯Mn state; (c) comparison of the Heisenberg JAB parameters derived from SF-SAPT with existing theoretical and experimental data; note that only in the study of Buchachenko et al.72 the JAB parameter is available as a function of R, whereas other data are provided for the respective equilibrium distances. For the experimental value, the distance at which JAB is marked corresponds to the CCSD(T) 1Σg+ minimum obtained in Ref. 72.

FIG. 7.

First-order exchange energy for the Mn⋯Mn complex: (a) diagonal and spin-flip contributions to the exchange energies compared with the high-spin exchange energy; the extremely small contribution of spin-flip exchange is evident; (b) the quality of the S2 approximation for the high-spin (undecaplet) Mn⋯Mn state; (c) comparison of the Heisenberg JAB parameters derived from SF-SAPT with existing theoretical and experimental data; note that only in the study of Buchachenko et al.72 the JAB parameter is available as a function of R, whereas other data are provided for the respective equilibrium distances. For the experimental value, the distance at which JAB is marked corresponds to the CCSD(T) 1Σg+ minimum obtained in Ref. 72.

Close modal

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

FIG. 8.

The singlet-triplet energy splitting for the staggered conformation of the phenalenyl dimer. The first-order SAPT calculations and the CASSCF(2,2) calculations use the aug-cc-pVDZ basis set. The MR-AQCC(2,2)/6-31G(d) splitting, computed in Ref. 31 for the minimum separation of the singlet complex, is included for comparison.

FIG. 8.

The singlet-triplet energy splitting for the staggered conformation of the phenalenyl dimer. The first-order SAPT calculations and the CASSCF(2,2) calculations use the aug-cc-pVDZ basis set. The MR-AQCC(2,2)/6-31G(d) splitting, computed in Ref. 31 for the minimum separation of the singlet complex, is included for comparison.

Close modal

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

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

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

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

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

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

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

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

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

c0=S(SASB)|SASASBSB,
(A1)
c1=S(SASB)|SA(SA1)SB(SB+1).
(A2)

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

We start by noting that the action of the spin-flip operator ŜAŜB+ on the initial product state |SASASBSB⟩ produces

ŜAŜB+|SASASBSB=SA(SA+1)SA(SA1)SB(SB+1)+SB(1SB)|SA(SA1)SB(SB+1)=2SASB|SA(SA1)SB(SB+1).
(A3)

Since Ŝ2=ŜA2+ŜB2+2(ŜA)z(ŜB)z+ŜAŜB++ŜA+ŜB, one might rewrite the action of the operator ŜAŜB+ as

ŜAŜB+|SASASBSB=Ŝ2(ŜA2+ŜB2+2(ŜA)z(ŜB)z+ŜA+ŜB)|SASASBSB=Ŝ2(SA(SA+1)+SB(SB+1)2SASB+0)|SASASBSB.
(A4)

The last term in the above equation yields zero since the reverse spin-flip operator ŜA+ŜB cannot raise any more spins on A or lower any more spins on B. Now, by projecting Eqs. (A3) and (A4) onto ⟨S (SASB)| and using the fact that ⟨S (SASB)|Ŝ2 = S(S + 1)⟨S (SASB)|, one finds that

2c1SASB=S(S+1)SA(SA+1)SB(SB+1)+2SASBc0,
(A5)

which, after a simple rearrangement, gives the formula (33) for Z(SA, SB, S).

1.
S.
Aloisio
and
J. S.
Francisco
,
Acc. Chem. Res.
33
,
825
(
2000
).
2.
A.
Galano
,
M.
Narciso-Lopez
, and
M.
Francisco-Marquez
,
J. Phys. Chem. A
114
,
5796
(
2010
).
3.
E.
Roueff
and
F.
Lique
,
Chem. Rev.
113
,
8906
(
2013
).
4.
I.
Ratera
and
J.
Veciana
,
Chem. Soc. Rev.
41
,
303
(
2012
).
5.
M.
Kirste
,
X.
Wang
,
H. C.
Schewe
,
G.
Meijer
,
K.
Liu
,
A.
van der Avoird
,
L. M.
Janssen
,
K. B.
Gubbels
,
G. C.
Groenenboom
, and
S. Y.
van de Meerakker
,
Science
338
,
1060
(
2012
).
6.
E.
Lavert-Ofir
,
Y.
Shagam
,
A. B.
Henson
,
S.
Gersten
,
J.
Kłos
,
P. S.
Żuchowski
,
J.
Narevicius
, and
E.
Narevicius
,
Nat. Chem.
6
,
332
(
2014
).
7.
S.
Chefdeville
,
Y.
Kalugina
,
S. Y. T.
van de Meerakker
,
C.
Naulin
,
F.
Lique
, and
M.
Costes
,
Science
341
,
1094
(
2013
).
8.
P. J.
Knowles
,
C.
Hampel
, and
H.-J.
Werner
,
J. Chem. Phys.
99
,
5219
(
1993
).
9.
J. D.
Watts
,
J.
Gauss
, and
R. J.
Bartlett
,
J. Chem. Phys.
98
,
8718
(
1993
).
10.
B.
Jeziorski
and
H. J.
Monkhorst
,
Phys. Rev. A
24
,
1668
(
1981
).
11.
R. J.
Bartlett
and
M.
Musiał
,
Rev. Mod. Phys.
79
,
291
(
2007
).
12.
13.
H.-J.
Werner
and
P. J.
Knowles
,
J. Chem. Phys.
89
,
5803
(
1988
).
14.
K.
Andersson
,
P.-A.
Malmqvist
,
B.
Roos
,
A.
Sadlej
, and
K.
Wolinski
,
J. Phys. Chem.
94
,
5483
(
1990
).
15.
J.
Finley
,
P.
Malmqvist
,
B. O.
Roos
, and
L.
Serrano-Andrés
,
Chem. Phys. Lett.
288
,
299
(
1998
).
16.
C.
Angeli
,
R.
Cimiraglia
, and
J.-P.
Malrieu
,
J. Chem. Phys.
117
,
9138
(
2002
).
17.
C.
Camacho
,
R.
Cimiraglia
, and
H. A.
Witek
,
Phys. Chem. Chem. Phys.
12
,
5058
(
2010
).
18.
A. I.
Krylov
,
Chem. Phys. Lett.
338
,
375
(
2001
).
19.
Y.
Shao
,
M.
Head-Gordon
, and
A. I.
Krylov
,
J. Chem. Phys.
118
,
4807
(
2003
).
20.
S. V.
Levchenko
and
A. I.
Krylov
,
J. Chem. Phys.
120
,
175
(
2004
).
21.
A. I.
Krylov
,
Acc. Chem. Res.
39
,
83
(
2006
).
22.
P. M.
Zimmerman
,
F.
Bell
,
M.
Goldey
,
A. T.
Bell
, and
M.
Head-Gordon
,
J. Chem. Phys.
137
,
164110
(
2012
).
23.
N. J.
Mayhall
and
M.
Head-Gordon
,
J. Phys. Chem. Lett.
6
,
1982
(
2015
).
24.
B.
Jeziorski
,
R.
Moszyński
, and
K.
Szalewicz
,
Chem. Rev.
94
,
1887
(
1994
).
25.
E. G.
Hohenstein
and
C. D.
Sherrill
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
304
(
2012
).
26.
P. S.
Żuchowski
,
R.
Podeszwa
,
R.
Moszyński
,
B.
Jeziorski
, and
K.
Szalewicz
,
J. Chem. Phys.
129
,
084101
(
2008
).
27.
M.
Hapka
,
P. S.
Żuchowski
,
M. M.
Szczȩśniak
, and
G.
Chałasiński
,
J. Chem. Phys.
137
,
164104
(
2012
).
28.
V.
Aquilanti
,
D.
Ascenzi
,
M.
Bartolomei
,
D.
Cappeletti
,
S.
Cavalli
,
M.
de Castro Vitores
, and
F.
Pirani
,
J. Am. Chem. Soc.
121
,
10794
(
1999
).
29.
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
).
30.
M.
Bartolomei
,
E.
Carmona-Novillo
,
M. I.
Hernández
,
J.
Campos-Martínez
, and
R.
Hernández-Lamoneda
,
J. Chem. Phys.
133
,
124311
(
2010
).
31.
Z.
Cui
,
H.
Lischka
,
H. Z.
Beneberu
, and
M.
Kertesz
,
J. Am. Chem. Soc.
136
,
5539
(
2014
).
32.
M. V.
Parkes
,
D. F.
Sava Gallis
,
J. A.
Greathouse
, and
T. M.
Nenoff
,
J. Phys. Chem. C
119
,
6556
(
2015
).
33.
M.
Bartenstein
,
A.
Altmeyer
,
S.
Riedl
,
R.
Geursen
,
S.
Jochim
,
C.
Chin
,
J.
Hecker Denschlag
,
R.
Grimm
,
A.
Simoni
,
E.
Tiesinga
 et al,
Phys. Rev. Lett.
94
,
103201
(
2005
).
34.
C.
Chin
,
V.
Vuletić
,
A. J.
Kerman
,
S.
Chu
,
E.
Tiesinga
,
P. J.
Leo
, and
C. J.
Williams
,
Phys. Rev. A
70
,
032701
(
2004
).
35.
H. L.
Williams
and
C. F.
Chabalowski
,
J. Phys. Chem. A
105
,
646
(
2001
).
36.
P. E. S.
Wormer
and
A.
van der Avoird
,
J. Chem. Phys.
81
,
1929
(
1984
).
37.
T.
Ćwiok
,
B.
Jeziorski
,
W.
Kołos
,
R.
Moszyński
, and
K.
Szalewicz
,
J. Chem. Phys.
97
,
7555
(
1992
).
38.
K.
Patkowski
,
T.
Korona
, and
B.
Jeziorski
,
J. Chem. Phys.
115
,
1137
(
2001
).
39.
K.
Patkowski
,
B.
Jeziorski
,
T.
Korona
, and
K.
Szalewicz
,
J. Chem. Phys.
117
,
5124
(
2002
).
40.
P.
Gniewek
and
B.
Jeziorski
,
Phys. Rev. A
90
,
022506
(
2014
).
41.
P.
Gniewek
and
B.
Jeziorski
,
J. Chem. Phys.
143
,
154106
(
2015
).
42.
P.
Gniewek
and
B.
Jeziorski
,
Phys. Rev. A
94
,
042708
(
2016
).
43.
B.
Jeziorski
,
G.
Chałasiński
, and
K.
Szalewicz
,
Int. J. Quantum Chem.
14
,
271
(
1978
).
44.
B.
Jeziorski
,
M.
Bulski
, and
L.
Piela
,
Int. J. Quantum Chem.
10
,
281
(
1976
).
45.
R.
Schäffer
and
G.
Jansen
,
Theor. Chem. Acc.
131
,
1235
(
2012
).
46.
R.
Schäffer
and
G.
Jansen
,
Mol. Phys.
111
,
2570
(
2013
).
47.
R.
Moszyński
,
B.
Jeziorski
,
S.
Rybak
,
K.
Szalewicz
, and
H. L.
Williams
,
J. Chem. Phys.
100
,
5080
(
1994
).
48.
K.
Patkowski
,
K.
Szalewicz
, and
B.
Jeziorski
,
J. Chem. Phys.
125
,
154107
(
2006
).
49.
F. A.
Matsen
,
D. J.
Klein
, and
D. C.
Foyt
,
J. Phys. Chem.
75
,
1866
(
1971
).
50.
A.
Hesselmann
,
G.
Jansen
, and
M.
Schütz
,
J. Chem. Phys.
122
,
014103
(
2005
).
51.
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
).
52.
H. L.
Williams
,
E. M.
Mas
,
K.
Szalewicz
, and
B.
Jeziorski
,
J. Chem. Phys.
103
,
7374
(
1995
).
53.
SAPT2012: An Ab Initio Program for Many-Body Symmetry-Adapted Perturbation Theory Calculations of Intermolecular Interaction Energies
, by
R.
Bukowski
,
W.
Cencek
,
P.
Jankowski
,
M.
Jeziorska
,
B.
Jeziorski
,
S. A.
Kucharski
,
V. F.
Lotrich
,
A. J.
Misquitta
,
R.
Moszyński
,
K.
Patkowski
,
R.
Podeszwa
,
F.
Rob
,
S.
Rybak
,
K.
Szalewicz
,
H. L.
Williams
,
R. J.
Wheatley
,
P. E. S.
Wormer
, and
P. S.
Żuchowski
(
University of Delaware and University of Warsaw
,
2012
), http://www.physics.udel.edu/∼szalewic/SAPT/SAPT.html.
54.
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, https://chemrxiv.org/articles/Psi4NumPy_An_Interactive_Quantum_Chemistry_Programming_Environment_for_Reference_Implementations_and_Rapid_Development/5746059,
2018
.
55.
F.
Weigend
,
A.
Köhn
, and
C.
Hättig
,
J. Chem. Phys.
116
,
3175
(
2002
).
56.
R. A.
Kendall
,
T. H.
Dunning
, Jr.
, and
R. J.
Harrison
,
J. Chem. Phys.
96
,
6796
(
1992
).
57.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
, and
M.
Schütz
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
242
(
2012
).
58.
W.
Kutzelnigg
,
J. Chem. Phys.
73
,
343
(
1980
).
59.
B.
Jeziorski
, private communication (
2018
).
60.
M.
Przybytek
,
K.
Patkowski
, and
B.
Jeziorski
,
Collect. Czech. Chem. Commun.
69
,
141
(
2004
).
61.
M.
Przybytek
and
B.
Jeziorski
,
J. Chem. Phys.
123
,
134315
(
2005
).
62.
M.
Semczuk
,
X.
Li
,
W.
Gunton
,
M.
Haw
,
N. S.
Dattani
,
J.
Witz
,
A. K.
Mills
,
D. J.
Jones
, and
K. W.
Madison
,
Phys. Rev. A
87
,
052505
(
2013
).
63.
T. V.
Tscherbul
,
J.
Kłos
,
A.
Dalgarno
,
B.
Zygelman
,
Z.
Pavlovic
,
M. T.
Hummon
,
H.-I.
Lu
,
E.
Tsikata
, and
J. M.
Doyle
,
Phys. Rev. A
82
,
042718
(
2010
).
64.
R.
Côté
,
A.
Dalgarno
, and
M. J.
Jamieson
,
Phys. Rev. A
50
,
399
(
1994
).
65.
R.
Ahlrichs
,
M.
Hoffmann-Ostenhof
,
T.
Hoffmann-Ostenhof
, and
J. D.
Morgan
,
Phys. Rev. A
23
,
2106
(
1981
).
66.
R. K.
Nesbet
,
Phys. Rev.
135
,
A460
(
1964
).
67.
Z.
Wang
,
R. R.
Lucchese
, and
J. W.
Bevan
,
J. Phys. Chem. A
108
,
2884
(
2004
).
68.
S.
Yamamoto
,
H.
Tatewaki
,
H.
Moriyama
, and
H.
Nakano
,
J. Chem. Phys.
124
,
124302
(
2006
).
69.
I.
Negodaev
,
C.
de Graaf
, and
R.
Caballol
,
Chem. Phys. Lett.
458
,
290
(
2008
).
70.
D.
Tzeli
,
U.
Miranda
,
I. G.
Kaplan
, and
A.
Mavridis
,
J. Chem. Phys.
129
,
154310
(
2008
).
71.
C.
Camacho
,
S.
Yamamoto
, and
H. A.
Witek
,
Phys. Chem. Chem. Phys.
10
,
5128
(
2008
).
72.
A. A.
Buchachenko
,
G.
Chałasiński
, and
M. M.
Szczȩśniak
,
J. Chem. Phys.
132
,
024312
(
2010
).
73.
M.
Cheeseman
,
R. J.
Van Zee
,
H. L.
Flanagan
, and
W.
Weltner
,
J. Chem. Phys.
92
,
1553
(
1990
).
75.
P. G.
Szalay
and
R. J.
Bartlett
,
Chem. Phys. Lett.
214
,
481
(
1993
).
76.
Z.
Mou
,
Y.-H.
Tian
, and
M.
Kertesz
,
Phys. Chem. Chem. Phys.
19
,
24761
(
2017
).
77.
A. J.
Misquitta
and
K.
Szalewicz
,
J. Chem. Phys.
122
,
214109
(
2005
).