Density-fitted open-shell symmetry-adapted perturbation theory and application to π-stacking in benzene dimer cation and ionized DNA base pair steps

The quantification and understanding of non-covalent interactions underlies the work of both chemists and biochemists. Although weak, non-covalent interactions play decisive roles in numerous processes:¹ formation of pre-reactive complexes, determination of preferred orientation in enantioselective reactions,^2,3 and conformation of proteins⁴ and macromolecular assemblies.^5,6 Small energy differences in non-covalent interactions govern the most stable crystalline phase and hence the bulk properties of numerous materials.⁷ 

The most widespread method to quantitatively compute non-covalent interaction energies is the so-called supermolecular method, where the energy difference between the full system and its monomers is computed with a suitable level of theory. Some of the most common methods to access supermolecular interactions are Density Functional Theory (DFT)⁸ and Møller-Plesset perturbation theory (MP2)⁹ or its empirically tuned variants like Spin Component Scaled MP2 (SCS-MP2)¹⁰ or SCS-MP2 reparametrized for Molecular Interactions, SCS-MI-MP2.¹¹ Excellent accuracy is available at a higher computational cost through coupled cluster methods, among which Coupled Cluster with Singles, Doubles, and perturbative Triples, CCSD(T),^12,13 extrapolated at the complete basis set limit is considered the gold standard of quantum chemistry.

Even very accurate interaction energies yield only limited insight into the underlying physical phenomena. To alleviate this shortcoming, numerous interaction energy decomposition analysis methods have been devised. Some of them rely on stepwise relaxation of constrained monomer wavefunctions within a non-covalent dimer to obtain electrostatic, exchange, polarization, and charge transfer components.^14–21 Other methods rationalize interaction energies in terms of deviations, computed through perturbation theory, from an ideal Lewis structure reference;^22–25 treat the system as a collection of interacting fragments;^26–32 partition the system in real space;^33–36 or exploit more original partitions.^37–43 

As an alternative to the supermolecular method, perturbation theory obtains interaction energies directly from monomer wavefunctions, without computing the dimer explicitly. Perturbation approaches are free of the Basis-Set Superposition Error (BSSE) by construction and offer a different way to analyze interaction energies. In particular, Symmetry-Adapted Perturbation Theory (SAPT)^44,45 naturally yields a decomposition of the total interaction energy in terms of physically meaningful components such as electrostatics, exchange, induction, and dispersion. Moreover, SAPT is systematically improvable through higher order perturbations and can account for intramonomer correlation effects.

All methods mentioned above were primarily developed and applied to closed-shell systems. However, the interaction of radicals with either open- or closed-shell molecules bears considerable importance to a number of essential processes and devices such as photosynthesis, DNA oxidative damage, dye-sensitized solar cells, or organic electronics. As a consequence, open-shell extensions of energy decomposition methods have been devised for the most common schemes.^46–48 SAPT was also extended to the open-shell case as soon as 1980 for the $H 2 +$ cation;⁴⁹ however, a general many-electron implementation of the first-order terms did not become available until 1995.⁵⁰ Only recently was a full implementation including second-order terms published, based either on Restricted Open-shell Hartree-Fock (ROHF) or Kohn-Sham (ROKS)⁵¹ or on their Unrestricted counterparts, UHF and UKS⁵² for both SAPT based on Density Functional Theory, SAPT(DFT), and SAPT based on Hartree-Fock, SAPT(HF) [the latter, when using uncoupled dispersion energies,⁵³ is also referred to as SAPT0 because the intramolecular correlation is treated at zeroth-order]. Both papers mention that the density fitting approximation could be used to access large molecules, but the actual open-shell density-fitted implementation was never published or made available to our knowledge.

Hence, open-shell SAPT has only been applied to small systems up to now.^54–58 The largest open-shell SAPT0 computation to our knowledge was performed on radical cation homodimers of furan and thiophene.⁵⁹ However, the interaction of large radicals attracts more and more interest in the literature.^60–64 A crucial process generating large radicals is the oxidative damage of DNA, which can potentially influence both base pairing and base pair stacking.^65–69 In the neutral case, stacking is dispersion-dominated,^70,71 but one might expect the introduction of a charge to increase the overall binding energy and make it electrostatically dominated through cation-π interactions. How do physical components of π-π interactions evolve as the system becomes open-shell?

