Effect of Electron Correlation on Intermolecular Interactions: A Pair

In the DLPNO-CCSD framework, the occupied orbitals of the reference HF ..... orbital theory as the beryllium atom has a closed-shell electronic struct...
0 downloads 0 Views 1MB Size
Subscriber access provided by University of Winnipeg Library

Quantum Electronic Structure

The Effect of Electron Correlation on Intermolecular Interactions: A Pair Natural Orbitals Coupled Cluster Based Local Energy Decomposition Study Ahmet Altun, Frank Neese, and Giovanni Bistoni J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00915 • Publication Date (Web): 29 Nov 2018 Downloaded from http://pubs.acs.org on November 30, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

The Effect of Electron Correlation on Intermolecular Interactions: A Pair Natural Orbitals Coupled Cluster Based Local Energy Decomposition Study Ahmet Altun, Frank Neese,* Giovanni Bistoni* Max-Planck-Institut für Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Mülheim an der Ruhr, Germany

ABSTRACT: The development of post-Hartree-Fock (post-HF) energy decomposition schemes that are able to decompose the HF and correlation components of the interaction energy into chemically meaningful contributions is a very active field of research. One of the challenges is to provide a clear-cut quantification to the elusive London dispersion component of the intermolecular interaction. London dispersion is well-known to be a pure correlation effect and as such it is not properly described by mean field theories. In this context, we have recently developed the Local Energy Decomposition (LED) analysis, which provides a chemically meaningful decomposition of the interaction energy between two or more fragments computed at the domain-based local pair natural orbitals coupled cluster (DLPNO-CCSD(T)) level of theory. In this work, this scheme is used in conjunction with other interpretation tools to study a series of molecular adducts held together by intermolecular interactions of different nature. The HF and correlation components of the interaction energy are thus decomposed into a series of chemically meaningful contributions. Emphasis is placed on discussing the physical effects associated with the inclusion of electron correlation. It is found that four distinct physical effects can contribute to the magnitude of the correlation part of intermolecular binding energies (𝐸𝐶𝑖𝑛𝑡): (i) London dispersion, (ii) the correlation correction to the reference induction energy, (iii) correlation correction to the electron sharing process, (iv) correlation correction to the permanent electrostatics. As expected, the largest contribution to the correlation binding energy of neutral, apolar molecules is London dispersion, as in the argon dimer case. In contrast, the correction for the HF induction energy dominates 𝐸𝐶𝑖𝑛𝑡 in systems in which an ACS Paragon Plus 1 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 39

apolar molecule interacts with charged or strongly polar species, as in Ar–Li+. This effect has its origin in the systematic underestimation of polarizabilities at the HF level of theory. For similar reasons, electron sharing largely contributes to the correlation binding energy of covalently bound molecules, as in the beryllium dimer case. Finally, the correction for HF permanent electrostatics significantly contributes to 𝐸𝐶𝑖𝑛𝑡 in molecules with strong dipoles, such as water and hydrogen fluoride dimers. This effect originates from the characteristic overestimation of dipole moments at the HF level of theory, leading in some cases to positive 𝐸𝐶𝑖𝑛𝑡 values. Our results are apparently in contrast to the widely accepted view that 𝐸𝐶𝑖𝑛𝑡 is typically dominated by London dispersion, at least, in the strongly interacting region. Clearly, post-HF energy decomposition schemes are very powerful tools to analyze, categorize, and understand the various contributions to the intermolecular interaction energy. Hopefully, this will eventually lead to insights that are helpful in designing systems with tailored properties. All analysis tools presented in this work will be available free of charge in the next release of the ORCA program package.

KEYWORDS: DLPNO-CCSD(T); local energy decomposition; London dispersion; interaction energy; correlation energy. 1. INTRODUCTION Highly correlated wave function based methods are widely used for the calculation of weak intermolecular interactions within a supramolecular approach. In particular, the coupled-cluster method with single, double, and perturbative treatment of triple excitations (CCSD(T))

is

considered the “gold standard” quantum mechanical method for the calculation of interaction energies.1,2 From a physical point of view, most of the key components of the intermolecular interaction energy, such as permanent electrostatics, induction, and orbital mixing, are approximately described by the Hartree-Fock (HF) reference. To unravel these contributions, energy decomposition analysis (EDA)3–5 approaches decompose the mean field interaction energy between two molecules into

ACS Paragon Plus 2 Environment

Page 3 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

several physically meaningful terms, thus providing a quantitative framework in which to discuss experimental observables related to the nature of the intermolecular interaction. Most of these schemes are based on the pioneering work of Morokuma.6,7 By contrast, London dispersion is well-known to be a pure electron-correlation effect that is not properly described by mean field theories. This component of the interaction plays a fundamental role in controlling the stability and reactivity of a wide range of systems.8–22 Hence, the correlation contribution to the interaction energy is usually taken as an estimate for the magnitude of London dispersion.23–26 However, London dispersion is not the only contribution to the intermolecular correlation energy, and thus previous SAPT27,28 and MP2 studies29 have questioned this assumption. In particular, the long-range correlation binding energy of molecules with strong permanent dipoles was found to be dominated by the correlation correction to the permanent electrostatics obtained at the HF level.30 It was thus suggested that any post-HF EDA scheme should include terms correcting the permanent electrostatics of the reference determinant.30 In this work, we quantitatively discuss the physical contributions affecting both the reference and the correlation binding energy in a series of prototype molecular adducts held together by intermolecular interactions of different nature. The selected systems (see Figure 1) can be grouped into four distinct categories depending on the dominant interaction type: (i) London dispersiondominated interactions, e.g., the Argon (Ar) dimer;31 (ii) Induction energy-dominated interactions, e.g., the interaction of Ar with a positively charged lithium atom (Li+);32 (iii) Electron-sharing dominated interactions, e.g., the beryllium (Be) dimer;33 (iv) electrostatically dominated interactions, such as water dimer (H2O)2 and hydrogen fluoride dimer (HF)2. The recently introduced Local Energy Decomposition (LED) analysis34 in the domain-based local pair natural orbital CCSD(T) (DLPNO-CCSD(T)) framework35–44 is used to provide a clearcut definition of London dispersion at the DLPNO-CCSD(T) level as well as of all other important components of the binding energy.34 This scheme has already been successfully applied to the study of H-bonding interactions,34,45 agostic interactions,46 and the interaction of frustrated Lewis pairs.47 Its results have been also compared with those of the popular Symmetry Adapted Perturbation ACS Paragon Plus 3 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 39

Theory (SAPT).34,45 For H-bonding species, it was shown45 that LED and DFT-SAPT converge to the same values for the London dispersion and electrostatic energy components of the interaction in the asymptotic limit.

Figure 1. The investigated molecules with the key optimized bond lengths.

In order to provide a thorough analysis of the nature of such intermolecular interactions, we implemented and applied new tools into the ORCA program package48,49 within the LED scheme. They allow for the separate quantification of permanent and induced electrostatics to the interaction energy as well as for an in-depth understanding of the mechanism through which orbital relaxation lowers the total energy of the system. Moreover, to assess the contributions of individual orbitals to such relaxation effects, the so called Natural Orbitals for Chemical Valence (NOCV) theory50 combined with Extended Transition State (ETS)51 method of Ziegler, abbreviated herein as NOCVETS,52 is also implemented in the ORCA package. ACS Paragon Plus 4 Environment

Page 5 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

The paper is organized as follows. First, the theoretical aspects of the LED scheme relevant to this study and of all the new functionalities implemented in ORCA are described. Following the computational details, the results of the extended DLPNO-CCSD(T)/LED analysis are discussed separately on each interaction type. The last section is devoted to conclusions of the work.

2. THEORETICAL ASPECTS AND COMPUTATIONAL DETAILS All calculations were carried out with a development version of ORCA program package.48,49 All new functionalities described in this article will be made available to the scientific community free of charge with the next release of ORCA. 2.1. General Considerations. Before discussing the specific details of the LED analysis, we want to clarify a few aspects of interaction energy calculations in the hope to avoid confusion between the various interaction terms defined in the LED scheme below. Consider two non-interacting fragments X and Y. Let them be described by exact wave functions Ψ𝐹𝐶𝐼,𝑋 and Ψ𝐹𝐶𝐼,𝑌 (for example, full-configuration interaction (FCI) wave functions at the basis set limit). The exact wave function of the super system is then given by the antisymmetrized product Ψ𝐹𝐶𝐼,𝑋𝑌 = 𝐴(Ψ𝐹𝐶𝐼,𝑋 ∙ Ψ𝐹𝐶𝐼,𝑌), where 𝐴 is the antisymmetry operator. The associated energies would be 𝐸𝑋, 𝐸𝑌, and 𝐸𝑋𝑌 = 𝐸𝑋 + 𝐸𝑌. They can be written in terms of the HF and correlation 𝐶 energies of the fragments, e.g., 𝐸𝑋 = 𝐸𝐻𝐹 𝑋 + 𝐸𝑋. Suppose that the two fragments are now brought into

a weak interaction regime. In a Gedankenexperiment, we could perform two calculations: (i) the exact calculation on the super system, and (ii) a calculation on the super system in which the electronic structure, e.g., the densities, at the fragments are frozen at their non-interacting values. 𝐹𝐶𝐼 This provide us with energies 𝐸𝐹𝐶𝐼 𝑋𝑌 ≠ 𝐸𝑋 + 𝐸𝑌 and a “frozen” energy 𝐸𝑋𝑌 . This frozen energy will

be different from the exact energy since terms are missing that describe the changes in the electronic structure due to the presence of the interacting second fragment. In particular, the emergence or the change in the multipole moments of the fragments will be missing from the frozen energy. For this reason, we refer to

ACS Paragon Plus 5 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

𝐹𝐶𝐼 ∆𝐸𝑖𝑛𝑑𝑢𝑐𝑡𝑖𝑜𝑛,𝑋𝑌 = 𝐸𝐹𝐶𝐼 𝑋𝑌 ― 𝐸𝑋𝑌

Page 6 of 39

(1)

as the “induction” energy.53 In accordance with the discussion above, it can be further decomposed 𝐶 into a HF (∆𝐸𝐻𝐹 𝑖𝑛𝑑𝑢𝑐𝑡𝑖𝑜𝑛,𝑋𝑌) and a correlation component (∆𝐸𝑖𝑛𝑑𝑢𝑐𝑡𝑖𝑜𝑛,𝑋𝑌). Its HF part will be

