Overcoming Artificial Multipoles in Intramolecular Symmetry-Adapted Perturbation Theory

Intramolecular symmetry-adapted perturbation theory (ISAPT) is a method to compute and decompose the noncovalent interaction energy between two molecular fragments A and B connected via a linker C. The existing ISAPT algorithm displays several issues for many fragmentation patterns, including an artificially repulsive electrostatic energy (even when the fragments are hydrogen-bonded) and very large and mutually cancelling induction and exchange-induction terms. We attribute those issues to the artificial dipole moments at the interfragment boundary, as the atoms of A and B directly connected to C are missing electrons on one of their hybrid orbitals. Therefore, we propose several new partitioning algorithms which reassign one electron, on a singly occupied link hybrid orbital, from C to each of A/B. Once the contributions from these link orbitals are added to fragment density matrices, the computation of ISAPT electrostatic, induction, and dispersion energies proceeds exactly as normal, and the exchange energy expressions need only minor modifications. Among the link partitioning algorithms introduced, the ISAPT(SIAO1) approach (in which the link orbital is obtained by a projection onto the intrinsic atomic orbitals of a given fragment followed by orthogonalization to this fragment's occupied space) leads to reasonable values of all ISAPT corrections for all fragmentation patterns and exhibits fast and systematic basis set convergence. This improvement is made possible by a significant reduction in magnitude (even though not a complete elimination) of the unphysical dipole moments at the interfragment boundaries. We demonstrate the utility of the improved ISAPT partitioning by examining intramolecular interactions in several pentanediol isomers, examples of linear and branched alkanes, and the open and closed conformations of a family of N-arylimide molecular torsion balances.


Introduction
Weak noncovalent interactions are ubiquitous in nature as they occur both between separate molecules and between nonbonded fragments of the same molecule. As a result, the knowledge and understanding of these interactions is indispensable in many areas of chemistry, physics, materials science, and biology. When the energy of a weak interaction is computed using quantum chemistry, it is highly desirable to use a computational protocol that is both quantitative and conceptual, that is, it provides an accurate overall value and its meaningful decomposition into individual terms that can be separately interpreted and understood. There are many approaches that satisfy one of these conditions but only very few methods can do both; among the latter, symmetry-adapted perturbation theory (SAPT) 1,2 is a particularly attractive choice thanks to its many robust variants and efficient computer implementations. The interaction energy in SAPT is decomposed into four major contributions: electrostatics (Coulomb interaction of unperturbed molecular charge densities, including charge penetration), exchange (the short-range repulsion stemming from the Pauli exclusion principle), induction (the polarization of one molecule by the interacting partner), and dispersion (the correlation between instantaneous charge density fluctuations).
For example, at the lowest (qualitatively accurate) SAPT0 level of theory, the interaction energy is approximated as a sum of the following corrections, grouped together to account for the four major contributions: exch + E (20) ind,resp + E (20) exch−ind,resp + δ HF + E (20) disp + E (20) exch−disp (1) In Eq. (1) and below, the consecutive superscripts denote orders of perturbation theory with respect to intermolecular interaction and intramolecular correlation (thus, there is no intramolecular correlation in SAPT0). Furthermore, the additional subscript "resp" denotes the relaxation (response) of each molecule's Hartree-Fock (HF) orbitals to the electric field of the interacting partner, and the δ HF term brings in higher-order induction and exchange-3 induction effects contained in the supermolecular HF interaction energy. SAPT has been employed to compute and interpret intermolecular interactions for numerous complexes of theoretical and experimental interest; many of these applications, together with major advances in the SAPT methodology, have been summarized in the recent reviews. [3][4][5][6][7] While there are many options to decompose the intermolecular interaction energy, and several ways to do the same thing for a covalent bond, [8][9][10] the physical decomposition of intramolecular noncovalent interactions is much less understood. While it is not at present possible to quantify the entire nonbonded interaction involving all molecular fragments, if one can identify two nonbonded fragments A and B connected to each other through a linker fragment C, there are a few options to define and partition the interaction energy between A and B in the presence of C. On one hand, one could remove C and cap the dangling bonds of A and B with hydrogen atoms. As the interaction in question is now intermolecular, one can use any of the standard SAPT variants to evaluate its contributions. 11,12 However, the fragmentation approach requires altering the molecular system and, if the linking fragment is small, artificial means to avoid an overly pronounced repulsion between the capping hydrogens. Two alternatives that do not require cutting and capping have been introduced by Gonthier and coworkers. [13][14][15] The first one involves localizing the occupied orbitals on the fragments A, B, C and constructing a non-Hermitian zeroth-order Hamiltonian 13 in which the A-B interactions are removed following the Chemical Hamiltonian approach introduced by Mayer. 16 These removed interactions are then brought back as the perturbation that gives rise to the noncovalent intramolecular interaction. 15 The need for a non-Hermitian, biorthogonal formalism makes the theory somewhat more complicated, and the resulting energy decomposition is a little different than in standard SAPT: the polarization corrections cannot be separated from their exchange counterparts (for example, the entire first-order energy is a single term), but the second-order energy includes an explicit delocalization (charge transfer) term in addition to the induction and dispersion ones.
The second approach 14 leads to the same set of corrections as intermolecular SAPT 4 (at its lowest, SAPT0 level). In this method, referred to as ISAPT, the occupied orbitals are localized according to the intrinsic bond orbital (IBO) scheme. 17 Subsequently, the HF wavefunctions for A and B are obtained by solving the HF equations for each fragment embedded in the HF potential of linker C (that is, the A-C and B-C interactions are included but the A-B one is not). The product of the A and B wavefunctions obtained in this way can then be used as the zeroth-order wavefunction for a standard SAPT expansion.
In any variant, the treatment of the interfragment boundaries is critical to the performance of the method. In the default ISAPT protocol, the doubly occupied IBOs describing the interfragment bonds are treated as part of the linker C (it is assumed that (A,C) and (B,C) are each connected by a single bond only). To maintain zero net charges on A and B (as long as the entire molecule is electrically neutral), one proton worth of the charge on the atom from A (B) that takes part in the interfragment bond is reassigned to C. Thus, no spurious charges are created at the interfragment boundaries (obviously, the presence of such charges would badly skew the SAPT electrostatic energy). However, the same cannot be said about spurious dipole moments, even when the interfragment bonds are nonpolar.
Imagine, for illustration, that the fragmentation is performed for the propane molecule to calculate the noncovalent interaction between two CH 3 -groups connected by the -CH 2linker. Then, each of the fragments A and B features a central carbon atom, with its nuclear charge reduced by one, and three hydrogen atoms located along the directions of three sp 3 hybrid orbitals. However, as the fourth hybrid orbital that would complete the tetrahedral symmetry is missing, such a fragment has a substantial dipole moment, and we might expect that the interaction between these spurious dipoles on A and B might dominate the SAPT electrostatic contribution. If the doubly occupied interfragment IBOs, together with the corresponding +1 nuclear charges, are instead reassigned to fragments A and B or equally shared between A/B and C, 14 the spurious dipole moments remain, and an additional problem appears when the linker is small and the resulting charge densities of A and B are too close together. For instance, in our propane example, such a reassignment is not acceptable because it would lead to assigning a +1 nuclear charge on the same central carbon atom to both A and B.
While the process of fragmentation and the definition of intramolecular SAPT corrections is not unique as it does not lead to any experimental observables, a useful approach should lead to ISAPT results that make physical sense. For example, an intramolecular hydrogen bond should show up as a favorable interaction both in terms of the total SAPT energy and its electrostatic component, similar to an intermolecular hydrogen bond. This is not always the case for any of the variants (link orbital assignments) of the original ISAPT method of Ref. 14. A good example is the 2,4-pentanediol molecule that involves an intramolecular hydrogen bond. In the fragmentation scheme of Fig. 1 (a), which was extensively studied in Ref. 14, the ISAPT results make sense and resemble the SAPT decomposition for the water dimer; in particular, the electrostatic energy is attractive (negative). However, in the fragmentation scheme of Fig. 1 (b), the electrostatic energy is strongly positive, which does not make physical sense as the addition of nonpolar hydrocarbon chains to A and B cannot offset the favorable dipole-dipole interaction between the -OH groups. We believe that the reason for this unphysical behavior are the artificial dipole moments resulting from the interfragment boundaries cutting through the (nonpolar) C-C bonds. tert-butyl derivative of hexaphenylethane, 18 proposed a solution to replace the embedding in the linking fragment (in this case, two central carbon atoms) by an embedding in a set of capping hydrogen atoms at very short C-H bond lengths to reduce the artificial dipoles. In most cases, this approach indeed led to a sensible, somewhat attractive electrostatic energy; however, such a correction is very system specific as it cannot be expected that the embedding in the capping hydrogens can replace an embedding in any linker fragment C. Therefore, the issue of unphysical dipoles appearing at the ISAPT interfragment boundary remains unsolved, limiting the applications of the method to highly favorable cases and/or cases where the errors cancel (which might be expected when investigating differences between similar systems).
In this work, we propose several algorithms to repartition the link bonds between fragments to reduce the unphysical fragment dipoles while preserving the embedding in the actual, unaltered chemical system. In a sense, we follow a direction stated, but not pursued,

Methodology
In the following discussion, the capital letters K, L, M, . . . will denote atomic orbital (AO) basis functions (the AO basis always covers the entire molecule). The fragments A and B interact with each other only noncovalently, but are covalently connected via the linker C.
We will assume that (A,C) and (B,C) are each connected by a single bond, with φ AC and φ BC denoting the connecting IBOs. 17 The initial set of IBOs is obtained by localizing the HF wavefunction for the entire molecule, and its orthogonal subspaces containing IBOs for the respective fragments will be denoted by B The connecting atom on A or B, besides missing one proton, is missing an electron on one of its hybrid orbitals, the one responsible for the covalent bond with C. Thus, even for a nonpolar environment such as an sp 3 carbon in a hydrocarbon molecule, an artificial dipole moment is created on the fragment A/B because one hybrid orbital, with its bonding electron, has been removed. We believe that this artificial dipole moment is the cause of often unphysical interfragment electrostatic energies in ISAPT.
The IBOs φ AC and φ BC can be viewed as linear combinations of "intrinsic hybrid orbitals" (IHOs) on the atoms connected by the linking bond -combinations of intrinsic atomic orbitals (IAOs) 17 resulting in a hybrid orbital pointing in the bond direction. We propose to partially undo the IAO→IBO transformation and identify IHOs χ x and χ y that, in some sense, constitute the A-fragment and B-fragment parts of φ AC and φ BC , respectively. More than one algorithm of determining these IHOs is possible and we will present two specific choices below. Normally, we do assume that χ x and χ y , like the whole link orbital space where they come from, are orthogonal to the occupied spaces B A and B B , respectively. Note that this orthogonalization will have to be explicitly enforced as χ x and χ y are not contained in the subspace B C ; we will also investigate what happens if one forgoes this orthogonalization.
In this process of bond orbital decomposition, we associate one electron with each of χ x and χ y and, for the purpose of updating density matrices, we will assume that half of this electron comes in with spin up and half with spin down. In our new link assignments, the single electron on χ x is reassigned from C to A together with the corresponding +1 nuclear charge on the A atom connected to C. Thus, the assignment of nuclear charges is now very simple -the full charges of all nuclei in A belong to A -but the assignment of electronic 9 charges becomes more complicated. An analogous reassignment takes place between C and B -the net result is that A and B gain one electron each relative to the default ISAPT link assignment. Our hope is that, by the addition of electrons on χ x and χ y , the fragment charges around the linking atoms become more spherically symmetrical, and any unphysical dipole moments at the fragment boundary should be substantially diminished.
At this point, the redefinition of the ISAPT first-order electrostatic energy is straightforward. The nuclear potentials of fragments A and B are modified to account for the entire charge +Ze − , not +(Z − 1)e − , on the linking atom. The HF-level fragment density matrices D A and D B are supplemented by the contributions of the (singly occupied) orbitals χ x and χ y , respectively. The updated nuclear potentials and density matrices are then used in standard SAPT expressions to determine the ISAPT E (10) elst term between fragments A and B.
The E (10) exch correction, within its customary S 2 approximation, 1 requires, in addition to the full density matrix, its partitioning into spin density matrices. It is natural to assume that the supplementary contributions of χ x and χ y are divided equally between the spin-up and spin-down density matrix (half spin-up and half spin-down). However, this assumption does not determine the exchange energy uniquely, as the interactions between the electrons on χ x and χ y depend on the mutual alignment of their spins. In the two limiting cases that both lead to equal spin-up and spin-down density matrices, the link spinorbitals ψ x and ψ y can have parallel spins, or perpendicular spins, Note that the parallel case corresponds to ψ x |ψ y = χ x |χ y and includes one electron 10 worth of exchange interaction between the link spinorbitals. The perpendicular case leads to ψ ⊥ x |ψ ⊥ y = 0; thus, no exchange interaction between the link spinorbitals exists. Both cases introduce an interdependence between the spins of the A-C and B-C link electrons which can be viewed as unphysical -the A-C and B-C partitionings should be independent of each other, leading to a random mutual orientation of spins. To model this random spin coupling, we propose to adopt as the final exchange energy the average of the parallel and perpendicular approaches: Following the standard derivation of MO-based exchange corrections within the densitymatrix formalism 21 (valid in both dimer and monomer basis sets) and recasting the resulting formulas to the AO basis, [22][23][24] we obtain the following formulas for the parallel-spin and perpendicular-spin variants: In the above equation, the upper plus/minus signs correspond to E (10) exch (S 2 ) and the lower signs to E (10)⊥ exch (S 2 ). Furthermore, the AO-basis density matrices are computed from the 12 LCAO MO coefficients C iK in the usual way: where the indices i and j run over occupied orbitals of A and B, respectively, and indices x, y correspond to the link orbitals χ x , χ y (the factors 1 2 result from the link orbital containing only half an electron of a given spin). Note that D A and D B already contain the linkorbital contributions D X and D Y , respectively, but the numerical factors multiplying the link-orbital terms need to be adjusted by the presence of the additional terms involving D X and D Y . The elements of the commonly used generalized Coulomb and exchange matrices are given by The case of the full nonapproximated E exch is more complicated and will be presented in detail in Appendix A.
A more subtle issue is defining a suitable induction energy E (20) Ind (following Ref. 25, the capitalized correction name denotes a sum of pure induction and exchange-induction terms) and, in a coherent manner, the δ HF correction that captures higher-order induction 13 and exchange-induction effects beyond that specific form of E (20) Ind . For these purposes, we note that the single electron on χ x (χ y ) does contribute (via the electron density) to the electrostatic potential of a given fragment that polarizes the other one. However, the orbitals χ x (χ y ) themselves are assumed to be frozen (unpolarizable). In other words, we consider the The coupled perturbed Hartree-Fock (CPHF) coefficients C i a /C j b for each fragment and the pure induction correction E (20) ind,resp can now be computed using standard SAPT0 formulas (note that a and b denote virtual orbitals of A and B, respectively). The exchange-induction correction E (20) exch−ind,resp , computed in ISAPT within the S 2 approximation, again depends on whether a parallel (Eq. (2)) or perpendicular (Eq. (3)) spin coupling of link electrons is assumed, and we will again opt for the average of the two as our final correction. The resulting AO-based formula for both couplings, derived via the density matrix formalism, 21 where C A is the matrix of fragment-A CPHF coefficients C i a backtransformed to the AO basis using the appropriate LCAO MO coefficients C iK and C aL : 26 An analogous formula for the other exchange induction term E Unfortunately, our proposed definition of interfragment induction energy does not lend itself to an infinite-order generalization that can be used to define δ HF . On the noninteracting fragment side, the induction energy accounts for a polarization of orbitals from B A in the electrostatic potential of B B ⊕ (1/2){χ y } and the other way around. Accordingly, one should optimize the noninteracting occupied orbitals for A in the orthogonal complement of B C , with the embedding potential including interactions with nuclei of C and the sum of electron densities coming from χ x and the rest of B C (recall that the latter has had two electrons removed and reassigned to A and B). Note that the embedding potential does not contain the interaction with the single-electron density from χ y , the IHO at the more distant fragment boundary, which has been reassigned to fragment B; thus, the entire A-B interaction has been temporarily switched off.
On the interacting side, one optimizes the occupied orbitals of AB also in the orthogonal complement of B C , with the nuclei of C and the Coulomb and exchange operators corresponding to the density of B C forming the HF embedding potential. Note that this potential is obtained from the same HF calculation for the entire molecule. Moreover, the reassignment of electrons on χ x and χ y does not matter this time because the HF system is embedded in the sum of densities coming from χ x , χ y , and the rest of B C . The orbitals χ x and χ y , while frozen in the HF optimization, do belong to the fragments A and B, so that their interaction with each other and with the nuclei of the other fragment needs to be included in the HF interaction energy (it is a constant term that does not depend on the orbitals being optimized). Unfortunately, in this algorithm, the treatment of the two electrons assigned to χ x and χ y is not consistent between the fragment calculations, where they are described by frozen orbitals, and the molecular one, where they belong to doubly occupied orbitals which are variationally optimized. As a result, we were not able to define a supermolecular δ HF correction that is consistent with the proposed reassignment of link bond electrons and the resulting E (10) elst , E (10) exch , E (20) ind,resp , and E (20) exch−ind,resp corrections. In the calculations below, this term will be taken from ISAPT calculations with the default assignment of the entire link orbitals to C.
We note at this point that the embedding defined in the fragment HF calculations described above (omitting the IHO at the distant fragment boundary), while perfectly consistent with our electrostatic and induction energy expressions, is not the same as the embedding used to compute the initial HF orbitals spanning B A and B B (which, by default, always takes into account the entire linker including both χ x and χ y ). If the linker is sufficiently large, the difference in the embedding density occurs far away from the fragment in question and its effect should be small, however, we have to consider whether an orbital reoptimization is in order before the computation of ISAPT corrections. On the other hand, χ x and χ y depend on the spaces B A and B B because they are orthogonalized to these spaces. This circular dependence suggests that iteration to self-consistency is the most rigorous way to implement the process: determine χ x and χ y which are orthogonal to the initial spaces B A and B B from the default link assignment algorithm, reoptimize B A in an embedding excluding χ y and B B in an embedding excluding χ x , update χ x /χ y by enforcing orthogonality to new B A /B B , and so on. However, a fully self-consistent algorithm is impractical -each iteration requires converged HF calculations for both fragments. On the other hand, as differences in embedding should be minor and occur relatively far away from a given fragment, this process should be very quickly convergent, and we can limit ourselves to zero, one, or two iterations. This leads to three variants of our approach: Zero iterations: the occupied spaces B A /B B are taken directly from the default ISAPT variant, that is, the HF orbitals have been obtained using an embedding in the entire B C space. These occupied spaces are used to construct χ x and χ y , and the ISAPT corrections are obtained from the resulting augmented density matrices and the polarization of the occupied spaces.
One iteration: after the IHOs χ x and χ y are constructed as above, the HF calculations for fragments are repeated, this time with the embedding potential for A excluding To complete the specification of our method, we have to define how the IHOs χ x and χ y are determined. We propose two algorithms for this, both of which proceed with a projection step (onto the fragment A or B) and an orthogonalization step (relative to all doubly occupied orbitals of this fragment). Once the specific IBO φ AC responsible for the A-C link bond is identified, its projection onto A can be performed in one of two ways.
In the first approach, inspired by the ALMO method, 27 φ AC is represented in the AO basis and projected onto the space of basis functions centered on atoms of A only (all other coefficients of this orbital are zeroed). This approach will be termed Splitting of Atomic Orbitals (SAO). In the second approach, φ AC is represented in the IAO basis, 17 projected onto the space of IAOs centered on atoms of A only by zeroing all other coefficients, and transformed back to the AO basis. This approach will be termed Splitting of Intrinsic Atomic Orbitals (SIAO). In either method, the final IHO χ x is obtained from the projected orbital by a Schmidt orthonormalization to all doubly occupied orbitals of A, that is, to the B A space. Note that the SIAO projection is not strictly an ALMO (it has nonzero coefficients of basis functions on the other fragments resulting from the tails of the IAOs), and the SAO projection ceases to be an ALMO after orthogonalization. However, in either case, χ x is still mostly localized on A. The same projection-then-orthogonalization procedure is performed at the other interfragment boundary to determine χ y .
The final specification of the new ISAPT algorithm includes both the SAO/SIAO recipe to determine link IHOs and the level of self-consistency between those IHOs and the occupied orbital space expressed by the number of iterations described earlier in this section. For example, the SAO0 method involves the SAO projection algorithm resulting in χ x /χ y orthogonalized to the original occupied spaces B A /B B (zero iterations), while the SIAO2 one uses the SIAO projection leading to χ x /χ y orthogonalized to the twice updated occupied spaces B A /B B . While the orthogonalization step is formally required to obtain a valid fragment density matrix when augmented by χ x /χ y , we will perform a limited set of numerical tests to investigate the practical consequences of not orthogonalizing (note that the lack of orthogonalization makes the algorithms with 1 and 2 iterations identical as χ X ≡ χ X ≡ χ X for X = x, y).
Finally, we note a common formal flaw of both selections of the IHOs χ x and χ y : they are not entirely contained in the occupied space B C . As a result, from the point of view of the linker C, the reassignment of one electron each to χ x and χ y involves subtracting some electron density that is not present in the system, leaving out a residual density matrix for C that is not positive definite (two negative eigenvalues). We do not see a simple way to avoid this flaw; however, it should be inconsequential in practice as only the density matrices of A and B are used to compute the ISAPT corrections and, for the purpose of recomputing the occupied orbital spaces B A and B B , the repartitioning only changes the embedding potential at the more distant interfragment boundary.

Results and discussion
Our improved ISAPT algorithms have been implemented in a development version of the Psi4 software package. 28 The modified Psi4 code is available on GitHub at https://github.com/konpat/psi4/tree/isapt. The interaction energy components within a number of representative molecules are analyzed and compared to the previously proposed ISAPT method by Parrish et al. 14 In particular, pentanediol isomers, n-heptane, 2,4dimethylpentane, and folded and unfolded bicyclic N-arylimide based molecular balances (halogens: Cl, Br, and I interact noncovalently with aromatic fragments: benzene, phenanthrene, pyrene, and ethylene) 29 are chosen for the study. All geometries are optimized at the MP2 level of theory employing the aug-cc-pVDZ basis set. To examine the accuracy of energy contributions, single-point ISAPT energy calculations are performed using three orbital bases, aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ (aDZ, aTZ, aQZ for short).
For bulky molecules, the N-arylimide balances, the systems are optimized using the B3LYP-D3(BJ) method 30 with the aDZ basis set (with cc-pVDZ-pp on the heavy iodine atom). The resulting ISAPT interaction energies are computed using the aDZ basis set, with cc-pVTZ centered on the iodine atom. Overall, each molecular system is tested with seven link as-     the fragments in P244 exhibit significantly stronger bonding than the fragments in P242 (Fig.   3), which is consistent with the addition of aliphatic chains that amplify the attractive dispersion and charge penetration effects. Both models can also be compared with an analogous hydrogen bonded configuration of the water dimer, involving the same OH· · · OH geometry with the remaining hydrogens added in the direction of the O-C bonds in 2,4-pentanediol at MP2/aDZ-optimized distances (Fig. 7). One can see that the attractive electrostatic energy obtained from the intramolecular SIAO1 method is in good agreement with the water dimer value computed with standard intermolecular SAPT0. This agreement is accidentally too good: both electrostatic terms primarily stem from the dipole-dipole interaction of the very similar polar fragments, but 2,4-pentanediol includes more charge penetration due to the hydrocarbon chains. Indeed, the presence of the additional aliphatic chains in P244 increases both the interfragment density overlap (thus leading to larger first-order exchange) and the polarizability (enhancing the induction and dispersion terms). can be obtained by examining the corresponding link orbitals (IHOs). An example comparison of the link IBOs with the IHOs resulting from both schemes, for the n-heptane molecule and its fragmentation pattern later referred to as C74, is presented in Fig. 8. For this example, we note that both schemes lead to IHOs that roughly resemble the carbon sp 3 hybrid orbitals pointing towards the other fragment, however, the SIAO1 orbitals appear to be more localized to the respective A/B fragments (note that all panels in Fig. 8 use the same isosurface value). This would suggest that the two IHOs χ x and χ y should have smaller overlap for SIAO1 than for SAO1 (as before, the primes signify that one iteration towards self-consistency of orbital spaces and link IHOs has been performed within the SAO1/SIAO1 schemes). Indeed, we observe that the χ x |χ y overlap integral tends to be much smaller in magnitude for SIAO1 than for SAO1. For example, in the 2,4-pentanediol system, the P244 fragmentation scheme, and the aDZ, aTZ, and aQZ basis sets, | χ x |χ y | amounts to 0.133, 0.0632, and 0.0724, respectively, for SAO1 and 3.02x10 −4 , 3.09x10 −4 , and 3.11x10 −4 , respectively, for SIAO1. Formally, neither SIAO1 nor SAO1 link hybrids are strictly localized on the A/B fragments: the IAOs have small tails on other atoms, and the orthogonalization does not strictly preserve the localization. The reason why the SIAO1 link hybrids are more strongly localized than the SAO1 ones is that the minimal basis used to construct IAOs 17 does not give them the flexibility to improve the basis set description on atoms other than their center. On the contrary, in a sufficiently large AO basis, the functions centered on one atom actively improve the basis set description around other atoms, although this improvement varies erratically from basis to basis. As a result, the SAO projection scheme leaves out larger tails than the SIAO one, but these tails are not stable with the basis set.
In our opinion, this behavior leads to both a larger | χ x |χ y | overlap and a larger basis set instability for the SAO1 link assignment scheme relative to the SIAO1 one.
As the OH· · · O interacting distance increases when going from 2,4-to 1,4-and 1,5pentanediols, the electrostatic contribution in a primarily dipole-dipole interaction decreases.
Among the 1,4-pentanediol models, the P142 partitioning reports the most favorable electrostatic energy (Fig. 5), which amounts to −3.04 kcal/mol for SIAO1 and −9.01 kcal/mol for SAO1 using the aDZ basis. A notable reduction in the induction term compared to the original ISAPT(C) variant occurs, leading to a net interaction energy of −4.50 kcal/mol and −0.38 kcal/mol in SIAO1 and SAO1, respectively. Interestingly, a similar trend is observed in the P156 partitioning of the 1,5-pentanediol system (Fig. 6), yet the electrostatic term can more likely be attributed to charge penetration in this case. The 1,5-pentanediol series shows that the fragments are mainly bound by the electrostatic force, along with a roughly equal contribution from the induction and dispersion terms that sum to the net attractive energy.

Alkanes
Due to the nonpolarity of the C-C bond and the negligible electronegativity difference in-  ) and (b)) and the reassigned link IHOs from ISAPT(SAO1) (panels (c) and (d)) and from ISAPT(SIAO1) (panels (e) and (f )). The system is n-heptane in the C74 fragmentation pattern, and the basis set is aDZ.
(1-propyl or 2-propyl). Thus, the A-B interaction misses the important 1,3-methyl-methyl stabilizing effects (protobranching 36 ) within each propyl group, but it misses the same number of such effects in both cases.
Both n-heptane and 2,4-dimethylpentane are nonpolar molecules, and their fragments obtained by cutting through C-C bonds should be nonpolar as well. However, as shown in  The ISAPT(SIAO1) data in Fig. 9 indicate that the overall attractive interaction energy is predominantly controlled by the dispersion contribution and secondarily by the electrostatic term, but the latter effect is more than counterbalanced by the repulsive first-order exchange. The 2,4-dimethylpentane system exhibits roughly twice the intramolecular electro-