In the present work, we attempt to answer this question. First, we introduce an efficient density-fitted open-shell implementation of SAPT based on a UHF reference. We then demonstrate the code reliability by reproducing results obtained on small systems by Hapka et al.⁵² and examine the accuracy of the density fitting approximation. The behavior of SAPT0(UHF) on open-shell systems is examined with a prototypical example of a challenging π-π system, the benzene dimer cation. We consider here three parallel and one T-shaped configuration at various intermonomer distances. Finally, we demonstrate the efficiency of the present implementation by applying it to 10 different DNA base pair steps and their radical cations.

First, our open-shell SAPT0 implementation was tested by comparison with the literature results using a set of small van der Waals complexes.^51,52 For consistency with the reference results, we used the aug-cc-pVTZ basis set for CN⋯Ne and OH⋯He and the aug-cc-pVQZ basis set for He⋯NH and NH⋯NH. All basis sets were supplemented with sets of midbond functions as previously described.⁵¹ The geometries of the van der Waals complexes and the Molpro^72,73 reference energies to full precision were provided by Hapka.¹¹⁴ 

Except for the above-mentioned van der Waals complexes, all SAPT0 computations on both neutral and cationic complexes were carried out in the jun-cc-pVDZ basis set,⁷⁴ known to provide maximal error cancellation for closed-shell systems.⁷⁵ To assess the performance of SAPT0(UHF) on a particularly challenging case, we chose the benzene radical cations previously studied by Krylov and co-workers.^76,77 The benzene monomers are kept frozen at their neutral geometries in four dimer configurations: sandwich, edge-displaced (or x-displaced), vertex-displaced (or y-displaced, see Figure 3), and T-shaped (see Figure 2). The geometries were extracted from the available supplementary material in Ref. 77 and used without reoptimization for both the radical cation and neutral computations.

The reference interaction energies for the radical systems were computed with DF-FNO-EOM-IP-CCSD^78–80/cc-pVTZ⁸¹ in Q-Chem.⁸² To improve computational efficiency, the density fitting^83,84 and the FNO approximations^85–88 were applied, with a criterion of 99.99% of the virtual occupation conserved for FNO, so that the resulting error was below 0.05 kcal/mol (see the supplementary material Table IV). Density fitting employed the RI-cc-pVTZ⁸⁹ auxiliary basis set. No BSSE correction was applied as it was found to be detrimental to the accuracy of our results (see the supplementary material Section III).

The convergence of the unrestricted Hartree-Fock wavefunctions proved particularly challenging, especially for the benzene dimer radical cation, with several stable low-lying SCF solutions. This difficulty may originate in the large number of low-lying excited states in these systems.⁷⁷ Accordingly, we explored different initial guesses for the SCF procedure and retained the one yielding the lowest energy. In some cases, orbitals from nearby points on the potential energy surface were used as guesses for other geometries to ensure the global SCF minimum was achieved across all points. The stability⁹⁰ of all UHF wavefunctions was verified.

Ten different DNA base pair steps at the average B-DNA geometry were generated from the Nucleic Acids Database^91,92 by Parker et al.⁷¹ and obtained from the corresponding supplementary material.⁷¹ Computations on these base steps made full use of density fitting, with the jun-cc-pVDZ-JKFIT^93,112 auxiliary basis set for all generalized (see below) Coulomb and exchange matrices in SCF and SAPT0(UHF), and the jun-cc-pVDZ-RI^89,94,113 basis to fit products of occupied-virtual orbitals in the dispersion energy. The base pairs are denoted WX:YZ, where WX and YZ form π-π stacked pairs while X:Y and W:Z form H-bonded base pairs. The geometries were kept the same for neutral and radical cation computations, and the convergence and stability of the open-shell Hartree-Fock wavefunctions were checked as for the benzene dimer radical cation case.