synonymous with the effects of orbital relaxation defined below. The correlation contribution will contain analogous contributions as well as a number of other effects that are explained below. We point this out, since in the original LED scheme, there is a correlation contribution that describes the correction to the exaggerated HF permanent electrostatics and to the underestimated HF induced electrostatics, as it will be detailed in the following. 2.2. Local Energy Decomposition (LED) Analysis. The DLPNO-CCSD(T) and LED methods have been described extensively elsewhere.35–42 In the following, we only recall the aspects relevant to this study. Within a supermolecular approach, the energy of a molecular adduct XY relative to the total energies of its noninteracting fragments (X and Y), i.e, the binding energy of fragments (E), can be expressed as:

E = Egeo−prep + Eint

(2)

where Egeo−prep is the geometric preparation energy (also called deformation energy). It is the energy required to bring the isolated fragments from their gas-phase geometries to their in-adduct geometries. Eint is the interaction energy of the fragments and can be divided into the HF (𝐸𝐻𝐹 𝑖𝑛𝑡 ) and correlation (𝐸𝐶𝑖𝑛𝑡) contributions: 𝐶 𝐸𝑖𝑛𝑡 = 𝐸𝐻𝐹 𝑖𝑛𝑡 + 𝐸𝑖𝑛𝑡

(3)

The decomposition of HF and correlation components of the interaction energy in the LED framework is described in the following sections. 2.2.1. LED of HF Interaction Energy. In the DLPNO-CCSD framework, the occupied orbitals of the reference HF determinant are initially localized through the Foster-Boys54 or the Pipek-Mezey scheme.55 In the LED scheme, these orbitals are then assigned to the fragment in which they are dominantly localized, and thus the HF energy of the XY adduct can be decomposed ACS Paragon Plus 6 Environment

Page 7 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

into intra- and intermolecular contributions.34,47 The intramolecular HF energy corresponds to the energy of the isolated X and Y fragments when their electronic structures are distorted to be optimal for the interaction. The intermolecular HF energy includes the electrostatic (Eelstat) and exchange (Eexch) interactions between the electron densities of the distorted fragments. Based on this partitioning, the HF binding energy 𝐸𝐻𝐹 𝑖𝑛𝑡 reads: 𝐻𝐹 𝐸𝐻𝐹 𝑖𝑛𝑡 = 𝐸𝑒𝑙 ― 𝑝𝑟𝑒𝑝 + 𝐸𝑒𝑙𝑠𝑡𝑎𝑡 + 𝐸𝑒𝑥𝑐ℎ

(4)

where the static electronic preparation energy 𝐸𝐻𝐹 𝑒𝑙 ― 𝑝𝑟𝑒𝑝 is the energy required to distort the electronic structures of the isolated fragments from their ground state to the one that is optimal for the interaction and constitutes the repulsive part of 𝐸𝐻𝐹 𝑖𝑛𝑡 . It is obtained by subtracting the sum of the ground state HF energies of isolated X and Y fragments in their in-adduct geometries from the decomposed intramolecular HF energy. The LED terms just discussed are extracted from the fully relaxed localized molecular orbitals of the adduct. Hence, they already incorporate all charge transfer (CT) and polarization effects occurring upon bond formation. However, in some chemical applications, it might be useful to provide a separate quantification of these effects. With this aim, one can define a “frozen state” (Ψ0𝐻𝐹) in which CT and polarization effects are not present,56 and estimate their magnitudes as the difference between the HF energy of the adduct and the energy of Ψ0𝐻𝐹. This strategy is commonly used in Morokuma-like decomposition schemes. Although various definitions5,57,58 are possible for Ψ0𝐻𝐹, the most common approach is to define Ψ0𝐻𝐹 as the antisymmetrized product of the wave functions of the isolated fragments frozen in their in-adduct geometry as in section 2.1: Ψ0𝐻𝐹 = 𝐴(Ψ𝐻𝐹,𝑋 ∙ Ψ𝐻𝐹,𝑌)

(5)

where Ψ𝐻𝐹,𝑋 (Ψ𝐻𝐹,𝑌) is the HF determinant for the X (Y) fragment. The overall interaction energy at the HF level can be thus written as: 𝐻𝐹,0 𝐻𝐹 𝐸𝐻𝐹 𝑖𝑛𝑡 = 𝐸𝑖𝑛𝑡 + 𝐸𝑖𝑛𝑡 ― 𝑜.𝑟.

(6)

ACS Paragon Plus 7 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 39

where 𝐸𝐻𝐹,0 𝑖𝑛𝑡 is the interaction energy of the fragments in the frozen state, i.e. the difference between the energy of Ψ0𝐻𝐹 and the HF energy of the isolated fragments frozen at the same geometry. Hence, 𝐸𝐻𝐹 𝑖𝑛𝑡 ― 𝑜.𝑟. represents the energy lowering associated with the orbital relaxation from the frozen state Ψ0𝐻𝐹 to the HF determinant (Ψ𝐻𝐹, also called “SCF state” in the following), which accounts for CT and polarization contributions. For instance, it includes a part of the induction energy, which is the energy of the interaction of permanent multipoles and/or charges of one fragment with the induced multipoles on the other. Note that a separate quantification for CT and polarization contributions is not given in this work and the 𝐸𝐻𝐹 𝑖𝑛𝑡 ― 𝑜.𝑟. term includes both effects. 0 Analogous to the decomposition of 𝐸𝐻𝐹 𝑖𝑛𝑡 , upon localization of the occupied orbitals of Ψ𝐻𝐹

and their assignment to fragments, 𝐸𝐻𝐹,0 𝑖𝑛𝑡 can be decomposed into: 𝐻𝐹,0 0 0 𝐸𝐻𝐹,0 𝑖𝑛𝑡 = 𝐸𝑒𝑙 ― 𝑝𝑟𝑒𝑝 + 𝐸𝑒𝑙𝑠𝑡𝑎𝑡 + 𝐸𝑒𝑥𝑐ℎ

(7)

where 𝐸𝐻𝐹,0 𝑒𝑙 ― 𝑝𝑟𝑒𝑝 is the energy required to distort the electronic structures of the isolated fragments from their ground states to the ones they have in the frozen state due to the effect of the antisymmetry operator. Hence, it is conceptually similar to the “Pauli repulsion” term commonly found in EDA schemes. Importantly, the LED scheme allows one to decompose the repulsive ( 0 𝐸𝐻𝐹,0 𝑒𝑙 ― 𝑝𝑟𝑒𝑝) and the attractive exchange (𝐸𝑒𝑥𝑐ℎ) terms in the frozen state, which are typically reported

together in EDA schemes (with one exception23). Finally, 𝐸0𝑒𝑙𝑠𝑡𝑎𝑡 represents the electrostatic interaction between the electron densities of the fragments in the frozen state, thus providing a quantification to the permanent electrostatic interaction (i.e., the interaction between charges and/or permanent multipoles on different fragments). Hence, the LED terms of 𝐸𝐻𝐹 𝑖𝑛𝑡 defined previously can be written as a sum of the frozen state and orbital relaxation correction terms: 𝐻𝐹,0 𝐻𝐹 𝐸𝐻𝐹 𝑒𝑙 ― 𝑝𝑟𝑒𝑝 = 𝐸𝑒𝑙 ― 𝑝𝑟𝑒𝑝 + 𝐸𝑒𝑙 ― 𝑝𝑟𝑒𝑝 ― 𝑜.𝑟.

(8)

𝐸𝑒𝑙𝑠𝑡𝑎𝑡 = 𝐸0𝑒𝑙𝑠𝑡𝑎𝑡 + 𝐸𝑒𝑙𝑠𝑡𝑎𝑡 ― 𝑜.𝑟.

(9) ACS Paragon Plus 8 Environment

Page 9 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

𝐸𝑒𝑥𝑐ℎ = 𝐸0𝑒𝑥𝑐ℎ + 𝐸𝑒𝑥𝑐ℎ ― 𝑜.𝑟.

(10)

Hence, this approach allows for a quantification of the CT and polarization energy (𝐸𝐻𝐹 𝑖𝑛𝑡 ― 𝑜.𝑟.) analogous to the one provided in many variational decomposition schemes, while, at the same time, providing information on the mechanism through which these components lower the total energy of the system. 2.2.2. LED of Correlation Interaction Energy. The DLPNO-CCSD(T) correlation energy can be expressed as a sum of electron pair correlation energies ε𝑖𝑗, where the subscripts label the occupied orbitals, plus the perturbative triples correction. All pair correlation energies are initially computed at the local MP2 level. Electron pairs with a contribution to the correlation energy higher than an energy threshold (TCutPairs) are called “strong pairs (SP)” and are treated at the CCSD level in the DLPNO-CCSD framework. The remaining pairs with expected negligible contribution are kept at the local MP2 level and are called “weak pairs” (WP). In the DLPNO-CCSD equations, the virtual space is spanned by a set of Pair Natural Orbitals (PNOs) that is specific for each electron pair. Only PNOs with an occupation number larger than a threshold (TCutPNO) are kept. If “TightPNO” settings are used,43,44 DLPNO-CCSD(T) typically recovers about 99.9% of the canonical CCSD(T) correlation energy. Hence, the DLPNO-CCSD(T) correlation contribution to the interaction energy (𝐸𝐶𝑖𝑛𝑡) can be expressed as:

𝐸𝐶𝑖𝑛𝑡 = 𝐸𝐶𝑖𝑛𝑡― 𝐶𝐶𝑆𝐷 + 𝐸𝐶𝑖𝑛𝑡― (𝑇) = (𝐸𝐶𝑖𝑛𝑡― 𝑆𝑃 + 𝐸𝐶𝑖𝑛𝑡― 𝑊𝑃) + 𝐸𝐶𝑖𝑛𝑡― (𝑇)

(11)

where 𝐸𝐶𝑖𝑛𝑡― 𝐶𝐶𝑆𝐷 is the DLPNO-CCSD correlation binding energy, which includes strong 𝐸𝐶𝑖𝑛𝑡― 𝑆𝑃 and weak pairs 𝐸𝐶𝑖𝑛𝑡― 𝑊𝑃 contributions, and 𝐸𝐶𝑖𝑛𝑡― (𝑇) is the triples correction contribution to the binding energy. Based on the localization of occupied orbitals, 𝐸𝐶𝑖𝑛𝑡― 𝑊𝑃 and 𝐸𝐶𝑖𝑛𝑡― (𝑇) can be decomposed into electronic preparation and intermolecular contributions.34 However, for the dominant 𝐸𝐶𝑖𝑛𝑡― 𝑆𝑃 contribution, a more sophisticated approach is used that takes the locality of both occupied and virtual orbitals into account.34 In fact, pair correlation energies can be expressed as a

