Watching the Molecules Move

On October 23-26, we were pleased to host a gathering of 185 scientists from around the world for the International Symposium on Computational Molecular Dynamics. These participants hailed from 14 countries and 27 states, and in less than 72 hours they bombarded the attendants with 25 plenary lectures and 107 poster papers. This article will summarize only a small fraction of the topics covered in these lectures and papers, but perhaps it will give a glimmer of the state of the field and the excitement we experienced seeing the new methods, the physical data generated, and many eyefuls of static and dynamic visualizations.

The talks covered a wide range of system sizes and time scales. Peter Lomdahl of Los Alamos National Laboratory discussed calculations on up to 250 million particles and in some calculations followed large-system motions up to nanosecond time intervals. The results were visualized by flying a computational airplane ("spy plane") over a barely visually resolvable bulk atomic landscape. In fascinating contrast, Craig Martens of the University of California at Irvine, presented calculations in which the motion of two iodine atoms in a rare gas crystal was followed in excruciating detail on the femtosecond time scale. The significance of every half oscillation of each of these two atoms was analyzed from several points of view. Other talks considered intermediate numbers of atoms (2 < N < 250,000,000) on intermediate time scales (fs < t < a few ns).

Materials Chemistry and Materials Physics

Uzi Landman of Georgia Tech presented a movie showing the properties of a layer of hexadecane, several molecules thick, on a gold surface interacting with a moving nickel tip. The role of molecular motion in the design of lubricants was elucidated in a most stimulating way. The wriggling hexadecane molecules groping for the receding Ni tip reminded one viewer of the sinners in Michelangelo's Last Judgment reaching for Jesus' soon-to-recede cloud. Whereas Michelangelo's picture seems to indicate that escape from purgatory was possible only with the help of angels, some of Uzis molecules escaped unassisted to the heavenly probe, for which they received a few cheers from the audience. Ted Davis and coworkers Dave Keffer, Dan Kroll, Alon McCormick, and Lianrui Zhang of the University of Minnesota along with Henry White of the University of Utah presented poster papers on molecular dynamics and Monte Carlo simulations of fluids confined in nanoporous and microporous materials. Enhanced layering and preferential dipole orientation were observed (showing features similar in some respects to those observed in Landman's simulation), and the dielectric constant of water in porous structures was found to be much smaller than that in bulk water.

Michael Gillan of the University of Keele (United Kingdom) discussed first principles calculations of the site dependence of catalytic oxygen exchange at oxide surfaces. These calculations revealed a strong site dependence for the energy to remove an O atom from a MgO crystal, e.g., 8.9 eV at a kink or 10.9 eV at a step, compared to 12.5 eV in the third layer of the bulk.

Michael Klein of the University of Pennsylvania presented a survey of his considerable computational progress in explaining conducting solutions of Li in ammonia, including the role of individual bipolarons, biplaron aggregates, and conduction by electron percolation through tubular networks. This is a beautiful example of computational molecular dynamics and visualization providing a clear theoretical framework for understanding in a field where years of analytic theory had yielded a considerably muddier picture.

Paulette Clancy of Cornell demonstrated the power of computational chemical dynamics on two problems of great practical importance. The first is the question of solute segregation at a growing solid/liquid interface. How slow must one keep the growth of the solid phase to be sure solute impurities are pushed ahead at the advancing interface rather than entrapped? To answer this, she simulated systems of 3000-8000 atoms for about a nanosecond. A second question she discussed is the mechanism of natural gas hydrate formation and dissolution in natural gas pipelines. The formation of hexagonal faces was identified as a kinetic bottleneck.

Several other presentations focused on properties of organic materials. James Bareman of Xerox Research Center of Canada discussed molecular dynamics simulations of the various crystalline forms of acridine, an example of the polymorph problem that has important consequences for the pharmaceutical, pigment, and electrophotographic industries. Martin Perry of the United States Naval Academy discussed molecular dynamics simulations of the friction that results when two reconstructed diamond surfaces are placed in sliding contact. Thomas Sewell of Los Alamos National Laboratory discussed Monte Carlo calculations of the properties of benzene crystals as functions of thermodynamic state parameters. Ilja Siepman of the University of Minnesota presented a successful calculation of the shear viscosity and vapor-liquid coexistence curve of n-decane. Mark Mondello of Exxon presented a calculation of the diffusion constants of n-decane and other large alkanes. Jutta Köhler of Wacker Chemie in Munich, Germany presented molecular dynamics simulations of the separation of chiral cyclodextrin complexes by enantioselective gas chromatography.