In symmetry-adapted perturbation theory,^44,45 the dimer Hamiltonian $H ˆ$ is partitioned as

where $H ˆ A$ and $H ˆ B$ are the Hamiltonians of isolated monomers A and B, respectively. $F ˆ$ is the Fock operator and $W ˆ$ is the fluctuation operator representing the electronic correlation contribution to the energy of the monomers. In the UHF case, $F ˆ X$ is the sum of the α and β Fock operators. In this version of SAPT(UHF), the dimer wavefunction does not contain any spin coupling terms between the monomers, thus only the high-spin state of the dimer is accessible. Note that if $F ˆ X$ is instead written as an effective ROHF Fock matrix, all the equations of SAPT0(UHF) are still valid, with the UHF orbital energies being replaced by the eigenvalues of the ROHF Fock matrix.⁵² This means that SAPT0(UHF) is compatible with both UHF and ROHF reference wavefunctions. $V ˆ$ is the intermolecular interaction operator, containing the repulsion between nuclei, the electron-electron repulsion, and the electron-nuclear attraction between monomers. The zeroth-order Hamiltonian in SAPT is $H ˆ 0 = F ˆ A + F ˆ B$ and the Schrödinger equation is perturbatively expanded in orders of $V ˆ$⁠, $W ˆ A$⁠, and $W ˆ B$⁠. The resulting equations are then symmetry-adapted to properly include exchange effects at low order.⁴⁴ The dimer interaction energy is

where i, j, and k are perturbation orders, respectively, in $V ˆ$⁠, $W ˆ A$⁠, and $W ˆ B$⁠, and l = j + k is the total intramonomer correlation order. Truncation of the energy expansion at different orders defines different levels of SAPT. The present paper is concerned with the efficient implementation of open-shell SAPT0, which contains the following terms:⁹⁵ 

The consecutive terms on the right-hand side of Equation (3) represent the electrostatic, exchange, induction, exchange-induction, dispersion, exchange-dispersion, and Hartree-Fock correction to the interaction energy, respectively. The induction terms include orbital relaxation if the corresponding amplitudes are obtained through Coupled-Perturbed Hartree-Fock (CPHF) equations. Note that our current implementation of CPHF does not support ROHF reference orbitals, thus only uncoupled ROHF induction and exchange-induction can be computed at present. In the remainder of this paper, we present UHF-based results that include orbital relaxation for induction terms.

The Hartree-Fock correction δ_HF,r is computed from the supermolecular Hartree-Fock interaction energy $E int HF$ as

This correction ensures that SAPT0 recovers the Hartree-Fock interaction energy, which is reasonably accurate for electrostatics, exchange, and induction, in addition to a second-order correction for dispersion and exchange-dispersion. The detailed formulae for the UHF-based SAPT energy components were previously reported in the literature,^51,52 albeit with some typographical mistakes. Thus, the formulae were verified thoroughly before implementation and are reproduced in the supplementary material in a form convenient for direct implementation.

Density fitting was introduced in the open-shell SAPT0 energy expressions using similar techniques as for the most recent implementation of closed-shell SAPT0 expressions in Psi4.^96–98 The SAPT0 energy terms contain contractions of two-electron integrals with various quantities such as overlap matrices, density matrices, induction, or dispersion amplitudes. The contractions found in the electrostatic, exchange, induction, and exchange-induction terms can be conveniently expressed as generalized Coulomb and exchange matrices defined, respectively, as⁵² 

where A and B are general matrices, i runs over the molecular orbitals of monomer X, the two-electron integrals $v μ λ ν σ$ are

and the indices μ, ν, λ, and σ run over the AO basis functions. For example, $E exch ( 10 ) ( S 2 )$ can be written as

