A new approach for computing the atom-in-molecule [quantum theory of atoms in molecule (QTAIM)] energies in Kohn–Sham density-functional theory is presented and tested by computing QTAIM energies for a set of representative molecules. In the new approach, the contribution for the correlation-kinetic energy (Tc) is computed using the density-functional theory virial relation. Based on our calculations, it is shown that the conventional approach where atomic energies are computed using only the noninteracting part of the kinetic energy might be in error by hundreds of kJ/mol.

Atomic properties defined within the quantum theory of atoms in molecules (QTAIM) are useful to describe and predict phenomena in molecules and materials in fields ranging from solid state physics and x-ray crystallography to drug design and biochemistry.1 The definition of QTAIM atomic properties is based on quantum mechanical principles and the topology of the electron density.2 At a stationary point on the potential energy surface, the QTAIM atomic energies can be computed from the N-electron correlated wave function, Ψ, using the simple equation,2 

EΩ=TΩN2ΩrΨ(r,r2,,rN)rΨ(r,r2,,rN)dr2,,drNdr=N4Ω[Ψ(r,r2,,rN)r2Ψ(r,r2,,rN)+Ψ(r,r2,,rN)r2Ψ(r,r2,,rN)]dr2,,drNdr.
(1)

(Away from the stationary point there are corrections due to the Feynman forces.2) Here, EΩ is the energy of an atom with volume Ω.

Equation (1) cannot be used directly in Kohn–Sham density-functional theory (KS-DFT) because there the correlated N-electron wave function is not known.3 In this paper we will develop alternative expressions, derived from the DFT virial expressions. These expressions were tested with calculations performed using the Amsterdam Density Functional (ADF) program.4–6 So the QTAIM atomic energies complement the set of atomic properties that can be computed with the fast QTAIM method we recently introduced.7–9 

The virial theorem for a molecule can be expressed as10,11

E=T[ρ]AXAEXA
(2a)
or

V[ρ]T[ρ]=2+1TAXAEXA.
(2b)

Here E is the total molecular electronic energy, T[ρ] is the (exact) kinetic energy functional, Xα are the nuclear coordinates, and V[n] is the potential energy functional.10 For a molecule at its equilibrium geometry Eq. (2) reduces to

E=T[ρ],
(3a)
V[ρ]T[ρ]=2.
(3b)
In the conventional KS formalism, T[ρ] is given by11–13 

T[ρ]=Ts+Tc.
(4)

Ts is the noninteracting kinetic energy which can be expressed by inserting the KS determinant in Eq. (1), or in terms of the KS orbitals ψi(r),

Ts=12i=1Nψi(r)2ψi(r)dr.
(5)

Tc is the correlation-kinetic energy, which can be computed as11 

Tc=Exc[ρ]ρ(r)rVxc(r)dr={εxc(r)+ρ(r)rVxc(r)}dr.
(6)

Exc, Vxc, and εxc are the exchange-correlation (xc) energy, the xc potential, and the xc energy density, respectively. To avoid the (numerical) evaluation of the gradient of Vxc, the last term of the right-hand side of Eq. (6) can be integrated by parts,14,

E={12i=1Nψi(r)2ψi(r)εxc(r)+(3ρ(r)+rρ(r))Vxc(r)}dr,
(7a)
E=κ(r)dr,
(7b)
κ{12i=1Nψi(r)2ψi(r)εxc(r)+(3ρ(r)+rρ(r)Vxc(r))}.
The potential energy functional is given by11,12

V[ρ]=VKS[ρ]Tc=ρ(r){AMZA|rXA|+12ρ(r)|rr|+εxc(r)}drTc.
(8)

From Eqs. (4)–(6), (7a), (7b), and (8), Eq. (3) can be rewritten as

E=(Ts+Tc)=κ(r)dr,
(9a)
γVKSTcTs+Tc=2.
(9b)
Equation (9) establishes the virial theorem in KS-DFT for a molecule at its equilibrium geometry. Now, since the QTAIM atomic volumes form a nonoverlapping partition of the space, R3=AΩA,1,2 the energy of atom A can be defined using Eq. (9a),

EΩAΩAκ(r)drTsΩATcΩAΩAi=1Nψi(r)(122)ψi(r)drΩA{εxc(r)+(3ρ(r)+rρ(r))Vxc(r)}dr,
(10)

and the molecular energy is the sum of the atomic energies,

E=AEΩA.
(11)

