Manuscript received 21 August 1995; manuscript accepted
17 October 1995.
Address correspondence to Dr. Lee G. Pedersen, MD A3-06,
NIEHS, P.O. Box 12233, Research Triangle Park, NC 27709. Telephone: (919)
541-4630. 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.
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 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:
Fx,i = mi dvx,i/dt
= -
E(xi,..., zN)/xi, i = 1 to 3N.
Here Fx,i is the x force of the ith particle, vx,i
is the corresponding velocity (vx,i = dxi/dt), 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 (A/r12-B/r6),
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 (qAqB/r)
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 (~1/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 Å is popular)
being truncated to make the calculations tractable. The truncation of the
charged interactions, however, has been shown to lead to artifactual (3-11) 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 Ewald (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 coordinates of the theory of X-ray crystallography, 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 N log(N) rather than N2
for the formal Ewald procedure (the procedure can actually be optimized
to behave as order N3/2) (N is 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).
Equation 1.
If a system of N-charged particles is repeated on a regular
cubic lattice l, the electrostatic energy is
where the charges are (qi, qj), rij
is the distance between charges, and the sum over the lattice is defined
by the integers lx, ly, lz where l = [lzL,
lyL, lzL] (l=0 is omitted for i=j). Ewald (13) showed
that this term could be written as
The second term in this equation is amenable to a fast Fourier
transform solution once the appropriate grid for evaluation is established.
The parameter ß is chosen to be large enough so that the minimum image
convention can be adopted.
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-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 time-dependent
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 Schrödinger 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).
Equation 2.
In conventional Schrödinger quantum mechanics, the
wave function can be found from the variational theorem by minimizing (with
respect to the parameters) the function (25)
where 
t is the trial wave function and
H is the electrostatic Hamiltonian. From this function, estimates can be
computed for any measurable property of the system. Most modern methodology
involves choosing the wave function as an antisymmetrized determinant of
molecular orbitals, which can hold a maximum of two electrons and which
themselves are linear expansions of atomic orbitals centered on the atoms
of the system. The solution of the Hartree-Fock equations gives the expansion
coefficients. The wave function can then be improved to account for the
dispersion energy (25) by configuration interaction or perturbation techniques.
In density functions theory (26), the orbital-density description
of the energy is given by:
There is no wave function in DFT; the electron density is everything.
The orbitals Oi that result from the solution are used to define
the electron density and do not have the same meaning as in Hartree-Fock
theory. Similarly, the orbital energy is not an approximation for the ionization
energy. Powerful interpretive concepts have arisen from the development
of DFT (26):
which is a chemical reativity index. Since p is not a smooth function
of N (when N is an integer), the derivative p /N is discontinuous. For integer N there exist three possible reactivity
indices. The left-hand derivative f-(r)
provides information about sites in a molecule susceptible to attack by
electrophiles. The right-hand derivative f+(r) is a reactivity
index for attack by a nucleophile, and the average of f+(r) and
f-(r) gives information about free radical
(homolytic) attack. Plots of the Fukui function about a molecule can show
positions susceptible to electrophilic or nucleophilic attack. The importance
and usefulness of the reactivity index is demonstrated in 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 10-8
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.
A
|
 |