ACS Paragon Plus 9 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 39

sum of double excitation contributions from pairs of occupied orbitals to the virtual space spanned by the PNOs. The PNOs can be localized using the Foster-Boys54 or the Pipek-Mezey scheme55 and assigned to the fragment in which they are dominantly localized based on Mulliken population analysis. This simple approach allows one to define several excitation types, which are shown graphically in Figure 2.

Figure 2. Excitation types of strong electron pairs in the DLPNO-CCSD/LED framework. For the sake of simplicity, only charge transfer excitations from X to Y are shown.

Accordingly, the overall 𝐸𝐶𝑖𝑛𝑡― 𝑆𝑃 contribution can be expressed as: 𝐶 ― 𝑆𝑃 𝐶 ― 𝑆𝑃 𝐶 ― 𝑆𝑃 𝐶𝑇(𝑋→𝑌) + 𝐸𝐶𝑇(𝑋←𝑌) 𝐸𝐶𝑖𝑛𝑡― 𝑆𝑃 = (𝐸𝐶𝑒𝑙――𝑆𝑃 𝑝𝑟𝑒𝑝 + 𝐸𝐶 ― 𝑆𝑃 𝐶 ― 𝑆𝑃 ) + 𝐸𝐷𝐼𝑆𝑃 = 𝐸𝑛𝑜 ― 𝑑𝑖𝑠𝑝 + 𝐸𝐷𝐼𝑆𝑃

(12)

where 𝐸𝐶𝑒𝑙――𝑆𝑃 𝑝𝑟𝑒𝑝 is the dynamic electronic preparation energy of the adduct XY, which is obtained by subtracting from the intramolecular contribution (see Figure 2) the strong pairs energies of the isolated fragments in their in-adduct geometry. It provides a repulsive contribution to the interaction 𝐶𝑇(𝑋←𝑌) energy. The 𝐸𝐶𝑇(𝑋→𝑌) terms arise from the double excitations that do not conserve the 𝐶 ― 𝑆𝑃 and 𝐸𝐶 ― 𝑆𝑃 𝐶𝑇(𝑋→𝑌) charge within each fragment. As the 𝐸𝐶𝑒𝑙――𝑆𝑃 terms typically + 𝐸𝐶𝑇(𝑋←𝑌) 𝑝𝑟𝑒𝑝 and the 𝐸𝐶 ― 𝑆𝑃 𝐶 ― 𝑆𝑃

compensate to each other,34,47 they can be combined in most applications to represent the nondispersive strong pairs correlation contribution to the interaction energy, i.e. 𝐸𝐶𝑛𝑜――𝑆𝑃𝑑𝑖𝑠𝑝. It mostly corrects for permanent and induced electrostatics approximately described by the HF method, ACS Paragon Plus 10 Environment

Page 11 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

― 𝑆𝑃 as it will be discussed in the following. Finally, 𝐸𝐶𝐷𝐼𝑆𝑃 is the London dispersion attraction between

the fragments. Note that R. Mata and coworkers recently developed the so called “Dispersion Interaction Density” (DID) function,59 which is a simple yet powerful tool for the analysis of the London ― 𝑆𝑃 dispersion at the local MP2 level. This scheme can be extended to the study of 𝐸𝐶𝐷𝐼𝑆𝑃 . The DID

function (r) (where r = (x, y, z) is the radius vector) is defined as: Γ(𝒓) =

1 2

∑∑𝜀 𝑖𝜖𝑋 𝑗𝜖𝑌

(𝑑𝑖𝑠𝑝) 𝑖𝑗

(

)

1 1 |𝜑𝑖(𝒓)|2 + |𝜑𝑗(𝒓)|2 𝑁𝑖 𝑁𝑗

(13)

where the sum i is restricted to occupied orbitals located on fragment “X” [𝜑𝑖(𝒓)] and the sum j is restricted to occupied orbitals located on fragment “Y” [𝜑𝑗(𝒓)]. Ni and Nj are the occupation number of 𝜑𝑖(𝒓) and 𝜑𝑗(𝒓), respectively (2 in the closed-shell case), and 𝜀(𝑑𝑖𝑠𝑝) is the component of the pair 𝑖𝑗 energy corresponding to the dispersion interaction in the LED scheme. As the weak-pair contribution typically describes the correlation energy of very distant pairs of ― 𝑆𝑃 electrons, it has mainly dispersive character. Thus, it can be summed with 𝐸𝐶𝐷𝐼𝑆𝑃 for obtaining the ― 𝐶𝐶𝑆𝐷 total dispersion contribution at the DLPNO-CCSD level, labeled as 𝐸𝐶𝑑𝑖𝑠𝑝 . The remaining part of

the DLPNO-CCSD correlation interaction energy, i.e., the sum of the 𝐸𝐶𝑛𝑜――𝑆𝑃𝑑𝑖𝑠𝑝 and intramolecular contribution of weak pairs, is labeled as 𝐸𝐶𝑛𝑜――𝐶𝐶𝑆𝐷 𝑑𝑖𝑠𝑝. Thus, 𝐶 ― 𝐶𝐶𝑆𝐷

― 𝐶𝐶𝑆𝐷 + 𝐸𝑛𝑜 ― 𝑑𝑖𝑠𝑝 𝐸𝐶𝑖𝑛𝑡― 𝐶𝐶𝑆𝐷 = 𝐸𝐶𝑑𝑖𝑠𝑝

(14)

In its simplest form, the DLPNO-CCSD(T) correlation interaction energy can be thus expressed as: ― 𝐶𝐶𝑆𝐷 𝐶 ― (𝑇) + 𝐸𝐶𝑛𝑜――𝐶𝐶𝑆𝐷 𝐸𝐶𝑖𝑛𝑡 = 𝐸𝐶𝑑𝑖𝑠𝑝 𝑑𝑖𝑠𝑝 + 𝐸𝑖𝑛𝑡

(15)

Of course, depending on the system of interest, one may want to analyze the components of 𝐶 ― 𝐶𝐶𝑆𝐷 𝐸𝑛𝑜 ― 𝑑𝑖𝑠𝑝 and 𝐸𝐶𝑖𝑛𝑡― (𝑇) described above. However, for the sake of simplicity, eq. 15 will be the

base of our correlation energy analysis in this study. For the sake of clarity, all the LED terms just discussed that constitute the total binding energy are summarized in Figure 3 schematically.

ACS Paragon Plus 11 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 39

Figure 3. Schematic overview of the energy terms constituting the binding energy in the DLPNOCCS(D)/LED framework.

2.3. London Dispersion Formula. In this study, the dispersion energy obtained within the LED scheme will be compared with its estimate obtained by using the simple equation of F. London introduced for explaining the stability of weakly interacting nonpolar molecules.60,61 The simple approximate expression for the London dispersion energy between two atoms (X and Y) is as follows: 𝐸(6) 𝑑𝑖𝑠𝑝,𝐿 = ―

𝐶6

3 ∙ 𝐼𝑋 ∙ 𝐼𝑌 ∙ 𝛼𝑋 ∙ 𝛼𝑌 = ― 𝑟6 2 ∙ (𝐼𝑋 +𝐼𝑌) ∙ 𝑟6

(16)

where C6 is the coefficient of the atom pairwise induce dipole-induce dipole interaction; IX and IY are the first ionization potential of X and Y; αX and αY are the polarizabilities of X and Y; and r is the distance between X and Y.

ACS Paragon Plus 12 Environment

Page 13 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

In some cases, higher order contributions to the London equation may be significant. A simple expression for the first higher order contribution obtained by deriving a recursion formula of the C6 and C8 coefficients is as follows:62,63 𝐸(8) 𝑑𝑖𝑠𝑝,𝐿 = ―

𝐶8

=― 𝑟8

(

45 ∙ 𝐼𝑋 ∙ 𝐼𝑌 ∙ 𝛼𝑋 ∙ 𝛼𝑌 𝐼𝑋 ∙ 𝛼𝑋 2𝐼𝑋 +𝐼𝑌

8 ∙ 𝑒2 ∙ 𝑟8

+

𝐼𝑌 ∙ 𝛼𝑌

)

𝐼𝑋 +2𝐼𝑌

(17)

where C8 is the coefficient of the atom pairwise dipole-quadrupole interaction; and e is the elementary charge. Then, with these two contributions, London dispersion energy estimate becomes (8) 𝐸𝑑𝑖𝑠𝑝,𝐿 = 𝐸(6) 𝑑𝑖𝑠𝑝,𝐿 + 𝐸𝑑𝑖𝑠𝑝,𝐿

(18)

In applying these formula for the systems involving Ar and Li+ atoms, we used the experimental64–67 polarizabilities of 1.642 and 0.028 Å3, and ionization potentials of 15.759 and 76.638 eV, respectively.

2.4. Natural Orbitals for Chemical Valence and Extended Transition State (NOCV-ETS) Method. Further insights into orbital relaxation effect discussed above in Section 2.2.1 can be gained via the NOCV-ETS scheme.52 In this scheme, the difference in the electron densities of the Ψ𝐻𝐹 and Ψ0𝐻𝐹 states 𝛥𝜌(𝒓) is expressed in terms of NOCVs, i.e. the eigenfunctions 𝜑k of the electron density difference operator: 𝑁

𝑽=

∑(|𝜑

0 0 ⟩⟨𝜑𝐻𝐹 𝑖 | ― |𝜑𝑖 ⟩⟨𝜑𝑖 |)

𝐻𝐹 𝑖

(19)

𝑖=1

0 where the set of 𝜑𝐻𝐹 𝑖 constitutes the occupied HF orbitals and the set of 𝜑𝑖 constitutes the occupied

orbitals of the frozen state. The NOCVs can be grouped in pairs of complementary orbitals (𝜑-k, 𝜑k) that have eigenvalues of the same magnitude but with opposite sign. (20)

𝑽𝜑 ± 𝑘 = ± 𝑣𝑘𝜑 ± 𝑘. In terms of NOCV pairs, 𝛥𝜌(𝒓) reads: 𝑁/2

𝛥𝜌(𝒓) =

𝑁/2

∑ Δ𝜌 (𝒓) = ∑ 𝑣 [ ― 𝜑 𝑘

𝑘=1

𝑘

2 ―𝑘(𝒓)

+ 𝜑2𝑘(𝒓)]

𝑘=1

ACS Paragon Plus 13 Environment

(21)

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 39