where $A ⋅ B$ denotes a scalar product between two matrices, S^AO is the atomic orbital overlap matrix, $P σ X = C σ X ( C σ X ) †$ with $C σ X$ containing occupied orbital coefficients, $K σ X$ is the usual exchange matrix, $ω X = V X + J α X + J β X$⁠, and $h σ X = ω X − K σ X$⁠. Within Psi4, a specialized object performs the evaluation of the generalized Coulomb and exchange matrices and handles efficiently both exact integrals and approximate density-fitted integrals. This feature allows us to write a single code for both the exact and the density-fitted open-shell SAPT terms that can be expressed through generalized Coulomb and exchange matrices. Orbital relaxation in the $E ind , r ( 20 )$ and $E exch − ind , r ( 20 )$ energies is obtained by solving the CPHF equations, which are also expressed through generalized Coulomb and exchange matrices and are thus available both in their exact and density-fitted versions. $E disp ( 20 )$ and $E exch − disp ( 20 )$ cannot be written entirely in terms of generalized Coulomb and exchange matrices; however, they benefit from a factorization written directly in terms of the density-fitted three-index tensors.⁹⁸ Accordingly, dispersion and exchange-dispersion can only be computed within the density-fitted approximation in the current implementation. All SAPT0 energy terms benefit from multi-threaded parallelization for multicore architectures.

The validity of the present implementation was assessed by comparison with published SAPT0(UHF) results from Molpro^72,73 for CN⋯Ne, OH⋯He, He⋯NH, and NH⋯NH.^51,52 Using exact integrals, we successfully reproduced the reference results for each individual energy term, with a maximum deviation of 4 ⋅ 10⁻⁶ kcal mol⁻¹ (see Table I in the supplementary material) and a maximum relative deviation of 0.0012%. The remaining difference can be attributed to minor variations between integral screening thresholds and SCF convergence criteria.

Dispersion and exchange-dispersion terms are excluded from the above consideration since our implementation of these terms always makes use of the density-fitting approximation, whereas Molpro’s implementation only uses exact integrals. The density fitting approximation is known to be very accurate for closed-shell SAPT,^99–101 and the observed maximum deviation with Molpro is less than 4 ⋅ 10⁻⁴ kcal mol⁻¹ for dispersion terms (see the supplementary material Table I), or a maximum relative deviation of 1.5%. To confirm the reliability of density-fitting in the open-shell case, we repeated SAPT0(UHF) computations on the 4 van der Waals complexes considered above with density fitting for all terms. We also computed exact and density-fitted SAPT0(UHF) energy components for the larger vertex-displaced benzene dimer radical cation at intermonomer separations d = 3.1, 5.0, and 20.0 Å (see Figure 3). Note that the Hartree-Fock monomer wavefunctions were density-fitted as well, thus the starting orbitals and orbital energies for DF-SAPT0(UHF) are slightly different from the exact computation. However, the maximum error is only 3.5 ⋅ 10⁻⁴ kcal mol⁻¹ (see the supplementary material Table II) and the maximum relative deviation 1.6% for the set of small molecules. Considering the larger benzene dimer radical cation, the maximum absolute error observed is 2.5 ⋅ 10⁻² kcal mol⁻¹, corresponding only to a relative error of 0.11%. The maximum relative deviation in the perturbation expansion is 0.5% for the exchange component at a distance d = 5.0 Å. At a distance d = 20.0 Å, the relative error on δ_HF,r is 46%, but this component is 3 orders of magnitude smaller than the total SAPT0 energy, which then deviates by less than 0.02% (see the supplementary material Table III).

To evaluate the performance of our multithreaded parallelized density-fitted code,⁹⁸ we timed a SAPT0/jun-cc-pVDZ computation for CA:TG on one of our compute nodes, equipped with an Intel i7-3930K processor with 6 compute cores clocked at 3.20 GHz and 64 Gb of DDR4 RAM. Psi4 was compiled on the node with Intel compiler version 15.0.3. The CA:TG computation was successively performed with 1, 2, 4, and 6 threads with no competing job on the node, both with exact integrals (except for the dispersion terms) and density fitting. Speed-ups were computed as the ratio of the wall time in seconds for N threads over the wall time in seconds for a single thread. Wall times only include the SAPT0(UHF) module and not the preceding Hartree-Fock computations (see Figure 1).