No approximations were made in Eqs. (2)–(11) so these equations are exact. For approximate KS calculations, however, Eqs. (3) [and thus Eqs. (9)] are not satisfied exactly because of the finite basis set, the imperfect geometry optimization, and (most importantly) the approximate xc functional. Table I shows our calculations of the virial ratio, γ [Eq. (9b)], for a set of representative molecules. All calculations reported here were performed with a development version of ADF using a TZP Slater basis at the equilibrium geometry. We used the Dirac exchange functional15 with the Vosko–Wilk–Nusair correlation functional16 for the local density approximation (LDA) calculations and the Perdew–Burke–Ernzerhof (PBE) functional17 for the generalized gradient approximation (GGA) calculations. (ADF default parameters were used unless otherwise stated.4–6) From Table I we can see that the value of γ is not strictly 2 but there is a small deviation due to the reasons mentioned above. Table I also shows the value of the ratio between VKS and Ts,

γVKSTs,
(12)

which is also close to 2 because Tc is much smaller (<1%) than either Ts and |VKS|. Unsurprisingly, the virial theorem is more nearly satisfied when the contribution from Tc is included; this suggests that Tc should be included in QTAIM computations of atomic energies.

Table I.

Values of the virial ratio for both LDA and GGA calculations.

MoleculeLDAGGA-PBE
γγγγ
H2O 2.006 02 2.000 40 2.003 95 2.000 38 
NH3 2.007 26 2.000 21 2.004 74 2.000 20 
CH2O 2.006 33 2.000 36 2.004 10 2.000 36 
C6H6 2.007 68 2.000 27 2.005 07 2.000 34 
Fe(C5H5)2 2.002 71 2.000 03 2.001 81 2.000 05 
Cr(CO)6 2.003 10 2.000 11 2.002 01 2.000 13 
MoleculeLDAGGA-PBE
γγγγ
H2O 2.006 02 2.000 40 2.003 95 2.000 38 
NH3 2.007 26 2.000 21 2.004 74 2.000 20 
CH2O 2.006 33 2.000 36 2.004 10 2.000 36 
C6H6 2.007 68 2.000 27 2.005 07 2.000 34 
Fe(C5H5)2 2.002 71 2.000 03 2.001 81 2.000 05 
Cr(CO)6 2.003 10 2.000 11 2.002 01 2.000 13 

Because there is a small error in Eq. (9a), there is also a small error in Eq. (11). In order to overcome this problem, it is conventional to correct the atomic energies by scaling them with the so-called virial factor η,2,18–24

EΩ=(1γ)(TsΩ+TcΩ)η(TsΩ+TcΩ).
(13)

Tables II and III show our calculations for the QTAIM scaled and unscaled energies. To perform the QTAIM atomic integrations we used the fast grid-based method recently introduced;7–9 we ensured that every atomic integration was accurate to at least LΩ103a.u..7–9 In Tables II and III, we can see that energies obtained using the atomic virial theorem without the Tc correction,

EΩ=(1γ)TsΩηTsΩ,
(14)

can be considered a first approximation to the atomic energies, Eq. (13). The biggest differences between the “first approximation” (EΩ) and the improved result from Eq. (13)(EΩ) occur for hydrogen atoms; the atomic energies of hydrogen atoms sometimes differ by more than 100%. The importance of these differences will depend on the particular application where QTAIM is used. For example, recently Vila and Mosquera25 defined the strain energy of oxiranes based on the (noncorrected) QTAIM atomic energies obtained from KS-DFT calculations [Eq. (14)]. Their strain energy definition is in terms of relative differences between the QTAIM energies of the constituent atoms. Thus some of their conclusions might be modified when the correct KS-DFT QTAIM energies are used [Eq. (13)]. (See Ref. 26 for another interesting chemical application of QTAIM energies that might be affected by the conclusions of the present article.)

Table II.

Values of Ts, Tc, and the QTAIM atomic energies within the LDA approximation. The SCF energy is shown in the last column. Units are hartree.

MoleculeTsTc(Ts+Tc)E=Ts×ηE=(Ts+Tc)×ηSCF energy
H2O 
74.692 27 −0.304 23 −74.388 03 −75.141 89 −74.418 02   
0.379 90 0.363 91 −0.743 81 −0.382 19 −0.744 11   
Total 75.452 07 0.423 59 −75.875 66 −75.906 26 −75.906 25 −75.906 25 
  
NH3 
54.309 17 −0.844 52 −53.464 65 −54.703 21 −53.475 97   
0.464 06 0.412 26 −0.876 32 −0.467 43 −0.876 51   
Total 55.701 36 0.392 26 −56.093 62 −56.105 50 −56.105 49 −56.105 50 
  
