{'fraction_non_alphanumeric': 0.040370554542608246, 'fraction_numerical': 0.01578356078604711, 'mean_word_length': 4.858653431160984, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Benchmarking the accuracy of seniority-zero wavefunction methods for non-covalent interactions

In this paper, we scrutinize the ability of seniority-zero wavefunction-based methods to model different types of non-covalent interactions, such as hydrogen bonds, dispersion, and mixed non-covalent interactions as well as prototypical model systems with various contributions of dynamic and static electron correlation effects. Specifically, we focus on the pair Coupled-Cluster Doubles (pCCD) ansatz combined with two different flavours of dynamic energy corrections, (i) based on a perturbation theory correction and (ii) on a linearized coupled-cluster ansatz on top of pCCD. We benchmark these approaches against the A24 data set [\v{R}ez{\'{a}}{\v{c}} and Hobza, J.~Chem.~Theory~Comput., 9, 2151-2155 (2013)] extrapolated to the basis set limit and some model non-covalent complexes that feature covalent bond breaking. By dissecting different types of interactions in the A24 data set within the Symmetry-Adapted Perturbation Theory (SAPT) framework, we demonstrate that pCCD can be classified as a dispersion-free method. Furthermore, we found that both flavours of post-pCCD approaches represent encouraging and computationally more efficient alternatives to standard electronic structure methods to model weakly-bound systems, resulting in small statistical errors. Finally, a linearized coupled cluster correction on top of pCCD proved to be most reliable for the majority of investigated systems, featuring smaller non-parallelity errors compared to perturbation-theory-based approaches.


I. INTRODUCTION
A reliable quantum chemical description of potential energy landscapes of chemical reactions is of crucial importance in chemistry, physics, biology, and astrophysics [1][2][3][4][5][6]. To predict the rates of chemical reactions, one has to know the dynamics of the nuclei, which is determined by the electronic energy at a given geometry. This involves, however, the ability to accurately calculate energy differences for molecular systems with stretched bonds and in multi-configurational states. Despite the significant progress in the development of multi-reference quantum chemistry methods in the past years [5,[7][8][9][10][11][12][13][14][15][16], it is still very challenging to tackle such problems. Conventional quantum chemistry methods that employ an active orbital space to describe static and non-dynamic correlation are difficult to use because of their complexity, nonblackbox character, and factorial scaling with respect to the system size. Hence, there is a great need for further development of electronic structure methods that allow us to accurately and reliably model bond-breaking situations without defining active spaces, but at the same time recover properly subtle effects such as dispersion interactions that originate from dynamic electron correlation.
Another promising geminal-based method is the Antisymmetric Product of 1-reference orbital Geminal (AP1roG) model [42,43], also known as the pair Coupled Cluster Doubles (pCCD) ansatz [43,44], which combined with an orbital optimization protocol [45][46][47] represents a versatile tool to model strongly-correlated closed-shell systems at low computational cost. When the molecular orbitals are variationally optimized, the pCCD model becomes size-consistent and the resulting molecular basis is appropriate to model stretched bonds in contrast to standard canonical Hartree-Fock orbitals. As a consequence, potential-energy surfaces with a correct description of dissociation processes are obtained at meanfield-like cost. Yet, the exponential form of the ansatz allows for the desired (linear) scaling with the number of electrons (size-extensivity). This size-consistent and size-extensive methodology has proven to be useful in a wide scope of applications ranging from purely theoretical models [44,45,48,49] to real chemical problems including bond-breaking processes [46,47,[50][51][52][53][54][55]. The missing fraction of dynamic correlation energy on top of pCCD can be added a posteriori, similar as for APSG, using density functional theory [56][57][58], many-body perturbation theory [24,59], and coupled cluster theory [21,22] corrections. Moreover, several extensions to pCCD have been proposed to model exited states [60][61][62].
Recently, some of us [24] developed a new family of dynamic energy corrections on top of a pCCD reference function, where the influence of a diagonal and off-diagonal zero-order Hamiltonian, a single-and multideterminant wavefunction as dual, and different excitation operators used to construct the projection manifold were carefully examined [24]. Benchmarking those models against 15 reaction energies composed of weakly interacting closed-shell molecules shows that using an offdiagonal zero-order Hamiltonian and the full quantumchemical Hamiltonian as perturbation operator reduces the error with respect to CR-CC(2,3) reference data to 2 kcal/mol. The application of a linearized coupled cluster singles and doubles correction on top of pCCD brings the results even closer to the reference values. Keeping in mind that these pCCD-based methods are also capable of modeling covalently stretched bonds in closed-shell systems, these findings motivate us to test their performance for more challenging weakly-interacting systems including covalent bond-breaking processes.
Despite being much weaker in nature than forces bonding molecules internally (ionic and covalent), van der Waals or non-covalent interactions are responsible for a wide range of physical, chemical, and biological phenomena ranging from small molecules to complex systems in condensed phase [63,64]. Van der Waals forces might, for example, significantly influence molecular structures, phase transitions, and melting points. A reliable quantum mechanical description of small and medium-sized van der Waals complexes is nicely achieved within Symmetry-Adapted Perturbation Theory (SAPT) [63][64][65][66]. Unfortunately, this beautiful methodology, as other methods designed specifically to model van der Waals complexes [67], has problems when static electron correlation effects have to be captured as, for instance, present in covalent bond breaking. Standard quantum chemistry methods that are good for static electron correlation, such as, the multi-reference configuration interaction (MRCI [5]) and complete active space second-order perturbation theory (CASPT2 [68,69]), are on the other hand, not reliable for van der Waals interactions due to lack of size-consistency. Thus, we believe that our seniority-zero-based methods and their extensions represent a promising alternative that can bridge this gap.
In this article, we will focus on non-covalent interactions and scrutinize (i) to what extend the pCCD method combined with an orbital optimization protocol can describe the interaction energy in weakly-bound systems and (ii) which of the dynamic energy corrections on top of pCCD is most suitable to account for the dispersion energy in these systems. Our starting point will be the A24 data set designed by Řezáč and Hobza [70] to test different types of non-covalent interactions. This molecular test set is composed of 24 non-covalent complexes with all sorts of non-covalent interactions such as hydrogen bonds, mixed electrostatics/dispersion, and dispersion-dominated interactions including π − π stacking. Furthermore, we investigate the performance of our newly developed methods for a few non-covalently interacting monomers such as the H 2 · · · H 2 , H 2 · · · LiH, and H + 3 · · · H 2 model complexes, in which the H−H, Li−H, and H−H bonds are stretched, respectively.
This work is organized as follows. In section II, we briefly recapitulate our theoretical methods. Section III contains all computational details. Our numerical results are presented in section IV. Finally, in section V we deliver our conclusions and outlook.

II. THEORY
In this section we will briefly review the pCCD model [42,44] combined with a variational orbital optimization procedure [45,47] and post-pCCD dynamic energy corrections [22,24,59] that have been investigated in this work. For a complete picture of the methodology used, we refer to refs. 22, 24, 42, 44, 45, 47, 59. Furthermore, we will scrutinize the definition of the interaction energy, its computational protocol as well as main issues concerning the accuracy of the obtained potential energy surfaces.

A. pCCD and orbital optimization
In the pCCD ansatz, the cluster operator is restricted to electron pair excitations T p , where a † p and a p (a † p and ap) are the electron creation and annihilation operators for α (β) electrons and |0 is some independent-particle wavefunction, for instance, the Hartree-Fock (HF) determinant.
In Eq. (1), {t a i } are the electron-pair amplitudes andT p = occ i=1 virt a=1 t a i a † a a † a aīa i is the electron-pair excitation operator that excites an electron pair from an occupied (iī) to a virtual orbital (aā) with respect to |0 . Although being size-extensive by construction, pCCD is not size-consistent and thus does not yield reliable potential energy surfaces. Size-consistency can be recovered if the one-particle basis functions are optimized. This optimization can be done in a fully variational manner, [45,47] analogous to orbital-optimized coupled clus-ter [71] (OCC), or using approximate seniority-based projection techniques. [46,47] In this work, the orbital basis has been optimized using a variational orbital optimization protocol. [45,47] For pCCD, the orbitals are then chosen to minimize the pCCD energy expression subject to the constraint that the pair amplitude equations are satisfied. The energy Lagrangian reads [45,47] where {λ a i } are the Lagrange multipliers. In the above equation, κ is the generator of orbital rotations, where (κ pq ) is a skew-symmetric matrix and transforms into a new orthogonal basis with a transformation U = e −κ and e κĤ e −κ is the Hamiltonian in the rotated basis. We should note that the above sum runs over all orbital indices. Thus, in orbital-optimized pCCD, the occupiedoccupied and virtual-virtual orbital rotations are nonredundant and have to be considered in the orbital optimization scheme. In conventional CC theory, only the occupied-virtual rotations are non-redundant [71]. The requirement that the partial derivative of L with respect to the Lagrange multipliers {λ a i } evaluated for the current orbitals is stationary results in the standard set of equations for the pCCD amplitudes {t a i }, ∂L ∂λ a i κ=0 while the stationary requirement of L with respect to the geminal coefficients, ∂L ∂t a i | κ=0 = 0, leads to a set of equations for the Lagrange multipliers, analogous to the Λ-equations in CC theory, The variational orbital gradient is the partial derivative of L with respect to the orbital rotation coefficients {κ pq }, To obtain the unitary transformation matrix U , we expand the energy Lagrangian as a function of κ up to second order where A is the molecular orbital Hessian. Thus, minimizing the Lagrangian with respect to {κ pq } leads to the well-known equation for the orbital rotation coefficients After the orbital gradient and Hessian are evaluated (we employed a diagonal approximation of the exact orbital Hessian), the matrix representation of κ can be determined from the above equation. The transformation matrix is then approximated to second order as U ≈ 1 − κ + 1 2 κ † κ and orthogonalized. This transformation matrix is used to transform the current orbitals into the new basis. The above steps (evaluation of eqs. (4), (5), (6), and (8)) are repeated until convergence, for instance, in the energy and orbital gradient, is reached. No convergence acceleration techniques, such as the Direct Inversion of the Iterative Subspace, were employed.

B. Dynamic correlation energy corrections on top of pCCD
The pCCD wavefunction allows us to reliably model non-dynamic/static electron correlation effects. [] It misses, however, a large fraction of the dynamic correlation energy, which can be attributed to excitations beyond electron pairs. In the following, we will briefly discuss corrections on top of pCCD that allow us to include dynamic correlation a posteriori, focusing on different flavours of perturbation theory models and a linearized coupled cluster correction. An in depth discussion of all investigated corrections can be found in refs. 22, 24, 59.

Perturbation theory models with a pCCD reference function
The perturbation theory models that have been developed for a pCCD reference function use Rayleigh-Schrödinger perturbation theory (RSPT) of second order. In all approaches, theĤ 0 Hamiltonian is chosen to be the inactive Fock matrix (with and without off-diagonal elements), where pi||qi are the two-electron integrals in physicists' notation containing the Coulomb pi|qi and exchange pi|iq terms. Note that we have taken the zero-order Hamiltonian in its normal-product form (see also Ref. 24 for more details). Based on the choice of (Ĥ 0 ) N (with or without off-diagonal elements), the perturbation is adjusted accordingly, In the above equation, the perturbation operator has been shifted by the pCCD correlation energy (indicated by the ) so that the first order correction to the energy vanishes, where Ψ | is the dual state. The second-order energy can be calculated from the well-known expression of RSPT containing the first-order correction to the wavefunction |Ψ (1) and the shifted normal-product perturbation Hamiltonian, In all PT2 models presented in Ref. 24, the first-order correction to the wavefunction is written as a linear expansion of Slater determinants, whereT is some excitation operator that creates states that are orthogonal to the pCCD reference wavefunction, that is, we have Ψ (1) |pCCD = 0. Specifically,T is restricted to double excitations without pair excitationŝ T 2 (indicated by the ) as well as singles and doubles excitations (again without electron-pair excitations)T 1 + T 2 . The PT2 amplitudes are determined by solving a set of equations, where the sum runs over all Slater determinants included in the expansion Eq. (13) and the bar over the bra-states indicates that the final equations are spin-summed. Based on the choice of the zero-order Hamiltonian, the excitation operatorT , and the dual state Ψ |, we can derive different PT2 corrections on top of a pCCD reference function. Various, recently presented PT2 models [24,59] are summarized in Table I and scrutinized in Ref. 24. Specifically, the nature of the dual state is abbreviated by SD (single determinant) or MD (multi determinants), while the shape of the zero-order Hamiltonian is indicated by a little letter (d: diagonalF N ; o: off-diagonalF N ). The choice of the excitation operator is given in parentheses after each acronym and covers double (d) and single and double (sd) excitations. We should note that in the original PT2b model, [59] pair excitations are included in the first-order correction to the wavefunction Eq. (13) and henceT 2 includes all double excitations. To distinguish betweenT 2 andT 2 in the PT2b method, the exclusion of pair excitations is explicitly mentioned in parentheses, e.g., (d\p) indicates that T 2 has been chosen as excitation operator.

A linearized coupled cluster correction with a pCCD reference function
A different approach to include dynamic correlation effects a posteriori is to use a coupled cluster ansatz with a pCCD reference function, Specifically for pCCD, the cluster operator is chosen to contain single excitationsT 1 and (non-pair) double exci-tationsT 2 with respect to the reference determinant of pCCD. In the linearized CC (LCC) correction, the cluster amplitudes t ν are determined by solving a linear set of coupled equations where the Baker-Campbell-Hausdorff expansion e −TĤ eT =Ĥ + [Ĥ,T ] + 1 2 [[Ĥ,T ],T ] + . . . has been truncated after the second term. If we exploit the exponential form of the pCCD wavefunction, we can rewrite the above equation as Note that, by construction,T contains excitations beyond electron pairs, which are accounted for inT p . Thus, the hybrid pCCD-LCCSD approach represents a simplification of CCSD and fpCCSD [21], where non-pair amplitudes are treated linearly. In contrast to conventional LCC, where the last term of Eq. (17) is missing, thê T p operator introduces non-linear terms, which explicitly couple the pair amplitudes with all non-pair amplitudes. The energy can be calculated by projecting against the reference determinant of |pCCD , or equivalently,

C. Computation of interaction energies
The main goal of the present work is to assess whether the post-pCCD methods are suitable to model systems bound by non-covalent forces. To this end, we will use the supermolecular approach that allows us to conveniently calculate interaction energies within any given quantum chemistry method, provided it is size-consistent. Specifically, in the supermolecular approach the interaction energy is defined as where the superscripts A and B denote the interacting subsystems, AB corresponds to the dimer molecule, and X indicates the employed method (e.g., E int HF corresponds to the interaction energy calculated within the supermolecular Hartree-Fock approach).
It is well-known [72,73] that such a subtraction must be performed with caution. To this end, total energies of the monomers and of the dimer had been obtained consistently, that is, within exactly the same methods, basis sets, geometries, and numerical thresholds of accuracy.
To probe the composition of the pCCD and post-pCCD interaction energies and to better understand the role of each component in forming different types of noncovalent compounds, we will utilize the interaction energy decomposition provided within the SAPT framework. In the many-electron formulation of SAPT [63,64,66,74] the interaction energy is obtained directly, with no subtraction whatsoever. Instead, it is obtained as the sum of the lowest corrections in the perturbation series, in which the system Hamiltonian is partitioned into the intermonomer interaction V AB and the Hamiltonians corresponding to the monomersĤ A andĤ B , TheV AB operator includes all Coulomb repulsions between the electrons of monomer A and B, the repulsion between the nuclei of A and B, and all attractive interactions between the electrons of monomer A and the nuclei of monomer B, and vice versa. Such a partitioning of the Hamiltonian defines the perturbation series in terms of the λ parameter. The zeroth-order wavefunction in SAPT is the product of the wavefunctions of the noninteracting monomers. Furthermore, in the SAPT decomposition of interaction energies, it is well known that the energy expressions of arbitrary order have to be modified to ensure the proper permutation symmetry of all electrons within the dimer. This procedure (see, for example, Ref. 75) results in the appearance of the so-called exchange corrections in the perturbation series, which is required for a proper electron exchange between the monomers. Nowadays, the many-electron SAPT approach is a well-establish theory and routinely applied. For practical reasons, however, only the lowest two orders inV AB are used to calculate the interaction energies, The first two corrections on the right hand side of the above equation represent the electrostatic interaction energy of molecular charge distributions and the correction for the exchange energy, which is the dominant repulsive effect. E (2) ind is the induction energy and originates from the mutual polarization of the monomers in the field of their corresponding partners, while E terms are the dispersion and exchange-dispersion energies [76][77][78][79] and are related to the correlated motion of electrons between the monomers.
In mean-field methods such as Hartree-Fock theory or the CASSCF approach with a minimal active space, dynamic electron correlation is missing, which often results in catastrophic errors when predicting interaction energies. Most strikingly in such calculations is that the dispersion energy is missing. [80][81][82][83] Thus in dispersion dominated systems, like systems with non-polar molecules or atoms, this deficiency might even lead to the wrong sign of the interaction energy. On the other hand, mean-field methods correctly reproduce the induction and electrostatic energies as well as higher-order induction effects. [80,84] It is, however, not entirely clear how our geminal-based methods perform for interaction energies of non-covalent systems. In order to reproduce the dispersion energy in van der Waals complexes, it is essential to include intersystem excitations with at least one single excitation on each subsystem. In the (restricted) pCCD ansatz (cf. Eq. (1)), theT p operator is by construction a double excitation operator of opposite spin pairs, but bearing the same spatial parts. Hence, it is not possible for the dimer |pCCD wavefunction to reproduce the dispersion energy. Recently Garza et al. 58 made a similar observation for the helium and neon dimers (for these systems pCCD yields repulsive interaction energies), which confirms the inability of pCCD to reproduce dispersion interactions.
On the other hand, orbital-optimized pCCD is exact for two-electron systems. Hence, the intramonomer correlation energy should be reproduced to some extent. A comparison between the supermolecular pCCD and dispersion-free SAPT (labeled from now on as SAPT2df) interaction energies and the supermolecular Hartree-Fock energy will shed some light on this puzzle. Throughout this paper, we will consider dispersion-free SAPT as a sum of the following terms where SAPT2 is defined as elst,resp + E (11) exch + E (12) with elst + E exch + E (20) ind,resp + E (20) exch−ind,resp + E (20) disp + E (20) exch−disp + δ HF .
Thus, dispersion-free SAPT accounts for intramonomer correlation through electrostatics and exchange in first order as well as induction in second order. It also includes infinite-order induction effects by adding the socalled δ HF term, which is the difference between E elst + E (10) exch + E (20) ind,resp + E (20) exch−ind,resp and the supermolecular Hartree-Fock interaction energy. The superscript t in Eq. (24) indicates the "true" correlation contribution. We refer the reader to Ref. 63 for more details concerning eqs. (24) and (25).
The supermolecular pCCD-LCCSD interaction energy includes double excitations that occupy different spatial orbitals between the systems. It is, thus, quite obvious that the dispersion energy should be properly recovered within this method. However, it is well known that dispersion can be very sensitive to the quality of the monomer wavefunction [74]. In particular, the contribution of triply excited configurations is often crucial and can contribute to even 30% of the total interaction energy [85][86][87]. In this paper, we will compare this LCCSD correction to the closed-shell supermolecular CEPA(0) (or equivalently linearized CCSD) and CCSD interaction energies for the A24 set to assess the performance of pCCD-LCCSD.
Finally, we expect that the pCCD-based perturbation theories can reproduce dispersion energies similar to Møller-Plesset perturbation theory of second order (MP2). As it is well-known, the MP2 interaction energy in the supermolecular approach includes the E (20) disp interaction energy (with exchange counterpart), where the second number in the superscript refers to the correction in terms of the intramonomer correlation [78,79,88]. Such an approximation of the dispersion energy is equivalent to a Hartree-Fock description of the monomers, which can be very poor in some cases, like the interaction of π-stacked aromatic molecules [89,90].
Since the corrections to pCCD are designed to describe both static and dynamic electron correlation, we will carefully examine the interaction energies in systems where the multi-reference character of the wavefunction directly affects the interaction energy. We expect pCCD to correctly reproduce static correlation associated with the multi-configurational description of the monomers, while the LCCSD and PT models-as mentioned aboveproperly recover the dispersion energy.

A. Basis sets and extrapolation to the basis set limit
For all atoms, we considered the augmented correlation consistent series of basis sets developed by Dunning and coworkers. Specifically for the A24 test set (see Figure 1), we used the aug-cc-pVDZ, aug-cc-pVTZ, and augcc-pVQZ basis sets [91]. The corresponding auxiliary basis sets [92,93] for density fitting were employed for calculating all SAPT energies using the PSI4 software package [94,95]. In all pCCD-based methods, pivoted incomplete Cholesky decomposition (CD) was used to approximate the electron repulsion integral tensors. The threshold cutoff (δ CD ) of the CD procedure was set to values ranging from 10 −8 to 10 −10 for all systems, which is more than sufficient to retrieve µE h accuracy. Furthermore, we performed an extrapolation to the basis set limit using a two-step procedure [96,97]. First, the basis set limit of the Hartree-Fock energy was obtained by fitting an exponential function of the form [96] to the Hartree-Fock energies (E SCF ) obtained with the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets. In the above equation, X indicates the cardinal number of the basis set (2 for D, 3 for T, etc.). Second, the correlation energies were extrapolated by a two-point fit to the function [96,97] E corr (X) = E corr where E corr (X) is defined as E corr (X) = E tot (X) − E SCF (X). For the correlation energy, we included only results obtained with the aug-cc-pVTZ and aug-cc-pVQZ basis sets in the fitting procedure. In all pCCD-based methods, the Hartree-Fock energy was extrapolated via Eq. (26), then the pair correlation energy was treated using the extrapolation scheme given in Eq. (27). All a posteriori correlation corrections on top of pCCD were extrapolated in the same fashion. All other calculations were performed using the augcc-pVDZ and aug-cc-pVTZ [91] basis sets if not stated otherwise. For all investigated molecules, we applied a counterpoise correction to minimize the basis set superposition error [73] (vide infra).
In the case of H + 3 · · · H 2 two types of basis sets were employed: aug-cc-pVDZ for 2-electron full configuration interaction(FCI)-based SAPT calculations and a modified aug-cc-pVTZ, where all diffuse d-type functions were neglected (see section S2 in the Supplementary Information for more details).

B. pCCD and pCCD-based corrections
All pCCD-based calculations were performed using our own software package Piernik [98]. In the variational orbital optimization pCCD calculations, the total energy convergence threshold was set to 10 −11 E h (or tighter) with a geminal coefficients threshold of 10 −12 for the absolute tolerance in the maximum norm for the residual (or tighter), while the absolute tolerance in the maximum norm of the orbital gradient was set to 10 −6 (or tighter). For the LCCSD and all PT2 corrections with a pCCD reference function (in the optimized orbital basis), the convergence threshold for the energy was set to 10 −10 E h (or tighter), while the absolute tolerance in the maximum norm of the amplitudes residue was set to 10 −9 (or tighter).
One of the most important issues in supermolecular calculations is to ensure that both the dimer and monomers converge to consistent states during the optimization procedure. To this end, we monitored the orbitals and their occupancies in both monomer and dimer calculations, while for the dissociating monomers we employed the orbitals from an adjacent point along the dissociation pathway as initial guess orbitals (specifically, moving towards decreasing monomer separations).

C. Reference methods
All CEPA(0), CCSD, and SAPT calculation for the A24 test set were performed with the Psi4 software package [94,95]. In all calculations, we applied canonical restricted Hartree-Fock (RHF) orbitals.
The SAPT calculations included intramonomer electronic correlation as described in refs. 63, 66, 74. The density-fitting employed in all SAPT calculations used JK-fit and RI-fit auxiliary basis sets corresponding to the aug-cc-pVXZ (X = 2, 3, 4) primary basis set. The SAPT2df interaction energy was obtained using Eq. (23).
For H + 3 · · · H 2 2-electron sub-systems, we performed SAPT calculations based on a FCI wavefunction. To this end, we used the SAPT2012 package [99] (with the Molpro program) to generate the molecular integrals, which was later interfaced with the computer code utilized in Refs. 100-102. We should emphasize that this code allows us to consider only small basis sets, not larger than approximately 80 basis functions.
For LiH· · · H 2 , we also performed CASSCF and multireference configuration interaction singles and doubles (MRCI-SD) calculations. In order to correctly reproduce the dissociation limit of the LiH molecule, we included only the lowest antibonding orbital into the active space so that the active space contains the 1s orbital of Li, the σ and σ * orbitals of LiH, and the σ(1s) orbital on H 2 . We also performed FCI calculations for LiH with a frozen 1s shell for the Li atom. The energy convergence thresholds for FCI, CASSCF, and MRCI-SD were all set to 10 −7 E h . Finally, we should note that we encountered convergence difficulties in CEPA(0) calculations for molecular systems with stretched bonds or near conical intersections. Thus, the corresponding CEPA(0) results are not shown.

IV. RESULTS AND DISCUSSION
A.
pulsive than the CCSD interaction energies, they feature tighter outer fences representing 95 % (see Figure 2) of the A24 test set with one outliner (system no. 9) around −0.57 kcal/mol. We should stress that pCCD-PT2b(d) outperforms both pCCD-LCCSD and CCSD, reducing the MAE and RMSE by roughly a factor of 2. However, the (averaged) good performance of pCCD-PT2b(d) comes at the price of a wider variation compared to all tested quantum chemistry methods (see also Figure 2). Specifically, while for some systems the pCCD-PT2b(d) interaction energy is almost identical to the CCSD(T) reference energy (e.g., system 24 in Table II), a PT2b correction underestimates the interaction energy by over 50 % for other test systems, like BH 3 · · · CH 4 (see also Table II). Finally, we should stress that the PT2b variants perform best compared to other perturbative approaches on top of pCCD. Similar observations have been made for reaction energies as reported in Ref. 24.
Having tested various CC and PT2 methods that include dynamic correlation effects on top of pCCD a pos-teriori, we will now focus on the pCCD supermolecular interaction itself. As we have pointed out in Sec. II C, we expect that the pCCD method does not include any dispersion interaction between two subsystems as its wavefunction ansatz allows only for pair excitations and the dispersion interaction of lowest order arises from double excitations, where each of the monomers contributes to at least one single excitation (so-called broken-pair excitations). Although dispersion contributions in the supermolecular pCCD energy are not expected, one can argue that correlated electrostatic and induction contributions are present since pCCD includes intramonomer correlation. In particular, pCCD with optimized orbitals is exact for two-electron systems [47]. For these reasons, we decided to compare the pCCD supermolecular interaction energy with the one obtained from supermolecular Hartree-Fock calculations, which include the uncorrelated electrostatics and induction contributions, but no dispersion, as well as from dispersion-free SAPT2 calculations (see Eq. 23), which include additional contributions of correlated electrostatics and induction. This analysis allows us to scrutinize the contributions present in supermolecular pCCD. The corresponding data is presented as inset in Figure 3. Statistically, the interaction predicted by pCCD (see inset of Figure 3) tends to be more repulsive than obtained within SAPT2df and HF (median of 0.55 and 0.73 kcal/mol respectively). For systems 14-24, which are generally dominated by dispersion, pCCD agrees with both theories (with differences of approximately 0.2 kcal/mol), except for BH 3 · · · NH 3 , where both SAPT2df and HF yield smaller interaction energies by about 0.5 kcal/mol. Interestingly, pCCD predicts a slightly positive interaction energy for systems 6-13, where both SAPT2df and supermolecular HF yield either marginally positive or negative interaction energies. The absolute differences between E int pCCD , E int HF , and E int SAPT2df for the hydrogen-bonded systems [70] 1-5 are larger than for the remaining systems in the A24 test set. This might suggest that the dispersion-free interaction energy is more sensitive for systems that contain more polar molecules. What is still missing in SAPT2df   and HF is the coupling of the intramonomer correlation with the higher-order induction energy. Clearly, if intramonomer correlation strongly affects the dipole moments or polarizabilities, the remainder of the induction energy terms higher than second order is also affected by intramonomer correlation effects. This may lead to larger discrepancies in interaction energies between pCCD and SAPT2df and pCCD and HF, respectively.

B. Systems with mixed static and dynamic correlation
In this subsection, we will investigate the interaction energies in small model systems, where one of the monomers contains stretched bonds and thus the whole system features both static and dynamic electron correlation. Specifically, our model test systems include (a) T-shaped H 2 · · · H 2 with one of the monomers being stretched, (b) H 2 · · · LiH with a stretched Li-H bond, and (c) the avoided crossing in H + 3 · · · H 2 . Based on our previous observations for the A24 set, we will focus on the most promising post-pCCD corrections investigated in this work, that is, PT2b(d) and LCCSD.

T-shaped H2 · · · H2
In the H 2 · · · H 2 model system, one of the H 2 molecules is stretched beyond its equilibrium distance and hence requires a multi-determinant description of its ground state. Simultaneously, the interaction energy between the subsystems is very sensitive to the choice of the dynamic electron correlation correction. These two features make this system an ideal candidate for testing new quantum chemistry methods designed to model static and dynamic correlation. Furthermore, due to its small size, FCI reference calculations are computationally feasible and can be used to assess the accuracy and reliability of any approximate electron correlation method. This particular model system was very recently studied by Pastorczak et al. [25] and Hapka et al. [26]. Specifically, they considered the T-shaped (C 2v point group symmetry) H 2 · · · H 2 system in which the hydrogen molecule positioned on the C 2 axis was being stretched. This model system was investigated with GVB-based methods and augmented with an in-depth SAPT analysis [78,79,105]. The authors concluded that the H 2 · · · H 2 system is bound by dispersion forces and represents a stringent test case to assess the coupling between dynamic and static correlation, whose importance increases in the dissociation limit of the stretched hydrogen molecule.
Our results for post-pCCD approaches are presented in Figure 4. They are further compared to EERPA-GVB and RCCSDT reference data [25]. pCCD-LCCSD captures about 90% of the total interaction energy around the equilibrium geometry (absolute difference of 0.012 kcal/mol) and in the dissociation limit (absolute difference of 0.01 kcal/mol at R = 15.0 a.u.) with respect to FCI. In the intermediate region, however, the pCCD-LCCSD interaction energy overestimates the FCI inter- action energy by about 0.01 kcal/mol. This overcorrelation results from the fact that pCCD (and hence pCCD-LCCSD) is exact in the numerical regime of monomers built from two electron systems, while it is not exact for the dimer containing four electrons (see Table S5 of the Supporting Information) and in the case of the dimer it gives lower total energies than FCI. What is missing in the pCCD-LCCSD method are triple excitations (related to connected and some disconnected clusters contributions). As we already mentioned in previous section, the contribution of triple excitations to the dispersion energy is always negative, thus, the interaction energy calculated according to Eq. (20) is underestimated for the equilibrium distances of H 2 molecules due to the fact that the dispersion interaction is dominant in this system. We should also note that pCCD-LCCSD is not variational (similar to other CC-based methods) and might predict more than 100% of the correlation energy, hence it possible that the interaction energy is overestimated in the intermediate regime. The results obtained from the perturbative PT2b(d) correction are rather poor for the T-shaped H 2 · · · H 2 model system compared to all other investigated methods (see Figure 4). For all considered geometries, pCCD-PT2b(d) consistently underestimates the interaction en- ergy by about 0.4 kcal/mol. The largest differences are observed around the equilibrium geometry and amount to 40%. Moreover, the interaction energy limit is not correctly reproduced by the pCCD-PT2b(d) approach as the H 2 +H interaction energy appears to be repulsive, as opposed to more accurate theories. Nonetheless, pCCD-PT2b(d) can qualitatively reproduce the overall shape of the interaction well up to an H· · · H distance of approximately 5 a.u.. Finally, the observed non-parallelity error (NPE) for pCCD-PT2b(d) is 0.020 kcal/mol and slightly worse than the corresponding NPE of pCCD-LCCSD (0.018 kcal/mol). The slightly larger NPE predicted by pCCD-PT2b(d) can be associated with the artificial bump located at a distance of 5-7 a.u.

H2 · · · LiH
The H 2 · · · LiH complex with the stretched LiH molecule represents another prototype system for testing the interplay between static and dynamic correlation. The interaction energy of this model system strongly depends on the proper dissociation limit of the LiH unit, but is also altered by the induction and electrostatic energy contributions. Any deviation from the neutral separation of the LiH molecule into its individual products will generate a nonphysically high value of the dipole moment, which strongly polarizes H 2 . As a result the description of the interaction energy might be completely wrong. In contrast to H 2 · · · H 2 , the dispersion interaction is weaker than electrostatics and induction. Thus, the overall performance of the tested pCCD and post-pCCD methods should be better than in the molecular hydrogen dimer (in terms of relative errors).
In our studies we adopted a model T-shaped structure, where the Li atom is facing the center of mass of the H 2 molecule at a distance of 4.0 a.u. The LiH bond is then stretched until its dissociation limit. Simultaneously, the H 2 molecule remains fixed in its equilibrium structure for every single point considered along the potential energy surface (see the Supporting Information for more details). Within a single-determinant description (as in Hartree-Fock theory) of the stretched LiH molecule, the corresponding dissociation limit would conform to Li + · · · H − and thus result in a very large dipole moment. Furthermore, the interaction between the H 2 · · · LiH subunits is dominated by a electrostatic quadrupole-dipole interaction, which decays as R −4 , where R is the separation between the centers of mass of H 2 and LiH, respectively. When the LiH fragment is stretched, the dipole moment predicted by HF theory linearly goes to infinity. Thus, the electrostatic interaction for E int HF decays as R −3 . For an exact theory, the LiH molecule dissociates into two neutral atoms. Hence, we expect that the interaction energy very quickly flattens out as soon as the Li and H atoms stop overlapping, while the interaction energy goes to the Li+H 2 limit and becomes repulsive.
Numerical results for the H 2 · · · LiH complex are shown in Figure 5. We compare CC-based interaction energies to MRCI-SD and CASSCF calculations as well as to FCI results. As predicted above, this model system has only little contributions of dispersion interaction given the rather small differences between pCCD/CASSCF and other methods that include dynamic correlation. As expected, RHF features a wrong asymptotic behaviour of the interaction energy, which clearly does not flatten out. If the active space is properly chosen, CASSCF yields an interaction energy that goes to the correct limit, while for large interatomic Li· · · H separations (>10 a.u.) it becomes almost flat. Furthermore, the pCCD interaction energy is very close to CASSCF results: it is slightly lower in the intermediate region (5-10 a.u.), while the limit lies slightly above the limit predicted by CASSCF. As expected, the pCCD-LCCSD interaction energy curve lies slightly above RCCSD(T) results (absolute error of 0.15 kcal/mol) in the region of small Li-H separations (4 a.u.), but is closer to the FCI curve (absolute error of 0.6 kcal/mol) in the dissociation limit (13 a.u.) than the CASSCF interaction energy (absolute error of 1.1 kcal.mol). MRCI-SD calculations improve the CASSCF description of the interaction energy and reduce the absolute error to 0.27 kcal/mol with respect to FCI. Surprisingly, pCCD-PT2b(d) again outperforms pCCD-LCCSD in both regions by 0.01 and 0.15 kcal/mol for small Li-H separations and in the vicinity of dissociation, respectively. This also translates into smaller NPEs for pCCD-PT2b(d), which yields an NPE of 0.482 kcal/mol compared to 0.667 kcal/mol predicted by pCCD-LCCSD.

H + 3 · · · H2
For our final test system, we have chosen the symmetric stretching of the H 2 moiety in the H + 3 · · · H 2 supermolecule. Specifically, the molecular geometry of the H + 3 unit can be arranged in such a way that the su- That is, we can pick a molecular geometry for the H + 3 monomer near which the supermolecular system exhibits an avoided crossing [106][107][108] between the two lowestlying configurations. These configurations correspond to the (H + · · · H 2 ) and (H· · · H + 2 ) dissociation channels. Due to the avoided crossing, the H + 3 · · · H 2 dissociation process features a sudden change in the chemical character and electronic wavefunction along the chosen reaction coordinate. In this specific setup (see also Figure 6), the interaction between the individual subsystems changes from (H + · · · H 2 and H 2 ) to (H· · · H + 2 and H 2 ). Thus, two interaction regimes are present: in the former case, dispersion forces dominate, while in the latter case induction takes over. Here, we will explore the crossover between those two regimes more closely. Specifically, we will scrutinize the transition between two chemically different electronic configurations of one of the monomers.
The reaction coordinate of the H + 3 · · · H 2 model system is constructed in a specific way: the H + 3 unit is strongly stretched to form an isosceles triangle such that the distant H atom is separated from the (center of mass of the) remaining H atoms by 9 a.u. The interatomic distance in the H 2 fragment of the H + 3 unit is then being varied, while the H 2 monomer is kept fixed in its equilibrium position throughout the reaction pathway (see also Figure 6) and is positioned at 6 a.u. with respect to the center of mass of the individual H 2 fragments. For interatomic distances smaller than 2.4 a.u., the H + 3 unit features a neutral H 2 molecule that is polarized by a proton, while for stretched H-H bond lengths a charge transfer occurs and the H + 3 system is composed of a neutral H atom that interacts with a H + 2 molecular ion. In the former case, the dominant interaction energy of the H + 3 unit with the H 2 molecule is the dispersion energy between the H 2 fragments, while in the latter case the induction energy significantly increases and becomes the dominant interaction contribution as the positive charge on the H + 2 ion strongly polarizes the H 2 molecule. To scrutinize the transition between the individual contributions to the interaction energy, we have performed second order SAPT (RS and Symmetrized Rayleigh-Schrödinger (SRS)) calculations based on FCI 2-electron wavefunctions of the H + 3 and H 2 monomers. Finitefield FCI calculations of the respective monomers allow us then to obtain the second order induction. Knowing that the second-order RS energy can be decomposed into dispersion and induction, we are able to retrieve second-order dispersion as the difference between the RS and induction energy. Moreover, the SRS counterpart will provide some insights into second order exchange effects. The corresponding numerical results are presented in Figure 7. For distances smaller than 2.4 a.u., E SRS + E (2) SRS perfectly agrees with the interaction energy obtained from supermolecular FCI. Nevertheless, dispersion dominates E (2) RS in this region contributing 75-85% of the total interaction energy (see inset of Figure 7), whereas induction is responsible for the remaining 15-25%. After the crossing point (r H-H > 2.50 a.u.), the situation changes completely. Corrections calculated from up to second-order SRS do no longer match with FCI results, highlighting the importance of higher-order contributions (e.g., third-order induction). Note the small kink on the plot of the dispersion energy, which is due to the small energy denominator close to the avoided crossing of H + 3 . Having understood the nature of the interaction and the progression of the interaction energy in the H + 3 · · · H 2 model system, we can now assess the quality of pCCD-based approaches. Figure 8 shows the interaction energy obtained from pCCD, pCCD-LCCSD, and pCCD-PT2b(d) compared to the FCI reference. As one would expect, pCCD follows the shape of the FCI curve. However, due to the lack of a proper description of dynamic correlation, it severely underestimates the attractive aspect of the interaction (difference up to 0.3 kcal/mol). Adding dynamic correlation using a posteriori corrections greatly improves the performance. In the dispersion-dominated region, differences with respect to the FCI results are reduced by one order of magnitude if dynamic correlation is included on  Figure. top of the pCCD wavefunction (errors amount to 0.03 kcal/mol). For the induction-dominated part of the reaction coordinate, the LCCSD correction slightly overcorrelates (see Table S9 of the Supporting Information) the total energy of the dimer resulting in a stronger attraction in pCCD-LCCSD compared to FCI results for interatomic separations of 2.5 to 3.5 a.u. In contrast to the H 2 · · · LiH model system, PT2b(d) does not outperform LCCSD, yielding difference of about 0.06 kcal/mol in the dispersion and about 0.12 kcal/mol in induction branch, respectively.
In this article, we have investigated the performance of some recently developed post-pCCD methods in describing interaction energies in non-covalent compounds and non-covalent complexes that feature covalent bond breaking. To reliably model chemical reactions, an (approximate) quantum chemistry method has to account for both types of electron correlation effects: static and dynamic. Only in this case, it is possible to correctly reproduce equilibrium bond lengths, dissociation limits, and possible bond-breaking/forming scenarios. Given that most of the reactive potentials used in molecular dynamics simulations are obtained from expensive and complicated MRCI or perturbation theory-based approaches, we seek new types of methods that are robust, inexpensive, and at the same time rigorously size-consistent.
In this work, we first scrutinized the performance of pCCD and post-pCCD-based methods in modeling dispersion-and hydrogen-bonded systems present in the A24 set. We showed that the pCCD interaction energies do not include any dispersion energy and that pCCD gives similar results as dispersion-free SAPT (up to the second order, augmented by the δ (2) HF term). The largest differences between these two approaches amount to 0.55 kcal/mol. The best performance of post-pCCD approaches was observed for pCCD-PT2b(d) and pCCD-LCCSD. Both methods correctly model the dispersion interaction effects. The absence of triply-excited contributions in the LCCSD correction (related to connected and some disconnected clusters contributions) leads to a systematic underestimation of the interaction energy (with a ME of −0.284 kcal/mol), which is similar to CCSD (with an error of about −0.125 kcal/mol) and CEPA(0) (with a ME of −0.125 kcal/mol) interaction energies. Surprisingly, the statistically best performance for the A24 data set was obtained with the PT2b(d) approach, resulting in the smallest ME of approximately −0.118 kcal/mol. The remaining PT2b models yield similar, albeit slightly worse results. All the other perturbative corrections on top of pCCD (PT2SDd, PT2MDd, PT2SDo, PT2MDo) give errors in the order of 1 kcal/mol or more (See Table  S1 of the Supporting Information).
In the second part of our study, we focused on modeling interaction energies in complex model systems featuring bond-breaking processes: the T-shaped H 2 dimer, the LiH· · · H 2 system, and the H + 3 · · · H 2 complex. The bare pCCD interaction energy in the T-shaped H 2 dimer is purely repulsive. The more sophisticated pCCD-LCCSD model properly reproduces the shape of the potential energy. However, the pCCD-LCCSD interaction energy is underestimated by about 0.1 kcal/mol around the minimum and in the vicinity of dissociation. In the intermediate region, it slightly overestimates the reference FCI results. Nonetheless, the overall differences are very small (< 0.1 kcal/mol). The less expensive pCCD-PT2b(d) approach properly reproduces the shape of the potential of the T-shaped H 2 dimer around its equilibrium structure.
Closer to the dissociation limit, the pCCD-PT2b(d) potential, however, deviates from the reference potential more significantly. Furthermore, the complete pCCD-PT2b(d) potential energy surface is shifted by about 0.4 kcal/mol compared to the FCI reference curve. Given the performance of post-pCCD methods on the A24 set, these result are understandable. We believe that inclusion of triple excitations in these models would improve the description of dispersion interactions and the overall performance.
The T-shaped LiH· · · H 2 complex is dominated rather by induction interaction than dispersion forces, as opposed to the H 2 dimer. Hence, the overall better agreement of pCCD-LCCSD with FCI for this system was to be expected. The errors in interaction energies, however, increase for larger Li-H bond distances (the error changes from 0.02 to 0.6 kcal/mol). pCCD-PT2b(d) yields similar errors with respect to FCI as pCCD-LCCSD, albeit its performance is slightly better throughout all distances (with errors up to 0.15 kcal/mol).
Finally, we assessed the performance of various post-pCCD methods in describing the avoided crossing in the H + 3 · · · H 2 system that features a slightly distorted H + 3 ion molecular structure. The contributions from dispersion and induction rapidly change as the H 2 fragment in the H + 3 unit is stretched and correspond to either H + · · · H 2 (dispersion branch) or H· · · H + 2 (induction branch). The pCCD method systematically underestimates the interaction energy by about 0.2 kcal/mol and 0.35 kcal/mol in the dispersion and induction branch, respectively. For the dispersion branch, pCCD-LCCSD slightly underestimates the FCI interaction energy, which can be explained by neglecting triple excitation in the cluster operator. The opposite is true for the induction branch, where pCCD-LCCSD slightly overestimate the interaction energy. The PT2b(d) correction improves the pCCD interaction energies in both branches, but its performance in the induction branch is slightly worse compared to the LCCSD correction.
To conclude, our study highlights that the pCCD-LCCSD method is very competitive in terms of accuracy compared to other methods that are currently used to model bond-breaking processes, such as MRCI-SD or CASPT2. The method is size-consistent and quite robust: we do not need to define active spaces and a set of linear equations is solved for the cluster amplitudes. Furthermore, such pCCD-based methods (linearized or not) can be systematically improved by considering higherorder excitation operators in the cluster operator, like triples, etc. Since pCCD is dispersion-free, the development of a fully-fledged SAPT methodology based on a pCCD description of the monomers is promising. Some initial work in this particular direction has been recently presented by Hapka et al. [26], who show that extended RPA might be a suitable framework for the derivation of second-order contributions.

SUPPORTING INFORMATION
The Supporting Information is available free of charge on the ACS Publications website at DOI:.
Coordinates and energies for the dissociation process of model systems with mixed dynamic and static correlation energy (H 2 · · · H 2 , H 2 · · · LiH, and H + 3 · · · H 2 ). A modified aug-cc-pVTZ basis set for H. A24 statistical errors for all variants of perturbation theory corrections on top of pCCD. ACKNOWLEDGEMENT F.B., P.T., and P.Sz.Ż thank an OPUS grant of the National Science Centre, Poland, (no. 2015/19/B/ST4/02707). P.T. thanks a POLONEZ 1 grant financed by Marie-Skłodowska-Curie COFUND. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie SkÅĆodowska-Curie grant agreement No 665778. K.B. acknowledges financial support from a Marie-Skłodowska-Curie Individual Fellowship project no. 702635-PCCDX and a scholarship for outstanding young scientists from the Ministry of Science and Higher Education.
Calculations have been carried out using resources provided by Wroclaw Centre for Networking and Supercomputing (http://wcss.pl), grant No. 411.