THE EMERGENCE OF COMPUTATIONAL CHEMISTRY
Computational chemistry has its roots in the early attempts by theoretical physicists, beginning in 1928, to solve the Schrödinger equation (see Box 2.1) using hand-cranked calculating machines. These calculations verified that solutions to the Schrödinger equation quantitatively reproduced experimentally observed features of simple systems such as the helium atom and the hydrogen molecule. Approximate solutions for larger systems and exact solutions to simple model problems allowed chemists and physicists to make qualitative explanations of spectra, structure, and reactivity of all types of matter.
During the Second World War, electronic computers were invented, and in the decade after the war these became available for general use by scientists. At the same time, physicists generally became more interested in nuclear structure and lost interest in the details of molecular structure and spectra. Hence, beginning in the mid-1950s, a new discipline was developed, primarily by chemists, in which serious attempts were made to obtain quantitative information about the behavior of molecules via numerical approximations to the solution of the Schrödinger equation, obtained by using a digital computer. The present success of this field has come largely from the enormous increase in speed, and decrease in cost, of computers, with significant improvements also attributable to many developments in algorithms and methodology. During the 1960s, three major developments in algorithms and methodology made quantum chemistry a useful tool: computationally feasible, accurate basis sets were developed; reasonably accurate approximate solutions to the electron correlation problem were demonstrated; and formulas for analytic derivatives of the energy with respect to nuclear position were derived. These developments were incorporated into several software packages that were made readily available to most chemists in the early 1970s, leading to an explosion in the literature of applications of computations to chemical problems. These programs are used to predict and explain the structure and reactivity of molecules and to complement the information obtained from many types of spectral measurements. Refinement of the program packages has, of course, continued, with emphasis on increased accuracy, increased size of molecules that can be studied, and adaptation to new computer hardware. The present methods have evolved from those that were used to study 1- and 2-atom systems in 1928 through those that were used to study 2- to 5-atom systems in 1970, to the present programs that produce useful quantitative results for molecules with up to 10 to 20 atoms. Much of the current research in new methods is aimed at developing methods that are feasible for even larger molecules.
A classic example of the power of the theoretical/computational approach is the work in the 1960s by W. Kolos and L. Wolniewicz. Explicitr12 calculations had been introduced for the hydrogen molecule in 1933 by James and Coolidge, and Kolos and Roothaan, working together in Mulliken's lab, improved these calculations in 1960. Subsequently, Kolos teamed up with Wolniewicz to author a sequence of papers of increasing accuracy. Their results diverged from the accepted (experimentally derived) dissociation energy of H2. When all known corrections were included, Kolos and Wolniewicz's best estimate of the discrepancy (in 1968) was 3.8 cm-1 Thus prodded, experimentalists reexamined the issue and in 1970 a new spectrum of better resolution and a new assignment of the vibrational quantum numbers of the upper electronic state were published. Both of these results were within experimental uncertainty of the best theoretical result.
While the emphasis of one aspect of computational chemistry has been on solving the many-body electronic structure problem, another group of chemists has focused on using the resulting potential energy surface for studying nuclear motion. This has led to a collection of programs for doing
BOX 2.1 The Schrödinger Wave Equation
The time evolution of a (nonrelativistic) quantum mechanical system is prescribed by the Schrödinger wave equation. For a particle with mass m and position r that is moving under the influence of a potential V(r), this equation reads:
where H is the linear Hermitian Hamiltonian operator:
Here is Planck's constant divided by 2π and has the very small value 1.05 x 10-27 erg-second characteristic of the submicroscopic quantum regime. The wavefunction ψ generally is complex; its amplitude squared ¦ψ¦2 provides the probability distribution for the position of the particle at time t. The linearity of the Schrödinger wave equation and the Hermiticity of H guarantees conservation of probability.
Time-independent solutions to the wave equation that have physical significance fall into two categories and require two distinct boundary conditions. The first (''bound states'') are square-integrable eigenfunctions of H; they are bounded, vanish at infinity, and provide a discrete spectrum of real energy eigenvalues for H. The second ("scattering states") are a continuum of solutions that lie in energy above the infinite-r limit of V(r) (if it exists), and the wavefunction remains nonzero (and oscillatory as a plane wave or sinusoidal function) in this asymptotic limit.
Quantum mechanical phenomena implied by the Schrödinger wave equation have a strongly counterintuitive flavor. The discreteness of bound-state energies contrasts starkly with the continuous energies available to a dynamical system in classical mechanics. The corresponding discreteness of the energy spectrum of transitions between pairs of quantum bound states has spawned the colloquially overused and misused phrase "quantum leap." In addition, the wavefunction can be nonvanishing in regions where V exceeds the total energy, leading to the well-known, but often mystifying, tunneling phenomenon; this permits particles to pass through potential energy barriers that in classical mechanics would provide absolute blockage. Less well known, but equally mystifying, is the reverse phenomenon originally pointed out by von Neumann and Wigner (1929): some potential functions V(r) have the capacity to trap particles in square-integrable eigenstates at energies above the absolute maximum of V, again in contradiction to classical mechanical behavior.
The Schrödinger wave equation adapts naturally to the description of many mutually interacting particles. The case most prominently considered in chemistry is that of two or more electrons, where V includes Coulombic interactions for electron-electron repulsions and electron-nucleus attractions. However, it is important to account for the fact that electrons in particular are endowed with intrinsic angular momentum, denoted by "spin," that is quantized to display only two allowed values, "up" and "down." In the strictly nonrelativistic regime (chemically, that which involves only atoms with low atomic numbers), these electron spins are invariants of motion and can be formally eliminated from the mathematical problem, provided that the remaining spatial wavefunction for the electrons satisfies appropriate symmetry conditions. In the case of two-electron systems (examples are the helium atom and the hydrogen molecule), the spatial wavefunction must either be symmetric under interchange of electron positions (if one "up" spin and one "down" spin are present) or antisymmetric (if both spins are "up" or both are "down"). Analogous, but more complicated, interchange conditions are applicable to the spatial wavefunction when more electrons are present.
von Neumann, J., and E. Wigner, 1929, Phys. Z. 30:465.
classical, semiclassical, and quantum calculations. Since 1980, use of these programs has become a routine tool for modeling molecules and gas-phase chemical reactions. These computations yield collision cross sections, both differential and integral, for elastic, inelastic, and reactive events. These approaches require, as with transition state theory, potential energy surface(s) obtained using quantum chemical methods of solution of the electronic Schrödinger equation. The Schrödinger equation for nuclear motion is solved subject to a scattering boundary condition, which takes the form of coupled differential, integro-differential, algebraic, or integral equation systems. The methods used to solve these coupled systems of equations are drawn from the applied mathematics literature as well as from algorithm improvements developed by computational chemists.
Meanwhile, simpler approximations have long been used by chemists to estimate the energy of molecules near their equilibrium geometry. In the molecular mechanics approach (see Box 2.2) the total energy of a chemical system is approximated by a sum of simple terms involving distances between atoms, bond angles, and dihedral angles. These terms involve estimated parameters that are assumed to have the same values as similar parameters obtained by data fitting for simpler molecules. (Chemists have long known that many structural and energetic features of molecules are nearly transferable between similar subfragments of molecules.) This representation of the energy has made possible the modeling of biological systems and rational drug design. It is also at the heart of the computational engine of many programs that produce three-dimensional computer graphics images of molecules. Molecular mechanics has become so prevalent that many chemists now equate it with computational chemistry. This approach has allowed the modeling of molecules with thousands of
BOX 2.2 Molecular Mechanics/Molecular Dynamics
Molecular mechanics and molecular dynamics refer to methods for computing certain molecular properties, particularly molecular structure and relative energy. They both typically use fairly simple potential energy functions that are derived from classical mechanics (e.g., a parabolic function to calculate the energy required to stretch or to compress a chemical bond). In addition, they both rely on parameters that are derived either from experiment (e.g., infrared spectroscopy and X-ray crystallography) or from quantum mechanics-based calculations (e.g., high-level ab initio molecular orbital calculations). A collection of potential energy functions and the associated parameters that are employed for molecular mechanics/molecular dynamics calculations is frequently referred to as a "force field"; thus, calculations that utilize the molecular mechanics or molecular dynamics approach are often referred to as empirical force field calculations.
The molecular mechanics method is generally employed to compute the relative energies of different geometries (conformations) of the same molecule that arise from rotations about chemical bonds as well as relative energies of intermolecular complexes. Often, energy minima are sought; thus, the molecular mechanics method is frequently coupled with optimization procedures. On the other hand, in the molecular dynamics method, Newton's equations of motion are solved by using the gradient of the above-mentioned potential energy function (force field) to compute the dynamic trajectory of a molecule or of an ensemble of molecules. Both the molecular mechanics and the molecular dynamics methods have found widespread use in the modeling of biomolecular systems, for which quantum mechanical calculations are simply not practical due to the overwhelming number of particles involved. Nonetheless, these methods are quite accurate for the estimation of certain molecular properties (i.e., those for which classical mechanics is appropriate), and they have been successfully employed to compute conformational energies (as described above), to estimate the binding affinity of small molecules bound to a macromolecular receptor, and as an adjunct for the refinement of structures derived from protein X-ray crystallography and protein nuclear magnetic resonance spectroscopy.
BOX 2.3 Chemist, Mathematician, or Physicist?
Douglas R. Hartree was an example of a truly interdisciplinary scientist. He invented the self-consistent field approach in 1927 to approximate the solution of the Schrödinger equation and did a number of calculations on atomic structure. His undergraduate studies had been interrupted by the First World War, during which he did ballistics computations. After completing his Ph.D., he held three chairs, first in Applied Mathematics at Manchester, then in Theoretical Physics at Manchester, then in Mathematical Physics at Cambridge. In 1934, he built the first differential analyzer in England. From 1939 until 1945, he used this in support of the war effort. After the war he was a consultant on the ENIAC computer and developed some of the first computer-based methods for solving partial differential equations. He taught numerical analysis at Cambridge University and wrote an insightful text on the subject, first published in 1952.
atoms. The practical disadvantage is that only structural types previously encountered in smaller molecules can be parameterized for larger molecules, so many parameters remain unknown. The conceptual disadvantage is that this is no longer a first principles theory and the connection to the Schrödinger equation is unclear. Hence, there can be no rigorous estimate of the potential errors in this approach and its success relies on chemical intuition for finding suitable molecules from which to develop the "transferable" parameters.
Another important thread in theoretical chemistry has been the study of many-particle systems such as liquids, solid materials, and biological macromolecules. The major framework for this study has been statistical mechanics—a subject with its formal roots in the nineteenth century. In the 1930s, the study by physical chemists of structure and thermodynamics accelerated with the advent of simple ideas about intramolecular and intermolecular forces. Equilibrium statistical mechanics has offered many questions of principle—for example, the question of the nature, and even definition, of phase transitions. These questions fostered a long-standing cross-fertilization between workers in both the mathematical and chemical communities (see Box 2.3). Similarly the study of phenomena away from equilibrium (e.g., the transport phenomena of hydrodynamics and the chemical rates) attracted the fundamental thinkers in statistical mechanics starting in the 1950s. Recently, corresponding deep questions of principle about disordered systems such as glasses have attracted workers from both communities.
Although a large part of statistical mechanics can be studied without computers, machine calculations for many-body simulation made an early impact in the 1950s and have grown to be the dominant mode of investigation. Monte Carlo methods, invented at the weapons laboratories by workers such as Fermi, Ulam, von Neumann, Metropolis, and Teller, were used immediately to address the many-body problems relevant to the thermodynamics of liquids. Such Monte Carlo approaches were adapted quickly to the study of polymers as well. The numerical solution of Newton's laws for many-particle models, so-called molecular dynamics (see Box 2.3), was also first carried out by theoretical chemists in the late 1950s and early 1960s. The application of molecular dynamics and Monte Carlo methods to proteins and other biomolecules in the 1970s has led to their widespread use throughout the theoretical and experimental chemical communities. Since significant advances in the efficiency of the algorithms used in molecular dynamics and Monte Carlo simulation are needed to address the forefront questions such as protein folding, a renewed contact of theoretical chemists with the numerical mathematics community has recently involved collaborative efforts of mathematicians, chemists, and physicists.
The advent of molecular quantum mechanics was followed by a very successful theory of chemical
reaction rates that modeled a reactive event as passage over a reaction barrier on a multidimensional potential energy surface representing the energy as a function of the internal coordinates of the reacting system. In its simplest form, the model corresponds to the system moving from reactants to transition states (the critical configuration), from which the system moves to reaction products. This conceptually simple model has remained the predominant approach for estimating rates of chemical reactions. Because of the multidimensionality of the reactive system, however, it is computationally difficult to implement rigorously. Over the years, efforts have focused on improving methods to estimate reaction barriers and properties of the reactants, and these have required better solutions of the electronic and nuclear transition states.
The roots of much of the mathematics now finding application to computational chemistry extend back at least to the eighteenth or nineteenth century, although, as illustrated in Chapters 3 and 4 of this report, the most up-to-date developments in the mathematical sciences can also be very natural tools. Group theory traces its origin to fundamental studies of geometries, but from it has come the theory of groups of motions, continuous groups, Lie groups, and Lie theory. The need to understand functions on the sphere and other surfaces led to the representation theory of groups and to various kinds of function theory. These theories grew up with the creation of quantum mechanics and fed, and were fed by, quantum mechanics. Much of operator theory and integral equations came from physics and engineering, as did the general theory of harmonic analysis. Numerical linear algebra and numerical analysis developed largely as tools for fluid mechanics and military applications, but their usefulness is vastly more widespread than that.
After World War II the mathematics community entered a period of intense development of its core, accelerating the growth of fields such as topology, number theory, algebraic geometry, and graph theory. Advances were largely motivated by questions generated by the internal structure of mathematics and not by contact with the outside world. In recent years, however, attention has once again turned outward, and the products of this intense period are now being applied widely in novel ways. The advent of modern computing capacity has enabled mathematicians to generate computational algorithms that yield answers—when combined with proper modeling techniques—to important practical problems. Success has been achieved in signal processing, sound and image compression, flow problems, and electromagnetic theory. Historically, mathematical scientists have worked more closely with engineers and physicists than with chemists, but recently many fields of mathematics such as numerical linear algebra, geometric topology, distance geometry, and symbolic computation have begun to play rules in chemical studies.
Before proceeding to accounts of past and potential contributions that mathematics can make to progress in chemistry, it should be emphasized that the challenge of interdisciplinary research is not one of scientific content alone, but also one of scientific process. Neither the chemist nor the mathematician is generally the optimal person to construct a mathematical model, as the model by its very nature lies at the interface between theory and observation. To build the model, an iterative process of refinement is required, in which mathematical considerations motivate approximations that need to be checked against reality, and in which key chemical insights necessarily force levels of mathematical complexity. It is exactly this need for iterative model construction that may motivate the collaboration of mathematicians and chemists, against the self-referential and conservative tendencies of each discipline. Focusing on this process of iterative model construction can help clarify the rules of the collaborators in interdisciplinary research, and by extension illustrate the goals for their respective disciplines as attempts are made to lower the hurdles to such collaborations. The model is both the interface between the disciplinary boundaries and the lingua franca between the cultures.