CH2O 
36.919 22 −0.739 37 −36.179 85 −37.152 82 −36.192 72   
74.863 36 0.212 44 −75.075 80 −75.337 05 −75.102 52   
0.572 04 0.600 52 −1.172 56 −0.575 66 −1.172 97   
Total 112.926 66 0.674 09 −113.600 76 −113.641 19 −113.641 19 −113.641 19 
  
C6H6 
37.488 90 −1.052 78 −36.436 12 −37.776 83 −36.445 95   
0.582 43 1.334 83 −1.917 26 −0.586 90 −1.917 78   
Total 228.427 98 1.692 28 −230.120 26 −230.182 42 −230.182 36 −230.182 39 
  
Fe(C5H5)2 
Fe 1259.6627 −0.8370 −1258.8257 −1263.0793 −1258.8666   
37.5263 −0.6798 −36.8465 −37.6281 −36.8477   
0.5677 1.2031 −1.7708 −0.5692 −1.7709   
Total 1640.6031 4.3963 −1644.9994 −1645.0529 −1645.0528 −1645.0529 
  
Cr(CO)6 
Cr 1040.8165 −1.3822 −1039.4343 −1044.0428 −1039.5506   
36.8912 −2.5203 −34.3709 −37.0055 −34.3748   
75.0018 3.6031 −78.6050 −75.2343 −78.6137   
Total 1712.1744 5.1151 −1717.2895 −1717.4817 −1717.4816 −1717.4816 
MoleculeTsTc(Ts+Tc)E=Ts×ηE=(Ts+Tc)×ηSCF energy
H2O 
74.692 27 −0.304 23 −74.388 03 −75.141 89 −74.418 02   
0.379 90 0.363 91 −0.743 81 −0.382 19 −0.744 11   
Total 75.452 07 0.423 59 −75.875 66 −75.906 26 −75.906 25 −75.906 25 
  
NH3 
54.309 17 −0.844 52 −53.464 65 −54.703 21 −53.475 97   
0.464 06 0.412 26 −0.876 32 −0.467 43 −0.876 51   
Total 55.701 36 0.392 26 −56.093 62 −56.105 50 −56.105 49 −56.105 50 
  
CH2O 
36.919 22 −0.739 37 −36.179 85 −37.152 82 −36.192 72   
74.863 36 0.212 44 −75.075 80 −75.337 05 −75.102 52   
0.572 04 0.600 52 −1.172 56 −0.575 66 −1.172 97   
Total 112.926 66 0.674 09 −113.600 76 −113.641 19 −113.641 19 −113.641 19 
  
C6H6 
37.488 90 −1.052 78 −36.436 12 −37.776 83 −36.445 95   
0.582 43 1.334 83 −1.917 26 −0.586 90 −1.917 78   
Total 228.427 98 1.692 28 −230.120 26 −230.182 42 −230.182 36 −230.182 39 
  
Fe(C5H5)2 
Fe 1259.6627 −0.8370 −1258.8257 −1263.0793 −1258.8666   
37.5263 −0.6798 −36.8465 −37.6281 −36.8477   
0.5677 1.2031 −1.7708 −0.5692 −1.7709   
Total 1640.6031 4.3963 −1644.9994 −1645.0529 −1645.0528 −1645.0529 
  
Cr(CO)6 
Cr 1040.8165 −1.3822 −1039.4343 −1044.0428 −1039.5506   
36.8912 −2.5203 −34.3709 −37.0055 −34.3748   
75.0018 3.6031 −78.6050 −75.2343 −78.6137   
Total 1712.1744 5.1151 −1717.2895 −1717.4817 −1717.4816 −1717.4816 
Table III.

Values of Ts, Tc, and the QTAIM atomic energies within the GGA (PBE) approximation. The SCF energy is shown in the last column. Units are hartree.

MoleculeTsTc(Ts+Tc)E=Ts×ηE=(Ts+Tc)×ηSCF Energy
H2O 
75.272 39 −0.487 56 −74.784 83 −75.569 76 −74.813 36   
0.404 95 0.379 49 −0.784 44 −0.406 55 −0.784 74   
Total 76.082 30 0.271 41 −76.353 71 −76.382 87 −76.382 84 −76.382 85 
  
NH3 
54.764 44 −1.016 31 −53.748 14 −55.024 07 −53.758 99   
0.494 18 0.423 85 −0.918 03 −0.496 52 −0.918 21   
Total 56.246 98 0.255 23 −56.502 22 −56.513 63 −56.513 62 −56.513 63 
  
