Chemistry With ADF

Chemistry with ADF
G. TE VELDE,1 F. M. BICKELHAUPT,2 E. J. BAERENDS,2

C. FONSECA GUERRA,2 S. J. A. VAN GISBERGEN,3
J. G. SNIJDERS,4 T. ZIEGLER5
1
Paragon Decision Technology BV, Julianastracit 30, Postbus 3277, NL-2001 DC Haarlem,
The Netherlands
2
Theoretical Chemistry, Vrije Universiteit, De Boelelaan 1083, NL-1081 HV Amsterdam,
The Netherlands
3
Scientific Computing & Modelling NV, Theoretical Chemistry, Vrije Universiteit,
De Boelelaan 1083, NL-1081 HV Amsterdam, The Netherlands
4
Theoretical Chemistry, Materials Science Centre, Rijksuniversitet Groningen, Nijenborgh 4,
NL-9747 AG Groningen, The Netherlands
5
Department of Chemistry, University of Calgary, 2500 University Drive NW, Calgary, Alberta,
Canada T2N 1N4
Received 16 June 2000; accepted 4 December 2000
ABSTRACT: We present the theoretical and technical foundations of the

Amsterdam Density Functional (ADF) program with a survey of the
characteristics of the code (numerical integration, density fitting for the Coulomb
potential, and STO basis functions). Recent developments enhance the efficiency
of ADF (e.g., parallelization, near order-N scaling, QM/MM) and its functionality
(e.g., NMR chemical shifts, COSMO solvent effects, ZORA relativistic method,
excitation energies, frequency-dependent (hyper)polarizabilities, atomic VDD
charges). In the Applications section we discuss the physical model of the
electronic structure and the chemical bond, i.e., the KohnSham molecular orbital
(MO) theory, and illustrate the power of the KohnSham MO model in
conjunction with the ADF-typical fragment approach to quantitatively
understand and predict chemical phenomena. We review the Activation-strain
TS interaction (ATS) model of chemical reactivity as a conceptual framework for
understanding how activation barriers of various types of (competing) reaction
mechanisms arise and how they may be controlled, for example, in organic
chemistry or homogeneous catalysis. Finally, we include a brief discussion of
exemplary applications in the field of biochemistry (structure and bonding of
DNA) and of time-dependent density functional theory (TDDFT) to indicate how
this development further reinforces the ADF tools for the analysis of chemical
phenomena. c 2001 John Wiley & Sons, Inc. J Comput Chem 22: 931967,
2001
Correspondence to: F. M. Bickelhaupt; e-mail: bickel@chem.

vu.nl; G. te Velde; e-mail: Bert.te.Velde@Paragon.nl
Contract/grant sponsors (to F.M.B): Deutsche Forschungs-
gemeinschaft, and Fonds der Chemischen Industrie
Journal of Computational Chemistry, Vol. 22, No. 9, 931967 (2001)

c 2001 John Wiley & Sons, Inc.
TE VELDE ET AL.
Keywords: ADF program; density functional theory; materials science; chemical

bond; reactivity
The potential VXC is the functional derivative

Introduction with respect to the density of EXC [], the
exchangeandcorrelation energy functional. The
THE SUCCESS OF DENSITY FUNCTIONAL one-electron molecular orbitals (MOs) i with cor-
THEORY (DFT) responding orbital energies i define the exact elec-
tronic charge density and give, in principle, ac-
K ohnSham DFT1 5 has in the recent past

emerged as an important first-principles com-
putational method6 15 to predict chemical prop-
cess to all properties because these are express-
ible as functional of the density, in particular the
energy. Moreover, they provide an intuitively ap-
erties accurately16 22 and to analyze and interpret pealing view of the system as being built from
these in convenient and simple chemical terms.23 31 independent-electron orbitals with all ensuing inter-
The reasons for its popularity and success are pretations.
easy to understand. In the first place, the DFT ap- The exact form of the exact energy density EXC (r),
proach is in principle exact. The exact exchange- representing and incorporating all exchange and
and-correlation (XC) functional is unknown, but correlation (XC) effects is unknown. From general
the currently available XC functionals provide in principles one can formulate conditions on what
most cases already a chemical accuracy of a few EXC (r) should look like, and several, more and
kcal/mol for binding energies.16, 22, 32, 33 Moreover, more advanced expressions have been advocated
the quest for more accurate ones based on a more for it in the literature. Their application to real
detailed understanding of their essential properties systems has been impressively successful, and it
is continuing.27, 28, 34 46 seems likely that we are on the right track, and
In the second place, it preserves at all lev- that a further increase of accuracy is a matter of
els of approximation the appealing one-electron time.
molecular orbital (MO) view on chemical reactions
and properties. The computed orbitals are suit- THE ADF PROGRAM
able for the typical MO-theoretical analyses and
interpretations.27, 29 31 The KS method effectively The Amsterdam Density Functional (ADF)
incorporates all correlation effects. program has been developed since the early
In the third place, it is a relatively efficient com- seventies6 at that time named HFS, for Hartree
putational method, and its fundamental scaling FockSlaterwith the explicit purpose to exploit
properties do not deteriorate when methodological the DFT computational advantages. With later
precision is increased, in particular, when a more (and still ongoing) additions and improvements,
accurate XC functional is applied. Recent research notably by the Theoretical Chemistry research
paves the way to implementations that scale only groups of Amsterdam6, 11 15, 62 73 and Calgary74 85
linearly with the system size.15, 47 60 This brings ADF has evolved into a state-of-the-art package for
within reach the treatment by fundamental quan- quantum-chemical research. This article is based
tum chemical methods of systems with hundreds, on the ADF2000/2001 releases. Documentation,
maybe even thousands of atoms. including user manuals and licensing and distri-
bution information, is available at the Web site
THE KOHNSHAM MO MODEL http://www.scm.com.
ADF supports a wide variety of exchange
The basic postulate in KohnSham DFT61 is that correlation (XC) functionals (but no exact-exchange
we can apply a one-electron formulation to the varieties that require the evaluation of multicen-
system of N interacting electrons by introducing a ter integrals; see the ExchangeCorrelation sec-
suitable local potential VXC (r), in addition to any ex- tion for details), and incorporates relativistic effects
ternal potentials Vext (r) and the Coulomb potential with the ZORA formalism, in the scalar approach
of the electron cloud VC (r), and solving or with spin-orbit terms included. Basis sets are
built with Slater-type orbital functions. The stan-
12 2 + Vext (r) + VC (r) + VXC (r) i (r) = i i (r) (1) dard sets in the ADF database range from minimal
932 VOL. 22, NO. 9

CHEMISTRY WITH ADF
to triple-zeta-plus-double-polarization quality (and

beyond, for some atoms). An easy-to-use utility Part I. Implementation:
program is provided to generate higher precision Technical Characteristics
basis sets as required. Energy gradients and sec-
NUMERICAL INTEGRATION
ond derivatives (computed analytically) allow the
computation of energy minima, transition states, All DFT implementations employ some kind of
reaction paths, and harmonic frequencies with IR numerical integration scheme, because the mathe-
intensities. Solvation and electric field environment matical expression for any relevant density func-
effects can be taken into account, and the com- tional makes an analytical integration of matrix el-
putable molecular properties include NMR chem- ements of the XC potential between basis functions
ical shifts, ESR, and various response proper- of whatever type impossible. ADF uses a Gaussian-
ties using the time-dependent density functional type quadrature method11, 12, 89 based on the parti-
theory (TDDFT): excitation energies, frequency- tioning of space in atomic cells. Each cell contains
dependent (hyper)polarizablities, RAMAN intensi- a core-region sphere, inside which the integrands
ties, and dispersion coefficients. Numerical preci- with cusps and singularities at the nucleus are han-
sion is controlled by a flexible numerical integration dled conveniently in spherical coordinates:
Z Z
scheme for the evaluation of integrals. QM/MM
f (r) dr = f (r) dr
and linear scaling techniques are available for the cell sphere
treatment of very large systems, while paralleliza- Z
tion speeds up the execution on multi-CPU sys- + f (r) dr (2)
cell-without-sphere
tems.
Z Z R Z 4
The first part of this articlesecond and third
f (r) dr = r2 dr f (r, ) d (3)
sectionsreviews the ADF implementation, regard- sphere 0 0
ing technical characteristics and functionality. The
package is large, and the number of features steadily The Jacobian factor r2 in eq. (3) removes the
Coulomb singularity of the integrand at the nuclear
expanding. The scope and size of this article limit
position. The remaining cell-minus-sphere region is
the amount of detail that can be shown in these sec-
treated, after suitable coordinate transformations,
tions, and some aspects will only be touched upon, by a Gaussian product formula in Cartesian-like co-
if at all. In particular, the article does not include ordinates.
a discussion of the solid-state counterpart of ADF The atomic cells fill all space inside the mole-
(BAND).12, 86 88 The second part of the article is de- cule. In the remaining outer region, farther away
voted to main themes in ADF from the outset, i.e., from the atoms, all integrands decay exponentially
the analysis of chemical bonds and reactions and with the distance from the molecule are smooth, and
the development of tools and concepts that help do not have cusps. An accurate evaluation by nu-
achieve a deeper and more detailed understanding merical sampling techniques over that part of space
of electronic structures and reaction mechanisms. is not very demanding.
We elaborate upon these subjects in sections four The ADF numerical integration module automat-
six on the chemical bond, chemical reactivity, and ically tunes the grid by varying the numbers of
an application in the field of theoretical biochem- points in the different regions (atomic spheres, cells
istry (structure and bonding of DNA), while the and the outer region), while monitoring a series
seventh section contains an exemplary application of test integrals such that these are evaluated with
an (input-adjustable) precision. This procedure has
of TDDFT.
proven to be robust and reliable, and is able to pro-
To date, there has been no scientific publication
duce any desired accuracy in the actual calculations.
that covers both the various methods and tech-
It is user-friendly, because one only has to specify
niques in ADF (and the reason for particular choices the precision level, by a single parameter. At the
made) and how they are applied, for all the many same time it is flexible, as one can select any value
features in ADF together. We have found that many from 0.5 (so inaccurate as to be almost meaningless)
users of the ADF software see this as a problem. to 12.0 (effectively close to machine precision). Of
By integrating methodology and chemical appli- course, a higher precision implies more points and a
cations in the present article, we hope to fill the greater computational effort. As an indication of the
gap. scaling of computing time: increasing the precision
JOURNAL OF COMPUTATIONAL CHEMISTRY 933

TE VELDE ET AL.
parameter by 1.0 will roughly double the number sets for relativistic calculations within the ZORA
of grid points, and hence, the computing time. This formalism.65, 66, 68, 69, 90 Easy-to-use tools are avail-
is a crude estimate; much depends on the actual able to generate basis sets with any required accu-
case. racy, allowing (in principle) to approach the basis set
ADF applies numerical integration not only to limit and eliminate any inaccuracies resulting from
those integrals that cannot be done analytically, basis set deficiencies. The implemented method is
such as the matrix elements of the XC potential, but based upon even-tempered STO basis sets devel-
to most other integrals as well. This approach has oped by Chong, similar to earlier work,91 but with
two major advantages. In the first place, simplic- a few modifications that, for molecular properties,
ity, and hence efficiency, as well as reduced chances seem to give somewhat better results.
of implementation errors. Numerical integrals, as-
suming that the grid itself is available, are easy for FROZEN CORE APPROXIMATION
the developer to implement and for the computer to
evaluate: straightforward do-loops and vector con- An important aspect of ADF is the (optional)
structs. In the second place, it creates the freedom use of a frozen core to trim down the computa-
to select the types of basis functions that are in- tional time by reducing the size of the variational
trinsically the most suitable for electronic structure basis set. Because deep-core atomic orbitals change
calculations, but which might cause difficulties in very little upon bond formation, the frozen core ap-
analytical integration schemes, for instance Slater- proximation saves time without sacrificing much
type orbitals (STOs) or numerical atomic orbitals. It in quality. Use of frozen cores is preferred over
makes it almost trivial to implement the computa- pseudopotentialswhich ADF does not support
tion of all kinds of properties, as long as these are because it is, in essence, an all-electron calculation.
expressible as integrals over space of an integrand With a frozen core calculation you obtain the total
that one can evaluate without too much difficulty at charge density and potential in the valence and in
an arbitrary point r. the core region, ignoring only the slight change in
the deep-core orbitals upon formation of the chemi-
BASIS SETS cal bond.
To ensure orthogonality between the frozen core
ADF exploits Slater-type functions (STOs) f (r) = orbitals and the valence basis functions, the valence
Ylm ()rn er as basis set elements, in terms of which basis set used in the SCF equations is explicitly or-
the molecular orbitals (MOs) are expanded. The thogonalized to the frozen core. This is achieved by
center of the function is at a nucleus, the Ylm are adding to the valence set a series of auxiliary core
spherical harmonics, and the exponential factor functions core , precisely one for each frozen core
determines the long-range decay of the function. orbital. Each valence function is replaced by a linear
In contrast to Gaussian-type orbitals (GTOs), used combination
in many other quantum-chemical codes, STOs have X
the advantage of possessing the required cusp be- valence valence + c core (4)

