Ab initio Calculation of Potential Energy Surfaces

1 The Concept of a Potential Energy Surface

This chapter is concerned with the calculation of potential energy surfaces by non-empirical methods, i.e. by obtaining solutions to Schrödinger's equation within the Born-Oppenheimer approximation. The concept of a potential energy surface is a consequence of the separation of the nuclear and electronic motions as proposed by Born and Oppenheimer. The nuclei may be considered to move under the influence of a potential determined by their mutual electrostatic repulsion and by the total energy of the electrons, an energy which is determined for every possible static configuration of the nuclei. The gain in conceptual simplicity afforded by the Born-Oppenheimer procedure is enormous and its use underlies many of our chemical concepts. For example, energies of activation or energy barriers in general, potential constants, the frequencies of vibrational and rotational motions, and bond lengths and bond angles as determined by an energy minimum, are all concepts defined in terms of the properties of a potential surface.

To obtain information concerning a system of n electrons and N nuclei in the absence of external fields, one must solve the time-independent Schrödinger equation

[FORMULA NOT REPRODUCIBLE IN ASCII] (1)

where the properties of the ith stationary state of the system of energy Ei are obtainable from the eigenfunctions Ψi(x,R) which in turn depend upon the space and spin co-ordinates of all the electrons and nuclei (whose collective co-ordinates are denoted by x and R, respectively) in the system. The total Hamiltonian operator is

[FORMULA NOT REPRODUCIBLE IN ASCII] (2)