Note that the single example here merely provides qualitative information on the scaling, and a more significant, quantitative study of the SAPT0(UHF) performance is outside the scope of the present work. We observe that density fitting (DF) exhibits slightly lower speed-ups than the exact integrals but is still close to ideal behaviour when the number of threads increases. This demonstrates that our use of density fitting does not impede the single-node parallelism of our code. The slightly inferior efficiency of DF may be attributed to the reduced amount of work that can be distributed between threads.

The benzene dimer is a typical π-π stacking system, and its radical cation is prototypical of radical aromatic interactions. As such, it has been the subject of a number of previous studies,^76,102–104 among which the work of Krylov and co-workers is most relevant to us since they provided high-level EOM-IP-CCSD/6-31+G* binding energies at various geometries.^76,77 In the present work, we apply their methodology with a larger cc-pVTZ basis set and controlled approximations to speed up the computation, as summarized in Section II.

For the four configurations of benzene dimer studied, we scanned the potential energy surface along the intermonomer coordinate d. For the T-shaped configuration, d is the closest distance between the stem monomer and the geometrical center of the base monomer (see Figure 2), and for the parallel configurations, d is the distance between the two planes containing the monomers (see Figure 3). The optimal distance d_eq for these four configurations was taken from the literature.^76,77 The potential energy curves were obtained by shifting d_eq by −0.6, − 0.4, − 0.2, 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.2, and +1.7 Å.

Among the configurations studied, the sandwich and both parallel-displaced ones exhibit important charge delocalization between the monomers (see the supplementary material Tables VI–VIII).⁷⁶ SAPT obtains interaction energies perturbatively from monomer wavefunctions, where the charge is completely localized on one benzene molecule. The extended electronic structure reorganization induced by charge transfer is particularly challenging for perturbation theory, which fundamentally relies on a weak perturbation. The T-shaped configuration exhibits much less charge transfer (see the supplementary material Table IX), thus it offers us an opportunity to discern any intrinsic difficulties associated with treating radical interactions through SAPT. Since the charge is mostly localized on the stem, our SAPT computations start from a cationic stem and a neutral base, unless otherwise noted. Accordingly, we begin by examining the T-shaped configuration, before turning to the parallel ones. We will then discuss the physical components of the energy and their evolution compared to neutral benzene dimers.

For the T-shaped configuration, the SAPT0(UHF)/jun-cc-pVDZ interaction energies are quite accurate, being between +1.2 kcal mol⁻¹ and −1.5 kcal mol⁻¹ of the reference values, well within the bounds of the observed accuracy for closed-shell systems.⁷⁵ The potential energy curve is closely reproduced, with SAPT0(UHF) slightly overbinding at large distances.

The results for the three different parallel configurations are summarized in Figure 3. Around equilibrium, SAPT0(UHF) underestimates the binding energy by 5.0–7.0 kcal mol⁻¹, depending on the conformer. In spite of this difference, the equilibrium distance is reasonably reproduced, as well as the shape of the potential energy curve around equilibrium. At larger interfragment separation, SAPT0(UHF) transitions from underestimating to overestimating the binding energy, albeit by a much smaller absolute amount of about 0.7 kcal mol⁻¹. At large distance, the accuracy is thus similar to the case of closed-shell systems⁷⁵ and to the observed behavior for the T-shaped configuration.

The data presented above indicate that the SAPT0(UHF) error is larger for all parallel configurations, which exhibit significant charge transfer, than for the charge-localized T-shaped configuration. Moreover, the binding energy is underestimated around the intermonomer equilibrium distance and slightly overestimated at long range. The error measured is qualitatively correlated with the charge transfer occurring in the dimers. At short distance, the error decreases (see the supplementary material Section IV), which indicates that the terms missing in SAPT0(UHF) become less important in the repulsive region. These results suggest that SAPT0(UHF) is suitable to treat quantitatively open-shell systems with limited charge transfer and provides qualitative insight for other cases. A more extensive study of SAPT0(UHF) accuracy is underway in our laboratory.