B
|
C
|
D
|
Figure 1. Reactivity
index f-(r) for azulene. A, a ball-and-tube model of azulene.
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 HOMO2 mapped onto the isosurface.
The resulting surface is colored from zero (blue) to most positive (red).
D, the Fukui function f-(r) is mapped onto the isosurface. The
resulting surface is colored from zero (blue) to most positive (red).
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 fiavor 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-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.
REFERENCES
1. Brooks CL III. Methodological advances in molecular
dynamics simulations of biological systems. Curr Opin Struct Biol 5:211-215
(1995).
2. Wiberg KB, Rablen PR. Comparison of atomic charges derived
via different procedures. J Comp Chem 14:1504-1518 (1993).
3. Foley CK, Pedersen LG, Charison PS, Darden TA, Wittinghofer
A, Pai EF, Anderson MW. Simulation of the solution structure of the H-ras
p21-GTP complex. Biochemistry 31:4951-4959 (1992).
4. Hamaguchi N, Charifson P, Darden T, Xiao L, Padmanabhan
K, Tulinsky A, Hiskey R, Pedersen L. Molecular dynamics simulation of bovine
prothrombin fragment 1 in the presence of calcium ions. Biochemistry 31:8840-8848
(1992).
5. York DM, Darden TA, Pedersen LG. The effect of long-range
electrostatic interactions in simulations of macromolecular crystals: a
comparison of the Ewald and truncated list methods. J Chem Phys 99:8345-8348
(1993).
6. Frech M, Darden TA, Pedersen LG, Foley CK, Charifson
PS, Anderson, MW, Wittinghofer A. Role of glutamine-61 in the hydrolysis
of GTP by p21 H-ras. An experimental and theoretical study. Biochemistry
33:3237-3244 (1994).
7. York DM, Wlodawer A, Pedersen LG, Darden TA. Atomic-level
accuracy in simulations of large protein crystal. Proc Natl Acad Sci USA
91:8715-8718 (1994).
8. Auffinger P, Beveridge D. A simple test for evaluating
the truncation effects in simulations of systems involving charged groups.
Chem Phys Lett 234:413-415 (1995).
9. Cheatham TE III, Miller JL, Fox T, Darden TA, Kollman
PA. Molecular dynamics simulations on solvated biomolecular systems: the
particle mesh Ewald method leads to stable trajectories of DNA, RNA and
proteins. J Am Chem Soc 117:4193-4194 (1995).
10. Lee H, Darden TA, Pedersen LG. Molecular dynamics simulation
studies of a high resolution Z-DNA crystal. J Chem Phys 102:3830-3834 (1995).
11. York DM, Yang W, Lee H, Darden T, Pedersen LG. Toward
the accurate modeling of DNA: the importance of long-range electrostatics.
J Am Chem Soc 117:5001-5002 (1995).
12. Ewald PP. Die Berechnung optischer und elektrostatischer
Gitterpotentiale. Ann Phys (Leipzig) 64:253-287 (1921).
13. Allen MP, Tildesley DJ. Computer Simulation of Liquids.
Oxford:Clarendon Press, 1989;156-164.
14. Darden T, York D, Pedersen L. Particle mesh Ewald:
an N log(N) method for Ewald sums. J Chem Phys 98:10089-10092 (1993).
15. Smith PE, Marlow GE, Pettitt BM. Peptides in ionic
solutions: a simulation study of a bis(penicillamine) enkaphalin in sodium
acetate solution. J Am Chem Soc 115:7493-7498 (1993).
16. Hauptman J, Klein ML. An Ewald summation method for
planar surfaces and interfaces. Mol Phys 75:379-395 (1992).
17. Perera L, Essmann U, Berkowitz ML. Effect of the treatment
of long-range forces on the dynamics of ions in aqueous solutions. J Chem
Phys 102:450-456 (1995).
18. Weerasinghe S, Smith PE, Mohan V, Cheng Y-K, Pettitt
BM. Nanosecond dynamics and structure of a model DNA triple helix in saltwater
solution. J Am Chem Soc 117:2147-2158 (1995).
19. Figueirido F, Del Buono GS, Levy RM. On finite size
effects in computer simulations using the Ewald potential. J Chem Phys 103:6133-6143
(1995).
20. Greengard L. The Rapid Evaluation of Potential Fields
in Particle Systems. Cambridge, MA:MIT Press, 1987.
21. Schmidt KE, Lee MA. Implementing the fast multipole
method in three dimensions. J Stat Phys 63:1223-1235 (1991).
22. Board JA Jr, Causey JW, Leathrum JF Jr, Windemuth A,
Schulten K. Accelerated molecular dynamics simulation with the parallel
fast multipole algorithm. Chem Phys Lett 198:89-94 (1992).
23. Ding H-Q, Karasawa N, Goddard WA III. Atomic level
simulations on a million particles: the cell multipole method for Coulomb
and London nonbond interactions. J Chem Phys 97:4309-4315 (1992).
24. Prestera T, Talalay P. Electrophile and antioxidant
regulation of enzymes that detoxify carcinogens. Proc Natl Acad Sci USA
92: 8965-8969 (1995).
25. Levine IN. Quantum Chemistry. 4th ed. Englewood Cliffs,
NJ:Prentice Hall, 1991.
26. Parr RG, Yang W. Density-functional Theory of Atoms
and Molecules. New York:Oxford Press, 1989.
27. Levine IN. The variation method. In: Quantum Chemistry.
4th ed. Chap 8. Englewood Cliffs, NJ:Prentice Hall, 1991;189220.
28. Kohn W, Sham LJ. Self-consistent equations including
exchange and correlation effects. Phys Rev A 140:1133-1143 (1965).
29. von Barth U, Hedin L. A local exchange-correlation
potential for the spin polarized case. J Phys. C 5: 1629-1637 (1972).
30. Hehre WJ, Radom L, von Schleyer PR, Pople JA. Ab Initio
Molecular Orbital Theory. New York:John Wiley and Sons, 1985.
31. Frisch MJ, Frisch A, Foresman JB. Gaussian 94 User's
Reference. Pittsburgh, PA:Gaussian, 1994.
32. Vaida V, Simon JD. The photoreactivity of chlorine
dioxide. Science 268:1443-1448 (1995).
33. Lee JY, Ohyeon H, Lee SJ, Choi HS, Hansu S, Mhin BJ,
Kim KS. Ab initio study of s-trans-1,3-butadiene using various levels of
basis set and electron correlation: force constants and exponentially scaled
vibrational frequencies. J Phys Chem 99:1913-1918 (1995).
34. Francisco JS, Sander SP, Stanley SP, Lee TJ, Rendell
AP. Structures, relative stabilities and spectra of isomers of HClO2.
J Phys Chem 98:5644-5649 (1994).
35. Francisco JS. Theoretical characterization of FOOCl:
implications for the atmospheric chemistry between FOOx and ClOx
species. J Chem Phys 98:5650-5652 (1994).
36. Suzzi G, Orrú R, Clementi E, Lagana A, Crocchianti
S. Rate coefficients for the N+O2 reaction computed on an ab
initio potential energy surface. J Chem Phys 102:2825-2832 (1995).
37. Lu DH, Truong TN, Melissas VS, Lynch GC, Liu YP, Garrett
BC, Steckler R, Isaacson AD, Rai SN, Hancock GC, Lauderdale JG, Joseph T,
Truhlar DG. Polrate 4: a computer program for the calculation of chemical
reaction rates for polyatomics. Comput Phys Commun 71: 235-263 (1992).
38. Fu Y, Lewis-Bevan W, Tyrrell J. An ab initio investigation
of the reaction of trifluoromethane with hydroxyl radical. J Phys Chem 99:630-633
(1995).
39. Zhang Q, Bell R, Truong TN. Ab initio and density functional
theory studies of proton transfer reactions in multiple hydrogen bond systems.
J Phys Chem 99:592-599 (1995).
40. Abrash SA, Zehner RW, Mains GJ, Raff LM. Theoretical
studies of the thermal gas-phase decompostion of vinyl bromide on the ground
state potential-surface. J Phys Chem 99:2959-2977 (1995).
41. Schatz GC. Influence of atomic fine structure on bimolecular
rate constants: the Cl+HCl reaction. J Phys Chem 99:7522-7529 (1995).
42. Raz T, Levine RD. Four-center reactions: a computational
study of collisional activation, concerted bond switching and collisional
stabilization in impact heated clusters. J Phys Chem 99:7495-7506 (1995).
43. Franciso JS. An ab initio calculation of the energetics
for the FO+HCl
HOF+Cl reaction. Mol
Phys 82:831-833 (1994).
44. Rayez M-T, Rayez J-C, Sawerysyn J-P. Ab initio studies
of the reactions of chlorine atoms with fluoro- and chloro-substituted methanes.
J Phys Chem 98:11342-11352 (1994).
45. Melius CF. Thermochemical modeling. I. Application
to ignition and combustion of energetic materials. In: Chemistry and Physics
of Energetic Materials. New York:Kluwer Academic, 1992;21-49.
46 Dobbs KD, Dixon DA. Ab initio prediction of the activation
energies for the abstraction and exchange reactions of H with CH4 and SiH4.
J Phys Chem 98:5290-5297 (1994).
47. Davidson MM, Hillier IH, Hall RJ, Burton NA. Modeling
the reaction OH-+CO2HCO3-
in the gas phase and in aqueous solution: a combined density functional
continuum approach. Mol Phys 83:327-333 (1994).
48. Bock CW, Trachtman M, Niki H, and Mains GJ. Ab initio
study of the abstraction reactions of CF3O. J Phys Chem 98:7976-7980
(1994).
49. Morokuma K, Muguruma C. Ab initio molecular orbital
study of the mechanism of the gas phase reaction SO3+H2O:
importance of the second water molecule. J Am Chem Soc 116:10316-10317 (1994).
50. Bartolotti LJ, Edney EO. Investigation of the correlation
between the energy of the highest occupied molecular orbital (HOMO) and
the logarithm of the OH rate constant of hydrofluorocarbons and hydrofluoroethers.
Int J Chem Kinet 26:913-920 (1994).
51. Langenaeker W, De Proft F, Geerlings P. Development
of local hardness related reactiviy indices: their application in a study
of the SE at monosubstituted benzenes with in the HSAB context. J Phys Chem
99:6424-6431 (1995).
52. Arnett EM, Ludwig RT. On the relevance of the Parr-Pearson
principle of absolute hardness to organic chemistry. J Am Chem Soc 117:6627-6628
(1995).
53. Matsubara T, Maseras F, Koga N, Morokuma K. Application
of the newly "integrated MO+MM" method to the organometallic reaction:
Pt (PR3)2+H2. J Phys Chem (in press).
54. Hutter J, Tuckerman M, Parrinello M. Integrating the
Car-Parrinello equations. III: Techniques for ultrasoft pseudopotentials.
J Chem Phys 102:859-871(1995).
55. Brabec CJ, Maiti A, Bernholc J. Growth of carbon nanotubes:
a molecular dynamics study. Chem Phys Lett 236:150-155 (1995).
56. Hartsough D, Merz K Jr. Dynamic force field models:
molecular dynamics simulations of human carbonic anhydrase II using a quantum
molecular mechanical coupled potential. J Phys Chem 99:11266-11275 (1995).
[
Table
of Contents ] [
Citation
in PubMed ]
Last Update: September 2, 1998
[EHIS Home] [Comment
on article] [Tech Assistance]
[Search EHP] [Single
Copy Order Form]