where [??]e and [??N are the kinetic energy operators of the electrons and nuclei, respectively, and where the potential energy operator in atomic units is,

[FORMULA NOT REPRODUCIBLE IN ASCII] (3)

One can distinguish three levels of approximation in obtaining solutions to equation (1):2 (i) the Born-Oppenheimer approximation, often referred to as the clamped-nucleus approximation; (ii) the adiabatic approximation; and (iii) the non-adiabatic approximation. The solutions to equation (1) provided by the first of these are entirely adequate for most problems of chemical interest. For example, use of this approximation in the calculation of the low-lying vibrational and rotational levels of H+2 yields results which are correct within experimental accuracy.

In the Born-Oppenheimer approximation, one first transforms to a centre-of-mass, molecule-fixed co-ordinate system, thereby yielding a kinetic energy operator for the electrons referred to nuclei that are fixed with respect to their centre of mass and a nuclear kinetic energy operator which refers only to their internal motions, i.e. their vibrations and rotations and the various couplings between them.* Because of the large disparity in the masses of the electrons and the nuclei, the average kinetic energy of the former will be many times that of the latter, the ratio being proportional to mk/me ~ 2000 in the least favourable case when mk refers to the mass of a proton. Thus, in the limiting classical case, one obtains a picture of the electrons undergoing very rapid motion relative to the nuclei and adjusting almost instantaneously to changes in the nuclear positions. This suggests that the electronic motions are determined to a good approximation by just the static field of the nuclei, i.e. the electronic motion is determined by where the nuclei are but not by how fast they are moving. The motions of the nuclei in this separation of co-ordinates are governed by a potential whose negative gradient is simply the electrostatic force of repulsion between the nuclei and the attractive force exerted on the nuclei by the electronic charge distribution, a distribution whose form is determined by the electronic wavefunction evaluated for each configuration of the nuclei.

Mathematically, this separation of the kinetic motions of the electrons and nuclei amounts to approximating the total wavefunction Ψi (x, R) by a simple product of electronic and nuclear wavefunctions

[FORMULA NOT REPRODUCIBLE IN ASCII] (4)

The electronic function y>i(x; R) is obtained by solving the electronic Schrödinger equation for a fixed nuclear configuration,

[FORMULA NOT REPRODUCIBLE IN ASCII] (5)

where

[FORMULA NOT REPRODUCIBLE IN ASCII] (6)

The use of the semi-colon in the notation Ψi(x; R) is to denote the explicit dependence of Ψi on the electronic space-spin co-ordinates x and its implicit dependence, along with Ei on the nuclear co-ordinates R. Clearly, the solution of equation (5) for all spatial arrangements of the nuclei will generate an energy (hyper)-surface, the potential energy surface, which governs the motion of the nuclei as obtained by solving the nuclear eigenvalue equation

[FORMULA NOT REPRODUCIBLE IN ASCII] (7)

Equation (7) and its implied assumption of the separability of the electronic and nuclear motions is called the Born-Oppenheimer approximation.

The exact Schrödinger equation of motion, equation (1), may be equivalently stated in a manner which shows the neglected terms arising from the assumption of the product form for the wavefunction, equation (4). The exact eigenstate Ψi(x, R) is expanded in terms of the complete orthonormal set of functions Ψi(x; R) obtained from the solutions of the electronic equation, equation (5), in which case the nuclear wavefunctions χi (R) appear as the coefficients in the expansion. This procedure yields the following infinite set of coupled equations for the χi(R)

[FORMULA NOT REPRODUCIBLE IN ASCII] (8)

From this formal, but exact statement of Schrödinger's equation, which corresponds to the non-adiabatic approximation referred to above, one sees that the approximation made in obtaining equation (7) is to assume that

[FORMULA NOT REPRODUCIBLE IN ASCII] (9)

for all i and j. In another order of approximation, with the expansion of Ψi(x, R) truncated to the single product term as given in equation (4), one obtains in place of equation (8),

[FORMULA NOT REPRODUCIBLE IN ASCII] (10)

for the nuclear eigenfunction χiv(R). The term <Ψi|[??]NΨi> is called the adiabatic correction, a term which, as argued above, should contribute only 1/2000 or less than does Ei(R) to the total energy of the system.

Ko[??]os and Wolniewicz have obtained extremely accurate solutions to Schrödinger's electronic equation, equation (5), for H2 in its ground state over a large range of intemuclear distances. They further obtained solutions to the nuclear motion problem at the adiabatic level of approximation, obtaining estimates of the vibrational quanta which are larger than the experimental values by 0.5 — 0.9 cm-1. The adiabatic correction to the Born-Oppenheimer estimate of the binding energy of [FORMULA NOT REPRODUCIBLE IN ASCII] was found to be 4.95 cm-1. The theoretical estimate of the binding energy of Ha at the adiabatic limit, corrected for relativistic and radiative effects (corrections which are smaller than the adiabatic correction) is in agreement with the most recent experimental value.

The Born-Oppenheimer approximation and the concept of separate nuclear motion on a potential surface breaks down in the case of degeneracy or near degeneracy of the electronic states Ψ i(x; R), resulting in the so-called Jahn-Teller or Renner effects. In these cases, the non-adiabatic terms, <Ψ i|[??]N|Ψj>, can become very large as the energies approach one another. The large mixing of the electronic levels by the nuclear kinetic energy operator implies a strong dependence of the electronic wavefunction on the nuclear motion and results in a complete breakdown of the Born-Oppenheimer approximation.

2 Scope of the Review

The theoretical study of any chemical problem which falls within the domain of the Born-Oppenheimer approximation involves two basic steps: (a) the solution of Schrödinger's electronic equation, equation (5), for one or more electronic states over the range of nuclear configurations demanded by the problem; and (b) solving the nuclear equations of motion on the potential energy surface, or surfaces, thus obtained.

The calculations in part (b) may be of two types; the determination of the nuclear energy levels for bound states of the system, i.e. the quantized vibrational and rotational levels of the system, or the study of the dynamics of the chemical changes described by the surface in terms of quantum reactive scattering or classical trajectory calculations.

It is the purpose of this review to discuss and illustrate the methods presently employed to obtain potential energy surfaces by approximate, but non-empirical solutions to Schrödinger's electronic equation. In addition to discussing the different levels of approximation employed in these ab initio calculations, we emphasize the type of chemical system (in terms of its electronic structure) to which each level of calculation may be expected to yield usable results, i.e. results with acceptable errors or with predictable bounds on the error. Our interest will be primarily in surfaces which have been determined for the prediction and understanding of chemical reactions. This will include a survey of those calculations which have concentrated on determining the reaction path, and the geometry and properties of the system at points on this path, as well as those in which an essentially complete surface has been determined. The latter type of calculation coupled with either classical or quanta! treatments of the nuclear motion on such a surface provides a total theoretical prediction of a chemical reaction. This ultimate objective has been achieved in the case of the H + H2 exchange reaction.

A recent review of the dynamics of bimolecular reactions by Polanyi and Schreiber provides a detailed account of the classical trajectory method for treating the motion of a reactive system on a potential energy surface. Such classical motion studies on three-dimensional surfaces are now commonplace and have yielded a great deal of information regarding the microscopic dynamical behaviour of reactive systems. The majority of these studies, those by Polanyi et al., Bunker et al., Grice et al., and Karplus et al., for example, have employed potential energy surfaces determined by semi-empirical methods. Most notable of these is the 'extended LEPS (London, Eyring, Polanyi, and Sato) function'. Other semi-empirical or purely empirical surfaces have been introduced and used in trajectory studies by Blais, Bunker, and Parr, and by Porter and Karplus. Roach and Child have developed a pseudopotential method for semi-empirical variational calculations of a potential surface.

Quantum calculations of reactive scattering on a potential surface are much more difficult than the corresponding classical trajectory studies and are not as far advanced. Progress in this field has been reviewed recently by Levine, Ross, Karplus, and Light. Fortunately, it appears from the limited comparisons of the quantal and classical studies and from a comparison of the classical results with experimental findings, that the classical approach to the study of reaction dynamics is sufficiently accurate to yield informative and useful results. The classical trajectory studies in particular have resulted in a cataloguing of the properties of potential surfaces with respect to the possible observed types of reaction dynamics. Thus, one desires a knowledge of the principal features of a surface which determine the distribution of energy in the products, whether it appears primarily as vibrational excitation of the newly formed bond, as rotational energy, or as relative translational motion of the product molecules. Information regarding the product energy distribution over vibrational and rotational states is obtained from the i.r. chemi-luminescence method developed by Polanyi (for a recent review, see ref. 28). In addition, crossed molecular beam studies, by measuring the angular distribution as well as the translational energy distribution of the scattered molecules, require that predictions be made regarding the relative probabilities of 'forward', 'back' and angularly uniform scattering of the product molecules. For a recent review of the crossed molecular beam method, the reader is referred to an article by Kinsey. This type of categorization of the properties of potential energy surfaces is well documented in the recent review by Polanyi and Schreiber and will be referred to in discussing the potential energy surfaces obtained for particular reactive systems.

Through the use of the Rayleigh-Ritz variational principle, approximate solutions to Schrödinger's electronic equation, equation (5), may be obtained to any desired degree of accuracy for a many-electron system. Because of the formidable computational problems involved in obtaining these solutions, chemical accuracy (generally taken to imply errors of a few kcal mol-1 or less) has been achieved in only a few two- and three-electron systems. Considerable progress is being made, however, by methods reviewed here, and complete potential surfaces with an accuracy acceptable for trajectory studies for systems with up to 30 electrons and four nuclei now lie within the present bounds of accomplishment. In addition, of course, many questions concerning the electronic mechanisms of particular reactions may be answered by calculations performed at lower levels of approximation.

Of paramount importance in this latter category is the Hartree-Fock approximation. The so-called 'Hartree-Fock limit' represents a well-defined plateau, in terms of its methematical and physical properties, in the hierarchy of approximate solutions to Schrödinger's electronic equation. In addition, the Hartree-Fock solution serves as the starting point for many of the presently employed methods whose ultimate goal is to achieve solutions to equation (5) of chemical accuracy. A discussion of the Hartree-Fock method and its associated concept of a self-consistent field thus provides a natural starting point for the discussion of the calculation of potential surfaces.

3 Hartree-Fock and Self-Consistent-Field Calculations

The orbital approach to the problem of electronic structure reduces a many-electron problem to a corresponding number of one-electron problems. The Hartree-Fock solution represents the best attainable description of the electronic structure of a many-electron system in terms of the one-electron orbital approach.

In the orbital description of electronic structure, the motion of each electron is determined by a spin-orbital, ui a simple product of a space function Φiand a spin function σi

ui(k) [equivalent to]i(rk) σ i(Sk)

where σi may denote either an α or β spin function.

By analogy with a system in which each electron moves independently of the others, the total wavefunction is taken to be a product of the one-electron spin orbitals, one for each electron in the system. To satisfy the additional restrictions of the Pauli principle, and simultaneously to allow for the indistinguishability of the electrons, an antisymmetrized sum of such product functions must be used to approximate the many-electron wavefunction for the system. This is most conveniently done by expressing the wavefunction in the form of a determinant,

[FORMULA NOT REPRODUCIBLE IN ASCII]

since on expansion, a determinant equals a sum of such product terms [as represented in equation (11) by its diagonal element] representing all possible permutations of the n electrons among the n spin-orbitals (n! in number), each term in the sum being multiplied by(+1) or (-1) depending on whether the number of permutations is even or odd. The factor (n!)-1/2 is required to normalize Φ x; R) to unity when the spin-orbitals are orthonormal, i.e.

