This chapter highlights some of the most prominent research challenges from theoretical/computational chemistry that appear to be amenable to attack with the help of reasonable advances in the mathematical sciences. Most of the sections have as their starting point a challenge facing the chemical sciences, which is described in terms that should be accessible to the nonspecialist. The remaining sections—"Molecular Dynamics Algorithms," "Multivariate Minimization in Computational Chemistry," and "Fast Algebraic Transformation Methods"—contain discussions of topics from the mathematical sciences with specific insight into their relevance to theoretical/computational chemistry. Expositions and references are meant to give the mathematical reader sufficient insight and direction to be able into subsequently investigate the topic via deeper reading and especially by discussions with colleagues from the chemical sciences. Note that the topics surveyed in Chapter 3 are also sources of continuing research opportunity, some of which are identified therein.

As an overview of this chapter and a guide for navigating through it, the matrix in Figure 4.1 displays a subjective assessment of the depth of potential cross-fertilization between major challenges from theoretical and computational chemistry and relevant topics in the mathematical sciences. This matrix is based to some extent on intuition because it is an assessment of future research opportunities, not past results. An "H" in the matrix implies an overlap that appears clearly promising, while an "M" suggests that some synergy between the areas is likely. The absence of an H or an M should not be taken to imply that some clever person will not find an application of that technique to that problem at some point.

Two topics, polymers and protein folding, are consciously underrepresented in this report because the mathematical research opportunities related to these topics have been surveyed very recently by other reports from the National Research Council. For a discussion of mathematical opportunities related to the frontiers of polymer science, see pages 153–168 of *Polymer Science and Engineering: The Shifting Research Frontiers* (National Research Council, 1994). A chapter devoted to mathematical research into protein folding may be found in *Calculating the Secrets of Life* (National Research Council, 1995).

National Research Council, 1994, *Polymer Science and Engineering: The Shifting Research Frontiers*, National Academy Press, Washington, D.C.

