In many cases, a computational chemist calculates the energy
change Δ*E* of an interesting process
(e.g., conformational transition or formation of an intermolecular complex),
hoping to identify the energetic "components" of distinct physical
origin. Contributions to such
"energetic decomposition analysis" (EDA) may include:

- "Electrical" (EL) component, considered to represent the classical-like Coulombic interactions between atomic charges, bond dipoles and the like;
- "Steric exchange" (EX) component, considered to represent Pauli exchange-type repulsions between filled orbitals (or the quasi-classical "Lennard-Jones repulsion" between hard-shell sphere atoms);
- "Charge transfer" (CT) component, considered to represent resonance-type "delocalization" interactions between filled orbitals of one subunit and unfilled orbitals of the other.

The EL component may be further divided into static (ES) and induced "polarization" (POL) components [with associated "self-energy" (SE) correction for the energetic penalty at each polarizing center], and post-HF "correlation" (CORR) or other components may be added in more refined treatments. However, the contentious literature of conflicting EDA schemes centers mainly on the magnitude of quantum-mechanical CT vs. quasi-classical EL + EX components in additive decompositions of the form

Although there is widespread agreement on the general form of Eq. (1) and
the quantum mechanical integrals contributing to overall Δ*E*,
sharp disagreements persist on how individual integrals should be
regrouped for labelling as "EL", "EX", "CT" contributions.

The pioneer EDA formulation of Morokuma and coworkers [1] (MEDA) was initially developed for interacting A, B monomers of a supramolecular A···B complex, but similar schemes were later adopted for torsional interactions of alkyl radicals about a connecting bond [2]. The common assumption of these formulations is that the quasi-classical EL + EX contributions can be evaluated as

Here, ψ_{A}, ψ_{B} are antisymmetrized
determinantal wavefunctions for isolated A, B subunits,
placed at equilibrium A···B separation in the Hamiltonian
integral. The EL component
is commonly evaluated as a classical Coulomb integral between the
corresponding electron densities ρ_{A}, ρ_{B}
at the two centers, viz.,

and remaining
EX^{(MEDA)},
CT^{(MEDA)} components are then determined
by difference from Eqs. (1), (2). In particular, the CT
component is evaluated as the "remainder" term

independent of how
EL^{(MEDA)} is evaluated.

The integral
<ψ_{A}|H|ψ_{B}>_{eq}
superficially resembles the perturbative interaction integral
between unperturbed
ψ_{A}, ψ_{B}
wavefunctions of the non-interacting long-range limit. Indeed, as
*R*_{AB} → ∞, this interpretation is valid
and Eqs. (2)-(3) become *exact*. In the same limit,
ψ_{A}, ψ_{B} are mutually orthogonal and
CT^{(MEDA)} vanishes, in accordance
with London's long-range perturbation theory [3].

However, at the actual equilibrium separation of
the A···B complex,
attempted physical justification of Eqs. (2)-(4) in perturbation-theoretic terms
is certainly invalid. At any finite equilibrium
separation, the supposed "unperturbed wavefunctions"
ψ_{A}, ψ_{B} are *non*orthogonal [4]

and *no possible* physical (Hermitian) "unperturbed Hamiltonian"
could underlie such interpretation (because non-degenerate
solutions of Hermitian
eigenvalue problems are necessarily orthogonal). For H-bonded complexes,
equilibrium *R*_{AB} separations are commonly found to be
0.5-1.0 Ångstrom *in*side the formal van der Waals contact distance
at which exchange-type overlap is expected to
increase *exponentially*, and for torsional interactions the
overlap issues are still more severe. MEDA-style components are therefore
seen to be inherently dependent on AO overlap
(required, e.g., for any non-zero
evaluation of CT^{(MEDA)}),
exemplifying "overlap-dependent" interpretation
(and labeling) of CT.

Why is AO overlap a critical issue for evaluation of CT? As shown
in Eq. (4), MEDA evaluation of CT originates not in physical identification
of donor-acceptor (filled-unfilled) orbital interactions, but in the
overlap-dependent "remainder" of a long-range identity that is applied improperly
in the region of significant exchange-type interactions. This
contrasts sharply with NBO-based EDA techniques, which are always formulated in the
*overlap-free* domain of NAOs, assuring strict probability
conservation and consistent natural charge assignments throughout
subsequent hybridization, covalent bond formation, and supramolecular
complex formation. Disagreements between
MEDA and NBO-based EDA assignments of CT can therefore be traced squarely
to disparate handling of AO overlap [5].