Comparison of different SAO and SIAO variants, with and without orthogonalization
The results presented so far led us to designate the ISAPT(SIAO1) variant as the most meaningful one for practical calculations. We also made extensive comparisons to the IS- exch , as well as their averages (AVG), are presented in Fig. 14 for four systems: 2,4-pentanediol (P244), 1,4-pentanediol (P142), n-heptane (C73), and 2,4-dimethylpentane (C7B24). Due to the larger overlap of the link orbitals, the ISAPT(SAO1) first-order exchange energy is not just more erratic with respect to basis set, but it is also much more sensitive to the spin coupling than ISAPT(SIAO1). As the overlap integral between the two IHOs is zero in the perpendicular mode, the parallel coupling leads to a larger exchange energy in the SAO1 variant. In contrast, the PAR and PERP E (10) exch values from ISAPT(SIAO1) are highly consistent for all systems studied. In second order (see the Supporting Information), both E (20) exch−ind,resp and E (20) exch−disp are practically independent on spin coupling: the absolute differences between the PAR and PERP values for the systems and basis sets included in Fig. 14 do not exceed 0.042 kcal/mol (SAO1) and 2.7x10 −4 kcal/mol (SIAO1) for Figure 13: Comparison of ISAPT energy components within n-heptane (C73) and 2,4dimethylpentane (C7B24) molecules with (data marked "ORTH") and without (marked "NONE") link orbital orthogonalization, computed using the SAO1 and SIAO1 variants with the aXZ bases, X=D, T, Q. E (20) exch−ind,resp , and 0.0032 kcal/mol (SAO1) and 7.4x10 −5 kcal/mol (SIAO1) for E (20) exch−disp . For all other ISAPT results in this work, the average of the parallel and perpendicular spin-coupling values is chosen as the total exchange contribution.