Combining this partitioning of 𝛥𝜌(𝒓) with the ETS method, one can decompose the 𝐸𝐻𝐹 𝑖𝑛𝑡 ― 𝑜.𝑟 term into a series of contributions associated with pairs of complementary orbitals: 𝑁/2

𝐸𝐻𝐹 𝑖𝑛𝑡 ― 𝑜.𝑟. =

∑ 𝐸

𝐻𝐹 𝑖𝑛𝑡 ― 𝑜.𝑟.,𝑘

𝑘=1

𝑁/2

=

∑𝑣 [ ― 𝐹 𝑘

𝑇𝑆 ―𝑘, ― 𝑘

+ 𝐹𝑇𝑆 𝑘,𝑘]

(22)

𝑘=1

𝑇𝑆 where 𝐹𝑇𝑆 ―𝑘, ― 𝑘 and 𝐹𝑘,𝑘 are the diagonal elements of the so called “transition state Fock operator”

𝑭𝑻𝑺 in the NOCV basis (i.e., the Fock operator constructed from the average of the densities of the Ψ𝐻𝐹 and Ψ0𝐻𝐹 states).52 Note that, in the ORCA implementation of the NOCV-ETS scheme, the above analysis is possible for any arbitrary definition of the reference determinant (Ψ0𝐻𝐹 in the present study).

2.5. Computational Details. The fully relaxed geometries and potential energy scans (PESs) of the diatomic molecules (Ar2, ArLi+, and Be2) were obtained at the DLPNO-CCSD(T) level within RIJK approach68,69 using aug-cc-pV5Z basis set and the corresponding auxiliary counterparts.70–73 The optimized and experimental equilibrium bond lengths of Ar2 (3.785 Å and 3.761 Å,74 respectively) and Be2 (2.435 Å and 2.454 Å,75 respectively) are reasonably close to each other. To be consistent with our previous study,45 the geometries of (H2O)2 and (HF)2 were obtained at MP2 level within RIJK approach68,69 using aug-cc-pVTZ basis set and the corresponding auxiliary counterparts.70–73 The optimized geometry parameters of (H2O)2 and (HF)2 at this level has been already shown45 to be quite consistent with the available experimental and CCSD(T) data.76–79 In all calculations, only the core orbitals (1s for oxygen and fluorine; 1s, 2s, and 2p for argon) were frozen while all valence electrons were included in the correlation treatment. The single-point DLPNO-CCSD(T) and LED calculations of the diatomic and the larger (H2O)2 and (HF)2 molecular systems were performed using “TightPNO”43,44 settings. For the S66 benchmark set for nonconvalent interactions, DLPNO-CCSD(T)/TightPNO provides a maximum absolute total energy deviation of only 0.10 kcal/mol compared with the parent canonical CCSD(T).43,44 As one of the aims of the present study is to assess the decay of the LED components of 𝐸𝐶𝑖𝑛𝑡― 𝐶𝐶𝑆𝐷 with the intermolecular distance (including the long range in which 𝐸𝐶𝑖𝑛𝑡― 𝐶𝐶𝑆𝐷 ACS Paragon Plus 14 Environment

Page 15 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

vanishes), we decided to use even more conservative thresholds. In particular, all electron pairs were included in the CCSD treatment. Therefore, the weak-pairs term only corrects for the PNO truncation and does not affect significantly the overall correlation energy. In the DLPNO-CCSD(T) calculations, aug-cc-pV5Z and aug-cc-pVQZ basis sets were used with their matching auxiliary counterparts.70–72 As shown previously,45 the DLPNO-CCSD(T) interaction energies of (H2O)2 and (HF)2 calculated with aug-cc-pVQZ agree with the basis set superposition error corrected energies upon complete basis set extrapolation within 0.05 kcal/mol. The localization of the occupied orbitals of Ψ0𝐻𝐹 and Ψ𝐻𝐹 was obtained with the Foster-Boys scheme.54 The Pipek-Mezey55 scheme was applied for localizing PNOs in the LED scheme. The RIJK approximation was used in the HF part.

3. RESULTS AND DISCUSSION In this section, the binding energies of the systems shown in Figure 1 are decomposed using the LED scheme. The computed E values, their comparison with previously available experimental or high level computational data, and their decomposition into HF and correlation contributions are reported in Table 1. The computed E values agree very well with the available experimental data. They range from –0.27 (Ar▪▪▪Ar) to –8.40 (Ar▪▪▪Li+) kcal/mol. The HF contribution to the interaction energy is positive for Ar▪▪▪Ar and Be▪▪▪Be, and negative for Ar▪▪▪Li+, H2O▪▪▪H2O, and HF▪▪▪HF. On the other hand, the correlation binding energy is always negative and ranges from –0.51 (Ar▪▪▪Ar) to –10.60 (Be▪▪▪Be) kcal/mol. Its contribution ranges from the 17% (HF▪▪▪HF) to the 369% (Be▪▪▪Be) of the overall Eint, consistent with the fact that these adducts are held together by interactions of completely different nature. In the following, the decay of the LED components of the binding energies with the intermolecular distance r is shown and discussed for all the systems studied in this work.

ACS Paragon Plus 15 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 39

Table 1 The calculated equilibrium binding energy E, and the corresponding interaction energy

Eint of the adducts in kcal/mol. E Molecules

Ea

E

Eint Egeo-prep

Eint 𝐸𝐻𝐹 𝑖𝑛𝑡

𝐸𝐶𝑖𝑛𝑡

𝐸𝐶𝑖𝑛𝑡/Eint (%)

Ar▪▪▪Ar -0.28b -0.27 -0.27 0.00 0.24 -0.51 189 + Ar▪▪▪Li -8.40 -8.40 0.00 -5.99 -2.40 29 c Be▪▪▪Be -2.66 -2.87 -2.87 0.00 7.73 -10.60 369 HF▪▪▪HF -4.66d -4.65 -4.74 0.09 -3.92 -0.82 18 e H2O▪▪▪H2O -4.88 -5.01 -5.07 0.06 -3.69 -1.38 27 a The values based on experimental or previously available high level computational data. Note that dissociation and binding energies are equal with opposite sign. b Based on –D value (–0.284 ± 0.003 kcal/mol) derived from its experimental rovibronic spectra.74 e c Based on –D value (–2.658 ± 0.006 kcal/mol) derived from its experimental rovibronic spectra e recorded by using stimulated emission pumping method.75 d Obtained by subtracting –D value found on an experimentally refined empirical potential (–3.036 0 ± 0.003 kcal/mol)79 from the sum of the CCSD(T)/aug-cc-pVQZ level of the harmonic (1.805 ± 0.017 kcal/mol) and VPT2/aug-cc-pVQZ level of anharmonic (–0.185 ± 0.03 kcal/mol) zero-point vibrational energy contributions.80 e Obtained by subtracting the experimental –D (–3.16 ± 0.03 kcal/mol)81,82 and zero-point 0 vibrational energies (1.72 kcal/mol).78

Interaction Type 1: Dispersion-Dominated Interactions. It is useful to start our discussion with the analysis of the argon dimer interaction energy, which is expected to be dominated by London dispersion to a large extent. As already mentioned, the computed DLPNO-CCSD(T)/aug-cc-pV5Z interaction energy of Ar▪▪▪Ar at its equilibrium distance is –0.27 kcal/mol, which is reasonably close to the one derived from the experimental rovibronic spectra (–0.284 ± 0.003 kcal/mol). As expected, Ar2 is not bound at the HF level, showing a 𝐸𝐻𝐹 𝑖𝑛𝑡 value of +0.24 kcal/mol. 𝐻𝐹,0 The orbital relaxation effect that corresponds to the difference between 𝐸𝐻𝐹 is 𝑖𝑛𝑡 and 𝐸𝑖𝑛𝑡

practically negligible for any value of interatomic distance r (Figure 4), indicating that the electron density of the individual argon atoms is only marginally perturbed by the interaction at the HF level.

ACS Paragon Plus 16 Environment

Page 17 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 4. The computed binding energy profiles of Ar▪▪▪Ar as a function of the interatomic distance. The values on the black dotted vertical line belong to the equilibrium structure.

Hence, dynamic electron correlation plays a major role in the Ar▪▪▪Ar binding. It counteracts the repulsive HF energy and stabilizes the adduct. At the equilibrium distance, the total correlation effect is –0.51 kcal/mol. The DLPNO-CCSD correlation constitutes the main part of this energy (–0.45 kcal/mol). The contribution of the perturbative triples is small but noticeable (–0.06 kcal/mol). The LED analysis of the CCSD correlation gives further insights on the nature of the Ar▪▪▪Ar bonding (Figure 5, left panel). As expected, the DLPNO-CCSD contribution to the interaction energy (red triangles) is almost identical to the dispersion energy in all distance ranges. Hence, 𝐸𝐶𝑖𝑛𝑡 is dominated by the London dispersion contribution. Interestingly, the 3D contour plots of the one electron density difference (EDD) between unrelaxed83 DLPNO-CCSD and SCF (HF) densities shows that electron correlation shifts a small amount of electron density from the atoms to the middle of the Ar–Ar bond.

ACS Paragon Plus 17 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 39

Figure 5. (Left) Decomposed DLPNO-CCSD correlation interaction energy terms of Ar2 as a function of the interatomic distance. The nearly linear relation of the long-range dispersion interaction in the log-log scale is given as insert on the graph. The values on the black dotted vertical line belong to the equilibrium structure. The difference between the electron densities computed at the DLPNO-CCSD and SCF levels is shown for the equilibrium structure with the density isosurface contour value of ±2∙10-5 e/Bohr3. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively. (Right) The relation between the dispersion energy of Ar2 calculated with LED and London equation at several interatomic distances.

It is now interesting to compare our results with the one obtained using the simple London equation (eqs. 16–18). The correlation between the London and LED dispersion energy estimates of Ar2 is given in the right panel of Figure 5. The C6-only dispersion energy estimate of London equation underestimates significantly the dispersion energy with respect to the accurate LED estimate. The agreement improves significantly when the C8 contribution that accounts for instantaneous dipole-quadruple interactions is added to the London energy, as indicated by the closeness of the blue data points to the black diagonal line showing ideal relation in the right panel of Figure 5. As expected, 𝐸𝐶𝑖𝑛𝑡 is thus a good estimate for the London dispersion contribution to the intermolecular interaction of neutral species with no permanent multiples. At the equilibrium, these are typically dominated by instantaneous dipole-dipole and dipole-quadrupole interactions. ACS Paragon Plus 18 Environment

