Mathematical Challenges from Theoretical/Computational Chemistry

Most of the recent progress in theoretical chemistry has come from computation of numerical approximations to the solutions of realistic model problems. This has largely replaced earlier approaches based on analytic solutions to simplified model proble ms. Although there are interesting intellectual problems involved with analytic properties of the exact solutions, progress to date with this approach has had little impact on chemistry.

In quantum chemistry, improved numerical approximations made possible by the greatly decreased cost (in dollars and personnel time) of computing have had an enormous impact. The Journal of the American Chemical Society, which once reluctantly accepted a limited number of articles reporting computational results, now has some computational aspect coupled with experimental results in more than half of the articles published.

The Schrödinger wave equation was introduced earlier, inBox 2.1. The primary method now used to obtain useful results may be described briefly as follows. The time-independent, nonrelativistic Sch rödinger equation for a system of N electrons in the Coulomb field generated by K nuclei may be written in the Born-Oppenheimer approximation (using reduced units to eliminate mass and Planck's constant) as

(1)

H = U

(2)

is the Hamiltonian operator with

(3)

= (r,r,...,r; R,R,...,R) ,

U = U (R,R,...,R)

where r is the position vector for the ith electron and R is the position vector for the ith nucleus with charge Z

This is a partial differential eigenvalue equation in 3N variables. The equation and its solutions are parameterized by the nuclear positions and charges. The eigenvalues U, viewed as functions of R, yield the potential energy surfaces for the molecule in its ground and excited states after including the repulsive internuclear interactions. The lowest eigenvalue, in particular, plays a fundamental role in underst anding the chemical properties of a chemical system.

Equation (1) may be solved subject to either bound state or scattering boundary conditions. For a bound state, the function should give finite values for the integrals < img src="images/symbols/psi.gif">drdr ... for = 1, T, or V, where the asterisk denotes complex conjugation. For a scattering state, the wavefunction may have an exp (i k.r) behavior for large |r|. The w avefunction is further subject to conditions of continuity everywhere and differentiability except at the singularities in V.

The Hamiltonian is symmetric under permutation of the electron coordinates. Consequently, the eigenfunctions can be chosen to form bases for irreducible representations of the permutation group. It has been found empirically that only a few irreducible representations actually correspond to physical reality and the rest are "excluded." This observation is summarized in the Pauli exclusion principle. In the standard labeling of the irreducible representations of the symmetric group by a partition of < !i>N into the sum of integers, only those partitions containing no integers greater than 2 are observed.

In quantum chemistry, this observation is usually implemented by introducing an additional discrete value, called the spin , for each electron. This variable takes only th e values +/-1/2. The permutation operator is then defined to interchange the spin and position of a pair of electrons in the extended function (r,r,...). The Pauli exclusion principle then is simply stated by observing that th e only that occur in nature change sign under every pairwise permutation of the electron coordinates.

Equation (1) has no closed-form solutions for nontrivial chemical systems. The most important problem in electronic structure theory is the rapid construction of useful approximate solutions to this equation. While chemists have made much progress on t his problem, there is always the possibility that a fresh approach by mathematicians would lead to novel approaches. There is less possibility that changing the implementation of the methods now in use would lead to great improvement.

The dominant approach for obtaining approximate solutions at present is based on multiple expansions. First, a basis set is selected in three-dimensional space. Then this is transformed by a similarity transform to an orthonormal set of functions call ed "molecular orbitals." These orbitals are then used to construct Slater determinants, which are the simplest functions in 3N dimensions that obey the Pauli condition of antisymmetry under permutation of the positions and spins of the electrons . Finally, the wavefunction is expanded as a linear combination of Slater determinants. Often, the transformation from basis functions to molecular orbitals is chosen to optimize an approximation to the wavefunct ion from a severely truncated expansion in Slater determinants. The coefficients in the Slater determinant expansion and the approximate eigenvalue U may be computed by the Rayleigh-Ritz variational principle, perturbation theory, cluster expans ions, or some combination of those methods.

All of these methods require formation of matrix elements of the Hamiltonian in the Slater determinant basis,

(4)

where is a Slater determinant and the integration is over 3N dimensions. For a finite basis of orthonormal molecular orbitals, it is possible to replace this equation by a different one that gives the same matrix element

(5)

This is accomplished by utilizing an operator algebra that facilitates manipulation of Slater determinants, in which a is an annihilation operator and a is a creation operator. A correspondence can be established between the function space formed by all possible linear combinations of Slater determinants and the vector space formed by the creation operators so that

(6)

(7)

where, in terms of the orbitals ,

(8)

A major bottleneck in the present approach is the calculation of a large number of six-dimensional integrals g. The choice of basis set is limited to those functions for which these integrals are readily computed.

The second major bottleneck is the actual calculation of U for this second quantized Hamiltonian. Although the values of |U| for various chemical systems vary from 0 to several thousand, the range of Ufor interesting variations in the parameters for a given system is typically only 10 to 10. Most programs determine U to a precision of 10, regardless of the absolute magnitude. Because of the multiple expansions involved, the error in U is often 10 to 10 , but the error is often constant over the interesting range of nuclear position parameters to within +/-10%.

A major advantage of some of these procedures for obtaining U is that first and second derivatives with respect to the nuclear coordinate parameters can be obtained analytically. This has greatly facilitated searching U in this parameter space for mini ma, saddlepoints, and so on. The major disadvantage has been the very steep dependence of computer time on the number of electrons. Evaluating U for one set of parameters typically takes computing time proportional to N, with k ranging from 3 for the least accurate methods to 7 for the best methods now in common use. Finding one eigenvalue of in the full vector space is even more costly, with computer time proportional t o M, where M is the number of basis functions. Consequently, accurate calculations are limited to about 40 electrons and an improvement of a factor of 10 in computer speed will not change this very much.

Although this approach is very productive, it is also limited to small chemical systems. Progress for somewhat larger systems can be made by use of the Hohenberg-Kohn and Kohn-Sham theorems to give a useful density functional theory. These assume that it is possible to find an effective one-body potential so that, by solving the one-electron Schrödinger equation

the eigenfunctions will yield

(9)

where is the exact charge density of the system as defined by Equation (13) below, using the exact wavefunction defined by Equation (1). As before,, , and are parameterized by the nuclear charges and positions. Here is a functional of . From , the lowest potential energy surface U can be formed. The difficulty is that there is only an existence proof for and no systematic constructive procedure is yet known. Nevertheless, mu ch progress has been made by choosing empirically so it will correctly reproduce the properties of simple model systems such as the uniform electron gas and free atoms. The attraction of this method is that the compu tational cost grows only as N. Effort is being made to improve this even further to obtain methods whose costs grow as N or N.

Progress with density functional theory has been rapid in recent years. Direct solution of the Kohn-Sham equation on a three-dimensional grid now is possible, although basis set expansions are still more commonly used. The major limitation in this fiel d is still the lack of a scheme for finding the effective potential that can guarantee the desired accuracy in U for all chemical systems of interest.

For even larger systems, solution of the Schrödinger equation seems hopeless. In this case, U is usually directly approximated by taking advantage of empirically observed near-transferability of parameters between similar chemical systems. This r equires some input from the user of these programs to decide which atoms are bonded. Then, near a local minimum in U, it is possible to approximate U as a sum,

(10)

U U (for displacements from normal bond lengths and bond angles)

+ U (for pairwise nonbonded interactions).

Each pair of nonbonded atoms can be assigned parameters transferred from experimental data to allow a calculation of a U(RAB) energy contribution. Similarly, each chemically distinct type of bond can be assigned an energy for displacement from its usual bond length. Functions constructed in this way have accuracy approaching the desired 10 level needed for chemical prediction. This approach is now used for widely different problems, such as rational drug design, structure of liquids, predicting the shape of moderate-sized organic molecules, and protein folding.

There are many problems in numerical analysis and data handling associated with the present methods. These include the generation and manipulation of large numbers of six-dimensional integrals, finding eigenvalues and eigenvectors of large matrices, an d searching a complicated function for global and local minima and saddle points. There are also important questions about the construction of optimum expansion functions for most rapid convergence. At present, there are no useful error bounds for the e nergy or other properties derived from the wavefunction.

NAS Home Page | NAP Home Page | Reading Room | Report Home Page