CH2O 
37.283 93 −0.875 97 −36.407 96 −37.436 82 −36.421 18   
75.452 95 0.069 99 −75.522 94 −75.762 35 −75.550 36   
0.608 99 0.615 87 −1.224 87 −0.611 49 −1.225 31   
Total 113.954 87 0.425 77 −114.380 63 −114.422 15 −114.422 16 −114.422 16 
  
C6H6 
37.859 49 −1.208 07 −36.651 42 −38.051 62 −36.663 97   
0.616 12 1.390 08 −2.006 19 −0.619 24 −2.006 88   
Total 230.853 62 1.092 05 −231.945 67 −232.025 17 −232.025 10 −232.025 14 
  
Fe(C5H5)2 
Fe 1262.2958 −1.1179 −1261.1779 −1264.5862 −1261.2471   
37.9025 −0.8611 −37.0415 −37.9713 −37.0435   
0.5952 1.2627 −1.8579 −0.5963 −1.8580   
Total 1647.2732 2.8982 −1650.1714 −1650.2621 −1650.2620 −1650.2621 
  
Cr(CO)6 
Cr 1043.2084 −1.6351 −1041.5733 −1045.3090 −1041.7058   
37.2918 −2.6980 −34.5937 −37.3669 −34.5981   
75.5902 3.5114 −79.1015 −75.7424 −79.1116   
Total 1720.5000 3.2451 −1723.7451 −1723.9643 −1723.9642 −1723.9643 
MoleculeTsTc(Ts+Tc)E=Ts×ηE=(Ts+Tc)×ηSCF Energy
H2O 
75.272 39 −0.487 56 −74.784 83 −75.569 76 −74.813 36   
0.404 95 0.379 49 −0.784 44 −0.406 55 −0.784 74   
Total 76.082 30 0.271 41 −76.353 71 −76.382 87 −76.382 84 −76.382 85 
  
NH3 
54.764 44 −1.016 31 −53.748 14 −55.024 07 −53.758 99   
0.494 18 0.423 85 −0.918 03 −0.496 52 −0.918 21   
Total 56.246 98 0.255 23 −56.502 22 −56.513 63 −56.513 62 −56.513 63 
  
CH2O 
37.283 93 −0.875 97 −36.407 96 −37.436 82 −36.421 18   
75.452 95 0.069 99 −75.522 94 −75.762 35 −75.550 36   
0.608 99 0.615 87 −1.224 87 −0.611 49 −1.225 31   
Total 113.954 87 0.425 77 −114.380 63 −114.422 15 −114.422 16 −114.422 16 
  
C6H6 
37.859 49 −1.208 07 −36.651 42 −38.051 62 −36.663 97   
0.616 12 1.390 08 −2.006 19 −0.619 24 −2.006 88   
Total 230.853 62 1.092 05 −231.945 67 −232.025 17 −232.025 10 −232.025 14 
  
Fe(C5H5)2 
Fe 1262.2958 −1.1179 −1261.1779 −1264.5862 −1261.2471   
37.9025 −0.8611 −37.0415 −37.9713 −37.0435   
0.5952 1.2627 −1.8579 −0.5963 −1.8580   
Total 1647.2732 2.8982 −1650.1714 −1650.2621 −1650.2620 −1650.2621 
  
Cr(CO)6 
Cr 1043.2084 −1.6351 −1041.5733 −1045.3090 −1041.7058   
37.2918 −2.6980 −34.5937 −37.3669 −34.5981   
75.5902 3.5114 −79.1015 −75.7424 −79.1116   
Total 1720.5000 3.2451 −1723.7451 −1723.9643 −1723.9642 −1723.9643 

The virial factor η(η) effectively corrects the “virial theorem error” from approximate KS calculations and allows the SCF molecular energy to be evaluated as the sum of the atomic energies. Notice that the atomic correlation-kinetic energy is sometimes quite negative, especially on heavier atoms. This is not an error: It is possible for the atomic Tc to be negative even though the global quantity will always be positive because of the enormous freedom one has in choosing a definition for the local kinetic energy.27–29 

In conclusion, in this letter we derived, for the first time, exact formulas for the QTAIM atomic energies within the KS-DFT formalism. Comparing results from the new formula, Eq. (13), to the conventional atomic energy formula for DFT-based QTAIM calculations, Eq. (14), reveals significant differences, especially for Hydrogen atoms. This suggests that it is important to include the contribution of the correlation-kinetic energy in KS-DFT-based calculations of QTAIM atomic energies.

