Decomposition of Intermolecular Interaction Energies within the Local

Article pubs.acs.org/JCTC

Decomposition of Intermolecular Interaction Energies within the Local Pair Natural Orbital Coupled Cluster Framework Wolfgang B. Schneider, Giovanni Bistoni, Manuel Sparta, Masaaki Saitow, Christoph Riplinger, Alexander A. Auer, and Frank Neese* Max Planck Institute for Chemical Energy Conversion, Stiftstrasse 34-36, D-45470 Mülheim an der Ruhr, Germany S Supporting Information *

ABSTRACT: The local coupled cluster method DLPNO-CCSD(T) allows calculations on systems containing hundreds of atoms to be performed while typically reproducing canonical CCSD(T) energies with chemical accuracy. In this work, we present a scheme for decomposing the DLPNO-CCSD(T) interaction energy between two molecules into physical meaningful contributions, providing a quantification of the most important components of the chemical interaction. The method, called Local Energy Decomposition (LED), is straightforward and requires negligible additional computing time. Both the Hartree−Fock and the correlation energy are decomposed into contributions from localized or pairs of localized occupied orbitals. Assigning these localized orbitals to fragments allows one to differentiate between intra- and intermolecular contributions to the interaction energy. Accordingly, the interaction energy can be decomposed into electronic promotion, electrostatic, exchange, dynamic charge polarization, and dispersion contributions. The LED scheme is applied to a number of test cases ranging from weakly, dispersively bound complexes to systems with strong ionic interactions. The dependence of the results on the one-particle basis set and various technical aspects, such as the localization scheme, are carefully studied in order to ensure that the results do not suffer from technical artifacts. A numerical comparison between the DLPNOCCSD(T)/LED and the popular symmetry adapted perturbation theory (DFT-SAPT) is made, and the limitations of the proposed scheme are discussed.

■

INTRODUCTION In recent years the importance of noncovalent intermolecular interactions in chemistry has been clearly recognized and has risen to the forefront of chemical research. At the same time computational chemistry is an indispensable tool for their study and rationalization.1 Different than covalent bonds with high directionality and binding energy due to shared electrons and different than ionic bonds without directionality but high binding energy due to electrostatic attraction between two oppositely charged ions, weak or noncovalent interactions are generally weak, with relatively low directionality. Typical families of weak interactions are pure dispersion interaction, hydrogen bonds and halogen bonds, differing in the degree of directionality and involved binding partners. The interaction energy of weakly bonded, small systems ranges usually between a couple of kcal/mol for dispersion bonded systems up to dozens of kcal/mol between molecules with high polarity. Two types of approaches − perturbative and supermolecular − are usually applied to the calculation of weak interaction energies.2 Within a perturbative approach, weak interactions are computed by partitioning the Hamiltonian into contributions of noninteracting fragment terms and a series of perturbating potentials.3 Among the perturbative approaches, the most successful method is the symmetry-adapted perturbation theory (SAPT).4 Within SAPT, one obtains not only the total interaction energy but also its constituent electrostatic, © 2016 American Chemical Society

induction, dispersion, and exchange-repulsion contributions. These quantities are of great value for the rationalization of the underlying mechanism that gives rise to the interaction. However, the most accurate implementations (SAPT including high order terms of the perturbation series computed with accurate ab initio methods) are computationally demanding (O(N7) scaling), and this limits the applicability of these methods to relatively small systems and benchmark studies.5 For the sake of completeness it should be mentioned that several techniques known to speed up ab initio calculations can be implemented in the framework of SAPT, for example density fitting, making the methods more efficient.6,7 One of the most common versions is SAPT in conjunction with Density Functional Theory (DFT).8 Within a supermolecular approach, the interaction energy is evaluated as the difference between the total energy of the interacting adduct and the energy of the separated and geometrically relaxed monomers.2 Virtually any quantum mechanical method can be employed for the computation of the individual energies; however, the accuracy and applicability of the quantum mechanical models need to be carefully addressed.2 For example, it is commonly known that the dominant theoretical methodology for quantum chemical Received: May 20, 2016 Published: August 26, 2016 4778

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation

prefactor to be routinely applied to molecules with a few hundred atoms (the largest system treated featured more than 20000 basis functions and more than 1000 atoms).54 Hence, DLPNO-CCSD(T) single point energies can be obtained in competitive turnaround times for essentially any system for which a single point DFT energy can be calculated (for comparison between DLPNO-CCSD(T) and DFT benchmarks see ref 55). The advent of local correlation theories has hugely increased the range of applicability of coupled-cluster methods; however, it is important to recognize that the calculation of accurate interaction energies addresses only one aspect of the problem. A deeper understanding of the underlying phenomena is needed to rationalize London dispersion interactions. In this respect SAPT retains a key advantage over supermolecular approaches: SAPT provides a physically motivated decomposition of the interaction energy into its electrostatic, exchange-repulsion, induction, and dispersion terms. Although such a decomposition does not have a counterpart within the supermolecular approaches, some options are available, starting from the pioneering Energy Decomposition Analysis (EDA) of Morokuma,56 that is, however, limited to the single determinant case. Another possibility is the use of Natural Energy Decomposition Analysis (NEDA)57 implemented in the popular framework of the Natural Bond Orbitals (NBO) of Weinhold and co-workers.58 The advantage of NBO and similar a posteriori analysis schemes is that it is based on the electron density and therefore refined ab initio correlated densities can be used as input. A different approach, in which the decomposition of the interaction energy was derived directly from the MP2 pair correlation energy, rather than from the total electron density, was proposed by Grimme and coworkers.59 This analysis specifically targets noncovalent interactions, and it is exclusively based on the nature and symmetry of the localized internal orbitals. This limits the insight offered by the decomposition, as it does not discriminate between charge transfer and dispersion contributions of the correlation energy.59 Conversely, the energy decomposition proposed by Schü tz and Werner takes advantage of the local character of both occupied and virtual orbitals to partition the contribution to the correlation energy by distinct excitation classes.60 This allows for addressing the relative contribution of electrostatic and exchange-dispersion terms. Another technique developed by Head-Gordon et al. relies on the use of absolutely localized molecular orbitals (ALMOs) for the decomposition of the CCSD61 or MP262 interaction energies into chemically meaningful contributions. To the best of our knowledge, this technique has only been applied to the study of weakly interacting model systems. In this work we present a new energy decomposition scheme for the entire interaction energy, including the self-consistent field reference energy. The technique is equally appropriate for weak and strong interactions of molecules and is of a nonperturbative nature. The scheme is based on the DLPNO-approach for CCSD(T), and we refer to it as Local Energy Decomposition (LED). A test set of weakly interacting molecules with different bonding characters is used to assess the merits of the approach. Furthermore, the basis set dependence as well as the dependence of the LED results on the localization scheme is carefully investigated in order to exclude the possibility that LED results are plagued by numerical artifacts. Finally, a comparison to the well-established DFT-SAPT approach is carried out.

applications, namely DFT, does not properly describe weak interactions.9,10 In recent years, the two most popular approaches to deal with this shortcoming have been advocated by Grimme and co-workers11−14 and Truhlar and co-workers.15−17 The procedure known as empirical dispersion correction (DFT-D) relies on adding empirical force field type terms to existing standard density functionals.11−13 The magnitude of such dispersion corrections can be used as an estimate for the dispersive energy contribution to the interaction, but results depend on the different parametrizations used for the damping functions. Truhlar’s approach consists of including the dispersion energy into the extensive parametrization of gradient corrected (GGA), hybrid and metaGGA functionals.18−20 In this case, the distinction between dispersive and nondispersive energy components is not possible. Extensive studies by Hobza and co-workers,1,21,22 Sherrill and co-workers,23,24 and Pulay and co-workers25,26 among many others demonstrated that accurate results are obtained with wave function based ab initio methods if one pushes the treatment of dynamic correlation far enough. These approaches have the desirable feature of converging toward the exact solution of the many-particle Schrödinger equation and hence do not only yield excellent energies but also feature coconvergence of all molecular properties toward their exact values. Hence, specialized or empirical treatments are unnecessary since all relevant interactions (exchange-repulsion, electrostatics, and dispersion) are automatically included in the calculation.8,27−29 The drawbacks of the ab initio approaches are evidently their exploding computational cost (conventional coupled-cluster theory with single-, double-, and perturbatively treated triple excitations,30 CCSD(T), scales as O(N7)) with respect to increasing system size. Second, it might be argued that there is a lack of transparency in the physical interpretation of the many electron wave function. Hence, it would be desirable to have a first-principle method available that is computationally efficient enough for application to large systems with hundreds of atoms while still allowing for a transparent and physically sound interpretation of the results. To improve computational efficiency, several groups have developed strategies and approximations that exploit the locality of the electron correlation to expedite the method while trying to retain its accuracy and reliability.31−39 In a series of papers we have contributed with the development and implementation of the Domain-based Local Pair Natural Orbital CCSD(T) (DLPNO-CCSD(T)).40−48 The DLPNO approach is built upon the pioneering contributions of Meyer, Kutzelnigg, Ahlrichs, Staemmler, and co-workers49−53 and may be characterized by a) the localization of the occupied orbitals: this allows the electron pairs screening based on their contribution to the total correlation energy since these fall off very quickly (R−6) with the distance between the localized occupied orbitals, b) a compact virtual space, that is constructed by using approximate Natural Orbitals for each of the surviving electron pairs (PNOs),40−48 and c) the expansion of the PNOs in terms of projected atomic orbitals (PAOs) within a large domain of atoms that is assigned to each electron pair. In its most recent realization54 it has been demonstrated that a) DLPNO-CCSD(T) is linear scaling with respect to system size, b) it typically recovers about 99.9% of the basis set correlation energy relative to canonical CCSD(T) which translates to an accuracy between 0.25−1 kcal/mol for medium sized systems, and c) it has a sufficiently small 4779

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation

1. THEORETICAL BACKGROUND 1.1. DLPNO-CCSD(T). We recently described the theory and linear scaling implementation of the domain-based pair natural orbital coupled cluster (DLPNO-CCSD(T)) method.47,48,54 The most important features of this technique are summarized below. In the DLPNO-CCSD(T) approach the occupied orbitals of the reference Hartree−Fock (HF) determinant are initially localized, usually by means of either the Foster-Boys63 or the Pipek-Mezey64 scheme. The correlation energy (EC) is then expressed as a sum of electron pair correlation energies εij, where i and j label the occupied HF orbitals. EC = 2 ∑ Fi , aiĩ taiiĩ + i , a∼ii

kcal/mol for LoosePNO, NormalPNO, and TightPNO, respectively.66 Consistent with our previous observations, we recommend TightPNO to be used in conjunction with the LED analysis to ensure that the electron pairs that dominate the interaction energy are treated at the CCSD level. 1.2. Local Energy Decomposition (LED). Within a supermolecular approach, the dissociation energy of two interacting molecules (X and Y) is computed as the difference between the total energy of the adduct (XY) and the sum of the energies of the infinitely separated X and Y. To account for the shortcomings of the finite basis set, this quantity can be corrected by the counterpoise correction of Boys and Bernardi.67 Hence, the total dissociation energy ΔE can be expressed as

∑ εij + E(T) (1)

i≥j

XY ΔE = EXY (XY ) − [EXX (X ) + EYY (Y )]

The first term represents the contribution from the single excitations and vanishes if the Brillouin’s theorem is satisfied. In the PNO approximation, the expression for the pair energy is

− [EXXY (XY ) − EXXY (X ) + EYXY (XY ) − EYXY (Y )]

Here, the following notation is used: EBA(C) denotes the energy of A calculated at the equilibrium geometry of B with the basis X Y set of system C. The first term EXY XY(XY) − [EX(X) + EY(Y)] provides the uncorrected dissociation energy and features the optimized geometries of the isolated molecules. The second XY XY XY term [EXY X (XY) − EX (X) + EY (XY) − EY (Y)] describes the counterpoise correction for basis set superposition error (BSSE) effects, where the monomers are calculated at the dimer equilibrium geometry with the dimer basis set. For the sake of the analysis, we will assume that all energy terms related to the monomers have been appropriately corrected. While the counterpoise correction is designed to compensate for basis set incompleteness effects (and hence vanishes at the basis set limit), it is still advisible to perform the calculations with an as large and flexible Gaussian basis set as possible in order to minimize possible artifacts.68 A rearrangement of the energy terms in eq 3 leads to the following partitioning of the dissociation energy:

εij = ∑ (iaij̃ |jbij̃ )τãij̃ b ̃

ij ij

∼ a∼ijbij

(2)

where ãij and bij̃ are PNOs that belong to pair ij, (iãij|jb̃ij) are the two electron integrals in Mulliken (11|22) notation, while τaij̃ b ̃ = taij̃ b ̃ + taiij̃ t bj̃ are the cluster amplitudes in the PNO basis ij ij

ij ij

ij

2 (2τaij̃ b ̃ 1 + δij ij ij ij ij ij j i (ta ̃ b ̃ , taij̃ and t b ̃ denote ij ij ij

and τãij̃ b ̃ =

(3)

− τaij̃ b ̃ ) are contravariant amplitudes ij ij

the doubles and the singles amplitudes

of the coupled cluster equations). In eq 1, the singles amplitudes are expanded by the singles PNOs which are essentially a set of diagonal PNOs truncated by a much tighter threshold than the doubles PNOs.47,65 It is readily shown that pair correlation energies decay very quickly with the distance Rij between the centroids of the localized occupied orbitals i and j (i.e., at least as R−6 ij ). This permits the exclusion of a large number of “weak pairs” with negligible contributions to the correlation energy from the CCSD treatment. In the DLPNO-CCSD scheme, these pairs are treated with second-order perturbation theory.54 Their contribution to the overall correlation energy is denoted here as EC−WP. The remaining pairs, in the following called “strong pairs”, are treated at the coupled cluster level, with the external correlation space spanned by a highly compact set of pair natural orbitals (PNOs) that are local and different for each electron pair. The energy contribution associated with the strong pairs EC−SP typically represents the dominant part of the overall correlation energy.47 The last term of the correlation energy is associated with the perturbative triples correction EC−(T), as described in ref 54. The occupied and virtual orbitals, to be included in the different steps of the correlation treatment, are selected on the basis of a series of truncation thresholds that are described in refs 54 and 66 and whose effect on the accuracy of the calculation has been already carefully investigated. To render the method black box, three sets of thresholds for these parameters have been established.66 In order of increasing accuracy and computational cost, they are referred to as “LoosePNO”, “NormalPNO”, and “TightPNO”. For the S66 benchmark set,21,22 especially devised for benchmarking noncovalent intermolecular interactions, the maximum absolute deviations in the evaluation of the total energy with respect to CCSD(T) reference have been found to be 0.31, 0.26, and 0.10

XY (XY ) − (EXXY (XY ) + EYXY (XY ))] ΔE = [EXY + [(EXXY (X ) − EXX (X )) + (EYXY (Y ) − EYY (Y ))]

≡ ΔEint + ΔEgeo − prep (4)

In this way, the counterpoise corrected dissociation energy decomposes into a genuine “electronic interaction” part (ΔEint) ̀̀ and a geometric preparation (ΔEgeo−prep, also called deformation energy’’). The latter contribution is the energy difference arising from the geometrical distortion of X and Y from their optimal monomer geometries to the geometry adopted at the interaction conformation. In the following, a physical meaningful decomposition of ΔEint within the DLPNO-CCSD(T) scheme is provided. As there is obviously an infinite number of possibilities to decompose the total energy, the decomposition presented here is by no means unique. However, in our opinion, it is physically sound and leads to results that are readily chemically interpretable. For this purpose, all energy contributions in eq 4 are decomposed into their constituent Hartree−Fock (HF) and correlation contributions. Accordingly, the overall electronic interaction can be written as HF C ΔEint = ΔEint + ΔEint

4780

(5) DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation ⎡ ZA(X )ZB(Y ) (X , Y ) EHF = PX ↔ Y ⎢ ∑ −2 ∑ (X ) (Y ) ⎢ iX , AY ⎣ AX > BY |RA − R B |

The first term ΔEHF int accounts for the (static) electrostatic interaction between the monomers, polarization effects, and donor−acceptor interactions. The second term ΔECint (dynamically) corrects all the above-mentioned components and also accounts for the dispersive attraction resulting from instantaneous dipole−dipole interactions, whose precise treatment is crucial for the appropriate treatment of weak interactions. ΔEHF int is decomposed by exploiting the localization of the occupied orbitals (unless otherwise specified, Foster-Boys localization has been used for the occupied orbitals). Assuming a closed-shell system, the HF energy can be decomposed in terms of the HF orbitals {ψ} as

ZAZB + |RA − R B|

∑ A BY |RA − R B |

(X ) (X ) (X ) (X ) = V NN + T (X ) + V eN + V (J X ) + V exch EHF

∑ AX > BX

−2

ZA(X )ZB(X ) (X ) |RA − R(BX )|

∑ iX , AX

+

∑ iX ≥ jX

iX

+

∑ ⟨iX | − ∇i2X |iX ⟩ iX

ZA(X ) |r(i X ) − R(AX )|

4 (iX iX |jX jX ) − 1 + δij

−2

iX

∑ iX ≥ jX

(9)

Note that the expression of the kinetic energy only appears in the intramolecular term. Obviously, in the limit of noninteracting molecules, the (X) and (Y) terms equal the quantities obtained for the isolated fragments and the intermolecular contribution E(X,Y) HF vanishes. Writing ΔEHF int into a sum of partitioned energies and regrouping based on the affiliation to either of the fragment leads to (6)

=

⎦

iX , jY

⎡ ⎤ PX ↔ Y ⎢4 ∑ (iX iX |jY jY )⎥ = 4 ∑ (iX iX |jY jY ) + 4 ∑ (iY iY |jX jX ) ⎢ ⎥ iX > jY iY > jX ⎣ iX > jY ⎦ (10)

4 (ii|jj) 1 + δij

XY (X ) (Y ) (X , Y ) EXY (XY )HF = EHF + EHF + EHF

iX

where PX↔Y represents the symmetrization operator in terms of X and Y. An example is given for the Coulomb integrals in eq 10.

i

+

− R(AY )|

∑ (iXiX|jY jY ) − 2 ∑ (iXjY |iXjY )⎥⎥

iX , jY

∑ ⟨i| − ∇i2 |i⟩

ZA i − 2∑ i r | − i RA| i,A 2 (ij|ij) −∑ 1 δij + i≤j

ZA(Y ) |r(i X )

⎤

+4

EHF = VNN + Te + VeN + VJ + Vexch =

iX

∑ iX , AY

2 (iX j |iX j ) 1 + δij X X

iX

ZA(Y ) |r(i X ) − R(AY )|

iX

+4

∑ iX > jY

⎤ (iX iX |jY jY )⎥ ⎥ ⎦ (12)

The first and third terms are always repulsive, while the second term is always attractive. Note that E(X,Y) elstat represents the electrostatic interaction between the distorted electron density distributions of the fragments. This term differs from analogues

(8)

and the intermolecular HF energy is 4781

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation of alternative energy decomposition schemes70,71 in which the electrostatic interaction of the unperturbed fragments is considered. We will come back to this point later in the paper. The exchange term is simply (X , Y ) Eexch

⎡ ⎤ ⎢ = PX ↔ Y −2 ∑ (iX jY |iX jY )⎥ ⎢ ⎥ ⎣ iX > jY ⎦

interfragment correlation energy of the weak pairs E(X,Y) C−WP. By subtracting the weak-pairs contribution of the isolated calculated fragments in the dimer geometry EXY X (XY)C−WP and EYXY(XY)C−WP from the intracontributions, the electronic preparation ΔEC−WP el−prep of the weak pairs is obtained. An analogous procedure is carried out for the triples correction:

(13)

C − (T ) ΔE int = EC(X−,(YT)) + (EC(X−)(T ) − EXXY (XY )C − (T )) + (EC(Y−)(T ) − EYXY (XY )C − (T ))

This intermolecular exchange is always negative and describes a stabilizing component of the interaction, lowering the repulsion between electrons of the same spin. Note that this exchange contribution differs from the one included in the SAPT decomposition that, by construction, is always positive and represents the repulsive component of the interaction (see section 3.5). As mentioned above, this decomposition is by no means unique. Different energy and charge decomposition techniques can be applied in conjunction with the proposed methodology in order to elucidate particular aspects of the chemical interaction. 56,72−76 For instance, in contrast to other approaches,76,77 our scheme does not provide a quantification for the classical induction (also called “polarization’’) and charge transfer contributions but instead provides a different interpretative framework in which the electronic structure of molecules can be discussed. Nevertheless, it is still possible to draw connections between the proposed scheme and previously defined concepts. For instance, in our scheme, a strong induction component would result in both a large electronic preparation (due to the polarization of the fragments) and electrostatic contribution (due to the interaction between permanent and induced multipoles of the distorted electron clouds of the fragments). This approach, besides providing a novel, physically sound, decomposition of ΔEHF int , is particularly useful for the study of noncovalent interactions since, by exploiting the locality of both the internal and the virtual space, it permits to decompose the corresponding correlation contribution ΔECint at the coupled cluster level of theory. As discussed in section 1.1, the correlation energy in DLPNO-CCSD(T) is comprised of three terms: the correlation energy for the strong pairs treated at the coupled cluster level (EC−SP), the correction to account for the truncation of the pair space and truncation of the virtual PNO space computed at the MP2 level, called weak pair energy (EC−WP), and the perturbative triples correction (EC−(T)). In analogy to the decomposition of the HF energy in eq 11, the weak pairs are divided into intra- and intermolecular contributions and the electronic preparation energy is determined by

= EC(X−,(YT)) + ΔEC(X−)(T ) + ΔEC(Y−)(T ) (T ) ≡ EC(X−,(YT)) + ΔEelC−−prep

(15)

If all three considered occupied orbitals are localized on the same fragment, they contribute to the intrafragment correlation (Y) energies E(X) C−(T) and EC−(T). Else they are counted to the (X,Y) interfragment correlation energy of the triples EC−(T) . Accordingly, the electronic preparation of the triples correction ΔEC−(T) el−prep is obtained. For the strong pairs contribution, that usually represents the largest part of the correlation energy, a more sophisticated approach is applied. It takes into account both the locality of the occupied and that of virtual orbitals. EC−SP can be written as a summation of pair energy terms as EC − SP =

∑ ∑ (iaij̃ |jbij̃ )τãij̃ b ̃ ∼ i ≥ j a∼ijbij

ij ij

(16)

Since PNOs are constructed for a distinct pair of localized occupied orbitals, it is reasonable to assume that they are located in the same region of space as the pair of the occupied orbitals they are constructed for. However, it is necessary to also localize the PNOs for each pair to ensure that they can be assigned to a given fragment. Therefore, each PNO is localized by using the Pipek-Mezey scheme and subsequently assigned to a fragment based on Mulliken charge analysis. Afterwards, a transformation of the amplitudes and integrals to the localized PNO-basis is carried out. This procedure is straightforward and analogous to the one followed for the occupied orbitals. Delocalization effects arising from less than perfectly localized PNOs will be discussed below. If both occupied orbitals of an electron pair are located on one fragment, the corresponding PNOs are expected to be localized on the same fragment, too. However, especially for strongly interacting molecules, we have observed localized PNOs, that are located on different fragments than the occupied orbitals they belong to. It is important to point out that the resulting contribution is not just a BSSE effect. If the latter was the case the contributions would vanish at the basis limit. However, we have consistently observed the opposite effect, and these contributions become larger with increasingly large basis sets and still persist at the basis set limit. Based on the localization of the occupied as well as the virtual orbitals, different families of contributions, displayed in Figure 1, can be defined: 1. Intrafragment excitations, with all holes and particles localized on the same fragment. 2. Dispersion contributions, with each fragment possessing one hole and one particle. They can be separated in A) genuine dispersion and B) exchange dispersion, as depictured in Figure 1. 3. Charge transfer contributions with a different number of holes and particles for each fragment. These are A) a dynamic