In spite of the cationic character of the dimer configurations, dispersion forces are as important as electrostatics on most of the potential energy curves (see Figures 2 and 3). As expected, dispersion is slightly less important for the T-shaped complex; however, it still accounts for about 50% of the interaction energy at equilibrium. The high-order induction terms contained in the δ_HF,r correction represent the dominant part of the total induction energy. Induction is expected to be very significant in a charged complex, and the importance of the δ_HF,r correction for the interaction energy of ionic species has recently been highlighted.¹⁰⁵ We note that a large δ_HF,r may indicate that the underlying perturbation theory is starting to break down, which may also explain the quantitative deviations observed, especially for the parallel configurations. In the T-shaped configuration, switching the initial monomer charge site from the stem to the less favored base decreases attractive electrostatics and increases repulsive exchange interactions (see the supplementary material Section V for a more detailed discussion).

Including δ_HF,r into the SAPT0(UHF) energy makes it equivalent to the Hartree-Fock supermolecular interaction energy supplemented by a second-order dispersion correction. The quantitative error observed for parallel configurations has thus two possible origins: missing intermolecular terms, namely, high-order dispersion and coupled dispersion-induction; and missing intramonomer correlation effects. Since the monomer wavefunctions undergo significant electronic rearrangement upon charge transfer, intramonomer correlation effects are likely to play a significant role. Computation of these terms requires high-order SAPT formulae for open-shell monomers or a detailed open-shell SAPT(DFT) study which are outside the scope of the current work.

The present data on cationic benzene dimer configurations indicate that dispersion interactions are as important as electrostatics and are essential for a proper description of the system, in spite of its charge. For a semi-quantitative picture, high-order induction corrections through δ_HF,r are essential.

In Secs. IV B 1 and IV B 2, we examined the physical mechanisms responsible for binding in the benzene dimer cation and the performance of SAPT0(UHF) for this case. Here, we discuss the evolution of the SAPT0(UHF) energy components upon ionization of the neutral benzene dimer.

To facilitate visualization, we represent directly the difference in each SAPT0(UHF) energy component upon ionizing the benzene dimer to the corresponding radical cation, i.e., we plot $Δ E ( i l ) = E cat. ( i l ) − E neut. ( i l )$ where E^(il) is an energy component, cat. stands for the radical cation, and neut. stands for the neutral dimer. The results for the parallel configurations are gathered in Figure 4, and we begin our analysis by considering only those configurations.

Intuitively, we may expect the creation of an electric charge by removal of an electron to strengthen the intermolecular interaction, which now involves an ion. This is indeed the case for all dimer configurations examined here. Electrostatics is the main contributor to this change around equilibrium and at larger distances; however, at close range electrostatics becomes slightly less attractive upon ionization. Perhaps surprising at first, this effect may be attributed to the decreased charge penetration upon contraction of the electron density. The electron removal thus reduces close-range favorable electrostatic interactions with the partner’s partially unshielded nuclei.

Unexpectedly, second-order induction does not play a major role in the energetic changes observed. Higher-order induction terms, contained in the δ_HF,r correction, vary as much as second-order induction at large distance but become strongly dominant at equilibrium (distance d ≈ 3.3 − 3.1 Å) and shorter distances. The total induction change, including second- and higher-order effects, is as important as the change in electrostatics at long distances. The increasing variation in high-order induction at short range for parallel configurations is likely linked to charge transfer and strong polarization effects that become more important as monomers get closer (see the supplementary material Tables VI–VIII).

Exchange repulsion strongly decreases when removing an electron from neutral benzene dimer. The π clouds of the monomers reduce their spatial extent, which in turn decreases the Pauli repulsion for a fixed geometry. The decrease in exchange repulsion is directly related to the overlap between the monomers. Considering only the most compressed geometry, exchange repulsion is reduced by 15.0 kcal mol⁻¹ in the sandwich configuration, where the two monomers are directly stacked on top of each other. The edge-displaced parallel configuration exhibits less overlap between monomer orbitals, and accordingly exchange repulsion only varies by 10.5 kcal mol⁻¹. Finally, exchange is only reduced by 8.7 kcal mol⁻¹ in the vertex-displaced parallel dimer, where a vertex of the first monomer is almost on top of the second monomer center, yielding the least overlap.