The study of glasses and the glass transition is an area of potentially great rewards because of the widespread importance of these subjects in industrial applications. Paul Madden of the University of Oxford (United Kingdom) discussed calculations of the fragility and workability of SiO_2 and how this is affected by the addition of Na_2O. Raymond Mountain of the National Institute of Standards and Technology discussed the glass transition in polymers by studying how molecules become localized and remain associated with a given group of neighbors.

Ilan Benjamin of the University of California at Santa Cruz and Woods Halley of the University of Minnesota presented molecular dynamics simulations of electron transfer at interfaces. Benjamin's lecture concerned a water/dichloroethane interface; Halley's study involved a metal/electrolyte interface. The study of interface phenomena is very challenging, but these researchers showed that simulations provided insights unobtainable in any other way.

Direct Dynamics

An exciting highlight of the symposium was the large number of papers on the emerging technique of direct dynamics, which has proven to be a powerful method both in many-body materials simulations and also for the study of few-body reactions. Trygve Helgaker of the University of Oslo presented an exciting overview of the subject, covering work from the 1970s to the present. He quoted the following definition: "Direct dynamics calculations [are] the calculation of rates or other dynamical observables directly from electronic structure information, without the intermediary of fitting the electronic energies in the form of a potential energy function" [Journal of Physical Chemistry, Vol. 95, page 4618 (1991)]. Trygve's talk climaxed in his discussion of his own trust region method for carrying out such calculations, with applications to unimolecular dissociation and proton transfer reactions.

François Gygi of Institut Romand de Recherche Numerique en Physique des Materiaux, Lausanne, Switzerland, presented an extension of the Car-Parrinello approach to direct dynamics in which generalized forces are used to adapt the coordinate system used for the electronic structure calculation such that it is finer where the action is and coarser where it is not.

Emily Carter of U.C.L.A. argued that for transition metals an orbital-based direct-inversion-in-the-iterative-subspace method, in which the electronic wave function is iterated to convergence at each dynamics step, is faster than Car-Parrinello Lagrangian procedures with the reversible RESPA algorithm. This debate may continue.

Jim Chelikowsky of the University of Minnesota argued for yet another alternative for direct dynamics calculations, namely the use of high-order finite differences to solve the electronic structure problem.

Quantum mechanical direct dynamics calculations showed up in the talk by Uzi Landmann, who presented path integral molecular dynamics calculations of H_3O^+ inversion, and in a poster paper by Wei-Ping Hu, Ivan Rossi, and Donald Truhlar of the University of Minnesota, in collaboration with Jose Corchado and Joaquin Espinosa-Garcia of the University of Extramadura in Spain. The latter presentation involved a multidimensional tunneling calculation of the OH + NH_3 -> H_2O + NH_2 reaction using dual-level techniques in which low-level electronic structure calculations are carried out at many points, and their accuracy is improved by interpolating corrections from a few points where high-level calculations are carried out. Both small-curvature and large-curvature semiclassical tunneling were included.

Paulette Clancys talk included a stimulating comparison of the advantages and costs of using semiempirical potential potentials vs. using tight-binding (i.e., extended Hückel) or ab initio direct dynamics for materials simulation. She emphasized the time-consuming aspect of finding accurate tight-binding parameters. A similar problem occurs in gas-phase dynamics if one tries to optimize parmeters for neglect-of-diatomic-differential-overlap (NDDO) electronic wave functions. Ivan Rossi and Donald Truhlar of the University of Minnesota presented a genetic algorithm approach for automating this process.