The intrinsic ambiguities of overlap-dependent charge assignments strongly affect how electrons are perceived to be distributed between one atomic center or another (as well known, e.g., in the pathologies of Mulliken population analysis). Morokuma-type evaluation of CT is intrinsically tied to ambiguous AO overlap terms that confuse assignments of atomic charge, resulting in profound disagreements over "transfer" between one monomer or the other. Related Ziegler-Rauk, Bickelhaupt-Baerends, BLW, and ALMO-EDA methods [6] handle overlap-density terms in slightly different ways, but all differ sharply from overlap-free NBO methods in this respect.

While NCE, STERIC, and $DEL keywords allow quantitation of NBO-based EL-, EX-, and CT-type energy contributions, they do not account for expected non-additive coupling terms between these components (corresponding to "higher-order effects" in perturbative parlance). Important coupling effects include:

- geometric and electronic "deformations" in complex formation;
*CT-induced*electronic shifts that significantly alter apparent EL/EX component values;*R*-dependent basis set superposition (BSSE) artifacts (often significant in non-augmented AO basis sets).

To incorporate such higher-order effects consistently in traditional EDA-style additive components, Glendening and coworkers [7] developed the "Natural Energy Decomposition Analysis" (NEDA) method. The NEDA keyword requires an intricate sequence of interactive calculations between NBO6 and host ESS, and is currently implemented only in GAMESS, FIREFLY (PC-GAMESS), and NWCHEM host systems.

Similar to Eq. (1), NEDA partitions the
usual expression for the binding energy Δ*E*
of the A···B complex

into a sum of "electrical" (EL), "core" (CORE), and "charge transfer" (CT) components which may be abbreviated as

All NEDA energy evaluations in Eq. (6) employ the *full* basis set of
the complex, so that Δ*E* automatically
incorporates the full Boys-Bernardi "counterpoise
correction" [8] for BSSE.

The key step of NEDA analysis is evaluation of perturbed ("deformed")
wavefunctions
(ψ_{A}^{(def)},
ψ_{B}^{(def)}) for each monomer
as local eigenvectors of the NBO Fock matrix
in the respective monomer blocks. (Each monomer wavefunction
ψ_{A}^{(def)}, ψ_{B}^{(def)}
thereby inherits the overlap-free character of the underlying NAOs.) The
antisymmetrized (determinantal) product of perturbed fragment wavefunctions

then defines the starting "localized" (ψ^{(loc)}) wavefunction
that underlies evaluation of each NEDA component.

Consistent with its physical description as "delocalization" from
each fragment into unfilled orbitals of the other, the
charge transfer (CT) component is defined *directly* as

rather than, as in Eq. (4), a "remainder" term. As usual, the energy difference
between ψ^{(loc)} and the perturbed fragment wavefunctions
ψ_{A}^{(def)}, ψ_{B}^{(def)} is
attributed to a sum of electrostatic (ES), polarization (POL) and
exchange (EX) contributions, viz.,

The "deformation" energy (DEF) of forming each perturbed
localized ψ_{A}^{(def)} from its starting fragment
wavefunction ψ_{A} is evaluated in a similarly direct
manner as

corresponding to combined electric field and quantum mechanical effects experienced by each fragment within the complex.

In accordance with classical electrodynamics, one can employ linear response theory to evaluate the self-energy (SE) correction ("polarization penalty") for each center, thereby expressing the total electrical component as

However, the self-energy must then be removed (subtracted) from the deformation (DEF) penalty to obtain the net repulsive CORE component, viz.,

By combining the definitions for EL [Eq. (12)], CORE [Eq. (13)],
and CT [Eq. (9)], one readily verifies that these components satisfy
the basic NEDA *ansatz*, Eq. (7).

In addition to the energies *E*(ψ_{A}),
*E*(ψ_{A}^{(def)}), NEDA also calculates
the dipole moments (μ_{A},
μ_{A}^{(def)})
for the wavefunctions ψ_{A},
ψ_{A}^{(def)} of each fragment. The difference
between these two vector quantities is a measure of the induced
dipole due to complex formation, which asymptotically reduces to
the corresponding classical expression in the long-range limit.

NEDA differs significantly from other NBO keyword options. The
"NEDA" keyword is *only* allowed within a $DEL keylist (not in the usual
$NBO keylist) and must always appear with accompanying
"END" in keylist-like "NEDA [fragment units] END"
form. For the sake of specificity in the tutorial instructions
below, we assume that GAMESS is the
host ESS program and that the user has some general
familiarity with GAMESS keyword groups and syntax.

