First-order symmetry-adapted perturbation theory for multiplet splittings

Patkowski, Konrad; Żuchowski, Piotr S.; Smith, Daniel G. A.

doi:10.1063/1.5021891

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

I. INTRODUCTION

The interactions of open-shell atoms, molecules, and ions are of great significance to many areas of molecular physics and chemistry. As a matter of fact, all chemical reactions at some stage have to experience bond breaking and rearrangement and to describe such processes one has to deal with interacting fragments of reactants or products. Many electronically excited states of molecules have open electronic shells, so in studies of the interactions of such molecules with background gas or solvent, similar problems arise. Interacting radicals are of key importance in atmosphere, where they undergo chain reactions as well as form metastable states and stable complexes.^1,2 Collisions of radicals (such as OH, CH₂, CN, and many more) with hydrogen (H and H₂) are a subject of study for astrochemistry³ as they provide information about the conditions which molecules experience in the interstellar clouds. Stabilized organic radicals are relevant in organic and materials chemistry,⁴ 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.⁸ and Watts et al.,⁹ 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,¹⁰ the proposed variations of MRCC are far from being black-box methods (see Refs. 11 and 12 for a review). At present, to account for both static and dynamical correlation effects in molecular interactions for systems exhibiting a non-unique electronic configuration, multireference configuration interaction¹³ (MRCI) or multireference perturbation theories^14–16 are often applied. In the case of the MRCI method, it is possible to reproduce the short-range part of the potential well when size-consistency corrections are applied. Such a posteriori corrections are not exact and non-unique, and the remaining size-consistency error might be on the order of the interaction energy. In the case of the perturbation methods, the same problem also arises in general; moreover, some of these methods suffer from the problem of intruder states.¹⁷ 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.¹⁸ The spin-flip equation of motion approach has been implemented within a variety of different theories, in particular, Hartree-Fock,¹⁸ density functional theory,¹⁹ configuration interaction doubles,¹⁸ and finally, using the CC framework.^20,21 The theory has been generalized to multiple spin flips;²² however, the resulting calculations can be quite demanding, and it is often advantageous to describe the exchange splittings between all multiplets by a single exchange coupling constant J_AB within the Heisenberg spin Hamiltonian. In order to determine J_AB, it is sufficient to perform a single spin-flip calculation starting from the high-spin, single-determinant reference state.²³

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.²⁶ and Hapka et al.²⁷ introduced the symmetry-adapted perturbation theory based on, respectively, restricted and unrestricted Slater determinants (with Hartree-Fock or Kohn-Sham orbitals). Present implementations of SAPT, however, are valid only for the single-reference high-spin case of the dimer in which the spin quantum number of the complex is equal to the sum of spin quantum numbers of monomers (S = S_A + S_B).

The interaction of two open-shell species with the total spin of the dimer smaller than the sum of spin quantum numbers of the monomers still poses a challenge. This is the case since low-spin configurations are typical examples of multireference systems. The knowledge of the so-called exchange splitting, which is defined as the difference between the highest and lowest possible spin state of a given dimer, is very important in interactions of small open-shell molecules (e.g., in the O₂ ${}^{3}Σ_{g}^{-}$ dimer^28–30) as well as for interactions involving stable organic radicals such as phenalenyl³¹ or unsaturated metal sites within metal organic frameworks in the O₂ adsorption process.³² The new approaches of calculating low-spin potential energy surfaces could be applied to model fragments of interaction potentials in reaction dynamics, for instance, near the channels which correspond to dissociation into two radicals or near the bond-breaking geometry of the reactive complex. Spin splitting also plays an important role in cold physics since it is the term that strongly couples hyperfine states of alkali-metal atoms and drives Feshbach resonances.^33,34 It should be stressed that in the SAPT approach exchange splitting does indeed come exclusively from exchange corrections: the SAPT terms arising from polarization theory²⁴ 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 ${\hat{S}}^{2}$ ⁠. The theory introduced in this paper is beneficial for several reasons: (i) as usually in SAPT, the interaction energy components are calculated directly, which contrasts with supermolecular calculations that always involve subtraction; (ii) there is no need for the active space selection for the complex: in fact, no multireference calculation is necessary (similar to the spin-flip electronic structure approaches); (iii) it is possible to calculate the energy in the monomer-centered basis set, and the overall scaling of the method [once transformed molecular-orbital (MO) integrals are available] is just o⁴, 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,²⁶ intramonomer correlation effects are not included except for the possible use of Kohn-Sham orbitals.³⁵ 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.¹⁸ However, while in the latter theories the z component of the total spin of the system (the M_S quantum number) is changed upon the spin flip, in SF-SAPT only the $M_{S_{A}}$ and $M_{S_{B}}$ numbers for individual molecules are changed: the overall M_S value is conserved.

To date, there are only a few applications of perturbation theory to low-spin complexes. In 1984, Wormer and van der Avoird used the Heitler-London formula with appropriately spin-adapted linear combinations of Slater determinants (in the specific case of 4 active electrons) to predict the splittings between different multiplets of the O₂⋯O₂ system.³⁶ 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.³⁷ 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 $H_{2}^{+}$ 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 $H_{2}^{+}$ and discussed a new approach to exchange energy via a variational volume-integral formula.

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

II. THEORY

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

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

E_{i n t} = \frac{⟨ Ψ_{A} Ψ_{B} | V A P_{S M_{S}} | Φ_{A B} ⟩}{⟨ Ψ_{A} Ψ_{B} | A P_{S M_{S}} | Φ_{A B} ⟩},

(1)

in powers of the intermolecular interaction operator $V = \sum_{i \in A} \sum_{j \in B} \tilde{v} (i j)$ and separating the polarization and exchange effects in each order. In Eq. (1), $A$ is the (N_A + N_B)-electron antisymmetrizer and Φ_AB is the polarization (RS) expansion of the wavefunction for the complex.²⁴ We will assume that Ψ_A and Ψ_B are ground-state high-spin open-shell determinants. The occupied orbitals in Ψ_A are $χ_{i_{1}}$ ⁠, $χ_{i_{2}}$ ⁠, …, $χ_{i_{k_{A}}}$ (inactive, doubly occupied) and $χ_{a_{1}}$ ⁠, $χ_{a_{2}}$ ⁠, …, $χ_{a_{2 S_{A}}}$ (active, occupied by α electrons only). The occupied orbitals in Ψ_B are $χ_{j_{1}}$ ⁠, $χ_{j_{2}}$ ⁠, …, $χ_{j_{k_{B}}}$ (inactive, doubly occupied) and $χ_{b_{1}}$ ⁠, $χ_{b_{2}}$ ⁠, …, $χ_{b_{2 S_{B}}}$ (active, occupied by β electrons only).

The spin projector $P_{S M_{S}}$ commutes with V and the zeroth-order Hamiltonian H₀ = H_A + H_B as the latter operators are spin-independent. Moreover, the operator $Ŝ^{2}$ ⁠, and thus also $P_{S M_{S}}$ ⁠, 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,

E_{i n t} = \frac{⟨ P_{S M_{S}} Ψ_{A} Ψ_{B} | V A | Φ_{A B}^{'} ⟩}{⟨ P_{S M_{S}} Ψ_{A} Ψ_{B} | A | Φ_{A B}^{'} ⟩} .

