вход по аккаунту


Molecular Structures from Density Functional Calculations with Simulated Annealing.

код для вставкиСкачать
Molecular Structures from Density Functional Calculations
with Simulated Annealing
By Robert 0. Jones*
The geometrical structure of any aggregate of atoms is one of its basic properties and, in
principle, straightforward to predict. One chooses a structure, determines the total energy E
of the system of electrons and ions, and repeats the calculation for all possible geometries. The
ground state structure is that with the lowest energy. A quantum mechanical calculation of the
exact wave function Y would lead to the total energy, but this is practicable only in very small
molecules. Furthermore, the number of local minima in the energy surface increases dramatically with increasing molecular size. While traditional ab initio methods have had many
impressive successes, these difficulties have meant that they have focused on systems with
relatively few local minima, or have used experiment or experience to limit the range of
geometries studied. On the other hand, calculations for much larger molecules and extended
systems are often forced to use simplifying assumptions about the interatomic forces that limit
their predictive capability. The approach described here avoids both of these extremes : Total
energies of predictive value are calculated without using semi-empirical force laws, and the
problem of multiple minima in the energy surface is addressed. The density functional formalism, with a local density approximation for the exchange-correlation energy, allows one to
calculate the total energy for a given geometry in an efficient, if approximate, manner. Calculations for heavier elements are not significantly more difficult than for those in the first row
and provide an ideal way to study bonding trends. When coupled with finite-temperature
molecular dynamics, this formalism can avoid many of the energetically unfavorable minima
in the energy surface. We show here that the method leads to surprising and exciting results.
1. Introduction
A knowledge of the exact wave function Y of an interacting system of electrons and ions would allow one to determine many quantities of interest, including the total energy
E (Eq. l), where 2 is the Hamiltonian of the system. Calcu-
lations of this type can be carried out for systems with very
few electrons, the most familiar example being the hydrogen
atom. However, approximations for Y are generally unavoidable, and are usually based on the variational principle
of RayLeigh and Ritz: If I@) is an approximate wave function, then Equation 2 is valid, where E,, is the exact energy
of the ground state.
The simplest form of a many-particle wave function is
that of Hurtree (1928), who represented @ as a product of
single-electron functions (Eq. 3), where each di satisfies a
Schrodinger equation whose potential term is given by the
mean field of all other electrons (Eq. 4). V,,, is the external
Dr R. 0.Jones
Institut fur Festkorperforschung
Forschungszentrum Julich, W-5170Julich (FRG)
Verlagsgeselisshafr mbH. W-6940 Weinheim. 1991
field of the nuclei, and the Coulomb potential cp is determined by solving Poisson’s equation (we use atomic units,
i.e., the unit of charge is the electron charge e, of distance the
Bohr radius uorand of mass the electron mass m), Replacing
the product by a single (Slater) determinant-the HartreeFock (HF) approximation-leads to an additional “exchange” potential, without changing the single-particle picture. Although this has been the basic method of atomic and
molecular physics for many years, a single determinant does
not generally lead to a satisfactory value of the total energy.
While a linear combination of determinants (“configuration
interaction”, CI) should lead in principle to an exact wave
function and energy, the numerical effort required increases
explosively with increasing electron number.
These arguments apply to a one geometry, that is, for a
particular set of internal coordinates of a system. If we have
N atoms in an aggregate ( N > 2), the total number of internal coordinates is 3N - 6, since the internal energy of the
aggregate is independent of any translational or rotational
motion. With so many independent interatomic degrees of
freedom, the scale of the problem is clear: 1) Finding the
energy minimum by calculating the exact energy for all possible geometries is impracticable, since it would require a
vast number of calculations in a system where one is difficult!
2 ) The ability to locate the closest minimum to a particular
geometry would simplify the calculations greatly, particularly if the number of possible isomers is small. In general,
however, it is difficult to predict the most stable structure,
since the number of local minima in the energy surface can
be very large indeed.
The enumeration of topologically different isomers consistent with a given chemical formula is one of the oldest problems in theoretical chemistry. In 1874, only a decade after the
0570-083319lj0606-0630$3.50+ 2510
Angew. Chem. Inf. Ed. Engl. 30 (1991) 630-640
chemical bond was first denoted by a line joining the nuclei,
Cayle-v[’*showed that the number of isomers could grow
rapidly with increasing N , and much work in the intervening
period has confirmed this. Hoare and McInnes,[zl for example, located all the minima in small aggregates with
simple, pairwise interatomic interactions, and found a rapid,
perhaps exponential, increase in their number with increasing N . In fact, for aggregates of identical atoms interacting
with a pairwise potential, Wille and Vennikt3]have shown
that there is no known algorithm for determining the ground
state energy and structure that grows with time as a power of
N . Such problems are classified as “NP-hard”[41 and are
termed intractable. It is indeed sobering when mathematicians derive such a result. Moreover, the parameters in force
laws such as those used by Hoare and Mclnnes are usually
derived from experiment, and one cannot be sure that the
results are not affected by the choice of parameters.
This discussion shows that there are two distinct problems
to be solved, and the difficulties are multiplicative in nature.
The approach we now describe addresses both aspects. First,
the density functional formalism provides a numerically eficient scheme for calculating the ground state energy of the
interacting system of ions and electrons without attempting
to calculate the wave function Y . Second, we use the strategy
of simulated annealing (SA) to determine low-lying energy
minima. We show that, even in small molecules, there are
energetically favorable structures that had not been found
on the basis of experiment or intuition. In Section 2 we outline the density functional (DF) formalism and discuss the
form of the exchange<orrelation energy Ex,and local density approximations to it. The results of local spin density
(LSD) calculations given in Section 3 illustrate the level of
accuracy to be expected, and we show in Section 4 how DF
calculations can be combined with molecular dynamics. We
present some applications in Section 5 and our concluding
remarks in Section 6.
2. Density Functional Formalism
2.1. Formal Basis
The basic theorems of the density functional formalism
’ ’ showed that:
were derived by Hohenberg and K ~ h n . ~They
1) Ground state (GS) properties of a system of electrons
and ions in an external field Vex,can be determined from a
knowledge of the electron density n(r)alone. The total energy
E is such a functional of the density, that is, E = E(n).
2) E(n) satisfies the variational principle E(n) 2 EGs,and
the density for which the equality holds is the ground state
density nGs.
A simple and general proof of these assertions has been given
by Levy.[61
The step to a practical scheme for energy calculations
resulted from the observation by Kohn and Sham[’] that the
minimization of E(n) can be simplified if we write E as in
Equation 5. Here To is the kinetic energy that a system with
density n would have in the absence of electron-electron
interactions, q ( r ) is the Coulomb potential, and the remainder Ex, is our definition of the exchangesorrelation energy.
The variational principle now yields Equation 6, where p is
the Lagrange multiplier associated with the requirement of
constant particle number.
Equation 6 applies to the interacting system of electrons
and ions. The solution of the corresponding equation for a
system of noninteracting particles, Equation 7, can be found
6E(n) 6To
6n(r) - &(r)
+ V ( r )= p
by solving the Schrodinger equation, Equation 8a. The density is given by Equation 8 b, where the sum runs over orbitals i with occupation numbersJ;. The problems (6) and (7)
are mathematically identical, provided that V(r) is defined
according to Equation 9. This condition can be satisfied in
practice by a “self-consistency’’ cycle: From a starting value
of n we compute V(r) and the n’ that results, and continue
this procedure until n = n’.
Robert 0 .Jones was born in 1941 in Kojonup, Western Australia. He graduated in 1963from the
University of Western Australia, Perth, and obtained his Ph.D. in 1967 from the University of
Cambridge, England, under K Heine. From 1967 to 1970 he worked as Research Associate,
Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, N X USA and since
1970 as Scientist at the Forschungszentrum Jiilich (KFA). In 1975 he was NORDITA Professor
at the Chalmers Technical University, Goteborg, Sweden and has been Fellow of the Australian
Institute of Physics since 1978.
Angen. Chem. lnr. Ed. Engl.
30 (1991) 630-640
The problem of determining the total energy of a system
of electrons and ions has therefore been reduced to the solution of a single-particle equation of Hartree form (cf. Eq. 4),
leading to the energy and electron density of the lowest state
and all quantities derivable from them. This method can be
extended to the lowest energy state of a given symmetry by
restricting the variation to densities compatible with the corresponding symmetry quantum numbers.['] In the form of
Equations 8 and 9, density functional theory is an orbital
theory and it shares the main advantage of the Hartree-Fock
method, namely the interpretation of the results in a singleparticle picture. In contrast to the H F potential, however,
the effective potential Y(r) has a local dependence on the
This situation may seem to be almost too good to be true.
Before celebrating, however, we should note the following:
1) The total energy E of the most stable state of a given
symmetry is found by solving equations describing a fictitious system of noninteracting "electrons", and the eigenvalues ci and eigenfunctions ll/i in Equation 8 can not be viewed
as excitation energies and wave functions of the interacting
system. The use of such fictitious systems leads to additional
complications, since it is not always possible to draw an
unambiguous connection between the states of the interacting and noninteracting systems. In the latter, the f ; will be
compatible with a single determinant-or at most a small
number of determinants-of the same symmetry.["
2) While in the energy expression ( 5 ) To, the electron-ion
and electron4ectron interaction energies can readily be
evaluated, Ex, is dtlfined only as the difference between the
sum of these terms and the exact energy. It is a small but far
from negligible part of the total energy of atoms and molecules. Moreover, if there is no simple expression for Ex,,how
do we determine in practice its functional derivative in Equation 6? Approximations to Ex, therefore play an essential
part in applications of D F theory, and we now study its form
in more detail.
2.2. Exchange-Correlation Energy
The connection between the noninteracting and interacting systems-both of density n ( r + c a n be drawn simply by
considering an electron-electron interaction 14 and increasing i.
slowly from 0 to 1. One can use this construction to
derive an exact expression (10) for Ex,,['] in which nxc is defined as in Equation l l . The pair-correlation function
g(r, J , A ) describes the fact that an electron at point r reduces
the probability of finding one at J, so that Ex,can be viewed
as the interaction between n(r) and the exchange-correlation
hole n,,(r). This picture is, of course, very similar to that of
the exchange hole introduced by Slarer for electrons of parallel spin."']
nxc(r,J - r)
= n(r')
dl(g(r, r', A)
- 1)
Ex, = 2 drn(r) $ dR RZR
- 1dSZn,,(r, R )
nxe(r,R), so that approximations for it can give an accurate
value even if the description of the nonspherical parts of n,,
is very poor. Furthermore, a sum rule following from the
definition of the pair-correlation function requires that the
exchange-correlation hole contain one electron, that is, for
all r , Equation 13 is valid. If we consider - nxc(r,r' - r) as
J dJn,,(r,
J - r) = - I
a normalized weight factor and define locally the radius of the
exchange-correlation hole (Eq. 14), Equation 15 follows.
Provided the sum rule (13) is satisfied, Ex,depends only weakly on the details of nxcand is determined by the first moment
of a function whose second moment is known exactly.
The experience of the last twenty years shows that the
exact exchange-correlation energy functional E,,(n(r)) must
be extremely complicated. The calculation of correlation
functions (see Eq. 11, for example) generally requires the
knowledge of the exact wave function. Nevertheless, the
above argument shows that simple, useful approximations to
Ex,may exist. The practical necessity of approximating Ex,is
an essential difference between the density functional and
configuration interaction approaches. While the latter seeks
an exact numerical solution of the Schrodinger equation and
precise values of related quantities, even the best solution of
the DF equations can only reflect the accuracy of the approximation used for Exc. Furthermore, while the CI approach can be improved systematically by increasing the
number of determinants and basis functions, a corresponding scheme for improving the accuracy of LD calculations is
presently not known.
2.3. Local Density Approximations for Ex,
Local density calculations have a long history in condensed matter physics, where a special role has been played
by the approximate exchange potential introduced by
Slater.["l It is sometimes overlooked that Slater's derivation
was based on the concept of the exchange-hole discussed
above, and not on the particular form of the exchange energy
density obtained for a homogeneous electron gas. The approximation for Ex,in most widespread use today is the local
spin density (LSD) approximation (16). Here c,,(n,, n i ) is
(1 1)
The isotropic nature of the Coulomb interaction suggests
the variable substitution R = J - r in (lo), leading to
Equation 12.['l Ex,depends only on the spherical average of
the exchange and correlation energy per particle of a homogeneous, spin-polarized electron gas with spin-up and spinAngew. Chem. Int. Ed. Engl. 30 (1991) 630-640
down densities n , and n , , respectively. This approximation
is free of adjustable parameters-&,, is taken from the best
available calculations for the uniform electron gas-and its
use can be justified in systems with small departures from
homogeneity. Density distributions are far from homogeneous in atoms, molecules, and most solids, and the usefulness of the approximation relies on arguments such as those
in Section 2.2. The approximation satisfies the sum rule (13);
more details are given in the review article by Jones and
Gunnarsson." '1
small for H, and Li,. In cases where it leads to an energy
minimum, the H F approximation usually leads to a small
underestimate of the equilibrium interatomic separations.
The values for r, calculated using the LSD approximations
are generally in good agreement with experiment. An overestimate of 1 - 2 % is common.['5. l61 Ground state vibration
frequencies are correspondingly low. The dipole moment of
CO and its variation with internuclear separation are given
significantly better by the LSD approximation than by H F
calculations.[' 7 J
3. DF Calculations with Local Density
3.2. Alkaline Earth Dimers Be2, Mg,, Ca,
Density functional calculations with the LSD approximation are common in extended systems, and the results show
some general features. Equilibrium geometries are given reliably in most materials, and energy variations in the neighborhood of local minima are described well, so that vibration
(phonon) frequencies are generally in reasonable agreement
with experiment. In this section, we select examples from the
applications to molecular systems[' 'I to indicate the level of
accuracy to be expected in such systems. They illustrate two
strengths of the method: the ability to study molecules containing elements from the whole periodic table, and the advantages of an interpretation based on a single-particle picture. They also demonstrate the tendency of LD calculations
to overestimate strengths of bonds involving sp orbitals.
3.1. First-row Diatomic Molecules
Diatomic molecules of first-row atoms Li-F have served
as a testing ground for numerous methods of electronic
structure calculation. In Table I we compare measured potential well depths['2] in first-row dimers (homonuclear di-
While the LSD and H F approximations generally lead to
similar bond lengths in first-row molecules, Be, is quite different. The lowest-lying state in the alkaline earth dimers
[I a,(fl> 1uu(lJ)]) has equal occupancies of bonding
and antibonding orbitals, and Hartree-Fock calculations['81
lead to repulsive curves. Bonds in these molecules have often
been referred to as Van der Waals bonds, in the expectation
that the bond strength increases monotonically Be, + Mg,
-+ Ca, . . . with increasing atomic polarizabilities. The earliest CI calculations supported this picture, since they predicted for Be, a very weak bond with an equilibrium separation
near 9 a, .['91
The local density approximation gives qualitatively different results and a different picture of bonding in this homologous series.[201The energy curves show minima in all
cases, and the equilibrium separations in Mg, and Ca, agree
well with measured values. The bond length in Be, (4.86 a,)
was found to be much shorter than in CI calculations. Most
striking, however, is the variation in bond energies shown in
Figure 1 . The bond in Be, is stronger than in Mg,, and the
Table 1. Experimental and calculated well depths [eV] for the experimental
ground states of the first-row dimers.
Expt [a]
LSD [h]
XO [b, cl
HF [dl
- 1.37
[a] For Be, from Ref. [25],otherwise Ref. (121. [b] Ref. [15]. [c] Local density
calculations with exchange only. [d] Total energies for experimental geometries
from Ref. [13]. [el HF calculations for Be, give a purely repulsive energy curve.
atomic molecules) with values calculated using HartreeFock,t'31 and local spin density approximations for Ex, and
Ex (exchange only).[7%
14,151 The H F approximation leads to
substantial underestimates of the binding energies, particularly for singlet ground states. The LSD values overestimate
the stability of these molecules, particularly beyond carbon
in this row, although the deviations from experiment are
Angew. Chrm. Int. Ed. Engl. 30 (1991) 630-640
He Be Mg Ca Sr Ba Ra
Fig. 1. Bond energies EB calculated for the
state of alkaline earth dimers
(full curve, left scale, Ref. [ZO]) and cohesive energies of bulk materials (broken
curve, right scale, ref. [22]). Experimental values Refs. [12] and [21]) are given
as crosses.
variation with atomic number parallels the irregular behavior observed in the bulk cohesive energies (Fig. 1, experiment: ref. [21]; LD calculations: ref. [22]). The bond energies in the molecules and solids are overestimated by the LD
approximation in those cases where experimental values are
known. The error is significant in Be,, as we show in Table 2,
where we include the results of more recent LD and CI calcul a t i o n ~ [241
~ ~and
* gas-phase experimental result^.[^^^ The D F
and experimental equilibrium separations are in very good
agreement, and the extensive CI calculations reproduce both
the well depth and equilibrium separation satisfactorily.
Table 2. Experimental and calculated spectroscopic constants of Be,
Experiment [25]
CI [24]
CI [23]
4.73 5 0.03
we [cm-'1
D, lev]
. ..
ca. 0.1 1
0.09 k 0.01
3.3. Diatomic Molecules of the Fourth Main Group:
C, ,Si,, Ge,
The absence of core p states, the strong attractive potential
experienced by the 2p states, and their relatively compact
nature apply to all first-row atoms. As a consequence, their
bonding properties differ qualitatively from those of atoms
in the remainder of the periodic table. The experimental
and the excitation
ground state in C, , for example, is
energy (T, = 0.09 eV) to the 3n,(Z0,n,3)
is very small.I'2]
The remarkable ease of the og-+ K, transfer is not found in
the heavier dimers of this group; these have '2; (2 0: 1 K:)
ground states, excitation energies of 1.0 to 1.5 eV to the
state (1 K,"), and substantially weaker n
[a] Estimated. [b] R. 0. Jones, unpublished results (1982).
To discuss the bonding in this homologous series, it is
useful to consider the valence orbitals (Fig. 2) of the atoms.
With the exception of He,, there is a substantial overlap
between the electron densities on the two atoms, particularly
in Be,. Furthermore, the radial extent of the orbitals does
not increase smoothly with core size, but shows a "secondary
periodicity". For example, the valence s orbital in Mg is not
Fig. 3. Radial valence orbital functions for atoms of the fourth main group: a)
s functions, b) p functions.
Fig. 2. Radial valence orbital functions for helium and alkaline earth atoms: a)
s functions for the ' S (ns') states, b) p functions for the 3P (ns'np') state. The
dashed curve is the 2p function of He (ls2p). The arrows denote half of the
calculated equilibrium separation (rJ2) in each dimer [20].
The radial functions of the valence orbitals of these elements (Fig. 3) show a similar behavior to that discussed
above. The valence functions of Si and Ge are remarkably
much more extended than that in Be. The 2p functions in the
Mg core have no repulsive (orthogonalization) effect on the
3s function and are sufficiently extended that the additional
core charge is incompletely screened. A similar orbital compression is evident in Sr and Ra, where 3d and 4f functions
enter the core for the first time. The unoccupied p orbitals
are generally more extended than the s valence functions, but
Figure 2 shows that Be is an exception. The similar extent of
the 2s and 2p orbitals suggests that sp mixing will be favored.
This picture is in accord with the pseudopotential theory of
cohesion in the bulk: a stronger sp hybridization (and larger
band gap) in Be than in Mg, and increasing sd hybridization
with increasing atomic number. It is encouraging that the
trends predicted by the DF calculations, particularly the relative strength of the bond in Be,, were confirmed by both
exhaustive CI calculations and by experiment.
-1.00 -
Fig. 4. Dependence of the self-consistent eigenvalues for 'Egstate ofC,, Si,,
and Ge, on the interatomic distance r. The valence eigenvalues for the isolated
atoms are given and the arrows on the r axis denote half of the calculated
equilibrium separation ( r J 2 ) in each case.
Angew. Chern. Inf. Ed. Engl. 30 (f991) 630-640
similar, since the relatively diffuse 3d core density in Ge
screens the additional nuclear charge imperfectly. This is
reflected in very similar bulk properties, such as ground state
geometries, and in the eigenvalues of Si, and Ge, shown in
Figure 4. The C, eigenvalues reflect the compact valence p
function and show why the multiplet structure is so different.
state obtained by placing four electrons in the nu
orbital has a very similar energy to the "C, state obtained by
transferring two electrons to the 2 crgorbital. The eigenvalues
in Figure 4 suggest that occupancy of the latter orbital will
be favored in the heavier dimers, resulting in weaker n
bonds. The consequences for the formation of simple saturated molecules have been discussed by Harris and Jones.[26J
3.4. Triatomic Molecules of the Sixth Main Group:
As a further example, we study the energy hypersurfaces
of 0, and SO,,[27Jwhich have ground states with 'A,(C,J
symmetry and similar bond angles (116.8" and 119.4", respectively). Although the absorption spectrum of ozone has
been studied in detail, particularly for ultraviolet radiation,
little experimental information is available on the excited
state energy surfaces, and theoretical studies have been essential. Attention has been focused in the past on the presence of two minima in the 'A, energy hypersurfaces and the
energy difference between them. The H F approximation
gives a qualitatively incorrect ordering of the low-lying states
of 0,, since there are two configurations with substantial
contributions to the ground state wave function.
In Figure 5 we show the energy hypersurfaces for lowlying states of ozone. D F calculations reproduce the ground
However, there are qualitative differences between the energy ordering of the states. In particular, the excitation from
the ground state to the low-lying triplet states requires substantially more energy for SO, than for O,, and the energy
separation between the two 'A, states (4.1 eV) is almost
three times as large for SO,.
These differences can be related to the larger radius of the
sulfur atom and the higher energy of the p 0rbita1.l'~~The
ground state (x' A,) has the valence configuration
1a: 3 b: 4a: 2 by, and the 2 'A, state corresponds to the excitation (3 b: -+ 2 b;). While the 1 a, and 3 b, orbitals have
significant contributions only from the outer (oxygen)
atoms, the higher-lying 4a, and 2b, orbitals have a strong
contribution from the central atom (0in O,, S in SO,). The
eigenvalue spread in SO, is therefore much larger, so that the
large excitation energy xi A, + 2 'A, in SO, reflects simply
the double excitation to the high-lying 2b, orbital.
The above molecules are only a small sample of those that
have been studied with the D F method. Nevertheless, they
illustrate some general features of the results. As in extended
systems, equilibrium geometries and vibration frequencies
are described satisfactorily, even in molecules-such as Be,
and 0,-where the Hartree-Fock method gives qualitatively
incorrect results. While energy variations near equilibrium
geometries are given satisfactorily, this is not true of all regions of configuration space. The errors in the dissociation
energies of first-row diatomic molecules and the ozone family (up to ca. 2 eV), for example, are unacceptably large and
have contributed to the motivation behind the search for
improved approximations.["] The single-particle picture
and the relative ease of performing calculations for elements
with large atomic numbers is very useful in understanding
the differences between first-row and heavier elements. A
discussion on similar lines was presented subsequently by
4. The Molecular Dynamics/Density Functional
EB lev1
Fig. 5. Energy surfaces for low-lying states of 0,. ec = bond angle. The values
on the curves give the bond lengths [a,] that optimize the energy for each
angle a.
state geometry well, and the excitation energy between the
two 'A, minima (1.4 eV) is in reasonable agreement with the
most recent CI calculations, which give results between
1.O eV and 1.4 eV.1271In both 0, and SO,, the bond lengths
increase in the order x' A, -+ 1 3B, -+ 1 3B, -+ 2 'A,, and the
bond angles in the order 2'A, + 1 3B, -+x'A, -+ I3B,.
Angen,. Chem. I n [ . Ed. Engl. 30 (1991) 630-640
While we cannot expect to be able to determine the exact
ground state if the energy hypersurface has many local minima, it is essential to develop methods for avoiding energetically unfavorable local minima. Kirkpatrick et al.1291noted
that this can be achieved by allowing the system to evolve at
finite temperature and implemented this "simulated annealing" strategy using a Monte Carlo sampling. Molecular dynamics (MD) provides an alternative, and Car and Parrinelshowed that it could be combined with the density
functional (DF) scheme for calculating total energies. The
use of the LSD approximation for the exchangexorrelation
energy avoids the parameterization of the interatomic forces
common in MD schemes,[311and should give a reliable description of energy variations for large regions of configuration space.
To perform geometry optimization using D F calculations,
we face two minimization problems: The total energy must
be minimized for each geometry by varying the density (the
solution of Eq. 7), and the geometry with the lowest energy
must be found. This is done conveniently by viewing E as a
function of two interdependent sets of degrees of freedom,
{I)~} and { R , } , and using standard MD techniques (Eq. 17,
where 2, are the ionic charges).[301
The energy minimization required for structure determination can be performed using standard MD techniques,[301
since the evolution of the system of electrons and ions is
reproduced exactly by the evolution of a dynamic system
defined by the Lagrangian (18) and the equation of motion 19a and 19b. M I and R, denote the ionic masses and
coordinates, pi are fictitious "masses" associated with the
electronic degrees of freedom, dots denote time derivatives,
and the Lagrangian multipliers Aij are introduced to satisfy
the orthonormality constraints on the tji(r, t ) . From these
orbitals and the resultant density n(r, t ) = C I$i(r, t ) l z we
5.1. Sulfur Clusters S,
The elements of the sixth main group provide some of the
best characterized of small atomic clusters. These elements
are unique in that many allotropes comprise regular arrays
of well-separated ring molecules, and X-ray structure analyses have been performed for S,, n = 6- 8,lO - 13,18,20, and
Sen, n = 6, 8. S, has been found as microcrystals, and S,,
exists in two distinct forms. Several mixed crystals of the
form Se,S, and a range of sulfur oxides and ions are also
known. The preparation and structure of these clusters has
been reviewed by S t e ~ d e l . ~ ~ ~ ]
The theoretical results for sulfur clusters S , to S, 3[341 show
that it is possible to determine low-lying energy minima even
if the starting geometries are far from correct: linear chains
in S, to S, and nearly planar rings in S, to S 1 3 .The final
structures agreed well with experiment in all cases where
X-ray data were available. For example, the D,,symmetry in
S,, (Fig. 6b) is reproduced very closely, and the structural
parameters (bond length d = 3.97a0, bond angle ci = 106",
dihedral angle y = 88') are in good agreement with measured values (3.88a0, OL = 106.2", y = 87.2'). The discrepancy in bond length is reduced by using a more accurate pseudopotential, as shown by calculations for the D,, form of S,,
where the structure ( d = 3.900a0, LY = 108.5", y = 98.3") is
very close to the experimental result ( d = 3.893 0.006a0,
ci = 108.0 f 0.7", y = 98.3 f 2.1").[3s]
use Equation 17 to evaluate the total energy E, which acts as
the classical potential energy in the Lagrangian function 18.
The artificial Newton's dynamics for the electronic degrees
of freedom, together with the assumption < M I , prevent
transfer of energy from the classical to the quantum degrees
of freedom, even for long simulations.
5. Applications of the MD/DF Scheme to Molecules
MD/DF calculations have been performed for several
families of molecules containing main-group elements, using
a large unit cell with periodic boundary conditions to guarantee a weak interaction between the molecules. We give
some examples here to show the value of the approach in
several contexts.[321The scattering properties of the individual atoms are described by nonlocal pseudopotentials, and
the eigenfunctions +i are expanded in a plane wave basis. We
set the velocities di and R equal to zero for an arbitrary
starting geometry, and use an efficient self-consistent iterative diagonalization technique to determine those iithat
minimize E. With the electrons initially in their ground state,
the dynamics (Eq. 19a and b) generate Born-Oppenheimer
(BO) trajectories without the need for additional diagonalization/self-consistency cycles for the electrons. In the present applications, the minimum of the potential energy surfaces can be found by alternating MD calculations at finite
temperature with steepest descents (or conjugate gradient)
determinations of the nearest local minimum. Further details
are given in the original papers.
Fig. 6 . Structures of a) S , and b) S , ,
Cases where the most stable isomers are not yet known are
of particular interest. There is no experimental information
on S, and S, ,for example, and S, has only been prepared as
crystallites that are too small for X-ray diffraction analysis.
The most stable S, isomer in our calculations was an unusual
rectangular structure made up of two weakly interacting
dimers. This minimum of the energy surfaces is very shallow
and it is easy to distort the molecule. The detailed ab initio
calculations of Rughavachuri et aI.[361support our assignment of the ground state. This is also true in s,, where we
found a ground state isomer with C, symmetry (Fig. 6a) as
well as several local minima with C, symmetry with energies
> 0.2 eV higher. The calculated ground state structure is
consistent with measured Raman spectra, which indicate
narrow ranges of bond lengths and bond angles, and dihedral angles in the range 70- 130".
5.2. Structural Change in S,O
A surprising example of structural optimization is provided
by the S,O molecule.[37]In the course of the simulation
Angew. Chem. Inr. Ed. Engl. 30 (1991) 630-64Q
shown in Figure 7-1200 time steps at T = 2000K-the
molecule changed from one with a stable local minimum and
a ring structure similar to that In S, (Fig. 7a) to a structure
with different topology: an oxygen atom outside the S,
structure (Fig. 79. The structure found on reducing the temperature SlOWlY to T = 0 K (Fig. 7J) agrees well with the
experimental ground ~tate.1~’’
The entire structural change
I1 21
il 31
I1 I 1
I1 51
Flg 8 Possible cyclic structures for Se,S, and Se,S, The numbers in parentheses indicate the positions of the two atoms of the minority element
Calculations have been carried out on seven- and eightmembered rings of the type Se,S,,’471 starting with the
ground state structures of S, (C,) or S, (D4J. The eightmembered rings are characterized by small deviations from
the crown-shaped (D4J structures, and the possible cyclic
structures in SezS6 and Se,S, are shown in Figure 8. The
results for these molecules show several remarkable features,
illustrated by the structural parameters for Se,S, given in
Table 3. In all eight structures, the bond lengths (ds-s,d j _ s , ,
Table 3. Molecular parameters d[u.]. a and y [ “ ] for all isomers of Se2S6(see
Fig. 8). Bond lengths between atoms i and j are denoted by d,,,bond angles at
atom i by a, and dihedral angles at bond ij by yi,.
Fig. 7. Evolution of S,O from ring (a) to ground state structure 6). The time
difference between successive frames is 1.1 x
(including cooling) took only about 8000 time steps
(< 10-’2s), although the energy barrier between ring and
ground state structures is ca. 5 eV and the energy difference
only about 0.2 eV. This shows just how efficient molecular
dynamics can be in generating geometrical configurations.
5.3. Isomers of Se,S,, Se,S,
Sulfur and selenium have much in common, and binary systems have been studied extensively in the vapor phaser3*,391
and as liquids and as solid solution^.^^^-^^] The different
species tend to crystallize together, leading to structures with
sulfur and selenium atoms distributed over the atomic
The possible complexity is evident in the eight-membered rings S,Se, -”, where thirty different crown-shaped
isomers can occur. Apart from the intrinsic interest in molecules of the form Se,S,, an understanding of their structures
could provide insight into the structures of related liquid and
amorphous materials, including Se itself.
The interconversion of various Se,S, isomers and ringchain equilibration processes involves reactions of the type
in (20). This reaction is endothermic in the
+ -Se-Se-
2 -Se-S-
diatomic S-Se molecules in the gas phase,[451the measured
enthalpy change is +7.9 f 6.7kJmol-’ and the sum of the
dissociation energies on the left side of (20) is greater than
the sum on the right. Theoretical predictions of the relative
stability of different isomers are difficult to make,’461 and it
is natural to ask whether DF calculations can describe reliably the small energy differences involved.
Angew. Chem. Inr. Ed. Engl. 30 (1991) 630-640
ds,-s,) are the same to within 0.02 a, (3.93 aor4.21 a,, 4.44a0),
and ds-s and dS,+, are the same as found previously for the
ring structures of S, and Se,. The (1,2) isomers are the most
stable in both Se,S, and Se,S,, and the energies of the other
three structures are ca. 0.08 eV higher and degenerate to
within 3 meV in each case. We now show that these findings can be understood in terms of a very simple model.
The transferability of the bond lengths and the known
connection between bond length and bond strengthrz6]suggest that we can associate with each bond type a contribution
to the binding energy, E,-,, Es-se, and Es,-s,. This leads to
(21) for the (1,2) structure in Se,S6 and to (22) for all remain-
ing structures. This argument predicts that the energies of
the last three structures will be equal, and will differ from the
energy of the (1,2) structure by A E (23).
If we apply the same argument to the Se,S, structures,
Equations 21 -23 still apply when S and Se are interchanged.
This explains why the energy orderings in Se,S, and Se,S,
are the same, and why the most stable (1,2) isomers are
separated from three structures with almost identical energies. The separation in the present calculations (ca. 0.08 eV,
7.4 kJmol- ’) is consistent with the endothermic reaction
(20) and in reasonable agreement with the measured heats of
reaction in the gas and liquid phases.
In spite of its simplicity, this model-based on a very
careful geometry optimization-gives results that agree with
trends found in recent measurements, particularly the relative abundance of structures in which the minority atoms are
adjacent. The symmetry between the energy ordering of the
isomers of Se-rich and S-rich molecules is perturbed, of
course, by differences between the elements such as the
atomic radius, but the essential features should be unchanged. The trend to segregation of minority components
could also occur in the disordered liquid and amorphous
states. Segregation in this model is a consequence of the sign
of B E in Equation 23, and other systems with a small value
of A E or with one of opposite sign would behave differently.
5.4. Phosphorus Clusters P,
The structural variety shown by phosphorus exceeds that
of all elements other than sulfur and perhaps boron, and
there have been many studies of its crystalline and amorphous
A problem of continuing interest in gasphase clusters has been the detection of only p,, p,, p,, and
(possibly) P, in the vapor phase at high temperatures.[501
There have been numerous speculations about the apparent
absence of heavier clusters,[5 and several calculations indicate that the cubic form of P, is less stable than two P,
tetrahedra.[52.5 3 1 Recently, Martint541was able to detect by
mass spectroscopy P, molecules up to n = 24 in a low-temperature beam. In spite of this, almost nothing is known
about the structure of clusters with n > 4.
MD/DF calculations have been performed on phosphorus
clusters up to PlO,[sJ1
and calculated and experimental geometries and vibration frequencies agree well in those cases (P,,
P4) where the latter are known. In the larger clusters there
were several unexpected findings: 1) Although the tetrahedral structure is energetically favored in P,, there is a large
“basin of attraction” for a D,, “roof” structure, that is, this
structure is the closest minimum for a large region of configuration space. 2) The roof structure is a prominent feature in
the low-lying isomers in P, to P,. The most stable isomers
found in P,, P,, and P, have a P,-roof with an additional
one, two, and three atoms, respectively.
The two forms of P, shown in Figure 9 are energetically
more stable than the cubic structure. The structure in Figure 9 a shows again the preference for the “roof’ structure
(two such units are bonded together), and the isomer shown
in Figure 9 b was the most stable we found. It can most easily
be understood as a (distorted) cube with one bond rotated
through 90” and is ca. 40 kcalmol-’ more stable than the
cubic form. This “wedge” or “cradle” structure is found as
a structural unit in violet (monoclinic, Hittorf) phosphoru~.[~~I
Fig. 9. Two isomers of P,.
There is a pronounced analogy between the structures of
the P, isomers and those of the corresponding isoelectronic
hydrocarbons C,H,. The cubic form of the latter (cubane)
has been prepared by Eaton and
and can be converted catalytically to the wedge-shaped form ~uneane.[~’]
third possible isomer of C,H, (analogous to Fig. 9a) was
estimated to have a strain energy intermediate between the
other two, in agreement with the energy ordering we found
in the P, isomers.15s1
The results for the P, clusters show two particularly interesting features: 1) The presence of a large basin of attraction
for the energetically unfavorable roof structure in P,. This
shows that the most favorable structures are not necessarily
surrounded by the most extensive minima in the energy surfaces. 2) The existence of two isomers of P, that are more
stable than the cubic form. Although attention had been
focused for many years exclusively on the last, our finitetemperature simulation led almost immediately to more
stable forms. The lowest energy for P, lies ca. 0.5 eV below
the energy of two P, tetrahedra in the LD calculations. It
would be interesting to see the energy difference found using
ab initio methods.
5.5. Aluminum Clusters Al,
The calculations for the sulfur and phosphorus clusters
reproduce known structures satisfactorily. This is also true
for small aluminum clusters, where spectroscopic data are
available for Al,, and ab initio calculations have been performed for a limited range of geometries up to Al, . [ 5 9 1 If we
can reproduce known structures of small clusters with our
methods, we may have confidence in our predictions for
larger clusters of these elements. Aluminum is a prototype
“simple metal” of the condensed matter theorist, which is
characterized by a weak electron-ion interaction and a valence electron charge distribution that is nearly uniform. It is
natural then that calculations of the electronic structure of
A1 clusters have often used a “spherical jellium” model,[601
where both the electron density and the positive-ion distribution are assumed to be uniform inside a sphere of appropriate size. The focus in such calculations has been on particularly stable electronic configurations with “magic numbers”
of electrons.
Angew. Chem. Int. Ed. Engl. 30 (1991) 630-640
MD/DF calculations up to n = 10 give structures for
n = 2-4 in good agreement with other calculations, and a
satisfactory description of variations in ionization energies.
There were, however, some unexpected results for larger
clusters,I6’] including the existence of a class of energetically
favorable structures that can by no means be viewed as
spherical. Figure 10 shows three such structures for Al,.
They can be viewed as distorted planar arrays of approximately equilateral triangles, connected so that adjacent triangles are either nearly coplanar (dihedral angle near zero)
or have dihedral angles ca. 35“-55”. These features are also
apparent in bulk aluminum (the face-centered-cubic structure comprises equilateral triangles with dihedral angles zero
or 54.7”) and in the t( form of the neighboring boron-group
element gallium, as shown in Figure 10d. There are, of
course, many possible structures of this type in molecules of
finite extent. It is surprising, nevertheless, that D F calculations show that they are of comparable stability to the more
compact “metallic” structures of Al,. The structure shown
in Figure 1Oc was, in fact, the most stable found in this
Fig. 10. a) -c) Three isomers of Al,. d) Crystal structure of a-Ga.
We have noted above how useful the D F approach has
been in finding trends in families of molecules, and in relating them to single-particle properties such as the valence
orbital. This is also true in the elements of the third main
group. We find that the bonds in gallium clusters are generally shorter than the corresponding bonds in aluminum clusters, and that this can again be correlated with the range of
the valence orbitals.[611
6. Concluding remarks
Theory should be able to contribute much to the determination of molecular structures. The calculation of the energy
surfaces and the location of low-lying minima should lead to
the most stable isomers, but there are severe problems: The
determination of the energy E from the exact wave function
is impracticable in all but the smallest molecules, and there
are usually many local minima in the energy surface.
The approach described here addresses both aspects. The
density functional formalism, with the local spin density approximation for the exchange-correlation energy, enables us
to carry out energy calculations in an efficient manner, and
the use of finite temperature molecular dynamics prevents
the molecule becoming trapped in unfavorable energy minima. The results have been encouraging. The structures of
small clusters of sulfur are described well and, in the case of
Angew. Clirm. Int. Ed. Engl. 30 (1991) 630-640
S,O, we can simulate structural changes involving a high
energy barrier. The final structure (an 0 atom outside an
S,-ring) is in good agreement with the structure determined
by X-ray diffraction. The heterocyclic ring molecules Se,S,
are characterized by structures with very small energy differences and energy barriers that are small on a thermal scale.
The calculations led in this case to a simple model of the
energy differences in mixed Se-S ring structures and predictions for molecules for which calculations have not yet been
performed. In the case of phosphorus clusters, we predict
structures with plausible trends and with an interesting analogy to isoelectronic hydrocarbons (CH), . In Al, clusters we
have found a class of energetically favorable structures that
was quite unexpected. There have been so many surprises, in
fact, that one can only assume that many calculations using
other methods have overlooked important minima in molecular energy surfaces. The number of minima and the fact that
many can have comparable stability show how difficult it is
to determine structures on the basis of “chemical intuition”
alone. There is little doubt that other methods of calculating
electronic structures would benefit from the implementation
of simulated annealing.
The successes of this approach should not lead us to overlook its limitations. Even the most careful calculation reflects the approximations used, and the LSD approximation
is an essential part of most of the calculations performed so
far. Overestimates in the binding energy are evident in all the
molecules discussed above, and there are indications that an
inadequate description of exchange energy differences is an
important factor.[”] Improved approximations would be
useful, if implemented in a systematic way, but attempts in
this direction have had mixed success. Moreover, there are
many local minima in the energy surfaces even for clusters
containing less than 10 atoms. The implementation of molecular dynamics techniques does not change this, and no method can guarantee finding the global energy minimum, or
even all of the most important minima. The development of
alternative optimization schemes and of simplified energy
functionals[621will be important in the future.
Density functional calculations are common in condensed
matter physics and it is not surprising that workers from this
field have played the major role in theoretical developments
and in molecular applications. In fact, density functional
methods have been slow to find favor amongst chemists. The
coupling to molecular dynamics has opened up perspectives
not currently available to some of the traditional methods of
theoretical chemistry, and-perhaps in conjunction with
these methods-the prospects of applications to many problems of chemical interest are very good.
I have benefited greatly from collaboration over the years
with 0. Gunnarsson, J. Harris, and D. Hohl, and from discussions with many colleagues. Special thanks are due to the
developers ojthe MDIDF method, R. Car and M . Purrinello.
and to G . Eilenberger for critical remarks on the manuscript.
Received: October 5. 1990 [A8381E]
German version: Angew. Chum. 103 (1991) 647
[l] A. Cayley, Philos. Mag. (4) 47 (1874) 444.
121 M . R. Hoare, J. A. McInnes, Adv. Phys. 32 (1983) 791.
[3] L. T. Wille, J. Vennik, J. Phys. A 18 (1985) L419, L1113.
[4] M. R. Garey, D. S. Johnson: Computers andlntractability: A Guide to the
Theory of NP-Completeness, Freeman, San Francisco 1979; NP-complete
means nondeterministic polynomial time complete [29].
[S] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864.
[6] M. Levy, Proc. Narl. Acad. Sci. USA 76 (1979) 6062.
[7] W. Kohn, L. J. Sham, Phys. Rev. 140 (1965) A1133.
[8] 0. Gunnarsson, B. 1. Lundqvist, Phys. Rev. ( E j 13 (1976) 4274; see also
J. Harris, R. 0. Jones, J. Phys. F4 (1974) 1170.
[9] J. Harris, P h p Rev. A 29 (1984) 1648.
[lo] J. C. Slater, Phys. Rev. 81 (1951) 385; ibid. 82 (1951) 538.
[ l l ] For a review on the formalism and some applications see R. 0. Jones, 0.
Gunnarsson, Rev. Mod. Phys. 61 (1989) 689.
[12] K. P. Huber, G. Herzberg: Molecular Spectra and Molecular Structure. I K
Consfants of Diatomic Molecules, Van Nostrand Reinhold, New York
[I31 P. E. Cade, A. C. Wahl, At. Data Nucl. Data Tables 13 (1974) 339.
[14] R. Gaspar, Acta Phys. Hung. 3 (1954) 263.
[IS] G. S . Painter, F. W. Averill, Phys. Rev. E 26 (1982) 1781.
[16] A. D. Becke, J. Chem. Phys. 84 (1986) 4524.
[I71 See, for example, 0. Gunnarsson, J. Harris, R. 0.Jones, J. Chem. Phys. 67
(1977) 3970.
[18] M. R. A. Blomberg, P. E. M. Siegbahn, B. 0.Roos, In!. J. Quantum Chem.
Quantum Chem. Symp. 14 (1980) 229.
[I91 M. R. A. Blomberg, P. E. M. Siegbahn, Int. J. Quantum Chem. 14 (1978)
[20] R. 0. Jones, J. Chem. Phys. 71 (1979) 1300.
[21] K. A. Gschneidner, Jr., Solid State Phys. 16 (1964) 276.
I221 V. L. Moruzzi, J. F. Janak, A. R. Williams: Calculaied Electronic Properties of Metals, Pergamon, New York 1978.
1231 R. J. Harrison, N. C. Handy, Chem. Phys. Lett. 98 (1983) 97.
[24] 9. H. Lengsfield 111, A. D. McLean, M. Yoshimine, 9. Liu, J. Chem. Phys.
79 (1983) 1891.
[25] V. E. Bondybey, J. H. English, J. Chem. Phys. 80 (1984) 568.
1261 J. Harris, R. 0. Jones, Phys. Rev. A 18 (1978) 2159; ibid. 19 (1979) 1813.
1271 R. 0. Jones, J. Chem. Phys. 82 (1985) 325.
[28] W. Kutzelnigg, Angew. Chem. 96 (1984) 262; Angew. Chem. Int. Ed. Engl.
23 (1984) 272.
I291 S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi. Science 220 (1983) 671.
[30] R. Car, M. Parrinello. Phys. Rev. Lett. 55 (1985) 2471.
[31] F. Stillinger, T. A. Weber, R. A. LaViolette, J. Chem. Phys. 85 (1986) 6460.
1321 See also R. Car, M. Parrinello, Phys. Rev. Lett. 60 (1988) 204; P. Ballone,
W. Andreoni, R. Car, M. Parrinello, ibid. 60 (1988) 271.
[33] R. Steudel, Stud. Inorg. Chem. 5 (1984) 3.
[34] D. Hohl, R. 0.Jones, R. Car, M. Parrinello, J. Chem. Phys. X9 (1988) 6823.
(351 R. 0. Jones, D. Hohl, Int. J. Quantum Chem. Quantum Chem. Symp. 24
(1990) 141.
[36] K . Raghavachari, C. M. Rohlfing, J. S. Binkley, J. Chem. Phys. 93 (1990)
[37] D. Hohl, R. 0. Jones, R. Car, M . Parrinello,J. Am. Chem. SOC.111 (1989)
[38] R. Cooper, J. V. Culka, J. Inorg. Nucl. Chem. 29 (1967) 1217.
[39] M. Schmidt, E. Wilhelm, Z . Nafurforsch. B 2 5 (1970) 1348.
[40] R. Steudel, E. M. Straws: The Chemistry oflnorganic Homo- and Heterocycles, Vol. 2, Academic, London 1987, p. 769.
[41] R. Steudel, R. Laitinen, Top. Curr. Chem. 102 (1982) 177.
[42] H. Bitterer (Ed.): Selenium, Gmelin Handbuch der Anorganischen Chemie,
Erganzungsband E2, Springer, Berlin 1984.
[43] R. Laitinen, L. Niinisto, R. Steudel, Acta Chem. Scand. Ser. A 33 (1979)
[44] T. Maekawa, T. Yokokawa, K. Niwa, Bull. Chem. SOC.Jpn. 46 (1973) 761.
[45] J. Drowart, S . Smoes, J. Chem. Soe. Faraday Trans. 2 73 (1977) 1755.
[46] R. Laitinen, T. Pakkanen, 1 Mol. Sfruct. THEUCHEM. 91 (1983) 337;
ibid. 124 (1985) 293.
[47] R. 0. Jones, D. Hohl, J. Am. Chem. Soe. I f 2 (1990) 2590.
1481 J. Donohue: The Structures of the Elements, Wiley, New York 1974,
Chap. 5 (boron group), 8 (nitrogen group), 9 (oxygen group).
(491 D. E. C. Corbridge: Phosphorus. An Outline of its Chemistry, Biochemistry
and Technology, Elsevier, Amsterdam 1985.
[SO] J. D. Carette, L. Kenuin, Can. J. Phys. 39 (1961) 1300.
[SI] H. Bock, H. Miiller, Inorg. Chem. 23 (1984) 4365.
1521 K. Raghavachari. R. C. Haddon, J. S. Binkley, Chem. Phys. Lett. 122
(1985) 219.
[S3] R. Ahlrichs, S. Brode, C. Ehrhardt, J. Am. Chem. SOC.107 (1985) 7260.
[54] T. P. Martin, 2. Phys. D 3 (1986) 221.
[55] P2 to P,: R. 0. Jones, D. Hohl, J. Chem. Phys. 92 (1990) 6710.
[56] H. Thurn, H. Krebs, Acta Crystallogr. Sect. E 25 (1969) 125.
[57] P. E. Eaton, T. W. Cole, Jr., J. Am. Chem. SOC.86 (1964) 962, 3157.
[58] L. Cassar, P. E. Eaton, J. Halpern, J. Am. Chem. SOC.92 (1970) 6366.
1591 T. H. Upton, J. Chem. Phys. 86 (1987) 7054; L. G. M. Petterson, C. W.
Bauschlicher, Jr., T. Halicioglu, J. Chem. Phys 87 (1987) 2205, and references cited therein.
[60] For a review, see W. A. deHeer, W. D. Knight, M. Y. Chou, M. L. Cohen,
Solid State Phys. 40 (1987) 94.
[61] R. 0. Jones, unpublished results.
1621 See, for example, J. Harris, D. Hohl, J Phys. Condens. Matter 2 (1990)
5161, and references cited therein.
Angew. Chem. Int. Ed. Engl. 30 (1991) 630-640
Без категории
Размер файла
1 160 Кб
structure, molecular, calculations, simulated, annealing, function, density
Пожаловаться на содержимое документа