Full details of "[fragment units]" specification will be described below. However, in the simplest and most common closed-shell application, when the desired NEDA fragments A, B,... are the same as the "molecular units" returned by default NBO analysis, default NEDA analysis could be run merely by appending the following lines to the GAMESS input deck:

$NBO $END $DEL NEDA END $END

For an open-shell species, the corresponding $DEL keylist could be of the form

$DEL ALPHA NEDA END BETA NEDA END $END

again implying that the desired fragments A, B,... are those of default
NBO search. (Note however that the program will halt if the molecular units
of the alpha spin set are different from those of the beta spin
set, so more detailed specification of NEDA fragment units is necessary
in this case.) One can also choose to perform NEDA with alternative
NAO rather than default NBO basis set, in which case
the default "units" are individual *atoms* rather than
bonded molecular units of the NBO search (see below).

NEDA fragment wavefunctions can be evaluated by any $SCF group options available in GAMESS input, including SOSCF, DIIS, damping, and direct energy minimization (DEM) methods to control convergence. As in other $DEL applications, use of symmetry should generally be disabled by setting NOSYM=1 in the GAMESS $CONTRL group and NOPK=1 in the $INTGRL group.

Several options are available for specifying the NEDA fragments, if different from default NBO molecular units. The most commonly useful options are:

**Parenthesized Unit-Lists.**In the special case that one wishes to*combine*multiple molecular units (as defined and numbered in the "NBO Summary" output section) into composite NEDA fragments, one can insert a parenthesized list of units to be included in each fragment in the NEDA ... END input. For example, if default NBO analysis returned six molecular units (numbered 1-6 in the NBO Summary listing) that one wished to consider as four fragments [A = (1,3), B = (2,6), C = (4), D = (5)] for NEDA purposes, one could use the NBO input lines shown below:

$NBO $END $DEL NEDA (1,3) (2,6) (4) (5) END $END

**$CHOOSE specification.**Often the most convenient option is to specify desired alternative molecular units by a $CHOOSE keylist, in which case the simple "NEDA END" again suffices. For example, if one wishes to analyze BH_{3}NH_{3}in terms of A = BH_{3}, B = NH_{3}fragments, but the default NBO search returns a single molecular unit with B1-N2, B1-H3, B1-H4, B1-H5, N2-H6, N2-H7, N2-H8 connectivity, the NBO input lines could be:

$NBO $END $CHOOSE LONE 2 1 END BOND S 1 3 S 1 4 S 1 5 S 2 6 S 2 7 S 2 8 END $END $DEL NEDA END $END

**NAO/atom-list specification.**An alternative option is to request that NEDA be performed in the NAO (rather than default NBO) basis, as specified by a command of the form "NAO NEDA [unit-lists] END". In this case the default "units" are individual atoms, and each parenthesized "unit-list" specifies the atoms contained in the desired NEDA fragment. This is often the simplest and most flexible option for non-standard molecular units. For example, the BH_{3}NH_{3}analysis described above could be requested alternatively by input lines$NBO $END $DEL NAO NEDA (1,3,4,5) (2,6,7,8) END $END

Other $DEL command keywords can be included as usual in the keylist, including alternative NEDA fragment analyses. Each such NEDA or non-NEDA command will be processed in turn for output display. See Sec. B.5 of the NBO Manual for further information on $DEL keylist syntax.

Here's a sample GAMESS input deck for RHF/4-31G calculation of the water dimer (far too low a theoretical level for current "real world" application):