N-arylimide molecular balances
Noncovalent intramolecular interactions of aromatic groups are crucial in chemical and biological processes, especially in rational drug design, new drug discovery, and the development are separately optimized; thus, their geometric difference involves some rearrangements in addition to the rotation of the substituted benzene ring. To separate the actual noncovalent intramolecular interactions from the geometry relaxation effects (the latter do not lend themselves to a SAPT-like decomposition), we also examine the "openR" structures, which are open-balance conformations obtained from the optimized closed structure directly by a 180-degree rotation of the halogen-substituted benzene ring (no geometry reoptimization is performed so the "openR" structure is not a local minimum). An analogous "closedR" structure, a closed conformation obtained by a 180-degree rotation of the open one, involves very short halogen-π distances and is strongly repulsive; therefore, we do not discuss it any further. In all structures, the noncovalent interaction energy of the halogen atom (fragment A) with the bulk of the molecule including the aromatic surface (fragment B) is investigated, with a single phenylene ring constituting the linker fragment C, as depicted in Fig. 15.

Cl-ben-closed
Br-ben-closed Br-ben-open I-phe-closed I-pyr-closed I-pyr-openR Figure 15: Representative open and closed structures of several molecular balances investigated in this work. The balances vary by the halogen atom (Cl, Br, or I) and its noncovalently interacting partner, the aromatic surface (ben = benzene, phe = phenanthrene, pyr = pyrene). A missing bond signifies a fragment boundary.
The ISAPT(C) and ISAPT(SIAO1) intramolecular interaction energy contributions for elst energy for closed conformations is usually overly attractive. While this has not happened before for the systems studied in this work, it is a manifestation of the very same issue of unphysical dipole moments at the interfragment boundary: after all, these dipole moments can be aligned unfavorably or favorably.
Second, the ISAPT(C) induction and exchange-induction terms for closed conformations are dramatically large. When the individual contributions are so overestimated, it is hard to expect their sufficiently complete cancellation: indeed, the sum E (20) ind,resp + E (20) exch−ind,resp is up to 12 times larger for ISAPT(C) than for ISAPT(SIAO1). The issues with the ISAPT(C) induction energy lead to sometimes unphysically large values of the δ HF term -see the Supporting Information. Therefore, the total ISAPT molecular balance interaction energies Interestingly, bromine-containing systems exhibit the strongest electrostatic attraction forces in the closed configuration, ranging from −0.71 to −2.99 kcal/mol compared with −1.14 to −2.31 kcal/mol in iodine systems. Furthermore, the benzene surface produces the strongest electrostatic energies for each halogen substituent, about 0.1 -0.2 kcal/mol greater than in phenanthrene and pyrene systems. The total ISAPT(SIAO1) interaction energy in closed conformations shows that binding increases with both the size of the aromatic surface and the atomic number of the halogen (Fig. 16). This relationship shows a good agreement with the experimental findings, 29 which quantified the halogen-π interactions to contribute about −1.2 kcal/mol for unimolecular systems. For the open state of the molecular balance, the "openR" structure (where the -Ph-X fragment is rotated 180 degrees without reoptimization) has a slightly larger dihedral angle (by about 6 degrees) between the aromatic surface and the imide group connecting to the linker fragment than the "open" one. With the SIAO1 method, a slightly more repulsive first-order exchange and somewhat stronger induction and dispersion energies are observed in the "openR" systems. As a result, the total interaction energy of the "open" systems is up to 0.4 kcal/mol less attractive than of the "openR" systems, as shown in Table 2. However, the overall differences in ISAPT0 interaction energy components between the "open" and "openR" structures are quite minor, showing that geometric relaxation is not a crucial factor for the performance of these molecular balances; the difference in intramolecular nonbonded interactions, which can be readily studied with ISAPT(SIAO1), is a much more important quantity. The increase in the electron density and polarizability of the halogen substituent from chlorine to iodine is correlated with their total binding energies. In the benzene based balances, the overall interaction energy increases from −0.43 kcal/mol for chlorine to −0.52 kcal/mol for bromine and to −0.98 kcal/mol for iodine in the "openR" configurations. The corresponding increase in the pyrene based balances is from −0.58 kcal/mol for Cl to −0.75 kcal/mol for Br and −1.46 kcal/mol for I. Figure 17: Decomposition of the folding energy (closed−openR) in representative molecular balance systems computed at the ISAPT0(SIAO1)/aDZ level of theory. The effects marked "Ind20" and "Disp20" represent the complete induction and dispersion components, including the respective exchange-induction and exchange-dispersion contributions.
As shown by the ISAPT0(SIAO1) decomposition, the closed conformations of N-arylimide molecular balances are mainly bound by the dispersion forces. The predicted relative folding energy can thus be computed as δE int total (closed − openR) as presented in Fig. 17. At equilibrium, due to a close halogen-π contact, the closed structure is strongly stabilized by dispersion forces, which grow with both the halogen atomic number (and thus its polarizability) and the size of the π surface. The closed conformation is secondarily stabilized, relative to the open one, by the electrostatic contribution which is, however, more than counterbalanced by the exchange repulsion term. The effect of induction energy on the relative conformational stability is negligible.