National Research Council, 1995, *Calculating the Secrets of Life*, National Academy Press, Washington, D.C., in press.

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 problems. Although there are interesting

Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.

Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
4
MATHEMATICAL RESEARCH OPPORTUNITIES FROM THEORETICAL/COMPUTATIONAL CHEMISTRY
Introduction
This chapter highlights some of the most prominent research challenges from theoretical/computational chemistry that appear to be amenable to attack with the help of reasonable advances in the mathematical sciences. Most of the sections have as their starting point a challenge facing the chemical sciences, which is described in terms that should be accessible to the nonspecialist. The remaining sections—"Molecular Dynamics Algorithms," "Multivariate Minimization in Computational Chemistry," and "Fast Algebraic Transformation Methods"—contain discussions of topics from the mathematical sciences with specific insight into their relevance to theoretical/computational chemistry. Expositions and references are meant to give the mathematical reader sufficient insight and direction to be able into subsequently investigate the topic via deeper reading and especially by discussions with colleagues from the chemical sciences. Note that the topics surveyed in Chapter 3 are also sources of continuing research opportunity, some of which are identified therein.
As an overview of this chapter and a guide for navigating through it, the matrix in Figure 4.1 displays a subjective assessment of the depth of potential cross-fertilization between major challenges from theoretical and computational chemistry and relevant topics in the mathematical sciences. This matrix is based to some extent on intuition because it is an assessment of future research opportunities, not past results. An "H" in the matrix implies an overlap that appears clearly promising, while an "M" suggests that some synergy between the areas is likely. The absence of an H or an M should not be taken to imply that some clever person will not find an application of that technique to that problem at some point.
Two topics, polymers and protein folding, are consciously underrepresented in this report because the mathematical research opportunities related to these topics have been surveyed very recently by other reports from the National Research Council. For a discussion of mathematical opportunities related to the frontiers of polymer science, see pages 153–168 of Polymer Science and Engineering: The Shifting Research Frontiers (National Research Council, 1994). A chapter devoted to mathematical research into protein folding may be found in Calculating the Secrets of Life (National Research Council, 1995).
References
National Research Council, 1994, Polymer Science and Engineering: The Shifting Research Frontiers, National Academy Press, Washington, D.C.
National Research Council, 1995, Calculating the Secrets of Life, National Academy Press, Washington, D.C., in press.
Numerical Methods for Electronic Structure Theory
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 problems. Although there are interesting

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
Adaptive and multiscale methods
H
M
M
M
M
M
M
M
Special bases
H
M
M
H
Differential geometry
M
M
M
M
H
H
Functional analysis
H
M
H
H
H
M
M
M
Graph theory
M
M
M
H
M
M
M
Group theory
H
M
M
M
M
M
M
Optimization
H
H
H
H
M
H
H
M
H
H
H
Numerical linear algebra
H
H
M
M
M
Number theory
M
M
M
Pattern recognition
M
M
H
H
M
H
H
H
Probability and statistics
H
H
M
H
M
H
H
H
Several complex variables
H
H
M
M
M
M
Topology
M
M
H
H
M
M
H
H
H
Dynamical systems
M
H
H
M
H
M
Quantum electronic structure
Molecular mechanics
Condensed-phase simulations
Density functionals
N-representability
Design of molecules
Construction of potential energy functions
Gas-phase dynamics
Polymers
Topography of potential energy surfaces
Biological macro-molecules (including protein folding)
FIGURE 4.1. Subjective assessment of depth of potential cross-fertilization between major areas of the mathematical sciences and theoretical and computational chemistry. An "H" implies an overlap that is clearly promising and an "M" suggests some synergy is likely.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
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, in Box 2.1. The primary method now used to obtain useful results may be described briefly as follows. The time-independent, nonrelativistic Schrö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
where
is the Hamiltonian operator with
where ri is the position vector for the ith electron and Ri is the position vector for the ith nucleus with charge Zi.
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 Uk, viewed as functions of Ri, 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 rule in understanding 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 ψ*qyd3r1d3r2 ... for θ = 1, T, or V, where the asterisk denotes complex conjugation. For a scattering state, the wavefunction may have an exp (ikr) behavior for large |r|. The wavefunction 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 N into the sum of integers, only those partitions containing no integers greater than 2 are observed.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
In quantum chemistry, this observation is usually implemented by introducing an additional discrete value, called the spin ξi, for each electron. This variable takes only the values ±1/2. The permutation operator is then defined to interchange the spin and position of a pair of electrons in the extended function ψ(r1,ξ1, r2ξ2, . . .). The Pauli exclusion principle then is simply stated by observing that the 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 this 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 called "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 wavefunction 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 expansions, or some combination of those methods.
Slater determinants are determinants whose elements are n distinct orbitals for n distinct electrons.
In the Hartree-Fock approximation an exact Wavefunction is replaced with an antisymmetrized product of single-particle orbitals (i.e., a Slater determinant).
All of these methods require formation of matrix elements of the Hamiltonian in the Slater determinant basis,
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
This is accomplished by utilizing an operator algebra that facilitates manipulation of Slater determinants, in which ai is an annihilation operator and ai 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
That is, the Slater determinant, ΦI, in which the orbital labels i1 . . . iN appear, corresponds to the abstract element |ΦI>. The operator , called the "second quantized Hamiltonian," takes the form

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
where, in terms of the orbitals φi,
A major bottleneck in the present approach is the calculation of a large number of six-dimensional integrals gijkl. 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 U for interesting variations in the parameters for a given system is typically only 10-1 to 10-3. Most programs determine U to a precision of 10-6, regardless of the absolute magnitude. Because of the multiple expansions involved, the error in U is often 10-1 to 10-2, 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 minima, 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 Nk, 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 to MN+4, 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
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 εi 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, much 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 computational cost grows only as N3. Effort is being made to improve this even further to obtain methods whose costs grow as N2 or N.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
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 field 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 requires 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,
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-3 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, and 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 energy or other properties derived from the wavefunction.
The N-and V-Representability Problems
Conventional approaches for computing the solution to the wavefunction have a strong dependence on the number of electrons. Therefore, searches are constantly under way for methods of comparable accuracy with a better scaling. In this regard, people have considered methods based on density matrices, density functionals, Monte Carlo diffusion equations, effective core potentials, and so on. In particular, because the energy of an atom or molecule is a linear function of the density matrix and the one-and two-body distribution functions derived from it, density matrix methods raise the hope that one could dispense with computing the associated 3N-dimensional Wavefunction and deal with simpler three-dimensional density functions. A number of mathematical issues are raised by the attempt to reformulate computational chemistry in terms of particle distribution functions instead of wavefunctions.
The N-body distribution function is given simply as the product of the wavefunction and its complex conjugate,
Here, x symbolizes the collection of coordinates x1, x2, . . ., xN, describing the positions ri and spins ξi of all N particles. In the chemical literature, this quantity is usually called the N-body density

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
matrix. This may be averaged over an ensemble with an arbitrary set of probabilities (weights) wi to give the ensemble distribution,
Because of permutational symmetry, all identical particles enter this distribution in an equivalent way. From the N-body distribution, one-and two-body position distribution functions can be obtained by integration over the other coordinates:
The distribution function P(x1,x2) is also called the two-body reduced density matrix. Similarly, the density matrix may first be Fourier transformed and then used to derive the one-body momentum distribution
For most cases of interest, knowledge of the one-and two-body position distribution functions and the one-body momentum distribution function suffices to determine the energy. For the Hamiltonian H with given external potential Vex(ri) and given two-body interaction potentials g(rij),
the energy is given by
There are a number of outstanding mathematical problems associated with these distributions. For many potentials, the behavior of ρ( r) is known near the singularities of the potential. Similarly, the behavior of γ(r12) is known for a Coulomb interaction near the singularities, and the form of π(p) is known for large and small p. The N-representability problem consists of finding the conditions on this set of three functions such that they could all come from the same N-body density matrix. Many inequalities are known, but no general solution has been obtained. The V-Representability problem consists of determining the further restriction imposed by considering only those distributions that can come from a wavefunction that is the eigenfunction of some H for a fixed g.
The two-body reduced density matrix Γ(2) can be integrated to give the one-body reduced density matrix

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
Density matrices may be treated as integral kernels and factored in terms of their eigenfunctions
where
These functions fi are known as natural orbitals. It is conjectured, but not proved, that they form a complete set when Γ(1)(1) is derived from an exact eigenfunction of H with Coulomb interactions. The complementary functions
are certainly not complete in the 3(N — 1) dimensional space. The extended Koopmans theorem claims, however, that if ψ is the exact eigenfunction of H for an N-electron system, then the ground state wavefunction of the N — 1 electron system may be expanded exactly in the set of Fi. Both ''proofs" and "disproofs" of this conjecture exist in the literature (see Morrison, 1993, and Sundholm and Olsen, 1993, for a recent exchange of opinions and a list of relevant earlier papers).
Because the energy is a linear function of the one-and two-body distribution functions, the variational minimum will lie on the boundary determined by the N-representability conditions. Unfortunately, only an incomplete set of necessary conditions are known, but these are already so complex that further work in this area has been abandoned by chemists.
In density functional theory, only the density ρ (r) is used. Hohenberg and Kohn showed that the energy is a functional of this density for a fixed two-body potential g(rij). Density functional theory has become of practical importance, but progress is hindered by lack of knowledge of the properties of the functional and lack of a systematic procedure for constructing a convergent sequence of functionals. Some past work in this field has been summarized in several monographs (Davidson, 1976; Dahl and Avery, 1984; Kryacho and Ludena, 1989; March, 1989; Parr and Yang, 1989; Sham and Schlinter, 1989; Gadre and Pathah, 1991).
References
Dahl, J.P., and J. Avery, 1984, Local Density Approximations in Quantum Chemistry and Solid State Physics, Plenum Press, New York.
Davidson, E.R., 1976, Reduced Density Matrices in Quantum Chemistry, Academic Press, New York.
Gadre and Pathah, 1991, Advances in Quantum Chemistry 22:211.
Kryacho, E.S., and E.V. Ludena, 1989, Density Functional Theory of Many-Electron Systems, Kluwer Press, Dordrecht.
March, N.H., 1989, Electron Density Theory of Atoms and Molecules, Academic Press, New York.
Morrison, R.C., 1993, The oxactness of the extended Koopmans Theorem—A numerical study—Comment, J. Chem. Phys. 99:6221.
Parr, R.G., and W. Yang, 1989, Density Functional Theory of Atoms and Molecules, Oxford University Press,

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
Oxford.
Sham, L.J., and M. Schlinter, 1989, Principles and Applications of Density Functional Theory, World Scientific, Teaneck, N.J.
Sundholm, D., and J. Olsen, 1993, The exactness of the extended Koopmans theorem—A numerical study—Response, J. Chem. Phys. 99:6222.
Melding of Quantum Mechanics with Simpler Models
A daunting challenge for the future is to accurately model chemical reactions in phases and at the active site of enzymes. An ability to do so would be of great importance in designing new biological catalysts as well as fully understanding the chemical mechanism of those that already exist. This would be of significant technological as well as scientific importance. One could imagine that many new molecules could be made and made much more efficiently by such catalysts.
Methods to do this in an approximate way have been available since 1976 (Warshel and Levitt, 1976) and have involved using the Schrödinger equation (Equation 1 or a variant or empiricized form of it) for the parts of the system where bonds are being made or broken and thus the electronic structure is changing, combined with representations such as Equation (10), which assume transferable electronic structure, for the remaining atoms of the system. Typically, the number of atoms for which Equation (1) must be solved is much smaller, of the order of 20 to 30, than the number in the whole system, which for the chemical reactions mentioned is typically more than a thousand.
In fact, some exciting results for simple reactions involving organic molecules in different solvents have been achieved (Blake and Jorgensen, 1991). In these cases, one has solved Equation (1) to high accuracy for a simple reactive pathway in vacuo and then, employing these energies, has used free energy calculation methods to evaluate the solvation free energy of different structures along the reactive pathway. This is in some sense a proof of concept for the combined application of Equations (1) and (10) because impressive agreement with experiment has been achieved in these simple, welldefined cases.
For more complex cases, such as enzyme reactions, the reaction pathway might involve many steps, and some of the reacting groups are chemically bonded to the protein, thus requiring some additional technical challenges in simulating the atoms at the junction between those that are participating in the chemical reactions and those that are not. In addition, one might have to consider many conformations of the enzyme and its substrate and accurately represent their relative energy by using the energy function of Equation (10), all the while considering the electronic energy (Equation 1) and the relative total free energy of the system.
As noted above, progress on this problem has been made when employing much simplified representations of the electronic structure of the system, which enable the solution of equations such as Equation (1) for the few "quantum mechanical" atoms as rapidly or more so than the classical molecular dynamical equations of motion, using Equation (10) as a potential energy (Field et al., 1990; Warshel, 1991).
These methods use semi-empirical or empirical valence bond approximations to solve Equation (1). Although these methods are not highly accurate, the use of non-empirical quantum mechanical

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
BOX 4.1 Electronic Phase Transitions
The electronic properties of extended bulk materials such as metals and superconductors have been studied extensively by physicists. However, as the field has progressed from simple elemental materials to more complex synthetic ones, computational chemistry has come to the fore in addressing these problems. The "exotic materials" include the organic superconductors such as Bechgaard salts, polyacetylene, the fullerenes, and the high-temperature superconductors based on copper oxides. In addition, many new amorphous materials are technologically interesting in their electronic properties. For instance, electronic conduction in polymers, such as polyvinylcarbazole, is essential to some xerographic processes. The chemical complexity of these systems puts a premium on understanding the fundamental physics in new ways that do not usually rely on the simple symmetries present in elemental materials. All of these systems exhibit electronic phase transitions as the chemical composition or doping is changed: their electronic states change qualitatively. Sometimes this transition is from being an insulator to a metal, sometimes from a metal to a superconductor or to some complex magnetically ordered structure.
One of the simplest electronic phase transitions is the transition between extended and localized states of a single electron moving in a random potential. Even though this problem is at the heart of the study of the electronic properties of any disordered material, the traditional methods of simply combining rough semiquantitative theories and experiment have been insufficiently powerful to resolve all of the important issues. One reason is that real materials have many physical influences in addition to disorder, such as interactions of vibrations and interactions between different electrons in the same material. Some mathematical work, such as rigorous theorems related to one-dimensional Anderson localization, has already helped in understanding this problem. On the other hand, loopholes in these theorems have been uncovered when the disorder is of a special, correlated kind. For instance, certain 1-D systems do not have only localized states as the simple theorems had indicated. These unusual sorts of extended states arise in systems with certain kinds of correlated disorder or with quasi-periodic Hamiltonians.
Interestingly, far from being a mathematical curiosity, these exceptions to the simple theorems about one-dimensional localization seem to be at the heart of understanding the behavior of materials such as polyacetylene. In two-or three-dimensional materials, experiment has amply demonstrated the existence of both extended and localized states. However, there are still numerous controversies about the applicability of simple phase transition ideas to these electronic phase transitions. Are they described by the usual scaling phenomenology of ordinary thermal phase transitions?
It has been argued that, in fact, such descriptions may break down because the wavefunctions at the transition are multifractal. Thus, the study of these electronic phase transitions has much in common with problems addressed in quantum chaos, where the structure of the wavefunctions needs to be understood in a statistical way. The interacting electron systems and their phase transitions also carry mathematical challenges. In one dimension, the interactions between electrons again can cause them to behave as if they were insulating. These one-dimensional many-electron problems lead to exactly solvable Schrödinger equations. The exact integrability of these classical models is connected with conformal invariance and the existence of solitons of nonlinear partial differential equations in one dimension. The question of whether such electronic interaction effects give rise to new phases for higher-dimensional systems doubtless has connections with the problem of solitons and exact integrable systems in higher dimensions—a problem that has attracted many in the area of partial differential equations.
Finally, a large number of interesting phase transitions occur in disordered systems that also have interactions. These include the Kondo effect in which isolated electronic sites behave as if they have spins that can be quenched, as well as the exotic spin glass phases that have proved useful at least as analogies in many other areas of chemistry and biology. At the moment, the most useful theory of these system is based largely on the use of the unrestricted Hartree-Fock approximation. There are numerous mathematical questions connected with the Hartree-Fock problem for such disordered systems. Many of the ideas invented by Lieb in his proof of the stability of matter may have practical use here. Again, the rigorous mathematical analysis can be of significant help in understanding whether or not certain approximations can be used confidently in elucidating the qualitative physics. In addition, the same issues will arise when quantitative calculations are contemplated for specific materials.
The study of electronic phase transitions in dusters will make these issues even more pronounced. The usual theory of phase transitions concentrates on infinite systems, and only dominant thermodynamic contributions are computed. For clusters, the finite size effects and the asymptotics of the approach to the infinite limit become important mathematical problems.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
methods for systems of 20 to 30 atoms (the traditional ab initio approach) requires 104 to 106 more computations than the semi-empirical or empirical approaches and much more time than it takes to solve the classical equations of motion for the rest of the system (Singh and Kollman, 1986).
To elaborate, calculation of the free energy of a complex chemical system by using classical molecular dynamics (Kollman, 1993) requires one to calculate the energy of the system and its gradient with respect to all the 3N coordinates. This can be done for noncovalent processes (those using only Equation 10) quite efficiently because the energy function in Equation (10) is very simple and its derivatives are quick and easy to evaluate. When one adds quantum mechanical (bond making or breaking) effects via Equation (1), in order to make the calculation of the free energy tractable, one must be able to evaluate the quantum mechanical energy and its gradient for the few quantum mechanical atoms as rapidly as the classical molecular mechanical energy and gradient for the thousands of atoms in the remainder of the system. This can be done by using simpler empirical (Warshel, 1991) and, to a reasonable approximation, semi-empirical (Field et al., 1990) quantum mechanical methods, but not with the first principle ab initio methods.
Density functional methods (Labanowski and Andzelm, 1991), particularly the divide-and-conquer strategy (Yang, 1991), show promise in leading to accurate and rapid solutions of Equation (9) for the electronic structure, but they are still a long way from being fully developed, so one cannot tell how efficient and useful they will be in this regard.
Thus, accurate simulation of chemical reactions at the active sites of macromolecules will likely require significant progress in the conformational search problem, even if one considers only the active site of the enzyme. As should be emphasized, the "conformational search problem" requires one not only to consider many conformations, but also to rank their relative free energy in solution. On top of this, one places the problem of accurate and very rapid electronic structure calculations. The above problems are very challenging conceptually, practically, and computationally.
References
Blake, J.F., and W. Jorgensen, 1991, Solvent effects on a Diels-Alder reaction from computer simulations, J. Am. Chem. Soc. 113:7430–7432.
Field, M.J., P.A. Bash, and M. Karplus, 1990, A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations, J. Comput. Chem. 11:700–733.
Kollman, P., 1993, Free energy Calculations—Applications to chemical and biochemical phenomena, Chem. Rev. 93:2395–2417.
Labanowski, J., and J.W. Andzelm, eds., 1991, Density Functional Methods in Chemistry, Springer-Verlag, New York.
Singh, U.C., and P.A. Kollman, 1986, A combined ab initio QM/MM method for carrying out simulations on complex systems: Application to the CH3Cl + Cl- exchange reaction and gas phase protonation of polyethers, J. Comput. Chem. 7:718–730.
Warshel, A., 1991, Computer Modeling of Chemical Reactions in Enzymes and Solutions, John Wiley, New York.
Warshel, A., and M. Levitt, 1976, Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme, J. Mol. Biol. 103:227–249.
Yang, W., 1991, Direct calculation of electron density in density-functional theory-implementation for benzene and a tetrapeptide, Phys. Rev. A. 44:7823–7826.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
Here Za and Zb are the nuclear charges, rai is the distance separating nucleus a from electron i, and R is the distance between nuclei a and b.
If one assumes that the electronic Schrödinger equation,
can be solved exactly for the complete set of eigenfunctions ψk (r,R) and eigenvalues Ek(R), then the total Schrödinger equation for nuclear and electronic motions,
can be solved by expanding Ψ as follows,
Here r represents all the coordinates of the electrons and k is the set of electronic quantum numbers. This leads to a set of equations for the functions φk(R) that determine the nuclear motion of the system,
An a sterisk denotes complex conjugate. The quantities defined by Equation (28) give rise to velocity-

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
dependent forces on the nuclei. For real ψl, however, the diagonal term vanishes.
In practice, Equation (27) is extremely difficult to solve and various approximations to it are introduced. The BO approximation corresponds to neglecting all of the coupling terms. Equation (27) then becomes a Schrödinger equation for nuclear motion,
In this approximation, the El(R) determined from Equation (24) become the potential energy for the nuclear motion.
The adiabatic approximation corresponds to neglecting all nondiagonal terms in Equation (27), which results in a Schrödinger-type equation for nuclear motion,
where the potential for nuclear motion is
The diagonal elements Ell'(R) are effectively a correction to the potential energy due to the coupling between the electronic and nuclear motions.
The adiabatic approximation gives the best potential energy function. As defined here the adiabatic approximation to the energy is an upper bound to the true energy since it can be expressed as the expectation value to the correct Hamiltonian for the molecule evaluated with an approximate wave function.
The nonadiabatic approximation corresponds to consideration of the nondiagonal as well as the diagonal elements Elk'(R). This is extremely difficult to carry out but has been accomplished for H2 by using variational basis set (Kolos and Wolniewicz, 1960) and quantum Monte Carlo approaches (Traynor et al., 1991; see also Ceperly and Alder, 1987).
For processes occurring in excited states, it appears clear that for molecules more complicated than H2 alternative approaches must be found. This is an area demanding fundamental improvements if breakthroughs are to be achieved.
References
Ceperley, D.M., and B.J. Alder, 1987, Phys. Rev. B 36:20–92.
Hirschfelder, J.O., and W.J. Meath, 1967, Adv. Chem. Phys. 12:1.
Kolos, W., and L. Wolniewicz, 1960, Rev. Mod. Phys. 35:13–23.
Traynor, C.A., J.B. Anderson, and B.M. Boghosian, 1991, A quantum Monte-Carlo calculation of the ground-state energy of the hydrogen molecule, J. Chem. Phys. 94:3657–3664.
Evaluation of Integrals with Highly Oscillatory Integrands: Quantum Dynamics with Path Integrals
The solutions of many quantum dynamics problems can be formulated in terms of path integrals. Indeed, there is a formulation of quantum mechanics (less familiar than the Schrödinger wave function formulation, the Heisenberg matrix formulation, or the Dirac Hilbert space formulation) based

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
entirely on the use of path integrals (Feynman and Hibbs, 1965). Path integrals are closely related to some of the mathematical approaches used to describe classical Brownian motion. A prototypical mathematical problem and a class of methods for solving it are described first; then a more general discussion of the problem is presented.
Prototypical Problem
This discussion of the prototypical problem follows closely the discussion of Doll et al. (1988). Consider the integral I(t) = dx ρ(x) eitf(x) where t is a real parameter corresponding to the time, x is the set of coordinates for a many-dimensional space, and ρ(x) and f(x) are functions in that space. The function ρ(x) is positive semidefinite and can without loss of generality be taken as normalized to unity. Thus, it can be regarded as a probability density in the space. The function f(x) is real (for this prototype problem).
For problems of interest, the dimensionality of the space can be very large (from 1 to a hundred to thousands), and the complexity of the function precludes analytic integration. For any particular x, ρ(x) and f(x) can be evaluated numerically, and typically there is much more computation required for the evaluation of ρ than f.
The range of the x integration may be bounded, but more often it is unbounded. For typical problems of interest, however, integrals of the form dx ρ(x) (f(x))n exist for all positive integers n.
The goal is usually to evaluate the integral for a specific value of t (or for a set of values of t) by making use of calculated values of ρ and f at appropriately chosen points and to estimate the accuracy of the answer. In some cases, the desired quantity is the area under the function
or perhaps the Fourier transform
Techniques for evaluation of such integrals of I would be extremely worthwhile, but this discussion focuses on the evaluation of I(t) for specific values of t.
Discussion of the Problem
A standard way of attacking this problem would be a Monte Carlo integration scheme in which points x are generated with a probability distribution of ρ(x) and f is evaluated at these points. Then
where N points are generated and xn is the nth point. For small values of t, this can be a practical procedure. For large values of t, however, this procedure is extremely inefficient because of the large amount of cancellation among the various terms.
If t is large and f is not a constant, it is clear that the integrand is highly oscillatory; thus, small regions of space where f is varying will tend to contribute small net amounts to the integrand for large t. Nevertheless, individual points in such a region will contribute to the sum an exponential whose magnitude is unity. Cancellation among the various terms will give an accurate and small answer only when each such small region is sampled many times in the Monte Carlo evaluation.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
The regions of the space in which df(x)/dx is zero or small will make the major contributions to the integral for large t. Thus, a scheme that concentrates the Monte Carlo sampling at and near such regions would be desirable. Such schemes are useally called "stationary-phase" methods, since at such points the phase of the exponential is stationary.
Stationary-Phase Monte Carlo Methods
A variety of sampling methods have been developed that focus on the stationary-phase points (Doll, 1984; Filinov, 1986; Doll et al., 1988). Here the approach of Doll et al. (1988), which is similar to the others, is outlined.
The integral of interest can be rewritten identically as
where
and P(y) is an arbitrary normalized probability distribution. (This result is exact if the range of integration is over all positive and negative values of the integration variables.) The function P( y) is typically chosen to be a Gaussian distribution centered at y = 0. D(x,t) is called a "damping function," because for such a typical choice of P(y), D(x,t) is largest where f is stationary and hence it damps the integrand in Equation (29) away from the stationary points.
An evaluation of D(x,t) poses the same problems as evaluation of I(t). However, for specific choices of P, various approximations for D can be constructed and used as the basis for a more efficient sampling method for the evaluation of I. For example, if P(y) is a narrow Gaussian function with a maximum at y = 0, then for large t an approximation can be constructed by replacing ρ(x-y)/ρ(x) by unity and performing a Taylor expansion of the exponent. One then obtains an approximation for D called the "first-order gradient approximation":
Then
This integral, which is the major contribution to the correct answer, is generally more amenable to evaluation by sampling, in this case by using ρ(x)D0(x,t) as the probability distribution. A complete method would require the development of procedures to estimate the correction terms that must be added on the right to make this an equality (see Doll et al., 1988, for additional details). The essential idea of these procedures is that use of a precise approximate formula for D automatically includes much of the cancellation responsible for decreasing the value of I.
The stationary-phase Monte Carlo methods provide some useful ideas for formulating tractable sampling methods for solving some problems of interest. It would be worthwhile to develop better versions of these methods both for I(t) and for its Fourier transform.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
Alternative Approaches to the Prototype Problem
In many cases of interest, the function I(t) is known to be an analytic function of t at the origin. This knowledge may be of use in developing approximate evaluation techniques based on analytic continuation. Most analytic continuation methods are based on a more detailed formulation of the problem of interest that is tied closely to the physical nature of the problem. But even for this simple formulation of the prototype problem, one might ask whether the fact that I(t) is an analytic function of (complex) t for some region that includes the origin suggests any alternate ways of calculating I(t) accurately.
For example, one might consider the use of Padé approximants (i.e., ratios of polynomials; Baker and Gammel, 1970), which are often used for approximate analytical continuation of analytic functions, especially meromorphic functions. Expanding the exponential in I(t) gives
where
The first few Taylor series coefficients can then be estimated by sampling on ρ(x) in the usual way, but the difficulty of making accurate estimates grows as n increases. These coefficients might then be used to fit I(t) with a Padé approximant.
This approach raises a number of mathematical questions.
Are there any sampling techniques that would be especially effective in evaluating the coefficients for large n?
Does the statistical error in the αn make it inappropriate to use Padé approximants to fit the function? Are there better alternatives?
If the Padé approximant method is valid, how is the uncertainty in I(t) related to the statistical uncertainty of the individual αn?
Does the fact that the integral expression for I(t) is dominated by stationary points in the x space give any insights into how to evaluate αn or how to construct appropriate Padé approximants?
Other Formulations and Solutions of the Basic Problem
The nature and dimensionality of the x space and the specific form of the ρ(x) function are determined by the nature of the quantum mechanical problem to be solved. In general, however, the coordinates of x are Cartesian (or other) coordinates for a mechanical system. Similarly, the form of the function f(x) is determined by the nature of the problem. Various analytic continuation techniques have been applied to attack specific special cases of the prototype problem. The question then arises as to whether these are isolated solutions of specific problems or special cases of an underlying, not yet discovered, general method for solving problems of this type. This is not the place to delve into the details of specific problems; instead, some of the relevant mathematical principles involved are highlighted. For additional discussion of related problems, the reader is referred to Doll (1984), Makri (1991), and Wolynes (1987).

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
Analytic Continuation in Time. In some problems, the goal is to evaluate a function of time t of the form
More generally, there is a function of a complex variable z of the form
that is analytic in a region including the positive real axis and part of the positive imaginary axis, and for which h(t) = lim 0+H( + it). The function g(z,x) is an analytic function of all its arguments, and f(x) is real (for real x]). In many cases, the formulation in Equation (31) is required to make the problem well defined, because g(it, x) is purely imaginary and the integrals in Equation (30) are not absolutely convergent (for an infinite range of integration).
The calculation of H(z) for real z presents a tractable sampling problem because g(z,x) is real for real z. Then, sampling a set of points distributed with a probability density proportional to g(z, x) and calculating f at those points allow H to be calculated for real z. In fact, as long as z has a positive real part, the problem of evaluating H can be converted to a well-defined sampling problem, because g(z,x) - for large positive and negative values of the integration variables. Then g(z,x) can be used as a probability distribution, and the right side of Equation (31) can be expressed as a ratio of two averages over this distribution:
For complex values of z, however, this sampling problem may become intractable because the quantities to be averaged, exp [ig(z,x)]f( x) and exp [ig(z,x)], are complex and oscillatory. This can pose the same difficulty as that discussed for the prototypical problem above. This difficulty becomes acute in the limit that z becomes purely imaginary, because the "probability distribution function" approaches unity and is not normalizable.
A variety of analytic continuation methods have been used to solve problems of this type. One strategy is to perform tractable sampling calculations to obtain H(z) for real z and then use Padé approximant techniques to estimate H as z approaches the imaginary axis (Miller et al., 1983; Thirumalai and Berne, 1983; Jacquet and Miller, 1985; Yamashita and Miller, 1985; Makri and Miller, 1989). These techniques motivate the same sorts of mathematical questions as those listed above for a different approach to analytic continuation.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
Deformation of the Integration Contours. Another technique has been used in certain situations in which the integrand is an analytic function of all the variables in x. The integration is from- to along the real line for each of these variable. Deformation of these contours of integration and evaluation of the resulting integrals by sampling methods can be an effective technique if the deformed contours come near or pass through stationary-phase points of the integrand (Chang and Miller, 1987; Mak and Chandler, 1991). In one situation in which this was used (Mak and Chandler, 1991), the contour was described by a set of parameters that were varied during the calculation to maximize some measure of the success of the sampling process being used to calculate the integral.
Most General Formulation. Perhaps the most general formulation of these problems is the following one. Given a set of complex variables z1, ..., zn and a function of these variables g(z1, ...., Zn) that is an analytic function of all these variables for a domain that includes the real axis for each of them, calculate the following complex function of m complex variables:
The method should work for highly oscillatory functions and should take into account the fact that each evaluation of g and the location of stationary-phase points are expensive. Here z1, . . . , zm are analogous to the t in previous examples, and zm+1, ..., zn are the integration variables.
References
Baker, Jr., G.A., and J.L. Gammel, 1970, The Padé Approximant in Theoretical Physics, Academic Press, New York.
Chang, J., and W.H. Miller, 1987, Monte-Carlo path integration in real-time via complex coordinates, J. Chem. Phys. 87:16–48.
Doll, J.D., 1984, Monte Carlo Fourier path integral methods in chemical dynamics, J. Chem. Phys. 81:35–36.
Doll, J.D., Thomas L. Beck, and David L. Freeman, 1988, Quantum Monte-Carlo dynamics—The stationary phase Monte-Carlo path integral calculation of finite temperature time correlation-functions, J. Chem. Phys. 89:5753–5763.
Feynman, R.P., and A.R. Hibbs, 1965, Quantum Mechanics and Path Integrals, McGraw-Hill, New York.
Filinov, V.S., 1986, Calculation of the Feynman integrals by means of the Monte-Carlo method, Nucl. Phys. B 271:717–725.
Jacquet, R., and W.H. Miller, 1985, J. Phys. Chem. 89:21–39.
Mak, C.H., and D. Chandler, 1991, Coherent-incoherent transition and relaxation in condensed-phase tunneling systems, Phys. Rev. 44:2352–2369.
Makri, N., 1991, Feynman path integration in quantum dynamics, Computer Physics Communications 63:389.
Makri, N., and W.H. Miller, 1987, Monte-Carlo integration with oscillatory integrands: implications for Feynman path integration in real time, Chem. Phys. Lett. 139:10–14.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
Makri, N., and W.H. Miller, 1989, Exponential power series expansion for the quantum time evolution operator, J. Chem. Phys. 90:904–911.
Miller, W.H., S.D. Schwartz, and J.D. Tromp, 1983, Quantum mechanical rate constants for bimolecular reactions, J. Chem. Phys. 79:488–990.
Thirumalai, D., and B.J. Berne, 1983, On the calculation of time correlation functions in quantum systems: path integral techniques, J. Chem. Phys. 79:502–933.
Wolynes, P.G., 1987, Imaginary time path integral Monte Carlo mute to rate coefficients for nonadiabatic barrier crossing, J. Chem. Phys. 87:65–59.
Yamashita, K., and W.H. Miller, 1985, "Direct" calculation of quantum mechanical rate constants via path integral methods: Application to the reaction path Hamiltonian, with numerical test for the H + H2 reaction in 3D, J. Chem. Phys. 82:5475–5484.
Fast Algebraic Transformation Methods
Numerous problems require rapid evaluation of a product of a matrix and a vector,
This arises in the transformation from one basis in a linear space to another. If computed as written, this transformation computes the n values of the transform in an amount of work proportional to n2. However, there are both exact and approximate methods that can compute such a transform much more efficiently in a variety of cases.
One of the most famous is the Fourier transform, which is frequently used to transform from physical space to frequency space. The fast Fourier transform (FFT) is a well-known algorithm for computing the n values of the transform in an amount of work proportional to n log2 n. The FFT is an exact method, not an approximation: in exact arithmetic, it computes the sum exactly by a reordering of the summands. Although this algorithm is remarkably successful, recent advances have been made by further exploiting its algebraic structure.
Efficient implementation of large multidimensional FFTs on modern computers is memory bound This can be seen from the following timings for a Sun Sparcstation 10 performing an FFT on a two-dimensional image of size 2K by 2K words using the standard FFT algorithm. When carried out on a machine with 64 megabytes of random access memory (RAM), the algorithm executes in 66 seconds. However, on a machine with only 32 megabytes of RAM, the algorithm takes more than two hours, due to the fact that the 2K by 2K dataset cannot fit in RAM and data must be moved in and out of memory during the calculation. Recently, J.R. Johnson and R.W. Johnson (private communication, June 1994) reported that a truly multidimensional version of the FFT algorithm has been developed that accomplishes the same computation in only 75 seconds on a 32-megabyte machine and 47 seconds on a 64-megabyte machine.
A recent example of approximate methods is the fast multipole expansion technique, which takes advantage of special properties of the matrix (aij). This technique has been successfully applied to evaluate Coulombic potentials in molecular dynamics simulations (Ding et al., 1992). This idea has been extended (Draghicescu, 1994) to a broad class of matrices using a general approximation procedure. These more general techniques would apply directly to the non-Coulombic potentials

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
arising in periodic problems through the use of Ewald sums.
A closely related algebraic operation occurs in the evaluation of multilinear forms such as those that appear in density functional methods for calculating ground electronic states in quantum electronic structures. A bilinear form B((zi), (xi)) can always be written in terms of a matrix:
As written, the bilinear form would require O(n2) operations to evaluate. It can be computed instead by using the intermediate vector (yi) defined above, calculation of which requires O(n) operations, as
This way of evaluating a bilinear form can greatly reduce overall computation time, depending on how efficiently yi can be computed.
In density functional methods for calculating ground electronic states in quantum electronic structures, it is desired to approximate multilinear forms including integrals of the form
where
(Here, μi denotes a multi-index, say (pi,qi, ri) so that
In particular, such a representation arises via an expansion
Hence,
The evaluation of ρ at a single point provides an example of the alternatives for evaluating a bilinear form. Using the representation involving the matrix (Cijrequires n2 operations, whereas the expression involving the χ's requires m·n operations. If m n, the latter approach will be more efficient.
The above integral can be viewed as a quadrilinear form involving the coefficients Multilinear forms can be evaluated in a variety of ways. It is tempting to represent them in terms of a precomputed tensor (a matrix for a bilinear form). Recently, Bagheri et al. (1994) have observed that it can be more efficient, in both time and memory, not to precompute expressions.

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
The original integral can be written
By using a suitable quadrature rule, this could be approximated via
Let Q = q2 denote the total number of quadrature points; the quadrature approximation can be calculated in an amount of work proportional to Qnm.
On the other hand, the original integral can be expressed as
where computation of the intermediate quantities
has been introduced, which adds only O(n2m) work. The summation over i, i', j, j' requires an amount of work and storage proportional to n4 to evaluate. It would be more efficient to use the former technique if the quadrature could be done sufficiently accurately with Q n3 /m points.
The integral form involving ρ can be written as
where K is an operator defined by
If one thinks of an integral simply as a very large summation, the analogy between the integral form and the original bilinear form is apparent. With a rapid way to evaluate K, there might be a faster way to evaluate the bilinear form. (Alternative ways to evaluate ρ, based on similar considerations, have already been discussed.)
Observe that K is simply the solution operator for Poisson's equation

OCR for page 43

Mathematical Challenges from Theoretical/Computational Chemistry
With this in mind, there are numerous ways to evaluate the action of K approximately, using discrete methods for approximating the solution of Poisson's equation. For example, the multigrid method can be used to do this in a very efficient way. The analogy with quadrature suggests using a grid with q points and solving Poisson's equation, with a multigrid method in O(q) work. Then the complete evaluation could be done in O(qmn) work. This potentially would be faster than direct evaluation of the integrals.
References
Bagheri, B., L.R. Scott, and S. Zhang, 1994, Implementing and using high-order finite element methods, Finite Elements in Analysis and Design 16:175–189.
Ding, Hong-Qiang, Naoki Karasawa, and William A. Goddard III, 1992, Atomic level simulations on a million particles: The cell multipole method for Coulomb and London nonbond interactions, J. Chem. Phys. 97:4309–4315.
Draghicescu, C.I., 1994, An efficient implementation of particle methods for the incompressible Euler equations, SIAM J. Numer. Anal. 31: 1090–1108.