17152
J. Phys. Chem. 1996, 100, 17152-17156
Natural Energy Decomposition Analysis: The Linear Response Electrical Self Energy Gregory K. Schenter*,† and Eric D. Glendening*,‡ EnVironmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory, Richland, Washington 99352, and Department of Chemistry, Indiana State UniVersity, Terre Haute, Indiana 47809 ReceiVed: May 7, 1996; In Final Form: July 17, 1996X
An extension of the natural energy decomposition analysis (NEDA) is described that leads to a reduced, three-component treatment of ab initio molecular interaction potentials. These components include the classical electrical (static and induced) and quantum mechanical core (σ-σ) and charge transfer (σ-σ*) interactions. The electrical component is, in turn, represented by three terms: the static interaction of unperturbed fragment charge densities, a contribution involving polarized charge densities, and an electrical self polarization energy. Extended basis set calculations demonstrate the high numerical stability of the three-component NEDA approach. Applications are presented for the Li+(H2O) and water dimer complexes. The long-range behavior of the interaction potentials of these complexes is entirely determined by the classical electrical interaction of the fragments. Core and charge transfer effects are only significant for molecular separations within roughly 1 Å of equilibrium.
I. Introduction Natural energy decomposition analysis (NEDA)1,2 was recently introduced as an efficient and numerically stable computational tool for analyzing molecular interactions potentials. Like the popular Morokuma analysis,3-5 NEDA evaluates contributions such as electrostatic interaction, polarization, and charge transfer, energy components that in part form the basis for our conceptual models of molecular interactions. Previous applications have demonstrated the high numerical stability of the method and reasonable behavior of its energy components.6-11 In its present form,2 NEDA evaluates electrostatic (ES), polarization (POL), charge transfer (CT), exchange (EX), and deformation (DEF) components. It is advantageous, however, to consider the reduction of these five components to three, one that arises solely from classical electrical interactions (both static and induced) and two others associated with quantum mechanical core (σ-σ) and charge transfer (σ-σ*) interactions. ES and POL clearly contribute to the electrical interaction, CT to the charge transfer interaction, and EX (the intermolecular exchange of like-spin electrons) to the core interaction. The repulsive DEF component is, on the other hand, somewhat problematic, since the source of deformation depends on the region of the potential under consideration. At short range, DEF arises from Pauli repulsions that prevent significant interpenetration of molecular charge densities. At long range, DEF is the electrical energy penalty (self energy) to polarize the interacting molecules. DEF thus appears to include both long-range electrical and short-range core contributions. In this work, we apply linear response theory to NEDA to separate the electrical self energy (SE) from the DEF component. In doing so, we reduce the molecular interaction potential to three components: electrical interaction (EL), charge transfer (CT), and core repulsions (CORE). Section II briefly describes the NEDA method and discusses the linear response treatment of the SE contribution. Section III then presents two applications to the Li+(H2O) and water dimer complexes. Extended † Pacific Northwest National Laboratory. Email address: gk_schenter@ pnl.gov. ‡ Indiana State University. Email address:
[email protected]. X Abstract published in AdVance ACS Abstracts, September 1, 1996.
S0022-3654(96)01299-3 CCC: $12.00
basis set calculations of the Li+(H2O) complex demonstrate the high numerical stability of the three-term analysis. We show that the long-range regions of the Li+(H2O) and (H2O)2 interaction potentials are entirely described by EL and that EL behaves in a manner consistent with traditional empirical treatments of electrical interaction that employ instantaneous polarization response. CT and CORE are only significant within about 1 Å of equilibrium separation, i.e., at geometries where there exists significant overlap of the molecular charge densities. II. Natural Energy Decomposition Analysis In its original form,1 NEDA partitioned the total interaction energy, ∆E, into ES, CT, and DEF contributions. More recently,2 the method was extended by separating the POL and EX contributions from ES so that ∆E is treated as a sum of five terms,
∆E ) ES + POL + CT + EX + DEF
(1)
We now identify the electrical SE and define the total electrical interaction EL,
EL ) ES + POL + SE
(2)
and the CORE contribution,
CORE ) EX + DEF - SE
(3)
such that the interaction energy is given by the three-term partition
∆E ) EL + CT + CORE
(4)
In doing so, we separate the charge density interactions of EL that are well described by classical electrodynamics from the CT and CORE contributions that arise from the quantum mechanical nature of the density. NEDA is based on Weinhold’s natural bond orbital (NBO) approach.12-14 NBO calculates a set of mutually orthogonal bond orbitals that, aside from small orthogonalization tails, are strictly localized either on single atoms (core orbitals, lone pairs) or on atom pairs (bonds, antibonds). The monomer units of the complex are naturally defined by the connectivity of the © 1996 American Chemical Society
Natural Energy Decomposition Analysis
J. Phys. Chem., Vol. 100, No. 43, 1996 17153
calculated bonds. Thus, NBO analysis of water dimer reveals two water monomers, each consisting of an oxygen 1s core, two lone pairs, and two OH bonds. The NEDA partitioning employs the SCF-converged Fock matrix for the full complex (supermolecule) expressed in the NBO basis. NEDA, at least in its present form, is thereby restricted to Fock matrix-based electronic structure methods. To simplify the presentation, we limit our discussion here to restricted Hartree-Fock (RHF) wave functions. A key element of the NEDA approach is the evaluation of the intermediate monomer wave functions, ΨAdef, one for each of the monomers that comprise the complex. The wave function for the intermediate, perturbed monomer A is given by
of Ftot(r). In terms of molecular orbitals, we write the following explicit expression for the charge density associated with the perturbed monomer A on A
FA(r) )
∑R
on A
ZRδ(r - RR) - 2 ∑ |φa(r)|2
where ZR is the nuclear charge at coordinate RR and the summations are over all nuclei and orbitals that comprise monomer A only. Similarly, we define FA0(r) as the unperturbed charge density of the isolated monomer that is associated with the wave function ΨA, on A
on A
FA0(r) )
on A
ΨAdef ) A[∏ φa]
(13)
a
(5)
|φa′(r)|2 ∑R ZRδ(r - RR) - 2 ∑ a′
(14)
a
where A is the antisymmetrizer and the molecular orbitals are the local block eigenvectors of the NBO Fock matrix of the full complex.1,2 This wave function is perturbed or “deformed” in the sense that the orbitals describe the redistribution of the charge density that arises from electric field and quantum mechanical effects experienced by the monomer within the complex. A localized wave function for the complex is then constructed as an antisymmetrized product of perturbed monomer wave functions
Ψloc ) A[∏ΨAdef]
In terms of charge densities, ES can be written as
ES )
1
FA0(r)FB0(r′)
∑ ∫ dr dr′
2 A,B
(15)
A*B
while the interaction of the perturbed monomer charge densities gives the ES + POL contribution,
ES + POL )
(6)
1
FA(r)FB(r′)
∑ ∫ dr dr′
2 A,B
|r - r′|
(16)
A*B
A
NEDA further employs the supermolecule wave function, Ψ, and the SCF-converged wave functions for the isolated monomers, ΨA. The latter are given by
At this point, to be consistent with a linear response description of the self polarization energy,15 we define the induced monomer charge density as
∆FA(r) ) FA(r) - FA0(r)
on A
ΨA ) A[∏ φa′]
(7)
a′
where the prime indicates that the molecular orbitals are fully optimized in the absence of all other monomers. Given these wave functions, the interaction energy and its NEDA components are obtained as follows:
∆E ) E(Ψ) - ∑E(ΨA)
(8)
CT ) E(Ψ) - E(Ψloc)
(9)
A
ES + POL + EX ) E(Ψ ) - ∑ loc
E(ΨAdef)
DEF ) ∑[E(ΨAdef) - E(ΨA)]
(11)
A
Further details of the evaluation of these components are given elsewhere.1,2 It is convenient to consider the representation of the total charge density of the interacting system, Ftot(r). The orthogonal orbitals that define the localized wave function provide a welldefined partitioning of the charge density of the system that is not associated with CT. We denote the charge density associated with Ψloc by F(r). From its construction, F(r) can be partitioned in terms of the charge densities of the perturbed monomers (ΨAdef) that we denote FA(r). Explicitly,
(12)
A
In cases where CT is negligible, F(r) is a good representation
(17)
and require that self energy be a quadratic functional of the form
SE )
1 2
∫dr dr′ ∆FA(r)KA(r - r′)∆FA(r′) ∑ A
(18)
where KA(r - r′) is the linear response kernel. This energy is the cost to form the induced charge density. Furthermore, we require that ∆FA(r) be determined from the condition
δ EL ) 0 δ∆FA(r)
(10)
A
F(r) ) ∑FA(r)
|r - r′|
(19)
representing an instantaneous linear response of the induced charge density to other charges. The following expression results,
δ
EL ) ∫dr′ KA(r - r′)∆FA(r′) +
δ∆FA(r)
FB(r′)
∑B ∫dr′ |r - r′| B*A
(20)
to give
SE ) ∑SEA
(21)
A
where SEA is the self energy of monomer A,
SEA ) -
1 2
∑ ∫ dr dr′
B B*A
∆FA(r)FB(r′) |r - r′|
(22)
17154 J. Phys. Chem., Vol. 100, No. 43, 1996
Schenter and Glendening
By collection of terms, EL may be written as
EL )
1
∑ ∫ dr dr′
FA0(r)FB(r′) (23)
|r - r′|
2 A,B
A*B
This treatment of the total electrical energy is consistent with the traditional approach used in molecular dynamics simulations that employ instantaneous polarization response. For the case of two monomers, A and B, the contributions to the electrical energy may be written explicitly using standard notation16 as
ZRZβ
on A on B
ES )
∑R ∑β
RRβ
on A on B
- 2∑
∑ b′
R
on A on B
2∑
∑β
a′
〈 |R |b′〉 ZR
b′
〈| |〉 Zβ
a′
Rβ
a′ + 4 ∑
R
∑ b′
〈| |〉 ZR
b′
RR
R′
∑β
a′
(a′a′|b′b′) - 2
a′
Zβ
on A on B
a
2
β
a
Zβ
on A on B
a
Rβ
β
(24)
∑R ∑b
RR
(aa|bb) (25)
a
b
(aa|bb) -
a
Zβ
β
Rβ
b
on A on B
a′ + 2
(a′a′|bb) (26)
a′
b
on A on B
b - 2∑ a
ZR
on A on B
b′
R
b -
on A on B
〈| |〉 ∑ ∑∑〈 | | 〉 ZR
b
RR
Figure 1. Basis set convergence of the interaction energy and NEDA components for Li+(H2O). The squares are values calculated with the cc-pVXZ family of basis sets (X ) D, T, Z). The circles are the corresponding values for the aug-cc-pVXZ basis sets. Exponential fits of the aug-cc-pVXZ data points are shown.
a′ -
on A on B
a′
a′
Rβ
b
a +4
a -2
on A on B
on A on B
b
R
b′
Rβ
ZR
on A on B
4
a
Zβ
on A on B
b′ + 2 ∑
on A on B
SEB )
(a′a′|b′b′) ∑ b′
a′
〈| |〉 ∑∑ ∑∑〈 | | 〉 ∑∑〈 | | 〉 ∑∑ ∑∑〈 | | 〉 ∑∑ ∑∑〈 | | 〉 ∑∑ on A on B
POL ) 2 ∑
SEA )
R on A on B
b′
RR
(aa|bb) -
b
on A on B
b′ + 2 ∑ a
(aa|b′b′) ∑ b′
(27)
Although we only focus in the present work on the HartreeFock implementation of NEDA, note that eqs 15, 16, and 22 represent general expressions that are equally applicable to correlated electronic structure methods. In fact, it is possible, in principle, to extend NEDA to post-Hartree-Fock methods provided one can evaluate the correlated monomer charge densities of eq 17 and localized wave function of eq 6. The three-term NEDA approach has been implemented in the GAMESS17 version of the NBO program.18 The method can be applied to RHF and UHF wave functions using either the conventional disk-based or direct two-electron integral techniques. All wave functions (including those of the monomers) are evaluated in the full atomic orbital basis set of the complex. In addition to being convenient (the two-electron integrals need only be evaluated once in a conventional calculation), the full basis set treatment corrects the calculated interaction energy and NEDA components for unphysical basis set superposition errors (BSSE) according to the counterpoise method of Boys and Bernardi.19 All calculations reported in the following section were performed with GAMESS. III. Applications A. Li+(H2O). As a first application, we consider the strong charge-dipole interaction of the Li+(H2O) complex. This
closed-shell system is well described at the RHF level of theory9,20 and was previously employed to demonstrate the numerical stability of the five-term NEDA approach.2 The high numerical stability of the three-term NEDA approach is likewise revealed by an extended basis set treatment of the Li+(H2O) complex. Calculations were performed on a C2V symmetry complex [R(Li-O) ) 1.85 Å, the Li-O collinear with the bisector of the HOH bond angle] with a constrained water geometry [R(OH) ) 0.9572 Å, ∠(HOH) ) 104.52°].21 Figure 1 shows the interaction energy and NEDA components evaluated with Dunning’s cc-pVXZ and aug-cc-pVXZ basis sets, where X ) D, T, and Q for double-, triple-, and quadruple-ζ quality, respectively).22-24 The augmented correlation consistent sets (aug-cc-pVXZ) differ from the standard sets (cc-pVXZ) by the addition of diffuse functions that are often required to treat quantitatively molecular interactions. The interaction energy and NEDA components of Li+(H2O) converge rapidly with increasing basis set quality. Previous studies with the correlation consistent sets have shown approximate exponential convergence of various properties, and similar convergence behavior is seen in Figure 1. By fitting the aug-cc-pVXZ data with a three-parameter function of the form
P(x) ) (P0 - PCBS)e-kx + PCBS
(28)
where x ) 2, 3, 4 for X ) D, T, Q, we estimate complete basis set limits (PCBS) of -35.8 (∆E), -47.0 (EL), -6.2 (CT), and 17.4 (CORE) kcal mol-1. Since the aug-cc-pVTZ results differ by less than 0.3 kcal mol-1 from these estimated limits, we use this basis set in the remainder of the calculations reported here. Figure 2 shows the three-term NEDA partition of the Li+(H2O) interaction potential. The separation of long-range electrical interaction from short-range quantum mechanical effects is evident. EL clearly dominates the long-range region of the potential through the charge-dipole interaction of Li+ with H2O. Indeed, the EL and ∆E curves differ by less than 0.04 kcal mol-1 for all Li-O distances greater than 2.75 Å (equilibrium separation is ∼1.85 Å). Overlap of the Li+ and H2O charge densities is weak in this region of the potential such that the CORE and CT contributions remain negligible. At shorter distances, where CT and CORE become significant, ∆E deviates strongly from EL. CT arises in the Li+(H2O) complex
Natural Energy Decomposition Analysis
Figure 2. Three-term NEDA (aug-cc-pVTZ) of the interaction potential for Li+(H2O). The dashed curve represents the classical electrical interaction assuming a point charge at Li and a point dipole and polarizability at the center-of-mass of H2O (cf. text).
Figure 3. Five-term NEDA (aug-cc-pVTZ) of the interaction potential for Li+(H2O).
from transfer of electron density out of an axial lone pair on O into the valence 2s orbital of Li. However, the effect remains relatively weak, only stabilizing the complex by 5.8 kcal mol-1 (16% of ∆E) at equilibrium. Figure 2 also compares the EL component of the Li+(H2O) potential with the classical electrical approximation, EL(cl). The latter is evaluated assuming instantaneous polarization response of a polarizable dipole interacting with a point charge. The water charge density is represented by a point dipole (µz ) 1.975 D) and polarizability (Rzz ) 8.427 a03) placed at its center of mass.25 The EL and EL(cl) curves differ by less than 0.2 kcal mol-1 for all Li-O distances greater than 3 Å. The curves diverge somewhat at shorter Li-O distances where the Li+ and H2O charge densities overlap significantly and the point chargedipole approximation fails. Comparison of the three- and five-term NEDA treatment of the Li+(H2O) potential reveals the advantage of the former. The five-term analysis is shown in Figure 3, where the ES + POL contribution is presented by a single curve. The ES + POL contribution is significantly stronger than ∆E, in sharp contrast to EL of Figure 2. For example, ES +POL (-23.60 kcal mol-1) is 3.50 kcal mol-1 stronger than ∆E (-20.10 kcal mol-1) at 2.75 Å, whereas EL (-20.13 kcal mol-1) is nearly identical with ∆E at that distance. Recall that ES + POL and EL differ by the SE contribution of the latter (eq 2). Instead, SE is retained in the DEF component of the five-term partition. Thus,
J. Phys. Chem., Vol. 100, No. 43, 1996 17155
Figure 4. Three-term NEDA (aug-cc-pVTZ) of the interaction potential for (H2O)2.
the long-range behavior of the interaction potential is reflected in three components (ES, POL, DEF) of the five-term partition whereas it is isolated in a single component (EL) of the threeterm partition. B. (H2O)2. We next turn to the weak hydrogen-bonding interaction of the water dimer. RHF/NEDA calculations were performed with the aug-cc-pVTZ basis set for the structure shown below:
The monomer geometries were constrained to experimentally determined parameters [R(OH) ) 0.9572 Å, ∠(HOH) ) 104.52°]21 with linear O‚‚‚H-O hydrogen bond. The orientation of the electron donor monomer relative to the acceptor was also constrained so that a 120° angle subtends the bisector of ∠(HOH) of the donor unit and the O-O line of centers. The interaction energy is -3.59 kcal mol-1 at equilibrium separation (∼3 Å), nearly identical with the estimated complete basis set limit (-3.55 kcal mol-1) for the RHF method.26 Figure 4 shows the three-term NEDA of the dimer interaction potential. The five-term analysis of the dimer is qualitatively similar to Figure 4 and will not be discussed here. The leading electrical term is the dipole-dipole contribution, so the potential is considerably more shallow than that of Li+(H2O). The longrange electrical nature of the dimer potential is again clearly reflected by the NEDA treatment. ∆E and EL differ by less than 0.05 kcal mol-1 for all O-O distances greater than 4.5 Å. For shorter distances, ∆E and EL diverge somewhat as the CORE and CT interactions become significant. Charge transfer and electrical interactions both contribute importantly in the equilibrium geometry. At an O-O separation of 3.0 Å, the CT term (-6.06 kcal mol-1) is considerably stronger than ∆E (-3.59 kcal mol-1) and 73% as strong as EL (-8.29 kcal mol-1). CT effects increase rapidly with decreasing O-O distance and, in fact, are stronger than EL at distances slightly less than the equilibrium distance. For example, at 2.8 Å, CT (-12.38 kcal mol-1) exceeds EL (-11.64 kcal mol-1) by about 0.7 kcal mol-1. The essential contribution of CT in hydrogen-bonding interactions and other weakly bonded systems has been stressed elsewhere.14 In the water dimer, CT arises from the transfer of electron density from an oxygen lone pair of the donor unit into the proximal O-H antibond of the acceptor. A perturbative estimate of the strength of this
17156 J. Phys. Chem., Vol. 100, No. 43, 1996 interaction (-5.36 kcal mol-1 at 3.0 Å) suggests that it is primarily responsible for the NEDA CT component (-6.06 kcal mol-1). IV. Summary Natural energy decomposition analysis (NEDA) is a numerically stable and efficient computational tool for analyzing ab initio interaction potentials. In the present work, we applied linear response theory to NEDA to evaluate the electrical self polarization energies of interacting molecules. These energies, together with the NEDA components described in previous work,1,2 facilitate the reduction of the potential into three terms: electrical interaction (EL), core repulsion (CORE), and charge transfer (CT). Applications to Li+(H2O) and water dimer demonstrate that EL, resulting from the classical static and induced interactions of the molecular charge densities, is entirely responsible for the long-range behavior of the interaction potential. Quantum mechanical CORE (σ-σ) and CT (σ-σ*) effects only contribute importantly at short-range (typically within about 1 Å of equilibrium separation) where charge densities overlap considerably. NEDA should be particularly valuable for guiding the development of empirical potentials and mixed quantum mechanics/molecular mechanics (QM/MM) methods when consistency with ab initio results is desired. By fitting potentials to reproduce the separate physical components of NEDA rather than the total interaction energy, it may be possible to circumvent ambiguities due to underdetermined and highly correlated parameters. Atomic charges and polarizabilities can be selected to reproduce the EL component. As shown previously,1,2,27 multipole moments of the monomer charge densities provide a useful decomposition of the EL energy. We have shown in the present work that EL is reasonably well approximated by a classical electrical treatment, assuming instantaneous polarization response. The CORE and CT terms could be effectively modeled through appropriately constructed exponentials that account for the overlap of σ-σ and σ-σ* orbitals on neighboring molecules. Acknowledgment. This work was supported by the Office of Basic Energy Sciences of the U.S. Department of Energy under Contract No. DE-AC06-76RLO1830. E.D.G. also acknowledges the support of Associated Western Universities, Inc.
Schenter and Glendening under Grant No. DE-FG06-89ER-75522 with the U.S. Department of Energy and the Camille and Henry Dreyfus Faculty Start-up Grant Program for Undergraduate Institutions. Pacific Northwest National Laboratory is operated by Battelle Memorial Institute. The authors gratefully acknowledge Dr. Rick A. Kendall for his thorough review of the manuscript. References and Notes (1) Glendening, E. D.; Streitwieser, A. J. Chem. Phys. 1994, 100, 2900. (2) Glendening, E. D. J. Am. Chem. Soc. 1996, 118, 2473. (3) Morokuma, K. J. Chem. Phys. 1971, 55, 1236. (4) Kitaura, K.; Morokuma, K. Int. J. Quantum Chem. 1976, 10, 325. (5) Morokuma, K. Acc. Chem. Res. 1977, 10, 294. (6) Glendening, E. D.; Feller, D.; Thompson, M. A. J. Am. Chem. Soc. 1994, 116, 10657. (7) Thompson, M. A.; Glendening, E. D.; Feller, D. J. Phys. Chem. 1994, 98, 10465. (8) More, M. B.; Glendening, E. D.; Ray, D.; Feller, D.; Armentrout, P. B. J. Phys. Chem. 1996, 100, 1605. (9) Glendening, E. D.; Feller, D. J. Phys. Chem. 1995, 99, 3060. (10) Glendening, E. D.; Feller, D. J. Phys. Chem. 1996, 100, 4790. (11) Glendening, E. D.; Feller, D. J. Am. Chem. Soc. 1996, 118, 6052. (12) Foster, J. P.; Weinhold, F. J. Am. Chem. Soc. 1980, 102, 7211. (13) Reed, A. E.; Weinstock, R. B.; Weinhold, F. J. Chem. Phys. 1985, 83, 735. (14) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. ReV. 1988, 88, 899. (15) Bottcher, C. J. F.; Bordewijk, P. Theory of Electric Polarization; Elsevier: Amsterdam, 1978. (16) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry; McGrawHill: New York, 1989. (17) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, J.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., Jr. J. Comput. Chem. 1993, 14, 1347. (18) Glendening, E. D.; Badenhoop, J. K.; Reed, A. E.; Carpenter, J. E.; Weinhold, F. NBO 4.0; Theoretical Chemistry Institute: Madison, WI, 1994. (19) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (20) Feller, D.; Glendening, E. D.; Kendall, R. A.; Peterson, K. A. J. Chem. Phys. 1994, 100, 4981. (21) Benedict, W. S.; Gailar, N.; Plyler, E. K. J. Chem. Phys. 1956, 24, 1139. (22) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (23) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796. (24) Woon, D. E.; Dunning, T. H., Jr. To be published. (25) Electrical properties of H2O were calculated at the RHF/aug-ccpVTZ level of theory for the constrained geometry described in the text. (26) Feller, D. J. Chem. Phys. 1992, 96, 6104. (27) In particular, see Figure 8 of ref 2 for the Li+(H2O) system and Figure 4 of ref 1 for the water dimer.
JP9612994