Finally, in all cases the dispersion interaction becomes less attractive upon ionization, consistent with the decrease in polarizability of the cation. Overall, the effect of ionization on π-π stacking interaction in the case of the parallel configurations can be divided into two contributions. The first is the electrostatic one, and intuitively induces an increased electrostatic attraction at equilibrium and large distances, while enhancing attractive induction contributions (especially at high order) at all distances. The second contribution regroups effects linked to the number of electrons in the system. Removal of an electron decreases exchange repulsion, charge penetration, and dispersion attractions. The effect on exchange dominates over the dispersion at short distances, resulting in a net attractive effect of the ionization. At long distances, the r⁻⁶ dependence of dispersion makes it dominant over the exponentially decreasing exchange. However, the ionization is still overall attractive because of long-range electrostatics.

The T-shaped configuration binds through a CH⋯π interaction instead of π-π stacking. Figure 5 reveals whether the above observations on parallel configurations still hold. The electrostatic contributions exhibit the same overall behaviour, except that they do not become less attractive at the shortest distance explored in this work. The second- and higher-order inductions behave similarly to the parallel cases, albeit with changes of a smaller magnitude overall. In particular, the high-order induction terms contained in δ_HF,r do not exhibit strong changes at short distances, consistent with less charge transfer and a less drastic charge reorganization. Electrostatics and total induction are the most affected components upon ionization of the T-shaped configuration. The variation of the exchange repulsion upon ionization is much smaller than in the parallel case, reaching only 4 kcal mol⁻¹ at the shortest distance. This is consistent with the limited overlap of the two T-shaped monomers. Finally, dispersion becomes again less attractive as an electron is removed from the system. Thus, in spite of quantitative differences, the trends in the T-shaped benzene dimer are in agreement with the observations made on the parallel configurations.

We now turn to a more biologically relevant system, DNA. Many factors can compromise the structural integrity, hence the function of this biologically essential molecule.^106–108 Among them, oxidative damage^109,110 can result in the photo- or chemo-induced removal of an electron from DNA. The resulting hole can be mobile on the double strand¹¹¹ and localize on a hole acceptor site. Here, we study the influence of oxidation on the stacking interaction in 10 DNA base pair steps at their average B-DNA geometries, as reported in the work of Parker and Sherrill.⁷¹ In the present work, a base pair step WX:YZ denotes a tetramer of nucleobases, where a hydrogen-bonded Watson-Crick nucleobase pair X:Y interacts via π-π stacking with another nucleobase pair W:Z.

Our results are summarized in Figure 6 where the base steps have been grouped by nucleobase composition for clarity. Note that we report again the difference in each SAPT0(UHF) energy component upon ionization, i.e., $ΔE(il)=Ecat.(il)−Eneut.(il)$⁠. For each base step, the charge was attributed to the base pair with the lowest Hartree-Fock ionization energy in the jun-cc-pVDZ basis (see the supplementary material Table X), and all open-shell wavefunctions were checked for instabilities.

Upon electron removal, the stacking interaction between nucleobase pairs always becomes more attractive, by 3.1–8.6 kcal mol⁻¹, which is significant compared to the neutral base step interaction energies, ranging from −8.4 to −17.8 kcal mol⁻¹. The radical cationic base step interaction energies span a smaller range, from −14.5 to −20.9 kcal mol⁻¹.

Qualitatively, the ionization has the same effect on the energy components in the base steps as in the benzene dimers: the π electron clouds become more compact, thus dispersion is always less attractive and exchange less repulsive. Their sum, which we might call “net dispersion,” exhibits a small change of less than 1 kcal mol⁻¹ in the base steps, which can be repulsive or attractive. Electrostatics and induction become more attractive since a charge is created. These last two components are generally the most affected by ionization. The δ_HF,r correction to the energy is smaller than second-order induction, thus the perturbation expansion seems well-behaved and the energy components are expected to be reasonably reliable. Overall, ionization affects mostly electrostatics, making it significantly more attractive, while reducing the magnitude of dispersion interactions.