[FORMULA NOT REPRODUCIBLE IN ASCII] (12)

where the orthonomality obtains for both space and spin functions separately.

The energy of such a determinantal wavefunction, Φ(χ;R), is obtained by the quantum mechanical averaging of the electronic Hamiltonian operator given in equation (6)

[FORMULA NOT REPRODUCIBLE IN ASCII] (13)

or

[FORMULA NOT REPRODUCIBLE IN ASCII] (14)

where

[FORMULA NOT REPRODUCIBLE IN ASCII] (15)

and

[FORMULA NOT REPRODUCIBLE IN ASCII] (16)

The kinetic energy of the electron described by ui together with its energy of interaction with all of the nuclei in the system is represented by hi. The sum over Jij and Kij in equation (14) represents the energy of repulsion between the electrons in the system, the Coulomb integral Jij representing the electrostatic repulsion between the two charge distributions u*i (k)ui(k)] and u*jj(l)uj(l) and the exchange integral Kij correcting this estimate of the repulsion in a manner dictated by the antisymmetry principle.

The forms of the orbitals ui are determined by demanding that the energy obtained in equation (13) be a minimum with respect to arbitrary variations in the orbitals. This procedure leads to a set of n coupled equations, the Hartree-Fock equations, of the form

[FORMULA NOT REPRODUCIBLE IN ASCII] (17)

The Hartree-Fock operator [??](k) is generally expressed as a sum of the two terms, [??](k) = [??] (k) + [??] (k) where [??](k) is the Hamiltonian for motion of the electron in orbital ui(k) in the field of the bare nuclei and where the operation of [??](k) on ui (k)is given by

[FORMULA NOT REPRODUCIBLE IN ASCII] (18)