Page 19 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Interaction Type 2: Induction Energy-Dominated Interactions. The binding energy of Ar–Li+ amounts to –8.4 kcal/mol at the DLPNO-CCSD(T) level of theory (see Table 1 and Figure 6a). Of course, the argon atom has a spherical electron density distribution and no permanent multipoles. Hence, the permanent electrostatic interaction between the positively charged Li atom and the neutral Ar atom is expected to be negligible. This can be better understood using a point charge model in which the Li+ is replaced by a positive charge. In this case, the attractive interactions between the point charge and the 18 electrons on the Ar atom roughly cancel with the repulsive interaction between the point charge and the 18 protons. Thus, the negative value of the binding energy is expected to arise from the interplay of attractive induction and London dispersion components.

Figure 6. (a) The computed binding energy profiles of Ar–Li+ as a function of the interatomic distance. The values on the black dotted vertical line belong to the equilibrium structure. (b) The electron density difference of the SCF and frozen states () and its largest decomposed NOCV contributions (1, 2, and 3) constructed with the density isosurface contour value of ±0.002 e/Bohr3 at the equilibrium structure. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively.

ACS Paragon Plus 19 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 39

Ar and Li+ are not bound in the frozen state (binding energy of +4.31 kcal/mol at the equilibrium) similar to the Ar2 case (Figure 6a). This is consistent with the absence of any significant (permanent) electrostatic interaction. However, at HF level, they are bound (–5.99 kcal/mol). Hence, in contrast to Ar2, orbital relaxation is very large (–10.30 kcal/mol at the equilibrium) and constitutes the main binding contribution. The contour plot of the difference between the SCF and the frozen one electron densities (Figure 6b) at the equilibrium geometry shows a significant Ar → Li+ polarization of the electron density on the Ar atom upon bond formation due to the positive charge of the Li atom. The NOCV-ETS decomposition shows that both σ and  orbitals are significantly polarized (Figure 6b).

Figure 7. Decomposed HF interaction energy terms of the SCF (top) and frozen (bottom) states of Ar–Li+ as a function of the interatomic distance. The nearly linear relations of the long range electrostatic in the log-log scale and of the long-range exchange interaction in the semi-log scale are given as inserts on the graphs. The values on the black dotted vertical line belong to the equilibrium structure. ACS Paragon Plus 20 Environment

Page 21 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

More detailed information on the mechanism through which orbital relaxation effects lower the total energy of the system can be obtained by comparing LED analysis of HF interaction energy 𝐻𝐹,0 𝐸𝐻𝐹 𝑖𝑛𝑡 (Figure 7, top panel) with that of the frozen state 𝐸𝑖𝑛𝑡 (Figure 7, bottom panel). The LED

analysis of 𝐸𝐻𝐹 𝑖𝑛𝑡 (Figure 7, top panel) shows a significantly large electrostatic interaction Eelstat between the distorted electron densities of the fragments. In contrast, the electrostatic interaction between the undistorted electron densities in the frozen state is almost negligible (𝐸0𝑒𝑙𝑠𝑡𝑎𝑡, bottom panel of Figure 7). Hence, the overall electrostatic interaction is, as expected, dominated by the induction component. It can be understood as the electrostatic interaction between the positive charge on Li+ and the induced dipole on the argon atom. The correlation energy contribution to the interaction energy 𝐸𝐶𝑖𝑛𝑡 amounts to –2.40 kcal/mol at the equilibrium, which includes only –0.06 kcal/mol contribution from perturbative triples. The LED analysis of 𝐸𝐶𝑖𝑛𝑡― 𝐶𝐶𝑆𝐷 (Figure 8, left panel) shows that the dispersion contribution is almost negligible, in striking contrast to the Ar2 case. Hence, 𝐸𝐶𝑖𝑛𝑡 is almost entirely dominated by the nondispersive term 𝐸𝐶𝑛𝑜――𝐶𝐶𝑆𝐷 𝑑𝑖𝑠𝑝. In this system, this attractive component provides a correction to  3 𝐸𝐻𝐹 𝑖𝑛𝑡 ― 𝑜.𝑟.. In particular, HF underestimates the polarizability of the Ar atom (HF: 1.588 Å 

DLPNO-CCSD(T): 1.642 Å3; experiment:64 1.642 Å3). Consequently, it also underestimates the extent of induction and CT contributions. This is made evident by the EDD plot between the DLPNO-CCSD and the SCF (HF) reference densities (Figure 8, left panel, bottom right), which shows that electron correlation significantly enhances the Ar → Li+ polarization. Again, the London equation provides results that are in qualitative agreement with the LED ones (Figure 8, right panel). Quantitatively speaking, the C6-only estimate of the London equation agrees very well with the accurate LED estimate of the dispersion energy. The inclusion of the C8 contribution somewhat deteriorates the agreement, with the overall dispersion still being much smaller than 𝐸𝐶𝑖𝑛𝑡.

ACS Paragon Plus 21 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 39

Figure 8. (Left) Decomposed DLPNO-CCSD correlation interaction energy terms of Ar–Li+ as a function of the interatomic distance. The values on the black dotted vertical line belong to the equilibrium structure. The electron density difference of the DLPNO-CCSD and SCF states is shown for the equilibrium structure with the density isosurface contour value of ±0.0001 e/Bohr3. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively. (Right) The relation between the dispersion energy of Ar–Li+ calculated with LED and simple London equation at several interatomic distances.

Summarizing the above reported results, when neutral species with no permanent multiples interact with positively charged or highly polar molecules, 𝐸𝐶𝑖𝑛𝑡 can be dominated by the correlation correction for orbital relaxation effects, which includes induction and CT energy contributions. Hence, the magnitude of 𝐸𝐶𝑖𝑛𝑡 can not be directly used for quantifying the London dispersion contribution in these situations. However, the LED scheme can be used to disentangle dispersive and non-dispersive contributions in 𝐸𝐶𝑖𝑛𝑡, providing a useful quantification of London dispersion that is consistent with the one obtained using the simple London dispersion formula. Interaction Type 3: Electron Sharing-Dominated Interactions. The Be2 dimer is introduced as an “impossible molecule” in several textbooks based on molecular orbital theory as the beryllium atom has a closed-shell electronic structure (1s22s2).33 However, Be has long been known to form a stable dimer in gas phase that was subjected to numerous experimental and theoretical studies.75 The ACS Paragon Plus 22 Environment

Page 23 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

experimental equilibrium interaction energy of Be2 (–2.658 ± 0.006 kcal/mol) derived from the rovibronic spectra recorded by using stimulated emission pumping method75 is reasonably close to the present DLPNO-CCSD(T) prediction of –2.87 kcal/mol.

Figure 9. (a) Binding energy curve of Be2 as a function of the interatomic distance. The values on the black dotted vertical line belong to the equilibrium structure. (b) The electron density difference of the SCF and frozen states () and its decomposed NOCV contributions ( and *) constructed with the density isosurface contour value of ±0.0015 e/Bohr3 at the equilibrium structure. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively.

The interaction between the two Be atoms is repulsive in the frozen state in all distance ranges (Figure 9a). Although orbital relaxation reduces the repulsion significantly (29.31  7.73 kcal/mol at the equilibrium; i.e., by –21.58 kcal/mol), the interaction is still repulsive at the SCF level (see Figure 9a). Note that in the Be2 dimer orbital relaxation describes a completely different physical effect than in Ar–Li+. The corresponding EDD plot of the SCF and frozen states (Figure 9b) shows that it causes significant electron density accumulation at the middle of the two Be atoms. NOCVETS decomposition of this EDD (Figure 9b) shows that it arises completely from the mixing of the two doubly occupied 2s orbitals on the Be atoms to form the HOMO–1 and HOMO molecular ACS Paragon Plus 23 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 39

orbitals (MOs) at the SCF level (reported in Figure 10 at the left). The corresponding - and *NOCV components of the overall density differences contribute –14.16 and –7.42 kcal/mol to the orbital relaxation effect. Thus, they fully account for the overall orbital relaxation effect.

Figure 10 (a) Occupied valence HF molecular orbitals in Be2, i.e., HOMO (top) and HOMO–1 (bottom). (b) The corresponding Foster Boys localized orbitals. (c) The corresponding Foster Boys localized orbitals in the frozen state.

Complementary information on the source of the huge orbital relaxation energy in Be2 can be obtained through the LED analysis (Figure 11). In this scheme, the MOs are localized into atomic orbitals (Figure 10), thus providing a simple valence bond like description of the electronic structure in Be2 at the HF level. In the frozen state (lower panel), the most significant term is the repulsive electronic preparation (34.07 kcal/mol at the equilibrium). Although the attractive intermolecular exchange due to the overlap of electron densities on the Be atoms is significant (–5.98 kcal/mol at the equilibrium), it is too small compared to the repulsive electronic preparation. The repulsive electrostatic interaction is negligible in the frozen state (1.23 kcal/mol at the equilibrium).

ACS Paragon Plus 24 Environment

Page 25 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 11. Decomposed HF interaction energy terms of the SCF (top) and frozen (bottom) states of Be2 as a function of the interatomic distance. The nearly linear relations of the long range electrostatic in the log-log scale and of the long-range exchange interaction in the semi-log scale are given as inserts on the graphs. The values on the black dotted vertical line belong to the equilibrium structure.

Orbital relaxation increases all these terms significantly. The resulting accumulation of electron charge in the middle of the bond (EDD plot of SCF and frozen states in the top insert of Figure 9) increases the electrostatic interaction from 1.23 to –45.6 kcal/mol. Moreover, it also increases the spatial overlap between the orbitals located in the two fragments, as shown by the comparison of the 2s orbitals in the frozen state with the localized HF ones (Figure 10). Consequently, the inter-fragment exchange interaction increases from –5.98 to –23 kcal/mol due to orbital relaxation effects. Compared with the other presently studied molecular adducts, the interfragment exchange interaction plays a significantly larger role in the stability of Be2.

ACS Paragon Plus 25 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 39