Summary
In For the corresponding exchange corrections, a choice needs to be made whether the single electrons on the A-C and B-C link orbitals undergo parallel, perpendicular, or averaged spin coupling. This is a purely technical issue that has a minor effect on the ISAPT energy contributions.
The selection of the link IHOs involves a projection of the bond orbital onto a suitably designed fragment A/B space. We propose two algorithms for this projection: in the SAO variant, the target space is spanned by the AO basis functions centered on A/B, and in the SIAO one, the target space is spanned by the intrinsic atomic orbitals 17 centered on A/B.
Each of these two choices leads to several minor variants depending on the level of selfconsistency between the noninteracting orbital space of, say, A (which should not include any interaction with the B link IHO) and the A link IHO (which should be orthogonal to the other occupied orbitals of A). We find out that a minimal level of this self-consistency is essential (the SAO0 and SIAO0 variants are not reliable), but one iteration towards selfconsistency, as in the SAO1 and SIAO1 methods, is entirely sufficient.

Associated Content
The Supporting Information is available free of charge at ...
• Derivation of the modified E (20) exch−disp expressions involving parallel and perpendicular spin couplings of link orbitals, and additional tables and figures (PDF) • Cartesian coordinates of all systems studied in this work (TXT) The modified Psi4 code including new ISAPT variants is available at https://github.com/konpat/psi4/tree/isapt.