havior as well as the appropriate long-range decay.
This allows the construction of high-quality ba- The condition that each such modified valence func-
sis sets with relatively small numbers of functions. tion be orthogonal to each frozen core orbital deter-
GTOs, in comparison, typically would need a factor mines exactly all the coefficients c . The frozen core
of 3 more functions for the same level of basis set orbitals are taken from very accurate single-atom
quality. (In practice, GTO schemes usually employ calculations with large STO basis sets. Note that the
contracted basis sets, which, of course, reduce again core functions do not add variational freedom to the
the computational effort that the uncontracted set (valence) basis set, and hence, do not increase the
would require.) size of the secular (Fock) matrix to be diagonalized.
The ADF package comprises a database with The ADF database offers for most atoms in the
a large number of basis sets for all elements of periodic table several basis sets with different lev-
the periodic table, ranging from minimum single- els of frozen core approximation, including all-
(SZ) to doubly polarised triple- (TZ2P) basis electron basis sets to avoid the frozen core approxi-
sets (http://www.scm.com). It incorporates spe- mation altogether. The exact definition of the frozen
cial extended sets with very diffuse functions core orbitals (expansion functions and coefficients)
for (hyper)polarizability calculations,67, 70 73 and used with a basis set is stored in the corresponding
others with very contracted deep-core function database file.
934 VOL. 22, NO. 9

CHEMISTRY WITH ADF
SYMMETRY VXC (r) and any other (external) potentials and terms
of the Fock operator and evaluate the Fock matrix
ADF exploits molecular symmetry for analysis elements by a single numerical integration, which
purposes and to enhance the computational effi- we have to carry out anyway, namely for VXC (r).
ciency. The latter applies in particular to numeri- X
1
cal integrals, where only the integration over the F wk (rk ) 2 + VCoulomb (rk )
symmetry-unique region in space has to be carried 2
k
out (assuming that the integrands have been pro-
jected onto their totally symmetric component). The + VXC (rk ) + (rk ) (5)
number of integration points is thereby reduced
by a factor equal to the number of symmetry op- The summation runs over all integration points rk .
erators. Almost all of the commonly encountered The only effort, apart from the trivial addition of
point group symmetries that have real irreducible VC (r) to VXC (r), lies therefore in computing the value
representations (irreps) have been implemented, in- of VC (rk ).
cluding nonabelian point groups with up to an
Density Fitting
eight-fold symmetry axis. The program recognizes
the symmetry automatically, but the user may spec- The density fit implies that the exact density,
ify to use a lower symmetry via input. Preparations expressed as one- and two-center products of ba-
have been made to extend the symmetry function- sis functions from the computed occupied MOs, is
ality to include all point group symmetries. replaced by an expansion in one-center functions
The KohnShamFock matrix is constructed on only.6 The auxiliary fit functions are (also) STO-
the basis of functions that transform according type functions, but different from the basis func-
to a particular irrep and subspecies of the point tions. In fact, for a single-atom calculation where all
group. The symmetry-adapted functions are com- products of basis functions are one-center products,
binations of the original functions, orthogonal- the perfect set of fit functions would include all prod-
ized on the frozen core as described earlier, and ucts of basis functions. These fit functions are then
Lwdin orthonormalized among themselves, result- not only different from the basis functions, but there
ing in an orthonormal and symmetry-adapted basis are also more of them.
set in which the Fock matrix is block-diagonal. In molecules, the exact representation of two-
Diagonalizationblock by blockyields the MOs center density terms by one-center expansions is
grouped by symmetry representation. practically impossible. A very accurate approxima-
tion is obtained, however, by taking a sufficiently
COULOMB POTENTIAL VIA CHARGE
large set of fit functions, and in particular, by adding
DENSITY FITTING
polarization fit functions, i.e., STO-type fit func-
The analytical calculation of Coulomb matrix el- tions with higher angular momentum quantum
ements would require two-electron integrals with numbers. Practice has shown that this is an ade-
its characteristic n4 scaling, leaving distance cutoffs quate and accurate approach. Moreover, the ADF
aside. The involved four-center integrals, in the case utility program to generate basis sets with any re-
of the physically appropriate STO basis functions, quired completeness accuracy can equally be ap-
are far from trivial. Numerical integration, in combi- plied to generate high-precision fit sets in case the
nation with density fitting, greatly reduces the effort. existing fit sets may be found to be inadequate.
The evaluation is simpler and the formal scaling is The density fit in ADF is performed on an atom-
only n3 rather than n4 . pair basis: the total density is partitioned in its dis-
tinct two-center (and one-center) components, and
Numerical Integration each atom-pair-density is fitted by fit functions on
these two atoms only. As a result, the computational
With the exception of the Coulomb potential, all weight of the fitting procedure scales principally
terms in the KohnSham operator of eq. (1) are lo- only as n2 instead of n3 (linear scaling techniques
cal in the sense that they can be evaluated readily further reduce this). The possible benefit from using
as functions of r. Rather than separately evaluating off-pair fit functions is heavily system dependent,
all Coulomb matrix elements, we first compute (see and the fit sets would intrinsically be less robust and
below) the value of the Coulomb potential VC (r) in accurate if we relied on this. If we want to make
each point r of the numerical integration grid. We sure that the fit quality is sustained independent
then add VC (r) to the exchangecorrelation potential of the actual system, the pair-wise fitting has to be

TE VELDE ET AL.
adequate in itself, and the use of off-pair functions new trial density matrix Pk+1 , etc., until Pk and Pk+1
becomes merely a computational burden. are identical, and we have reached self-consistency.
The fit equation for the expansion coefficients ci Simply plugging the new Pk+1 into the iterative
is procedure will often fail to converge and result in
X oscillatory behavior. Various methods exist to make
exact (r) ci fifit (r) (6)
the SCF converge, the simplest one being simple
i
damping, in which Pk and Pk+1 are averaged to deter-
This is solved in a least-squared sense under the mine a new P as trial density matrix. ADF employs
constraint that the number of electrons in the fit the DIIS (direct inversion of the iterative subspace)
representation is preserved. Having obtained the method to obtain convergence.92 95 This is a gen-
approximate fit density approximation [eq. (6)], the eralization of simple damping, taking the results of
corresponding Coulomb potential is several previous cycles to construct a guess for the
X next cycle. With DIIS, it typically takes between 5
VCfit (r) = ci fic (r) (7)
and 30 cycles to converge the SCF equations.
i
For an STO-type fit function f (r) = Ylm ()rn er the LARGE MOLECULES
corresponding Coulomb potential function f c (r) is
given by DFT methods have, in principle, the same scal-
Z ing behavior as HartreeFock (HF), namely n4 when
f (r0 )
f c (r) dr0 using two-electron integrals for the Coulomb poten-
|r r0 |
Z tial matrix elements. Most implementations, how-
0 X
l
4 0 r< ever, scale formally as n3 as a result of using either
= dr0 Ylm ()(r0 )n er 9
0 0 ( )
0 0
2l0 + 1 l m rl+1
> numerical integration or density fitting (or both)
lm
Z
for handling the Coulomb potential. As a conse-
4 rl< 0
= 9lm () (r0 )2+n er dr0 quence, they are intrinsically much more efficient
2l + 1 0 rl+1
> than correlated ab initio formalisms of higher-than-
(R r HF methodological accuracy (MP2, for instance).
0 n+l+2 r0
4 0 (r ) e dr0
= 9lm () Nevertheless, even a (formal) n3 scaling behavior
2l + 1 rl+1
makes the treatment of larger molecules a formi-
Z ) dable task. Two approaches have emerged to deal

0 nl+1 r0
+ rl (r ) e dr0 (8) with this problem, and most present-day codes ex-
0 ploit these to at least some extent, and display
At large distances r a long-range multipole term better scaling behavior as a result. In the first place,
dominates, which decays as q/rl+1 , where l is the there has been an increasing effort, still continuing,
angular momentum quantum number of the STO to develop intrinsically order-n algorithms (linear
fit function, and q is the strength. For the l = 0 scaling with system size). The idea is to exploit
monopole function, q would be the total charge distance effects based on the notion that the inter-
contained in the fit function, i.e., the integral over actions of a particular atom with other atoms in
all space of the fit function itself. For small r, the the molecule can be ignored, without significant
multipole term is modified by an exponentially de- loss of accuracy, for all atoms except those in its
caying small-range term; at r = 0, the function is neighborhood. If we can somehow carry this con-
zero. sideration over to all aspects of the calculation we
The database files in the ADF package contain should arrive at a method where the effort per atom
for each atom not only the basis functions (and the is limited by the maximum size of environment it
frozen core information) but also the fit functions chemically/physically sees. This is intrinsically
centered on that atom. not system size dependent, so the total effort should
scale linearly with the number of atoms.
SELF-CONSISTENT FIELD (SCF) PROCEDURE In the second place, one may realize that in
calculations on large systems we are usually only
The self-consistent solution of the KohnSham interested in a small active site where, say, the
equations is obtained in an iterative SCF procedure. chemical bond is formed. The remaining part of
From a trial density matrix Pk we compute the Fock the system must be included in the calculation, be-
matrix Fk and diagonalize it to obtain the eigenvec- cause leaving it out would remove relevant bound-
tors, i.e., the MOs. The occupied MOs then yield a ary conditionstypically: steric and electrostatic
936 VOL. 22, NO. 9

CHEMISTRY WITH ADF
effectsand thereby affect the outcomes too much. tween (basis, fit) functions are only computed when
The accurate treatment of those outer parts, how- the functions ranges overlap, while the actual eval-
ever, has often not much impact on the computed uation of these integrals only addresses those blocks
strength of the chemical bond, the energy profile of of points that fall in the range of both functions si-
the reaction, etc. These outer parts of the system can, multaneously. Many of the dominating aspects of a
hence, be treated by much faster methods of lower calculation are in this way reduced to algorithms
methodological accuracy (molecular mechanics), in that scale linearly. Of course, the system must be
combination (QM/MM) with the accurate but more large enough to actually see the linear scaling be-
demanding quantum-mechanical methods for the havior, and for some aspects the onset of the linear
crucial centre of the system. regime may become apparent only for fairly large
molecules.
Linear Scaling Analyzing all the steps in a calculation we find
a large number of them that scale worse than lin-
The basic elements in a calculation that deter- early, and we would have to treat all of them to
mine the computational effort and that increase formally claim an order-n program. From a practical
with system size are (1) the numbers of basis and point of view, however, we need to deal only with
fit functions, and (2) the number of numerical inte- the dominant aspects in calculations on systems that
gration grid points. Both aspects scale linearly with are of interest now and in the very near future,
the system size. The number of fit functions is larger lets say up to a few hundreds, or maybe a thou-
than the number of basis functions; up to a factor sand atoms. The ADF development projectstill
of 2 is quite common. The number of grid points is ongoingto establish linear scaling has focussed
also proportional to the size of the system (number first on the SCF procedure, which is concentrated
of atoms), but with an even much larger proportion- around the construction and processing of the Fock
ality factor. For the analysis of scaling behavior such matrix, and the TDDFT modules where the compu-
differences are irrelevant. We will, therefore, refer to tational structure has many similar features. Other
these items only by the scaling properties, and de- CPU-intensive steps, in the preparation before the
note each of the scaling numbers as n. SCF procedure can start, and in the postprocessing
To obtain better-scaling algorithms, ADF applies and computation of properties, have structural ele-
locality properties15 based on the common chemical ments that are also very similar, and that will, hence,
reasoning that interatomic interactions over large be amenable to the same treatment developed for
distances in a molecule can be ignored, which is, the Fock matrix.
in a sense, already reflected in the use of expo- Within the SCF procedure two relatively impor-
nentially decaying atomcentered basis functions. To tant aspects have not been treated completely yet.
each basis and fit function one assigns a range, be- First, the diagonalization of the Fock matrix scales
yond which it is assumed negligible. The grid points as n3 in the standard approach. At the moment,
are grouped together in small subsets (blocks of we have not yet tackled this aspect because the
points), typically a few tens to a few hundred, that total time spent in diagonalization has not been
all lie in the same region of space. Each block is then dominant so far in ADF calculations. Second, the
associated with a certain location (center and ex- evaluation in the grid points of the multipole long-
tension) in space, and the functions range and the range components of the Coulomb potential is still
location of the block determine whether the func- an order-n2 step because the 1/r type decaynot
tion needs to be evaluated in that particular block. exponentialprohibits a short-range cutoff. Fortu-
Application of a yes-or-no cutoff, as an alterna- nately, processing of the multipole terms is not very
tive to standard screening techniques, entails the demanding, but with increasing sizes of molecules
risk of discontinuities in the energy surface. How- this will sooner or later become significant, and this
ever, with stringent criteria this effect appears (so aspect is, therefore, a prime target for future efforts.
far) to be minor. The advantage of the cutoff ap- Obviously, the structure of the multipole potential is
proach is simplicity and a somewhat better perfor- suitable for tree-based fast multipole methods,96 as
mance improvement: the relevant data structure all components are based on the fit functions expan-
the lists of atoms that are in each others sphere sions, which can be organized and grouped together
of influencecan be computed in advance, and ap- per atom and per (l, m)-value. In ADF, we only need
plied wherever needed. to evaluate the Coulomb potential in the integra-
Carrying these considerations through has the tion points, instead of calculating Coulomb matrix
effect that (overlap, Fock, etc.) matrix elements be- elements directly. This will make the use of hierar-

TE VELDE ET AL.
FIGURE 1. Scaling behavior of the fit procedure, that

FIGURE 3. Scaling behavior of the buildup of the Fock
is, the calculation of fit coefficients, with and without
matrix, with and without distance effects.
distance effects.
chical multipole methods even simpler than when be just an effective average of underlying terms
they have to be applied to the direct calculation of with different fundamental scaling behavior. The fit
Coulomb matrix elements. procedure (Fig. 1) does not show perfect linear scal-
Figures 1, 2, and 3 display timing resultsin the ing because of memory-cache and I/O effects that
original and new implementationfor a systematic will have to be further analyzed. The result for the
sequence of semi-one-dimensional alkanes, where (Coulomb) potential in Figure 2, a scaling order of
the carbon atoms lie in a plane to form a regular 1.5 instead of 1.0, is caused by the long-range mul-
seesaw chain and the hydrogen atoms are below tipolar term, which has not yet been treated, and
and above the plane. In these systems the scaling hence, implies a term with scaling order 2. The Fock
improvements are exhibited most clearly due to matrix setupnot including the computation of the
their geometrical form where relatively many atoms potentialsscales perfectly linear.
are far apart. The scaling orders have been deter- With molecules of a few tens of atoms significant
mined by fitting the results over the data range, speed-ups are obtained by the implemented scaling
so they represent only the de facto scaling over the techniques. The truly linear regime is seen in sys-
range of used sizes, and will, generally speaking, tems with hundred atoms and more.
QM/MM
The quantum mechanicsmolecular mechanics

(QM/MM) approach combines a normal (DFT)
calculation on certain parts of the molecule or sys-
tem with an approximate, force field-based method
for the rest.84, 85, 97 Typically, this is applied when
one can partition a large system into a region
that requires an accurate quantum mechanical
descriptionthe active siteand a remainder that
acts only as a perturbation and can, therefore, be
dealt with by approximate methods. The calcula-
tion proceeds in an iterative fashion. A pure-MM
module optimizes the MM subsystem, taking the
QM parts into account in the given QM geome-
try and considering the QM atoms as (frozen) MM
FIGURE 2. Scaling behavior of the evaluation of the atoms without information on the wave functions.
(Coulomb and XC) potential in the numerical integration Next, a QM calculation on the QM subsystem is
grid, with and without distance effects. performed including the (steric and electrostatic) ef-
938 VOL. 22, NO. 9

CHEMISTRY WITH ADF
fects from the MM part in the SCF procedure and PARALLELIZATION

the MM forces on the QM atoms in the geometry
optimization step. With the new QM geometry the The implementation in ADF of parallelization13 is
MM module takes control again, until convergence. based on the distributed-memory concept, using the
single-program-multiple-data (SPMD) model. Al-
The MM module uses a force field that is speci-
though this may ignore some of the advantages of
fied on a separate data file. Currently available in
shared-memory architectures, it has the advantage
the ADF QM/MM procedure are AMBER9598 and
that it can be applied generally. Shared-memory
UFF99 force fields. The setup is very modular, and
machines are treated as being of the distributed-
other force fields can easily be inserted as additional memory type. For the parallel communication, PVM
options. and MPI are both supported in ADF. Where we have
When the boundary between the QM and MM to be specific in the following, PVM will be assumed
parts cuts one or more covalent bonds and the two for the sake of being concrete.
subsystems are not well separated, the QM sub- ADF is started on one of the nodes of the (vir-
system is defined with so-called capping atoms. tual) parallel machine; this process is called the
These are usually hydrogen atoms, and represent parent. The parent communicates with the system
the MM part in the QM cycles of the computation. to create the child processes; it then reads the in-
The real MM-atom that is bonded to the QM atom put and broadcasts this to the children. So, after
(and, hence, replaced by the capping atom in the the startup phase, a copy of ADF runs at each of
QM calculations) is called the link atom. the processes with exactly the same input data.
Various schemes have been advocated in the lit- In the serial (nonparallelized) parts of the program
erature to combine the QM and MM calculations. all these copies perform the same calculation, du-
They differ in their prescription for how the true plicating each others work. In the parallel parts,
bond between the QM-atom and the MM link atom each copy handles its share of the task and combines
is related to the pseudobond between the QM atom the results with the others. After all copies have
and the capping (hydrogen) atom. Most important finished the whole calculation, the child processes
approaches assume either or not a fixed direction terminate, their private files are destroyed, and we
are left with the results (output and other files) from
from the QM system to the capping atom, and use
the parent. A major advantage of the SPMD method
either a fixed distance to the capping atom or a fixed
is that the structure of the parallel program is virtu-
ratio between the distance to the capping atom and
ally identical to the structure of the serial program.
the distance to the true MM atom. The ADF pack- This makes it easier to maintain and develop both
age supports three schemes: (1) the one originally versions simultaneously. In fact, technically speak-
proposed by Singh and Kollman;100 (2) the approach ing, there are no separate parallel and serial versions
of Maseras and Morokuma;101 (3) a variation on at all.
the latter, using Cartesian rather than internal (Z-
matrix) coordinates.84 General Approach
In comparison with pure-QM calculations, even
when these are highly optimized by linear-scaling Communication between different processes is
techniques, the QM/MM approach is much more ef- relatively time-consuming, and easily becomes a
ficient. The MM contribution to the QM calculation bottleneck in massively parallel runs. Minimization
of communication is, therefore, crucial, even if at
cycles implies a small extra effort and the MM cal-
the expense of a slightly reduced balance of the
culation cycles are very fastdepending somewhat
workload distribution. Coarse-grained static load
on the chosen options. This is paid for by a reduced
balancing keeps the communication times small: the
accuracy, or more appropriately formulated, by an
data is kept local, and is not redistributed (static bal-
additional uncertainty in the results. One needs to ancing) while communication is infrequent (coarse
be sure that the partitioning has been adequate, and grain).
that not too much has been put into the MM approx- A general library has been developed for par-
imation. Until sufficient expertise has been built up, allelization operations. All calling sequences to
this may imply that one has to perform additional the procedures in the library are general, i.e., not
validation calculations. Nevertheless, the QM/MM specifically PVM or MPI or otherwise. The proce-
development clearly has paved the way to the treat- dure interfaces are wrappers that contain PVM and
ment of much larger systems than were considered MPI versions inside; the appropriate version is se-
doable in a not too distant past.102 105 lected at compile time, when the human installer

TE VELDE ET AL.
chooses the message-passing method. This ensures

that whatever is specific for the message-passing
method is isolated in the library, and when new
versions of PVM or MPI, or even completely new
parallel facilities are to be used, changes have to be
made only in this library.
Parallelization in ADF
The numerical integrations (matrix elements of

the Fock operator, in particular) consume a major
part of the CPU time in ADF. The number of in-
tegration points is typically of the order of a few
thousands per atom. The points are partitioned in
blocks of points, each block comprising a few tens
up to a few hundred points. The blocks are distrib- FIGURE 4. Parallel performance of ADF. Displayed
are the actual performance (obtained on an IBM SP2),
uted over the different processes, and each process
the hypothetical performance without startup effort and
stores on file only information related to its own the theoretically optimal performance given the fraction
blocks. Every numerical integral is implemented as of the computation that runs in parallel mode
an outer loop over the blocks, inside which the (Amdahls law).
integral-specific work is performed, including an
inner (vector-) loop over the grid points in that par-
ticular block. At the end of the loop-over-the-blocks ts and tp are the (elapsed) times of the serial and par-
the results from all processes are added. This ap- allel parts of the code, respectively, when executed
proach has parallelized in a coarse-grain manner, on a single processor. Tn is the time in an n-node run,
yet very adequately, the numerical integrals and all assuming perfect load balancing and zero commu-
related activities. nication overhead.
A second significant computational aspect of The actual performance suffers not only from
ADF is the determination of the density fit coef- Amdahls law. Imperfect load balancing, communi-
ficients at each cycle of the SCF procedure, as a cation overhead, and startup timesthe time spent
preliminary to the construction of the Coulomb po- in spawning the copies of the program before we
tential. The procedure to compute the coefficients can start computingfurther degrade it. Startup
involves a loop over (all) pairs of atoms. This has times are sometimes ignored when presenting the
been parallelized by distribution of the atom pairs successes of parallelization, but for the end-user
over the parallel processes. only the true elapsed times count. Figure 4 presents
the speedup for a (single-point) calculation on a
Pt(P(Ph)3 )3 CO complex. It also displays the predic-
Performance
tion by Amdahls law, based on timings performed
Not only the numerical integrals and the fit pro- in the ADF code to determine the serial and parallel
cedure have been parallelized but other, less CPU- parts of the calculation.
intensive tasks as well. This has not been a redun- The fair correlation between Amdahls law and
the actual results implies that communication over-
dant and only-academic effort. At first sight it may
head and load imbalance are not yet severe at these
sound impressive when 90% or even 99% of the
numbers of parallel processes. Figure 4 also shows
code runs perfectly in parallel, but such figures im-
that a large part of the deterioration when more
ply that the theoretical (maximum) speed-up is at
nodes are used is due to startup effects: routine
most a factor of 10, respectively 100, even for ex-
PPINIT does little more than spawning the child
tremely large numbers of nodes. This is expressed
processes.
by Amdahls law, which boils down to the simple
fact that however many processes we assign to a
Recent Developments
task, we cannot reduce the time further than to what
is spent in serial mode. The very good load balance, implicit in the results
tp shown in Figure 4 (from the parallelization project
Tn = t s + (9) a few years back), relied on the fact that in the
n
940 VOL. 22, NO. 9

CHEMISTRY WITH ADF
only a transformation of the basis set. If there are

symmetry-equivalent fragments, for example, the
six CO molecules in octahedral Cr(CO)6 , the pro-
gram generates symmetry combinations of the FOs
and uses the symmetrized fragment orbitals (SFOs)
as basis functions. The SFOs transform as the ir-
reducible representations (irreps) of the molecule,
allowing a symmetry-driven analysis of the results.
In absence of any symmetry the SFOs are identical
to the FOs.
The fragment approach offers considerable ad-
vantages. It enhances the interpretative power of
ADF as it leads to a more transparent picture of
bonding, which reduces from a complicated mix-
FIGURE 5. Parallel performance of ADF. Displayed are ing of many primitive basis functions (possessing
the parallelization speedups for the TDDFT module, and little physical relevance) to a few key interactions
separately for the construction of the Fock matrix and the between meaningful fragment (frontier) orbitals.
TDDFT matrix elements. This was already convincingly demonstrated in Ex-
tended Hckel studies,107 and it also holds true for
the KohnSham DFT approach in ADF. The frag-
ADF version of that time each block of integration ment approach also improves the numerical preci-
points carried almost exactly the same workload. sion. In ADF, energies are calculated directly, with re-
Later developments changed that. The performance spect to the fragments, by one single numerical
per block improved by linear scaling techniques, but P integral
of the difference energy density [, r] A A [, r]
not for each block to the same extent. A work-load between the overall molecule and the constituting
measurement procedure, based on counting the fragments.
numbers of basis and fit functions that are relevant Z X
for each individual block of points, has been applied 1E[] = [, r] A [, r] dr (10)
to restore the load balance. Preliminary results for A
a few computationally dominant parts of the code
In other words, we evaluate
R the energy of the over-
(1-cycle SCF and 1-cycle TDDFT) demonstrate that
all molecule, E[] = [, r] dr, and the energies
this approach is successful: Figure 5 is based on tim-
of each of the fragments, say the atomsR that con-
ings with ADF2000 on a cluster of Pentiums.106
stitute the overall molecule, EA = A [, r] dr, in
the same numerical integration grid. This provides
more accurate relative energies than subtracting to-
Functionality tal energies from separate calculations, because the
same relative numerical integration error applies to
ELECTRONIC STRUCTURE AND
BONDING ANALYSIS a much smaller quantity, yielding, in turn, a much
smaller absolute error.
Fragment Approach Note that the user has the freedom to make his
own choice of fragments. This is, however, not a
ADF builds a molecule from user-defined frag- matter of plain arbitrariness, and it does not make
ments, which may be single atoms or larger moi- the analysis tools less meaningful. On the contrary,
eties, for example, ligands, functional groups, or this freedom simply reflects the many perspectives
complete molecules in a donoracceptor complex. from which a particular chemical phenomenon can
In practice, this means that the results of the ADF be viewed.
calculation on a fragment are saved on a file and In practice, many calculations are performed us-
that the fragment files are then used in setting up ing as fragments the so-called basic atoms, which
the calculation on the overall system. The frag- are the smallest possible building blocks in ADF.
ment orbitals (FOs), i.e., the MOs from the calcu- The basic atoms are not necessarily physically real-
lations on the fragments, are employed as basis istic objectsindeed, usually they are not, as they
functions in the new calculation. This does not im- must be spin-restricted and spherically symmetric.
ply a basis set truncation or contraction because The computed (bonding) energy w.r.t. basic atoms,
the virtual FOs are included: the FOs constitute then, does not yield quantities that can be compared

TE VELDE ET AL.
to experimental data directly. Rather, one must cor- prefer to reserve the designation steric interaction
rect for the true ground state of the isolated single or repulsion for 1EPauli because that is, as already
atoms.108 mentioned, the only source of net repulsive inter-
actions between molecular fragments.31 Finally, the
Bond Energy Analysis wavefunction is allowed to relax from 9 0 to the
fully converged wave function 9. The associated or-
In the framework of KohnSham MO theory and bital interaction energy 1Eoi accounts for electron
in conjunction with the fragment approach, one can pair bonding, charge transfer (e.g., HOMOLUMO
decompose the bond energy between the fragments interactions) and polarization (empty/occupied or-
of a molecular systemsay, a base and a substrate
bital mixing on one fragment due to the presence
for E2 eliminationinto contributions associated
of another fragment). This can be further decom-
with the various orbital and electrostatic interac-
posed into the contributions from the distinct ir-
tions. In ADF, we follow a Mokokuma-type en-
reducible representations 0 of the interacting sys-
ergy decomposition method.23, 24, 31, 109, 110 The over-
tem [eq. (13)] using the extended transition state
all bond energy 1E is divided into two major com-
ponents [eq. (11)]. In the first place, the preparation method.112 In systems with a clear / separation,
energy 1Eprep corresponding to the amount of en- this symmetry partitioning proves to be very infor-
ergy required to deform the separated fragments, mative.
A and B say, from their equilibrium structure to X
1Eoi = 1Eoi,0 (13)
the geometry they acquire in the overall molecule
0
(1Eprep,geo ), and to excite them to their valence elec-
tronic configuration (1Eprep,el). In the second place, An extensive discussion of the physical meaning of
the interaction energy 1Eint between the prepared all the terms in the energy decomposition is given
fragments. in ref. 31.
1E = 1Eprep + 1Eint = 1Eprep,geo + 1Eprep,el + 1Eint
(11) Atomic Charges
In the following step, the interaction energy 1Eint is
further decomposed into three physically meaning- Additional insight into bonding and reactivity
ful terms, which are printed in the ADF output file. can be gained if, complementary to the investiga-
1Eint = 1Velst + 1EPauli + 1Eoi = 1E0 + 1Eoi (12) tion of the wave function, one analyses the elec-
tronic charge density distribution. A convenient
The term 1Velst corresponds to the classical electro- way to do so is by computing effective charges
static interaction between the unperturbed charge on the atoms in a molecule. In ADF, three differ-
distributions of the prepared fragments as they are ent approaches are available for routinely comput-
brought together at their final positions, giving rise
ing atomic charges: (1) the Voronoi deformation
to an overall density that is simply a superposi-
density (VDD) method,89, 113, 114 (2) the Hirshfeld
tion of fragment densities A + B . (Note that we
scheme,115, 116 and (3) the well-known Mulliken pop-
use the convention that energy terms containing
ulation analysis.117, 118 Practice has shown that, de-
potential energy only, kinetic energy only, or both
spite its popularity, the latter is strongly basis set
kinetic and potential energy are indicated by V, T,
and E, respectively.) For neutral fragments, 1Velst is dependent. Hence, Mulliken charges are rather un-
usually attractive.31, 111 The Pauli repulsion 1EPauli reliable. We recommend the VDD and Hirshfeld
arises as the energy change associated with going atomic charges instead.
from A +B the wave function 9 0 = NA[9A 9B ] that The VDD method is based on the deformation
properly obeys the Pauli principle through explicit density and a rigorous partitioning of space into
antisymmetrization (A operator) and renormaliza- nonoverlapping atomic areas, the so-called Voronoi
tion (N constant) of the product of fragment wave cells. The Voronoi cell of an atom A is the region in
functions. It comprises the destabilizing interactions space closer to nucleus A than to any other nucleus
between occupied orbitals, and is responsible for (cf. WignerSeitz cells in crystals). The VDD charge
any steric repulsion. In case of neutral fragments, of an atom A monitors the flow of charge into, or
it can be useful to combine 1Velst and 1EPauli in out of the atomic Voronoi cell as a result of turning
a term 1E0 [eq. (12)] which, in the past, has been on the chemical interactions between the atoms.
conceived as the steric interaction. However, we The change upon going from the promolecule, in
942 VOL. 22, NO. 9

CHEMISTRY WITH ADF
which the density is the superposition of the initial, [eq. (16b)] or the overlapping-initial-charge-
unperturbed atomic densities, to the converged SCF distribution Hirshfeld fuzzy cell [eq. (16c)],
density of the molecule, is measured by: respectively [in the latter case, one must R only
Z X realize that ZA in eq. (15) is equal to A (r) dr].
QA = Voronoi (r)
VDD
B (r) dr The VDD and Hirshfeld charges are evaluated by
cell-A B numerical integration in ADF, and both appear to
Z
be rather reliable. They match each other reasonably
= Voronoi 1(r)def dr (14)
cell-A well, are not very sensitive to basis set effects, and
their outcomes conform much better to chemical
The VDD method summarises the three-dimen-
intuition than Mulliken charges do. 113, 114
sional deformation density 1def on a per-atom
It may be noted here that, although the VDD
basis [eq. (14)]. It is conceptually simple, and affords
and Hirshfeld atomic charges conform in practice
a transparent interpretation based on the plausi-
rather well to common chemical reasoning, atomic
ble notion of charge redistribution due to chemical
charges in general, by whatever method defined,
bonding, i.e., the gain or loss of charge in well-
have an element of arbitrariness. Atomic charges
defined geometrical compartments of space.89, 114, 119
are meaningful primarily in the comparison of
See the Highly Polar Electron-Pair Bonding sec-
trends where any arbitrariness in the definition of
tion for an application of the VDD method and a
the charges is eliminated or at least strongly re-
brief description of two recent extensions made in
duced.
refs. 114 and 119.
The Hirshfeld charge of an atom A is defined
EXCHANGECORRELATION (XC) FUNCTIONALS
as the integral of the SCF charge density weighted
by the relative contribution, in each point in space, To calculate the self-consistent solutions of
of the initial unperturbed charge density of that the KohnSham equation the exchangecorrelation
atom in comparison to the density of the promole- (XC) potential VXC is derived from an approximate
cule already mentioned above, i.e., the superposi- expression EXC for the exact XC energy. Various ap-
tion of the initial, unperturbed atomic densities: proximations have been implemented in ADF, at
Z
A (r) the level of the local density approximation (LDA),
QHirshfeld
A = ZA (r) P dr (15) and with corrections from the generalized gradient
B B (r)
approximation (GGA) included. At the LDA level
One may interpret this as picking from the SCF pure-exchange, or Slaters X parameterization can
density the parts that belong (according to the be chosen, or the VoskoWilkNusair (VWN) para-
unperturbed densities) to the atom at hand. The cur- meterization of electron gas data.120, 121 The latter ef-
rent ADF implementation supports Hirshfeld (de- fectively includes correlation effects. Optionally, this
formation) charges only with respect to the initial may further be modified with the correction for the
fragments, rather than to single atoms. The principle correlation between like-spin electrons as proposed
is the same, and in calculations that use single atoms by Stoll.122 At the GGA level, the first implemented
as building fragments they are identical, of course. options included the formulas of Becke (exchange),
The Voronoi deformation density (VDD) and Hir- Perdew (correlation), PerdewWang (both), Lee
shfeld atomic charge of an atom A can both be writ- YangParr (correlation).123 126
ten as the integral over all space of the deformation The fundamentals of density functionals and
density times a weight function wMethod
A (Method = their essential characteristics are a subject of ongo-
VDD, Hirshfeld). ing research. Development of the Van Leeuwen
Z
Baerends potential (LB94)34 targeted the long-range
QMethod
A = 1def (r)wMethod
A (r) dr (16a)
all
space behavior of the XC potential and indeed im-
proved properties that critically depend on the
wVDD (r) = 1 for r Voronoi cell A long-range behavior such as (hyper)polarizabilities
A (16b)
=0 for r
/ Voronoi cell A and high-lying excitation energies. This has been
further improved upon by the orbital-dependent
A (r) SAOP (statistical average of orbital potentials),
wHirshfeld
A (r) = P (16c) which have yielded very good results for response
B B (r)
properties.127, 128 Research on this subject is continu-
The atom-dependent weight function wMethod
A is ing, and new XC potentials are being developed.129
either the hard-boundary Voronoi atomic cell Recently published energy XC functionals devel-

TE VELDE ET AL.
oped elsewhere, such as the meta-GGA by Perdew pect in a calculation. The ZORA formalism is the
and coworkers,40, 41 have been implemented as recommended method in ADF.
well. The relativistic calculation is very similar in
The available functionals do not include B3LYP structure to a nonrelativistic one, having only a few
and other exact exchange functionals when these additional terms in the Hamiltonian.
require multicenter integrals, because their compu-
tation would affect the computational advantage Pauli Formalism
of the DFT technology too much. The newer (pure
In the Pauli formalism relativistic effects are in-
DFT) functionals seem to be at least competitive in
cluded to first order in 1/c2 , where c is the velocity
accuracy.
of light. The Pauli Hamiltonian has relativistic terms
given by
RELATIVISTIC EFFECTS p4 2V 1
HPauli = 3 2
+ 2
+ (V p)
8m c 8mc 4mc2
ADF supports the incorporation of relativistic
HMV + HD + HSO (17)
effects by the quasi-relativistic method based on
the Pauli Hamiltonian63, 64, 130 137 or the more re- The mass-velocity term HMV describes the rel-
cently developed zero-order regular approximation ativistic increase of electron mass with veloc-
(ZORA).65, 66, 68, 69, 90 The Pauli Hamiltonian con- ity, and can be considered a correction to the
taining the first-order relativistic correction terms nonrelativistic kinetic energy. Its matrix elements
(Darwin, mass-velocity and spin-orbit coupling) are
presents a fundamental problem as it is not 2
2 2
hi |HMV |j i = i j (18)
bounded from below. A large part of the rela- 8
tivistic effects can be accounted for by the first- The Darwin term HD is the nonclassical contri-
order perturbation theoretic scheme if the first-order bution of the so-called Zitterbewegung, i.e., the
corrections to the density are taken into account quantum-mechanical effect associated with rapid
(coupled).63, 64 In frozen core calculations one ob- movements of electrons around their average po-
serves a considerable improvement over the cou- sitions, and can be considered the correction to
pled first-order perturbation treatment when the the nonrelativistic potential. The matrix elements
Pauli Hamiltonian is diagonalized in a restricted are
space consisting of just the valence (occupied and
virtual) orbitals from conventionally used double- 2
2
2
hi |HD |j i = i V j + i V j
zeta and triple-zeta STO basis sets. This holds es- 8
X
pecially for very heavy elements like actinides.135
+2 hk i |V|k j i (19)
[See ref. 136 for an extensive discussion and doc-
k
umentation, in particular for metalligand bond
energies.] The diagonalization of the Pauli Hamil- = e /(40 hc) is the fine structure constant,
2
tonian is usually denoted as the quasi-relativistic HMV and HD are the scalar relativistic terms. Their
method.138 For the quasi-relativistic method it is es- matrix elements have the same symmetry as the
sential that the frozen core approximation be used nonrelativistic Fock operator. The extra effort they
and that the valence basis is limited, to avoid a bring along is the calculation of the gradients of the
variational collapse. Indeed, it is sometimes found, orbitals in the integration grid; the Laplacian has to
to the dismay of users that are not fully aware be calculated anyway for the nonrelativistic kinetic
of this aspect, that their Pauli relativistic calcula- energy.
tions give worse results when they apply bigger
Spin-Orbit Coupling
(more accurate) basis sets and/or relax the frozen
core. The spin-orbit term HSO results from the inter-
The ZORA method is generally found to give a action of the electron magnetic moments with the
more accurate description of the relativistic effects. magnetic field generated by their own orbital mo-
It also avoids the problems in the Pauli formalism tion. The computation becomes more involved now
arising from the singular nature of the Coulomb than for the scalar relativistic terms because the
potential, thereby eliminating the above-mentioned spin-orbit term does not have the ordinary mole-
counterintuitive considerations of the basis set as- cular symmetry. Instead, it is invariant under the
944 VOL. 22, NO. 9

CHEMISTRY WITH ADF
molecular point double group, which affects both The ZORA Hamiltonian reads
spin and spatial coordinates. From the orbitals (or c2
basis functions) that are adapted to the irreps of the HZORA = p p+V
2mc2 V
ordinary point group we obtain the spin orbitals X c2
adapted to the double group by coupling the space = pi pi
and spin parts with the (double group) Clebsch 2mc2 V
i
Gordan coefficients. mc2
X + (V p) + V (25)
(2mc2 V)2
1 (r, s) = 0 (r)p (s)h00s p|1i (20)
p Again, this can be written as a sum of terms that
have the (single group) symmetry of the molecule
p (s) are the primitive spin functions (1 = , and a spin-orbit term that is invariant under the
2 = ) and 0s is the double group representa- double group.
tion spanned by these spin functions. The coupled
functions are now adapted to the irrep 1 of the dou- SOLVENT EFFECTS: COSMO
ble group. The spin-orbit matrix elements for these
functions are The conductor-like screening model (COS-
02 X MO)139 141 has been implemented in ADF83 to
01
1r HSO 1s = Bi (1|01 02 ) include solvent effects for reactions and material
i properties in solution rather than in the gas phase
h01 r|i |02 si (21) or the solid state. In this method the environment
effects are treated in an approximate fashion.
We have separated the space and spin part here, The solvent is handled as a dielectric medium
given by that induces charge polarisation on a suitably
2 X defined surface around the moleculethe cavity
h01 r|i |02 s i = i ijk hj 01 r|V|k 02 s i immersed in the solution. Any detailed properties
4
jk of the solvent are reduced to an assumed rigid
(22)
sphere size of the solvent molecules, which aspect
X is only used to determine how far the solvent can
Bi (1|01 02 ) = h1|01 0s pi(i )pq h1|02 0s qi penetrate into the region occupied by the molecule.
pq
This thereby plays a role to determine the effective
(23)
molecular surface, the solvent accessible surface
ijk is the totally antisymmetric LeviCivita tensor; (SAS).
(i )pq are the Pauli spin matrices. The Clebsch Technically the SAS is defined as the envelope
Gordan coefficients h1|01 0s pi can be calculated of a series of spheres.139, 142, 143 These are atom-
from the irreps of the double group. centered spheres with Van der Waals-type radii
Spin-orbit effects can usually be ignored for (a bit larger, usually, than standard Van der Waals
bonding energies and geometry minimizations of radii), augmented by auxiliary spheres, if neces-
closed shell molecules. The complete calculation can sary, to describe the effective cavity that contains
then be done with the ordinary (single) point group the molecule and that is not accessible for the sol-
using only the scalar relativistic terms. The current vent.
implementation does not support geometry opti- The surface-describing spheres are partitioned,
mization with spin-orbit terms included. typically in 60 surface triangles, and each triangle
is assigned a point charge, the strength of which
ZORA is adjustable. Because the spheres defining the SAS
may intersect, only those triangles are used that are
The ZORA method is obtained by rewriting the actually on the surface of the molecule and not,
energy expression and expanding in the parameter for instance, inside one of the other spheres. In the
E/(2mc2 V), which remains small even close to the SCF procedure, the surface-point charges are deter-
nucleus. Retaining only the zeroth order term we mined from the charge density and corresponding
get potential of the molecule itself, and the electrosta-
tic equations assuming, for the moment, a perfect
p 2 c2
E= +V (24) conductor solvent, i.e., a vanishing potential on the
2mc2 V cavity surface. This particular assumption makes

TE VELDE ET AL.
the equations a lot easier and faster to solve. The in a minimization usually keeps us in the local
actual dielectric value of the medium defines a cor- bowl containing the minimum so that the gra-
rection factor for computed properties such as the dients at any new geometry still point more or
energy. less towards the minimum. With transition states,
however, it is much more common that a too-large
FIRST-ENERGY DERIVATIVES step brings us to a location where the energy sur-
face has little resemblance anymore to the surface
Geometry Optimization close to the TS, and the procedure easily loses sight
The (first) derivatives of the energy with respect completely of where the TS might be. The aspects
to nuclear displacements are calculated analytically that discriminate TS search implementations from
for a particular geometry at the end of the SCF a simple minimization are largely related to these
procedure.74, 144, 145 They are used to identify sta- considerations. In the first place, a different Hessian
tionary points in the energy surface, notably for the update scheme is applied: obviously we have no
automatic optimization of the molecular structure. reason to apply a method that favors a positive
This is carried out by a Newton-type iterative pro- definite Hessian. Second, control by the program
cedure, based on a quadratic local approximation of, for instance, a maximum step length to take
of the energy surface and an initial estimate of the at each optimization step is tighter. Third, checks
Hessianbased on a force-field approximation in- are performed to verify that the Hessian has the
corporated in the program or read-in from a restart right structure (one negative eigenvalue) and, if it
file, for instance, obtained from an explicit ADF cal- is found otherwise, adjustments are applied, all this
culation of harmonic frequencies.75, 146 At each step to avoid that the procedure wanders away from the
in the optimization the estimate of the true Hessian TS and loses track.148 Despite such precautions, a
is improved by an updating procedure using the dif- successful search for a TS often requires that the ini-
ference of current and previous gradients in relation tial Hessian is fairly accurate. A force field-based
to the difference in geometries. Various Hessian up- estimate may then be inadequate, and one has to
date schemes are available in ADF. The default (for perform a calculation of frequencies first and feed
normal geometry optimization) is the BFGS (Broy- the result into the TS run, which makes the whole
den Fletcher Goldfarb Shanno) method,147 which procedure more cumbersome and time-consuming.
contains a bias to keep the Hessian positive-definite, A reasonable first estimate of the TS, crucial for
which should be the case near an energy mini- the successful application of the TS search algo-
mum. rithm, can, in many cases, be identified by a linear
transit procedure as described in the next section.
Transition States148
Reaction Path by a Linear Transit
Transition states are of high chemical interest
but computationally often hard to pinpoint. More- One of the options in ADF is to scan a path on
over, the accuracy of the results may be very criti- the energy surface. The path is defined by a sys-
cal when, for instance, the energy barrier between tematic variation of one or more nuclear coordinates
reactants and products is required at kcal preci- between user-specified initial and final values. The
sion. At first sight, the determination of a transition coordinates that are so varied are denoted the lin-
state (TS) is very similar to a minimization. A sta- ear transit (LT) parameters. If there is more than one
tionary point is searched where the gradient vector such parameter, all of them are varied simultane-
vanishes. The Hessian is updated in an iterative ously; the linear transit is therefore often called a
procedure, after some kind of initial estimate. The synchronous linear transit. For each LT point on the
difference between a TS search and a minimization scan, defined by a particular value of the path pa-
is only that at the TS the Hessian has a negative rameter(s), the remaining nuclear coordinates (or a
eigenvalue. subset of them), i.e., those that are not LT parame-
When we are close enough to the TS, the stan- ters, may be optimized. From a technical point of
dard procedure used for minimization will produce view, a LT run is simply a sequence of constrained
the TS. However, the energy surface around the TS optimizations.
often displays sizeable anharmonicity. Neglecting One of the most useful applications of a LT is to
third- and higher-order terms, as inherent to the get a reasonable estimate of a transition state. A clas-
Newton approach, is therefore likely to be less ac- sical example is the HCN to CNH transit, in which
curate. Moreover, taking wrong or too large steps the hydrogen atom travels from one end of CN to
946 VOL. 22, NO. 9

CHEMISTRY WITH ADF

the other. The LT parameter would here be the an- qi = x(1) m(1) , q2 = y(1) m(1) ,
(27)
gle (HCN), changing from 180 degrees down to q3 = z(1) m(1) , q4 = x(2) m(2) ,
zero in 10 steps, say. For each particular -value, the
HC and CN distances are optimized. By gathering Diagonalization of F, after projecting out the rigid
the energy values at the LT points and interpolating translations and rotations, gives the normal modes
them in a curve we obtain a fair view of the path, as eigenvectors. The eigenvalues are directly related
and hence, the location of the TS. to the harmonic frequencies.

k
Intrinsic Reaction Coordinate (IRC) vk = (28)
2
For insight in the properties of a full reaction path When building the matrix of force constants, one
the LT functionality provides only very crude in- can use Cartesian or internal (Z-matrix) displace-
formation. An alternative and more advanced pro- ments. In both cases the final matrix is transformed
cedure to trace the path of a chemical reaction is to the representation in mass-weighted Cartesians
provided by the intrinsic reaction coordinate (IRC) and then diagonalized to obtain the frequencies and
method149 as implemented in ADF.76, 150 The basic normal modes.
idea151 is to start at a transition state and slide It is possible to run a calculation of frequen-
down the hill towards the adjacent local minimum, cies where only a subset of the nuclear coordinates
at either side of the transition state, by following is displaced, thereby neglecting any coupling with
the line of steepest decent in the metric of mass- the other coordinates. This is a fast way to study
weighted coordinates. A full determination of the a few well-selected normal modes if one knows
path does not only produce the height of the bar- a priori that the omitted displacements will hardly
rier in a reaction, but also provides information contribute to the mode(s) one is interested in.
for thermodynamic properties, the reaction rate and In the computation of second energy derivatives
dynamics;152 the details of such analysis have not by numerical differentiation of first derivatives in
yet been implemented in ADF. slightly differing geometries, accuracy of the gradi-
The calculation starts by taking steps of con- ents is crucial. This puts high demands on numeri-
trolled length (in the mass-weighted metric) to- cal integration precision, SCF convergence, etc. Fur-
wards the minimum by determining the minimum thermore, for each of the displaced geometries a full
on the local sphere around a given point along the SCF has to be carried out to determine the gradients.
path, with as radius the chosen step length. The cal- This makes the computation of frequencies by nu-
culation, therefore, consists, as in the linear transit merical differentiation a demanding task, although
method, of a sequence of constrained optimizations, in the ADF program use is made of symmetry148, 153
and terminates when the adjacent minima to the TS so that only symmetry unique atoms are displaced.
have been identified. There are indications that ADF frequencies are
slightly superior to frequencies computed with pro-
grams that use Gaussian-type basis of comparable
SECOND DERIVATIVES: FREQUENCIES
AND IR SPECTRA size.154
Frequencies Analytical Second Derivatives
The original (and still available) implementation Use of analytical second derivatives avoids the
of harmonic frequencies in ADF75, 153 computes sec- numerical differentiation of gradients, which is
ond derivatives of the energy with respect to nu- rather critical for precision of the results, and it re-
clear displacements by numerical differentiation of moves the necessity to perform many SCF cycles.
the analytically computed first derivatives (gradi- At the other hand, one needs to solve the cou-
ents) in the geometry at hand and a series of slightly pled perturbed KohnSham equations where both
modified geometries. Going over all possible dis- derivatives are taken with respect to the nuclear
placements, we obtain the matrix of force constants, coordinates.155, 156 Analytical second derivatives are
which we can express in mass-weighted Cartesian faster and more robust and accurate, and have re-
coordinates qi . cently become available (in the export version of the
package) as a separate program that takes the nor-
2E mal (TAPE21) result file from an ADF calculation as
Fij = (26)
qi qj input.

TE VELDE ET AL.
Infrared (IR) Intensities75, 153 This symmetric 33 tensor relates the magnetic field
that is externally applied in the NMR experiment,
As a by-product in the calculation of second en- Bext , to the effective magnetic field that is felt at the
ergy derivativesharmonic frequenciesone ob- nucleus, Beff
tains the first derivatives of the dipole moments,
yielding the infrared (IR) intensities. Beff = Bext (31)

2 The isotropic shielding is the average of the di-
Ai = (974.86)gi (29)
Qi agonal elements of the tensor. The difference in
Qi is the normal mode (mass-weighted Cartesians isotropic shielding of a reference material and the
in atomic units, masses in atomic mass units) sample gives the NMR shift of a standard NMR
with degeneracy gi ; Ai is the absorption intensity experiment. The program calculates the absolute
in [km/mol] of normal mode Qi , and is the dipole shielding, as it has no knowledge of reference ma-
moment in atomic units. IR intensities are obtained terials.
both with the analytic and with the numeric ap- In the full equation for the tensor, we can discrim-
proach to second derivatives. inate three terms: the paramagnetic and diamag-
Raman intensities are calculated similarly from netic contributions, and a spin contribution. The
the derivatives of the polarisability tensor with re- latter is the Fermi-contact term if the Pauli spin-
spect to the normal modes (see also TDDFT section). orbit Hamiltonian is used, or the spin-orbit coupling
contribution in the ZORA approach. The three con-
ELECTROMAGNETIC PROPERTIES tributions to the tensor are calculated separately
(the spin-orbit term only in a spin-orbit calculation).
NMR Chemical Shifts
The GIAO formalism makes it possible to analyze
The ADF package allows the calculation of the the contributions to the shielding in terms of orbital
shielding tensor of nuclear magnetic resonance contributions.158 163
(NMR) spectroscopy.77, 78, 80 82, 157 A separate NMR
module that reads information about the molecule NMR Spin-Spin Coupling
from the standard TAPE21 output file of an ADF run
calculates the shielding data. The method has been Implementations are under way of methods
formulated in the framework of gauge including that will calculate NMR spinspin coupling
atomic orbitals (GIAO) and includes relativistic ef- constants.164 167
fects (ZORA or Pauli, scalar or including spin-orbit
effects).
The calculation of NMR parameters constitutes ESR
a theoretically challenging task, in particular when
applied to heavy elements. An earlier investiga- ADF supports the computation of Electron
tion suggested that the frozen core approximation Spin Resonance (ESR) properties in systems with
attractive because of its efficiency and necessary for one unpaired electron: the g-tensor, A-tensor (nu-
relativistic calculations with the Pauli formalism clear magnetic dipole Hyperfine interaction), and
would not be appropriate for shielding calculations. Q-tensor (nuclear electric quadrupole hyperfine
However, when terms that contain the core orthog- interaction).168 172 Work on ESR calculations for
onalization coefficients to first order in the magnetic systems with more than one unpaired electron is in
field are taken into account, as in the ADF imple- progress.
mentation, the differences between frozen core and
all-electron calculations are found to be small. One TIME-DEPENDENT DFT
has to take care, however, that the valence includes
at least the highest atomic shells and the next to Time-dependent DFT (TD-DFT)173 has been im-
highest (p- and d-shells). plemented in ADF174, 175 for the research of response
The NMR shielding tensor is the second deriva- properties, such as excitation energies, frequency-
tive of the energy with respect to the magnetic field dependent (hyper)polarizabilities and related prop-
Bk in each direction k = x, y, z, and the magnetic mo- erties of (at the moment) closed-shell molecules. The
ment of the nucleus t in each direction t = x, y, z implementation exploits the molecular symmetry
2E and typically requires similar amounts of CPU time
kt = (30) as an ordinary SCF calculation in ADF.
Bk t
948 VOL. 22, NO. 9

CHEMISTRY WITH ADF
The time-dependent KohnSham equations are equation becomes

similar to the time-independent ones
X
occ X
virt
0
2 s (r, r , ) = a (r)i (r)a (r0 )i (r0 )
i i (r, t) = Hi (r, t) = + V[](r, t) i (r, t) i a
t 2 !
(32) 2(i a )
X 2 (36)
(r, t) = nj j (r, t) (33) (i a )2 + 2
j
(Note: if we replace +2 by 2 , we obtain the
The potential V comprises the Coulomb poten- polarizabilities at imaginary frequencies, which are
tial due to the electron cloud, the nuclear and used for the determination of Van der Waals disper-
any other external potentials, and an effective sion coefficients.)
exchangecorrelation (XC) potential, all of them In the TDDFT case, one has to consider the re-
time-dependent. In principle, the effective poten- sponse of the noninteracting KS system (described
tial is different from the one defined in the KS by the single-particle KS response function) to an ef-
formulation of ground-state DFT. However, it is a fective field, schematically written as:
good approximation (the adiabatic approximation)
in many cases to assume that the functional depen- = s Veff (37)
dence of the time-dependent effective potential on
The change in the effective potential Veff in eq. (37)
the time-dependent density is the same as the de-
is split in its Coulomb, external, and XC compo-
pendence of the time-independent potential on the
nents.
time-independent density. Z
The objective is to determine the response of the (r0 , )
Veff = dr0 + Vext (r, ) + VXC (r, ) (38)
systemthe change in the densityto a change in |r r0 |
the external field. Using a perturbative approach,
The latter is related to the so-called XC kernel
the density is expanded in a functional Taylor series
fXC , which is the functional derivative of the time-
(in the perturbing potential). Many interesting prop-
dependent XC potential with respect to the time-
erties, for instance, the polarizability and excitation
dependent density.
energies, depend only on the first-order cnange in
Z
the density.
VXC (r, ) = dr0 fXC (r, r0 ; )(r0, ) (39)
The first-order density depends on the (exact) lin-
ear response function . After a Fourier transforma- d2
tion we can rewrite the time-dependent equations in fXC = 2
XC [] = (40)
d 0
frequency (energy)-dependent form
Z Here XC is the XC energy functional.
1 (r, ) = dr0 (r, r0 , )Vext (r0 , ) (34) In ADF, the adiabatic (local density) approxima-
tion (ALDA) is used for the kernel. Although this
In the KohnSham formulation of time-dependent kernel, which is then the functional derivative of the
DFT the response function is replaced by the re- ordinary LDA XC potential, is local in both space
sponse function s of the noninteracting Kohn and time, it has performed quite well in practical
Sham system, while the external potential is re- calculations. There are strong indications128, 129, 176
placed by the effective KohnSham potential. The that the approximation chosen for the usual XC po-
response function s can be expressed in the unper- tential in the SCF part of the calculation has much
turbed KS orbitals more impact on the accuracy of calculated response
properties.
X j (r)k (r)j (r0 )k (r0 ) Equations (37)(40) are solved self-consistently
s (r, r0 , ) = (nk nj ) (35)
(j k ) + i in an iterative fashion, using the external potential
jk
as initial guess for the effective potential; in other
If we consider the real density response in an ap- words, in the first iteration one neglects all XC and
plied electric field, we can choose the complex KS Coulomb screening effects. Rapid convergence is
orbitals to be real-valued and the infinitesimal can obtained if DIIS techniques are used in the iterative
be set to zero. Assuming, furthermore, that the oc- update procedure, similar as for the ordinary SCF
cupation numbers nj are all either 1 or zero, the cycles in ADF.

TE VELDE ET AL.
Linear Response Properties that uses the result files of response calculations on
the separate monomers.177
The linear response of a molecule to a frequency
dependent electric perturbation is given by the lin- Excitation Energies and Oscillator Strengths
ear polarizability tensor ij ().
X TDDFT offers a rigorous theoretical framework
i = 0i + ij ()Ej () (41) for the calculation of excitation energies in DFT, and
j is preferable in this respect to 1SCF approaches to
(i 0i ) is the change in the dipole moment of the reach the same goal.175, 178, 179 The basic equation for
molecule in direction i. Having solved the SCF equa- the determination of excitation energies and oscilla-
tions for the first-order change in the density, we tor strengths is the following eigenvalue equation
obtain the polarizability by Fi = i2 Fi (47)
Z
2 The linear dimension of the matrix is the prod-
ij () = dr i (r, )rj (42)
E uct of the numbers of virtual and occupied orbitals,
The poles of the polarizability tensor are the (verti- respectively. Because of its large size, this matrix
cal) excitation energies, while the strengths of these cannot be diagonalized directly. Instead, the iter-
poles are given by the oscillator strengths. The rela- ative Davidson procedure is invoked to obtain a
tion between the oscillator strengths, excitation en- few of the lowest excitation energies and oscilla-
ergies, and average dipole polarizabilities is given tor strengths of each irrep. The matrix contains
by: the same Coulomb and XC terms as are needed
for the calculation of the polarizability tensor. In
X fi
a () = (43) fact, the average dipole polarizability and its fre-
i
i
2
2 quency dependence (determined by the so-called
Cauchy coefficients) can also be determined if a di-
The polarizability also leads directly to the
rect diagonalization of the matrix is possible.
(isotropic) Van der Waals dispersion coefficient C6
The matrix is blocked by symmetry, and can
that determines the long-range dipoledipole dis-
be split in singlet and triplet parts, each of which
persion interaction between two systems, A and B.
is handled in turn in the implementation. Oscilla-
(This applies to atoms; molecules have also an
tor strengths, which determine the magnitudes of
anisotropic contribution to the R-6 interaction.)
the peaks in an absorption spectrum, are calculated
C6 (A, B) from the solution vectors Fi . An example of a chemi-
Edisp = (44)
R6AB cal application of this part of the code is given in the
Z second part of this article.
3
C6 (A, B) = A (i)B (i) d (45)
0 Nonlinear Response
The isotropic polarizability is the average of the
diagonal elements of the tensor. Higher order optical response is governed by the
Instead of a perturbative dipole field Er cos(t), hyperpolarizability tensors , , etc.
we may apply a more general multipole electric X j 1 X j
field. This yields the general multipolemultipole i = 0i + ij E0 + ijk E0 Ek0
2!
polarizability mm0 , which describes the response of j jk
a general molecular multipole moment qlm to such a 1 X j

+ ijkl E0 Ek0 El0 + (48)
field. 3!
X 0 0
jkl
qlm = qlm
0 + mm0 El m (46) is the first hyperpolarizability tensor (which van-
l0 m0
ishes for molecules with an inversion center), is
The l = 2 (quadrupole) and l = 3 (octupole) cases the second hyperpolarizability tensor, etc. If the ex-
have been implemented in ADF.177 In fact, one can ternal field is static, these tensors can be calculated
obtain all components mm0 up to l = 3 for a range with finite field techniques. In the experiment one
of frequencies in a single calculation. Similar to the always measures frequency-dependent hyperpolar-
isotropic dispersion coefficients of eq. (45), one can izabilities. For a direct comparison to experiment it
obtain all isotropic and anisotropic dispersion coef- is, therefore, important to include this frequency de-
ficients between two molecules from a calculation pendence as implemented in ADF.180
950 VOL. 22, NO. 9

CHEMISTRY WITH ADF
Derivatives w.r.t. Nuclear Coordinates
The variation of spectroscopic properties upon

nuclear displacement gives access to several prop-
erties. To compute this geometric dependency for
small displacements we need the geometric deriv-
atives at the equilibrium geometry. These are ob-
tained as in a normal calculation of harmonic FIGURE 6. Lewis structures of linear CN (1) and
frequencies, namely by numerical differentiation, CP (2) dimers.
using the data at the equilibrium geometry and at
slightly displaced geometries (some of the appli-
cations have yet to be implemented). Derivatives solid state and chemistry in solution. In the re-
of the dipole moment give rise to the IR intensi- mainder of this article, we elaborate first on two of
ties. The derivatives of the polarizability tensor chemistrys core issues: the chemical bond (this sec-
likewise yield the intensity and depolarization ra- tion) and chemical reactions (next section). There-
tio of Raman scattering (as a function of the laser after, an application in the field of biochemistry is
frequency in our case). Derivatives of the hyperpo- briefly discussed, namely, the structure and bond-
larizability tensors govern the hyper-Raman effect, ing of deoxyribonucleic acid, DNA. The emphasis
while derivatives of the excitation energies can be is on gaining insight, on being able to tell how a
used to estimate the intensity of the resonance Ra- certain overall effect is brought about by the under-
man effect. So far, only the IR and Raman intensities lying electronic mechanisms. More than once this
have been implemented in ADF. reveals surprising differences between seemingly
Geometrical derivatives are not only interesting similar phenomena, and vice versa. This is of great
because they give rise to additional properties, but practical importance. Not only does it yield a deeper
they also determine the vibrational average of a understanding of the investigated molecule, but it
quantity, which is usually what is measured experi- also allows a more rational tuning of the system to
mentally. achieve the desired properties.
In this chapter, we look into the elementary con-
Combination with Relativistic Options cepts of the electron-pair bond, the three-electron
bond, and the donoracceptor (or charge transfer)
All functionality of the RESPONSE module can bond in various chemical situations, and we eval-
be combined with the scalar relativistic options of uate these concepts within the framework of the
ADF, such as the scalar Pauli and scalar ZORA op- quantitative KohnSham MO model. These bond-
tions. It is not yet possible to include spin-orbit ing modes may interfere with and can be substan-
effects on response properties. If scalar relativis- tially affected by Pauli repulsion with closed shells
tic effects are taken into account, the relativis- (or same-spin electrons in open shells) and by clas-
tic orbital energies and orbitals are used in the sical electrostatic interactions.
post-SCF excitation energy or polarizability calcu-
lation. No further relativistic effects are included THE ELECTRON-PAIR BOND
in the calculation of the density response. By de-
fault the so-called scaled ZORA orbital energies An interesting example of similar trends that
are used (see ref. 181 for a discussion and an appli- originate from very different electronic mechanisms
cation). is provided by the two valence iso-electronic series
of CN and CP dimers,182, 183 schematically shown in
Figure 6, series 1ac and 2ac, respectively.
Part II. Applications: We have studied these systems at the BP86/TZ2P
The Chemical Bond level. All three CN dimers have linear equilibrium
structures and it is only the most stable CP dimer,
The scope of applicability of modern DFT codes PCCP (2a) that acquires such an unbent geometry.
is rapidly expanding as a result of extended func- In fact, CPCP and CPPC actually tend toward non-
tionality and increased efficiency of the implemen- linear geometries, the linear structures 2b and 2c
tations, reaching such diverse fields as organic, being second- and fourth-order saddle points. This
organometallic, and inorganic chemistry and even has been investigated and analyzed in detail.183 In
biochemistry, covering the gas phase as well as the the following, we focus on the linear species, and

TE VELDE ET AL.
in particular, on the central electron-pair bond that

holds together the monomers in 1 and 2.
The strength of the central bond decreases in both
series of CX (X = N, P) dimers if we go from cou-
FIGURE 8. The CN and CP dimers can be considered
pling via carbon to coupling via the heteroatom. ambident radicals.
The bond dissociation enthalpies (BDE) for 298 K
are 136.6, 113.5, and 68.2 kcal/mol along 1ac, and
amplitude at both ends of the diatomic. The HOMO
152.2, 68.7, and 7.5 kcal/mol along 2ac. In the CN
has somewhat more amplitude on the heteroatom
dimers, the decrease in bond strength goes with a
(it is the N or P lone pair orbital), whereas the
shortening of the bond (1.373, 1.305, and 1.274
SOMO has a slight bias toward the carbon atom.
along 1ac), whereas in the CP dimers it is accom-
In other words, CN and CP are ambident radicals
panied by bond elongation (1.336, 1.709, and 2.216
(see Fig. 8) that may, in principle, form electron-pair
along 2ac). The relatively short CC bond dis-
bonds (Fig. 9) either via carbon or via the hetero-
tance in 1a and 2a is indicative of a partial multiple-
atom, leading to 1ac and 2ac.
bond character, which is confirmed by our analysis.
As the SOMO has a bigger lobe on carbon, one
Thus, the as-yet unisolated PCCP (2a) is thermody-
may expect that the strength of the electron-pair
namically stable. However, as follows from further
bond (Fig. 9) and the stability of the dimer decreases
explorations, 2a may tend to polymerise. Synthetic
in the order XCCX > CXCX > CXXC (X =
strategies toward this molecule have to cope with
N, P).
this instability, for example, through coordination
At first sight this is nicely confirmed by the de-
of the product to transition metals.183 The analy- creasing bond strengths along the series 1ac and
sis, furthermore, shows that the contraction along 2ac (Fig. 8). However, a more quantitative analysis
1ac is due to the combined effect of the smaller ef- shows that, below the surface of the overall bond
fective size of nitrogen compared with carbon, and strength, these similar trends have quite different
to its higher electronegativity. The singly occupied origins. In the CN dimers, the reduction in bond
molecular orbital (SOMO) of CN, which is weakly strength is caused by an enormous increase in the
antibonding, has a lower amplitude on nitrogen. 1E0 repulsion (see the Bonding Energy Analysis
These effects lead to an onset of both repulsive and section) from 144.3 to 371.3 kcal/mol, which is not
bonding interactions at shorter bond length if N in- fully compensated by a sizeable increase in 1Eoi
stead of C is involved in the central bond. The fact from 288.2 to 446.5 kcal/mol. In contrast, in the
that the central bond along 2ac nevertheless elon- CP dimers, the bond strength decreases because of a
gates when CP binds via P instead of C is simply due weakening of the bonding orbital interactions 1Eoi
to the more diffuse character and larger effective (from 320.9 to 104.1 kcal/mol), especially those
size of phosphorus compared to carbon or nitro- associated with the bond (1E ) despite an op-
gen. posite trend of the repulsion 1E0 , which actually
We may wonder what cause underlies the similar decreases (from 161.5 to 91.0 kcal/mol).
trends of decreasing BDE for the CN and CP dimer The increase in Pauli-repulsion along the CN
sequences 1 and 2. The relevant frontier orbitals of dimers 1ac is easy to understand. It is caused by the
the CX monomers, i.e., the HOMO and SOMO in increase in overlap between the closed-shell HOMO
symmetry schematically represented in Figure 7 (for orbitals (Fig. 7) and between the HOMO orbitals (not
contour plots, see ref. 183), are essentially CX non- shown), which have higher amplitudes on nitro-
bonding orbitals.
The former provides the axial N or P lone pair,
whereas the latter carries the unpaired electron. An
important feature of the CN and CP frontier or-
bitals is their delocalized nature, with significant
FIGURE 9. electron-pair (left) and donoracceptor

(right) orbital interaction diagram for CN (1) and CP
FIGURE 7. Frontier orbitals of CX radicals, X = C, P. (2) dimers.
952 VOL. 22, NO. 9

CHEMISTRY WITH ADF
gen. The increase in stabilizing orbital interactions

along 1ac, however, does not correlate with the
bond overlap hSOMO |SOMO i that, as qualitatively
predicted, decreases from 0.46 to 0.21. Here we see
a subtle interplay of secondary SOMO HOMO in-
teractions that interfere with the pure electron-pair
bond. To first order, the primary SOMO + SOMO in-
teraction decreases in strength along 1ac (Fig. 6) FIGURE 10. Simple VB resonance schemes for linear
CN (1) and CP (2) dimers.
and the SOMO + SOMO in NCCN (1a) descends
strongly, and comes out close to (in fact, even below)
the HOMO + HOMO . This gives a strong mutual re-
corresponding set of unoccupied LUMO s on the
pulsion and a corresponding weakening of the elec- other CX monomer (Fig. 9D). In CNNC (1c), the
tron pair bonding interaction contained in 1Eoi . In orbital interactions (85.8 kcal/mol) are even larger
contrast, in CNNC (1c) the occupied HOMO HOMO than the net bond enthalpy of 68.2 kcal/mol (the
is strongly destabilized, due to a sizeable overlap orbital interactions in 1c are 360.7 kcal/mol).
and repulsion between the N lone pair orbitals, Therefore, the central bond in the CX dimers is best
and it approaches in first order the unoccupied conceived as having partial triple bond character, as
SOMO SOMO from below. The resulting donor shown in E in terms of a resonance between simple
acceptor interaction is about twice as strong as that VB structures (Fig. 10).
in 1a and 1Eoi in 1c is, hence, significantly more
bonding.182 Without this secondary interference of THE THREE-ELECTRON BOND AND ITS
the pure electron-pair bond with other closed shell ONE-ELECTRON BONDING COMPONENT
orbitals the differences in bond strength along 1ac
would have been much larger, i.e., NCCN would The traditional picture in MO theory of the two-
have been more and CNNC less stable. center three-electron (2c3e) bond between a closed-
It remains to understand why the repulsive and shell and a radical fragment, A and B say, is that of
bonding orbital interactions in the three linear CP a doubly occupied bonding MO ( ) and a singly
dimers follow the opposite trend with respect to occupied antibonding MO ( ) of the composite
those in the CN dimers, i.e., why both 1E0 and 1Eoi molecule AB (Fig. 11).
decrease along 2ac. In the first place, the CP orbitals It follows from Figure 11 that the 2c3e bond may
are more extended and diffuse at the P side. So, be viewed as composed of an electron-pair bond
the overlaps are smaller and achieve a lower max- counteracted by a destabilizing antibonding elec-
imum at a longer bond distance when phosphorus tron leading formally to a bond order of 1/2. The
gets involved in the central bond. Consequently, the above analysis considers the overall molecule and
Pauli-repulsion contained in 1E0 and, even more so, the properties of its MOs.
the bonding orbital interactions 1Eoi decrease along Recently, we have proposed a different
2ac instead of increasing as they do along 1ac. The although equivalentapproach to the 2c3e bond,
other important point is that the primary or first- namely from the perspective of the fragments
order electron pair bond is much less affected by and their mutual interaction.184 We take the
secondary interactions with other orbitals: the de- sulphursulphur three-electron bond in H2 SSH2 +
crease in SOMO + SOMO electron-pair bonding, if as an example (H2 SSH2 + is the only minimum
we go from 2a to 2c, is no longer compensated by a on the MP2/6-31G potential energy surface of
strongly stabilizing HOMO /SOMO interaction. This H4 S2 + 185, 186 ). This species is the archetype of
is simply due to the fact that the HOMO HOMO
does no longer come close enough to the empty
SOMO SOMO because of the smaller HOMO
HOMO splitting in the CP dimer, and because of the
larger HOMO /SOMO gap in the CP monomer (3.5 eV
in CP compared to 2.4 eV in CN).
Finally, a decomposition of the bonding orbital
interactions shows a smaller but significant contri-
bution of donoracceptor-type -bonding between
a (doubly degenerate) set of occupied HOMO or- FIGURE 11. The two-center three-electron
bitals on one CX monomer (X = N, P) and the (2c3e) bond.

TE VELDE ET AL.
dialkylsulfide dimer radical cations (R2 SSR2 + ),

a class of systems for which there is also much
interest in gas-phase experimental studies.187 190
At the BP86/TZ2P level we find the H2 SSH2 +
system bound by a sulphursulphur bond of 2.886
with a bond enthalpy for 298 K of 40.7 kcal/mol.
This bond enthalpy is some 10 kcal/mol stronger
than the best ab initio estimate [31.9 kcal/mol
at CCSD(T)]. This relatively large error for the
nonlocal functionals can be attributed to a well-
known deficiency of the existing exchange function-
als to properly cancel the self-interaction part of the FIGURE 12. Decomposition of the two-center
three-electron bond.
Coulomb energy in case of delocalized ionization
out of symmetry equivalent weakly (or non-) over-
lapping orbitals.191, 192 To asses if this deficiency of
bond energy. This places the physical model on a
the functionals might affect the physical model of
more quantitative basis. In particular, it becomes
the 2e3e bond, we have carried out the quantitative
clear that the three-electron bond is, in principle,
analysis at various levels of DFT, local (X-VWN) as
weaker than the corresponding one-electron bond-
well as nonlocal (BP86 and PW91). It appears that,
ing component. In our H2 SSH2 + model system, for
although numerically different, the relative propor-
example, the three-electron bond (1E2c3e = 25.9
tions of the different physical terms (1Velst , 1EPauli ,
kcal/mol) is about twice as weak as the correspond-
1Eoi ) in the SS interaction are very similar for all
ing one-electron bonding component (1E2c1e =
three levels, yielding the same physical picture. The
1Eoi = 51.3 kcal/mol) owing to the Pauli repul-
subsequent discussion is based on the BP86/TZ2P
sive term (1EPauli = 25.4 kcal/mol).
calculations, which give the least overbinding.
The sulphursulphur bond in H2 SSH2 + is THE ROLE OF PAULI STERIC REPULSION
mainly provided by the three-electron orbital inter- IN BONDING MODELS
actions between the 1b1 orbitals of the two frag-
ments, i.e., the 3px lone-pair of H2 S and the 3px So far, we have looked at different modes of
SOMO of H2 S+ (e.g., Fig. 11). Combining the re- bonding and how Pauli-repulsive orbital interactions
pulsive and bonding orbital interactions into one may either influence them (e.g., the electron-pair
term, 1E2c3e = 1EPauli + 1Eoi , one arives at a three- bond) or be an essential part (2c3e bond). In this
electron interaction 1E2c3e of 25.9 kcal/mol. This section, we examine a different role of Pauli re-
is about 60% of the net interaction 1Eint . However, pulsion, namely its role as the source of steric
the electrostatic interaction 1Velst although clearly repulsion between nonbonded groups. An interest-
smaller than 1E2c3e , is an important component ing example is the methyl radical and its heavier
also: with 18.4 kcal/mol it still contributes some- Group-14 congeners, AH3 (A = C, Si, Ge, Sn).
what more than 40% to the net sulphursulphur The methyl radical is known to be planar (D3h ),
interaction (1Eint = 44.3 kcal/mol). whereas the heavier AH3 radicals adopt pyrami-
Figure 12 illustrates how the three-electron dal (C3v ) structures. This is nicely confirmed in a
bond (G) can be interpreted in the fragment ap- recent BP86/TZ2P study,193 which shows that the
proach. The repulsive term 1EPauli contained in HAH angle in AH3 decreases ( = 120, 112.7,
1E2c3e is mainly due to the destabilizing interac- 112.4, and 110.6), whereas the inversion barrier in-
tion H between the unpaired 3p electron on H S+
x 2 creases along A = C, Si, Ge, and Sn (0.0, 3.7, 3.8, and
and the same-spin 3px electron of the lone-pair on 7.0 kcal/mol).
H2 S (the excess spin is arbitrarily chosen ). The On the basis of the qualitative MO theory, this
bonding orbital interaction 1Eoi is simply provided was traditionally explained through the operation
by a one-electron bond I between the 3px elec- of a second-order JahnTeller effect (Fig. 13).194, 195
tron of the lone-pair on H2 S and the empty 3px The mixing between the nonbonding npz SOMO
orbital on H2 S+, i.e., 1Eoi ' 1E2c1e . In this way, and the AH antibonding LUMO stabilizes and
the three-electron bond is naturally linked to the pyramidalizes AH3 . This effect is more pronounced
one-electron bond plus a characteristic Pauli repul- for the heavier (more electropositive and diffuse)
sive term. Important also is that we are able to central atoms A. The SOMOLUMO gap becomes
quantify the different energy terms in the overall smaller due to the higher energy of the npz SOMO
954 VOL. 22, NO. 9

CHEMISTRY WITH ADF
FIGURE 13. Walsh diagram for the pyramidalization FIGURE 14. Schematic orbital interaction diagram
of AH3 . for a planar (D3h , left) and a pyramidal (C3v , right)
AH3 radical.
and the less AH antibonding nature of the LUMO.

The JahnTeller effect is opposed by the increas- AH bond and HH distance in CH3 , which
ing energy of the 1e1 orbitals, which is ascribed leads to a larger HH h1s|1si overlap (0.32 com-
to the loss of AH bonding overlap. Thus, only pared to, e.g., 0.16 in SiH3 ). The HH repulsion is
CH3 remains planar because the JahnTeller effect slightly relieved and partly hidden when the AH
is not strong enough in this case to outweigh the bond is allowed to adapt, by elongation, to the pyra-
1e1 destabilization. Thus, the emphasis in this tra- midal geometry.
ditional model is on the trend in electronic effects, The AH interaction is provided by the forma-
i.e., bonding orbital interactions, that are not strong tion of three electron-pair bonds, one in A sym-
enough to pyramidalize CH3 . metry, and two degenerate ones in E symmetry,
A careful KohnSham MO analysis (at BP86/ as shown in Figure 14. Only if AH3 is pyrami-
TZ2P) leads to a substantial revision of this pic- dal (Fig. 14, to the right), the singly occupied npz
ture. To understand the trends in pyramidalization AO of the central atom can be stabilized by a
of the four AH3 radicals, we have analyzed how bonding overlap with the 1a01 orbital of the H3
the HH and AH interactions in AH3 , 1Eint (H3 ) fragment. The bonding AH orbital interactions
and 1Eint (AH3 ), change if we go from the op- in A symmetry hence favor pyramidalization. The
timized planar structure to a pyramidal geometry. overall AH bonding interaction 1Eint (AH3 ) is
The AH bond distance is kept fixed to its value largest for CH3 , but the additional stabilization
in the planar radical. The HAH angle is bent 1 1Eint (AH3 ) upon pyramidalization is the weak-
to its value in the optimized pyramidal structure est, namely 0.4, 5.3, 4.7, and 6.3 kcal/mol
(for CH3 the optimum HAH angle of SiH3 along CH3 , SiH3 , GeH3 , and SnH3 . Note that the
was used, because there is no stationary point cor- bonding orbital interactions in the methyl radical
responding to a pyramidal methyl radical). still favor pyramidalization, i.e., 1 1Eoi (CH3 ) =
As shown in Figure 14 (to the left and the right), 1.9 kcal/mol; only the 1 1E0 (CH3 ) contribution
the H3 fragment in AH3 has a (1a01)1 (1e01 )2 valence causes the overall change in CH interaction to
state, i.e., one electron in the bonding 1a1 combina- slightly disfavor a pyramidal methyl radical. The
tion of the three hydrogen 1s AOs and one electron trend in 1 1Eint (AH3 ) along A = C to Sn follows
in each of the degenerate antibonding 1e1 combi- that of the increasing gain in AH3 overlap and sta-
nations. This gives rise to steric HH repulsion bilizing orbital interaction 1 1EA between the npz
that can be conceived as Pauli repulsion between AO of the central atom and the 1a01 combination of
three same-spin electrons. There is a striking differ- the three hydrogen atom 1s AOs (i.e., 1s + 1s + 1s).
ence between CH3 and the heavier homologs. In The interactions in E symmetry are relatively insen-
CH3 the HH steric repulsion 1E0 (H3 ) and, hence, sitive to the pyramidalization process.
1Eint (H3 ) is significantly stronger, and increases In conclusion, our KohnSham MO approach
much more upon pyramidalization: the HH in- shows that the CH3 radical is planar mainly be-
teraction becomes more repulsive by 4.0 kcal/mol cause of the steric repulsion between the hydrogen
for CH3 and only by 1.1 kcal/mol or less for SiH3 , ligands. The steric HH repulsion is much weaker
GeH3 , and SnH3 . The reason is simply the shorter for the heavier central atom homologs in which the

TE VELDE ET AL.
ligands are farther removed from each other. Elec-

tronic effects (i.e., electron pair bonding between
central atom and hydrogen ligands) always favor
a pyramidal structure, although only slightly so for
the methyl radical. As a result, there is an increas-
ing degree of pyramidalization if we go from CH3
down to SnH3 . This differs from the classical expla-
nation for the trend in AH3 geometry and inversion
barrier as sketched in the introduction of this sec-
tion. The difference in interpretation is that the main
opposing factor to pyramidalization, i.e., the reason
why 1e1 MOs go up in energy in the Walsh dia-
gram of Figure 13, is the increase in repulsive HH
h1s|1si overlap and not the loss in hnpxy |1e01 i bond- FIGURE 16. Schematic orbital interaction diagram for
ing overlap. CH3 Li (orbital energies in eV).
The preceding discussion illustrates the rele-
vance of a quantitatively accurate description of the
various terms in the interaction. It allows us, for in- (RLi)n can be conceived as salt-like aggregates of
stance, to accurately predict the geometry and the lithium cations and carbanions, bound by electro-
inversion barrier of AH3 . It also yields the correct static forces.
qualitative MO model that adequately accounts for However, a detailed reexamination within the
the often-delicate balance of cooperating and coun- framework of KohnSham MO theory leads to a
teracting mechanisms behind the observed trends. principally different conception. The CLi bond in
CH3 Li may very well be envisaged as an electron-
pair bond (K, Fig. 15), even though it is certainly
HIGHLY POLAR ELECTRON-PAIR BONDING polar. Our analysis is based on nonlocal DFT com-
putations at the BP86/TZ2P level (see ref. 113 for
In the Electron-Pair Bond section, we have an- computational details).
alyzed electron pair bonding between equal frag- We begin with the CLi bonding mechanism in
ments. But what happens if the difference in elec- monomeric methyllithium. This can be analyzed in
tronegativity becomes very large, such as in the two ways: (1) homolytically, as an interaction be-
strongly polar carbonlithium bond.113 The picture tween CH3 and Licorresponding to a covalent
of the bond that has evolved from numerous ab initio view (K)and (2) heterolytically, as an interaction
studies is that of a predominantly ionic interac- between CH3 and Li+ this corresponds to an
tion (J, Fig. 15).196 198 This is mainly based on the ionic view (J). It is instructive to compare the two
results of a number of advanced schemes for com- approaches. We begin with the homolytic approach.
puting atomic charges. We mention the natural pop- The methyl 2a1 -SOMO and the 3 eV higher en-
ulation analysis of Weinhold,199 201 the integrated ergy lithium 2s-SOMO have a sizeable overlap of
projected population (IPP) of Streitwieser,202 and 0.33, and enter into a strongly polar electron pair
the atoms-in-molecules (AIM) approach of Bader.203 bond (see Fig. 16) with substantial 2a1 2s mix-
They yield strongly positive lithium atomic charges ing. The latter, together with a sizeable admixture
ranging from +0.75 to +0.90 electrons.204 208 Based of lithium 2pz , is responsible for the covalent char-
on these numbers it was concluded that the CLi acter of the CLi bond. The increased population
bond is ionic, and that organolithium oligomers of the methyl 2a1 of 1.40 electrons reflects polar-
ity. Note, however, that the lithium 2s still keeps
a population of 0.50 and that the lithium 2pz ae-
quires a population of 0.19 electrons! Interestingly,
the associated orbital interaction of 62.9 kcal/mol
is twice as large as the classical electrostatic inter-
action 1Velst of only 32.1 kcal/mol. The overall
homolytic bond enthalpy (for 298 K) amounts to
43.7 kcal/mol.
FIGURE 15. Ionic (J) and covalent (K) representation One might also analyze the interaction from
of the CLi bond. an ionic or heterolytic point of view (J). Not
956 VOL. 22, NO. 9

CHEMISTRY WITH ADF
unexpectedly, due to the charge separation that

occurs, this leads to an enormous increase in
the electrostatic interaction 1Velst from 32.1 to
198.9 kcal/mol. However, this choice also im-
plies that one pays a high price in terms of en-
ergy. The fragments have to be electronically ex-
cited by transferring an electron from lithium to
methyl. Subsequently, more than half of this elec-
tron is transferred back to lithium through donor
acceptor interaction of 15 kcal/mol between the
methyl 2a1 (now the lone-pair on CH3 ) and the
lithium 2s (LUMO of Li+ ). After all, this trans-
lates into the heterolytic dissociation (1Ehetero =
174.2 kcal/mol) being four times more endother-
mic, i.e., less favorable, than the homolytic disso-
ciation (1E = 45.5 kcal/mol). Furthermore, the
heterolytic approach provides an unbalanced start-
ing point for the oligomers: excessive electrostatic
repulsion, especially between lithium cations (e.g.,
+843 kcal/mol for Li4 4+ ), have to be compensated
by even larger donor/acceptor and electrostatic in- FIGURE 17. Schematic orbital interaction diagram for
teractions between the methyl anion cage and the (CH3 Li)4 (orbital energies in eV).
lithium cation cluster. This does not provide a trans-
parent description of CLi bonding, and it also
complicates understanding the cohesion within the
central lithium cluster. On the basis of these and effect is further reinforced by a strong 2s2p mixing.
other arguments113 we conclude that a homolytic This has an important consequence for the CLi
approach is in many respects the more natural electron pair bond 2a1 + 1a1 occurring in the A1
one. symmetry. The energy of the tetralithium 1a1 SOMO
Next, we discuss the CLi bond in the methyl- (6.2 eV) is so much lowered with respect to that of
lithium tetramer solely in terms of the homolytic the isolated-lithium 2s AO (3.2 eV) that it begins
approach. The computed overall bond energy 1E to approach the energy of the 2a1 partner SOMO
for the formation of (CH3 Li)4 from 4 CH3 + 4 Li (7.0 eV) on the more electronegative carbon frag-
is 308.6 kcal/mol and the methyllithium tetramer-
ment (Fig. 17). The electronegativity of lithium is
ization energy is 126.6 kcal/mol. The carbon
effectively increased in the cluster as regards the
lithium bond in (CH3 Li)4 is formed between an
interactions in the A1 symmetry. This leads to a
outer tetramethyl cage, (CH3 )4 with dCC = 3.6 ,
virtually covalent electron-pair bond in which the
and an inner tetralithium cluster, (Li)4 with dLi,Li =
respective partner SOMOs more or less keep their
2.4 (see Fig. 17). These fragments have quin-
original population of one electron! The 2a1 + 1a1
tet electronic valence states with one electron in
a CC (or LiLi) bonding orbital of A1 symme- electron-pair bond (1EA1 = 85.8 kcal/mol) is also
try and three electrons in a set of triply degenerate stronger than that in the methyllithium monomer
CC (or LiLi) antibonding orbitals of T2 symme- (1EA1 = 62.9 kcal/mol), among others, because
try. These frontier orbitals stem from the interacting of a more efficient overlap (0.55 compared to 0.33).
SOMOs of the respective monomeric fragments, i.e., The transfer of electrons out of the lithium clus-
methyl 2a1 and lithium 2s (Fig. 16). The tetram- ter into the methyl cage is caused nearly exclu-
ethyl fragment has a relatively small energy gap of sively by the three more polar 3t2 + 1t2 electron-pair
1.2 eV between the CC bonding 2a1 and the CC bonds in T2 symmetry, which together contribute
antibonding 3t2 SOMOs because the 4 methyl 2a1 387.0 kcal/mol to the carbonlithium interaction.
SOMOs are quite far away from each other. In the in- Note that this bonding mechanism does not only
ner cluster of lithium atoms, an interesting phenom- provide CLi bonding, but also increases the cohe-
enon occurs. Due to the enormous overlap of 0.65 sion within the Li4 cluster! This is because there is
between the Li 2s orbitals, the bonding 1a1 combi- a selective depopulation of the LiLi antibonding
nation SOMOlow is significantly stabilized, and this 1t2 orbitals (which hold a population of only 0.65

TE VELDE ET AL.
electrons each) while the LiLi bonding 1a1 orbital

keeps its population. Chemical Reactivity
Altogether, the KohnSham MO analysis clearly
shows that the polarity of the carbonlithium bond THE ACTIVATION STRAIN MODEL
decreases on oligomerization. This conclusion fully
agrees with the trends we obtain from VDD and In the previous sections we have focussed on
Hirshfeld analyses of the charge density distribu- bonding and structure. However, ADF is also highly
tion for the methyllithium monomer, dimer (left suitable for tackling chemical reactions. The recently
out in the discussion of the electronic structures) introduced30 activation strain (or activation strain
and tetramer. The VDD lithium charges are rel- TS interaction, ATS) model of chemical reactivity
atively small and they decrease from +0.38 via enables a transparent description and a detailed
+0.26 to +0.13 electrons along CH3 Li, (CH3 Li)2 , and understanding of activation barriers and relative ef-
(CH3 Li)4 , showing that the shift of electron density ficiencies of competing reaction mechanisms, and
how they may be affected by modifying the re-
from lithium to methyl decreases upon oligomeriza-
actants or reaction conditions (e.g., solvation). Ba-
tion. Whereas the values of the Hirshfeld lithium
sically, the activation strainTS interaction (ATS)
charges are somewhat larger, they display the same
model evolves from a fragment approach (a key
trend, again a decrease, namely from +0.49 via
concept in ADF) to analyzing transition states and
+0.42 down to +0.30 electrons along the same
activation barriers. In the ATS model, the activa-
series. The NPA lithium charges of +0.85, +0.88
tion energy 1E 6= of a reaction is decomposed into
and +0.86 electrons (MP4SDQ/6-31+G) are signif- 6=
icantly larger.113 Strikingly, they are rather constant, two terms. The first is the activation strain 1Estrain ,
and do not monitor the reduction of Li CH3 which is the energy necessary to deform the sep-
charge transfer upon oligomerization. We consider arate reactants to the geometry they adopt in the
activated complex. The second, the TS interaction
these three methods as well as other schemes (e.g., 6=
Baders AIM) as satisfactory approaches for the 1Eint , is the interaction energy between the de-
computation of basis set independent and chemi- formed reactants in the transition state (TS). This
cally meaningful atomic charges.113, 116, 209 Yet, nu- is illustrated for an SN 2 reaction in Figure 18. The
6=
merically, the methods yield rather different val- TS interaction 1Eint between the strained reactants
ues for the atomic charges. This shows that atomic can be further decomposed, for example, into the
charges are no absolute quantities; they become classical electrostatic interaction, Pauli repulsive or-
meaningful only by comparing a trend along differ- bital interaction and bonding orbital interactions
6=
ent molecules. The label 90% ionic for the CLi (1Eint = 1Velst + 1EPauli + 1Eoi ) using the extended
bond based on a computed lithium charge of +0.9 transition state (ETS) method developed by Ziegler
electrons is, therefore, not justified. and Rauk,30, 211 in the same way as 1Eint between
The VDD method has meanwhile been extended fragments of a stable minimum energy structure
in two ways.114, 119, 210 First, it can be used to study (see the Bond Energy Analysis section). Note, how-
the change in atomic charges in polyatomic frag- ever, that the activation strain (or ATS) model of
ments caused by the chemical bond between these chemical reactivity is, in principle, completely in-
fragments, for example, if two methyl radicals com-
bine to form ethane or two DNA bases associate
to become a WatsonCrick pair. Apart from pro-
viding additional insight, this approach eliminates
artefacts caused by the so-called front-atom prob-
lem that, in principle, all methods for computing
atomic charges suffer from.119 Second, the changes
in atomic VDD charges that are caused by the
chemical interactions can be decomposed into the
contributions from different irreducible representa-
tions. This enables one, for example, to examine
in a straightforward manner the charge redistrib-
ution in the and the electron systems sepa-
rately as long as they occur in different irreducible FIGURE 18. The activation strainTS interaction (ATS)
representations.119, 210 model of chemical reactivity.
958 VOL. 22, NO. 9

CHEMISTRY WITH ADF
dependent of any model or decomposition of the

chemical bond as such.
The important point about the ATS concepts is
that they provide intuitive and computable quanti-
ties that allow for a direct comparison of structurally
unrelated and, in general, nonisolobal transition
states. Although the model is not confined to DFT (it
can equally well be applied in ab initio and semiem-
pirical approaches), its combination with the Kohn
Sham MO theory and the ADF fragment approach
leads to a particularly powerful interpretative tool.
This can be applied, for example, to optimize a cata-
lyst, that is, to minimize the activation barrier of the
rate-determining step of the catalytic mechanism.
This can be done by tuning the catalysts electronic
structure, by choosing suitable ligands, such that
FIGURE 20. ATS approach to anti- vs. syn-E2
(for a given activation strain) the TS interaction with
stereochemistry.
the substrate is maximized (see refs. 30 and 211 for
applications to organometallic reactions and homo-
geneous catalysis). But the activation strain model is
also useful for understanding and directing elemen- The conclusions of ref. 212 are confirmed by the
tary organic reactions. In the following, we discuss preliminary results of recent BP86/TZ2P calcula-
two such examples: E2 elimination and SN 2 substi- tions.
tution. The reactions proceed through abstraction of a
-proton by the base under formation of ethyl-
ANTI- VS. SYN-ELIMINATION ene, the expelled leaving group F , and the con-
jugate acid HF. The latter two form a complex
Base-induced 1,2-elimination (E2, Fig. 19 upper FHF , which makes the gas-phase reaction exother-
branch) and nucleophilic substitution (SN 2, Fig. 19 mic by 40.4 kcal/mol. The syn-E2 transition state
lower branch) are two fundamental types of chem- (1E 6= = 0.5 kcal/mol) is 9 kcal/mol higher in
ical reactions. In principle, E2 elimination is always energy than that for anti-E2 elimination (1E 6= =
in competition with SN 2 substitution, and the two 9.5 kcal/mol); this is represented by the two lower
pathways may occur as unwanted side reactions of energy pathways in Figure 20. A traditional and
each other. A thorough understanding of what de- plausible explanation for this has been the increased
termines these processes is important for the design repulsion associated with the eclipsed conforma-
of efficient syntheses. tion of the substrate. But, this argument turns out
In this section, we focus on the stereochemistry to be erroneous. In the first place, eclipsed fluo-
of E2 eliminations, and wonder why they pro- roethane is only 2.3 kcal/mol higher in energy than
ceed preferentially via a coplanar H C C L the staggered conformer is; this does not account for
arrangement of the substrate, and why the anti- the difference of 9 kcal/mol in activation energies.
elimination generally prevails over syn-elimination. But, more importantly, it is the difference in energy
We consider here the anti-E2 and syn-E2 pathways between the deformed substrates in the respective
of the simple model reaction between a fluoride transition states that counts in this context, and not
anion as base and fluoroethane as substrate (F + the relative energy of the staggered and eclipsed
C2 H5 F, see Fig. 20). This has been investigated at the conformers of the isolated fluoroethane reactant.
BVWN-Stol1/TZ2P//X/DZP level of theory.30, 212 This brings us naturally to the concept of activa-
6=
tion strain (1Estrain). The syn-elimination reaches the
TS at a point where the C2 H5 F substrate is much less
deformed than in the anti-E2 transition state. This
is simply due to the fact that in the syn-E2 TS the
conjugate acid and the leaving group are located
directly in front of each other. A slight internal ro-
tation and a moderate elongation of the C H
FIGURE 19. E2 elimination and SN 2 substitution. and C F bonds (by 39 and 11%, respectively, at

TE VELDE ET AL.
where the base is a selective catalyst that favors the

high-energy process (Fig. 20). However, one restric-
tion applies to this analogy: the base is not recovered
after elimination has taken place; the catalyst is poi-
soned, so to say, with a product molecule (HF), and
is needed in stoichiometric amounts. The tendency
to adopt a coplanar arrangement of the C H and
C F bond, either anti or syn, can be understood
FIGURE 21. Schematic representation of the LUMO on the basis of the CC -bonding character of the
of the CH3 CH2 F substrate in the trans (L) and cis (M)
substrate LUMO. In that case, the overlap is op-
conformation that lead to anti- and cis-E2 elimination,
respectively. timal leading to a low LUMO energy and thus a
strong TS interaction with the base HOMO.
X/DZP) enable the abstracted proton and the leav- SUBSTITUTION VS. ELIMINATION
ing group to interact and stabilize the system. The
energy does not rise any further, and the transition In this section, we examine how base-induc-
state is reached. Such a stabilizing interaction is not ed anti-elimination competes with nucleophilic SN 2
possible in an early stage of the anti-elimination, substitution and the influence of (micro)solva-
and the anti-E2 TS is reached somewhat later, af- tion.30, 212, 215, 216 Figure 22 shows the respective re-
ter a much stronger elongation of the C H and action energy profiles for our model system F +
C F bonds (by 67 and 51%). When we consider C2 H5 F.
the activation strain in the substrate, i.e., the de- Starting from the same reactant complex (F ,
formation energy associated with bringing C2 H5 F C2 H5 F) from which E2 elimination may proceed,
into the geometry of the TS, the syn-elimination with the base can also move toward the backside of
83 kcal/mol should be the more favorable process the -methyl group and perform a symmetric
and not the anti-elimination with 126 kcal/mol. The and, therefore, thermoneutral nucleophilic substi-
two higher energy pathways in Figure 20 reflect this tution. The activation energy for SN 2 substitution
(UTS refers to uncatalyzed TS, i.e., with the base (0.5 kcal/mol) is again higher than that for E2
removed). The activation strain may be conceived elimination (9.5 kcal/mol, see Table I), in line
as a good estimate of the barrier for the correspond- with the experimental observation that (anionic)
ing thermal elimination (i.e., without the assistance eliminations strongly dominate substitution in the
of a base) of HF from fluoroethane which occurs in- gas phase.217 224 On the other hand, it is well
deed in a syn-fashion.213, 214 known that in the condensed phase, this relative
So it is the TS interaction 1E#int with the base that order is reversed or, at least, shifted in favor of
leads to the preference for anti-elimination. To un- substitution,225 227 and that the absolute rate of both
derstand this, we take a closer look at the substrates processes decreases drastically, by up to 20 orders of
LUMO, the acceptor orbital in the interaction with magnitude.228
the HOMO of the base. It is the -bonding combi-
nation between the lowest unoccupied ex orbitals
of the CH3 (left side of L and M in Fig. 21) and
the CH2 F fragments (right side of L and M in
Fig. 21) that has CH and CF antibonding char-
acter (see also the Role of the Pauli Steric Repulsion
section).193
Therefore, the antibonding character of the
LUMO is reduced and it decreases in energy, as the
C H and C F bonds elongate. This is the key
to the answer of our question. The LUMO of the
strongly deformed substrate anti-TS is much lower
in energy, by 2.2 eV, than that in the syn-TS. It
is, therefore, closer in energy to the HOMO of the
base and enters into a significantly more stabilizing
donoracceptor interaction.30, 212 This may be con- FIGURE 22. Reaction energy profiles for the SN 2 vs.
ceived as a textbook example of a catalytic process E2 competition in the gas phase.
960 VOL. 22, NO. 9

CHEMISTRY WITH ADF
TABLE I. new C F bond is already being formed to quite an

Activation Energies for E2 and SN 2 Reactions extent (dCF = 1.826 ), whereas the old one is elon-
of F + C2 H5 F.a gated only moderately (by 26%). Accordingly, the
activation strain of the substrate in TS(SN 2) is rather
Reaction E2b SN 2
low, approximately 31 kcal/mol (see Fig. 23). Yet,
the process is energetically unfavorable because of
Gas phase 9.5 0.5
Monosolvated basec 31.6 23.0
the significantly weaker and less stabilizing (donor
Disolvated basec 54.4 34.4 acceptor) TS interaction between the HOMO of the
F nucleophile and the high-energy LUMO of the
a In kcal/mol. BVWN-Stoll/TZ2P//X /DZP. little deformed substrate.
b Anti-elimination. Upon solvation, however, the F HOMO is stabi-
c With hydrogen fluoride model solvent molecules. lized (both by electrostatic and charge transfer inter-
actions with the solvent molecules) resulting in an
increase of the HOMOLUMO gap and, therefore,
a decrease of the TS interaction with the substrate
We have tried to model these solvent effects LUMO. Whereas this HOMOLUMO TS interaction
theoretically by microsolvation, that is, by intro- as such remains still the strongest for the E2 transi-
ducing only one or two HF molecules as models tion state, the weakening in this stabilizing interac-
for protic solvent molecules. Many solvation modes tion on solvation of the base is, in absolute terms,
are conceivable.216 Here, we examine the effect of the largest for this TS (Fig. 23). The high activation
mono- and disolvation of the F base. Surprisingly, strain, a characteristic of the elimination process,
already the introduction of one single HF solvent makes the barrier of the solvated E2 elimination
molecule causes a number of striking changes in the higher than that of the solvated SN 2 substitution
reaction energy profile: (1) The pronounced wells of (Fig. 23).
the reactant and product complexes become very We stress the equivalence of this view on solva-
shallow, and (2) the activation energies that were tion effects and the one based on differential sol-
negative are now becoming positive (Table I). The vation. The ATS approach analyzes the interaction
trend continues as the second solvent molecule is between base or nucleophile and the deformed sub-
brought into play. This can be understood in terms strate instead of that between solvent and reaction
of differential solvation of reactantswhich are system. The advantage is that the same reasoning
strongly stabilizedand transition stateswhich used for understanding the effect of solvation may
are weakly stabilized. Only one or two solvent also shed some light on the effect of varying the base
molecules bring us from the double-well potential, itself. For example, assuming that other factors do
typical for the gas phase, to a nearly unimodal
reaction energy profile that is characteristic for re-
actions in the condensed phase. Of course, we are
still far removed from a real condensed-phase reac-
tion that takes place in bulk solvent. But this result
shows that some very important effects occur in the
first solvation shell. Interestingly, we also find the
reversal of activation energies: on disolvation, for
example, the SN 2 barrier is 20 kcal/mol lower than
the E2 barrier (Table I).
Finally, we address the question why the intrin-
sic reactivity of F + C2 H5 F is in favor of elimi-
nation and why this solvation shifts the reactivity
toward substitution.30 Again, following the ATS ap-
proach, we begin by inspecting the activation strain
of the substrate in the anti-E2 and SN 2 transition
states. As pointed out before, the anti-E2 transi-
tion state is strongly deformed, and much energy
(126 kcal/mol) is needed to distort the substrate to FIGURE 23. Activation strain (1E#strain ) and
the associated geometry. The SN 2 transition state, on TS interaction 1E#int in E2 and SN 2 transition states
the other hand, is much less deformed and tighter: a for an unsolvated (gas) and a solvated (solv) base.

TE VELDE ET AL.
not change and leaving out overlap considerations, are long-long-short (2.91, 2.95, and 2.86 , 3d). It
it is easily seen that a stronger basewith a HOMO appears that the disagreement is not due to an in-
at high energytends to E2 elimination, whereas adequacy of the quantum chemical method, in our
a weaker basewith a HOMO at low energy case DFT at BP86 with the TZ2P basis (for which the
inclines toward SN 2 substitution. BSSE is less than 1 kcal/mol). Instead, it is caused by
a deficiency of the model systems used so far in the-
oretical computations (e.g., 3a for GC), namely, the
Structure and Bonding of DNA absence of the molecular environment (i.e., water,
sugar OH groups, counterions) that the base pairs
Very recently, we have investigated the struc- experience in the crystals studied experimentally.
ture and nature of the deoxyribonucleic acid (DNA) If we incorporate the major elements of this envi-
molecule, making use of the linearization and par-
ronment into our model, simulating them by wa-
allelization of the ADF code.119, 210, 229 The purpose
ter molecules and a Na+ ion, we achieve excellent
has been to provide a foundation for an accurate,
agreement with experiment at BP86/TZ2P (com-
quantitative description of the structure and en-
pare, for example, our 3c with the experimental 3d).
ergetics of DNA, the influence of the molecular
Note that our best simulation of the environment in-
environment (such as, e.g., water molecules and
cludes all close contacts of the GC base pair with
counterions) and, in particular, a better understand-
crystal water, sugar OH groups and counterions in
ing of the nature and behavior of this molecule that
the crystal. One may expect further improvements
is of crucial importance for the existence of all life.
if a second or even higher order environment shells
One of the key results (ref. 229, and references
are introduced but the effects are likely to be much
therein) is that we have unravelled a hitherto
less pronounced than those of the first shell. On the
unresolved discrepancy between theoretical230, 231
other hand, whether plain nucleic bases (3a) or more
and experimental232 234 hydrogen bond lengths in
realistic models for the nucleotides (3b) are used is
WatsonCrick base pairs. At first instance, we also
much less important. This finding has an enormous
arrive at the striking disagreement with experiment
scope because it shows that the already existing
as shown in Figure 24 for the three hydrogen bonds
density functionals are capable of adequately de-
in GC, for which we find a bond length pattern that
scribing biological molecules containing hydrogen
is short-long-long (2.73, 2.88, and 2.87 , 3a) at sig-
bonds.
nificant variance with the experimental values that
TDDFT: Absorption Spectrum

of Porphin
The electronic absorption spectrum of the free
base porphin (FBP) molecule has recently been
studied by various ab initio methods. We have re-
cently completed a similar study using the ADF-
RESPONSE module235 with satisfactory results,
supporting the conclusion from a CASPT2 study
(see ref. 236 and references therein). Our results for
the most important low-lying excited states of this
molecule are compared to experiment237 in Table II.
The FBP molecule is interesting because there
is an important role for interconfiguration mixing
in the low-lying excitations, first described success-
fully at the semiempirical level by Gouterman238
FIGURE 24. Hydrogen bond distances (in ) in plain
guaninecytosine (GC, 3a), GC with inclusion of the with his famous four-orbital model. This leads to the
backbone (3b) and one of our GC model systems so-called B and Q bands in the FBP molecule as well
containing water molecules and a sodium as in other porphyrins. One can indeed observe the
counterion (3c) from BP86/TZ2P (3a, 3c) and expected mixing by looking at the solution vectors
BP86/DZP (3b) computations and in the crystal of of the excitation energy eigenvalue equation. In this
sodium guanylyl-30 ,50 -cytidine nonahydrate from X-ray way one sees the weight of various orbital replace-
diffraction (3d).234 ments to a certain excitation. This can, furthermore,
962 VOL. 22, NO. 9

CHEMISTRY WITH ADF
TABLE II.
TDDFT (BP86/ALDA) and Experimental Excitation Energies (eV) and Oscillator Strengths (f) for Free Base
Porphin (FBP).
TDDFTa TDDFTa Experimentb Experimentb

Excitation Oscillator Excitation Oscillator
State Energy (eV) Strength f Energy Strength
1 1 B3u 2.16 0.01 1.982.02 (Qx ) 0.01

1 1 B2u 2.29 0.0005 2.332.42 (Qy ) 0.06
2 1 B3u 3.01 0.04 3.133.33 (B) 1.15
2 1 B2u 2.98 0.1338
3 1 B2u 3.41 0.8962 3.65 (N) <0.1
3 1 B3u 3.47 0.7293
4 1 B3u 3.77 0.1688 4.25 (L) 0.1
4 1 B2u 3.76 0.1272
a Ref. 235.
b Ref. 237 and other references in ref. 235.
be used to study differences between related sys- ferences. The research on the porphyrin spectra
tems. The effect of adding a metal in the centre is ongoing, and also includes the nonlinar optical
of the porphyrin or the effect of using a different properties.241
macrocycle immediately shows up in the positions Currently a fragment orbital analysis is being im-
of the excitation energies, the magnitude of the os- plemented to study the adsorption spectra of these
cillator strengths, and, perhaps most interestingly, and molecules in even greater detail.
in the weights of the solution vectors. From such
information, one can relate the positions and mag-
nitudes of peaks in the absorption spectrum to the
geometrical and electronic structure of the molecule. Outlook
In this way, the ADF-RESPONSE module can help
to understand differences between spectra of vari-
As we have discussed, ADF is a package for de-
ous molecules.
tailed and quantitative analysis of reactivity, struc-
An important feature of the TDDFT results in
Table II is the paired structure of the 1 B2u and the tures, and various properties in the gas phase and in
1
B3u excitations. All low-lying experimental bands solution and for complicated and large-scale calcu-
are assigned to a pair of 1 B2u and 1 B3u excitations, lations on a large variety of systems. Being primarily
which is in full agreement with the CASPT2 re- academic software, ADF will continue to be devel-
sults, and challenges previous coupled cluster and oped as the result of ongoing scientific research.
CI calculations.239, 240 Focus and emphasis in code development naturally
Although the TDDFT excitation energies are a adjust to trends in scientific interests. This will target
bit low with respect to experiment, the assignments both functionality and performance. On a short-
are unambiguous. In the experimental spectrum, term further work on XC potentials and energy
the B and N bands strongly overlap, so that one functionals, extension of the linear-scaling technol-
usually speaks of the BN band system in FBP. ogy and improved parallelization behavior are at
The TDDFT results give large oscillator strengths hand. With the systematic and rapid improvement
in the N band. For metal porphyrins on the other of both hardware and software the applicability of
hand, the TDDFT results agree with the experi- ADF will cover more and more fields and increas-
mental result that the B band is clearly separated ingly more complex problems in chemistry, biology,
from and much stronger than the N band. There- and physics.
fore, the role that the B and N bands play in FBP Online information about ADF, including user
differs substantially from the situation in closely re- manuals, other documentation, and information
lated systems. Analysis of the TDDFT results can about licenses and distribution, can be found at the
be a very useful aid in understanding such dif- internet address http://www.scm.com.

TE VELDE ET AL.
28. Gritsenko, O. V.; Schipper, P. R. T.; Baerends, E. J. J Chem

References Phys 1997, 107, 5007.
29. Stowasser, R.; Hoffmann, R. J Am Chem Soc 1999, 121,
1. Parr, R. G.; Yang, W. Density-Functional Theory of Atoms 3414.
and Molecules; Oxford University Press: New York, 30. Bickelhaupt, F. M. J Comput Chem 1999, 20, 114.
1989.
31. Bickelhaupt, F. M.; Baerends, E. J. In Reviews of Compu-
2. Dreizler, R. M.; Gross, E. K. U. Density Functional The- tational Chemistry; Boyd, D. B.; Lipkowitz, K. B., Eds.;
ory, an Approach to the Quantum Many-Body Problem; Wiley-VCH: New York, 2000, p. 1, vol. 15.
Springer Verlag: Berlin, 1990.
32. Redfern, P. C.; Blaudeau, J.-P.; Curtiss, L. A. J Phys Chem
3. Ellis, D. E., Ed. Electronic Density Functional Theory of A 1997, 101, 8701.
Molecules, Clusters and Solids; Kluwer Academic Publish-
33. Curtiss, L. A.; Redfern, P. C.; Raghavachari, K.; Pople, J. A.
ers: Dordrecht, 1995.
J Chem Phys 1998, 109, 42.
4. Dunitz, J. D.; Hafner, K.; Houk, K. N.; Ito, S.; Lehn, J.-M.;
34. van Leeuwen, R.; Baerends, E. J. Phys Rev A 1994, 49, 2421.
Raymond, K. N.; Rees, C. W.; Thiem, J.; Vgtle, F., Eds. Top-
ics in Current Chemistry; Springer: Berlin, 1996, vol. 183. 35. Proynov, E. I.; Sirois, S.; Salahub, D. R. Int J Quantum
Chem 1997, 64, 427.
5. Seminario, J. M., Ed. Recent Developments and Appli-
cations of Modern Density Functional Theory; Elsevier: 36. Schipper, P. R. T.; Gritsenko, O. V.; Baerends, E. J. Phys
Amsterdam, 1996. Rev A 1998, 57, 1729.
6. Baerends, E. J.; Ellis, D. E.; Ros, P. Chem Phys 1973, 2, 37. van Voorhis, T.; Scuseria, G. E. J Chem Phys 1998, 109, 400.
41. 38. Filatov, M.; Thiel, W. Phys Rev A 1998, 57, 189.
7. Sambe, J.; Felton, R. H. J Chem Phys 1975, 62, 1122. 39. Krieger, J. B.; Chen, J.; Iafrate, G. J.; Savin, A. In Electron
8. Ellis, D. E.; Adachi, H.; Averill, F. W. Surf Sci 1976, 58, Correlations and Materials Properties; Gonis, A.; Kioussis,
497. N., Eds.; Plenum: New York, 1999.
9. Dunlap, B. I.; Corolly, J. W. P.; Sabin, J. R. J Chem Phys 1979, 40. Perdew, J. P.; Kurth, S.; Zupan, A.; Blaha, P. Phys Rev Lett
71, 3396. 1999, 82, 2544.
10. Delley, B.; Ellis, D. E. J Chem Phys 1982, 76, 1949. 41. Perdew, J. P.; Kurth, S.; Zupan, A.; Blaha, P. Phys Rev Lett
1999, 82, 5179.
11. (a) Boerrigter, P. M.; te Velde, G.; Baerends, E. J. Int J Quan-
tum Chem 1988, 33, 87; (b) te Velde, G.; Baerends, E. J. 42. Tsuneda, T.; Suzumura, T.; Hirao, K. J Chem Phys 1999,
J Comput Phys 1992, 99, 84. 110, 10664.
12. te Velde, G.; Baerends, E. J. Phys Rev B 1991, 44, 7888. 43. Becke, A. D. J Chem Phys 2000, 112, 4020.
13. Fonseca Guerra, C.; Visser, O.; Snijders, J. G.; te Velde, G.; 44. Boese, A. D.; Doltsinis, N. L.; Handy, N. C.; Sprik, M.
Baerends, E. J. In Methods and Techniques in Computa- J Chem Phys 2000, 112, 1670.
tional Chemistry; Clementi, E.; Corongiu, G., Eds.; Stef: 45. Proynov, E. I.; Chermette, H.; Salahub, D. R. J Chem Phys
Cagliari, 1995, p. 305. 2000, 113, 10013.
14. Andzelm, J. W.; Nguyen, D. T.; Eggenberger, R.; Salahub, 46. Boese, A. D.; Doltsinis, N. L.; Handy, N. C.; Sprik, M.
D. R.; Hagler, A. T. Comput Chem 1995, 19, 145. J Chem Phys 2000, 112, 1670.
15. Fonseca Guerra, C.; Snijders, J. G.; te Velde, G.; Baerends, 47. Rohklin, V. J Comput Phys 1985, 60, 187.
E. J. Theor Chem Acc 1998, 99, 391. 48. Greengard, L.; Rohklin, V. J Comput Phys 1987, 73, 325.
16. Ziegler, T. Chem Rev 1991, 91, 651. 49. White, C. A.; Johnson, B. G.; Gill, P. M. W.; Head-Gordon,
17. Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; M. Chem Phys Lett 1994, 230, 8.
Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys Rev B 1992, 50. White, C. A.; Head-Gordon, M. J Chem Phys 1994, 101,
46, 6671. 6593.
18. Becke, A. D. J Chem Phys 1992, 96, 2155. 51. Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Science 1995, 271,
19. Becke, A. D. J Chem Phys 1993, 98, 1372. 51.
20. Becke, A. D. J Chem Phys 1993, 98, 5648. 52. Kutteh, R.; Apra, E.; Nichols, J. Chem Phys Lett 1995, 238,
173.
21. Becke, A. D. J Chem Phys 1993, 97, 9173.
22. Ziegler, T. Can J Chem 1995, 73, 743. 53. Adamson, R. D.; Dombroski, J. P.; Taylor, S. W.; Gill,
P. M. W. J Phys Chem 1996, 100, 6272.
23. Ziegler, T.; Rauk, A. Inorg Chem 1979, 18, 1558.
54. Challacombe, M.; Schwegler, E.; Almlf, J. J Chem Phys
24. Ziegler, T.; Rauk, A. Inorg Chem 1979, 18, 1755. 1996, 104, 4685.
25. Baerends, E. J.; Gritsenko, O. V.; van Leeuwen, R. In Chem- 55. Gill, P. M. W. Chem Phys Lett 1997, 270, 193.
ical Applications in Density Functional Theory; Laird,
56. Schwegler, E.; Challacombe, M.; Head-Gordon, M. J Chem
B. B.; Ross, R. B.; Ziegler, T., Eds.; American Chemical So-
Phys 1997, 106, 5526.
ciety: Washington, DC, 1996, p. 20.
57. Millam, J. M.; Scuseria, G. E. J Chem Phys 1997, 106, 5569.
26. Baerends, E. J.; Gritsenko, O. V.; van Leeuwen, R. In New
Methods in Quantum Theory; Tsipis, C. A.; Popov, V. S.; 58. Daniels, A. D.; Millam, J. M.; Scuseria, G. E. J Chem Phys
Herschbach, D. R.; Avery, J. A., Eds.; Kluwer Academic 1997, 107, 425.
Publishers: Dordrecht, 1996, vol. 8. 59. Stratman, R. E.; Scuseria, G. E.; Frisch, M. J. Chem Phys
27. Baerends, E. J.; Gritsenko, O. V. J Phys Chem A 1997, 101, Lett 1997, 257, 213.
5383. 60. Kudin, K. N.; Scuseria, G. E. Chem Phys Lett 1998, 283, 61.
964 VOL. 22, NO. 9

CHEMISTRY WITH ADF
61. Kohn, W.; Sham, L. Phys Rev A 1965, 140, 1133. 96. Burant, J. C.; Strain, M. C.; Scuseria, G. E.; Frisch, M. J.
62. Baerends, E. J.; Ros, P. Int J Quantum Chem 1978, S12, Chem Phys Lett 1996, 248, 43.
169. 97. Warshel, A.; Levitt, M. J Mol Biol 1976, 103, 227.
63. Snijders, J. G.; Baerends, E. J. Mol Phys 1978, 36, 1789. 98. Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz,
64. Snijders, J. G.; Baerends, E. J.; Ros, P. Mol Phys 1979, 38, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Cold-
1909. well, J. W.; Kollman, P. A. J Am Chem Soc 1995, 117, 5179.
65. van Lenthe, E.; Baerends, E. J.; Snijders, J. G. J Chem Phys 99. Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Skiff, W. M.
1993, 99, 4597. J Am Chem Soc 1992, 114, 10024.
66. van Lenthe, E.; Baerends, E. J.; Snijders, J. G. J Chem Phys 100. Singh, U. C.; Kollman, P. A. J Comput Chem 1986, 7, 718.
1994, 101, 9783. 101. Maseras, F.; Morokuma, K. J Comput Chem 1995, 16, 1170.
67. van Gisbergen, S. J. A.; Snijders, J. G.; Baerends, E. J. 102. Woo, T. K.; Margl, P. M.; Blchl, P. E.; Ziegler, T. J Phys
J Chem Phys 1995, 103, 9347. Chem B 1997, 101, 7877.
68. van Lenthe, E.; Snijders, J. G.; Baerends, E. J. J Chem Phys 103. Cavallo, L.; Woo, T. K.; Ziegler, T. Can J Chem 1998, 76,
1996, 105, 6505. 1457.
69. van Lenthe, E.; van Leeuwen, R.; Baerends, E. J.; Snijders, 104. Deng, L.; Woo, T. K.; Cavallo, L.; Margl, P.; Ziegler, T. J Am
J. G. Int J Quantum Chem 1996, 57, 281. Chem Soc 1997, 119, 6177.
70. van Gisbergen, S. J. A.; Osinga, V. P.; Gritsenko, O. V.; van 105. Deng, L.; Ziegler, T.; Woo, T. K.; Margl, P.; Fan, L.
Leeuwen, R.; Snijders, J. G.; Baerends, E. J. J Chem Phys Organometallics 1998, 17, 3240.
1996, 105, 3142. 106. van Gisbergen, S. J. A.; Fonseca Guerra, C.; Baerends, E. J.
71. van Gisbergen, S. J. A.; Snijders, J. G.; Baerends, E. J. J Comput Chem 2000, 21, 1511.
J Chem Phys Lett 1996, 259, 599. 107. Hoffmann, R. Angew Chem Int Ed 1982, 21, 711.
72. van Gisbergen, S. J. A.; Snijders, J. G.; Baerends, E. J. 108. Baerends, E. J.; Branchadell, V.; Sodupe, M. Chem Phys
J Chem Phys 1998, 109, 10644. Lett 1997, 265, 481.
73. van Gisbergen, S. J. A.; Snijders, J. G.; Baerends, E. J. 109. Morokuma, K. J Chem Phys 1971, 55, 1236.
J Chem Phys 1998, 109, 10657.
110. Kitaura, K.; Morokuma, K. Int J Quantum Chem 1976, 10,
74. Versluis, L.; Ziegler, T. J Chem Phys 1988, 322, 88. 325.
75. Fan, L.; Ziegler, T. J Chem Phys 1992, 96, 9005. 111. Baerends, E. J. In Cluster Models for Surface and Bulk
76. Deng, L.; Ziegler, T. Int J Quantum Chem 1994, 52, 731. Phenomena; Nato Asi Series B; Pacchioni, G.; Bagus, P. S.;
77. Schreckenbach, G.; Ziegler, T. J Phys Chem 1995, 99, 606. Parmigiani, F., Eds.; 1992, p. 189, vol. 283.
78. Schreckenbach, G.; Ziegler, T. Int J Quantum Chem 1996, 112. Ziegler, T.; Rauk, A. Theoret Chim Acta 1977, 46, 1.
60, 753. 113. Bickelhaupt, F. M.; van Eikema Hommes, N. J. R.; Fonseca
79. Schreckenbach, G. Ph.D. thesis, University of Calgary, Cal- Guerra, C.; Baerends, E. J. Organometallics 1996, 15, 2923.
gary, 1996. 114. Bickelhaupt, F. M.; Fonseca Guerra, C.; Handgraaf, J. W.;
80. Schreckenbach, G.; Ziegler, T. Int J Quantum Chem 1997, Baerends, E. J., in preparation.
61, 899. 115. Hirshfeld, F. L. Theoret Chim Acta 1977, 44, 129.
81. Wolff, S. K.; Ziegler, T. J Chem Phys 1998, 109, 895. 116. Wiberg, K. B.; Rablen, R. B. J Comput Chem 1993, 14, 1504.
82. Wolff, S. K.; Ziegler, T.; van Lenthe, E.; Baerends, E. J. 117. Mulliken, R. S. J Chem Phys 1955, 23, 1833.
J Chem Phys 1999, 110, 7689. 118. Mulliken, R. S. J Chem Phys 1955, 23, 2343.
83. Pye, C. C.; Ziegler, T. Theor Chem Acc 1999, 101, 396. 119. Fonseca Guerra, C.; Bickelhaupt, F. M.; Snijders, J. G.;
84. Woo, T. K.; Cavallo, L.; Ziegler, T. Theor Chem Acc 1998, Baerends, E. J. Chem Eur J 1999, 5, 3581.
100, 307. 120. Ceperley, D. M.; Alder, B. J. Phys Rev Lett 1980, 45, 566.
85. Woo, T. K.; Patchkovskii, S.; Ziegler, T. Comput Sci Eng 121. Vosko, S. H.; Wilk, L.; Nusair, M. Can J Phys 1980, 58,
2000, 2, 28. 1200.
86. Philipsen, P. H. T.; van Lenthe, E.; Snijders, J. G.; Baerends, 122. Stoll, H.; Pavlidou, C. M. E.; Preuss, H. Theoret Chim Acta
E. J. Phys Rev B 1997, 56, 13556. 1978, 49, 143.
87. Olsen, R. A.; Kroes, G. J.; Baerends, E. J. J Chem Phys 1999, 123. Perdew, J. P.; Wang, Y. Phys Rev B 1986, 33, 8800.
111, 11155.
124. Perdew, J. P. Phys Rev B 1986, 33, 8822.
88. Philipsen, P. H. T.; Baerends, E. J. Phys Rev B 2000, 61, 1773.
125. Becke, A. D. Phys Rev A 1988, 38, 3098.
89. te Velde, G. Ph.D. thesis, Vrije Universiteit, Amsterdam,
126. Lee, C.; Yang, W.; Parr, R. G. Phys Rev B 1988, 37, 785.
1990.
127. Gritsenko, O. V.; Schipper, P. R. T.; Baerends, E. J. Chem
90. van Lenthe, E.; Ehlers, A. E.; Baerends, E. J. J Chem Phys
Phys Lett 1999, 302, 199.
1999, 110, 8943.
128. Schipper, P. R. T.; Gritsenko, O. V.; van Gisbergen, S. J. A.;
91. Raffenetti, R. C. J Chem Phys 1973, 59, 5936.
Baerends, E. J. J Comput Phys 2000, 112, 1344.
92. Pulay, P. Chem Phys Lett 1980, 73, 393.
129. Grning, M.; Gritsenko, O. V.; van Gisbergen, S. J. A.;
93. Pulay, P. J Comput Chem 1982, 3, 556. Baerends, E. J. J Chem Phys 2001, 114, 652.
94. Hamilton, T. P.; Pulay, P. J Chem Phys 1986, 84, 5728. 130. Ziegler, T.; Snijders, J. G.; Baerends, E. J. J Chem Phys 1981,
95. Fischer, T. H.; Almlf, J. J Phys Chem 1992, 96, 9768. 74, 1271.

TE VELDE ET AL.
131. DeKock, R. L.; Baerends, E. J.; Boerrigter, P. M.; Snijders, 166. Khandogin, J.; Ziegler, T. J Phys Chem A 2000, 104, 113.
J. G. Chem Phys Lett 1984, 105, 308. 167. Autschbach, J.; Ziegler, T. J Chem Phys 2000, 113, 936.
132. DeKock, R. L.; Baerends, E. J.; Boerrigter, P. M.; Hengel- 168. Schreckenbach, G.; Ziegler, T. J Phys Chem A 1997, 101,
molen, R. J Am Chem Soc 1984, 106, 33387. 3388.
133. Boerrigter, P. M. Ph.D. thesis, Vrije Universiteit, Amster- 169. van Lenthe, E.; van der Avoird, A.; Wormer, P. E. S. J Chem
dam, 1987. Phys 1997, 107, 2488.
134. Boerrigter, P. M.; Buijse, M. A.; Snijders, J. G. Chem Phys 170. van Lenthe, E.; van der Avoird, A.; Wormer, P. E. S. J Chem
1987, 111, 47. Phys 1998, 108, 4783.
135. Boerrigter, P. M.; Baerends, E. J.; Snijders, J. G. Chem Phys 171. van Lenthe, E.; van der Avoird, A.; Hagen, W. R.; Reijerse,
1988, 122, 357. E. J. J Phys Chem A 2000, 104, 2070.
136. Ziegler, T.; Tschinke, V.; Baerends, E. J.; Snijders, J. G.; 172. van Lenthe, E.; Baerends, E. J. J Chem Phys 2000, 112,
Ravenek, W. J Phys Chem 1989, 93, 3050. 8279.
137. Li, J.; Scnreckenbach, G.; Ziegler, T. J Am Chem Soc 1995, 173. Gross, E. K. U.; Dobson, J. F.; Petersilka, M. In Density
117, 486. Functional Theory; Nalewajski, R. F., Ed.; Springer: Hei-
138. Bersuker, I. B.; Budnikov, S. S.; Leizerov, B. A. Int J Quan- delberg, 1996.
tum Chem 1977, 11, 543. 174. van Gisbergen, S. J. A. Ph.D. thesis, Vrije Universiteit, Am-
139. Pascual-Ahuir, J. L.; Silla, E.; Tomasi, J.; Bonaccorsi, R. sterdam, 1998.
J Comput Chem 1987, 8, 778. 175. van Gisbergen, S. J. A.; Snijders, J. G.; Baerends, E. J. Com-
140. Klamt, A.; Schrmann, G. J Chem Soc Perkin Trans 1993, put Phys Commun 1999, 118, 119.
2, 799. 176. van Gisbergen, S. J. A.; Kootstra, F.; Schipper, P. R. T.; Grit-
141. Klamt, A. J Phys Chem 1995, 99, 2224. senko, O. V.; Snijders, J. G.; Baerends, E. J. Phys Rev A 1998,
142. Pascual-Ahuir, J. L.; Silla, E. J Comput Chem 1990, 11, 1047. 57, 2556.
177. Osinga, V. P.; van Gisbergen, S. J. A.; Snijders, J. G.;
143. Silla, E.; Tunn, L.; Pascual-Ahuir, J. L. J Comput Chem
Baerends, E. J. J Chem Phys 1997, 106, 5091.
1991, 12, 1077.
178. Casida, M. E. In Recent Advances in Density-Functional
144. Versluis, L. Ph.D. thesis, University of Calgary, 1989.
Methods; Chong, D. P., Ed.; World Scientific: Singapore,
145. Fan, L.; Ziegler, T. J Chem Phys 1991, 95, 7401. 1995, p. 155.
146. Fan, L.; Versluis, L.; Ziegler, T.; Baerends, E. J.; Ravenek, W. 179. Bauernschmitt, R.; Ahlrichs, R. Chem Phys Lett 1996, 256,
Int J Quantum Chem Symp 1988, S22, 173. 454.
147. Broyden, C. G. J Inst Math Applic 1970, 6, 76. 180. van Gisbergen, S. J. A.; Snijders, J. G.; Baerends, E. J. Phys
148. Fan, L.; Ziegler, T. J Am Chem Soc 1992, 114, 10890. Rev Lett 1997, 78, 3097.
149. Gonzalez, C.; Schlegel, H. B. J Phys Chem 1990, 94, 5523. 181. Rosa, A.; Baerends, E. J.; van Gisbergen, S. J. A.; van
150. Deng, L.; Ziegler, T. J Chem Phys 1993, 99, 3823. Lenthe, E.; Groeneveld, J. A.; Snijders, J. G. J Am Chem
151. Fukui, K. Acc Chem Res 1981, 14, 363. Soc 1999, 121, 10356.
152. Kraka, E.; Dunning, T. H., Jr. In Advances in Molecular 182. Bickelhaupt, F. M.; Nibbering, N. M. M.; van Wezenbeek,
Electronic Structure Theory; Dunning, T. H., Jr., Ed.; JAI E. M.; Baerends, E. J. J Phys Chem 1992, 96, 4864.
Press: Greenwich, CT, 1990, vol. 1. 183. Bickelhaupt, F. M.; Bickelhaupt, F. Chem Eur J 1999, 96,
153. Fan, L.; Ziegler, T. J Phys Chem 1992, 96, 6937. 4864.
184. Bickelhaupt, F. M.; Diefenbach, A.; de Visser, S. V.; de Kon-
154. Torrent, M.; Gili, P.; Duran, M.; Sola, M. J Chem Phys 1996,
ing, L. J.; Nibbering, N. M. M. J Phys Chem A 1998, 102,
104, 9499.
9549.
155. Brces, A.; Dickson, R. M.; Fan, L.; Jacobsen, H.; Swerhone,
185. Clark, T. J. J Am Chem Soc 1988, 110, 1672.
D.; Ziegler, T. Comput Phys Commun 1997, 100, 247.
186. Gill, P. M. W.; Radom, L. J. J Am Chem Soc 1988, 110,
156. Jacobsen, H.; Brces, A.; Swerhone, D. P.; Ziegler, T. Com-
4931.
put Phys Commun 1997, 100, 263.
187. de Visser, S. V.; de Koning, L. J.; Nibbering, N. M. M. Int J
157. Schreckenbach, G.; Ziegler, T. Theor Chem Acc 1998, 99, 71.
Mass Spectrom Ion Processes 1996, 157/158, 283.
158. Ruiz-Morales, Y.; Ziegler, T. J Phys Chem A 1998, 102, 3970.
188. de Visser, S. V.; Bickelhaupt, F. M.; de Koning, L. J.; Nib-
159. Ruiz-Morales, Y.; Schreckenbach, G.; Ziegler, T. J Phys bering, N. M. M. Int J Mass Spectrom Ion Processes 1998,
Chem 1997, 101, 4121. 179/180, 43.
160. Ruiz-Morales, Y.; Schreckenbach, G.; Ziegler, T. Organ- 189. James, M. A.; McKee, M. L.; Illies, A. J. J. J Am Chem Soc
ometallics 1996, 15, 3920. 1996, 118, 7836.
161. Ruiz-Morales, Y.; Schreckenbach, G.; Ziegler, T. J Phys 190. Deng, Y.; Illies, A. J.; James, M. A.; McKee, M. L.; Peschke,
Chem 1996, 100, 3359. M. J. J Am Chem Soc 1995, 117, 420.
162. Schreckenbach, G.; Ruiz-Morales, Y.; Ziegler, T. J Chem 191. Sodupe, M.; Bertran, J.; Rodrguez-Santiago, L.; Baerends,
Phys 1996, 104, 8605. E. J. J Phys Chem A 1999, 103, 166.
163. Ehlers, A. W.; Ruiz-Morales, Y.; Baerends, E. J.; Ziegler, T. 192. Noodleman, L.; Post, D.; Baerends, E. J. Chem Phys 1982,
Inorg Chem 1997, 36, 5031. 64, 159.
164. Dickson, R. M.; Ziegler, T. J Phys Chem 1996, 100, 5286. 193. Bickelhaupt, F. M.; Ziegler, T.; von Ragu Schleyer, P.
165. Khandogin, J.; Ziegler, T. Spectr Acta A 1999, 55, 607. Organometallics 1996, 15, 1477.
966 VOL. 22, NO. 9

CHEMISTRY WITH ADF
194. Albright, T. A.; Burdett, J. K.; Whangbo, M.-H., Orbital 217. Nibbering, N. M. M. Acc Chem Res 1990, 23, 279.
Interactions in Chemistry; WileyInterscience: New York, 218. Bickelhaupt, F. M.; Buisman, G. J. H.; de Koning, L. J.; Nib-
1985. bering, N. M. M.; Baerends, E. J. J Am Chem Soc 1995, 117,
195. Gimarc, B. M. Molecular Structure and Bonding; Academic 9889.
Press: New York, 1979. 219. Bickelhaupt, F. M.; de Koning, L. J.; Nibbering, N. M. M.
196. Streitwieser, A., Jr.; Bachrach, S. M.; Dorigo, A. E.; von J Org Chem 1993, 58, 2436.
Ragu Schleyer, P. In Lithium Chemistry; Sapse, A. M.; von 220. DeThuri, V. F.; Hintz, P. A.; Ervin, K. M. J Phys Chem A
Ragu Schleyer, P., Eds.; Wiley & Sons: New York, 1995. 1997, 101, 5969.
197. Lambert, C.; von Ragu Schleyer, P. In Carbanionen. Meth- 221. Li, C.; Ross, P.; Szulejko, J. E.; McMahon, T. B. J Am Chem
oden der Organischen Chemie (Houben-Weyl); Hanack, Soc 1996, 118, 9360.
M., Ed.; Thieme: Stuttgart, 1993, vol. E19d.
222. Lum, R. C.; Grabowski, J. J. J Am Chem Soc 1992, 114, 9663.
198. Lambert, C.; von Ragu Schleyer, P. Angew Chem 1994,
223. DePuy, C. H.; Gronert, S.; Mullin, A.; Bierbaum, V. M. J Am
106, 1187.
Chem Soc 1990, 112, 8650.
199. Reed, A. E.; Weinstock, R. B.; Weinhold, F. J Chem Phys
224. Jones, M. E.; Ellison, G. B. J Am Chem Soc 1989, 111, 1645.
1985, 83, 735.
225. March, J. Advanced Organic Chemistry; Wiley-Inter-
200. Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem Rev 1988, 88,
science: New York, 1997, 4th ed.
899.
226. Reichardt, C. Solvents and Solvent Effects in Organic
201. Reed, A. E.; von Ragu Schleyer, P. J Am Chem Soc 1990,
Chemistry; VCH: Weinheim, 1988, 2nd ed.
112, 1434.
202. Collins, J. B.; Streitwieser, A., Jr. J Comput Chem 1980, 1, 227. Isaacs, N. S. Physical Organic Chemistry; Longman: Har-
81. low, 1995, 2nd ed.
203. Bader, R. W. F. Atoms in Molecules, a Quantum Theory; 228. Riveros, J. M.; Jos, S. M.; Takashima, K. Adv Phys Org
Clarendon Press: Oxford, 1990. Chem 1985, 21, 197.
204. Streitwieser, A., Jr.; Williams, J. W.; Alexandratos, S.; 229. Fonseca Guerra, C.; Bickelhaupt, F. M.; Snijders, J. G.;
McKelvey, J. M. J Am Chem Soc 1976, 98, 4778. Baerends, E. J. J Am Chem Soc 2000, 122, 4117.
205. Ritchie, J. P.; Bachrach, S. M. J Am Chem Soc 1987, 109, 230. Bertran, J.; Oliva, A.; Rodrguez-Santiago, L.; Sodupe, M.
5909. J Am Chem Soc 1998, 120, 8159.
206. Hiberty, P. C.; Cooper, D. L. J Mol Struct (Theochem) 1988, 231. (a) Sponer, J.; Leszczynski, J.; Hobza, P. J Phys Chem 1996,
169, 437. 100, 1965; (b) Brameld, K.; Dasgupta, S.; Goddard, W. A.,
III. J Phys Chem 1997, 101, 4851.
207. Cioslowski, J. J Am Chem Soc 1989, 111, 8333.
232. Saenger, W. Principles of Nucleic Acid Structure; Springer
208. Kaufmann, E.; Raghavachari, K.; Reed, A. E.; von Ragu
Verlag: New York, 1984.
Schleyer, P. Organometallics 1988, 7, 1597.
233. Seeman, N. C.; Rosenberg, J. M.; Suddath, F. L.; Kirn, J. J. P.;
209. Meister, J.; Schwarz, W. H. E. J Phys Chem 1994, 98,
Rich, A. J Mol Biol 1976, 104, 109.
8245.
234. Rosenberg, J. M.; Seeman, N. C.; Day, R. O.; Rich, A. J Mol
210. Fonseca Guerra, C.; Bickelhaupt, F. M. Angew Chem 1999,
Biol 1976, 104, 145.
111, 3120; Angew Chem Int Ed 1999, 38, 2942.
211. Bickelhaupt, F. M.; Ziegler, T.; von Ragu Schleyer, P. 235. van Gisbergen, S. J. A.; Rosa, A.; Ricciardi, G.; Baerends,
Organometallics 1995, 14, 2288. E. J. J Chem Phys 1999, 111, 2499.
212. Bickelhaupt, F. M.; Baerends, E. J.; Nibbering, N. M. M.; 236. Serrano-Andrs, L.; Merchn, M.; Rubio, M.; Roos, B. O.
Ziegler, T. J Am Chem Soc 1993, 115, 9160. Chem Phys Lett 1998, 295, 195.
213. Ferguson, H. A.; Ferguson, J. D.; Holmes, B. E. J Phys Chem 237. Edwards, L.; Dolphin, D. H.; Gouterman, M.; Adler, A. D.
1998, 102, 5393. J Mol Spectr 1971, 38, 16.
214. Kato, S.; Morokuma, K. J Chem Phys 1980, 73, 3900. 238. Gouterman, M. J Mol Spectr 1961, 6, 138.
215. Bertran, J. In New Theoretical Concepts for Understand- 239. Gwaltney, S. R.; Bartlett, R. J. J Chem Phys 1998, 108, 6790.
ing Organic Reactions; Bertran, J.; Csizmadia, I. G., Eds.; 240. Tokita, Y.; Hasegawa, J.; Nakatsuji, H. J Phys Chem A 1998,
Kluwer Academic Publishers: Dordrecht, 1989, p. 231. 102, 1843.
216. Bickelhaupt, F. M.; Baerends, E. J.; Nibbering, N. M. M. 241. Ricciardi, G.; Rosa, A.; van Gisbergen, S. J. A.; Baerends,
Chem Eur J 1996, 2, 196. E. J. J Phys Chem A 2000, 104, 635.

Chemistry With ADF

Uploaded by

Copyright:

Available Formats

Chemistry With ADF

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Chemistry With ADF

Uploaded by

Copyright:

Available Formats

Chemistry with ADF

G. TE VELDE,1 F. M. BICKELHAUPT,2 E. J. BAERENDS,2

Received 16 June 2000; accepted 4 December 2000

ABSTRACT: We present the theoretical and technical foundations of the

Correspondence to: F. M. Bickelhaupt; e-mail: bickel@chem.

Journal of Computational Chemistry, Vol. 22, No. 9, 931967 (2001)

Keywords: ADF program; density functional theory; materials science; chemical

The potential VXC is the functional derivative

K ohnSham DFT1 5 has in the recent past

932 VOL. 22, NO. 9

to triple-zeta-plus-double-polarization quality (and

JOURNAL OF COMPUTATIONAL CHEMISTRY 933

934 VOL. 22, NO. 9

JOURNAL OF COMPUTATIONAL CHEMISTRY 935

936 VOL. 22, NO. 9

JOURNAL OF COMPUTATIONAL CHEMISTRY 937

FIGURE 1. Scaling behavior of the fit procedure, that

The quantum mechanicsmolecular mechanics

938 VOL. 22, NO. 9

fects from the MM part in the SCF procedure and PARALLELIZATION

JOURNAL OF COMPUTATIONAL CHEMISTRY 939

chooses the message-passing method. This ensures

The numerical integrations (matrix elements of

940 VOL. 22, NO. 9

only a transformation of the basis set. If there are

JOURNAL OF COMPUTATIONAL CHEMISTRY 941

942 VOL. 22, NO. 9

JOURNAL OF COMPUTATIONAL CHEMISTRY 943

944 VOL. 22, NO. 9

JOURNAL OF COMPUTATIONAL CHEMISTRY 945

946 VOL. 22, NO. 9

Frequencies Analytical Second Derivatives

JOURNAL OF COMPUTATIONAL CHEMISTRY 947

948 VOL. 22, NO. 9

The time-dependent KohnSham equations are equation becomes

JOURNAL OF COMPUTATIONAL CHEMISTRY 949

a general molecular multipole moment qlm to such a 1 X j

950 VOL. 22, NO. 9

Derivatives w.r.t. Nuclear Coordinates

The variation of spectroscopic properties upon

JOURNAL OF COMPUTATIONAL CHEMISTRY 951

in particular, on the central electron-pair bond that

FIGURE 9. electron-pair (left) and donoracceptor

952 VOL. 22, NO. 9

gen. The increase in stabilizing orbital interactions

JOURNAL OF COMPUTATIONAL CHEMISTRY 953

dialkylsulfide dimer radical cations (R2 SSR2 + ),

954 VOL. 22, NO. 9

and the less AH antibonding nature of the LUMO.

JOURNAL OF COMPUTATIONAL CHEMISTRY 955

ligands are farther removed from each other. Elec-

956 VOL. 22, NO. 9

unexpectedly, due to the charge separation that

JOURNAL OF COMPUTATIONAL CHEMISTRY 957

electrons each) while the LiLi bonding 1a1 orbital

958 VOL. 22, NO. 9

dependent of any model or decomposition of the

JOURNAL OF COMPUTATIONAL CHEMISTRY 959