Our data outline general trends for the DNA base steps; however, some particular points do not fit with the general picture. For example, the component most affected by ionization for AC:GT is exchange and not electrostatics. A closer examination reveals that the neutral AC:GT base step exhibits the largest exchange repulsion of the 10 configurations studied (see the supplementary material Table XXII) and that the two base pairs have a very strong overlap. For AT:AT, ionization affects induction more than electrostatics, which may be caused by an increased polarizability of the base step compared to the other systems.

Turning to the energy components themselves instead of their variation, we see that electrostatics is more important than induction for both the radical cations (see Figure 7) and the neutral case. The only exception is the neutral GG:CC base step, where the two dipoles of the base pairs are aligned almost on top of each other, generating a repulsive electrostatics interaction that is dominated by the attractive induction component. In spite of the introduction of a charge, the cationic complexes are still clearly dispersion-dominated. This highlights again the importance of π-π stacking interactions even for ions, where electrostatics may intuitively seem more important.

In the present work, we have reported an efficient implementation of open-shell SAPT0 based on density-fitted electron repulsion integrals, using libraries provided by the Psi4 program package. The efficiency of this code on multicore architectures was demonstrated on the case of a CA:TG base pair step. The SAPT0(UHF) code will soon be released in the public version of Psi4.

The accuracy of SAPT0(UHF) for radical cation interactions was tested with the challenging benzene dimer radical cation, for which parallel configurations exhibit important charge transfer. In spite of this extensive charge reorganization, SAPT0(UHF) is still able to yield at least qualitatively correct results and allows us to investigate the relative importance of different energy components. The T-shaped configuration of the benzene dimer exhibits much less charge transfer and is more readily amenable to a perturbation theory treatment. Accordingly, the error relative to EOM-IP-CCSD reference interaction energies does not exceed the bounds of SAPT0 accuracy on closed-shell systems, and the energy components obtained are expected to be accurate. In all cases, and despite the creation of a charge, the dispersion interaction still represents a significant part of the benzene dimer radical cation interaction.

Ten B-DNA base pair steps were investigated in their neutral and radical cationic states to probe how π-stacking in nucleobases is affected by ionization. The neutral base pairs are bound mainly by dispersion, electrostatics generally being the next main attractive force. Upon ionization, the energy components are qualitatively affected in the same way as for the benzene dimer: exchange becomes less repulsive, dispersion less attractive, electrostatics and induction become more attractive. Electrostatics is generally most affected, followed closely by induction. However, the radical cation base steps are still mainly dispersion-bound.

The Psi4 implementation of SAPT0(UHF) presented here allows us to investigate large radical systems. Preliminary tests indicate that radicals are not inherently problematic for this perturbation theory and that SAPT0(UHF) accuracy does not seem to deteriorate compared to closed-shell systems. Charge transfer, or more generally extensive charge reorganization upon interaction, affects the quantitative aspects of SAPT0(UHF) but still preserves qualitative features of the systems studied. Our applications to π-π stacking of radical cations point to larger electrostatic and induction contributions than in their neutral counterparts; however, dispersion still plays a major role that becomes even more important with system size. A more extensive study of SAPT0(UHF) accuracy is currently underway in our laboratory.

See the supplementary material for Cartesian geometries, definition of employed basis sets, comparison to reference data, benchmark energies for benzene dimer radical cations, correlation of errors with charge transfer, ionization potentials of DNA base pairs, detailed energy component analysis, SAPT0(UHF) formulae, and selected total energies.

M. Hapka is acknowledged for graciously providing reference data for open-shell SAPT0. This work was supported by a postdoctoral fellowship from the Swiss NSF to J.F.G. (Grant No. P2ELP2_155351) and research grants from the United States National Science Foundation (Grant Nos. CHE-1566192 and ACI-1147843).