Selected new developments in computational chemistry.

Molecular dynamics is a general technique for simulating the time-dependent properties of molecules and their environments. Quantum mechanics, as applied to molecules or clusters of molecules, provides a prescription for predicting properties exactly (in principle). It is reasonable to expect that both will have a profound effect on our understanding of environmental chemistry in the future. In this review, we consider several recent advances and applications in computational chemistry. ImagesFigure 1.


Introduction
Due largely to concurrent advances in theory, computer technology, numerical algorithm development, and the rapid influx of experimental information about the structure of molecules, computational chemistry stands poised to make significant contributions to our understanding of biological, and ultimately environmentally significant, processes. In this review we will consider some of the recent advances that seem most promising for the near future. We initially outline the recently developed particle mesh Ewald (PME) method. PME provides a mathematical approach for accurate evaluation of the Coulomb interactions in macromolecules. The development now makes possible the dynamical study of large environmentally important molecules by the well-established rules of classical molecules. We then discuss density functional theory (DFT). DFT is a relatively new methodology for investigating the next level of complexity-the electronic structure of Manuscript received 21 August 1995; manuscript accepted 17 October 1995. Address correspondence to Dr. Lee  Fax: (919) 541-7887. E-mail:pedersen@niehs.nih.gov Abbreviations used: PME, particle mesh Ewald; DFT, density functional theory; PDB, protein databank; HOMO, highest occupied molecular orbital. molecules. As with the PME method, significantly larger molecules of environmental interest can now be studied at the electronic level. Although both methods are global with respect to the expected impact on physical chemistry, it is in the area of environmental chemistry that profound applications may ultimately be found. Accurate Representation of Long-range Coulomb Forces The major source of information for the three-dimensional structures of biomolecules is that derived from X-ray crystallographic or from nuclear magnetic resonance (NMR) studies. The resulting information is normally tabulated as atomic Cartesian coordinates and isotropic thermal B factors. The latter are measures of the degree of thermal motion and are directly proportional to the mean square deviation from the average structure. For instance, a B factor of 5.0 for an atom indicates little motion, whereas a value of 75.0 indicates such large motion that the value of the coordinates can be taken to be quite uncertain. The structure files are ultimately (in most cases) deposited in databases such as the Brookhaven Protein Databank (PDB). These structures, universally displayed on two-dimensional graphics monitors, are static, unmoving representations of the molecules that do little to reveal the underlying motion of the atoms. The individual atomic motions can, however, be simulated by the technique of molecular dynamics (1). The underlying idea is that one assumes a potential energy function, or force field, and then solves Newton's equations of motion for the system. If there are N atoms in the system, then there are 3N equations, one for each Cartesian coordinate: Here FXi is the x force of the ith particle, vx,i is the corresponding velocity (vx,i = dxildt), and E(xi,..., ZN) is the potential energy function from which the force derives. The potential energy function is taken to approximate the more rigorous quantum mechanical solution. The function in its simplest form is a sum of quadratic terms in the displacement from classical equilibrium of bond stretches (A-B), angle bends (A-B-C), improper torsion angle contraints, a truncated Fourier expansion for the proper (bonded) torsion angles (A-B-C-D), Lennard-Jones terms for the nonbonded pairwise interactions (Alr'2-Blr6), which accounts for short-range repulsion of the outer electrons in a pair of interacting nonbonded atoms and for the long-range attractive dispersion interaction, and finally, a Coulomb's law term (qAqBlr) for the interactions of charges qA and qB. Its most important terms are those that affect space exclusion (atoms cannot be in the same place at the same time) and the Coulomb interactions that are long range (-l/r). The potential energy function is parameterized by establishing values for the constants in the function for the defined atom types. Parameterization derives from spectroscopic and thermodynamic measurements or from accurate quantum mechanical calculations, usually on small model compounds. In most current molecular dynamics programs, for instance, atoms are assigned static charges determined from accurate quantum mechanical calculations (2). The interaction between these charges is normally handled by Coulomb's Law, with interactions larger than a particular distance (8 A is popular) being truncated to make the calculations tractable. The truncation of the charged interactions, however, has been shown to lead to artifactual (3)(4)(5)(6)(7)(8)(9)(10)(11) Environmental Health Perspectives -Vol 04, Supplement a March 996 behavior for simulations of proteins and nucleic acids. In this section we describe recent progress in one approach for accurately accommodating long-range interactions between charges. Consider a crystal that is made of identical unit cells. Each unit cell holds a finite number of symmetry-related molecules. If we assign each atom in the molecule a charge, then we can determine a given atom's electrostatic interaction energy with all other charged atoms in the crystal by summing over Coulomb interactions in the unit cell and subsequently summing over all unit cells in the crystal. These summations are very slowly converging. In 1921, the Prussian physicist Paul Ewa!d (12) showed that the sums could be rewritten as a set of alternate summations, a direct space sum and an indirect space sum, that had much improved convergence properties (13). The indirect space sum, which has the natural co rdinates of the theory of X-ray crystallo raphy, was recently effectively solved by an implementation of fast Fourier transforms (5,7,14). The key idea was to approximate the charges as if placed on a regular grid using interpolation or spline formulae. The procedure, the PME, leads to a numerical algorithm with a time requirement of Nelog(N) rather than N2 for the formal Ewald procedure (the procedure can actually be optimized to behave as order N312) (Nis the number of atoms in the unit cell). For large macromolecules, the time speedup with PME is impressive, and thus we can now perform accurate calculations on systems that were previously unreachable. Although numerical, the error in the procedure is well defined in terms of the parameters of the method (14). The method also recently has been tested for periodic box images of single molecules surrounded by solvent and counterions (9) (Equation 1).
Applications of the Ewald sum to simulation problems in other laboratories include studies of salt solutions (15), planar interfaces (16), mobility of ions in solution (17), and DNA triple helices (18).
Theoretical details of the finite size effect that occur in Ewald-based simulations have been considered by Figueirido et al. (19).
An alternate approach for performing the long-range Coulomb summation, the fast multipole method (20)(21)(22), proceeds by first dividing a volume into approximate M subvolumes. An external potential for each subvolume is determined that then interacts directly with particles in nearby boxes but with a Taylor series expansion for particles in distant boxes. The overall time requirement can be shown to go as N, the number of charges, with the economy based on representing the interactions between volumes distantly spaced by truncated multipole expansions. A version of the fast multipole method has recently been applied to large macromolecular systems (23).
It is clear that our understanding of how to accurately simulate the motions of macromolecules, or clusters thereof, is increasing at a rapid pace. The application of these new techniques to problems in environmental chemistry should soon follow. It is now known (24), for instance, that exposure of animal cells to low levels of toxic compounds such as arsenicals, mercury derivatives, dimercaptans, peroxides, quinones, or diphenols leads to elevation of glutathione levels and induction of detoxification enzymes, e.g., glutathione S-transferases, quinone reductase, epoxide hydroxylase, and UDP-glucuronosyltransferase. The application of the newly emerging techniques should have a significant effect on our abilities to provide accurate representations of the timedependent interactions of drugs or toxic molecules with proteins such as those just described (once their three-dimensional structures are known), nucleic acids, or membranes.