(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 $P_{S M_{S}} Ψ_{A} Ψ_{B}$ ⁠, with all resulting perturbation corrections (in the expansion of $Φ_{A B}^{'}$ ⁠) 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 $P_{S M_{S}}$ reduces the zeroth-order space of (2S_A + 1)(2S_B + 1) asymptotically degenerate product functions corresponding to a pair of interacting multiplets to a single spin-adapted combination with the requested (S, M_S). For Eq. (1), any of the (2S_A + 1)(2S_B + 1) product functions, including Ψ_AΨ_B, is a valid zeroth-order state for nondegenerate RS perturbation theory as neither V nor H₀ mix states with different $M_{S_{A}}$ or $M_{S_{B}}$ ⁠. We choose Eq. (1) as the starting point for the derivation below.

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

A \sim (1 + P) A_{A} A_{B},

(3)

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

P = - \sum_{i \in A} \sum_{j \in B} P_{i j}

(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|X|Ψ_AΨ_B⟩, we obtain (note that $A$ ⁠, $A_{A}$ ⁠, and $A_{B}$ all commute with $P_{S M_{S}}$ ⁠)

\begin{matrix} (⟨ V ⟩ + E_{e x c h}^{(10)}) (⟨ P_{S M_{S}} ⟩ + ⟨ P P_{S M_{S}} ⟩) = ⟨ V P_{S M_{S}} ⟩ + ⟨ V P P_{S M_{S}} ⟩ . \end{matrix}

(5)

Note that we have separated the electrostatic contribution $E_{e l s t}^{(10)} = ⟨ V ⟩$ from the exchange one, and we introduced the zero into the E⁽¹⁰⁾ 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 $P_{S M_{S}}$ on the product Ψ_AΨ_B. The ROHF determinants Ψ_A and Ψ_B are pure spin functions corresponding to the quantum numbers $(S_{A}, M_{S_{A}} = S_{A})$ and $(S_{B}, M_{S_{B}} = - S_{B})$ ⁠, respectively. In the zeroth-order space spanned by (2S_A + 1)(2S_B + 1) product functions that correspond to total spin S_A for molecule A and total spin S_B for molecule B (and any combination of $M_{S_{A}}, M_{S_{B}}$ ⁠), there exists exactly one function $Ψ_{S, M_{S}}$ with a total spin S ∈ {|S_A − S_B|, …, S_A + S_B} and its projection M_S ∈ {−S, …, S}. This function is a linear combination of only those product functions that correspond to $M_{S_{A}} + M_{S_{B}} = M_{S}$ ⁠. Therefore, the projection $P_{S M_{S}} Ψ_{A} Ψ_{B}$ picks, up to normalization, this very function $Ψ_{S, M_{S}}$ out of the zeroth-order space. In other words, $P_{S M_{S}} Ψ_{A} Ψ_{B}$ produces a linear combination of products of functions of A and B with the same values S_A, S_B but different $M_{S_{A}}, M_{S_{B}}$ (however, the latter two numbers add up to the same total M_S). The coefficients in this linear combination are the Clebsch-Gordan coefficients $⟨ S M_{S} | S_{A} M_{S_{A}} S_{B} M_{S_{B}} ⟩$ —we will denote the coefficients ⟨S (S_A − S_B)|S_A S_A S_B − S_B⟩, ⟨S (S_A − S_B)|S_A (S_A − 1) S_B (−S_B + 1)⟩, and ⟨S (S_A − S_B)|S_A (S_A − 2) S_B (−S_B + 2)⟩ by c₀, c₁, and c₂, respectively. Specifically, the spin projection can be expressed as follows:

P_{S M_{S}} Ψ_{A} Ψ_{B} = c_{0} Ψ_{A} Ψ_{B} + c_{1} Ψ_{A}^{↓} Ψ_{B}^{↑} + c_{2} Ψ_{A}^{↓ ↓} Ψ_{B}^{↑ ↑} + \dots,

(6)

where, for example, $Ψ_{A}^{↓}$ is a normalized wavefunction (linear combination of determinants) corresponding to the spin quantum numbers (S_A, S_A − 1). Up to a constant, this function can be obtained from Ψ_A by the action of the spin-lowering operator $Ŝ_{-}$ ⁠.

Equation (6) implies that

⟨ P_{S M_{S}} ⟩ = c_{0},

(7)

⟨ V P_{S M_{S}} ⟩ = c_{0} ⟨ V ⟩ .

(8)

Therefore, Eq. (5) becomes

⟨ V ⟩ ⟨ P P_{S M_{S}} ⟩ + c_{0} E_{e x c h}^{(10)} + E_{e x c h}^{(10)} ⟨ P P_{S M_{S}} ⟩ = ⟨ V P P_{S M_{S}} ⟩ .

(9)

Now, the last term on the l.h.s. is neglected as it requires at least two electron exchanges (one in $E_{e x c h}^{(10)}$ and one in $⟨ P P_{S M_{S}} ⟩$ ⁠), and we arrive at the final formula for the SRS-like first-order exchange energy,

E_{e x c h}^{(10)} = c_{0}^{- 1} (⟨ V P P_{S M_{S}} ⟩ - ⟨ V ⟩ ⟨ P P_{S M_{S}} ⟩) .

(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 $Ŝ_{-} = \sum_{k = 1}^{N_{A}} Ŝ_{-} (k)$ applies the conventional spin-lowering operator $Ŝ_{-}$ (k) = $Ŝ_{x}$ (k) − $i Ŝ_{y}$ (k) to each of the N_A spins. It can be shown that when Ψ_A is a restricted high-spin determinant like in this work, the lowering operator acts on it as follows:

\begin{matrix} Ŝ_{-} Ψ_{A} & = & \sum_{n = 1}^{2 S_{A}} |χ_{i_{1}} α χ_{i_{1}} β \dots χ_{i_{k_{A}}} α χ_{i_{k_{A}}} β χ_{a_{1}} α \dots) \\ (\times χ_{a_{n - 1}} α χ_{a_{n}} β χ_{a_{n + 1}} α \dots χ_{a_{2 S_{A}}} α|, \end{matrix}

(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 2S_A determinants entering the linear combination are clearly orthonormal. Therefore, $Ψ_{A}^{↓} = (1 / \sqrt{2 S_{A}}) Ŝ_{-} Ψ_{A}$ is normalized.

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

\begin{matrix} E_{e x c h}^{(10)} & = & ⟨ V P ⟩ + c_{1} / c_{0} ⟨ Ψ_{A} Ψ_{B} | V P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩ \\ + c_{2} / c_{0} ⟨ Ψ_{A} Ψ_{B} | V P | Ψ_{A}^{↓ ↓} Ψ_{B}^{↑ ↑} ⟩ + \dots + \\ - ⟨ V ⟩ (⟨ P ⟩ + c_{1} / c_{0} ⟨ Ψ_{A} Ψ_{B} | P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩) \\ + (c_{2} / c_{0} ⟨ Ψ_{A} Ψ_{B} | P | Ψ_{A}^{↓ ↓} Ψ_{B}^{↑ ↑} ⟩ + \dots) . \end{matrix}

(12)

The first interesting observation from Eq. (12) is that the standard, closed-shell like contribution (the “spin-diagonal” term) $⟨ V P ⟩ - ⟨ V ⟩ ⟨ P ⟩$ 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.,⁴⁷

\begin{matrix} ⟨ V P ⟩ - ⟨ V ⟩ ⟨ P ⟩ & = & - [{\tilde{v}}_{α β}^{β α} + S_{α^{'}}^{β} ({\tilde{v}}_{α β}^{α α^{'}} - {\tilde{v}}_{α β}^{α^{'} α}) + S_{β^{'}}^{α} ({\tilde{v}}_{α β}^{β^{'} β} - {\tilde{v}}_{α β}^{β β^{'}})) \\ - (S_{α^{'}}^{β} S_{β^{'}}^{α^{'}} {\tilde{v}}_{α β}^{α β^{'}} - S_{α^{'}}^{β^{'}} S_{β^{'}}^{α} {\tilde{v}}_{α β}^{α^{'} β} + S_{α^{'}}^{β} S_{β^{'}}^{α} {\tilde{v}}_{α β}^{α^{'} β^{'}}] . \end{matrix}

(13)

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

\begin{matrix} {\tilde{v}}_{μ ν}^{λ κ} = (λ μ | κ ν) + {(v_{A})}_{ν}^{κ} S_{μ}^{λ} / N_{A} + {(v_{B})}_{μ}^{λ} S_{ν}^{κ} / N_{B} + V_{0} S_{μ}^{λ} S_{ν}^{κ} / N_{A} N_{B}, \end{matrix}

(14)

$(λ μ | κ ν)$ is a two-electron integral in the (11|22) notation, ${(v_{X})}_{ν}^{κ}$ is the nuclear attraction integral with all nuclei of X = A,B, $V_{0}$ 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

\begin{matrix} E_{e x c h, d i a g}^{(10)} = ⟨ V P ⟩ - ⟨ V ⟩ ⟨ P ⟩ & = & - {\tilde{v}}_{a j}^{j a} - {\tilde{v}}_{i b}^{b i} - 2 {\tilde{v}}_{i j}^{j i} + {\tilde{v}}_{a j}^{a^{'} a} S_{a^{'}}^{j} - {\tilde{v}}_{a j}^{a a^{'}} S_{a^{'}}^{j} + {\tilde{v}}_{a j}^{j j^{'}} S_{j^{'}}^{a} - 2 {\tilde{v}}_{a j}^{j^{'} j} S_{j^{'}}^{a} + {\tilde{v}}_{i j}^{a i} S_{a}^{j} - 2 {\tilde{v}}_{i j}^{i a} S_{a}^{j} - {\tilde{v}}_{a b}^{j b} S_{j}^{a} + {\tilde{v}}_{i b}^{b b^{'}} S_{b^{'}}^{i} \\ - {\tilde{v}}_{i b}^{b^{'} b} S_{b^{'}}^{i} + {\tilde{v}}_{i b}^{i^{'} i} S_{i^{'}}^{b} - 2 {\tilde{v}}_{i b}^{i i^{'}} S_{i^{'}}^{b} + {\tilde{v}}_{i j}^{j b} S_{b}^{i} - 2 {\tilde{v}}_{i j}^{b j} S_{b}^{i} - {\tilde{v}}_{a b}^{a i} S_{i}^{b} + {\tilde{v}}_{a j}^{i a} S_{i}^{j} + {\tilde{v}}_{i b}^{b j} S_{j}^{i} - 2 {\tilde{v}}_{a j}^{a i} S_{i}^{j} - 2 {\tilde{v}}_{i b}^{j b} S_{j}^{i} \\ + 2 {\tilde{v}}_{i j}^{i^{'} i} S_{i^{'}}^{j} + 2 {\tilde{v}}_{i j}^{j j^{'}} S_{j^{'}}^{i} - 4 {\tilde{v}}_{i j}^{i i^{'}} S_{i^{'}}^{j} - 4 {\tilde{v}}_{i j}^{j^{'} j} S_{j^{'}}^{i} + {\tilde{v}}_{a b}^{a^{'} b} S_{a^{'}}^{j} S_{j}^{a} + {\tilde{v}}_{a b}^{a b^{'}} S_{i}^{b} S_{b^{'}}^{i} + 2 {\tilde{v}}_{i b}^{a b} S_{a}^{j} S_{j}^{i} + 2 {\tilde{v}}_{a j}^{a b} S_{i}^{j} S_{b}^{i} \\ + 2 {\tilde{v}}_{a j}^{a^{'} j} S_{a^{'}}^{j^{'}} S_{j^{'}}^{a} + 2 {\tilde{v}}_{a j}^{a j^{'}} S_{i}^{j} S_{j^{'}}^{i} + {\tilde{v}}_{a j}^{a j^{'}} S_{a^{'}}^{j} S_{j^{'}}^{a^{'}} - {\tilde{v}}_{a j}^{a^{'} j^{'}} S_{a^{'}}^{j} S_{j^{'}}^{a} + 2 {\tilde{v}}_{i b}^{i b^{'}} S_{i^{'}}^{b} S_{b^{'}}^{i^{'}} + 2 {\tilde{v}}_{i b}^{i^{'} b} S_{i^{'}}^{j} S_{j}^{i} + {\tilde{v}}_{i b}^{i^{'} b} S_{i^{'}}^{b^{'}} S_{b^{'}}^{i} \\ - {\tilde{v}}_{i b}^{i^{'} b^{'}} S_{i^{'}}^{b} S_{b^{'}}^{i} + 4 {\tilde{v}}_{i j}^{a j} S_{a}^{j^{'}} S_{j^{'}}^{i} - 2 {\tilde{v}}_{i j}^{a j^{'}} S_{a}^{j} S_{j^{'}}^{i} + 4 {\tilde{v}}_{i j}^{i b} S_{i^{'}}^{j} S_{b}^{i^{'}} - 2 {\tilde{v}}_{i j}^{i^{'} b} S_{i^{'}}^{j} S_{b}^{i} + 2 {\tilde{v}}_{i j}^{i^{'} j} S_{i^{'}}^{b} S_{b}^{i} + 4 {\tilde{v}}_{i j}^{i^{'} j} S_{i^{'}}^{j^{'}} S_{j^{'}}^{i} \\ + 2 {\tilde{v}}_{i j}^{i j^{'}} S_{a}^{j} S_{j^{'}}^{a} + 4 {\tilde{v}}_{i j}^{i j^{'}} S_{i^{'}}^{j} S_{j^{'}}^{i^{'}} - 2 {\tilde{v}}_{i j}^{i^{'} j^{'}} S_{i^{'}}^{j} S_{j^{'}}^{i}, \end{matrix}

(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 ${\tilde{v}}_{a b}^{b a}$ 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} | V P | Ψ_{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} | V P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩ = \int \tilde{v} (12) ρ_{int}^{↓ ↑} (12) d τ_{12},

(16)

⟨ Ψ_{A} Ψ_{B} | P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩ = \frac{1}{N_{A} N_{B}} \int ρ_{int}^{↓ ↑} (12) d τ_{12},

(17)

where

ρ_{int}^{↓ ↑} (12) = N_{A} N_{B} \int Ψ_{A} Ψ_{B} P Ψ_{A}^{↓} Ψ_{B}^{↑} d τ_{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,

\begin{matrix} ρ_{int}^{↓ ↑} (12) & = & - ρ_{A}^{↓} (1 | 2) ρ_{B}^{↑} (2 | 1) - \int ρ_{A}^{↓} (1 | 4) Γ_{B}^{↑} (24 | 21) d τ_{4} \\ - \int Γ_{A}^{↓} (13 | 12) ρ_{B}^{↑} (2 | 3) d τ_{3} \\ - \int Γ_{A}^{↓} (13 | 14) Γ_{B}^{↑} (24 | 23) d τ_{3} d τ_{4} . \end{matrix}

(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^{'}) = N_{A} \int Ψ_{A}^{*} (1,2, \dots, N_{A}) Ψ_{A}^{↓} (1^{'}, 2, \dots, N_{A}) d τ_{1}^{'},

(20)

\begin{matrix} Γ_{A}^{↓} (12 | 1^{'} 2^{'}) & = & N_{A} (N_{A} - 1) \int Ψ_{A}^{*} (1,2,3, \dots, N_{A}) \\ \times Ψ_{A}^{↓} (1^{'}, 2^{'}, 3, \dots, N_{A}) d τ_{12}^{'} . \end{matrix}

(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^{'}) = \frac{1}{\sqrt{2 S_{A}}} \sum_{n = 1}^{2 S_{A}} χ_{a_{n}} α (1) χ_{a_{n}} β (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

\begin{matrix} Γ_{A}^{↓} (12 | 1^{'} 2^{'}) & = & \frac{1}{\sqrt{2 S_{A}}} \sum_{n = 1}^{2 S_{A}} [\sum_{m = 1}^{k_{A}} (χ_{a_{n}} α (1) χ_{i_{m}} α (2) χ_{a_{n}} β (1^{'}) χ_{i_{m}} α (2^{'}) - χ_{i_{m}} α (1) χ_{a_{n}} α (2) χ_{a_{n}} β (1^{'}) χ_{i_{m}} α (2^{'}))) \\ ((- χ_{a_{n}} α (1) χ_{i_{m}} α (2) χ_{i_{m}} α (1^{'}) χ_{a_{n}} β (2^{'}) + χ_{i_{m}} α (1) χ_{a_{n}} α (2) χ_{i_{m}} α (1^{'}) χ_{a_{n}} β (2^{'}))) \\ (+ \sum_{m = 1}^{k_{A}} (χ_{a_{n}} α (1) χ_{i_{m}} β (2) χ_{a_{n}} β (1^{'}) χ_{i_{m}} β (2^{'}) - χ_{i_{m}} β (1) χ_{a_{n}} α (2) χ_{a_{n}} β (1^{'}) χ_{i_{m}} β (2^{'}))) \\ ((- χ_{a_{n}} α (1) χ_{i_{m}} β (2) χ_{i_{m}} β (1^{'}) χ_{a_{n}} β (2^{'}) + χ_{i_{m}} β (1) χ_{a_{n}} α (2) χ_{i_{m}} β (1^{'}) χ_{a_{n}} β (2^{'}))) \\ (+ \sum_{m = 1}^{2 S_{A}} (χ_{a_{n}} α (1) χ_{a_{m}} α (2) χ_{a_{n}} β (1^{'}) χ_{a_{m}} α (2^{'}) - χ_{a_{m}} α (1) χ_{a_{n}} α (2) χ_{a_{n}} β (1^{'}) χ_{a_{m}} α (2^{'}))) \\ ((- χ_{a_{n}} α (1) χ_{a_{m}} α (2) χ_{a_{m}} α (1^{'}) χ_{a_{n}} β (2^{'}) + χ_{a_{m}} α (1) χ_{a_{n}} α (2) χ_{a_{m}} α (1^{'}) χ_{a_{n}} β (2^{'}))] . \end{matrix}

(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:

\begin{matrix} ⟨ Ψ_{A} Ψ_{B} | V P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩ & = & \frac{1}{2 \sqrt{S_{A} S_{B}}} [- {\tilde{v}}_{a b}^{b a}) + {\tilde{v}}_{a b}^{b j} S_{j}^{a} + {\tilde{v}}_{a j}^{j b} S_{b}^{a} - 2 {\tilde{v}}_{a j}^{b j} S_{b}^{a} + {\tilde{v}}_{a b}^{b b^{'}} S_{b^{'}}^{a} - {\tilde{v}}_{a b^{'}}^{b b^{'}} S_{b}^{a} + {\tilde{v}}_{a b}^{i a} S_{i}^{b} + {\tilde{v}}_{i b}^{a i} S_{a}^{b} - 2 {\tilde{v}}_{i a}^{i b} S_{b}^{a} + {\tilde{v}}_{a b}^{a^{'} a} S_{a^{'}}^{b} \\ - {\tilde{v}}_{a^{'} a}^{a^{'} b} S_{b}^{a} - 2 {\tilde{v}}_{a b}^{i j} S_{i}^{b} S_{j}^{a} - 2 {\tilde{v}}_{a j}^{i b} S_{i}^{j} S_{b}^{a} + 4 {\tilde{v}}_{a j}^{i j} S_{i}^{b} S_{b}^{a} + 4 {\tilde{v}}_{i b}^{i j} S_{a}^{b} S_{j}^{a} - 4 {\tilde{v}}_{i j}^{i j} S_{a}^{b} S_{b}^{a} + 2 {\tilde{v}}_{i b}^{i b^{'}} S_{a}^{b} S_{b^{'}}^{a} - 2 {\tilde{v}}_{a b}^{i b^{'}} S_{i}^{b} S_{b^{'}}^{a} \\ - 2 {\tilde{v}}_{i b^{'}}^{i b^{'}} S_{a}^{b} S_{b}^{a} + 2 {\tilde{v}}_{i b^{'}}^{a b^{'}} S_{a}^{b} S_{b}^{i} + 2 {\tilde{v}}_{a j}^{a^{'} j} S_{a^{'}}^{b} S_{b}^{a} - 2 {\tilde{v}}_{a b}^{a^{'} j} S_{a^{'}}^{b} S_{j}^{a} - 2 {\tilde{v}}_{a^{'} j}^{a^{'} j} S_{a}^{b} S_{b}^{a} + 2 {\tilde{v}}_{a^{'} b}^{a^{'} j} S_{a}^{b} S_{j}^{a} + {\tilde{v}}_{a^{'} b}^{a^{'} b^{'}} S_{a}^{b} S_{b^{'}}^{a} \\ (- {\tilde{v}}_{a^{'} b^{'}}^{a^{'} b^{'}} S_{a}^{b} S_{b}^{a} - {\tilde{v}}_{a b}^{a^{'} b^{'}} S_{a^{'}}^{b} S_{b^{'}}^{a} + {\tilde{v}}_{a b^{'}}^{a^{'} b^{'}} S_{a^{'}}^{b} S_{b}^{a}] . \end{matrix}

(24)

For the renormalization term, Eq. (17), we find [for example, by reanalyzing Eq. (24), replacing each $\tilde{v}$ by $\frac{1}{N_{A} N_{B}}$ times the appropriate product of overlap integrals] that

⟨ Ψ_{A} Ψ_{B} | P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩ = - \frac{1}{2 \sqrt{S_{A} S_{B}}} S_{a}^{b} S_{b}^{a} .

(25)

The ROHF electrostatic energy is equal to

E_{e l s t}^{(10)} = ⟨ V ⟩ = 4 {\tilde{v}}_{i j}^{i j} + 2 {\tilde{v}}_{a j}^{a j} + 2 {\tilde{v}}_{i b}^{i b} + {\tilde{v}}_{a b}^{a b} .

(26)

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

\begin{matrix} \frac{1}{2 \sqrt{S_{A} S_{B}}} E_{e x c h, 1 fl i p}^{(10)} & = & ⟨ Ψ_{A} Ψ_{B} | V P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩ - ⟨ V ⟩ ⟨ Ψ_{A} Ψ_{B} | P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩ \\ = & \frac{1}{2 \sqrt{S_{A} S_{B}}} [- {\tilde{v}}_{a b}^{b a}) + {\tilde{v}}_{a b}^{b j} S_{j}^{a} + {\tilde{v}}_{a j}^{j b} S_{b}^{a} - 2 {\tilde{v}}_{a j}^{b j} S_{b}^{a} + {\tilde{v}}_{a b}^{b b^{'}} S_{b^{'}}^{a} - {\tilde{v}}_{a b^{'}}^{b b^{'}} S_{b}^{a} + {\tilde{v}}_{a b}^{i a} S_{i}^{b} + {\tilde{v}}_{i b}^{a i} S_{a}^{b} - 2 {\tilde{v}}_{i a}^{i b} S_{b}^{a} + {\tilde{v}}_{a b}^{a^{'} a} S_{a^{'}}^{b} \\ - {\tilde{v}}_{a^{'} a}^{a^{'} b} S_{b}^{a} - 2 {\tilde{v}}_{a b}^{i j} S_{i}^{b} S_{j}^{a} - 2 {\tilde{v}}_{a j}^{i b} S_{i}^{j} S_{b}^{a} + 4 {\tilde{v}}_{a j}^{i j} S_{i}^{b} S_{b}^{a} + 4 {\tilde{v}}_{i b}^{i j} S_{a}^{b} S_{j}^{a} + 2 {\tilde{v}}_{i b}^{i b^{'}} S_{a}^{b} S_{b^{'}}^{a} - 2 {\tilde{v}}_{a b}^{i b^{'}} S_{i}^{b} S_{b^{'}}^{a} + 2 {\tilde{v}}_{i b^{'}}^{a b^{'}} S_{a}^{b} S_{b}^{i} \\ (+ 2 {\tilde{v}}_{a j}^{a^{'} j} S_{a^{'}}^{b} S_{b}^{a} - 2 {\tilde{v}}_{a b}^{a^{'} j} S_{a^{'}}^{b} S_{j}^{a} + 2 {\tilde{v}}_{a^{'} b}^{a^{'} j} S_{a}^{b} S_{j}^{a} + {\tilde{v}}_{a^{'} b}^{a^{'} b^{'}} S_{a}^{b} S_{b^{'}}^{a} - {\tilde{v}}_{a b}^{a^{'} b^{'}} S_{a^{'}}^{b} S_{b^{'}}^{a} + {\tilde{v}}_{a b^{'}}^{a^{'} b^{'}} S_{a^{'}}^{b} S_{b}^{a}] . \end{matrix}

(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} | V P | Ψ_{A}^{↓ ↓} Ψ_{B}^{↑ ↑} ⟩ = 0$ ⁠, where $Ψ_{A}^{↓ ↓}$ and $Ψ_{B}^{↑ ↑}$ are proportional to $Ŝ_{-} Ŝ_{-}$ Ψ_A and $Ŝ_{+} Ŝ_{+}$ Ψ_B, respectively. Similar to the previous term,

⟨ Ψ_{A} Ψ_{B} | V P | Ψ_{A}^{↓ ↓} Ψ_{B}^{↑ ↑} ⟩ = \int \tilde{v} (12) ρ_{i n t}^{↓ ↓ ↑ ↑} (12) d τ_{12},

(28)

where the double spin-flip interaction density matrix $ρ_{i n t}^{↓ ↓ ↑ ↑} (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

\begin{matrix} Γ_{A}^{↓ ↓} (12 | 1^{'} 2^{'}) & \sim & \sum_{m, n = 1}^{2 S_{A}} [χ_{a_{n}} α (1) χ_{a_{m}} α (2) χ_{a_{n}} β (1^{'}) χ_{a_{m}} β (2^{'}) - χ_{a_{m}} α (1) χ_{a_{n}} α (2) χ_{a_{n}} β (1^{'}) χ_{a_{m}} β (2^{'})) \\ (- χ_{a_{n}} α (1) χ_{a_{m}} α (2) χ_{a_{m}} β (1^{'}) χ_{a_{n}} β (2^{'}) + χ_{a_{m}} α (1) χ_{a_{n}} α (2) χ_{a_{m}} β (1^{'}) χ_{a_{n}} β (2^{'})] . \end{matrix}

(30)

According to Eq. (29), only the last term in the expression for $ρ_{i n t}^{↓ ↓ ↑ ↑}$ survives,

ρ_{i n t}^{↓ ↓ ↑ ↑} (12) = - \int Γ_{A}^{↓ ↓} (13 | 14) Γ_{B}^{↑ ↑} (24 | 23) d τ_{3} d τ_{4} .

(31)

Such an expression vanishes upon the integration over dτ₁ and dτ₂, either alone or with a spin-independent operator such as $\tilde{v}$ (12) [as in Eq. (28)]. In particular, in the integration over dτ₁, 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 S_A and S_B, in an arbitrary dimer spin state |S_A − S_B| ≤ S ≤ S_A + S_B, is given by the formula which combines the orbital expressions for $E_{e x c h, d i a g}^{(10)}$ and $E_{e x c h, 1 fl i p}^{(10)}$ defined in Eqs. (15) and (27), respectively. By simplifying the $c_{1} / (2 c_{0} \sqrt{S_{A} S_{B}})$ coefficient, as outlined in the Appendix, we obtain

\begin{matrix} E_{e x c h}^{(10)} & = & ⟨ V P ⟩ - ⟨ V ⟩ ⟨ P ⟩ + c_{1} / c_{0} [⟨ Ψ_{A} Ψ_{B} | V P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩ - ⟨ V ⟩ ⟨ Ψ_{A} Ψ_{B} | P | Ψ_{A}^{↓} Ψ_{B}^{↑} ⟩] \\ = & E_{e x c h, d i a g}^{(10)} + Z (S_{A}, S_{B}, S) E_{e x c h, 1 fl i p}^{(10)}, \end{matrix}

(32)

where

\begin{matrix} Z (S_{A}, S_{B}, S) = \frac{S (S + 1) + 2 S_{A} S_{B} - S_{A} (S_{A} + 1) - S_{B} (S_{B} + 1)}{4 S_{A} S_{B}} . \end{matrix}

(33)

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

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

In the following, we will denote AO indices by capital letters assuming, for simplicity, that the same AO basis set has been used to expand molecular orbitals of A and B (according to the so-called dimer-centered basis set formalism⁵²). 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, $C_{i α}^{K} = C_{i β}^{K}$ ⁠. We will use boldface letters for matrices and denote by $A \cdot B = \sum_{K L} A_{K L} B_{K L}$ the scalar product of two matrices. The inactive and active parts of density matrices for monomer A are $P_{K L}^{i A} = C_{i α}^{K} C_{i α}^{L} = C_{i β}^{K} C_{i β}^{L}$ ⁠, $P_{K L}^{a A} = C_{a α}^{K} C_{a α}^{L}$ ⁠, and similarly for monomer B (note that the active electrons on A all have spin α, the active electrons on B all have spin β). Therefore, the total density matrix is P^A = 2P^iA + P^aA. Now, the generalized Coulomb and exchange matrices are defined as

J {[X]}_{K L} = (K L | M N) X_{M N},

(34)

K {[X]}_{K L} = (K M | N L) X_{M N}

(35)

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

ω_{A} = v_{A} + 2 J^{i A} + J^{a A},

(36)

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

h_{A}^{α} = v_{A} + 2 J^{i A} + J^{a A} - K^{i A} - K^{a A},

(37)

h_{A}^{β} = v_{A} + 2 J^{i A} + J^{a A} - K^{i A},

(38)

h_{B}^{α} = v_{B} + 2 J^{i B} + J^{a B} - K^{i B},

(39)

h_{B}^{β} = v_{B} + 2 J^{i B} + J^{a B} - K^{i B} - K^{a B} .

(40)

With this notation, the diagonal and single spin-flip terms of the ROHF-based $E_{e x c h}^{(10)}$ formula, Eqs. (15) and (27), respectively, can be rewritten as follows:

\begin{matrix} E_{e x c h, d i a g}^{(10)} & = & - P^{i B} \cdot (2 K^{i A} + K^{a A}) - P^{a B} \cdot K^{i A} - (P^{i A} S^{A O} P^{i B}) \cdot (h_{A}^{α} + h_{A}^{β} + h_{B}^{α} + h_{B}^{β}) - (P^{a A} S^{A O} P^{i B}) \cdot (h_{A}^{α} + h_{B}^{α}) \\ - (P^{i A} S^{A O} P^{a B}) \cdot (h_{A}^{β} + h_{B}^{β}) + 2 (P^{i B} S^{A O} P^{i A} S^{A O} P^{a B}) \cdot ω_{A} + 2 (P^{i B} S^{A O} P^{i A} S^{A O} P^{i B}) \cdot ω_{A} \\ + (P^{a B} S^{A O} P^{i A} S^{A O} P^{a B}) \cdot ω_{A} + (P^{i B} S^{A O} P^{a A} S^{A O} P^{i B}) \cdot ω_{A} + 2 (P^{i A} S^{A O} P^{i B} S^{A O} P^{i A}) \cdot ω_{B} \\ + 2 (P^{i A} S^{A O} P^{i B} S^{A O} P^{a A}) \cdot ω_{B} + (P^{i A} S^{A O} P^{a B} S^{A O} P^{i A}) \cdot ω_{B} + (P^{a A} S^{A O} P^{i B} S^{A O} P^{a A}) \cdot ω_{B} \\ - 2 (P^{i A} S^{A O} P^{i B}) \cdot K [P^{i A} S^{A O} P^{i B}] - 2 (P^{a A} S^{A O} P^{i B}) \cdot K [P^{i A} S^{A O} P^{i B}] - 2 (P^{i A} S^{A O} P^{a B}) \cdot K [P^{i A} S^{A O} P^{i B}] \\ - (P^{a A} S^{A O} P^{i B}) \cdot K [P^{a A} S^{A O} P^{i B}] - (P^{i A} S^{A O} P^{a B}) \cdot K [P^{i A} S^{A O} P^{a B}], \end{matrix}

(41)

\begin{matrix} E_{e x c h, 1 fl i p}^{(10)} & = & - P^{a B} \cdot K^{a A} - (P^{a A} S^{A O} P^{a B}) \cdot (h_{A}^{α} + h_{B}^{β}) + (P^{a A} S^{A O} P^{i B}) \cdot K^{a B} + (P^{i A} S^{A O} P^{a B}) \cdot K^{a A} + 2 (P^{i B} S^{A O} P^{a A} S^{A O} P^{a B}) \cdot ω_{A} \\ + (P^{a B} S^{A O} P^{a A} S^{A O} P^{a B}) \cdot ω_{A} + 2 (P^{i A} S^{A O} P^{a B} S^{A O} P^{a A}) \cdot ω_{B} + (P^{a A} S^{A O} P^{a B} S^{A O} P^{a A}) \cdot ω_{B} \\ - 2 (P^{a A} S^{A O} P^{a B}) \cdot K [P^{i A} S^{A O} P^{i B}] - 2 (P^{a A} S^{A O} P^{a B}) \cdot K [P^{i A} S^{A O} P^{a B}] - 2 (P^{a A} S^{A O} P^{i B}) \cdot K [P^{a A} S^{A O} P^{a B}] \\ - 2 (P^{a A} S^{A O} P^{i B}) \cdot K [P^{i A} S^{A O} P^{a B}] - (P^{a A} S^{A O} P^{a B}) \cdot K [P^{a A} S^{A O} P^{a B}], \end{matrix}

(42)

where S^AO is the overlap matrix in the AO basis.

III. COMPUTATIONAL DETAILS

The MO-based formulas for the arbitrary-spin first-order exchange energy, Eqs. (15) and (27), have been implemented in two versions to verify the numerical correctness of the codes: into the ROHF-based SAPT code of Ref. 26 forming a part of the SAPT2012 package⁵³ and into the development version of the psi4 package⁵¹ using the straightforward psi4numpy framework⁵⁴ 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 $E_{e x c h, d i a g}^{(10)}$ and $E_{e x c h, 1 fl i p}^{(10)}$ [Eqs. (15) and (27), respectively] add up exactly to the formula for $E_{e x c h}^{(10)}$ in the high-spin (triplet) state in accordance with Eq. (33) for S_A = S_B = 1/2 and S = 1. The AO-based first-order exchange expressions, Eqs. (41) and (42), have been implemented 2into the development version of psi4, making use of its efficient evaluation of generalized Coulomb and exchange matrices, both with and without the DF approximation. Whenever the latter approximation was applied, the standard aug-cc-pVXZ/JKFIT set⁵⁵ was used as the auxiliary basis accompanying the orbital set aug-cc-pVXZ.⁵⁶ As shown in Sec. IV E, the DF approximation to $E_{e x c h}^{(10)}$ performs very well, and the density fitted AO-based code can be applied to much larger systems than the ones studied here. Unless stated otherwise, we used the aug-cc-pVTZ basis set.

IV. NUMERICAL RESULTS

In this section, we study the applications of theory developed in Sec. II to several test systems. For the first one, the Li⋯H complex, it is possible to compare the spin-flip theory with exact SAPT calculations of Patkowski et al.³⁸ Then, we focus on other atom-atom complexes including the lithium dimer, Li⋯N, and N⋯N systems. We also examine the Mn⋯Mn system, very extensively studied in the past, for which the spin-exchange splitting is surprisingly low. Finally, we also show a test of exchange splittings for molecular cases: the O₂(⁠ ${}^{3}Σ_{g}^{-}$ ⁠) 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.⁵⁷ The exchange splitting ΔE between the interaction energies for the highest (E_int(S_A + S_B)) and lowest (E_int(|S_A − S_B|)) multiplicities from ab initio calculations can be compared with the spin-flip SAPT term, Eq. (27), by multiplying the latter by a following factor:

Δ E^{(10)} = (1 + \frac{1}{2 max (S_{A}, S_{B})}) E_{e x c h, 1 fl i p}^{(10)},

(43)

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

A. Li⋯H

The first system we investigate is the Li⋯H complex, where the interaction between two doublets gives rise to a triplet state and a singlet state. For both of them, full configuration interaction (FCI) interaction energies as well as SAPT corrections to high order (in several variants including SRS) have been obtained before.³⁸ 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 $E_{e x c h}^{(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.

View large Download slide

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.

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 $E_{e x c h, d i a g}^{(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 $E_{e l s t}^{(10)}$ [Eq. (26)] up to terms that vanish as the fourth or higher powers of intermolecular overlap integrals.⁵⁹ The apparent Coulomb energy does not have this property and the difference between the two (that is, $E_{e x c h, d i a g}^{(10)}$ ⁠) does not vanish in the S² approximation. The only exceptions are the interactions involving one- and two-electron systems, such as H⋯H,³⁷ He (³S)⋯H,⁶⁰ and He (³S)⋯He (³S),⁶¹ when the mathematical Coulomb energy coincides with the mean energy of physical states. In these cases, $E_{e x c h, d i a g}^{(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 $E_{e x c h, d i a g}^{(10)}$ term is much smaller than $E_{e x c h, 1 fl i p}^{(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, $E_{e x c h}^{(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 $E_{e l s t}^{(10)}$ that is nearly negligible for this dispersion-bound complex. As a result, the sum $E_{e l s t}^{(10)} + E_{e x c h, d i a g}^{(10)}$ recovers only a small fraction (8% at the triplet van der Waals minimum distance of 11.5a₀) of the FCI apparent Coulomb energy, while $E_{e x c h, 1 fl i p}^{(10)}$ recovers 71% of the FCI exchange energy.

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

The total spin for the Li⋯Li complex can be either singlet or triplet, for Li⋯N—triplet or quintet, and for two interacting nitrogen atoms it can take multiplicities from singlet to septet. For these systems, unlike for the Li⋯H interaction, the benchmark SAPT corrections with exact monomer wavefunctions are not available. In Fig. 2, we show the comparison of exchange splittings from SF-SAPT, CASSCF, and MRCI (Davidson corrected). As expected, the first-order result fails to recover the exchange splitting around the chemical minima, which are at 5a₀ for Li⋯Li, 3.5a₀ for Li⋯N, and 2.1a₀ 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.9a₀ for the lithium dimer,⁶² 10.2a₀ 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.2a₀ for the nitrogen dimer.⁶³ 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 $E_{e l s t}^{(10)}$ and $E_{e x c h, d i a g}^{(10)}$ ⁠, show very little basis set dependence: our Li⋯Li tests at R = 7.9a₀ 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.

View large Download slide

Comparison of the ΔE⁽¹⁰⁾ 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.⁶⁴ were also plotted.

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.⁶⁴ 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 S² approximation by inspecting the ratio of $E_{e x c h}^{(10)} (S^{2})$ and $E_{e x c h}^{(10)}$ for the high-spin state (the nonapproximated values were computed using the high-spin code of Ref. 26). From Fig. 3 it is clear that while for the lithium dimer the single exchange approximation is quite drastic, even in the van der Waals minimum region, it works very well for the nitrogen dimer and the Li⋯N system. Such behavior can be attributed to a very small ionization potential I_p of the Li atom (4.2 eV) compared to the N atom (14.5 eV) which directly affects the decay rates of the wavefunctions [at long range, this decay rate is proportional to $\exp (- \sqrt{2 I_{p}} r)$ ⁶⁵].

FIG. 3.

View large Download slide

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 $E_{e x c h}^{(10)} (S^{2}) / E_{e x c h}^{(10)}$ ⁠. Color-coded vertical dashed lines correspond to the positions of minima for the high-spin complexes.

C. O₂⋯O₂

Oxygen dimer is one of the best studied van der Waals complexes with relevance to the chemistry of atmosphere.²⁸ Hence, this was one of the very first systems for which the exchange interaction was studied. Wormer and van der Avoird calculated the exchange splitting within the Heisenberg model and the J_AB parameter was obtained with variation-perturbation theory.³⁶ 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.²⁹ 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.28a₀).

In Fig. 4 we compare the splitting between the highest and lowest spin states of the O₂⋯O₂ complex for four basic angular geometries: H-shape, linear, T-shape, and X-shape with previous studies of Wormer and van der Avoird³⁶ and Bartolomei and co-workers,²⁹ 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.5a₀ but it is not clearly seen in Fig. 4. This node can, however, be seen on the logarithmic plot of the diagonal and spin-flip parts of the exchange splitting, Fig. 5, as a dip in the absolute value of the splitting. Note that the diagonal part of the exchange energy is very large in the oxygen dimer for all configurations including X-shape, almost two orders of magnitude larger than the spin-flip part. Such a big domination of the diagonal part is responsible for the fact that the oxygen molecules strongly repel each other at short range and form only weakly bound complexes for all spin states of the dimer. Finally, let us remark that the single exchange approximation works very well for the oxygen dimer, as shown by the comparison of full and S²-approximated $E_{e x c h}^{(10)}$ values for the high-spin (quintet) complex displayed in Fig. 6. This behavior is consistent with the observation made by Wormer and van der Avoird³⁶ that the single-parameter Heisenberg model recovers the quintet-triplet-singlet splittings accurately.

FIG. 4.

View large Download slide

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 (³ $Σ_{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. 5.

View large Download slide

Diagonal and spin-flip parts of the first-order exchange energy for the dimer of ground-state oxygen molecules (³ $Σ_{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 7a₀ corresponding to a change of sign from positive (for small R) to negative.

FIG. 6.

View large Download slide

The quality of the S² approximation for the high-spin quintet state of the molecular oxygen dimer. The oxygen molecules are in their ground ³ $Σ_{g}^{-}$ state and their mutual orientation corresponds to four basic geometric configurations.

D. Mn⋯Mn

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 4s² shell and the resulting spin exchange splitting is very small, on the order of 10 cm⁻¹ in the minimum of the potential energy curve, which is 2 orders of magnitude less than the binding energy. The first ab initio study of exchange interaction in the Mn₂ system was initiated by Nesbet⁶⁶ who used the Heisenberg model and estimated the J_AB parameter to be small (up to 62 K at an interatomic separation of 4.5a₀). 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.⁷¹ Clearly, a development of new methods for systems similar to the manganese dimer is warranted.

The ROHF method produces correct orbitals for the Mn atom (with a correct degeneracy of d orbitals) and can be straightforwardly applied to first-order SF-SAPT. In Fig. 7, the spin-flip SAPT exchange term is very small, nearly two orders of magnitude smaller than the diagonal exchange energy. Similar to the oxygen dimer, this causes a very small exchange splitting and very small differences between the potential well depths and equilibrium distances for different multiplets. It is also worthwhile to inspect the quality of the S² approximation by comparing the high-spin exchange energies. Near the van der Waals minimum of undecaplet Mn⋯Mn, the single exchange approximation reproduces about 90% of the full exchange. In order to assess the performance of first-order SF-SAPT, we have computed the J_AB parameter and compared it with previous literature data in panel (c) of Fig. 7. In the literature, the J_AB parameter is usually given at the minimum separation for the high-spin complex (in the case of the work of Buchachenko et al.,⁷² 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.⁷³ (although we realize this might be somewhat fortuitous). Except for the value by Negodaev et al.,⁶⁹ all presented values of J_AB are smaller than the result derived from SF-SAPT.

FIG. 7.

View large Download slide

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 S² approximation for the high-spin (undecaplet) Mn⋯Mn state; (c) comparison of the Heisenberg J_AB parameters derived from SF-SAPT with existing theoretical and experimental data; note that only in the study of Buchachenko et al.⁷² the J_AB 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 J_AB is marked corresponds to the CCSD(T) ¹ $Σ_{g}^{+}$ minimum obtained in Ref. 72.

E. Phenalenyl dimer

The dimer of the doublet phenalenyl radical is an example of “pancake bonding,” a strong interaction between π-stacked radicals that has gathered significant interest in recent years.⁷⁴ 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.³¹ At the interplanar separation of 3.104 Å corresponding to the pancake bonded minimum depicted in Fig. 8, the reference singlet-triplet splitting, computed using the high-level multireference averaged quadratic coupled cluster (MR-AQCC) approach⁷⁵ 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.⁷⁶

FIG. 8.

View large Download slide

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.

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 $E_{e x c h}^{(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 $E_{e l s t}^{(10)}$ ⁠, $E_{e x c h, d i a g}^{(10)}$ ⁠, and $E_{e x c h, 1 fl i p}^{(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⁽¹⁰⁾ 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 $E_{e x c h, d i a g}^{(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,⁵¹ any improvements made to these core psi4 objects will also be extended to the SF-SAPT calculations.

V. SUMMARY

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

We investigated the behavior of the diagonal and spin-flip components of the first-order SAPT exchange correction on a number of interatomic and intermolecular complexes. For the Li⋯H system, we compared the exchange energies with existing FCI-based SAPT calculations³⁸ 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 $E_{e x c h}^{(10)}$ for the highest and lowest spin state to the splittings obtained with the CASSCF and MRCI methods. The $E_{e x c h}^{(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 $E_{e x c h}^{(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.²⁹ A particularly interesting test of our new theory was the manganese dimer, which poses a great challenge for supermolecular calculations with a variety of multireference methods. The Heisenberg exchange coupling constant J_AB derived from $E_{e x c h}^{(10)}$ compares qualitatively well with literature data and predicts small spin splitting due to screening from the outermost doubly occupied (4s²) shell. However, it should be noted that for such a challenging system the overall agreement between existing supermolecular data is far from satisfactory. Interestingly, first-order spin-flip SAPT predicts the value of J_AB within the experimental bounds. Finally, a highly efficient AO-based implementation (with density fitting) of $E_{e x c h}^{(10)}$ allowed us to apply SF-SAPT to a larger system, the phenalenyl dimer. In that case, we obtained a qualitative agreement with reference CASSCF calculations. The agreement improves at larger separations which might be attributed to the importance of terms beyond the S² approximation. The new first-order exchange corrections have been implemented into the development version of the psi4 code⁵¹ as well as the sapt2012 package.⁵³ 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.⁷⁷ The ideas presented here can be generalized to the second-order exchange-induction and exchange-dispersion corrections as well as to a multireference, CASSCF-based description of the noninteracting wavefunctions. Moreover, while the single exchange approximation implies that double- and higher-spin-flip terms vanish, it is a different approximation than merely neglecting multiple spin flips. In fact, going beyond the S² approach (following the ideas of Refs. 44–46) might be useful even for two interacting doublets at short range (such as in the phenalenyl dimer) where there is only a single active spin on each monomer to be flipped. The work in all of these directions is in progress in our groups.

ACKNOWLEDGMENTS

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

APPENDIX: PROOF OF EQ. (33)

In this appendix, we calculate the coefficient $Z (S_{A}, S_{B}, S) = \frac{c_{1}}{2 c_{0} \sqrt{S_{A} S_{B}}}$ that appears when inserting Eq. (27) into Eq. (12), proving Eq. (33). Recall that

c_{0} = ⟨ S (S_{A} - S_{B}) | S_{A} S_{A} S_{B} - S_{B} ⟩,

(A1)

c_{1} = ⟨ S (S_{A} - S_{B}) | S_{A} (S_{A} - 1) S_{B} (- S_{B} + 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 |S_A S_A S_B − S_B⟩ produces

\begin{matrix} Ŝ_{A -} Ŝ_{B +} | S_{A} S_{A} S_{B} - S_{B} ⟩) & = & \sqrt{S_{A} (S_{A} + 1) - S_{A} (S_{A} - 1)} \sqrt{S_{B} (S_{B} + 1) + S_{B} (1 - S_{B})} | S_{A} (S_{A} - 1) S_{B} (- S_{B} + 1) ⟩) \\ = & 2 \sqrt{S_{A} S_{B}} | S_{A} (S_{A} - 1) S_{B} (- S_{B} + 1) ⟩) . \end{matrix}

(A3)

Since $Ŝ^{2} = Ŝ_{A}^{2} + Ŝ_{B}^{2} + 2 {(Ŝ_{A})}_{z} {(Ŝ_{B})}_{z} + Ŝ_{A -} Ŝ_{B +} + Ŝ_{A +} Ŝ_{B -}$ ⁠, one might rewrite the action of the operator $Ŝ_{A -} Ŝ_{B +}$ as

\begin{matrix} Ŝ_{A -} Ŝ_{B +} | S_{A} S_{A} S_{B} - S_{B} ⟩) & = & [Ŝ^{2} - (Ŝ_{A}^{2} + Ŝ_{B}^{2} + 2 {(Ŝ_{A})}_{z} {(Ŝ_{B})}_{z} + Ŝ_{A +} Ŝ_{B -})] | S_{A} S_{A} S_{B} - S_{B} ⟩) \\ = & [Ŝ^{2} - (S_{A} (S_{A} + 1) + S_{B} (S_{B} + 1) - 2 S_{A} S_{B} + 0)] | S_{A} S_{A} S_{B} - S_{B} ⟩) . \end{matrix}

(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 (S_A − S_B)| and using the fact that ⟨S (S_A − S_B)| $Ŝ^{2}$ = S(S + 1)⟨S (S_A − S_B)|, one finds that

\begin{matrix} 2 c_{1} \sqrt{S_{A} S_{B}} = [S (S + 1) - S_{A} (S_{A} + 1) - S_{B} (S_{B} + 1) + 2 S_{A} S_{B}] c_{0}, \end{matrix}

(A5)

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

REFERENCES

1.

S.

Aloisio

and

J. S.

Francisco

,

Acc. Chem. Res.

33

,

825

(

2000

).

https://doi.org/10.1021/ar000097u

Google Scholar

Crossref

PubMed

2.

A.

Galano

,

M.

Narciso-Lopez

, and

M.

Francisco-Marquez

,

J. Phys. Chem. A

114

,

5796

(

2010

).

https://doi.org/10.1021/jp101157b

Google Scholar

Crossref

PubMed

3.

E.

Roueff

and

F.

Lique

,

Chem. Rev.

113

,

8906

(

2013

).

https://doi.org/10.1021/cr400145a

Google Scholar

Crossref

PubMed

4.

I.

Ratera

and

J.

Veciana

,

Chem. Soc. Rev.

41

,

303

(

2012

).

https://doi.org/10.1039/c1cs15165g

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1126/science.1229549

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1038/nchem.1857

Google Scholar

Crossref

PubMed

7.

S.

Chefdeville

,

Y.

Kalugina

,

S. Y. T.

van de Meerakker

,

C.

Naulin

,

F.

Lique

, and

M.

Costes

,

Science

341

,

1094

(

2013

).

https://doi.org/10.1126/science.1241395

Google Scholar

Crossref

PubMed

8.

P. J.

Knowles

,

C.

Hampel

, and

H.-J.

Werner

,

J. Chem. Phys.

99

,

5219

(

1993

).

https://doi.org/10.1063/1.465990

Google Scholar

Crossref

9.

J. D.

Watts

,

J.

Gauss

, and

R. J.

Bartlett

,

J. Chem. Phys.

98

,

8718

(

1993

).

https://doi.org/10.1063/1.464480

Google Scholar

Crossref

10.

B.

Jeziorski

and

H. J.

Monkhorst

,

Phys. Rev. A

24

,

1668

(

1981

).

https://doi.org/10.1103/physreva.24.1668

Google Scholar

Crossref

11.

R. J.

Bartlett

and

M.

Musiał

,

Rev. Mod. Phys.

79

,

291

(

2007

).

https://doi.org/10.1103/revmodphys.79.291

Google Scholar

Crossref

12.

B.

Jeziorski

,

Mol. Phys.

108

,

3043

(

2010

).

https://doi.org/10.1080/00268976.2010.524169

Google Scholar

Crossref

13.

H.-J.

Werner

and

P. J.

Knowles

,

J. Chem. Phys.

89

,

5803

(

1988

).

https://doi.org/10.1063/1.455556

Google Scholar

Crossref

14.

K.

Andersson

,

P.-A.

Malmqvist

,

B.

Roos

,

A.

Sadlej

, and

K.

Wolinski

,

J. Phys. Chem.

94

,

5483

(

1990

).

https://doi.org/10.1021/j100377a012

Google Scholar

Crossref

15.

J.

Finley

,

P.

Malmqvist

,

B. O.

Roos

, and

L.

Serrano-Andrés

,

Chem. Phys. Lett.

288

,

299

(

1998

).

https://doi.org/10.1016/s0009-2614(98)00252-8

Google Scholar

Crossref

16.

C.

Angeli

,

R.

Cimiraglia

, and

J.-P.

Malrieu

,

J. Chem. Phys.

117

,

9138

(

2002

).

https://doi.org/10.1063/1.1515317

Google Scholar

Crossref

17.

C.

Camacho

,

R.

Cimiraglia

, and

H. A.

Witek

,

Phys. Chem. Chem. Phys.

12

,

5058

(

2010

).

https://doi.org/10.1039/c000196a

Google Scholar

Crossref

PubMed

18.

A. I.

Krylov

,

Chem. Phys. Lett.

338

,

375

(

2001

).

https://doi.org/10.1016/s0009-2614(01)00287-1

Google Scholar

Crossref

19.

Y.

Shao

,

M.

Head-Gordon

, and

A. I.

Krylov

,

J. Chem. Phys.

118

,

4807

(

2003

).

https://doi.org/10.1063/1.1545679

Google Scholar

Crossref

20.

S. V.

Levchenko

and

A. I.

Krylov

,

J. Chem. Phys.

120

,

175

(

2004

).

https://doi.org/10.1063/1.1630018

Google Scholar

Crossref

PubMed

21.

A. I.

Krylov

,

Acc. Chem. Res.

39

,

83

(

2006

).

https://doi.org/10.1021/ar0402006

Google Scholar

Crossref

PubMed

22.

P. M.

Zimmerman

,

F.

Bell

,

M.

Goldey

,

A. T.

Bell

, and

M.

Head-Gordon

,

J. Chem. Phys.

137

,

164110

(

2012

).

https://doi.org/10.1063/1.4759076

Google Scholar

Crossref

PubMed

23.

N. J.

Mayhall

and

M.

Head-Gordon

,

J. Phys. Chem. Lett.

6

,

1982

(

2015

).

https://doi.org/10.1021/acs.jpclett.5b00733

Google Scholar

Crossref

PubMed

24.

B.

Jeziorski

,

R.

Moszyński

, and

K.

Szalewicz

,

Chem. Rev.

94

,

1887

(

1994

).

https://doi.org/10.1021/cr00031a008

Google Scholar

Crossref

25.

E. G.

Hohenstein

and

C. D.

Sherrill

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

2

,

304

(

2012

).

https://doi.org/10.1002/wcms.84

Google Scholar

Crossref

26.

P. S.

Żuchowski

,

R.

Podeszwa

,

R.

Moszyński

,

B.

Jeziorski

, and

K.

Szalewicz

,

J. Chem. Phys.

129

,

084101

(

2008

).

https://doi.org/10.1063/1.2968556

Google Scholar

Crossref

PubMed

27.

M.

Hapka

,

P. S.

Żuchowski

,

M. M.

Szczȩśniak

, and

G.

Chałasiński

,

J. Chem. Phys.

137

,

164104

(

2012

).

https://doi.org/10.1063/1.4758455

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1021/ja9917215

Google Scholar

Crossref

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

).

https://doi.org/10.1039/b803555e

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1063/1.3479395

Google Scholar

Crossref

PubMed

31.

Z.

Cui

,

H.

Lischka

,

H. Z.

Beneberu

, and

M.

Kertesz

,

J. Am. Chem. Soc.

136

,

5539

(

2014

).

https://doi.org/10.1021/ja412862n

Google Scholar

Crossref

PubMed

32.

M. V.

Parkes

,

D. F.

Sava Gallis

,

J. A.

Greathouse

, and

T. M.

Nenoff

,

J. Phys. Chem. C

119

,

6556

(

2015

).

https://doi.org/10.1021/jp511789g

Google Scholar

Crossref

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

).

https://doi.org/10.1103/physrevlett.94.103201

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1103/physreva.70.032701

Google Scholar

Crossref

35.

H. L.

Williams

and

C. F.

Chabalowski

,

J. Phys. Chem. A

105

,

646

(

2001

).

https://doi.org/10.1021/jp003883p

Google Scholar

Crossref

36.

P. E. S.

Wormer

and

A.

van der Avoird

,

J. Chem. Phys.

81

,

1929

(

1984

).

https://doi.org/10.1063/1.447867

Google Scholar

Crossref

37.

T.

Ćwiok

,

B.

Jeziorski

,

W.

Kołos

,

R.

Moszyński

, and

K.

Szalewicz

,

J. Chem. Phys.

97

,

7555

(

1992

).

https://doi.org/10.1063/1.463475

Google Scholar

Crossref

38.

K.

Patkowski

,

T.

Korona

, and

B.

Jeziorski

,

J. Chem. Phys.

115

,

1137

(

2001

).

https://doi.org/10.1063/1.1379330

Google Scholar

Crossref

39.

K.

Patkowski

,

B.

Jeziorski

,

T.

Korona

, and

K.

Szalewicz

,

J. Chem. Phys.

117

,

5124

(

2002

).

https://doi.org/10.1063/1.1499488

Google Scholar

Crossref

40.

P.

Gniewek

and

B.

Jeziorski

,

Phys. Rev. A

90

,

022506

(

2014

).

https://doi.org/10.1103/physreva.90.022506

Google Scholar

Crossref

41.

P.

Gniewek

and

B.

Jeziorski

,

J. Chem. Phys.

143

,

154106

(

2015

).

https://doi.org/10.1063/1.4931809

Google Scholar

Crossref

PubMed

42.

P.

Gniewek

and

B.

Jeziorski

,

Phys. Rev. A

94

,

042708

(

2016

).

https://doi.org/10.1103/physreva.94.042708

Google Scholar

Crossref

43.

B.

Jeziorski

,

G.

Chałasiński

, and

K.

Szalewicz

,

Int. J. Quantum Chem.

14

,

271

(

1978

).

https://doi.org/10.1002/qua.560140306

Google Scholar

Crossref

44.

B.

Jeziorski

,

M.

Bulski

, and

L.

Piela

,

Int. J. Quantum Chem.

10

,

281

(

1976

).

https://doi.org/10.1002/qua.560100208

Google Scholar

Crossref

45.

R.

Schäffer

and

G.

Jansen

,

Theor. Chem. Acc.

131

,

1235

(

2012

).

https://doi.org/10.1007/s00214-012-1235-6

Google Scholar

Crossref

46.

R.

Schäffer

and

G.

Jansen

,

Mol. Phys.

111

,

2570

(

2013

).

https://doi.org/10.1080/00268976.2013.827253

Google Scholar

Crossref

47.

R.

Moszyński

,

B.

Jeziorski

,

S.

Rybak

,

K.

Szalewicz

, and

H. L.

Williams

,

J. Chem. Phys.

100

,

5080

(

1994

).

https://doi.org/10.1063/1.467225

Google Scholar

Crossref

48.

K.

Patkowski

,

K.

Szalewicz

, and

B.

Jeziorski

,

J. Chem. Phys.

125

,

154107

(

2006

).

https://doi.org/10.1063/1.2358353

Google Scholar

Crossref

PubMed

49.

F. A.

Matsen

,

D. J.

Klein

, and

D. C.

Foyt

,

J. Phys. Chem.

75

,

1866

(

1971

).

https://doi.org/10.1021/j100681a017

Google Scholar

Crossref

50.

A.

Hesselmann

,

G.

Jansen

, and

M.

Schütz

,

J. Chem. Phys.

122

,

014103

(

2005

).

https://doi.org/10.1063/1.1824898

Google Scholar

Crossref

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

).

https://doi.org/10.1021/acs.jctc.7b00174

Google Scholar

Crossref

PubMed

52.

H. L.

Williams

,

E. M.

Mas

,

K.

Szalewicz

, and

B.

Jeziorski

,

J. Chem. Phys.

103

,

7374

(

1995

).

https://doi.org/10.1063/1.470309

Google Scholar

Crossref

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.

Google Scholar

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

).

https://doi.org/10.1063/1.1445115

Google Scholar

Crossref

56.

R. A.

Kendall

,

T. H.

Dunning

, Jr., and

R. J.

Harrison

,

J. Chem. Phys.

96

,

6796

(

1992

).

https://doi.org/10.1063/1.462569

Google Scholar

Crossref

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

).

https://doi.org/10.1002/wcms.82

Google Scholar

Crossref

58.

W.

Kutzelnigg

,

J. Chem. Phys.

73

,

343

(

1980

).

https://doi.org/10.1063/1.439880

Google Scholar

Crossref

59.

B.

Jeziorski

, private communication (

2018

).

60.

M.

Przybytek

,

K.

Patkowski

, and

B.

Jeziorski

,

Collect. Czech. Chem. Commun.

69

,

141

(

2004

).

https://doi.org/10.1135/cccc20040141

Google Scholar

Crossref

61.

M.

Przybytek

and

B.

Jeziorski

,

J. Chem. Phys.

123

,

134315

(

2005

).

https://doi.org/10.1063/1.2042453

Google Scholar

Crossref

PubMed

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

).

https://doi.org/10.1103/physreva.87.052505

Google Scholar

Crossref

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

).

https://doi.org/10.1103/physreva.82.042718

Google Scholar

Crossref

64.

R.

Côté

,

A.

Dalgarno

, and

M. J.

Jamieson

,

Phys. Rev. A

50

,

399

(

1994

).

https://doi.org/10.1103/physreva.50.399

Google Scholar

Crossref

PubMed

65.

R.

Ahlrichs

,

M.

Hoffmann-Ostenhof

,

T.

Hoffmann-Ostenhof

, and

J. D.

Morgan

,

Phys. Rev. A

23

,

2106

(

1981

).

https://doi.org/10.1103/physreva.23.2106

Google Scholar

Crossref

66.

R. K.

Nesbet

,

Phys. Rev.

135

,

A460

(

1964

).

https://doi.org/10.1103/physrev.135.a460

Google Scholar

Crossref

67.

Z.

Wang

,

R. R.

Lucchese

, and

J. W.

Bevan

,

J. Phys. Chem. A

108

,

2884

(

2004

).

https://doi.org/10.1021/jp036603p

Google Scholar

Crossref

68.

S.

Yamamoto

,

H.

Tatewaki

,

H.

Moriyama

, and

H.

Nakano

,

J. Chem. Phys.

124

,

124302

(

2006

).

https://doi.org/10.1063/1.2178798

Google Scholar

Crossref

PubMed

69.

I.

Negodaev

,

C.

de Graaf

, and

R.

Caballol

,

Chem. Phys. Lett.

458

,

290

(

2008

).

https://doi.org/10.1016/j.cplett.2008.05.011

Google Scholar

Crossref

70.

D.

Tzeli

,

U.

Miranda

,

I. G.

Kaplan

, and

A.

Mavridis

,

J. Chem. Phys.

129

,

154310

(

2008

).

https://doi.org/10.1063/1.2993750

Google Scholar

Crossref

PubMed

71.

C.

Camacho

,

S.

Yamamoto

, and

H. A.

Witek

,

Phys. Chem. Chem. Phys.

10

,

5128

(

2008

).

https://doi.org/10.1039/b805125a

Google Scholar

Crossref

PubMed

72.

A. A.

Buchachenko

,

G.

Chałasiński

, and

M. M.

Szczȩśniak

,

J. Chem. Phys.

132

,

024312

(

2010

).

https://doi.org/10.1063/1.3292572

Google Scholar

Crossref

PubMed

73.

M.

Cheeseman

,

R. J.

Van Zee

,

H. L.

Flanagan

, and

W.

Weltner

,

J. Chem. Phys.

92

,

1553

(

1990

).

https://doi.org/10.1063/1.458086

Google Scholar

Crossref

74.

K. E.

Preuss

,

Polyhedron

79

,

1

(

2014

).

https://doi.org/10.1016/j.poly.2014.04.005

Google Scholar

Crossref

75.

P. G.

Szalay

and

R. J.

Bartlett

,

Chem. Phys. Lett.

214

,

481

(

1993

).

https://doi.org/10.1016/0009-2614(93)85670-j

Google Scholar

Crossref

76.

Z.

Mou

,

Y.-H.

Tian

, and

M.

Kertesz

,

Phys. Chem. Chem. Phys.

19

,

24761

(

2017

).

https://doi.org/10.1039/c7cp04637e

Google Scholar

Crossref

PubMed

77.

A. J.

Misquitta

and

K.

Szalewicz

,

J. Chem. Phys.

122

,

214109

(

2005

).

https://doi.org/10.1063/1.1924593

Google Scholar

Crossref

PubMed

2018

Author(s)

First-order symmetry-adapted perturbation theory for multiplet splittings

I. INTRODUCTION

II. THEORY

III. COMPUTATIONAL DETAILS

IV. NUMERICAL RESULTS

A. Li⋯H

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

C. O₂⋯O₂

D. Mn⋯Mn

E. Phenalenyl dimer

V. SUMMARY

ACKNOWLEDGMENTS

APPENDIX: PROOF OF EQ. (33)

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

First-order symmetry-adapted perturbation theory for multiplet splittings Free

I. INTRODUCTION

II. THEORY

III. COMPUTATIONAL DETAILS

IV. NUMERICAL RESULTS

A. Li⋯H

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

C. O2⋯O2

D. Mn⋯Mn

E. Phenalenyl dimer

V. SUMMARY

ACKNOWLEDGMENTS

APPENDIX: PROOF OF EQ. (33)

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

First-order symmetry-adapted perturbation theory for multiplet splittings

C. O₂⋯O₂