C − WP XY (XY )C − WP − EXXY (XY )C − WP − EYXY (XY )C − WP ΔE int = EXY Y) = EC(X−,WP + (EC(X−)WP − EXXY (XY )C − WP ) + (EC(Y−)WP − EYXY (XY )C − WP ) Y) = EC(X−,WP + ΔEC(X−)WP + ΔEC(Y−)WP Y) WP ≡ EC(X−,WP + ΔEelC−−prep

(14)

The contribution of the weak pairs of the adduct to the correlation energy EXY XY(XY)C−WP is decomposed depending on which fragment they are localized on. If both occupied orbitals of a weak pair are localized on the same fragment the corresponding energy is accounted to the intracontributions (Y) E(X) C−WP and EC−WP. If the occupied orbitals of the pair are localized on different fragments, they contribute to the 4782

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation ECCT−(SPX → Y ) = + +

∑ ∑ (iXaX̃ |jX bỸ )τãi ̃ jb ̃

XX

∼ iX > jX a∼X bY

X Y

∑ ∑ (iXaỸ |jX bỸ )τãi ̃ jb ̃

XX

∼ iX > jX a∼Y bY

+

Y Y

∑ ∑ (iXaỸ |jY bỸ )τãi ̃ jb ̃

XY

∼ iX > jY a∼Y bY

Y Y

∑ ∑ (iY aỸ |jX bỸ )τãi ̃ jb ̃

YX

∼ iY > jX a∼Y bY

Y Y

(20)

and X ← Y Charge Transfer: ECCT−(SPX ← Y ) = + +

charge polarization with both holes on one fragment and one particle per fragment; B) a double dynamic charge polarization, with both holes on one fragment and both particles on the other; C) another dynamic charge polarization, with one hole per fragment and both particles on the same fragment. For the sake of simplicity, Figure 1 shows only charge transfer excitation from X to Y. The corresponding energy expressions of those contributions are given in the following. Note that the subscript ij is omitted for clarity and that the indices X and Y denote the fragment on which the orbital is located.

(Y ) ECIntra = − SP

∑ ∑ (iXaX̃ |jX bX̃ )τãi ̃ jb ̃

XX

X X

∼ iX > jX a∼X bX

(17)

YY Y Y

+

XY

X Y

XY

Y X

∼ iY > jY a∼X bX

+

X X

∑ ∑ (iY aX̃ |jX bX̃ )τãi ̃ jb ̃

YX

∼ iY > jX a∼X bX

X X

∑ ∑ (iY aX̃ |jY bX̃ )τãi ̃ jb ̃

XY

∼ iX > jY a∼X bX

X X

(21)

Eq 22 is the basis of our Local Energy Decomposition (LED) technique, and its suitability will be illustrated in section 3 with some illustrative examples. C To summarize, ΔEHF el−prep and ΔEel−prep describe the HF and the correlation contributions to the electronic preparation of each of the molecules that are interacting. They are obtained from the difference in the electronic structure of the fragments (X,Y) in the dimer and isolated monomers. The terms E(X,Y) elstat and Eexch describe the intermolecular electrostatic and exchange interactions obtained for HF. The terms EDISP(X,Y) , ECT(X→Y) , and C−SP C−SP CT(X←Y) E C−SP describe the dispersion and charge transfer components of the correlation contribution of the strong (X,Y) pairs to the interaction. The terms E(X,Y) C−WP and EC−(T) provide corrections to the interaction energy for the weak pairs and the triples correction. The final term ΔEgeo−prep accounts for the geometric preparation of the system from the dimer geometry to the separated monomer geometries.

∑ ∑ (iXaỸ |jY bX̃ )τãi ̃ jb ̃ ]

∼ iX > jY a∼Y bX

YY

(22)

∑ (iXaX̃ |jY bỸ )τãi ̃ jb ̃

∼ iX > jY a∼X bY

∑ ∑ (iY aX̃ |jY bX̃ )τãi ̃ jb ̃