Thanh Truong of the University of Utah presented an application of direct dynamics to the reaction of H with H_2O, using variational transition state theory with small-curvature semiclassical tunneling. The energy along the minimum-energy path was treated by a dual-level technique, but frequencies were evaluated at a single level. These calculations lead to the conclusion that the widely used Walch-Dunning-Schatz-Elgersma potential energy surface allows too much tunneling for this prototype system. Other papers on direct dynamics included carbon cluster dynamics by Satoshi Itoh of the University of Illinois, Madhu Menon of the University of Kentucky, and Shep Smithline of Cray Research, CH_5^+ vibrational dynamics by Z.F. Liu and Y.T. Lee of the University of California at Berkeley, and mode-specific reactive decomposition of an ion-dipole complex in a poster paper by Gilles Peslherbe, Wei Chen, H.B. Schlegel, and Bill Hase of Wayne State University.

Electronic Structure

Advances and limitations in electronic structure techniques are often intimately connected to what is and is not doable in computational molecular dynamics. Steve Walch of NASA Ames Lab presented multi-reference configuration interaction (MRCI) calculations on various electronic problems with up to 10^8 configurations. The bond energy of N_2 was shown to be very slowly convergent with respect to basis set size in such large-CI calculations:

					Number of 
					bases
Basis set		D_e(kcal/mol)	functions

Cc-pVDZ	        	200		28
cc-pVTZ	        	217		60
cc-pVQZ	        	223		110
cc-pV5Z	        	225		182
Experiment		228

This provides a dramatic illustration of the need for dual-level approaches to dynamics. Walch illustrated the accuracy of ab initio calculations for small systems by calculations of the barrier height for unimolecular decomposition of the metastable HN_2 molecule. Here it is possible to extrapolate MRCI calculations to the basis set limit. Next minimum energy paths were calculated for the N+O_2 surface using both complete active space self-consistent field (CASSCF) and MRCI methods, and CASSCF gradient methods were used to map out potential energy surfaces for reactions involved in the formation of NO_x and soot in combustion of hydrocarbon fuels.

Jan Almlöf and coworkers of the University of Minnesota presented three technique-oriented electronic structure posters. One, with Angela Wilson, discussed new techniques for using localized orbitals in correlation-energy calculations. Another, with Yu Cheng Zheng, discussed an approach to nonlocal density functional calculations based on grid-free representations of the density. The third, with Matt Challacombe and Eric Schwegler, described techniques (that could be used in the fast multipole method or independently) for rapid assembly of the Coulomb matrix.

Four-body Problem

Exciting progress on the quantum mechanical four-body problem was reported by several workers. Daniel Neuhauser of U.C.L.A. and John Zhang of N.Y.U. presented accurate dynamics calculations of the full-dimensional OH + H_2 reaction. Evelyn Goldfield of Wayne State University presented a massively parallel study of a reduced-dimensionality model of the OH+CO reaction, a harder problem since it has two atoms heavier than hydrogen.

Bill Necoechea and Donald Truhlar of the University of Minnesota presented converged full-dimensional calculations of tunneling splittings and multi-quantum vibrational excitation energies of (HF)_2 that agreed excellently with earlier calculations by John Zhang, Zlatko Bacic, and their coworkers at N.Y.U. for two potential energy surfaces. The tunneling motion involves a concerted hydrogen bond switch: HF...HF -> FH...FH. The effect of monomer stretches on the tunneling motion is found to be small. Nadine Halberstadt of Université Paul Sabatier in France presented a quantum mechanical simulation of the spectrum of the very floppy He_2Cl_2 cluster which agreed well with an experiment by Ken Janda of the University of Pittsburgh.

Quantum Scattering and Dynamics

Claude Leforestier of the University of Paris South presented an interesting approach to quantum mechanical kinetics in which the recursive Lanczos scheme was used to compute the largest eigenvalues of the flux matrix for the reaction of H with O_2. The high-temperature rates were computed by a separable-rotation-like approximation. Leforestier emphasized the sensitivity to the assumptions about the transition state geometry used to compute effective moments of inertia.