Interestingly, an alternative decomposition of the orbital relaxation contribution into kinetic (T) and potential energy components (nuclei-electron attraction V, exchange K, and electronelectron repulsion J)34 reveals that the kinetic energy decreases due to electron sharing (T = – 78.95 kcal/mol) whilst all other components increase (V = 5.65 kcal/mol; K = 10.21 kcal/mol;

J = 41.52 kcal/mol), consistent with the Rüdenberg model of chemical bonding of a one-electron system H2+.84 Although the orbital relaxation effect is very large and stabilizes the Be2 dimer significantly, it is smaller in absolute value than the repulsive contributions in the HF binding energy. Hence, electron correlation is an essential part of the interaction. The DLPNO-CCSD correlation contributes to the interaction by –7.88 kcal/mol. Consequently, the DLPNO-CCSD interaction energy of Be2 is – 0.14 kcal/mol, which is smaller than that of Ar2. As previously pointed out,85 the inclusion of the triples contribution (–2.73 kcal/mol) significantly improves the agreement between theory (–2.87 kcal/mol) and experiment (–2.658 ± 0.006 kcal/mol).75 As seen from Figure 12, both dispersive and nondispersive correlation components of the CCSD correlation interaction are attractive and very large. The nondispersive term represents a correction to the orbital relaxation described at the HF reference. In fact, the EDD plot associated with the effect of electron correlation (see Figure 12, top) resembles the one associated to the orbital relaxation at the HF level (Figure 9). It shows that electron correlation shifts electrons from the Be atoms to the middle of the bond. However, it is worth mentioning that, in contrast to the Ar–Li+ case, London dispersion term is still the largest component of the correlation energy. In summary, the chemical bond in Be2 is dominated by electron sharing. It leads to a significant charge accumulation of electron density at the middle of the bond that increases the interfragment electrostatic and exchange attraction between the atoms while at the same time lowering the total kinetic energy. This physical effect is only partially captured at the HF level. Electron correlation significantly increases the amount of charge accumulated in the middle of the bond and further stabilizes the adduct through strong dispersion interactions.

ACS Paragon Plus 26 Environment

Page 27 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 12. Decomposed DLPNO-CCSD correlation interaction energy terms of Be2 as a function of the interatomic distance. The nearly linear relation of the long-range dispersion interaction in the log-log scale is given as insert on the graph. The values on the black dotted vertical line belong to the equilibrium structure. The electron density difference of the DLPNO-CCSD and SCF states is shown for the equilibrium structure with the density isosurface contour value of ±0.0002 e/Bohr3. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively.

Interaction Type 4: Electrostatically-Dominated Interactions. The water and hydrogen fluoride dimers [(H2O)2 and (HF)2] are prototypical examples of H-bond interactions. They are typically chosen to test newly developed methods. Recently, the stability of a series of conformers of these dimers was investigated with the LED scheme.34,45 Since interaction energies and the decomposed LED terms of the ground-state conformers of (H2O)2 and (HF)2 are quite similar to each other, herein we discuss the results on (HF)2 and include those on (H2O)2 in the Supporting Information (see Section S2). The binding energy of (HF)2 calculated presently at the DLPNO-CCSD(T)/aug-cc-pVQZ level (–4.65 kcal/mol) for its most stable equilibrium structure is almost the same as the previous estimate (–4.66 kcal/mol) obtained by subtracting –D0 value found on an experimentally refined empirical potential (–3.036 ± 0.003 kcal/mol)79 from the sum of the CCSD(T)/aug-cc-pVQZ level of the harmonic (1.805 ± 0.017 kcal/mol) and VPT2/aug-cc-pVQZ level of anharmonic (–0.185 ± 0.03 kcal/mol) zero-point vibrational energy contributions.80 ACS Paragon Plus 27 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 39

Figure 13. (a) The computed binding energy profiles of (HF)2as a function of the H-bond distance. The values on the black dotted vertical line belong to the equilibrium structure. (b) The electron density difference of the SCF and frozen states () and its largest decomposed NOCV contribution () constructed both with the density isosurface contour value of ±0.0025 e/Bohr3 at the equilibrium structure. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively.

The main contribution to the binding energy is 𝐸𝐻𝐹 𝑖𝑛𝑡 (–3.55 kcal/mol, see Figure 13a). The interaction in the frozen state 𝐸𝐻𝐹,0 𝑖𝑛𝑡 amounts to –1.29 kcal/mol. Hence, the orbital relaxation effect significantly contribute to the interaction (𝐸𝐻𝐹 𝑖𝑛𝑡 ― 𝑜.𝑟.= –2.26 kcal/mol). The EDD plot of the SCF and frozen states and its largest NOCV component (Figure 13b) demonstrate that the main contribution to the orbital relaxation in (HF)2 comes from orbitals of local -symmetry. Significant charge depletion regions of local -symmetry are located at the H-bond acceptor (F atom). This is accompanied by areas of local -symmetry charge accumulation in the Hbond region and by the subsequent polarization of the H–F bond in the H-bond donor site. Hence, orbital relaxation shifts electron density from the H-bond acceptor to the H-bond region as a result of the increased delocalization of the lone pair located on the H-bond acceptor towards the proton. ACS Paragon Plus 28 Environment

Page 29 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 14. Decomposed HF interaction energy terms of the SCF (top) and frozen (bottom) states of (HF)2 as a function of the H-bond distance. The nearly linear relations of the long-range electrostatic in the log-log scale and of the long-range exchange interaction in the semi-log scale are given as inserts on the graphs. The values on the black dotted vertical line belong to the equilibrium structure.

Further in-depth analysis of the interaction is possible via the LED analysis of SCF and frozen states (Figure 14). In general, the energy profiles of all contributions are qualitatively similar in both states, with a large and positive electronic preparation that is largely counteracted by an electrostatic interaction, and with a small exchange contribution. Quantitatively speaking, the overall 𝐸𝐻𝐹,0 𝑖𝑛𝑡 at 0 0 the equilibrium is –1.66 kcal/mol, with its constituting 𝐸𝐻𝐹,0 𝑒𝑙 ― 𝑝𝑟𝑒𝑝, 𝐸𝑒𝑙𝑠𝑡𝑎𝑡, and 𝐸𝑒𝑥𝑐ℎ contributions

that are 17.37, –16.36, and –2.67 kcal/mol, respectively. 𝐸𝑒𝑙𝑠𝑡𝑎𝑡 is 4.83 kcal/mol larger in absolute value than 𝐸0𝑒𝑙𝑠𝑡𝑎𝑡. This can be taken as a measure of induced electrostatics. Analogously, 𝐸𝐻𝐹 𝑒𝑙 ― 𝑝𝑟𝑒𝑝

ACS Paragon Plus 29 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 39

is 3.09 kcal/mol larger than 𝐸𝐻𝐹,0 𝑒𝑙 ― 𝑝𝑟𝑒𝑝, whilst exchange is only slightly (–0.52 kcal/mol) affected by orbital relaxation. The correlation energy provides also a notable stabilizing contribution to the depth of the potential energy well (see Figure 13). At equilibrium, it amounts to –1.10 kcal/mol (the corresponding interaction energy that excludes geometric preparation: –0.82 kcal/mol). The triples contributes to the interaction energy only by –0.13 kcal/mol while the DLPNO-CCSD contribution is –0.69 kcal/mol. The decomposition of the major part (DLPNO-CCSD) of correlation interaction ― 𝑆𝑃 energy into its dispersive (𝐸𝐶𝐷𝐼𝑆𝑃 ) and nondispersive (𝐸𝐶𝑛𝑜――𝑆𝑃𝑑𝑖𝑠𝑝) parts is shown in Figure 15.

Figure 15. Decomposed DLPNO-CCSD correlation interaction energy terms of (HF)2 as a function of the H-bond distance. The nearly linear relation of the long-range dispersion interaction in the loglog scale is given as insert on the graph. The values on the black dotted vertical line belong to the equilibrium structure. The electron density difference of the DLPNO-CCSD and SCF states is shown at the top right corner for the equilibrium structure with the density isosurface contour value of ±0.002 e/Bohr3. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively.

The major part of the correlation contribution is of dispersive character (–0.96 kcal/mol at the equilibrium). As expected, it decays as r-6 in the long range (see the insert in Figure 15). Note that the nondispersive part of the CCSD interaction energy (𝐸𝐶𝑛𝑜――𝐶𝐶𝑆𝐷 𝑑𝑖𝑠𝑝) is very small and positive. It is ACS Paragon Plus 30 Environment

Page 31 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

mostly a correction to the permanent electrostatic interaction calculated at the HF level. This effect originates from the characteristic overestimation of dipole moments at the HF level. In fact, the correlation component of the electron density (see Figure 15, top right) shows a significant electron density shifts in the F → H direction, which reduces the partial charges on the corresponding atoms. 𝐶 30 Note that the 𝐸𝐶𝑛𝑜――𝐶𝐶𝑆𝐷 𝑑𝑖𝑠𝑝 term dominates in the long range, leading to positive ∆𝐸𝑖𝑛𝑡 values.

4. CONCLUSIONS The LED scheme in the DLPNO-CCSD(T) framework provides a chemically meaningful decomposition of the interaction energy between two or more fragments. In this study, both the reference HF and the DLPNO-CCSD(T) correlation contributions to the binding energy of a series of molecular adducts held together by different interaction types were decomposed using an extension of the LED scheme. Particular emphasis was placed on discussing the physical effects incorporated into the electron correlation contribution to the interaction energy (𝐸𝐶𝑖𝑛𝑡). In contrast to the widely accepted view that considers 𝐸𝐶𝑖𝑛𝑡 as purely dispersive, it has been found that 𝐸𝐶𝑖𝑛𝑡 can be dominated by a number of different physical contributions, depending on the properties of the interacting fragments. The correlation binding energy of neutral, apolar molecules or atoms, such as Ar▪▪▪Ar, is indeed dominated by dispersion. By contrast, when a neutral apolar molecule interacts with charged or strongly polar species, as in Ar–Li+, the correlation binding energy is dominated by the induced electrostatics (induction energy) correction to the HF electrostatics. This effect originates from the characteristic underestimation of polarizabilities at the HF level. Analogously, orbital relaxation effects contribute significantly to the correlation binding energy of the Be▪▪▪Be interaction, which is of a substantial covalent nature. Finally, the correlation binding energy of molecules with strong dipoles, such as (H2O)2 and (HF)2, incorporates a correction to the HF permanent electrostatics. This effect arises from the characteristic overestimation of dipole moments at the HF level. These results clearly suggest that a meaningful decomposition of the correlation binding energy is necessary for quantifying the dispersion energy properly. The LED scheme provides a ACS Paragon Plus 31 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 39

clear-cut definition of dispersive and non-dispersive components of the correlation energy at the DLPNO-CCSD(T) level of theory.

ASSOCIATED CONTENT Supporting Information LED terms of all studied adducts at their equilibrium geometries. The distance dependence of LED terms of water dimer. This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected]; [email protected]. Notes The authors declare no competing financial interest.