Y) + ECCT−(SPX → Y ) + ECCT−(SPX ← Y ) + EC(X−,WP + EC(X−,(YT)) + ΔEgeo − prep

(18)

Consistent with the definition for HF, the weak pairs, and the triples (eq 11, 14, and 15), the electronic preparation energy (Y) ΔE(X) C−SP + ΔEC−SP is obtained by subtracting the strong pairs contribution of the isolated calculated fragments in the dimer geometry from the intrafragment contributions. By summing up the electronic preparation energies of the strong pairs, the weak pairs, and the triples, the electronic preparation of the C−WP C−(T) correlation energy ΔECel−prep = ΔEC−SP el−prep + ΔEel−prep + ΔEel−prep is obtained. The remaining contributions are: Genuine and exchange dispersion: (X , Y ) ECDISP = PX ↔ Y [ ∑ − SP

Y X

DISP(X , Y ) C (X , Y ) (X , Y ) ΔE = ΔEelHF − prep + ΔEel − prep + Eelstat + Eexch + EC − SP

∑ ∑ (iY aỸ |jY bỸ )τãi ̃ jb ̃

∼ iY > jY a∼Y bY

YY

∼ iY > jY a∼Y bX

The Charge Transfer (CT) might be viewed as a correction for the CT and polarization already accounted for at the HF level and corrects the overpolarization that is characteristic of the HF approach. Alternatively, the contributions can be interpreted as a ‘dynamic charge polarization’ given the fact that they arise only at the correlated level, where excitations can be viewed as instantaneous electron jumps. The fundamental mechanism by which these contributions enter to the interaction energy is that the CT event creates a cation/ anion pair that will have an attractive Coulombic interaction between them. In the time dependent interpretation, these dynamic charge polarization events correspond to instantaneous ion pair formations that occur with a small probability that decays exponentially with distance. Hence, the overall distance dependence is not 1/R as would be suggested by the Coulombic attraction. Clearly, excitations of this type are also part of the BSSE for calculations with finite basis sets. However, as the numerical results collected below will demonstrate, these contributions persist at the basis set limit, and hence they are physically meaningful and represent important contributions to the overall interaction energy. Collecting all the terms for both the HF and the correlation energy contributions we obtain for the dissociation energy

Figure 1. Graphical representation of the different groups of excitations for the strong pairs in the DLPNO-CCSD(T) approach. For 3. just charge transfer contributions X → Y are shown, while the corresponding contributions X ← Y are omitted.

(X ) ECIntra = − SP

∑ ∑ (iY aỸ |jY bX̃ )τãi ̃ jb ̃

(19)

X → Y Charge Transfer: 4783

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation

2. COMPUTATIONAL DETAILS All DLPNO-CCSD(T) calculations as well as the geometry optimizations were performed with a development version (based on version 3.0.3) of the ORCA suite of programs.78 All new functionalities described in this article will be made available to the scientific community free of charge with the upcoming (4.0) release of the software package. To calculate interaction energies and geometries the following procedure was applied: 1. The isolated molecules as well as the entire system were optimized using the PBE-exchange correlation functional with the D3-dispersion correction of Grimme79 and the resolution of identity approximation.80 Here, the def2-QZVPD basis set with matching auxiliary basis set was used. 2. DLPNO-CCSD(T) calculations of the isolated molecules and the entire system were carried out with TightPNO settings as described by Liakos et al.,66 using the cc-pVXZ basis sets and the corresponding auxiliary basis sets (X = D,T,Q,5). For the localization of the MOs, the Foster-Boys localization scheme has been applied. The PNOs were localized by using the PipekMezey-localization scheme, which in numerous numerical experiments has been shown to converge much faster than the Foster-Boys scheme. 3. The LED was carried out for the interacting molecules, and the energy terms were subsumed as described in section 1.2 Based on the same geometries that have been used for the DLPNO-CSD(T) calculations, DFT-SAPT has been performed using the Molpro program package, version 2012.1.81 The PBE0AC exchange correlation functional, as implemented in Molpro has been applied. This functional, suggested by Jansen,8 is the PBE082 functional with the asymptotic correction defined in a way that the long-range part of the functional contains 0.75 LB94-exchange.83 The applied shift parameter for the bulk potential within this correction was calculated as Δxc = εHOMO + IP. Here, εHOMO is the energy of the highest occupied molecular orbital of the monomer, obtained from a geometry optimization at the PBE0/aug-cc-pVQZ level of theory. IP is the ionization potential taken from the NIST Standard Reference Database.84

Figure 2. Models of the investigated systems.

are expected in addition to the dispersion contributions. In analogy to the findings of Ř ezáč et al.,22 who, based on SAPT, investigated the dispersion/electrostatics ratio for different types of weak bonds, the relative contribution of the dispersion contribution to the entire interaction energy is expected to decrease with increasing charge transfer and electrostatic contributions. As can be seen in Table 1, the dissociation energy ΔE for the chosen systems ranges from −0.37 kcal/mol for the methane dimer to −7.77 kcal/mol calculated for the ammonia−chlorine fluoride adduct, if the calculations are performed using DLPNO-CCSD(T). These values differ less than 0.50 kcal/ mol from the dissociation energies calculated by canonical CCSD(T), which is in line with previous benchmarks on weakly interacting systems.65,66 Consistent with the binding energy, an increased distortion of the interacting molecules from their monomer equilibrium geometries is obtained for the dimers in their optimized geometry and of course is reflected in an increasing ΔEgeo−prep. Thus, for the methane dimer the geometric preparation energy is almost zero, and it does not exceed 1 kcal/mol for the benzene dimer and the adducts containing water. However, for the strongly interacting NH3−FCl it is a rather large contribution with 4.40 kcal/mol which is counteracting about one-third of the intrinsically high interaction energy of −12.16 kcal/mol. Regarding the partitioning of the interaction energy ΔEint into their Hartree−Fock (HF) and correlation contributions, the dispersion dominated methane and benzene dimers display a repulsive interaction at the HF level, which is counterbalanced by the attractive correlation interaction. For the other systems, the interaction calculated at the HF level already provides the largest contribution to the interaction energy, while the correlation interaction provides a non-negligible contribution to the interaction energy (about 25−40%). Similar to the geometric distortion, the increasing interaction energy goes along with an increase of the electronic distortion at the HF level, ranging from 1.56 to 309.53 kcal/mol. This increase in the monomer HF energy is mostly compensated by (X,Y) the intermolecular HF contributions E(X,Y) elstat and Eexch . At this point it is important to recapitulate that our evaluation of the electrostatic contribution differs from other procedures. While in alternative decomposition methods70,71 the electrostatic interaction of the unperturbed fragments is considered, E(X,Y) elstat describes the electrostatic interaction of the fragments in their

3. RESULTS AND DISCUSSION In this section, the Local Energy Decomposition (LED) is applied to a test set of molecular adducts with prototypical intermolecular interactions. The basis set dependence of the decomposition as well as the influence of the localization scheme of the MOs will be discussed. On a strongly interacting system, the distance behavior of the LED is investigated. Furthermore, the LED is compared to DFT-SAPT. 3.1. LED for Different Types of Interactions. In this section the LED analysis is applied to the interaction of different molecular adducts, which are the methane dimer, water dimer, benzene dimer, ammonia−water, and ammonia− chlorine fluoride in the geometries as shown in Figure 2. These systems have been chosen as representatives of typical weak interactions. The interactions within the methane and benzene dimers are purely of dispersion type, the water dimer and the water−ammonia adduct are representatives for typical hydrogen bonded systems, and the ammonia−chlorine fluoride system exhibits a strong halogen bond.79 The first two systems are expected to display the weakest interactions among the investigated systems. For hydrogen bonded systems as well as in the case of the halogen bond, charge transfer contributions 4784

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

perturbed geometric and electronic structures of the adduct. Nevertheless, the trend observed here nicely follows the chemical intuition. Accordingly, the weakest interaction is reported for the methane dimer, whose constituting monomers feature only a small permanent quadrupole, while the strongest interaction is calculated between ammonia and chlorine fluoride with a strong dipolar interaction and a distinct halogen bond.85 The trend of the overall interaction energy is also reflected in the increasing contribution of the dynamic charge transfer (CT) component to the intermolecular correlation energy, and ECT(X←Y) . As discussed above, these quantities ECT(X→Y) C−SP C−SP may be interpreted as a ‘dynamic charge polarization’ arising from instantaneous intermolecular particle-hole formation. These dynamic CT terms typically present large values if compared with the estimated static CT contributions obtained from other energy decomposition techniques. However, we should stress here that static and dynamic charge transfer contributions have different physical origins and cannot be directly compared. Moreover, it is interesting to note that, for the systems studied here, the sum of the dynamic CT terms is almost equal (in absolute value) to ΔECel−prep. Consequently, the variations observed in the overall correlation contribution along the series correlate with the dispersion energy contribution. Having this in mind, the comparative analysis of ECT(X→Y) and C−SP among the different systems can still be used for ECT(X←Y) C−SP obtaining insights into the different nature of the interaction. For the methane and benzene dimers, these contributions are below 1 kcal/mol and are almost equal for both ways of excitations. The small deviation of 0.03 kcal/mol in the case of the methane dimer can be explained by uncertainties in the localization of the PNOs, which will be discussed in section 3.3. Conversely, for the other investigated systems, a remarkable deviation between the two components is observed, indicating a preferred direction for the dynamic charge transfer contribution (note that these results are independent of the particular TCutPNO threshold used in the calculation, as shown in Table S1 in the Supporting Information). This is always consistent with the direction of the actual charge transfer, as defined by Mulliken charge population. This is remarkable, as the partial charge on each atom is decreased, consistently with the counteracting of the HF overpolarization by dynamic charge polarization.86 However, the amount of charge transfer from one fragment to the other is increased. For illustration, the Mulliken charges per atom and fragment of the water−ammonia adduct are given in Table 2, computed at both the HF and the coupled cluster level of theory. The Mulliken charges of the latter are based on the unrelaxed DLPNO-CCSD density.87 The last component of the interaction that our analysis permits to put into a quantitative comparison is dispersion. In the case of the methane and the benzene dimers, that are not bound at the HF level, this contribution dominates over all correlation contributions and, as a result, can be regarded as the most important component of the interaction. In the remaining systems, despite the fact that charge transfer contribution dominates the intermolecular correlation contributions, dispersion is still far from being negligible and amounts to −1.24 kcal/mol, −1.61 kcal/mol, and −6.25 kcal/mol, respectively for water−water, water−ammonia, and ammonia−chlorine fluoride. This corresponds to one quarter of the total interaction energy for the hydrogen bonded adducts and half of the total

a Dissociation energy calculated by canonical CCSD(T). bDissociation energy calculated by DLPNO-CCSD(T). cThe labels for different contributions are as defined in section 1.2. For the applied geometries, see Figure 2. The values, that are obtained by using the cc-pVQZ basis set and the corresponding auxiliary basis, are given in kcal/mol.

−0.11 −1.16 −0.41 −0.58 −3.15 −0.02 −0.38 −0.06 −0.07 −0.18 −0.18 −0.98 −0.29 −0.64 −6.66 −0.21 −0.98 −8.05 −8.88 −24.57 −0.71 −4.73 −1.24 −1.61 −6.25 −0.40 −3.31 −4.30 −6.63 −45.65 −0.78 −6.37 −25.57 −38.93 −270.92 0.47 2.52 8.89 9.93 35.68 1.56 13.22 26.30 41.23 309.53 −0.76 −5.70 −1.16 −1.84 −5.13 0.38 3.53 −3.57 −4.33 −7.04 −0.38 −2.17 −4.73 −6.17 −12.16 0.00 0.13 0.23 0.50 4.40 −0.37 −2.04 −4.50 −5.66 −7.77 −0.39 −2.15 −4.57 −5.74 −8.14 CH4−CH4 C6H6−C6H6 H2O−H2O NH3−H2O NH3−FCl

ΔEint ΔEgeo−prep ΔEb ΔEa

Table 1. DLPNO-CCSD(T)-LED for Different Moleculesc

ΔEHF int

ΔECint

HF ΔEel−prep

C ΔEel−prep

E(X,Y) elstat

E(X,Y) exch

EDISP(X,Y) C−SP

ECT(X→Y) C−SP

ECT(X←Y) C−SP

E(X,Y) C−WP

(X,Y) EC−(T)

Journal of Chemical Theory and Computation

4785

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation and

Table 2. Mulliken Charges of the Water−Ammonia Adduct Calculated at the HF and the DLPNO-CCSD(T) Level of Theory Using the cc-pVQZ Basis N H H H O H H NH3 H2O

HF

DLPNO-CCSD(T)

−0.462 0.175 0.161 0.161 −0.617 0.332 0.249 0.035 −0.035

−0.457 0.173 0.161 0.161 −0.564 0.298 0.226 0.039 −0.039

(X , Y ) Y) Edisp = ECDISP + EC(X−,WP − SP

The dissociation energy of the above-discussed model systems, according to the suggested summation of eq 23, is shown in Table 3. In addition, the ratio of the dispersion energy Edisp to the entire interaction energy ΔEint (not including the geometric preparation, see eq 4) is given in the last column. Through this ratio, the character of different interactions can be readily characterized regarding the importance of dispersion. A ratio exceeding 1 unambiguously identifies a dispersion dominated interaction. Furthermore, a ratio between 0.5 and 1 points to an interaction with balanced contributions of dispersive and nondispersive character, while a ratio smaller than 0.5 suggests an interaction rather dominated by nondispersive contributions. In summary, the decomposition of the interaction energy for the investigated systems is consistent with the chemical intuition and provides a quantitative framework in which to analyze the interaction energies obtained at the DLPNOCCSD(T) level. The results presented in this section are consistent with chemical intuition for the analyzed interactions. Note that the results are in accordance with former investigations of various authors like Ř ezáč et al.22 or Schütz et al.,60 who report comparable ratios of dispersion interactions, electrostatics, or charge transfer contributions. 3.2. Basis Set Convergence. In the formulation presented in this work, LED is carried out considering counterpoise corrected energies. Due to their definition, introduced in section 1.2, the intermolecular terms are corrected for BSSE. Nevertheless, the LED contributions may be affected by a Basis Set Incompleteness Error (BSIE). To investigate how the BSIE affects the different energy contributions, the convergence of the LED-terms with the size of the basis set has been investigated for the methane dimer and the hydrogen bonded ammonia−water complex. Table 4 shows the basis set convergence behavior of the HFcontributions for these systems. As expected,60,88 the HF interaction energy and its constituting energy components rapidly converge with the basis set, and the qualitative picture emerging from our analysis is the same for all basis sets tested. For the weakly interacting methane dimer the maximal change from double to the quintuple basis of all considered values is 0.03 kcal/mol for E(X,Y) exch . For the strongly interacting water− ammonia system, a slightly worse convergence behavior is obtained for both ΔEHF int and its constituting components, but the results are still very stable. In particular, the maximal change is 0.25 kcal/mol for the ΔEHF el−prep term. Convergence for (X,Y) ΔEHF and E is reached using the quadruple-ζ basis set, el−prep exch while the electrostatic term is still not entirely converged even with the largest basis set employed.

interaction energy for the halogen bonded ammonia−chlorine fluoride system. (X,Y) The components E(X,Y) C−WP and EC−(T) irregularly increase with the binding strength and, if anything, correlate best with the dispersion interaction. At this point it is important to recall that the decomposition for these terms is, unlike for the strong pairs, just based on the occupied orbitals, and hence they are not further decomposed with respect to charge transfer and dispersion. However, in the case of E(X,Y) C−WP the electron pairs are usually far from each other, which reduces the importance of CT interactions. By this and the obtained correlation of E(X,Y) C−WP with the dispersion interaction, it is reasonable to assume that the intermolecular component of the weak pairs energy consists mainly of small dispersion type contributions. While the correlation with the size of dispersion contribution is also calculated for the triples correction, the same assumption cannot be made in general. As the triples are constructed from the same occupied orbitals being involved in the strong pairs, it is reasonable to expect a similar ratio of dispersion and charge transfer components as for the strong pairs. In the interest of simplicity it is reasonable to combine several contributions of the LED scheme. The electronic preparation energies at the HF level and coupled cluster level can be combined. Furthermore, one can define an “electronic interaction” contribution Eelint as the combination of the electrostatic and the CT terms. Finally, the solely dispersive character of the intermolecular weak-pair contribution permits its inclusion into the dispersion contribution. Hence, the equation for the dissociation energy can be written in a more compact form as (X , Y ) ΔE = ΔEgeo − prep + ΔEel − prep + Eelint + Eexch + Edisp + EC(X−,(YT))

(23)

with C ΔEel − prep = ΔEelHF − prep + ΔEel − prep

(24)

(X , Y ) Eelint = Eelstat + ECCT−(SPX → Y ) + ECCT−(SPX ← Y )

(25)

(26)

Table 3. Summary of the Contributions of the LED at the DLPNO-CCSD(T) Level of Theory

CH4−CH4 C6H6−C6H6 H2O−H2O NH3−H2O NH3−FCl

ΔE

ΔEgeo−prep

ΔEel−prep

Eelint

E(X,Y) exch

Edisp

E(X,Y) C−(T)

−0.37 −2.04 −4.50 −5.66 −7.77

0.00 0.13 0.23 0.50 4.40

2.03 15.74 35.19 51.16 345.20

−1.17 −8.33 −33.91 −48.45 −302.13

−0.40 −3.31 −4.30 −6.63 −45.65

−0.73 −5.11 −1.30 −1.67 −6.44

−0.11 −1.16 −0.41 −0.58 −3.15

4786

Edisp ΔEint

1.92 2.35 0.27 0.27 0.53

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation

We conclude that, although the results show pronounced basis set dependence in some cases, the overall trends are robust and the results with respect to the relative importance of dispersion versus charge transfer and HF versus correlation contributions are fairly stable. There are some fluctuations in the relative numbers obtained for each term that arise mostly from imperfect localization of the PNOs (the same behavior is obtained with NormalPNO settings, see Table S2 in the Supporting Information). Despite the fact that different localization schemes are available and can be applied, this behavior is, in our opinion, unavoidable. Nevertheless, these fluctuations are never large enough to change the fundamental conclusions. We emphasize again that there is no unique way to decompose an observable quantity into nonobservable contributions. Hence, no decomposition scheme can claim absolute results. We will revisit the issue of dependence of the results on the localization procedure in the following section. 3.3. Effect of Localization Scheme. By default, the DLPNO procedure relies on the localization of the occupied orbitals by means of the Foster-Boys procedure.63 Alternatively, the Pipek-Mezey scheme64 is also available. As the Pipek-Mezey orbital localization scheme is known to produce less stable results if augmented basis sets are used, the Foster-Boys procedure is suggested for the LED and applied in most of the studies presented in this work. However, as far as chemical interpretation is concerned, the Pipek-Mezey algorithm has the advantage of maintaining the separation of the orbitals according to their σ- and π-symmetry, whereas the former mixes them and returns the well-known “banana bonds”. For this reason, we present in this section a comparative study of how these two localization schemes affect the results obtained by means of our LED analysis. We base the discussion on the strongly interacting water− ammonia and on the benzene dimer case. The latter system is particularly interesting because of the already mentioned differences in the localization of σ and π occupied orbitals, which may affect the results of the decomposition. In Table 6, the Hartree−Fock terms of the LED analysis obtained by employing the two localization schemes are shown. Remarkably enough, both schemes give very similar results for the water− ammonia case (see Tables S3 and S4 in the Supporting Information for the basis set convergence behavior of the different terms with the Pipek-Mezey scheme). Conversely, the differences in the energy terms in the benzene dimer case are comparatively much larger. This is of course a special case due to the fact that the localized orbitals obtained with the two schemes qualitatively differ.

Table 4. Basis Set Dependence of the HF Contributions to the Interaction Energya basis CH4-CH4 cc-pVDZ cc-pVTZ cc-pVQZ cc-pV5Z NH3-H2O cc-pVDZ cc-pVTZ cc-pVQZ cc-pV5Z a

ΔEHF int

ΔEHF el−prep

E(X,Y) elstat

E(X,Y) exch

0.39 0.38 0.38 0.38

1.54 1.57 1.56 1.56

−0.77 −0.79 −0.78 −0.79

−0.37 −0.40 −0.40 −0.40

−4.23 −4.21 −4.33 −4.41

41.49 41.34 41.23 41.24

−38.96 −38.90 −38.93 −39.02

−6.75 −6.65 −6.63 −6.63

The energies are given in kcal/mol.

As expected,88 the correlation terms exhibit a slower convergence behavior compared to the HF components (see Table 5). For the methane dimer the difference between the correlation interaction calculated with the quadruple- and quintuple-ζ basis sets amounts to 0.04 kcal/mol and thus can be considered as converged. For the water−ammonia system this difference is 0.21 kcal/mol, and convergence is not yet achieved. It is remarkable that the convergence of the individual contributions with the basis size is not as slow as might have (X,Y) been anticipated. In fact, EDISP(X,Y) , E(X,Y) C−SP C−WP, and EC−(T) all show smooth convergence in all cases, being converged with quadruple-ζ quality basis sets. Having in mind that the difference of the dispersion contribution between triple-ζ and quintuple-ζ quality basis is less than 0.5 kcal/mol, it is possible to obtain a good estimate of the dispersion interaction already using small basis sets. On the other hand, the preparation and charge transfer terms exhibit a rather poor convergence in both cases. For the methane dimer, a fluctuation in the value of ΔECel−prep, ECT(X→Y) , C−SP and ECT(X←Y) is calculated by changing the cardinal number of C−SP the basis set. This behavior is explained by considering the limitation of the localization procedure applied to the PNOs in relation to the small value that these contributions play in this weakly interacting system. Since the localized PNOs can feature a significant degree of delocalization, their assignment to a given fragment is prone to uncertainties, and hence some fluctuations should be expected in the individual terms, especially for larger basis sets containing diffuse functions. In the strongly interacting NH3−H2O system, such oscillations are observed only for the small ECT(X←Y) term. C−SP

Table 5. Basis Set Dependence of the Correlation Contributions to the Interaction Energya basis CH4-CH4 cc-pVDZ cc-pVTZ cc-pVQZ cc-pV5Z NH3-H2O cc-pVDZ cc-pVTZ cc-pVQZ cc-pV5Z a

ΔECint

C ΔEel−prep

EDISP(X,Y) C−SP

ECT(X→Y) C−SP

ECT(X←Y) C−SP

E(X,Y) C−WP

(X,Y) EC−(T)

−0.42 −0.65 −0.76 −0.80

0.29 0.21 0.47 0.27

−0.45 −0.62 −0.71 −0.72

−0.09 −0.08 −0.21 −0.12

−0.10 −0.06 −0.18 −0.09

−0.01 −0.02 −0.02 −0.02

−0.05 −0.09 −0.11 −0.12

−0.29 −1.37 −1.84 −2.05

6.36 8.49 9.93 11.89

−1.03 −1.42 −1.61 −1.57

−4.67 −6.56 −8.88 −10.76

−0.60 −1.33 −0.64 −0.91

−0.01 −0.04 −0.07 −0.09

−0.32 −0.51 −0.58 −0.61

The energies are given in kcal/mol. 4787

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation

The intermolecular HF contributions show a much stronger dependence on the interaction energy than the intermolecular correlation contributions. Among the latter ones, the charge 3→FCl) transfer contributions ECT(NH changes stronger with the C−SP 3←FCl) and the dispersion distance, while the influence on ECT(NH C−SP contribution is comparably weak. Furthermore, all contributions show a smooth decay behavior with the distance. Regarding the physical interpretation of the results, E(X,Y) elstat and (X,Y) Eexch follow the expected exponential and polynomial decay behavior, shown by the linearity in the double logarithmic and logarithmic scales, respectively (see the insets of Figure 3). For both charge transfer contributions, the decay is approximately exponential, however, interestingly with slightly different decay rates. By plotting the dispersion contribution against the distance on a double logarithmic scale a linear dependency is obtained, with a slope close to six, indicating the expected R−6 decay of dispersion interactions. 3.5. Comparison to DFT-SAPT. A well-established decomposition method especially for weakly interacting systems, which has been applied extensively using correlated methods, is Symmetry Adapted Perturbation Theory (SAPT).5,90,91 If based on DFT, SAPT has the performance to treat the interaction of large systems at sufficient accuracy. Due to this, DFT-SAPT is the preferred method for a comparison of the results of LED for the investigated test cases. It is important to note that SAPT is inherently a molecular approach introducing the interaction of two molecules as a perturbation on the isolated monomers. On the other hand, in the LED-scheme, the interaction energy is calculated explicitly from the energy difference between the dimer and the isolated monomers, and, furthermore, the local energy contributions are obtained by decomposing the energy of the dimer. Thus, the LED is based on orthogonal orbitals throughout, whereas SAPT uses orbitals that are nonorthogonal between the monomers. Hence, contributions as obtained from DFT-SAPT are inherently difficult to compare to results from the LED. For example, as already mentioned, the exchange contribution to the interaction energy is inherently positive in the SAPT approach, while it is inherently negative in the LED approach. Despite these differences, it is still interesting to compare whether both approaches lead to a similar partitioning of dispersion and nondispersion terms in the interaction energy. In order to put the similarities and differences of the two schemes into perspective, a short introduction to DFT-SAPT in given in the following. In the framework of DFT-SAPT the interaction energy is defined as8,83

Table 6. Influence of the Localization Scheme on the Decomposition of the HF Interaction Energy of NH3−H2O and C6H6−C6H6 localization scheme NH3-H2O Foster-Boys Pipek-Mezey C6H6-C6H6 Foster-Boys Pipek-Mezey

ΔEHF int

HF ΔEel−prep

E(X,Y) elstat

E(X,Y) exch

−4.33 −4.33

41.23 38.16

−38.93 −35.63

−6.63 −6.86

3.53 3.53

13.22 19.25

−6.37 −9.64

−3.31 −6.08

Moving to the correlation components of the interaction energy (Table 7), which are much smaller than the ones calculated for the HF term, only minor differences have been found between the two localization schemes. This probably occurs because the HF terms are directly related to the localization of the occupied orbitals, while the correlation contributions also depend on the localization of the virtual orbitals, which are not changed in their localization scheme. Remarkably enough, the LED dispersion is almost independent of the localization scheme, even for the benzene dimer case. A stronger deviation is obtained for the charge transfer components. However, these high relative deviations appear for the CT contributions in the benzene dimer case, which are of small absolute value. On the basis of the discussion above, a small dependence of the decomposition on the applied localization scheme is to be expected, especially at the Hartee Fock level for symmetric systems. We also computed the LED contributions using Intrinsic Atomic Orbitals,89 finding results that are very close to the ones obtained through Pipek-Mizey localization. While a dependence of the results on the localization procedure can be expected, the results show that it is small in most cases and does not change the results of LED qualitatively. 3.4. Distance Dependence of LED Contributions. Any physically sound measure of dispersion interaction is expected to decay as R−6 with the intermolecular distance between the interacting fragments. Furthermore, the electrostatic interaction at the HF level E(X,Y) elstat is expected show a polynomial decay, (X,Y) while the HF exchange interaction Eexch should decay exponentially. In order to test whether the quantities extracted from LED follow this behavior, the decay of these contributions, with respect to the N−Cl distance in NH3− FCl, is presented in this section (see Figure 3). This system has been chosen as a prototype case, because different components, like dispersion and charge transfer, contribute significantly to the interaction energy. For completeness, also the decay of the corresponding ECT(X→Y) and ECT(X←Y) is discussed. Here, a C−SP C−SP simple prediction for the decay behavior is difficult as it provides just a correction for the charge transfer at the HF level.

(1) (2) (2) (2) (2) Eint = E(1) pol + Eexch + Eind + Eexch − ind + Edisp + Eexch − disp + δHF

(27)

Table 7. Influence of the Localization Scheme on the Decomposition of the Correlation Interaction Energy of NH3−H2O and C6H6−C6H6 localization scheme NH3-H2O Foster-Boys Pipek-Mezey C6H6-C6H6 Foster-Boys Pipek-Mezey

ΔECint

C ΔEel−prep

EDISP(X,Y) C−SP

ECT(X→Y) C−SP

ECT(X←Y) C−SP

E(X,Y) C−WP

(X,Y) EC−(T)

−1.84 −1.83

9.93 9.79

−1.61 −1.59

−8.88 −8.73

−0.64 −0.67

−0.07 −0.06

−0.58 −0.57

−5.70 −5.72

2.52 3.82

−4.73 −4.70

−0.98 −1.53

−0.98 −1.54

−0.38 −0.45

−1.16 −1.31

4788

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation

Figure 3. LED energy components decay with the N−Cl distance for NH3−FCl. The dashed gray line indicates the equilibrium distance. The (X,Y) distance dependence of E(X,Y) elstat (blue squares) and Eexch (green triangles) is plotted in the upper panel. The lower part shows the distance dependence (green squares) and the charge transfer contributions of the dissociation energy (blue crosses) beside the one of the dispersion interaction EDISP(X,Y) C−SP 3←FCl) (red circles) and the reverse NH3 to FCl (cyan triangles). Furthermore, two inserts per graph show the following: the decay FCl to NH3 ECT(NH C−SP of the electrostatic and the dispersion contributions on a double logarithmic scale (left) and decay of the Hartree−Fock exchange and of the charge transfer contributions on a logarithmic scale (right).

It includes the attractive first-order electrostatic interaction energy of the unperturbed monomers E(1) pol and the respective (1) repulsive exchange Eexch that arises from the antisymmetry conditions. Furthermore, the second-order induction energy (2) (2) E(2) ind is given, that can be split into Eind (X ← Y) and Eind (X → Y). The first term describes the interaction between the monomers X and Y, with X perturbed by the electrostatic field of Y and the latter the corresponding interaction with Y perturbed by X. E(2) disp is the second-order dispersion interaction energy originating from simultaneous excitations in both monomers. Also for the second-order terms the repulsive (2) exchange terms E(2) exch−ind and Eexch−disp originating from the antisymmetry are included. Interaction terms of higher than second order are approximated in δHF which is calculated by

HF − SAPT δHF = ΔEint (HF ) − Eint (1) (2) = ΔEint (HF ) − [E(1) pol (HF ) + Eexch(HF ) + Eind (HF ) (2) + Eexch − ind (HF )]

(28)

with Eint(HF) as the HF interaction energy of the entire system and EHF−SAPT as the SAPT interaction energy calculated at the int HF-level but omitting the dispersion contributions. It should be noted that the entire system, in the form of a dimer calculation, only enters at this point into DFT-SAPT. Furthermore, for DFT-SAPT, dispersion effects are only included in the intermolecular interaction, while intramolecular correlation is only included as far as the underlying DFT treatment accounts for it. The LED-decomposition scheme includes geometric distortion of the interacting monomers. In most SAPT4789

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation Table 8. Comparison of DLPNO-CCSD(T)/LED to DFT-SAPTa

CH4−CH4 C6H6−C6H6 H2O−H2O NH3−H2O NH3−FCl

ΔEint

EDFT−SAPT int

Edisp

(2) E(2) disp + Eexch−disp

ΔEHF int

(1) (2) (2) E(1) pol + Eexch + Eind + Eexch−ind

δHF

−0.38 −2.17 −4.73 −6.17 −12.16

−0.34 −2.50 −4.54 −6.05 −11.71

−0.73 −5.11 −1.30 −1.67 −6.44

−0.88 −5.57 −2.45 −3.38 −10.16

0.38 3.53 −3.57 −4.33 −7.04

0.57 3.43 −0.99 −0.58 11.75

−0.03 −0.36 −1.11 −2.09 −13.31

Edisp

(2) (2) Edisp + Eexch − disp

ΔEint

DFT − SAPT Eint

1.92 2.35 0.27 0.27 0.53

2.59 2.27 0.54 0.56 0.87

a

The geometric distortion is neglected for a better comparison. All energies are given in kcal/mol. Columns 1, 3, 5, and 8 refer to DLPNOCCSD(T)-LED values, columns 2, 4, 6, 7, and 9 refer to DFT-SAPT. See the text for further details.

interacting systems and weak polar bonds. However, the dispersion interaction as calculated by LED and the perturbatively evaluated dispersion contribution of the DFTSAPT differ by up to a factor of 2 with DFT-SAPT providing systematically larger numbers than DLPNO-CCSD(T)/LED. The difference of both methods in the calculated dispersion contribution at comparable interaction energies leads to a slightly different interpretation of the importance of dispersion contribution. The last two columns of Table 8 show that its relative importance is higher if calculated by DFT-SAPT compared to DLPNO-CCSD(T)/LED. However, the general trends from purely dispersion driven systems to polar bonded adducts is reflected in both theories. While DFT-SAPT can only be applied in the weak interaction limit, DLPNOCCSD(T) can provide accurate results for all interaction strengths while its energy components still remain interpretable. 3.6. Benzene Dimer. In this section, we present a detailed study on the π-interaction of the benzene dimer. For this purpose, the structure of the benzene dimer as described above is used - with the monomer arranged in an offset π−π stacking at a distance of approximately 3.6 Å. For the sake of a better chemical interpretation, Pipek-Mezey-localization has been applied in this section, as it is known to preserve the π/σseparation of the orbitals. The dissociation energy calculated by using the Pipek-Mezey (PM) localization scheme amounts to −2.07 kcal/mol and, hence, just slightly differs from the results obtained in section 3.1 that are obtained by using Foster-Boys localization (FB). The same holds for the outcome of the LED analysis (ΔEHF el−prep (X,Y) = 19.25 kcal/mol; E(X,Y) elstat = −9.64 kcal/mol; Eexch = −6.08 kcal/ mol; ΔECel−prep = 3.82 kcal/mol; EDISP(X,Y) = −4.70 kcal/mol; C−SP ECT(X→Y) = −1.53 kcal/mol; ECT(X←Y) = −1.54 kcal/mol; E(X,Y) C−SP C−SP C−WP + E(X,Y) C−(T) = −1.76 kcal/mol). Exploiting that the PM procedure retains σ- and π-separation, it is interesting to assign symmetry labels to the orbitals of pairs with i and j located on different fragments, as just these pair energies εiXjY contribute to the dispersion energy. With this information it is possible to analyze the dispersion energy EDISP(X,Y) regarding the character of the involved pairs. In Figure C−SP 4, the pair energies of all strong pairs εiXjY are shown, with i and j located on different fragments. For each benzene monomer there are 12 correlated orbitals of σ-symmetry (localized as 6 C−C and 6 C−H bonds) and 3 of π-type. The combinations of these orbitals would lead to 9 π−π pairs, 72 π−σ (and σ−π), and 144 σ−σ pairs. In our case, 9 strong π−π pairs and 72 strong π−σ pairs are selected as strong pairs for the CC level, but only 83 strong σ−σ pairs are obtained due to the pair energy cutoff parameters of the DLPNO-CCSD(T) procedure. While the first class (π−π pairs) shows the largest individual

approaches, this is not part of the DFT-SAPT energy decomposition. Hence, in the following we compare just the ΔEint values of the DLPNO-CCSD(T) interaction energy with results of the DFT-SAPT calculations. Table 8 shows the total interaction energies calculated with DFT-SAPT in comparison to the DLPNO-CCSD(T) values. In accordance with other studies,22 the difference of the two methods in the total interaction energy is small and below 0.5 kcal/mol. Considering the different contributions to the DFTSAPT interaction energy, the changing relevance of δHF for the different types of interaction is striking. Supposed to be a corrective for higher order terms, δHF gains more and more importance with increasing interaction energy. In the case of the weakly interacting and dispersion dominated methane and benzene dimers, it is a rather small contribution to the entire interaction energy (less than 20%). For the stronger interacting hydrogen bonded ammonia−water and water−water adducts, it gains more importance, and, in the case of the strongly interacting ammonia−chlorine fluoride adduct, the δHF correction even exceeds the entire interaction energy. As by δHF the interaction terms of higher order are estimated, it becomes clear that these terms increase for strongly interacting fragments (like for short intermolecular distances). Nevertheless, even in the case of ammonia−chlorine fluoride with δHF being comparable in magnitude to the overall interaction energy, DFT-SAPT is still in good agreement with the coupled cluster interaction energy. Regarding the dispersion contribution to the interactions, DFT-SAPT shows the same order in strength as LED for the investigated test systems: it increases in the order methane dimer, water dimer, ammonia−water adduct, benzene dimer, and ammonia−chlorine fluoride. However, compared to LED, the DFT-SAPT estimate of the dispersion energy is up to a factor of 2 larger. Note that this difference decreases by increasing the intermolecular distance, as shown in the Supporting Information (Table S5) in the water dimer case. Besides the methodological differences, an overestimation of dispersion might be caused by the perturbative treatment of dispersion within DFT-SAPT, as perturbation theory is known for overestimating weak binding energies.92 Furthermore, it is striking that, the larger the δHF contribution to the DFT-SAPT gets, the more LED and SAPT results deviate. For the nondispersive terms, despite δHF, an attraction is obtained for the water dimer and the water−ammonia system, which is dominated by dipole interactions and charge transfer. For the ammonia−chlorine flouride system a strong repulsion is obtained, indicating the strong distortion of the electron densities of the fragments in the electronic structure of the dimer. In summary, DFT-SAPT and DLPNO-CCSD(T) show good agreement for the entire interaction energy of weakly 4790

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation

cules. Accordingly, for the methane and benzene dimers, the dispersion energy is the dominating contribution of the interaction. We have also compared the DLPNO-CCSD(T)/LED scheme to DFT-SAPT as a well-established method for calculating and interpreting weak interactions of large molecules. The total interaction energies predicted by both approaches are in fairly good agreement, while the methods differ by up to a factor of 2 in their estimate of the dispersion contributions. However, qualitatively, both schemes yield comparable results with respect to the interpretation of dispersion dominated versus charge transfer dominated interactions. Besides giving insight in the fundamental contribution to an interaction, LED allows an interpretation of the underlying mechanisms. For example, the interaction between two benzene molecules in a π-stacking conformation was shown to only consist of about 30% genuine π−π interactions with the remaining contributions coming from the smaller but far more numerous σ−π and σ−σ interactions. In summary, DLPNO−CCDS(T)/LED is an efficient, accurate, and nonperturbative approach to calculate and interpret weak intermolecular interactions. It can be applied to virtually any system for which a single-point DFT calculation is still computationally feasible. Actual applications of the method to actual chemical problems will be reported in due course.

Figure 4. Decomposition of the intermolecular strong pair contributions for the benzene dimer case. The orbitals are separated according to their σ and π local symmetry. The type of the excitation is green, ECT(X→Y) , and ECT(X←Y) blue and violet, color-coded: EDISP(X,Y) C−SP C−SP C−SP respectively.

contribution to the dispersion interaction (summing up to 30% of the total interaction energy), the majority of the dispersion interaction (50%) is due to π−σ (and σ−π) pairs, while σ−σ pairs account for the remaining 20%. Hence, this very classical ‘π−π interaction’ case is not really dominated by genuine π−π interactions but rather by the interaction of the π-orbitals with the nearby σ-bonds. This coincides with the previous calculations of Grimme et al.,59 who, based on local MP2 calculations, investigated the intermolecular correlation energy of benzene dimers in various geometrical configurations.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00523. Tables S1−S5 (PDF)

4. CONCLUSIONS The DLPNO-CCSD(T) total energy can be decomposed into physically meaningful contributions that allow for a chemical interpretation of the coupled cluster results. The computational cost of the new procedure, called Local-Energy-Decomposition (LED), is negligible and all that is required is a user defined division of the system into chemically meaningful fragments. A correction for the basis set superposition error is inherent in the method. We have demonstrated that the genuine dispersion energy can be cleanly extracted from the calculations, converges smoothly with respect to basis set extension, and shows the proper decay with respect to intermolecular distance. Other than the dispersion energy contributions, the decomposition contains a dynamic charge polarization mechanism that is counteracting the overpolarization of the charge cloud obtained at the HF level. Furthermore, the interaction energy contains an exchange contribution that is inherently negative and falls off exponentially with distance. The physical picture for the interaction process involves a geometric and an electronic “preparation” of the system to the geometric and electronic structure that the interacting fragments have. In general, an increase of the kinetic energy of the fragments is counteracted by the favorable intermolecular interaction that leads to a decrease of the overall potential energy in just the same way Rüdenberg has described it for the formation of a covalent bond.69 Through application to a set of prototypical interaction types, such as dispersion dominated interactions as well as hydrogenand halogen-bonds, the LED was shown to be in accordance with chemical intuition and is hence a powerful tool in interpreting weak intermolecular interactions between mole-

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS W. B. Schneider and G. Bistoni gratefully acknowledge the financial support of the SPP 1807 “Control of London dispersion interactions in molecular chemistry” of the DFG. The authors thank Dipayan Datta for the implementation of the unrelaxed DLPNO-CCSD density into development version of ORCA.

■

REFERENCES

(1) Riley, K. E.; Hobza, P. Wires Comput. Mol. Sci. 2011, 1, 3. (2) Hohenstein, E. G.; Sherrill, C. D. Wires Comput. Mol. Sci. 2012, 2, 304. (3) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Chem. Rev. 1994, 94, 1887. (4) Szalewicz, K. Wires Comput. Mol. Sci. 2012, 2, 254. (5) Parker, T. M.; Burns, L. A.; Parrish, R. M.; Ryno, A. G.; Sherrill, C. D. J. Chem. Phys. 2014, 140, 094106. (6) Maurer, S. A.; Beer, M.; Lambrecht, D. S.; Ochsenfeld, C. J. Chem. Phys. 2013, 139, 184104. (7) Parrish, R. M.; Hohenstein, E. G.; Sherrill, C. D. J. Chem. Phys. 2013, 139, 174102. (8) Jansen, G. Wires Comput. Mol. Sci. 2014, 4, 127. (9) Perezjorda, J. M.; Becke, A. D. Chem. Phys. Lett. 1995, 233, 134. 4791

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792

Article

Journal of Chemical Theory and Computation (10) Kristyan, S.; Pulay, P. Chem. Phys. Lett. 1994, 229, 175. (11) Antony, J.; Grimme, S. Phys. Chem. Chem. Phys. 2006, 8, 5287. (12) Grimme, S. Wires Comput. Mol. Sci. 2011, 1, 211. (13) Ehrlich, S.; Moellmann, J.; Grimme, S. Acc. Chem. Res. 2013, 46, 916. (14) Grimme, S.; Ehrlich, S.; Goerigk, L. J. Comput. Chem. 2011, 32, 1456. (15) Cramer, C. J.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2009, 11, 10757. (16) Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2005, 1, 415. (17) Pu, J. Z.; Truhlar, D. G. J. Phys. Chem. A 2005, 109, 773. (18) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2006, 110, 5121. (19) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2004, 108, 6908. (20) Zhao, Y.; Truhlar, D. G. Acc. Chem. Res. 2008, 41, 157. (21) Jurecka, P.; Sponer, J.; Cerny, J.; Hobza, P. Phys. Chem. Chem. Phys. 2006, 8, 1985. (22) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. J. Chem. Theory Comput. 2011, 7, 2427. (23) Jaeger, H. M.; Schaefer, H. F.; Hohenstein, E. G.; Sherrill, C. D. Comput. Theor. Chem. 2011, 973, 47. (24) Flick, J. C.; Kosenkov, D.; Hohenstein, E. G.; Sherrill, C. D.; Slipchenko, L. V. J. Chem. Theory Comput. 2012, 8, 2835. (25) Pitonak, M.; Janowski, T.; Neogrady, P.; Pulay, P.; Hobza, P. J. Chem. Theory Comput. 2009, 5, 1761. (26) Janowski, T.; Ford, A. R.; Pulay, P. Mol. Phys. 2010, 108, 249. (27) Sedlak, R.; Riley, K. E.; Rezac, J.; Pitonak, M.; Hobza, P. ChemPhysChem 2013, 14, 698. (28) Pitonak, M.; Riley, K. E.; Neogrady, P.; Hobza, P. ChemPhysChem 2008, 9, 1636. (29) Reinhardt, P.; Piquemal, J. P. Int. J. Quantum Chem. 2009, 109, 3259. (30) Bartlett, R. J.; Musial, M. Rev. Mod. Phys. 2007, 79, 291. (31) Schutz, M.; Yang, J.; Chan, G. K.; Manby, F. R.; Werner, H. J. J. Chem. Phys. 2013, 138, 054109. (32) Schutz, M.; Werner, H. J. J. Chem. Phys. 2001, 114, 661. (33) Schutz, M.; Werner, H. J. Chem. Phys. Lett. 2000, 318, 370. (34) Schutz, M.; Hetzer, G.; Werner, H. J. J. Chem. Phys. 1999, 111, 5691. (35) Hampel, C.; Werner, H.-J. J. Chem. Phys. 1996, 104, 6286. (36) Saebø, S.; Pulay, P. Chem. Phys. Lett. 1985, 113, 13. (37) Saebø, S.; Pulay, P. J. Chem. Phys. 1987, 86, 914. (38) Saebø, S.; Pulay, P. J. Chem. Phys. 1988, 88, 1884. (39) Saebø, S.; Pulay, P. Annu. Rev. Phys. Chem. 1993, 44, 213. (40) Neese, F.; Hansen, A.; Liakos, D. G. J. Chem. Phys. 2009, 131, 064103. (41) Neese, F.; Hansen, A.; Wennmohs, F.; Grimme, S. Acc. Chem. Res. 2009, 42, 641. (42) Neese, F.; Wennmohs, F.; Hansen, A. J. Chem. Phys. 2009, 130, 114108. (43) Hansen, A.; Liakos, D. G.; Neese, F. J. Chem. Phys. 2011, 135, 214102. (44) Liakos, D. G.; Hansen, A.; Neese, F. J. Chem. Theory Comput. 2011, 7, 76. (45) Hansen, A. Ph.D. Thesis, Rheinischen Friedrich-WilhelmsUniversität Bonn, 2012. http://hss.ulb.uni-bonn.de/2012/2976/2976. htm (accessed Sept 1, 2016). (46) Huntington, L. M.; Hansen, A.; Neese, F.; Nooijen, M. J. Chem. Phys. 2012, 136, 064101. (47) Riplinger, C.; Neese, F. J. Chem. Phys. 2013, 138, 034106. (48) Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. J. Chem. Phys. 2013, 139, 134101. (49) Ahlrichs, R.; Kutzelnigg, W. Theor. Chim. Acta. 1968, 10, 377. (50) Meyer, W. Theor. Chim. Acta. 1974, 35, 277. (51) Ahlrichs, R.; Driessler, F.; Lischka, H.; Staemmler, V.; Kutzelnigg, W. J. Chem. Phys. 1975, 62, 1235. (52) Kutzelnigg, W.; Staemmler, V.; Hoheisel, C. Chem. Phys. 1973, 1, 27. (53) Staemmler, V. Theor. Chim. Acta 1973, 31, 49.

(54) Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. J. Chem. Phys. 2016, 144, 024109. (55) Liakos, D. G.; Neese, F. J. Chem. Theory Comput. 2015, 11, 4054. (56) Morokuma, K. Acc. Chem. Res. 1977, 10, 294. (57) Glendening, E. D.; Streitwieser, A. J. Chem. Phys. 1994, 100, 2900. (58) Glendening, E. D.; Landis, C. R.; Weinhold, F. Wires Comput. Mol. Sci. 2012, 2, 1. (59) Grimme, S.; Muck-Lichtenfeld, C.; Antony, J. Phys. Chem. Chem. Phys. 2008, 10, 3327. (60) Schutz, M.; Rauhut, G.; Werner, H. J. J. Phys. Chem. A 1998, 102, 5997. (61) Azar, R. J.; Head-Gordon, M. J. Chem. Phys. 2012, 136, 024103. (62) Thirman, J.; Head-Gordon, M. J. Chem. Phys. 2015, 143, 084124. (63) Boys, S. F. Rev. Mod. Phys. 1960, 32, 296. (64) Pipek, J.; Mezey, P. G. J. Chem. Phys. 1989, 90, 4916. (65) Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. J. Chem. Phys. 2016, 144, 024109. (66) Liakos, D. G.; Sparta, M.; Kesharwani, M. K.; Martin, J. M. L.; Neese, F. J. Chem. Theory Comput. 2015, 11, 1525. (67) Boys, S. F.; Bernardi, F. d. Mol. Phys. 1970, 19, 553. (68) Mentel, Ł. M.; Baerends, E. J. J. Chem. Theory Comput. 2014, 10, 252. (69) Ruedenberg, K. Rev. Mod. Phys. 1962, 34, 326. (70) Heßelmann, A.; Jansen, G.; Schütz, M. J. Chem. Phys. 2005, 122, 244108. (71) Azar, R. J.; Head-Gordon, M. J. Chem. Phys. 2012, 136, 024103. (72) Ziegler, T.; Rauk, A. Theor. Chim. Acta. 1977, 46, 1. (73) Mitoraj, M. P.; Michalak, A.; Ziegler, T. J. Chem. Theory Comput. 2009, 5, 962. (74) Mo, Y.; Gao, J.; Peyerimhoff, S. D. J. Chem. Phys. 2000, 112, 5530. (75) Bistoni, G.; Belpassi, L.; Tarantelli, F. J. Chem. Theory Comput. 2016, 12, 1236. (76) Lao, K. U.; Herbert, J. M. J. Chem. Theory Comput. 2016, 12, 2569. (77) Khaliullin, R. Z.; Cobar, E. A.; Lochan, R. C.; Bell, A. T.; HeadGordon, M. J. Phys. Chem. A 2007, 111, 8753. (78) Neese, F. Wires Comput. Mol. Sci. 2012, 2, 73. (79) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (80) Eichkorn, K.; Treutler, O.; Ohm, H.; Haser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240, 283. (81) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Wiley Interdisciplinary Reviews: Computational Molecular Science 2012, 2, 242. (82) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (83) Heßelmann, A.; Jansen, G. Chem. Phys. Lett. 2003, 367, 778. (84) Lia, S. G. Ionization Energy data from the NIST Standard Reference Database Number 69. http://webbook.nist.gov/chemistry (accessed Sept 1, 2016). (85) Wang, C.; Guan, L.; Danovich, D.; Shaik, S.; Mo, Y. J. Comput. Chem. 2015. (86) Cramer, C. J. Essentials of computational chemistry: theories and models; John Wiley & Sons: 2004. (87) In The unrelaxed DLPNO-CCSD density is implemented in a development version (based on version 3.0.3) of the ORCA suite of programs. (88) Van Duijneveldt-van De Rijdt, J. G. C. M.; Van Duijneveldt, F. B. J. Chem. Phys. 1992, 97, 5019. (89) Knizia, G. J. Chem. Theory Comput. 2013, 9, 4834. (90) Parker, T. M.; Sherrill, C. D. J. Chem. Theory Comput. 2015, 11, 4197. (91) Parrish, R. M.; Sherrill, C. D. J. Chem. Phys. 2014, 141, 044115. (92) Shirkov, L.; Makarewicz, J. J. Chem. Phys. 2015, 142, 064102.

4792

DOI: 10.1021/acs.jctc.6b00523 J. Chem. Theory Comput. 2016, 12, 4778−4792