A very recent development is the ability to carry out accurate quantum dynamics calculations for electronically nonadiabatic chemical reactions. First examples of such calculations were presented in posters by George Schatz of Northwestern University and by Steve Mielke, Tom Allison, Greg Tawa, and Donald Truhlar of the University of Minnesota and David Schwenke of NASA Ames Lab. Electronically nonadiabatic dynamics in solution were featured in a paper by Peter Rossky of the University of Texas, who discussed the ultrafast dynamics of hydrated electrons, including comparison with experimental work by Paul Barbara of the University of Minnesota. George Hagedorn of Virginia Polytechnic Institute and State University presented wave packet studies of adiabatic and nonadiabatic propagation through various types of avoided crossing regions, and Bill Lester of Lawrence Berkeley Laboratory presented reorientation cross sections for electronically excited H_2 colliding with He.

Of the quantal studies mentioned in the previous section and in this section so far, six used time-independent quantum mechanics and five used time-dependent quantum mechanics. Three papers presented new time-dependent quantum mechanical techniques, and two presented new techniques for time-independent computations. Daniel Neuhauser presented a time-dependent bound-state method in which short-time filters are used to get a basis for diagonalization. Zeki Kuruoglu of Bilkent University in Turkey presented a coupled-arrangement-channel wave packet approach to three-body collisions. Ronnie Kosloff of the Hebrew University in Jerusalem presented a stimulating analysis of how to get the asymptotically best attainable performance for quantum mechanical time-dependent treatments of very large problems.

Antonio Lagan of the University of Perugia in Italy presented a PVM-based parallel approach to the time-independent quantum scattering equations. B. "Ramu" Ramachandran of Louisianna Tech University, in a paper coauthored by Xudong Wu of Ohio State and Bob Wyatt of Texas, presented a new absorbing-potential, time-independent, variational method for chemical reactions.

Solvation

Quantum mechanical techniques for including solute polarization by the reaction field of a solvent are advancing rapidly and were represented at the symposium by several papers. Kenneth Merz of Penn State University and Anders Broo of Chalmers Technical University in Sweden presented hybrid quantum mechanical/molecular mechanical calculations of an S_N1 reaction and of absorption spectra, respectively. Ernest Mehler of CUNY presented a coulomb screening model for use in molecular dynamics calculations. Jayendran Rasaiah of the University of Maine presented a molecular dynamics simulation of the ionic mobility of alkali cations in water. A poster paper by Candee Chambers, Chris Cramer, and Donald Truhlar of the University of Minnesota and Edet Archibong, Ali Jabalameli, and Richard Sullivan of Jackson State University presented solvated studies of the transition states for internal rotation of dimethylthiourea.

David Ferguson of the University of Minnesota presented a new flexible water model and examples of its application to study the thermodynamic, dielectric, and dynamical properties of liquid water. The calculated solvation energies of water (-7.0 kcal) and neon (+2.7 kcal) both agree well with experiment (-6.3 and +2.7 kcal, respectively).

Jim Rustad of Pacific Northwest Labs and Woods Halley of the University of Minnesota presented molecular dynamics simulations of H_4SiO_4 and Fe(III)(OH)_6 deprotonation in aqueous solution as analogs of heterogeneous minerals in contaminated soils.

Susan Tucker of the University of California at Davis presented molecular dynamics simulations of the structure of a simple fluid at temperatures above its critical point. Such supercritical fluids are of technological importance because they exhibit novel--but poorly understood--solvating properties. Tucker's simulation results show that under certain thermodynamic conditions these fluids are not structurally homogeneous. Rather, they form alternating regions of high density ("clusters") and low density ("cavities"). Tucker showed that this inhomogeneous structure is an equilibrium property, even though each individual cluster persists only for a finite length of time. She also discussed a number of new statistical functions that she and coworker R. Ravi have developed to characterize these equilibrium structures.

Clusters

Jim Chelikowsky presented calculations of the photoemission spectrum of Si_n^- with n = 4 - 7 at 1500 K, at which temperature he obtained good agreement with experiment. Lower temperatures are harder to simulate because of strong vibrational-electronic coupling.

Julius Jellinek of Argonne National Laboratory presented a theoretical study of the fragmentation of Ni_n clusters with n = 12 - 14. The calculated preferred channel for fragmentation is ejection of a single atom, in good agreement with experiment.

Molecular Biology

