Subscriber access provided by University of Sussex Library
Article
Intricacies of van der Waals interactions in systems with elongated bonds revealed by electron-groups embedding and high-level coupled-cluster approaches. Ewa Pastorczak, Jun Shen, Michal Hapka, Piotr Piecuch, and Katarzyna Pernal J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00797 • Publication Date (Web): 18 Sep 2017 Downloaded from http://pubs.acs.org on September 19, 2017
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 free 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 accessible to all readers and 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.
Journal of Chemical Theory and Computation 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 63
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
Intricacies of van der Waals interactions in systems with elongated bonds revealed by electron-groups embedding and high-level coupled-cluster approaches Ewa Pastorczak,† Jun Shen,‡ Michal Hapka,¶ Piotr Piecuch,∗,‡,§ and Katarzyna Pernal∗,† Institute of Physics, Lodz University of Technology, ul. Wolczanska 219, 90-924 Lodz, Poland, Department of Chemistry, Michigan State University, East Lansing, Michigan 48824, United States, Faculty of Chemistry, University of Warsaw, ul. Pasteura 1, 02-093 Warsaw, Poland, and Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, United States E-mail:
[email protected];
[email protected] Abstract Non-covalent interactions between molecules with stretched intramonomer covalent bonds are a fascinating, yet little studied area. This shortage of information stems largely from the inability of most of the commonly used computational quantum ∗
To whom correspondence should be addressed Institute of Physics, Lodz University of Technology, ul. Wolczanska 219, 90-924 Lodz, Poland ‡ Department of Chemistry, Michigan State University, East Lansing, Michigan 48824, United States ¶ Faculty of Chemistry, University of Warsaw, ul. Pasteura 1, 02-093 Warsaw, Poland § Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, United States †
1
ACS Paragon Plus 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
chemistry methods to accurately describe weak long-range and strong non-dynamic correlations at the same time. In this work, we propose a geminal-based approach, abbreviated as EERPA-GVB, capable of describing such systems in a robust manner using relatively inexpensive computational steps. By examining a few van der Waals complexes, we demonstrate that the elongation of one or more intramolecular covalent bonds leads to an enhanced attraction between the monomers. We show that this rise in attraction occurs as the electron density characterizing intramolecular covalent bonds depletes and migrates toward the region between the monomers. As the covalent intramonomer bonds continue stretching, the intermolecular interaction potential passes through a minimum and eventually goes up. The findings resulting from our EERPA-GVB calculations are supported and further elucidated by the symmetryadapted perturbation theory and coupled-cluster (CC) computations using methods that are as sophisticated as the CC approach with a full treatment of singly, doubly, and triply excited clusters.
1
Introduction
Non-covalent or van der Waals interactions occur ubiquitously in nature. Although much weaker than the covalent bonds, they often have a significant impact on molecular geometries, barrier heights, phase transitions, and various properties and phenomena in condensed phases. The role of van der Waals forces in facilitating catalytic reactions, 1 rational design of materials, 2 molecular recognition, 3 and photoinduced chemical dynamics, 4 to name a few examples, is being appreciated more and more as the increasingly sophisticated experimental and computational techniques become available. Still, in spite of numerous advances in quantum chemistry, certain classes of non-covalent interactions remain challenging for the majority of the existing computational methods, especially when larger systems, which prevent us from using the highest ab initio levels of electronic structure theory, are to be examined. For example, both the widely used density functional theory (DFT) 5 and multi-configuration
2
ACS Paragon Plus Environment
Page 2 of 63
Page 3 of 63
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
self-consistent field (MC-SCF) 6 approaches experience difficulties with capturing the London dispersion forces, mainly because of their dynamic long-range non-local character. 7 The DFT methods can alleviate this problem through the introduction of parametrized, empirical, semi-classical dispersion corrections and the construction of non-local forms of the exchange-correlation functionals. 8 The wavefunction, MC-SCF-type, approaches can capture the long-range dispersion contributions through perturbative corrections, 6,9 but the resulting intermolecular potentials are often less accurate than desired, so much of the challenge remains. Having examined the role of van der Waals interactions between closed-shell monomers near their equilibrium geometries, one would naturally want to turn to weakly bound complexes attempting reactions through the activation or stretching of the intramonomer covalent bonds, but this is beyond reach of commonly used methods. Indeed, the widely used DFT methodology cannot, at least in its present form, correctly describe molecules far from equilibrium, especially those undergoing covalent bond breaking. This deficiency obviously hinders DFT’s ability to handle reactive potential energy surfaces, including those involving weakly bound van der Waals precursors and significant rearrangements of covalent bonds, where one has to consider long-range intermolecular and strong non-dynamic correlations at the same time. Popular ab initio wavefunction approaches face difficulties too. The widely used multireference methods based on MC-SCF, such as the second-order complete active space SCF (CASPT2) approach, 6,9 break chemical bonds, but being low-order approaches in the treatment of dynamic correlations they are often insufficient in providing an accurate description of weak intermolecular forces characterizing van der Waals precursor complexes or local minima and transition states having substantial contributions from non-covalent interactions. Furthermore, multireference methods require a great deal of expertise and become prohibitively expensive when the numbers of active orbitals and active electrons used in the calculations become larger. 6,10 Weakly bound intermolecular complexes can accurately be described by the single-reference methods based on coupled-cluster (CC) theory, 11–13
3
ACS Paragon Plus 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
including, for example, the widely used CCSD(T) approximation, 14 in which one adds a quasi-perturbative non-iterative correction due to connected triples to the energy obtained with the CC singles and doubles (CCSD) approach, 15 but they are not robust enough to handle covalent bond breaking unless one incorporates higher-than-doubly excited clusters fully, as in the CC approach with singles, doubles, and triples (CCSDT) 16,17 and its higher-order extensions, which are usually prohibitively expensive. 18 There has been significant progress in recent years in formulating alternatives to CCSD(T), which can handle both bond breaking and non-covalent interactions with computational costs that are considerably lower than those of CCSDT, including the left eigenstate completely renormalized CC approaches, 19–21 such as the CR-CC(2,3) triples correction to CCSD, 19,20,58 the active-space CC theories, 22 such as the CC method with singles, doubles, and the active-space treatment of triples, abbreviated as CCSDt, 59–62 and their merger through the so-called CC(P ;Q) formalism 21,23 represented in this work by the CC(t;3) correction to CCSDt. 21,23,39 Yet, in spite of these advances not much theoretical work has been done in the field of van der Waals interactions involving stretched intramonomer covalent bonds. The impact of noncovalent interactions on chemical reaction pathways and dynamics, including van der Waals wells characterizing precursor complexes and intermediates and saddle points having significant contributions from weak intermolecular interactions, has been highlighted for a number of elementary reactions. 24–27 Unfortunately, in some cases, the computations pointing out the importance of van der Waals interactions fail to include dispersion (see, e.g., ref 28). It seems fair to say that the struggle between the need for an accurate description of covalent bond dissociations and London dispersion forces, while keeping the computational costs at the relatively low level, continues. Related problems, where the covalent intramolecular bonds are stretched, but not broken, are simulations of ro-vibrational spectra of van der Waals complexes involving vibrationally excited monomers. It is well-known 29 that it is hard to obtain correct spectra of this type when the conventional perturbative approaches are employed, and yet the aforementioned
4
ACS Paragon Plus Environment
Page 4 of 63
Page 5 of 63
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
CCSD(T) approach is often suggested as the method of choice, in spite of the fact that, as already alluded to above, CCSD(T) struggles when covalent bonds are stretched. We show in this work that the use of CCSD(T) in calculations involving van der Waals complexes with stretched intramonomer bonds can indeed lead to severe problems. Apart from the molecules with stretched or broken bonds, another area where multireference effects and van der Waals interactions meet are strongly correlated systems, including compounds containing transition metals. 30 Such systems are ubiquitous in the aforementioned catalytic processes and adsorption phenomena, 31 but their computational description is often limited to the use of DFT, since the more suitable ab initio wavefunction methods are usually too expensive. This may lead to an erratic description (see, e.g., ref 32). The above problems can be addressed by exploiting expensive high-level ab initio methods, such as the CCSDT approach mentioned above, but the question arises if one can develop a more affordable methodology that works equally well for van der Waals interactions and intramonomer bond stretching, while being applicable to larger systems. We address this question in the present work by proposing a novel computational approach that combines the extended random-phase approximation (ERPA) 33 and generalized valence bond (GVB) 34,35 formalisms in the framework of electron embedding. 36 The resulting theory, Embedding Extended Random Phase Approximation GVB, abbreviated as EERPAGVB, improves upon the recently proposed ERPA-GVB method 36–38 and offers an excellent description of intermolecular interactions involving monomers with stretched covalent bonds, rivaling in accuracy some of the most sophisticated CC methods, such as CCSDT, at a small fraction of their computational costs. We employ the EERPA-GVB method to a few weakly bound systems involving stretched intramonomer bonds, ranging from the smaller H2 · · · H2 and NH3 · · · F2 species, in which the H–H bond in one of the H2 molecules and the F–F bond in F2 undergo dissociation, to the strongly correlated H12 · · · H12 dimer, in which all H–H bonds in the linear H12 chains are simultaneously stretched, undergoing a metal–insulator transition. The perfor-
5
ACS Paragon Plus 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
mance of EERPA-GVB is compared against the high-level CC methods, with a focus on the full CCSDT approach and the CC(t;3) approximation to it based on the aforementioned CC(P ;Q) framework. 21,23,39 We also use the symmetry-adapted perturbation theory (SAPT) 40 to provide insights into the behavior of the studied complexes upon intramonomer bond stretching.
2
EERPA-GVB: A novel theoretical approach based on electron-groups embedding
The challenge in description of intermolecular interactions in systems comprising strongly correlatefd monomers lies in the necessity of capturing substantial non-dynamic and longrange dynamic correlation effects with sufficient accuracy and on equal footing. Typically, the former correlations result from stretching covalent bonds in the interacting monomers, whereas the latter correlation effects are responsible for emergence of van der Waals interactions between molecules. The recently introduced ERPA approach based on the APSG reference wavefunction, abbreviated as ERPA-APSG, 36,37 where APSG stands for the antisymmetrized product of strongly orthogonal geminals, 34,35,41–43 and its simplified ERPAGVB variant 38 account for the strong correlation effects and are, therefore, suitable for describing interactions in multireference systems, but they underestimate the interaction energies. In this study, we demonstrate that the accuracy of ERPA-GVB calculations can be substantially improved using the EERPA-GVB approach, in which we correct the correlation energies within the interacting monomers by embedding electrons of one monomer in the field of the other monomer and correlating the embedded electron groups. The purpose of this section is to present a general EERPA formalism based on the generalized product function (GPF) ansatz for the wavefunction along with a correlation energy correction derived from ERPA. Already in the 1950s McWeeny formulated a GPF ansatz that provides a suitable zeroth6
ACS Paragon Plus Environment
Page 6 of 63
Page 7 of 63
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
order wavefunction for weakly interacting subsystems. 44 GPF for N electrons partitioned into M groups consisting of NI , I = 1, . . . , M , electrons is defined as
Ψ
GP F
= Aˆ
M Y
ΨI ,
(1)
I=1
where each ΨI is an antisymmetric NI -electron wavefunction. The antisymmetry of the ˆ which includes a proper total wavefunction is assured by the antisymmetrizing operator A, normalization factor. To reduce the complexity of the optimization problem, McWeeny imposed a strong orthogonality constraint on each pair of the group wavefunctions, i.e.,
∀I6=J
ˆ
ΨI (x1 , . . . , xNI )∗ ΨJ (x1 , . . . , xNJ ) dx1 = 0
(2)
(here and elsewhere in the present paper x combines spatial and spin coordinates of a single electron). Strong orthogonality leads to a simple expression for the two-electron reduced density matrix (2-RDM) of the total N -electron system, namely,
Γ
GP F
(x1 , x2 ; x′1 , x′2 )
= +
M X
I=1 M X I6=J
ΓI (x1 , x2 ; x′1 , x′2 ) (3)
γ I (x1 , x′1 )γ J (x2 , x′2 ) − γ I (x1 , x′2 )γ J (x2 , x′1 ) ,
where ΓI and γ I stand for, respectively, two- and one-electron reduced density matrices pertaining to ΨI . The total one-electron reduced density matrix (1-RDM) is a sum of group 1-RDMs, γ GP F (x1 , x′1 ) =
M X
γ I (x1 , x′1 ),
(4)
I=1
with a similar relation holding for a total electron density function, ρGP F =
PM
I=1
ρI . Writing
the 2-RDM of ΨGP F in the representation of the natural spinorbitals {ϕp (x)} of γ GP F and
7
ACS Paragon Plus 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 63
exploiting the fact that the strong orthogonality condition given by eq 2 implies 45
∀I6=J ∀p
ˆ
I
′
∗
′
′
Ψ (x , . . . , xNI ) ϕp (x )dx 6= 0
ˆ
⇒
ΨJ (x′ , . . . , xNJ )∗ ϕp (x′ )dx′ = 0
(5)
leads to the following expression for the matrix elements of 2-RDM in the representation of the natural spinorbitals:
F ΓGP pqrs
ΓIpqrs
GP F † † GP F = = Ψ |ˆ ar a ˆs a ˆq a ˆp |Ψ n n (1 − δ p q
where
and
if Ip = Iq = Ir = Is = I Ip Iq )(δpr δqs
(6)
− δps δqr ) otherwise
ΓIpqrs = ΨI |ˆ a†r a ˆ†s a ˆq a ˆp |ΨI
(7)
I ˆp |ΨI = np δpq , γpq = ΨI |ˆ a†q a
(8)
with a ˆ†p (ˆ ap ) representing the usual creation (annihilation) operators and np designating the occupation number associated with spinorbital ϕp (each group density matrix γ I is diagonal in the natural spinorbitals of the total γ GP F ). The relation given in eq 5 means that the natural spinorbitals are partitioned into disjoint subsets in which group wavefunctions are expanded. Symbol Ip in eq 6 indicates a subspace to which spinorbital ϕp belongs. The structure of the two-electron density matrix, eq 3 or 6, implies that the interaction of electrons assigned to different groups I are treated only at the Coulomb and exchange level. Indeed, the electronic energy expression obtained for the GPF ansatz reads E GP F =
X
EI +
I
X
IJ ECoul−Exch ,
(9)
I>J
where each group energy term, E I , accounts for the energy of the electrons assigned to group
8
ACS Paragon Plus Environment
Page 9 of 63
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
I, namely, EI =
X
δIp I hpp +
p
*
N + I X 1 I I Ψ Ψ , rij
(10)
i