ACKNOWLEDGMENTS We gratefully acknowledge the Priority Program “Control of Dispersion Interactions in Molecular Chemistry” (SPP 1807) of the Deutsche Forschungsgemeinschaft for financial support.

REFERENCES (1)

Ramabhadran, R. O.; Raghavachari, K. Extrapolation to the Gold-Standard in Quantum Chemistry: Computationally Efficient and Accurate CCSD(T) Energies for Large Molecules Using an Automated Thermochemical Hierarchy. J. Chem. Theory Comput. 2013, 9, 3986– 3994.

(2)

Řezáč, J.; Hobza, P. Extrapolation and Scaling of the DFT-SAPT Interaction Energies toward the Basis Set Limit. J. Chem. Theory Comput. 2011, 7, 685–689.

(3)

Pastorczak, E.; Corminboeuf, C. Perspective: Found in Translation: Quantum Chemical Tools for Grasping Non-Covalent Interactions. J. Chem. Phys. 2017, 146, 120901.

(4)

von Hopffgarten, M.; Frenking, G. Energy Decomposition Analysis. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 43–62.

ACS Paragon Plus 32 Environment

Page 33 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(5)

Journal of Chemical Theory and Computation

Phipps, M. J. S.; Fox, T.; Tautermann, C. S.; Skylaris, C.-K. Energy Decomposition Analysis Approaches and Their Evaluation on Prototypical Protein–drug Interaction Patterns. Chem. Soc. Rev. 2015, 44, 3177–3211.

(6)

Morokuma, K. Molecular Orbital Studies of Hydrogen Bonds. III. C=O···H–O Hydrogen Bond in H2CO···H2O and H2CO···2H2O. J. Chem. Phys. 1971, 55, 1236–1244.

(7)

Kitaura, K.; Morokuma, K. A New Energy Decomposition Scheme for Molecular Interactions within the Hartree-Fock Approximation. Int. J. Quantum Chem. 1976, 10, 325–340.

(8)

Wagner, J. P.; Schreiner, P. R. London Dispersion in Molecular Chemistry-Reconsidering Steric Effects. Angew. Chemie Int. Ed. 2015, 54, 12274–12296.

(9)

Liptrot, D. J.; Power, P. P. London Dispersion Forces in Sterically Crowded Inorganic and Organometallic Molecules. Nat. Rev. Chem. 2017, 1, 4.

(10)

Rösel, S.; Quanz, H.; Logemann, C.; Becker, J.; Mossou, E.; Cañadillas-Delgado, L.; Caldeweyher, E.; Grimme, S.; Schreiner, P. R. London Dispersion Enables the Shortest Intermolecular Hydrocarbon H···H Contact. J. Am. Chem. Soc. 2017, 139, 7428–7431.

(11)

Tasinato, N.; Grimme, S. Unveiling the Non-Covalent Interactions of Molecular Homodimers by Dispersion-Corrected DFT Calculations and Collision-Induced Broadening of RoVibrational Transitions: Application to (CH2F2)2 and (SO2)2. Phys. Chem. Chem. Phys. 2015, 17, 5659–5669.

(12)

Hansen, A.; Bannwarth, C.; Grimme, S.; Petrović, P.; Werlé, C.; Djukic, J.-P. The Thermochemistry of London Dispersion-Driven Transition Metal Reactions: Getting the “Right Answer for the Right Reason.” ChemistryOpen 2014, 3, 177–189.

(13)

Ehrlich, S.; Bettinger, H. F.; Grimme, S. Dispersion-Driven Conformational Isomerism in σBonded Dimers of Larger Acenes. Angew. Chemie Int. Ed. 2013, 52, 10892–10895.

(14)

Fokin, A. A.; Zhuk, T. S.; Blomeyer, S.; Pérez, C.; Chernish, L. V.; Pashenko, A. E.; Antony, J.; Vishnevskiy, Y. V.; Berger, R. J. F.; Grimme, S.; Logemann, C.; Schnell, M.; Mitzel, N. W.; Schreiner, P. R. Intramolecular London Dispersion Interaction Effects on Gas-Phase and Solid-State Structures of Diamondoid Dimers. J. Am. Chem. Soc. 2017, 139, 16696–16707.

(15)

Grimme, S.; Steinmetz, M. Effects of London Dispersion Correction in Density Functional Theory on the Structures of Organic Molecules in the Gas Phase. Phys. Chem. Chem. Phys. 2013, 15, 16031–16042.

(16)

Rösel, S.; Balestrieri, C.; Schreiner, P. R. Sizing the Role of London Dispersion in the Dissociation of All-Meta Tert-Butyl Hexaphenylethane. Chem. Sci. 2017, 8, 405–410.

(17)

Schreiner, P. R.; Chernish, L. V.; Gunchenko, P. A.; Tikhonchuk, E. Y.; Hausmann, H.; Serafin, M.; Schlecht, S.; Dahl, J. E. P.; Carlson, R. M. K.; Fokin, A. A. Overcoming Lability of Extremely Long Alkane Carbon–carbon Bonds through Dispersion Forces. Nature 2011, 477, 308–311. ACS Paragon Plus 33 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(18)

Page 34 of 39

Grimme, S.; Djukic, J.-P. Cation−Cation “Attraction”: When London Dispersion Attraction Wins over Coulomb Repulsion. Inorg. Chem. 2011, 50, 2619–2628.

(19)

Fanfrlík, J.; Pecina, A.; Řezáč, J.; Sedlak, R.; Hnyk, D.; Lepšík, M.; Hobza, P. B–H---π: A Nonclassical Hydrogen Bond or Dispersion Contact? Phys. Chem. Chem. Phys. 2017, 19, 18194–18200.

(20)

Wagner, J. P.; Schreiner, P. R. London Dispersion Decisively Contributes to the Thermodynamic Stability of Bulky NHC-Coordinated Main Group Compounds. J. Chem. Theory Comput. 2016, 12, 231–237.

(21)

Lyngvi, E.; Sanhueza, I. A.; Schoenebeck, F. Dispersion Makes the Difference: Bisligated Transition States Found for the Oxidative Addition of Pd(PtBu3)2 to Ar-OSO2R and Dispersion-Controlled Chemoselectivity in Reactions with Pd[P(iPr)(tBu2)]2. Organometallics 2015, 34, 805–812.

(22)

Guo, J.-D.; Nagase, S.; Power, P. P. Dispersion Force Effects on the Dissociation of “Jack-inthe-Box” Diphosphanes and Diarsanes. Organometallics 2015, 34, 2028–2033.

(23)

Su, P.; Li, H. Energy Decomposition Analysis of Covalent Bonds and Intermolecular Interactions. J. Chem. Phys. 2009, 131, 14102.

(24)

Casalz-Sainz, J. L.; Guevara-Vela, J. M.; Francisco, E.; Rocha-Rinza, T.; Martín Pendás, Á. Where Does Electron Correlation Lie? Some Answers from a Real Space Partition. ChemPhysChem 2017, 18, 3553–3561.

(25)

Szczȩśniak, M. M.; Scheiner, S. Studies of Dispersion Energy in Hydrogen‐bonded Systems. H2O–HOH, H2O–HF, H3N–HF, HF–HF. J. Chem. Phys. 1984, 80, 1535–1542.

(26)

He, X.; Fusti-Molnar, L.; Cui, G.; Merz, K. M. Importance of Dispersion and Electron Correlation in Ab Initio Protein Folding. J. Phys. Chem. B 2009, 113, 5290–5300.

(27)

Reinhardt, P.; Piquemal, J.-P. New Intermolecular Benchmark Calculations on the Water Dimer: SAPT and Supermolecular Post-Hartree-Fock Approaches. Int. J. Quantum Chem. 2009, 109, 3259–3267.

(28)

Riley, K. E.; Hobza, P. The Relative Roles of Electrostatics and Dispersion in the Stabilization of Halogen Bonds. Phys. Chem. Chem. Phys. 2013, 15, 17742–17751.

(29)

Gonthier, J. . F.; Thirman, J.; Head-Gordon, M. Understanding Non-Covalent Interactions: Correlated Energy Decomposition Analysis and Applications to Halogen Bonding. Chimia (Aarau). 2017, 72, 193–198.

(30)

Thirman, J.; H.-Gordon, M. Electrostatic Domination of the Effect of Electron Correlation in Intermolecular Interactions. J. Phys. Chem. Lett. 2014, 5, 1380–1385.

(31)

Podeszwa, R.; Szalewicz, K. Accurate Interaction Energies for Argon, Krypton, and Benzene Dimers from Perturbation Theory Based on the Kohn–Sham Model. Chem. Phys. Lett. 2005, 412, 488–493. ACS Paragon Plus 34 Environment

Page 35 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(32)

Journal of Chemical Theory and Computation

Grandinetti, F. Gas-Phase Ion Chemistry of the Noble Gases: Recent Advances and Future Perspectives. Eur. J. Mass Spectrom. 2011, 17, 423–463.

(33)

El Khatib, M.; Bendazzoli, G. L.; Evangelisti, S.; Helal, W.; Leininger, T.; Tenti, L.; Angeli, C. Beryllium Dimer: A Bond Based on Non-Dynamical Correlation. J. Phys. Chem. A 2014, 118, 6664–6673.

(34)

Schneider, W. B.; Bistoni, G.; Sparta, M.; Saitow, M.; Riplinger, C.; Auer, A. A.; Neese, F. Decomposition of Intermolecular Interaction Energies within the Local Pair Natural Orbital Coupled Cluster Framework. J. Chem. Theory Comput. 2016, 12, 4778–4792.

(35)

Neese, F.; Hansen, A.; Liakos, D. G. Efficient and Accurate Approximations to the Local Coupled Cluster Singles Doubles Method Using a Truncated Pair Natural Orbital Basis. J. Chem. Phys. 2009, 131, 64103.

(36)

Neese, F.; Hansen, A.; Wennmohs, F.; Grimme, S. Accurate Theoretical Chemistry with Coupled Pair Models. Acc. Chem. Res. 2009, 42, 641–648.

(37)

Neese, F.; Wennmohs, F.; Hansen, A. Efficient and Accurate Local Approximations to Coupled-Electron Pair Approaches: An Attempt to Revive the Pair Natural Orbital Method. J. Chem. Phys. 2009, 130, 114108.

(38)

Hansen, A.; Liakos, D. G.; Neese, F. Efficient and Accurate Local Single Reference Correlation Methods for High-Spin Open-Shell Molecules Using Pair Natural Orbitals. J. Chem. Phys. 2011, 135, 214102.