There were several exciting presentations in the computational biochemistry area. Three of these involved DNA or dinucleotides. Tamar Schlick of N.Y.U. presented a Debye-Hückel model of the salt effect on DNA supercooling. Alexei Stuchebrukhov of the University of California at Davis calculated the electron transfer rates between two transition metal complexes bridged by a DNA molecule. Tom Darden of the National Institute of Environmental Health Science reported on nanosecond simulations of dinucleotides.

Terry Stouch of Bristol Meyers Squibb presented a very ambitious study of molecular permeation of model cell membranes, a subject that is basic to the question of the bioavailability of drugs. He found, for example, that methane diffuses four to six times faster through the center region of phospholipid bilayers than it does near head groups. Stouch also discussed a new project involving a simulation of transmembrane proteins.

Jan Hermans of the University of North Carolina discussed the free energy changes in conformational transitions of peptides and outlined the advantages of "soft" simulations vis a vis "hard" experiments for studying peptide structure. Among the advantages of the computational approach that he listed are availability of nonexperimentally accessible detail and the ability to study metastable states, poorly soluble molecules, and arbitrary pathways.

There were many talks on the molecular dynamics of proteins. Dave Pearlman of Vertex Pharmaceuticals asked the question "Can Free Energy Perturbations Be Made More Reliable?" and presented the answer: "Yes--in favorable cases." The answer was illustrated by case studies. He also explored the question of single vs. dual topologics for mutations and found that while single topology (as used in AMBER) is sometimes much better, dual toplogy (as used in CHARMM) is sometimes required.

Lawrence Ho of Yale University presented a hybrid quantum mechanical and molecular mechanical study of structural aspects of malate dehydrogenase redox enzyme catalysis. In particular he presented a map of the complete minimum energy reaction path. He implicated low-barrier hydrogen bonds in the protein transfer process.

A poster paper by John Kremer and Bill Gleason of the University of Minnesota presented a simulation of the conformational changes in acidic fibroblast growth factor that are induced by its binding to the drug sucrose octasulfate in aqueous solution.

Minoro Saito of the Protein Engineering Research Institute in Japan presented calculations of a large pK_a shift of an active site residue as compared to the monomer. The calculated shift of 2.1 units for Asp10 in RNaseHI is in good agreement with experiment (2.4). Solvent effects were discussed in terms of decomposition analysis of the free energy.

Eric Boczko of Carnegie Mellon University and Julian Tirado-Rives presented unfolding simulations--the former based on a 48-residue three-helix bundle and the latter on 110-residue barnase. Progress in understanding the structural factors and visualizing the pathways is dramatic.

Jean Durup of Université Paul Sabatier in France presented the "egg model." This is a powerful simulation protocol that he demonstrated by its application to a 36,999-atom calculation of the dissociation of a protein-protein complex. As the proteins separate, 370 confined waters move into the space between them. Durup analyzed the consequences of three possible dissociation mechanisms for opening up the space for these waters: hinge motion, shear motion, and unflattening of touching surfaces.

Monte Pettitt of the University of Houston discussed new techniques for analyzing the hydration waters near the surface of myoglobin surrounded by 4706 water molecules. An impressive statistical analysis of these waters was presented, showing how computer experiments help resolve experimental discrepancies in the placement of water around proteins. In addition, he showed one particularly interesting snapshot, in which his techniques allowed us to "see" the channel to the binding pocket open and "visualize" the water molecules entering the pocket. Although Monte emphasized that this particular picture was not necessarily statistically significant ("This is an anecdote--not statistical mechanics"), we all enjoyed watching the molecules move.

The Supercomputer Institute is grateful to Evelyn Goldfield, M. Katharine Holloway of Merck Research Labs, William L. Jorgensen of Yale University, Peter Rossky, and George C. Schatz for their unselfish service on the organizing committee. The local members of the organizing committee were Jan Almlöf and Donald G. Truhlar. The symposium was sponsored by the Computers in Chemistry and Physical Chemistry Divisions of the American Chemical Society and the Division of Computational Physics of the American Chemical Society, and we are grateful to Paul Barbara, University of Minnesota, Andrew DePristo, Iowa State University, Richard DeVore, University of South Carolina, Gary Doolen, Los Alamos National Laboratory, Karen Rappaport, Hoechst Celanese Corporation, and Michael Scanlon, American Physical Society, for their help with these cosponsorships.