$CONTRL SCFTYP=RHF RUNTYP=ENERGY COORD=HINT NOSYM=1 $END< $INTGRL NOPK=1 $END $BASIS GBASIS=N31 NGAUSS=4 $END $NBO $END $DEL NEDA END $END $DATA Water dimer...(rhf/4-31g//from Umeyama and Morokuma, JACS 99, 1316 (1977) Cs Oxygen 8.0 LC 0.0000 0.0000 0.0000 + O I Hydrogen 1.0 TCT 0.9560 52.6000 90.0000 + O I J Oxygen 8.0 TCT 2.8800 120.0000 0.0000 + O I J Hydrogen 1.0 TCT 0.9560 105.2000 180.0000 + 3 O J Hydrogen 1.0 TCT 0.9560 105.2000 0.0000 + 3 4 O $END

As shown in the $DEL keylist, NEDA will perform default analysis using the two molecular units (water monomers) of the NBO search.

The first segment of NEDA output is shown below, essentially equivalent to usual $DEL output for deletion of all intermolecular delocalizations between the two units (deletion type 8; cf. Sec. B.5 of the NBO Manual):

NEDA: Natural Energy Decomposition Analysis Deletion of the NBO Fock matrix elements between orbitals: 1 3 4 7 8 11 12 15 16 17 18 19 20 and orbitals: 2 5 6 9 10 13 14 21 22 23 24 25 26 Orbital occupancies: Orbital No deletions This deletion Change ------------------------------------------------------------------------------ 1. CR ( 1) O 1 1.99982 1.99983 0.00001 2. CR ( 1) O 4 1.99979 1.99980 0.00000 3. LP ( 1) O 1 1.99920 1.99930 0.00010 4. LP ( 2) O 1 1.98029 1.99914 0.01885 5. LP ( 1) O 4 1.99994 2.00000 0.00006 6. LP ( 2) O 4 1.99922 1.99932 0.00010 7. BD ( 1) O 1- H 2 1.99912 1.99942 0.00030 8. BD ( 1) O 1- H 3 1.99912 1.99942 0.00030 9. BD ( 1) O 4- H 5 1.99960 1.99961 0.00000 10. BD ( 1) O 4- H 6 1.99919 1.99943 0.00024 11. BD*( 1) O 1- H 2 0.00032 0.00016 -0.00016 12. BD*( 1) O 1- H 3 0.00032 0.00016 -0.00016 13. BD*( 1) O 4- H 5 0.00027 0.00008 -0.00019 14. BD*( 1) O 4- H 6 0.01827 0.00014 -0.01814 15. RY ( 1) O 1 0.00004 0.00001 -0.00003 16. RY ( 2) O 1 0.00002 0.00000 -0.00002 17. RY ( 3) O 1 0.00001 0.00000 -0.00001 18. RY ( 4) O 1 0.00000 0.00001 0.00001 19. RY ( 1) H 2 0.00124 0.00127 0.00003 20. RY ( 1) H 3 0.00124 0.00127 0.00003 21. RY ( 1) O 4 0.00009 0.00000 -0.00009 22. RY ( 2) O 4 0.00002 0.00000 -0.00002 23. RY ( 3) O 4 0.00001 0.00000 0.00000 24. RY ( 4) O 4 0.00000 0.00000 0.00000 25. RY ( 1) H 5 0.00127 0.00120 -0.00007 26. RY ( 1) H 6 0.00159 0.00042 -0.00117 NEXT STEP: Evaluate the energy of the new density matrix that has been constructed from the deleted NBO Fock matrix by doing one SCF cycle. ------------------------------------------------------------------------------ -------------------------- RHF SCF CALCULATION -------------------------- DENSITY MATRIX CONV= 1.00E-05 ITER EX DEM TOTAL ENERGY E CHANGE DENSITY CHANGE ORB. GRAD 1 0 0 -151.8065596360 -151.8065596360 0.061186977 0.000000000 SCF IS UNCONVERGED, TOO MANY ITERATIONS TIME TO FORM FOCK OPERATORS= 0.0 SECONDS ( 0.0 SEC/ITER) TIME TO SOLVE SCF EQUATIONS= 0.0 SECONDS ( 0.0 SEC/ITER) FINAL RHF ENERGY IS 0.0000000000 AFTER 1 ITERATIONS ...... END OF RHF CALCULATION ...... STEP CPU TIME = 0.15 TOTAL CPU TIME = 0.4 ( 0.0 MIN) TOTAL WALL CLOCK TIME= 0.5 SECONDS, CPU UTILIZATION IS 93.48% ------------------------------------------------------------------------------ Energy of deletion : -151.806559636 Total SCF energy : -151.827680724 ------------------- Energy change : 0.021121 a.u., 13.254 kcal/mol ------------------------------------------------------------------------------

As shown in the final line, the calculated energy change for
this deletion is 13.254 kcal/mol. This corresponds to the
variational *loss* of stabilizing charge-transfer delocalizations,
so the CT component is simply the negative of this quantity, viz.,

This $DEL estimate of CT is in sensible qualitative agreement
with the corresponding
perturbative estimate (Δ*E*^{(2)} = -15.21 kcal/mol) for the
dominant *n*_{O(1)} → σ*_{O(4)-H(6)}
interaction (as usual, Δ*E*^{(2)} somewhat overestimates
the variational $DEL energy change, due to neglect of higher-order
coupling terms). From the $DEL occupancy changes, one can
also see that localizing
the wavefunction back-transfers about 0.018 electrons from NBO 14
(the σ*_{O(4)-H(6)} antibond acceptor) to NBO 8 (the
*n*_{O(2)} lone pair donor), effectively suppressing
the leading CT interaction of the dimer.

[Note that the warning message *"SCF is unconverged, too many iterations"*
should be ignored, because the single-pass $DEL
evaluation method is based on interrupting SCF
iterations after the first cycle to evaluate the energy for the input
$DEL-density, not the converged SCF density.]

The next segment of output reports on calculation of
the perturbed (ψ_{A}^{(def)})
and variationally optimized (ψ_{A}) wavefunctions
for the first fragment, the H(2)-O(1)-H(3) monomer:

-------------- Fragment 1: -------------- ...... END OF ONE-ELECTRON INTEGRALS ...... STEP CPU TIME = 0.03 TOTAL CPU TIME = 0.5 ( 0.0 MIN) TOTAL WALL CLOCK TIME= 0.5 SECONDS, CPU UTILIZATION IS 93.88% -------------------------- RHF SCF CALCULATION -------------------------- DENSITY MATRIX CONV= 1.00E-05 ITER EX DEM TOTAL ENERGY E CHANGE DENSITY CHANGE ORB. GRAD 1 0 0 -75.8820038771 -75.8820038771 0.051585286 0.000000000 ---------------START SECOND ORDER SCF--------------- 2 1 0 -75.9081413953 -0.0261375181 0.018452073 0.009561181 3 2 0 -75.9092479375 -0.0011065422 0.009032407 0.003595518 4 3 0 -75.9094058121 -0.0001578747 0.000764754 0.000898516 5 4 0 -75.9094124510 -0.0000066389 0.000263894 0.000128113 6 5 0 -75.9094127399 -0.0000002890 0.000084067 0.000046344 7 6 0 -75.9094127514 -0.0000000114 0.000018443 0.000008466 8 7 0 -75.9094127518 -0.0000000005 0.000003845 0.000002277 9 8 0 -75.9094127518 0.0000000000 0.000001199 0.000000691 ----------------- DENSITY CONVERGED ----------------- TIME TO FORM FOCK OPERATORS= 0.0 SECONDS ( 0.0 SEC/ITER) TIME TO SOLVE SCF EQUATIONS= 0.0 SECONDS ( 0.0 SEC/ITER) FINAL RHF ENERGY IS -75.9094127518 AFTER 9 ITERATIONS ...... END OF RHF CALCULATION ...... STEP CPU TIME = 0.05 TOTAL CPU TIME = 0.5 ( 0.0 MIN) TOTAL WALL CLOCK TIME= 0.5 SECONDS, CPU UTILIZATION IS 94.44% ...... END OF ONE-ELECTRON INTEGRALS ...... STEP CPU TIME = 0.00 TOTAL CPU TIME = 0.5 ( 0.0 MIN) TOTAL WALL CLOCK TIME= 0.5 SECONDS, CPU UTILIZATION IS 94.44% Dipole (def): 2.6394(x), 0.0793(y), 0.0000(z); 2.6406(tot) Debye Dipole (cp): 2.6541(x), -0.0947(y), 0.0000(z); 2.6557(tot) Debye Dipole (ind): -0.0146(x), 0.1740(y), 0.0000(z); 0.1747(tot) Debye

The calculated expectation
energy of localized ψ_{A}^{(def)} (-75.88200 a.u.)
is obtained by single-pass SCF evaluation
with the localized monomer MOs of the previous step, whereas that of
optimized ψ_{A} (-75.90941 a.u.) is obtained after complete SCF
convergence. Both wavefunctions are calculated in the *full* dimer basis
at the fixed monomer geometry in the complex, thereby incorporating the complete
counterpoise (cp) correction for this rather small basis.

The final lines of this segment display the calculated dipole moment for the perturbed (def) and optimized (cp) wavefunctions, showing the relatively small shift (induced "ind" dipole = 0.1747 Debye) from 2.6406 in the isolated monomer to 2.6557 Debye in the presence of the other monomer. In effect, the strict orthogonality required of the two monomer wavefunctions prevents the electron distribution of the donor from significantly penetrating the acceptor.

A similar report then follows for the acceptor fragment H(5)-O(4)-H(6):

-------------- Fragment 2: -------------- ...... END OF ONE-ELECTRON INTEGRALS ...... STEP CPU TIME = 0.04 TOTAL CPU TIME = 0.5 ( 0.0 MIN) TOTAL WALL CLOCK TIME= 0.6 SECONDS, CPU UTILIZATION IS 96.49% -------------------------- RHF SCF CALCULATION -------------------------- DENSITY MATRIX CONV= 1.00E-05 ITER EX DEM TOTAL ENERGY E CHANGE DENSITY CHANGE ORB. GRAD 1 0 0 -75.8962673806 -75.8962673806 0.016278993 0.000000000 ---------------START SECOND ORDER SCF--------------- 2 1 0 -75.9079696846 -0.0117023040 0.007570847 0.005067234 3 2 0 -75.9083141028 -0.0003444182 0.005463348 0.001684405 4 3 0 -75.9083555888 -0.0000414859 0.001012338 0.000672210 5 4 0 -75.9083571171 -0.0000015284 0.000352290 0.000297195 6 5 0 -75.9083574089 -0.0000002918 0.000092465 0.000043124 7 6 0 -75.9083574187 -0.0000000097 0.000014019 0.000006844 8 7 0 -75.9083574192 -0.0000000005 0.000004722 0.000002775 9 8 0 -75.9083574192 -0.0000000001 0.000000684 0.000000532 ----------------- DENSITY CONVERGED ----------------- TIME TO FORM FOCK OPERATORS= 0.0 SECONDS ( 0.0 SEC/ITER) TIME TO SOLVE SCF EQUATIONS= 0.0 SECONDS ( 0.0 SEC/ITER) FINAL RHF ENERGY IS -75.9083574192 AFTER 9 ITERATIONS ...... END OF RHF CALCULATION ...... STEP CPU TIME = 0.04 TOTAL CPU TIME = 0.6 ( 0.0 MIN) TOTAL WALL CLOCK TIME= 0.6 SECONDS, CPU UTILIZATION IS 95.16% ...... END OF ONE-ELECTRON INTEGRALS ...... STEP CPU TIME = 0.01 TOTAL CPU TIME = 0.6 ( 0.0 MIN) TOTAL WALL CLOCK TIME= 0.6 SECONDS, CPU UTILIZATION IS 96.77% Dipole (def): -0.8286(x), -2.6607(y), 0.0000(z); 2.7867(tot) Debye Dipole (cp): -1.0045(x), -2.3894(y), 0.0000(z); 2.5919(tot) Debye Dipole (ind): 0.1759(x), -0.2713(y), 0.0000(z); 0.3234(tot) Debye

In this case, the induced dipole for the acceptor is seen to be
rather sizable
(0.3234 Debye), reflecting the significant polarization of
σ_{O(4)H(6)} away from the donor
(thereby better exposing σ*_{O(4)H(6)} to nucleophilic
attack), as seen in the NBO polarization coefficients. Nevertheless,
as will be seen in the energy summary below, such polarization (induction)
changes make only secondary contributions to H-bond energy.

The final segment of NEDA output summarizes the final CT, ES, POL, EX component values for the complex, as well as the deformation (DEF) and self-energy (SE) values for each monomer:

Natural Energy Decomposition Analysis (Summary): Component Energy(wfn) Energy(wfn) (kcal/mol) ------------------------------------------------------------------------------ H4O2 -151.8276807(scf) -151.8065596(loc) CT = -13.25 ES = -12.24 POL = -3.89 EX = -1.62 1. H2O -75.8820039(def) -75.9094128(cp) DEF(SE) = 17.20( 0.10) 2. H2O -75.8962674(def) -75.9083574(cp) DEF(SE) = 7.59( 1.78) --------- E = -6.22 Electrical (ES+POL+SE) : -14.25 Charge Transfer (CT) : -13.25 Core (EX+DEF-SE) : 21.29 ------------ Total Interaction (E) : -6.22 ......END OF NBO ANALYSIS......

As shown in the final lines, the combined electrical (EL: -14.25 kcal/mol), charge transfer (CT: -13.25 kcal/mol), and core (CORE: +21.29 kcal/mol) contributions add up to the calculated total interaction energy, -6.22 kcal/mol.

In this case, CT (-13.25 kcal/mol) and ES (-12.24 kcal/mol) are the largest individual contributors to binding, with POL (-3.89 kcal/mol) playing an important secondary role. These strong attractions are necessary to overcome the strong CORE repulsions (+21.29 kcal/mol) at the equilibrium geometry of H-bonding (far inside van der Waals contact), resulting in the final net H-bond energy of -6.22 kcal/mol.

From the introductory discussion, we should generally expect strong differences
between NEDA (overlap-free) and MEDA-type (overlap-dependent) decompositions, and
the water dimer proves to be no exception. Through the hidden AO overlap terms [5]
that enter the <ψ_{A}|H|ψ_{B}> integral
of Eq. (2), the critical *n*_{O(1)} → σ*_{O(4)H(6)}
CT interaction is largely absorbed into "EL + EX" terms in MEDA-type
partitioning. CT^{(MEDA)} is then obtained [Eq. (4)] as a
feeble "remainder" term, evaluated by Umeyama and Morokuma [1] as
CT^{(MEDA)} = -2.1 kcal/mol for RHF/4-31G water dimer
(in slightly different geometry).

The near *order of magnitude* disagreement between what
NEDA vs. MEDA partisans consider
to be the accurately evaluated "CT" naturally inspires impassioned debate! (The
"CT vs. electrostatics" debate
in the H-bond case is mirrored
by corresponding "CT vs. sterics" debate in the rotation barrier case.) Still
greater disagreements separate *both* methods from the "Intermolecular
Perturbation Theory" of Stone and coworkers [9]. According to the latter, CT
is not merely of secondary importance, but "ill-defined"
and "part of the induction (polarization)" component
that "vanishes in the limit of a complete basis" [10].

Despite such points of disagreement with older methods, one can see
that NEDA clearly recognizes the importance of EL interactions. Indeed,
from NEDA results presented above, one might conclude superficially that
the largest contributions to H-bonding (ES + POL + SE: -14.25 kcal/mol)
are "electrical" in nature, consistent with the presumption expressed
in many textbooks. However, from
NCE analysis or similar considerations,
one can judge that much of the
quasi-classical EL interaction is *CT-induced*, amplified by the
deep inter-penetration of van der Waals radii that only resonance-type
CT interactions could enable.

Geometry optimization with deletions can be performed more straightforwardly with Gaussian than with GAMESS. [Such optimizations are possible with GAMESS (see the deletions tutorial for sample input) but require the use of RUNTYP=TRUDGE and Hildebrandt-style coordinate input.] Below is a sample Gaussian input file for performing $DEL-type optimization of the water dimer, using constrained monomer geometries. NEDA is not implemented in Gaussian, so the calculation instead "zeroes delocalizations" between the two monomers in the $DEL keylist. The zeroing of delocalizations between monomers is effectively equivalent to the CT localization of NEDA.

#rhf/4-31g pop=nbo6del opt nosymm Water dimer ...(rhf/4-31g//from Umeyana and Mokokuma, JACS 99, 1316 (1977)) CT-deleted structure: E($DEL) = -151.8222285 (dE = -4.36 kcal/mol) 0 1 O H 1 0.956 H 1 0.956 2 105.2 O 1 rhb 2 a214 3 d3214 H 4 0.956 1 a145 2 d2145 H 4 0.956 5 105.2 1 d1546 rhb 2.88 a214 107.68 d3214 114.64 a145 105.20 d2145 123.52 d1546 0.00 $nbo print=0 nbosum $end $del zero 2 deloc from 1 to 2 from 2 to 1 $end

[Note that such $DEL-type optimizations are performed by non-analytic derivative methods that are quite time-consuming, hence restricted to a limited number of optimization variables in practical applications.]

The results are shown below in graphical comparisons of the
Umeyama-Morokuma H-bonded complex (upper) with the no-CT counterpart
(lower). From visual comparison, one can see that the no-CT complex
differs strongly in both radial and angular characteristics from the
actual H-bonded dimer. The no-CT structure reverts to the
C_{2v}-symmetric structure expected for a classical
"dipole-dipole complex", with *no* tendency toward linear or
near-linear O···H-O alignment. One can also see
that H_{2}O···H_{2}O separation
increases by more than 0.3 Å to approximate van der Waals
contact distance.

Consistent with
*R*_{O···O} elongation,
the net binding energy of the no-CT complex weakens
by about 3.4 kcal/mol, to less than *half* the value of realistic
H-bonding. Qualitatively similar differences between the geometry and
energetics of full vs. no-CT complexes
persist at all higher levels of theory.

At the (unphysical) no-CT geometry, the neglected CT interaction energy
(ca. 0.001% of
*E*_{total}) is virtually negligible, and only the EL +
CORE components remain to exert significant forces on molecular
geometry. Although EL attraction appeared only slightly less
important than CT in realistic H-bond geometry, the $DEL calculations
show clearly that electrostatic forces alone are *un*able to
bring the monomers to the characteristic short distances, near-linear
O···H-O alignments, and robust binding energies
of authentic H-bond complexes. The resonance-type CT interaction is
thereby identified as the unique "smoking gun" necessary to transform
feeble dipole-dipole complexes to recognizable H-bonded species.
Similar results are obtained [11] across the entire gamut of H-bond
energies and types.

[1]
K. Kitaura and K. Morokuma, Int. J. Quantum Chem. **10**, 325 (1976);
H. Umeyama and K. Morokuma, J. Am. Chem. Soc. **99**, 1316 (1977);
K. Morokuma, Acc. Chem. Res. **10**, 294 (1977);
K. Morokuma and K. Kitaura, in H. Ratajczak, W. J. Orville-Thomas (eds.),
*Molecular Interactions* (John Wiley, New York, 1980), Vol. 1, pp. 21-66.

[2]
K. Morokuma and H. Umeyama, Chem. Phys. Lett. **49**, 333 (1977);
F. M. Bickelhaupt and E. J. Baerends, Angew. Chem. Int. Ed. **42**,
4183 (2003).

[3]
F. London, Z. Physik. Chem. **B11**, 222 (1930);
R. Eisenschitz and F. London, Z. Phys. **60**, 491 (1930);
F. London, Trans. Faraday Soc. **33**, 8 (1937);
cf. F. Weinhold and C. R. Landis, *Valency and Bonding* (Cambridge
U. Press, New York, 2005), p. 585ff.

[4]
It may be questioned how ψ_{A}, ψ_{B}
can be nonorthogonal when they are expanded in AOs on completely
disjoint sets of atoms. The answer is that basis AOs on one atom
commonly have overlap with AOs on *all* other atoms at finite
separation, so electron density associated with AO-overlap regions
cannot be uniquely assigned as being either "on A" or "on B". For
further details on the relationship of nonorthogonal AOs to
orthonormal NAOs, see the
"NAOs in the AO basis" sample output
from the MATRIX keyword.

[5]
C. T. Corcoran and F. Weinhold, J. Chem. Phys. **72**, 2866 (1980);
A. E. Reed, L. A. Curtiss, and F. Weinhold, Chem. Rev. **88**, 899 (1988);
F. Weinhold and J. E. Carpenter, J. Mol. Struct. (Theochem) **165**,
189 (1988);
F. Weinhold, Angew. Chem. Int. Ed. **42**, 4188 (2003).

[6]
Ziegler-Rauk: T. Ziegler and A. Rauk, Inorg. Chem.**18**, 1558 (1979);

Bickelhaupt-Baerends: F. M. Bickelhaupt and E. J. Baerends,
Rev. Comput. Chem. **15**, 1 (2000);
C. Esterhuysen and G. Frenking, Theor. Chem. Acc. **111**, 381 (2004);

BLW: Y. Mo and S. D. Peyerimhoff, J. Chem. Phys. **112**, 5530 (2000);
Y.Mo, L. Song, and Y. Lin, J. Phys. Chem. A **111**, 8291 (2007);
K Nakashima, X. Zhang, M. Xiang, Y. Lin, M. Lin, and Y. Mo,
J. Theo. Chem. Comput. **7**, 639 (2008);

ALMO-EDA: R. Z. Khaliullin, R. Lochan, E. Cobar, A. T. Bell, and
M. Head-Gordon, J. Phys. Chem. A **111**, 8753 (2007);
R. Z. Khaliullin, A. T. Bell, and M.Head-Gordon, Chem. Eur. J.
**15**, 851 (2008).

[7]
E. D. Glendening and A. Streitwieser Jr., J. Chem. Phys. **100**, 2900 (1994);
G. K. Schenter and E. D. Glendening, J. Phys. Chem. **100**, 17152 (1996);
E. D. Glendening, J. Am. Chem.Soc. **118**, 2473 (1996).

[8]
S. F. Boys and F. Bernardi, Mol. Phys. **19**, 553 (1970).

[9]
I. C. Hayes and A. J. Stone, Mol. Phys. **53**, 69 (1984);
A. D. Buckingham, P. W. Fowler, and A. J. Stone, Int. Rev. Phys. Chem.
**5**, 107 (1986), and refs. therein.

[10]
A. J. Stone, Chem. Phys. Lett. **211**, 101 (1993).

[11] A. E. Reed, F. Weinhold, L. A. Curtiss, and D. J. Pochatko,
J. Chem. Phys. **84**, 5687 (1985); F. Weinhold and C. R. Landis,
*Valency and Bonding* (Cambridge U. Press, New York, 2005), Sec. 5.2.