Applications of Quantum Mechanics to Environmental Chemistry
Since its inception in 1926, quantum mechanics has held great promise for all areas of chemistry. Again, it has been the evolution of high-speed computing machines with massive storage devices and creative algorithm development that has led to the belief that this promise will be kept (25). The last few years have seen the implementation of the alternate density functional approach (DFT) (26) to the Schrodinger wave function method. In the latter, one (in principle) writes down a wave function for a chemical system that is composed of atomic orbitals of the individual atoms and then finds by application of the variation theorem (27) the lowest energy (best) wave function for computation of the properties (energy, geometry, vibrational frequencies) of the system. In the former (26), the energy is determined by the electron density (not a wave function); i.e., if one knows the electron density, the energy (and other measurables) can be calculated.
The DFT approach appears to be suited for large systems, since its computer time requirements scale with the number of particles N as N3, while the wave function approach scales as N4 or higher (Equation 2). Ewald (13) showed that this term could be written as  Figure 1D.
Here we display f(r), which identifies sites in a molecule that are susceptible to attack by soft electrophiles. This is in contrast to hard electrophiles in which the principal interaction proceeds under charge control (electrostatic interaction). The density functional computer program DMOL (Biosym Technologies, Inc., San Diego, CA) was used to determine the HOMO (highest occupied molecular orbital) and the electron density. This computer package uses numerical techniques to solve the Kohn-Sham orbital-density equation (28). The local exchange-correlation functional developed by von Barth and Hedin (29) is used in conjunction with a triple numerical plus double polarization basis set. The self-consistent field solutions are required to converge to 108 hartree and the geometries are optimized until the gradient is near 0.001 hartree/bohr. The topographical format for visualizing the possible reactivity sites is generated by displaying an isosurface of the electron density (an isodensity value of 0.002 in atomic units) that encloses the van der Waals volumes of the individual atoms. On this surface, we have mapped the positive contribution of the reactivity index onto the isodensity surface and have color coded the surface from red (most positive) to blue (zero). All pictures were generated with AVS (Advanced Visual Systems, Inc., Waltham, MD). Figure 1D shows that for azulene, the most likely position of attack by an electrophile is the carbon atom in positions 3 and 5 of the 5-membered ring. This is apparent from the topographical display of f(r), and, in this case, the display of the HOMO2, in agreement with the experiment.
There is no wave function in DFT; the electron density is every-  Carbon atoms are colored green and hydrogen atoms are white; B, the electrostatic potential of azulene mapped onto the isosurface of the electron density that just encloses the van der Waals volumes of the atoms, an isodensity value of 0.002. The resulting surface is colored from most negative (blue) to most positive (red). C, the HOM02 mapped onto the isosurface. The resulting surface is colored from zero (blue) to most positive (red). D, the Fukui function L(r) is mapped onto the isosurface. The resulting surface is colored from zero (blue) to most positive (red).  The structure, vibrational spectra, and force constants. The HOOCI isomer is found to be the lowest HF/MP2 energy form on the general potential energy surface using very high-level QM calculations; heats of formation, vibrational frequencies and intensities are also computed.
The stability of FOOCI with respect to the FO + CIO, FOO + Cl, and CIOO + F dissociation limits is examined by very high-level computations on FOOCI. A novel gridding technique is defined for computing the potential energy surface. The Polyrate code (37) is then used to compute the thermal dependence of the rate constants with good comparison to experiment. Classical transition state barriers were computed for the reaction and then the reaction rate estimated with tunneling corrections using an Eckart barrier; excellent agreement with experiment is observed. The transition state for tautomerization of formamidine is studied as a function of hydration, conventional basis/method and methods within DFT; the BH&H-LYP method is found to be roughly equivalent to the Mp2 calculations. Classical trajectory calculations on a potential energy surface parameterized to fit the quantum mechanical starting, transition state, and end orientations. The reaction is studied on two PESs reactive (one excited) to determine the fine structure effect; the overall rate constants are found to be determined by the ground-state path. Parameterized London-Eyring-Polanyi surfaces are used to describe the reactive surfaces; the reactants are embedded in cold inert gas clusters. The cluster plays a significant role in activation of the reactants on impact. Thermodynamic and activation barriers are computed for these potentially important atmospheric reactions. Activation barriers are computed for a number of hydrogen abstraction reactions. A simple method is proposed for correcting the energies and frequencies; transition state calculations give reasonable agreement with theory once tunneling corrections are made. Similar to previous work, except more extensive basis used; reasonable agreement found at T> 300 K for the calculated rates (transition state theory) and experimental rates. The reaction is studied with several different methods and with solvation corrections made with two different continuum methods. Hydrogen abstraction by CF30 from CH4, C2H6, and H20 is studied; the transition state barriers and thermodynamic properties for the gas phase reactions are found. The activation barrier for the formation of sulfuric acid is much lower with twc water molecules in the transition state complex than with one; in agreement with experiment.  In this section we will consider new developments in the application of quantum mechanics to problems of immediate or potential environmental interest in the areas of single molecule properties and/or reactivity. The goal of Table 1 is to give the flavor of modern theoretical work which has as its lofty aim the computation of accurate reaction rates from knowledge of only the system molecules. The literature cited is not exhaustive but illustrative.

Conclusion
The future appears bright for the application of sophisticated computational chemistry techniques such as those considered in this review to enviromentally related questions that have a chemical basis. The union of molecular mechanics/dynamics and quantum mechanics is being intensely studied in several laboratories (53)(54)(55)(56). These techniques make possible the introduction of time-dependent phenomena such as polarization and charge exchange and are thus more realistic in accommodating actual electronic behavior; new general tools that make available quantum mechanical dynamics should appear soon. The direct application to molecular problems of environmental interest will follow in short order.