(39)

Huntington, L. M. J.; Hansen, A.; Neese, F.; Nooijen, M. Accurate Thermochemistry from a Parameterized Coupled-Cluster Singles and Doubles Model and a Local Pair Natural Orbital Based Implementation for Applications to Larger Systems. J. Chem. Phys. 2012, 136, 64101.

(40)

Riplinger, C.; Neese, F. An Efficient and near Linear Scaling Pair Natural Orbital Based Local Coupled Cluster Method. J. Chem. Phys. 2013, 138, 34106.

(41)

Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural Triple Excitations in Local Coupled Cluster Calculations with Pair Natural Orbitals. J. Chem. Phys. 2013, 139, 134101.

(42)

Liakos, D. G.; Hansen, A.; Neese, F. Weak Molecular Interactions Studied with Parallel Implementations of the Local Pair Natural Orbital Coupled Pair and Coupled Cluster Methods. J. Chem. Theory Comput. 2011, 7, 76–87.

(43)

Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse maps—A Systematic Infrastructure for Reduced-Scaling Electronic Structure Methods. II. Linear Scaling Domain Based Pair Natural Orbital Coupled Cluster Theory. J. Chem. Phys. 2016, 144, 24109.

(44)

Liakos, D. G.; Sparta, M.; Kesharwani, M. K.; Martin, J. M. L.; Neese, F. Exploring the Accuracy Limits of Local Pair Natural Orbital Coupled-Cluster Theory. J. Chem. Theory Comput. 2015, 11, 1525–1539.

ACS Paragon Plus 35 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(45)

Page 36 of 39

Altun, A.; Neese, F.; Bistoni, G. Local Energy Decomposition Analysis of Hydrogen-Bonded Dimers within a Domain-Based Pair Natural Orbital Coupled Cluster Study. Beilstein J. Org. Chem. 2018, 14, 919–929.

(46)

Lu, Q.; Neese, F.; Bistoni, G. London Dispersion Drives the Formation of Agostic Structures. Angew. Chemie Int. Ed. 2018, 57, 4760–4764.

(47)

Bistoni, G.; Auer, A. A.; Neese, F. Understanding the Role of Dispersion in Frustrated Lewis Pairs and Classical Lewis Adducts: A Domain-Based Local Pair Natural Orbital Coupled Cluster Study. Chem. - A Eur. J. 2017, 23, 865–873.

(48)

Neese, F. The ORCA Program System. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 73–78.

(49)

Neese, F. Software Update: The ORCA Program System, Version 4.0. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2018, 8, e1327.

(50)

Mitoraj, M.; Michalak, A. Natural Orbitals for Chemical Valence as Descriptors of Chemical Bonding in Transition Metal Complexes. J. Mol. Model. 2007, 13, 347–355.

(51)

Ziegler, T.; Rauk, A. On the Calculation of Bonding Energies by the Hartree Fock Slater Method. Theor. Chim. Acta 1977, 46, 1–10.

(52)

Mitoraj, M. P.; Michalak, A.; Ziegler, T. A Combined Charge and Energy Decomposition Scheme for Bond Analysis. J. Chem. Theory Comput. 2009, 5, 962–975.

(53)

Non-Covalent Interactions in Quantum Chemistry and Physics: Theory and Applications; de la Roza, A. O., DiLabio, G., Eds.; Elsevier, 2017.

(54)

Boys, S. F. Construction of Some Molecular Orbitals to Be Approximately Invariant for Changes from One Molecule to Another. Rev. Mod. Phys. 1960, 32, 296–299.

(55)

Pipek, J.; Mezey, P. G. A Fast Intrinsic Localization Procedure Applicable for Ab Initio and Semiempirical Linear Combination of Atomic Orbital Wave Functions. J. Chem. Phys. 1989, 90, 4916–4926.

(56)

Bistoni, G.; Rampino, S.; Tarantelli, F.; Belpassi, L. Charge-Displacement Analysis via Natural Orbitals for Chemical Valence: Charge Transfer Effects in Coordination Chemistry. J. Chem. Phys. 2015, 142, 84112.

(57)

Horn, P. R.; H.-Gordon, M. Alternative Definitions of the Frozen Energy in Energy Decomposition Analysis of Density Functional Theory Calculations. J. Chem. Phys. 2016, 144, 84118.

(58)

Mao, Y.; Ge, Q.; Horn, P. R.; Head-Gordon, M. On the Computational Characterization of Charge-Transfer Effects in Noncovalently Bound Molecular Complexes. J. Chem. Theory Comput. 2018, 14, 2401–2417.

(59)

Wuttke, A.; Mata, R. A. Visualizing Dispersion Interactions through the Use of Local Orbital Spaces. J. Comput. Chem. 2017, 38, 15–23. ACS Paragon Plus 36 Environment

Page 37 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(60)

London, F. The General Theory of Molecular Forces. Trans. Faraday Soc. 1937, 33, 8b–26.

(61)

London, F. Zur Theorie Und Systematik Der Molekularkräfte. Zeitschrift für Phys. 1930, 63, 245–279.

(62)

Davison, W. D. Atomic Dipole-Quadrupole Dispersion Forces. J. Phys. B At. Mol. Phys. 1968, 1, 139–144.

(63)

Hornig, J. F.; Hirschfelder, J. O. London Dispersion Forces between Unlike Molecules. J. Chem. Phys. 1952, 20, 1812–1812.

(64)

Newell, A. C.; Baird, R. C. Absolute Determination of Refractive Indices of Gases at 47.7 Gigahertz. J. Appl. Phys. 1965, 36, 3751–3759.

(65)

Pouchan, C.; Bishop, D. M. Static Dipole Polarizability of the Lithium Atom, Cation, and Anion. Phys. Rev. A 1984, 29, 1–5.

(66)

Moore, C. E. Ionization Potentials and Ionization Limits Derived from the Analyses of Optical Spectra. NSRDS-NBS 34, US Department of Commerce, Washington, DC, USA. 1970.

(67)

Miller, M. K. Atom Probe Tomography : Analysis at the Atomic Level; Kluwer Academic / Plenum Publishers, 2000.

(68)

Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: Optimized Auxiliary Basis Sets and Demonstration of Efficiency. Chem. Phys. Lett. 1998, 294, 143–152.

(69)

Eichkorn, K.; Treutler, O.; Öhm, H.; Häser, M.; Ahlrichs, R. Auxiliary Basis Sets to Approximate Coulomb Potentials. Chem. Phys. Lett. 1995, 240, 283–290.

(70)

Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007–1023.

(71)

Balabanov, N. B.; Peterson, K. A. Systematically Convergent Basis Sets for Transition Metals. I. All-Electron Correlation Consistent Basis Sets for the 3d Elements Sc–Zn. J. Chem. Phys. 2005, 123, 64107.

(72)

Peterson, K. A.; Dunning, T. H. Accurate Correlation Consistent Basis Sets for Molecular Core–valence Correlation Effects: The Second Row Atoms Al–Ar, and the First Row Atoms B–Ne Revisited. J. Chem. Phys. 2002, 117, 10548–10560.

(73)

Woon, D. E.; Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. IV. Calculation of Static Electrical Response Properties. J. Chem. Phys. 1994, 100, 2975–2988.

(74)

Herman, P. R.; LaRocque, P. E.; Stoicheff, B. P. Vacuum Ultraviolet Laser Spectroscopy. V. Rovibronic Spectra of Ar2 and Constants of the Ground and Excited States. J. Chem. Phys. 1988, 89, 4535–4549.

(75)

Merritt, J. M.; Bondybey, V. E.; Heaven, M. C. Beryllium Dimer-Caught in the Act of Bonding. Science 2009, 324, 1548–1551. ACS Paragon Plus 37 Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(76)

Page 38 of 39

Odutola, J. A.; Dyke, T. R. Partially Deuterated Water Dimers: Microwave Spectra and Structure. J. Chem. Phys. 1980, 72, 5062–5070.

(77)

Klopper, W.; M. van Duijneveldt-van de Rijdt, J. G. C.; van Duijneveldt, F. B. Computational Determination of Equilibrium Geometry and Dissociation Energy of the Water Dimer. Phys. Chem. Chem. Phys. 2000, 2, 2227–2234.

(78)

Howard, J. C.; Gray, J. L.; Hardwick, A. J.; Nguyen, L. T.; Tschumper, G. S. Getting down to the Fundamentals of Hydrogen Bonding: Anharmonic Vibrational Frequencies of (HF)2 and (H2O)2 from Ab Initio Electronic Structure Computations. J. Chem. Theory Comput. 2014, 10, 5426–5435.

(79)

Klopper, W.; Quack, M.; Suhm, M. A. HF Dimer: Empirically Refined Analytical Potential Energy and Dipole Hypersurfaces from Ab Initio Calculations. J. Chem. Phys. 1998, 108, 10096–10115.

(80)

Řezáč, J.; Hobza, P. Ab Initio Quantum Mechanical Description of Noncovalent Interactions at Its Limits: Approaching the Experimental Dissociation Energy of the HF Dimer. J. Chem. Theory Comput. 2014, 10, 3066–3073.

(81)

Samanta, A. K.; Czakó, G.; Wang, Y.; Mancini, J. S.; Bowman, J. M.; Reisler, H. Experimental and Theoretical Investigations of Energy Transfer and Hydrogen-Bond Breaking in Small Water and HCl Clusters. Acc. Chem. Res. 2014, 47, 2700–2709.

(82)

R.-Casterline, B. E.; Ch’ng, L. C.; Mollner, A. K.; Reisler, H. Determination of the Bond Dissociation Energy (D0) of the Water Dimer, (H2O)2, by Velocity Map Imaging. J. Chem. Phys. 2011, 134, 211101.

(83)

Datta, D.; Kossmann, S.; Neese, F. Analytic Energy Derivatives for the Calculation of the First-Order Molecular Properties Using the Domain-Based Local Pair-Natural Orbital Coupled-Cluster Theory. J. Chem. Phys. 2016, 145, 114101.

(84)

Ruedenberg, K. The Physical Nature of the Chemical Bond. Rev. Mod. Phys. 1962, 34, 326– 376.

(85)

Sosa, C.; Noga, J.; Bartlett, R. J. A Study of the Be2 Potential Curve Using the Full (CCSDT) Coupled-Cluster Method: The Importance of T4 Clusters. J. Chem. Phys. 1988, 88, 5974– 5976.

ACS Paragon Plus 38 Environment

Page 39 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

TOC

ACS Paragon Plus 39 Environment