Appendix: Nonapproximated first-order exchange energy in new ISAPT
In this Appendix, we will rederive the expression for the full, nonapproximated E Following the formulation in Ref. 25, the complete nonexpanded first-order SAPT0 energy E (10) = E (10) elst + E (10) exch is expressed as In Eq. (11), i runs over the occupied spinorbitals of A, j runs over occupied spinorbitals of B, and r, s run over occupied spinorbitals of both fragments. Furthermore, W AB is the constant intermolecular nuclear repulsion term, B ir = ψ i |V B |ψ r and A jr = ψ j |V A |ψ r are the matrix elements of the nuclear attraction potential of B/A, and ij||rs = ij|rs − ij|sr is the antisymmetrized two-electron integral in the physicists' notation. Finally, D rs are the elements of the inverse of the overlap matrix S rs of all occupied spinorbitals, D = S −1 . It is the lack of the block diagonal character of S, and thus also D, due to spin-up and spin-down coupling that necessitates a somewhat different treatment of Eq. (11) in the ISAPT/SAO and ISAPT/SIAO cases.
Once again, the explicit form of S and D depends to some extent on the spin coupling between the electrons occupying the spinorbitals ψ x and ψ y reassigned to fragments A and B in the SAO and SIAO approaches. We will consider both the parallel (Eq. (2)) and perpendicular (Eq. (3)) spin coupling of ψ x and ψ y . Let us assume that the occupied A and B spinorbitals ψ r are ordered as (ψ i (| ↑ ), ψ i (| ↓ ), ψ x , ψ j (| ↑ ), ψ j (| ↓ ), ψ y ). In this basis, the matrix S has the following explicit block form for each spin coupling: One can see that, in either case, S does not have a block-diagonal character and thus needs to be inverted as a whole (therefore, a (2N occ same-spin and opposite-spin terms, respectively. Now, going back to Eq. (11), we note that the B ir and A jr matrices are spin-diagonal so that only the same-spin block of D contributes to these terms. To despin the last term in Eq. (11), we have to break up the antisymmetrized two-electron integral and note that the ij|rs term requires the same spins within the (i, r) and (j, s) pairs while the ij|sr term imposes the same spins within the (i, s) and (j, r) pairs. As a result, the following non-zero terms contribute to the last sum in Eq. (11)  The one-and two-electron integrals involving the link orbital ψ y also differ between the parallel and perpendicular spin coupling. In the parallel case, the ir B ir D ri term has the following contributions from the link orbitals, when the index i falls on x and/or the index r falls on x or y: In the perpendicular case, the same term becomes with the B xy term canceling out in this case. The other one-electron term is completely analogous for the parallel spin coupling, and for the perpendicular one, once again, the term containing A yx cancels (note that, numerically, the parallel and perpendicular formulas are different in every term because the matrices S and D are different). We see that every link index reduces the term prefactor by a factor of √ 2, in accordance with the weight of the spinup and spin-down contributions to the overall link spinorbital. The same pattern is observed for the two-electron terms in Eq. (14), as will be illustrated by a complete decomposition of one of the terms involving both same-spin and opposite-spin contributions. Omitting the orbitals on all fragments. The coefficient C rK represents the weight of basis function K in the molecular orbital ψ r , that is, ψ r = K C rK φ AO K . Further, we define the back-transformed inverse-overlap matrix blocks D vw,spin , where v, w are orbital type indices (i, j, r, or s) and the spin case spin is either 'ss' or 'os', as in the following example: (D ri,ss ) KL = " ri C rK D ss ri C iL " = i i A second example shows the differences between the parallel spin coupling (upper signs) and perpendicular spin coupling (lower signs) in the opposite-spin case (the same-spin formulas are identical for both couplings): (D sj,os ) KL = " sj C sK D os sj C jL " = ij C iK D os ij C jL + (additionally, as stated before, for the perpendicular spin coupling D xy = 0). The intermediate expression in quotes would be appropriate in a closed-shell case where each occupied orbital of A and B holds two electrons, and it has been used before. 28,41 However, in our case we need to break down summations over individual-fragment indices into doubly occupied orbitals and the link orbital (i (i, x), j (j, y)), and the summations over both-fragment indices into all four above parts (r, s (i, x, j, y)). To accommodate the correct numerical prefactors for terms involving link orbitals in Eqs. (15)- (17), each term in Eqs. (19)- (20) is reduced by a factor of √ 2 for every link orbital x, y present. One should note that v, w in D vw,spin are not matrix indices but only denote the specific molecular orbital subblock over which the summation in Eq. (19) extends. Regardless of the v, w type, all these matrices are of the same size N AO × N AO . Finally, as implied in Eqs. makes use of the fact that the indices r and s are equivalent (have the same summation range), so after summing over the MO indices, for example, D ri,ss ≡ D si,ss . Equation (21) has the same form for the parallel and perpendicular spin coupling, but the D vw,spin matrix blocks, as explained above, are defined and computed differently.