J.I.R. gratefully acknowledges the Vrije Universiteit Amsterdam theoretical chemistry group and all SCM people for hosting him during the development of this work and for sharing their computer facilities. J.I.R. also thanks Miss Silvia Dimitrova for her insightful discussions and for THDs. J.I.R. and P.W.A. acknowledge research support from NSERC and insightful and encouraging discussions with Richard Bader. A.W.G. would like to thank Professor L. Visscher for generous support and is grateful for funding by NWO in the form of a VICI grant for L. Visscher.

1.
The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design
, edited by
C. F.
Matta
and
R. J.
Boyd
(
Wiley-VCH
,
Weinheim
,
2007
).
2.
R. F. W.
Bader
,
Atoms in Molecules: A Quantum Theory
(
Oxford University Press
,
New York
,
1990
).
3.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
4.
G.
te Velde
,
F. M.
Bickelhaupt
,
S. J. A.
van Gisbergen
,
C.
Fonseca Guerra
,
E. J.
Baerends
,
J. G.
Snijders
, and
T.
Ziegler
,
J. Comput. Chem.
22
,
931
(
2001
).
5.
C.
Fonseca Guerra
,
J. G.
Snijders
,
G.
te Velde
, and
E. J.
Baerends
,
Theor. Chem. Acc.
99
,
391
(
1998
).
6.
ADF2008, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, http://www.scm.com.
7.
J. I.
Rodríguez
, Ph.D. dissertation,
McMaster University
,
2008
.
8.
J. I.
Rodríguez
,
A. M.
Köster
,
P. W.
Ayers
,
A.
Santos-Valle
,
A.
Vela
, and
G.
Merino
,
J. Comput. Chem.
30
,
1082
(
2009
).
9.
J. I.
Rodríguez
,
R. F. W.
Bader
,
P. W.
Ayers
,
C.
Michel
,
A. W.
Götz
, and
C.
Bo
,
Chem. Phys. Lett.
472
,
149
(
2009
).
10.
S. K.
Ghosh
and
R. G.
Parr
,
J. Chem. Phys.
82
,
3307
(
1985
).
11.
M.
Levy
and
J. P.
Perdew
,
Phys. Rev. A
32
,
2010
(
1985
).
12.
W.
Yang
and
R. G.
Parr
,
Density Functional Theory of Atoms and Molecules
(
Oxford University Press
,
New York
,
1989
).
13.
S. K.
Ghosh
and
V. A.
Singh
,
J. Phys.: Condens. Matter
1
,
1971
(
1989
).
14.
H.
Ou-Yang
and
M.
Levy
,
Phys. Rev. Lett.
65
,
1036
(
1990
).
15.
P. A. M.
Dirac
,
Proc. Cambridge Philos. Soc.
26
,
376
(
1930
).
16.
S. H.
Vosko
,
L.
Wilk
, and
M.
Nusair
,
Can. J. Phys.
58
,
1200
(
1980
).
17.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
18.
C. F.
Matta
and
R. J.
Boyd
, in
The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design
, edited by
C. F.
Matta
and
R. J.
Boyd
(
Wiley-VCH
,
Weinheim
,
2007
).
19.
T. A.
Keith
, AIMALL (Dev. Version),
2008
, aim.tkgristmill.com.
20.
T.
Keith
, in
The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design
, edited by
C. F.
Matta
and
R. J.
Boyd
(
Wiley-VCH
,
Weinheim
,
2007
).
21.
MORPHY, a program for an automatic atoms in molecules analysis;
P. L. A.
Popelier
,
Comput. Phys. Commun.
93
,
212
(
1996
).
22.
C.
Gatti
, TOPOND-98: An Electron Density Topological Program for Systems Periodic in N(N=03) Dimensions, User’s Manual (
CNR_ISTM
,
Milano
,
1999
).
23.
C.
Gatti
,
V. R.
Saunders
, and
C.
Roetti
,
J. Chem. Phys.
101
,
10686
(
1994
).
24.
C.
Gatti
and
F.
Cargnoni
,
Proceedings III Convenzo Nazionale di Informatica Chimica
, Napoli, Italy,
1997
(unpublished).
25.
A.
Vila
and
R. A.
Mosquera
,
Chem. Phys.
287
,
125
(
2003
).
26.
P. M.
Quiñonez
,
A.
Vila
,
A. M.
Graña
, and
R. A.
Mosquera
,
Chem. Phys.
287
,
227
(
2003
).
27.
L.
Cohen
,
J. Chem. Phys.
80
,
4277
(
1984
).
28.
L.
Cohen
,
J. Chem. Phys.
70
,
788
(
1979
).
29.
P. W.
Ayers
,
R. G.
Parr
, and
A.
Nagy
,
Int. J. Quantum Chem.
90
,
309
(
2002
).