Historical Perspective and Current Outlook for Molecular Dynamics As

Ind. Eng. Chem. Res. 2010, 49, 3059–3078

3059

Historical Perspective and Current Outlook for Molecular Dynamics As a Chemical Engineering Tool 1. Introduction th

On the occasion of the AIChE’s 100 anniversary, we review the status of molecular dynamics simulation. The development of molecular dynamics has been highly interdisciplinary, but the roles of chemical engineers have been significant. The role of chemical engineers is evolving as we move toward tools for routine applications like optimized force fields, coarse grained connections with the scale of finite elements, and physical property estimation. The trends with time lead to expectations for future developments. Obstacles to broader application and possible remedies are also mentioned. This commentary is an outgrowth of a presentation made at the Centennial session on molecular modeling at the 2008 Fall American Institute of Chemical Engineers conference. It explores the origins of MD, its role as a useful tool for chemical engineers, and the outlook for future roles of MD in chemical engineering research. It is not intended as a comprehensive review of the subject. Citations are selected to illustrate key points in the authors’ perspective on MD. From our perspective, MD is more than a technology of the future. It is the foundation of many present macroscopic models. In several cases, MD obviates the need for macroscopic models entirely. We expect that this trend will continue, first inferring qualitative behavior then achieving accurate inference directly from the force fields themselves. As we approach an age of molecular self-assembly and materials engineered at the atomistic level, the importance of this tool in everyday practice can only increase. 2. What is Molecular Dynamics?

These complex force fields may include many-body effects and electronic polarizability. In so-called ab initio molecular dynamics (AIMD), energetic interactions between atoms are computed using quantum mechanical methods, thereby eliminating the need for a classical force field altogether, albeit at a much greater computational expense. This trade-off between fundamental precision and computational feasibility permeates any consideration of MD. MD simulation should not be confused with Monte Carlo (MC) simulation (MC methods are reviewed by D. N. Theodorou elsewhere in this volume). Both methods provide atomistic details of fluids and solids at equilibrium, but MC is based on moving atoms according to rules that guarantee equal and opposite fluctuations around the equilibrium. These movements do not satisfy the laws of classical mechanics, and the distinction can be quite considerable. For example, equilibrating a polymer melt by MD can be complicated by the tangling of long chains and the incidentally slow separation of the entangled atoms. In MC, it is feasible to establish atomic movements that swap the identity of chains to which the atoms belong, instantly eliminating the tangle. Clearly, it would not be reasonable to estimate transport properties from such an MC method. This distinction indicates different roles for MD and MC. MC can be the best approach when equilibrium properties are the primary interest. MD is most appropriate when rigorous estimates are desired for the dynamical trajectories of all the atoms.1 Equation 1 shows a pairwise force field commonly used in the contemporary literature

∑ k (r - r ) + ∑ k (θ - θ ) + ∑ k [1 + cos(nχ - δ)] + ∑ k (ψ - ψ )

φtot )

Molecular dynamics (MD) is a computational technique in which atoms are treated as particles moving under the influence of classical mechanics. It is a deceptively simple method, yet it is one of the most powerful tools in engineering and science to understand the behavior of fluids and materials at the atomistic level. Historically, MD has played a key role in corroborating theories of the fluid state. In recent years, it has also become an instrument for directly evaluating physical properties of interest. Going forward, we expect this interplay between theory and simulation to continue, leveraging the benefits of each to expand the understanding of matter at ever increasing length and time scales. 2.1. Fundamentals of Molecular Dynamics. The underlying assumption of MD is that atoms can be treated as classical particles within the Born-Oppenheimer assumption. The Hamiltonian of the system depends only on the position and momenta of the collection of atoms. Energetic interactions between atoms are modeled using potential functions (“force fields”) of varying complexity. As discussed below, these force fields can be relatively simple, capturing only the repulsive interactions between atoms. On the other end of the spectrum, force fields can be quite complex, capturing intermolecular and intramolecular interactions between atoms with a high degree of fidelity. * To whom correspondence should be addressed. E-mail: jelliott@ uakron.edu. † University of Notre Dame. ‡ University of Akron.

2

b

2

o

θ

bonds

o

angles

χ

dihedrals

N-1

N

∑∑ i)1

j>i

improper 12

{ [( ) ( ) ] 4εij

σij rij

-

σij rij

2

o

ψ

6

+

qiqj rij

+

}

(1)

In this equation, the total potential energy, φtot, is represented by a series of terms that capture different types of energetic contributions. The first four terms account for intramolecular energy due to bond stretching, bond angle bending, rotation around bonds of four adjacent atoms (“dihedral angle” motion), and so-called “improper angle” bending typically associated with out of plane motion of atoms in rings. The parameters in these terms (kb, kθ, kχ, kψ) and the nominal values of the intramolecular coordinates (r0, θ0, χ, ψ) are typically fit to reproduce geometries and energy profiles obtained from ab initio calculations of the isolated molecule in question. The parameters n and δ in the dihedral expression account for periodicity and phase shift of the dihedral potential over the full range of rotation. The last two terms in eq 1 accounts for so-called “nonbonded” interactions between atoms in different molecules or atoms in the same molecule but separated by (usually) three or more bonds. In eq 1, a simple 12-6 Lennard-Jones function is used to represent van der Waals interactions, with σ and ε representing the atomic collision diameter and energy well-depth, respectively. The Lennard-Jones function is only one form that may be used. Partial charges q are placed on atomic centers to mimic

10.1021/ie901898k  2010 American Chemical Society Published on Web 03/01/2010

3060

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

electrostatic interactions, modeled here with Coulomb’s law. Depending on the system under study, several of the terms in eq 1 may be absent. For example, simple nonpolar molecules can usually be modeled effectively without resorting to partial charges. There is no need of intramolecular terms for atomic species like argon, and even small molecular species like methane can often be modeled accurately as spheres. As described below, even van der Waals-type interactions can be modeled in a simplified way by using repulsive or discontinuous potential functions. Given a force field for the potential energy, one can then write the Hamiltonian of a system of N atoms as H(rN, pN) ) φ(rN) + K(pN)

(2)

where it is assumed that the kinetic energy, K, depends only on the momenta of the atoms (p) and is thus separable from the potential energy, φ, which depends only on atomic positions r. For an isolated system, the Hamiltonian will be a constant. The positions and momenta of each atom i are governed by Hamilton’s equations of motion2 pi ∂H ) ) r˙i ∂pi mi

(3)

∂H ) -p˙i ∂ri

(4)

These equations of motion can also be written equivalently in the Newtonian formulation by differentiating the momentum and substituting it in eq 4. The result is ∂H ) -mrï ∂ri

(5)

Recognizing that -

∂H ∂φ )) Fi ∂ri ∂ri

(6)

the conventional Newtonian equations of motion are recovered Fi ) mirï

(7)

Equations 3 and 4 or eq 7 can be integrated forward in time using finite difference techniques, given a suitable initial condition r(0), p(0). These equations will conserve the Hamiltonian, and the resulting trajectory will be consistent with the microcanonical (constant number of molecules, volume, and energy) ensemble. Various finite difference techniques may be used, as discussed below. On the surface, therefore, MD appears to be a quite simple technique. A suitable force field is required to capture the energetics of the system of interest. A system of 6N first order differential equations (or 3N second order differential equations) must then be solved numerically via finite difference techniques to obtain a trajectory of positions and momenta as a function of time. Given such a trajectory, thermodynamic and dynamic properties are then computed using well-known statistical mechanical analysis methods. As with anything, however, the devil is in the details and a number of advances have been necessary to make MD a useful tool for chemical engineering analysis. 2.2. Perspective. The fascinating story of MD is one of multidisciplinary and international contributions, including contributions from chemistry and physics as well as chemical

engineering. We emphasize this global aspect by highlighting the affiliations of many of the primary contributors. The history of MD is also a story of expedience and efficiency in the face of computational limitations. Within the fundamental framework of eqs 1-7, many simplifications can be made to address the trade-off between computational feasibility and faithfulness to the details of the atomic interactions. Historically, extreme simplification was necessary in light of primitive computing resources. Naturally, greater detail was included as computational resources improved, shedding light on the significance of various contributions in particular applications. Yet even today, simplified potential models play important roles in analyzing the fundamentals of larger and slower systems. Furthermore, successful theories of the interaction dynamics can accelerate the computations. At some point, simplifications in the algorithm or potential model can lead to a disconnection between the integrated dynamics and the time scales of actual atomic movements. Simplified models of this variety are useful for studying mesoscale behavior, like self-assembly. Nevertheless, we draw a distinction between these models and MD, similar to the distinction between MD and Monte Carlo (MC) methods. The history of MD and the perpetual drive for computational efficiency are so interwoven with those of mesoscale and MC models that the distinction is occasionally blurred. The distinctive role of MD, however, is to provide reliable trajectories of atoms in conjunction with classical mechanics. We defer to other authors for perspectives on MC1 and mesoscale3 methods. This article is organized generally along historical lines, but with occasional digressions or extrapolations in time. It is often convenient to outline how a classical technique is being adapted for present and future applications while discussing that technique in its historical context. Moreover, two approaches are pursued roughly in parallel: MD with continuous or discontinuous potential models. The algorithms for the two methods are quite different. Clearly, the discontinuous potentials are simpler, but they can be made to resemble the continuous models when desired, while maintaining the capacity to be simplified when computational speed would be benefited. It is helpful to present both approaches to provide a general perspective. Altogether, one perceives the continuing progress of MD through global contributions deriving from familiarity with computational expertise, statistical mechanical theory, and discernment of the essential physics while maintaining faithfulness to the trajectories of the modeled systems. We conclude with our perspective on MD’s role in chemical engineering and obstacles to its more prevalent use as a chemical engineering tool. 3. Early Developments Because the ability to carry out meaningful MD simulations requires access to computational resources, it is not surprising that the pioneers in the field were those who first had access to computers. In the 1950s when MD was “invented”, this meant researchers at national laboratories in the U.S. It is difficult to say when the first MD simulation was conducted, but certainly one of the first instances of something like an MD simulation in which the positions and momenta of particles were integrated forward in time occurred around 1953 when Enrico Fermi, John Pasta, and Stanislaw Ulam (with the vital but seldom-appreciated programming support of Mary Tsingou) simulated a onedimensional chain with both linear and weak nonlinear interactions.4 The simulations, run on the MANIAC I computer at Los Alamos National Laboratory, were an “experimental” test of


the way in which equipartition of energy of various modes is approached in such systems. Because of the nonlinearity of the system, analytic methods could not be used, but the problem was the perfect test of the new computers at Los Alamos. It is interesting that Fermi and co-workers actually called their simulations “experimental” tests of theory. Since the work was clearly not theory and there was not yet a concept of computational research, Fermi and co-workers resorted to calling it “experiment”. This foreshadowed the way in which the entire computational field would straddle the traditional fields of theory and experimentation, and become a new third leg of science. Fermi and co-workers believed that given a small perturbation, their system would eventually thermalize, showing ergodic behavior such that all the modes would be equally excited after some time. Instead, the calculations showed that the system exhibited an unexpected quasi-periodic behavior. The work was published as a Los Alamos technical report in 19555 and touched off a new field dedicated to the computational study of complex nonlinear systems. The first instance of an MD simulation of a fluid was carried out by Berni Alder and Tom Wainwright and published in 1957.6 These calculations were conducted at the (then named) University of California Radiation Laboratory in Livermore, California, and contain most of the elements of a modern MD simulation, although obviously system sizes and simulation times were modest due to computational limitations. The system under study was the hard sphere fluid, and Alder and Wainwright found evidence of a solid-fluid phase transition that had not been observed with earlier Monte Carlo simulations. Systems of 32 and 96 particles were simulated in a rectangular box with periodic boundary conditions. The pressure was computed using the virial theorem as well as radial distribution functions. Alder would go on to use MD simulations in a series of 18 papers spanning 21 years and, in the process, answer many fundamental questions regarding fluid behavior. As it happens, Alder was no stranger to chemical engineering. He obtained his BS in chemistry from the University of California at Berkeley and his MS in chemical engineering at Berkeley. He later went to the California Institute of Technology where he worked with Kirkwood. Berkeley and CalTech are two schools that place chemistry and chemical engineering in the same academic unit, and so it is not a stretch to say that MD has its roots as much in chemical engineering as it does in any other discipline. It is interesting that Alder may have also conducted one of the first Monte Carlo simulations. In a 1997 interview, Alder mentioned that he used a Monte Carlo procedure to equilibrate a system of hard spheres in England during the summer of 1950, running on the Manchester Mark I computer.7 This computer was built by Ferranti Inc., in conjunction with the University of Manchester.8 His early work, along with that of Kirkwood, Lewinson, and Frankel, is in fact acknowledged in the classic Monte Carlo paper by Metropolis and co-workers.9 The subject of hard sphere simulations falls in the general category of discontinuous potential MD (DMD). This methodology continued to develop, largely through Alder’s efforts but also with contributions from Rapaport,10 Hall,11 and Chapela.12 We review that development briefly. Shortly after the development of DMD, however, MD with continuous potentials began to advance. Continuous potentials facilitate the inclusion of detailed aspects of the atomic interactions. Continuous potentials represent the most common form of MD in current practice and that methodology is described at greater length.

3061

3.1. Molecular Dynamics With Discontinuous Potentials (DMD). Alder and Wainwright’s simulations appeared four years after the first Monte Carlo (MC) simulation paper, and it is fair to say that MD has always lagged behind MC in treatment of equilibrium properties. The MC method is more efficient computationally and offers more opportunities for finding “shortcuts” to the equilibrium properties. In MD, the trajectories are sacrosanct so shortcuts are relatively few, placing a high value on effective techniques. The beauty of the DMD method is its simplicity. It essentially reduces all molecular interactions into discrete collision events. These events can be attractive or repulsive in nature. The repulsive events are equivalent to familiar “billiard ball” physics. When an attraction occurs between nonbonded atoms, we call that an association. When the attraction occurs between bonded atoms, we can call that a bounce or a vibrating covalent interaction. In all cases, the collision equations are basically the same. At that point, the key features of the algorithm focus on efficiently sorting the list of events and collecting ensemble averages that characterize the various properties. 3.1.1. Simplicity of DMD. In general, we leave the detailed descriptions of MD fundamentals to standard references like the texts of Allen and Tildesley13 or Frenkel and Smit.14 As an introduction to DMD, however, the method is so simple that it is useful to present the key equations. The interaction of two hard spheres suffices to describe the DMD algorithm for all molecules. Through a reference frame argument, we can assume that one of the spheres is stationary near the x-axis and the other is traveling along the x-axis on a collision path. Courses in introductory physics typically derive the post collision dynamics in terms of the speeds of the two particles and their direction angles.15 Unfortunately, the vast majority of these introductions stop short of determining all four postcollision quantities, applying only conservation of energy and x- and y-momentum. The fourth constraint for smooth hard spheres comes from noting that the direction of the stationary particle is given by the position on its surface where impact occurs. Then, the angle of reaction of the stationary sphere must be along the vector between the spheres’ centers. Applying the conservation of momentum and energy, we have: ∆ui -∆uj -2(bcij /σ2)rcij ) ) mj (mi + mj) mi

(8)

where mi is the mass of the ith particle, ∆ui is the change in velocity, σ is the hard sphere diameter, rijc is the vector between spheres’ centers at the time of collision, bijc ) rijc · uij. The other principal equation is the time until the next collision. tcij ) {-bij - (b2ij - uij2(rij2 - σij2))1/2}/uij2

(9)

This shows one way that DMD differs significantly from continuous molecular dynamics (CMD). CMD may use 1-3 time step multiples. For example, the time step for intramolecular forces might be 1 fs (10-15 s), while the general time step is 5 fs. All atoms would be moved each 1 fs, and the intramolecular forces (which number relatively few but vary stiffly) would be recomputed each 1 fs, but the intermolecular forces (which may include many atoms) would be recomputed less often. DMD is the ultimate multiple time step algorithm. The time step is as long as the time until the next collision. This distinction might be significant in certain specialized applications. For example, a study of methyl radicals in a forest of growing carbon nanotubes might benefit from this approach.

3062


The simplicity of these equations and their familiarity from introductory physics courses can be beneficial when introducing molecular simulation, and thermodynamics in general. Students can compute the dynamics of the first collision themselves to demonstrate their comprehension, then let the computer perform the next million collisions. This is the mindset behind the DMD module at etomica.org.16 Students are led through a series of computer “experiments” that clarify molecular interactions, the nature of temperature, pressure, energy, corresponding states, and a number of other thermodynamic principles. The simulations are visual and interactive, appealing to a number of learning styles. 3.1.2. Alder’s Studies of Equilibrium and Transport Properties. From the beginning, Alder led the way in demonstrating the capabilities of MD simulation. In articles numbered I-XVIII, he examined the phase diagrams and transport properties of 2D and 3D square well and hard spheres, including applications to large enough systems to study hydrodynamics. These articles demonstrated the feasibility of using equilibrium MD (EMD) to characterize transport coefficients and the challenges of EMD for a property like viscosity. Alder also explored the transitions from fluid to solid and between solid phases as well as the phase equilibria of binary liquid mixtures. Alder limited his interests to spheres until relatively recently. Alder’s work exemplifies the power of MD to explore the phenomenology underlying fundamental behavior like phase transitions and hydrodynamics. His recent work also exemplifies the way that MD can contribute to efforts that span length and time scales. MD cannot replace continuum theories like the Navier-Stokes equations, but it can illuminate their deficiencies and point the way for making corrections. Using MD for this kind of leverage is the key to meaningful applications of MD in the future. 3.1.3. Extension to Polyatomic Molecules. In 1978, D. C. Rapaport from the Physics Department of Bar-Ilan University in Israel carried out the first MD simulation of a “polymer” in which freely jointed hard spheres were connected in chains as long as 50 units.10 The end-to-end distance and radius of gyration were computed for this model and compared with previous Monte Carlo simulations. Rapaport also refined and extended the DMD algorithm for hard spheres.17,18 Rapaport’s hard sphere algorithm takes advantage of a binary tree sorting method and extensive use of linked lists. This makes the computation time proportional to Nn log2 N where Nn is the number of spheres in neighboring cells and N is the number of spheres in the simulation. The basic idea is that spheres cannot collide unless they are in neighboring cells, as long as the cell dimensions are larger than the range of the potential. Otherwise, a cell crossing event must occur before the spheres can be close enough to collide. This approach is especially advantageous when the range of the simulated potential can be limited to the repulsive diameter with no attractive interactions. Attractive dispersion forces, as in a Lennard-Jones or square well potential, are longer range, increasing the number of neighbors in the range and slowing the calculations by an order of magnitude. Coulombic forces are even longer range and tend to slow the simulations by an additional order of magnitude. Since the attractive potentials vary relatively slowly, it is often feasible to delay evaluation of their forces with a multiple time step approach. Nevertheless, explicitly including attractive forces in the simulation usually slows the simulation by at least 1 order of magnitude. This is the motivation for seeking combinations of theory and simulation that permit the impact of attractive forces to be inferred quantitatively without simulating them

explicitly. Rapaport realized this and focused his efforts on the hard sphere algorithm. Rapaport’s key innovation for simulating polyatomic molecules was to realize that spheres can be tethered with potential wells that allow the bonds to vibrate but not break. The spheres inside their bonded wells travel like unbonded spheres. Therefore, the hard sphere algorithm can be applied directly. This approach was adopted by Carol Hall and co-workers in the Chemical Engineering Department at North Carolina State University to simulate thermodynamic and transport properties of tangent hard sphere chains, including simulations beyond the entanglement threshold of polymers.11,19,20 They identified fundamental trends in the diffusivity, viscosity, and thermal conductivity. Similarly, Chapela and co-workers in Mexico adopted the approach to study molecular fluids with more realistic descriptions of the bond lengths and multistep potentials to describe the attractive dispersion forces.12,21 Elliott and co-workers developed an approach similar to Chapela’s, but they demonstrated quantitative agreement with simulations when thermodynamic perturbation theory (TPT) was applied in lieu of explicitly simulating the attractive dispersion forces.22,23 Treating the attractive forces as perturbations enables refinement of their characterization without need of resimulating every trial potential. The discrete reference potential in DMD facilitates the application of perturbation theory because a soft reference potential like a shifted Lennard-Jones potential would require description in terms of both temperature and density. Nevertheless, discretized potentials are capable of reproducing the macroscopic results of a continuous force field potential (such as the transferable potentials for phase equilibria of “TraPPE” potential24) while maintaining the efficiency of DMD/ TPT, thereby providing the best of both worlds.21,25 This combination of DMD/TPT is referred to as SPEADMD, for step potential equilibria and discontinuous molecular dynamics. SPEADMD provides one example of the kind of leverage that can be achieved by combining theory and simulation. Applications of this approach are illustrated later in the context of physical property predictions. The ability of DMD to give accurate dynamics has been demonstrated in at least two ways. Liu and Elliott demonstrated the manner in which the lifetime of a hydrogen bond could be correlated with spectroscopic measurements of bonding decay dynamics.26 More recently, Gerek and Elliott have used DMD as the basis for correlating self-diffusivity measurements.27 Selfdiffusivity represents the simplest transport property that can be treated by MD and paves the way for extensions to thermal conductivity and viscosity. Naturally, details of the dynamics are correlated with details of the potential models and much work remains to be done to develop potential models that provide the greatest precision for the broadest range of properties. The challenges of developing potential models to characterize physical properties are discussed in section 4.4. 3.2. Continuous Potentials. Simulations involving continuous potentials followed closely on the heels of the hard sphere simulation work. The first instance of a continuous potential being used in an MD simulation of a fluid appears to have been carried out by Gibson, Goland, Milgram, and Vineyard at Brookhaven National Laboratory in 1960.28 This work utilized a Born-Mayer repulsive potential to simulate radiation damage in copper crystals. Systems of 500 atoms took about 1 min per time step to simulate on an IBM 704 computer. Aside from being the first paper to use a continuous potential, this amazingly comprehensive paper was also one of the earliest to use


animations, in which the positions of atoms were displayed as dots on cathode ray tubes. Another influential paper utilizing continuous potentialssthis time for the Lennard-Jones potentialswas carried out in 1964 by Aneesur Rahman at Argonne National Laboratory.29 This paper studied the behavior of liquid argon and remains a classic in the MD literature. In fact, one can easily find contemporary simulation studies that hew close to the path blazed by Rahman. In this paper, 864 Lennard-Jones atoms were simulated using periodic boundary conditions. Rahman showed how to compute most of the properties that are of interest in an MD simulation study, including radial distribution functions, self-diffusivities (via velocity autocorrelation functions and mean square displacements), van Hove-like correlation functions, and nonGaussian dynamical behavior parameters. Another MD pioneer whose name is associated with an integration scheme as well as a method for speeding up potential evaluations is Loup Verlet. In 1967 while working at Yeshiva University,30 Verlet introduced the concept of a “neighbor list” (now called a Verlet neighbor list) in which one keeps a list in memory of all the particles within a distance equal to the potential cutoff distance plus a “skin” thickness of each atom in the system. Rapaport’s linked list method is in the same spirit as a neighbor list. For simulations with a large number of atoms, the number of “neighbors” of any given atom is a small fraction of the total number of atoms, and only this list must be searched when computing the interaction energy for the central atom. As a result, the expensive loops involved in evaluating energies and forces are reduced significantly by only looping over those atoms that have a chance of being within the cutoff distance. The neighbor list is updated periodically to reflect the relative motion between atoms. This was one of the first of many subsequent algorithmic advances that have enabled MD to be applied to the kind of large, complex systems of interest to chemical engineers. In the same classic paper, Verlet also implemented an integration scheme that has come to be known as the “Verlet method”. Atomic positions are advanced at each time increment δt using the following equation ri(t + δt) ) -ri(t - δt) + 2ri(t) + ai(t)δt2

(10)

where ai(t) is the acceleration of atom i at time t. The method does not readily lend itself to the determination of velocities at a given time, so Verlet used a central difference method to estimate the velocities as Vi(t) ≈

ri(t + δt) - ri(t - δt) 2δt

(11)

The more commonly used and velocity-explicit “velocity Verlet” algorithm was introduced several years later by Swope, Andersen, Berens, and Wilson31 and takes the form 1 ri(t + δt) ) ri(t) + Vi(t)δt + ai(t)δt2 2 Vi(t + δt) ) Vi(t) +

ai(t) + ai(t + δt) δt 2

(12)

(13)

It is important to realize that despite Verlet’s pioneering work with MD, he was not a one-method scientist. In 1969, he worked with J.-P. Hansen at the Laboratoire de Physique Theórique et Hautes Energies in Orsay, France, to carry out an extensive series of Monte Carlo simulations in which the phase diagram of the Lennard-Jones system was mapped out.32 At the time (and still today), Monte Carlo was the superior technique for

3063

tracing coexistence curves, and so this was the technique he used. This is a valuable lesson that is sometimes lost on current practitioners who use only one method to compute all properties. Sometimes MD is not the most appropriate method. 3.2.1. Methods Improve. The first MD algorithms all sampled the microcanonical ensemble, in which the number of atoms, volume, and total energy are fixed. This is the natural ensemble of MD and is sufficient for many purposes given the equivalence of all ensembles in the thermodynamic limit. The microcanonical ensemble can be unwieldy on a practical level, however, especially when trying to simulate specific state points or make direct comparison with experimental systems where intensive variables like temperature and pressure are fixed. Monte Carlo simulations were far more convenient to use when simulations at predetermined state points were desired because different ensembles could be simulated easily by simply changing the probability density expression in the acceptance rules. Doing this for MD in a rigorous manner required development of new equations of motion that enabled the appropriate ensemble to be simulated. One of the earliest constant temperature MD simulations was conducted by Leslie Woodcock at the University of London in 1970.33 He used simple velocity rescaling to achieve constant set point temperatures for molten alkali chloride systems. A limitation of simple velocity rescaling, however, is that it does not properly sample the canonical ensemble probability distribution. Ten years later, Hans Andersen of Stanford University formulated an elegant means for “extending” the conventional Hamiltonian equations of motion (see eqs 3 and 4) in such a way so as to rigorously sample the canonical, isoenthalpic-isobaric, and isothermal-isobaric ensembles.34 Interestingly, it was chemical engineer J. W. Haile and physicist H. W. Graben who actually implemented the chemist Andersen’s algorithms and demonstrated the utility of these extended system methods for simulating in ensembles other than the microcanonical.35 This development provides a great example of the way in which molecular dynamics simulation had blossomed into a vibrant multidisciplinary field in which innovation was being driven by the diverse needs of different intellectual traditions. Haile was one of a long line of chemical engineering academics in this field trained under the direction of Prof. Keith Gubbins. Other methods for maintaining constant pressure and temperature were developed around the same time, including the Nose´-Hoover method in which external thermostat and barostat variables were introduced to maintain constant temperature and pressure.36,37 In this approach, a Lagrangian (rather than Hamiltonian) formulation to the equations of motion was used, with an extra variable having a fictitious “mass” introduced to control the temperature and pressure. This method also guarantees proper sampling of the canonical ensemble and gives the dynamics correctly. Several other methodological advances were developed during this time that enabled MD to become a practical tool for chemical engineering calculations. One particularly important one has to do with the way in which interactions are calculated between charges or dipoles. To enable polar or charged species to be modeled accurately, point charges typically are added to atomic centers (see eq 1). These electrostatic interactions are much longer ranged than the short-range van der Waals forces, and so, simple cutoffs do not work for these forces. In fact, any spatial interaction that decays no faster than r-3 (in threedimensional systems) is long ranged and simple cutoffs require huge box sizes. Practically, spherical cutoffs for charged interactions are an especially bad idea, since at any instant the

3064


charges within the cutoff may not cancel, and so unbalanced charging creates an artificial force. This problem was solved by de Leeuw, Perram, and Smith38 by exploiting old ideas proposed originally by Madelung and Ewald. The method is now generically referred to as the “Ewald sum” method. The basic idea is to recast the nonconvergent Coulombic sum into a combination of a real space and Fourier space sum, each of which are convergent. The method is applicable to point charges as well as dipolar systems.13 Many modifications have been proposed to the basic Ewald sum to speed the calculation, because it is a very time-consuming portion of any calculation. The most significant advances have involved using various “mesh” techniques that enable parallelization of the charge-charge calculations.39 In addition to these methodological advances, several other developments were needed to bring MD up to the point where it could be used to simulate realistic models of materials of interest to chemical engineers. A few of the major advances are summarized below. Rahman and Frank Stillinger conducted the earliest simulations of liquid water in 1971.40 These authors properly recognized the inherent difficulties of simulating water (difficulties that still confound modelers to this day), but nevertheless were able to compute a range of static and dynamic properties that provided insight into the behavior of liquid water. They found evidence for hydrogen bonding networks in their system and observed that diffusion took place through continuous cooperative motion between neighbors and not via a “hopping” mechanism. Shortly after in 1974, Rykaert and Bellemans simulated liquid butane using a realistic “united atom” model.41 They computed self-diffusion coefficients and saw very different dynamical behavior than what had been observed with simple spherical liquids, giving evidence that conformational motion of flexible molecules plays a role in their dynamical motion. In the late 1970s, two methods were developed that contributed to the efficient simulation of polyatomic molecules. Murad, Evans, Gubbins, Street, and Tildesley simulated fluid methane using a 5-site exponential-6 potential.42 In this work, quaternions were used to handle the rigid methane molecules, which enabled efficient simulations to take place. Evans43 and Evans and Murad44 had described this method a couple years earlier. Using a different approach, Ryckaert, Ciccotti, and Berendsen developed the “SHAKE” algorithm for constraining bond lengths.45 In both cases, the fast vibrational modes of stiff atomic bonds could be “ignored”, thereby enabling longer time steps to be taken. This in turn resulted in a great increase in the efficiency of MD simulations for multiatom molecular species. These and other algorithmic developments were coupled with ever increasing advances in computing speed and memory, such that by the mid-1980s large-scale simulations of realistic materials were possible. One particular advance of relevance to chemical engineering is the development of free energy perturbation methods for computing solubility and phase behavior, especially in aqueous systems. In 1982, Postma, Berendsen, and Haak published a paper in which isothermalisobaric MD was used to compute the free energy of cavity formation in liquid water.46 This work examined the hydrophobic effect of water and also served as a test of scaled particle theory. This was an early example of the dual role MD has continued to play to this day; it serves as both a probe of “real” systems and a test of theory. In 1987, Singh, Brown, Bash, and Kollman used free energy perturbation within an MD framework to study free energy differences between amino acids in aqueous

systems.47 Together, these two papers were among the first examples of the use of MD to compute free energies for systems of great practical relevance. 3.2.2. More Complex Systems. Model systems continued to grow in size and complexity, including mixtures and more detailed treatment of polar interactions. One reflection of this growing complexity was the approach toward polymeric materials. About ten years after Rapaport’s simulation of freely jointed hard sphere chains, Rigby and Roe of the Department of Materials Science and Engineering at the University of Cincinnati modeled polymers of the same chain length as Rapaport, but this time with potentials that accounted for bond stretching, angle bending, dihedral angle rotation, and site-site interactions through a Lennard-Jones potential.48 Their study uncovered a wealth of microscopic details, including evidence of a liquidto-glass transition. Interestingly, these authors specifically noted that they exploited both the vector processing capabilities of the then state-of-the-art Cray X-MP computer and the Verlet neighbor list technique. By this time, the twin developments taking place in algorithms and computer speed were opening up the field of molecular dynamics to materials scientists and chemical engineers, who were pushing the technique into unexplored territories. Two examples from this time period show how MD methods were penetrating the research community. In 1990, Kurt Kremer from the Physics Department at the Universita¨t Mainz and Gary Grest of Exxon Research and Engineering used MD to simulate entangled polymer melts.49 In 1991, Kevin Mansfield and Doros Theodorou of Berkeley’s Chemical Engineering Department studied glassy polymer surfaces with MD.50 No longer content with simple toy models, researchers in engineering fields and even industrial research laboratories began developing more realistic (if still simplified, due to computational limitations) models of materials and subjecting these models to in-depth MD simulations. As discussed below, these efforts have continued and have been driven by the applications areas, so that now very large-scale materials MD simulations are possible. 3.3. Software Developments. Most of the software used to conduct the research described above was written “in-house.” The papers from the earliest days reflect the fact that the MD community was relatively small and close knit. In those preinternet days, algorithms and methods were implemented on the basis of descriptions in publications while “trade secrets” were discussed at workshops and conferences. Organizations such as the UK’s Collaborative Computational Project 5 (CCP5) and Europe’s Centre Europeén de Calcul Atomique et Molećulaire (CECAM) organization greatly facilitated the exchange of information that helped lay the foundations of the entire field. Several excellent monographs on molecular dynamics and molecular simulation in general came out in the late 1980s and early 1990s, with perhaps the best-known and most influential being Computer Simulation of Liquids by Mike Allen and Dominic Tildesley.13 Published in 1987, this book contained an almost encyclopedic collection of simulation methods used up until that time. It explained how algorithms worked, who made the developments, and how results should be properly analyzed. Perhaps most importantly, it came with 37 example routines written in Fortran 77 that showed exactly how to implement the methods described in the book. It is probably not an exaggeration that these routines have been the place where the vast majority of molecular modeling practitioners first learned how to write simulation code. Whether it is a Brownian dynamics code for the Lennard-Jones fluid (program F.33) or a


routine for treating periodic boundary conditions (F.1), most of the basics can be found in these routines. There are probably hundreds of MD codes used and developed in research groups all over the world. A noticeable shift has occurred, however, toward research groups becoming users of a few well-established MD codes instead of deVelopers of their own local codes. This is perhaps inevitable and follows the trend of the electronic structure community. As new methods were developed, MD codes that incorporated these methods became faster and more complex. Writing a new code from scratch that incorporated all of these new developmentss especially efficient parallelizationsbecame a daunting task. As a result, the drive to abandon locally developed serial codes in favor of highly optimized parallel codes developed elsewhere became great. A large fractionsperhaps the majoritysof published MD work now uses one of a handful of codes such as LAMMPS,51 NAMD,52 DL_POLY,53 Tinker,54 CHARMM,55 Amber,56 or Gromacs.57 Each of these codes has an active user and development community who serve to add features, fix bugs, and improve the software. There are also commercially available MD codes, most notably those marketed by Accelrys58 and HyperChem.59 These tend to be easier to use and have better graphical front ends than the former “academic” codes, but their performance and functionality has tended to lag well behind those developed as open source codes. When commercial codes were first being developed, consortia of major corporations were formed to bring methods and software developed in academia into use by industry. Enthusiasm for such consortia has waned, and it is unclear whether this type of development model can be sustained.60 On the other hand, many of the open source MD codes have been developed and maintained through the generosity and expertise of a disparate user community. The open source or free software model has been successful in many different areas, and it appears to be a powerful driving force for molecular modeling software as well. The prospects of these open source MD codes hinges on the ability of future developers to have time and availability to continue their support efforts. It is not clear that the objectives of companies wanting to use molecular simulation and the network of open source developers and maintainers are always in line with one another. Mechanisms need to be found to support and expand molecular modeling software in the face of declining fundamental research investments by corporations. Funding agencies could help ensure these developments are made by supporting researchers who propose to add functionality to open source codes and by encouraging the support and maintenance of the main open source programs. Companies who see value in molecular simulation need to find ways to sponsor development of these tools as well. In general, the move toward MD users being less involved in code development and more involved in using codes has been beneficial, because it has allowed more people to conduct simulations. Figure 1 shows an example of the tremendous growth in the relative fraction of the science and engineering literature that utilizes MD simulations. It is clear that MD simulations are playing an ever more important role in the research community. The one cautionary note in all this, however, is that, as the barriers for conducting MD simulations have lowered, the risk that inexperienced researchers will misuse the technique and generate meaningless (at best) or wrong (at worst) results has increased. As with any computational technique, a healthy dose of skepticism and vigilance is required among the broader community to ensure high standards are

3065

Figure 1. Number of publications per 1000 indexed by Web of Science having the keyword “molecular dynamics” as a function of time. Since its inception in the late 1950s, there has been a steady growth in the relative percentage of the literature that deals with MD, growing to 0.5% of the total by 2008. In more specialized journals such as Physical Review Letters, the fraction is 3.5%. Data obtained by Carol Brach, University of Notre Dame.

maintained. This is one more motivation for teaching MD as part of the curriculum. 4. Current Methods for Expediting MD Simulations The theme of continued extension to more precise properties of more complex systems is ever present in the field of MD. More recent extensions reflect the approach to the continuum scale for homogeneous systems of relatively small molecules. This observation emphasizes that mesoscale analysis can be bypassed for many systems of interest. We cite three examples here of how increasingly sophisticated theories of macroscopic behavior are being interfaced with the microscopic details of MD. In nonequilibrium MD (NEMD), gradients in molecular velocity, energy and concentration are artificially imposed and their evolution is interpreted through continuum theory to infer transport coefficients. In transition state theory, key transitions are identified and characterized at the molecular scale enabling extrapolation to macroscopic length and time scales. In coarsegraining, attention is focused on the molecular interactions that matter, while critically evaluating continuum scale assumptions like the no-slip boundary condition. These extensions reflect a necessary familiarity with the fundamental strengths and limitations of each theory at every length and time scale. Within engineering, this familiarity comes most naturally to chemical engineers. 4.1. Nonequilibrium Molecular Dynamics. The MD methods described in section 1 have all assumed that the fluid under study was sampling a series of configurations consistent with thermodynamic equilibrium in the ensemble under study. Transport properties can be obtained for a system at equilibrium by computing the integral of the appropriate time correlation function. The general formula is γ)

∫

∞

0

dt〈ξ˙ (t)ξ˙ (0)〉

(14)

where γ is the transport coefficient, ξ is the perturbation in the Hamiltonian associated with the particular transport property under consideration and ξ˙ signifies a time derivative. Integrals of the form given by eq 14 are known as “Green-Kubo” integrals.13 The integrated form of eq 14 results in an “Einstein” formula for γ of the following form

3066


2tγ ) 〈(ξ(t) - ξ(0))2〉

(15)

For self-diffusivity, ξ is the Cartesian atom position and the time correlation function in eq 14 is of the molecular velocities. For the shear viscosity, the integral in eq 14 is of the time correlation of the off-diagonal elements of the stress tensor. For the thermal conductivity, the integral is over the energy current, and for the electrical conductivity, the integral is over the electric current.13 An important implicit assumption in eqs 14 and 15 is that the time oVer which these expressions are eValuated is much larger than the correlation time of the Variable ξ. This assumption is often satisfied easily for simple liquids, where relaxation times are fast. For systems with sluggish dynamics, however, application of eqs 14 and 15 can result in incorrect results. Very early on it was appreciated that MD simulations could be used to mimic real experiments, and that transport properties could be computed by constructing computational analogues of experimental transport measurements. For example, attempts were made to compute viscosities by shearing fluids between isothermal “walls”,61 although results were not always satisfactory due to system size effects. Evans and Morris were among the first to develop formalisms for conducting homogeneous nonequilibrium molecular dynamics (NEMD) simulations.62 These methods have since been refined and now comprise an important and widely used tool in MD simulations. Among chemical engineers, Peter Cummings was one of the early adopters of NEMD simulations. In 1986, he and A. D. Simmons carried out NEMD simulations of dense fluid methane.63 In an influential paper coauthored with Denis Evans, a review of the basic algorithms was provided in Industrial and Engineering Chemistry Research; this paper effectively introduced nonequilibrium methods to a much wider audience.64 Since then, NEMD simulation techniques have been used to model a wide range of organic liquids,65,66 water,67 mixtures,68 liquid metals,69 and gases in nanopores,70,71 to name just a few applications. Some of the more recent applications have involved large molecules representative of lubricants,72-75 ionic liquids,76 and traditional alkali halide molten salts.77 A number of recent developments have been made to the original NEMD algorithms that make implementation easier. Transient nonequilibrium methods have been used to enable rapid estimation of transport properties such as thermal conductivity,78 viscosity,79 and diffusivity.70 Hybrid methods that combine Monte Carlo and MD have been used to study diffusion processes,80,81 although the reliability of these methods has been called into question.82 The design and implementation of these and other algorithms has greatly increased the number of systems and processes amenable to simulation with MD. 4.2. Reversible Coarse Graining. It is not unusual to abbreviate the potential model when seeking to extend MD to longer time scales. Indeed, the early hard sphere (HS) work could be seen in this light. More recently, Rapaport adopted a similar approach to study hydrodynamics using a truncated Lennard-Jones potential designed according to the Weeks, Chandler, and Andersen (WCA) procedure.83 The WCA potential is soft like the Lennard-Jones potential but shortens the range from 2.5σ to roughly σ. Since the computation time is proportional to the number of neighbors, and that is proportional to the volume, shortening the range of the potential accelerates the simulation by about 2.53, roughly an order of magnitude. Kremer and Grest used a similar approach to study the dynamics of tangent sphere polymer chains.49

Ultimately, however, we would like to go beyond studies of phenomenology and apply our best estimates of the intermolecular interactions whenever we apply MD. Even so, it would be legitimate to apply a coarse grained method as long as the parameters of the coarse grained potential were reversibly determined from the high resolution model. We refer to these methods as “reversible” because they are designed to provide a mapping of the high resolution coordinates at any point in time. In Kremer’s latest work, several atoms are lumped into a single site for integrating the dynamics.84 Nevertheless, their motions are constrained such that the overall molecules move to satisfy the behavior of the high resolution atoms and a well-defined algorithm relates the positions of the lumped atoms at all times. The SPEADMD method is reversible in the sense that it strips away the attractive interactions but retains the information to reapply them at any time. Quantitative accuracy in thermodynamic properties has been demonstrated for this approach through the application of thermodynamic perturbation theory (TPT). Unfortunately, there is no theory for transport properties comparable to TPT at this time. It is feasible to apply a perturbation series as an empirical tool, however.85 4.3. Coarse Grained Solvents. Other applications of MD that can approach the mesoscale are what we refer to as direct coarse graining. In these applications, MD simulations are performed in their usual way, but the potential models are selected such that solvent interactions are dramatically simplified (hybrid solvent) or ignored entirely (implicit solvent). The most common solvent of interest is water, of course, and applications focus on self-assembly of relatively large but dilute molecules. Explicitly integrating the dynamics of the water molecules (comprising roughly 90% of the solution) in addition to those of the solute, is very slow. It can be argued that the primary impact of the solvent in these applications is to serve as a continuous field that alters the interactions between the solutes. Hydrophobicity could then be represented as an attractive interaction between solute segments for example, rather than being inferred from the water-solute interactions. We illustrate these approaches with a few examples. 4.3.1. Implicit Solvent. Carol Hall and co-workers have recently focused their MD efforts on self-assembly of amyloids.86 Amyloid plaques have been implicated in Alzheimer’s disease. The amyloid proteins are treated by DMD with a united atom model that maintains realistic bond lengths and bond angles. The interactions of most protein segments are treated as purely repulsive. Hydrophobic segments are treated with square well attractions. Hydrogen bonding interactions are treated as cones of square-well energy at specific orientations. If a pair of segments approaches such that the cones are aligned, then the usual square well attraction is applied. If the bonded cones subsequently become misaligned by a change in orientation, however, the hydrogen bond is allowed to dissociate, ignoring the impact of the dissociation energy on the dynamics. Technically, this violates the conservation of energy principle, but the simulations are run at constant temperature using an Andersen thermostat, so energy is not exactly conserved anyway. (A more rigorous, albeit slower, DMD treatment of hydrogen bonding has been described by Liu and Elliott.26) The Andersen thermostat controls the temperature by altering the velocity of randomly selected sites such that the average temperature approaches the set point. With this approach, Hall and coworkers can study self-assembly of polyalanine that takes days to occur in the lab. If the method can predict how alterations in protein structure or drug treatments can disrupt plaque formation,


it may lead to radical breakthroughs in the treatment of Alzheimer’s disease. The dynamics of micelles and microemulsions have been of interest for many years.87,88 Micelles and microemulsions consist of relatively few molecules each, but these molecules are suspended in a solution that is mostly water. The interface is dynamic, but molecules enter or leave the microstructure relatively slowly once it is formed. The ability to observe the self-assembly of these microstructures could help to design surfactants and cosurfactants. Hamad and co-workers have adapted a DMD method with nonadditive interactions to drive the self-assembly.89 They observe dynamics that occur in real time over a range of milliseconds. Microemulsions form just one example of a self-assembling system, but these applications demonstrate the feasibility of MD and the range of length and time scales accessible for studying self-assembly in general. At this time, there is no direct connection between the effective potentials applied in these methods and the high resolution potential models that are consistent with measurements of fluid structure, vapor pressure, transport properties, and so on. In other words, these effective potentials are not reversible, so the connection is lost between the simulated dynamics and the real dynamics. If that remains the case, then there is little advantage of MD over MC methods. It may be possible to construct a reversible connection, but efforts at present are focused on the phenomenology. Another problem with implicit solvents is that the thermodynamics and dynamics may be qualitatively inconsistent with the results of explicit solvent simulations. Yethiraj and co-workers have shown that poor choices for the parameters of the effective potential can lead to such inconsistencies, while more sophisticated choices provide qualitative consistency.90 Efforts to make coarse graining procedures more rigorous continue to evolve. Robert Evans and co-workers developed a rigorous mapping between the partition function of a coarse grained model of colloids and the underlying detailed molecular model.91 Greg Voth and co-workers have worked on force matching algorithms.92 Martin Schoen and co-workers in Berlin have had success with coarse grained dynamics.93 Altogether, these methods show promise but have yet to achieve universal acceptance. The companion article in this centennial issue by Gubbins and Moore briefly describes a number of related procedures like dissipative particle dynamics and Langevin equations.94 4.3.2. Semi-implicit Solvent. To overcome the potential inconsistencies of an implicit solvent, one approach used by Malevanets and Kapral is to explicitly simulate the solvent molecules and their interactions with solutes, but to ignore the interactions between the solvent molecules with themselves.95 This modification is straightforward, requiring no change to the MD algorithm itself. Only the input files for the intermolecular potential need to be modified. Yet Padding and Louis were able to reproduce experimentally observed sedimentation behavior that could not be explained with an implicit solvent.96 Briefly, Brownian motion of the sediment particle leaves behind a void that cannot be sensed by an implicit solvent. On the other hand, the explicit solvent molecules rush to fill this void, colliding with the sediment and inducing coupled motions between solvent and solute. The solvent-solvent interactions are not necessary and eliminating their computation reduces simulation times by 1-2 orders of magnitude. Similar to Malevanets and Kapral’s method, Alder’s recent MD work has extended his methodology to tangent sphere chains in the presence of coarse-grained solvent molecules.97

3067

The solvent molecules interact with the chain molecule in the usual way, but they interact with each other by a direct simulation Monte Carlo algorithm. This algorithm assigns positions and conducts collisions stochastically instead of computing the collisions by rigorous DMD. It is found that explicit particles (and collisions) must be treated when Knudsen numbers exceed 0.1.98 Under these conditions, continuum methods and the no-slip boundary condition become unreliable. The stochastic approach overcomes the limitations of ignoring solvent influence while cutting the computational load by 99%. On the other hand, the combination of a stochastic algorithm with DMD goes slightly beyond direct coarse graining. Therefore, we draw the line here between methods related to MD vs those related to mesoscale modeling. Similarly, we relegate dissipative particle dynamics and lattice Boltzmann techniques to the discussion of mesoscale methods.99 4.4. Transition State Theory. What can be done when dynamical processes are much slower than the time scales accessible to MD, even when time-saving methods like implicit solvents are used? At present, standard or implicit solvent MD methods can only be used to probe dynamical processes occurring on time scales up to about 10-6 s. Recent progress is pushing time scales toward 10-4 s for somewhat specialized systems using clusters of highly vectorized graphics processor units (ie. GPU’s instead of CPU’s).100 Perhaps 1 order of magnitude more can be achieved with great effort, but processes occurring on time scales of seconds or hours will be outside the realm of standard MD for the foreseeable future. These processes, which we characterize as rare event processes, can be treated within the framework of activated process theory or more specifically transition state theory. In transition state theory (TST), the system or species moves between two states separated by a free energy barrier that is much larger than the thermal energy (thereby making such transitions rare events). To make matters concrete, assume that TST is used to track a particle diffusing over a large energy barrier. (TST can also be used to examine conformational changes in molecules, adsorptiondesorption kinetics, and reaction kinetics. One of the key requirements, however, is that a reaction coordinate must be defined, which is not always an easy task). TST assumes that the particle is coupled to a heat bath and that a well-defined dividing surface separates the particle from its initial and final states. In classic TST, it is assumed that when the particle reaches the stationary point separating the two states, it has sufficient kinetic energy to overcome the barrier, and the crossing event is successful. If the coupling between the molecule and the thermal bath are too weak or strong, however, the particle may recross the barrier (fail to thermalize in the destination state) or fail to cross the barrier at all. Such “dynamical recrossing” events lead to smaller rate constants than predicted if these events are ignored. MD can be used to compute the free energy barriers separating the initial and final states, and can also be used to compute dynamical recrossing probabilities. Within the chemical engineering community, we can mention of few studies that used TST combined with MD to study diffusion processes. June and Theodorou used TST to estimate diffusion coefficients of small molecules in zeolite pores at infinite dilution.101 Tunca and Ford extended these ideas to study diffusion in porous materials at finite loading.102 Smit and co-workers have applied dynamically corrected TST to a wide range of sorbates and zeolites.103 Their studies provided a fundamental understanding of the nature of diffusion in these materials. TST has been combined with other coarse-graining methods to simulate the

3068


Figure 2. Arrhenius plots of simulated desorption rates for n-alkanes ranging from pentane to hexadecane. (Reprinted with permission from ref 107. Copyright 2006 American Institue of Physics.)

Figure 3. Simulated and experimental desorption energies for a series of n-alkanes on graphite. Experimental results labeled 17 and 22 are from refs 158 and 159, respectively. (Reprinted with permission from ref 107. Copyright 2006 American Institue of Physics.)

dynamics of very large molecules inside zeolites.104 Fichthorn and co-workers used MD and TST to examine the thermal desorption process of alkanes from a gold surface.105 Finally, when recrossing probabilities are found to be non-negligible, methods reviewed by Santiso and Gubbins can be considered.106 These and many other applications demonstrate the versatility of MD but also point out how important it is to consider combining MD with other techniques. Increasingly, advances in understanding long-time dynamical processes will involve the judicious use of MD coupled with other techniques such as TST to obtain meaningful results. This coupling of methods is difficult and requires care on the part of the researcher. 5. Chemical Engineering Applications Having discussed the history of MD and the development of the methods that enable realistic systems to be simulated, this section covers some examples of the application of MD to systems of direct interest to chemical engineering. The literature is vast and it is impossible to discuss even a small fraction of the work. The particular topics chosen for emphasis are thus intended to be simply a representative sample of the field. Our intent is to articulate how MD is a useful engineering tool at present as well as developing tool for the future. 5.1. Organics, Mixtures. In section 2.2.1, work was described in which simple hydrocarbon liquids were simulated. Recent work has focused on much more complex hydrocarbon systems and has shed light on the behavior of these systems at surfaces as well as in the bulk. As an example, Becker and Fichthorn used an accelerated MD method to examine the desorption of alkanes ranging from pentane to n-hexadecane from a graphite surface.107 This technique builds on the TSTinspired work of Fichthorn described earlier. The process of alkane desorption from a graphite surface is slow, such that, even with advanced computers and the latest methods, standard MD is incapable of probing the phenomenon. To overcome this, Fichthorn and her group use accelerated MD methods combined with a powerful technique called umbrella sampling and TST to examine rare event processes. They computed Arrhenius plots for the desorption of eight different alkanes ranging from pentane to n-hexadecane, as shown in Figure 2. Figure 3 shows a comparison of simulated and experimental desorption energies as a function of chain length. The agreement with experiment is excellent.108,109 Similarly impressive simulations have been conducted on bulk hydrocarbon systems of great complexity. For example, Zhang and Greenfield carried out MD studies of a model asphalt composed of a mixture of aliphatic and aromatic species.110

Figure 4. Comparison of experimental and simulated viscosities for various asphalt materials. Open and filled squares are simulations of two different model asphalts. Filled circles and triangles are different experimental viscosity measurements. (Used with permission from ref 158.)

Viscosities were computed at 443 K via Green-Kubo and Einstein methods, while at lower temperatures where the dynamics are much more sluggish, estimates were made by computing rotational relaxation times and then estimating the low temperature viscosity via a Debye-Stokes-Einstein approach. As shown in Figure 4, both the trend with temperature and the absolute values of the viscosity were in good agreement with experiment. 5.2. Diffusion in Nanoporous Materials. Nanoporous materials such as zeolites, porous carbons, nanotubes, and metal organic framework materials see wide application in separations and catalysis. Not surprisingly, these systems have come under intense study by chemical engineering modelers. Previously, we mentioned some of the work done on diffusion in these systems using TST, but the earliest work focused on calculating the nature of the fluid in the pores. One of the earliest MD studies of these systems was carried out by Magda, Tirrell, and Davis in 1985.111 A simple Lennard-Jones fluid was confined between structureless walls, and MD was used to compute local densities, solvation forces, interfacial tension and diffusion coefficients. Figure 5 shows how the fluid orders in two layers near the walls for a given wall spacing. Multiple layering was observed for pores having wider spacing. Later work has focused on adding atomistic details to specific porous systems. Studies have been conducted on zeolites112 and other porous materials.113 Comparisons have been made between self-diffusivities calculated from MD and those measured using NMR,114 neutron


Figure 5. Density profile between pore walls separated by three liquid diameters. The fluid shows enhanced density in two layers near the walls. “Reference 9” in the legend refers to the following: Snook, K.; van Megen, W. J. Chem. Phys. 1980, 72, 2907. (Reprinted with permission from ref 111. Copyright 1986 American Institute of Physics.)

and have opened up the possibility of using these materials in selective membrane separation systems. The effect of pore entrance and exit effects have also been determined from MD simulations,71,118 and predictions of anomalous “single file diffusion” have also been made using MD.119 These are just a small fraction of the MD studies done on nanoporous materials, but even this small sampling demonstrates that simulations are providing both quantitative predictions of transport properties for specific systems along with important mechanistic insights into the nature of molecular motion in these systems. Simulations have also played a key role in clarifying the difference between the types of diffusion coefficients that people measure in experiments. Much of the confusion and controversy among different experimental measurements can be traced to differences in the “diffusion coefficients” that were measured. Self-diffusivities, transport diffusivities, and corrected diffusivities can all have different dependencies on concentration, temperature, and system parameters. Since these quantities can be defined and calculated exactly in a simulation, these differences have been demonstrated by the use of MD simulations and greater clarity has arisen in the community as a result.117 5.3. Molten Salts and Room Temperature Ionic Liquids. Molten saltsstypically alkali halides at elevated temperaturess have been studied for a long time using molecular simulations, as evidenced by the pioneering work of Woodcock described in section 2.2.1.33 Simulations afford obvious advantages when the system of interest is at a high temperature state difficult to reach experimentally. The original force field functional forms used to model alkali halides is due to Huggins and Mayer and has the following form120 Uij(rij) )

Figure 6. Self-diffusivity for single-component CH4 and H2 and equimolar CH4/H2 mixtures inside a (10,10) single wall nanotube at room temperature, computed from equilibrium MD. Diffusion coefficient from EMD simulations of single-component CH4 in the zeolite ZSM-12. (b) Same as part a but showing the transport diffusivity. (Reprinted from ref 118. Copyright 2004 American Chemical Society.)

scattering,115 and conventional macroscopic techniques.116 Simulations predicted that diffusion in carbon nanotubes would be extremely high, much higher than that observed in porous materials having “rougher” surfaces.117 Figure 6 shows that diffusion coefficients of methane and hydrogen in nanotubes can be orders of magnitude greater than that in zeolites.117 These predictions have been subsequently confirmed experimentally

3069

qiqj Cij Dij + Bij exp(-Rijrij) - 6 - 8 rij rij rij

(16)

where i and j can be either a cation or anion and Uij is the potential energy of species i due to species j which are separated from each other by a distance rij. This is similar to the form of eq 1, although the repulsive part is exponential and there is an additional term in the attractive van der Waals terms. These latter terms were parametrized by Mayer,121 while Tosi and Fumi122 developed parameters for the repulsive part of this function to reproduce the solid phase properties of alkali halides having the rock salt structure. Many authors subsequently used this parametrization (or variations of it) to simulate a wide range of alkali halides in both the solid and molten states. Sangster and Dixon have reviewed the early work in this field.123 Work by Galamba and co-workers suggests that shear viscosities and thermal conductivities of molten NaCl and KCl are over predicted by 10-20% with this force field.124 While conventional alkali halides have been simulated for over thirty years, a new type of molten salt has recently emerged in the past decade that has sparked tremendous interest among modelers. “Room temperature” molten salts or “ionic liquids” are terms given to pure salts that tend to be liquids at temperatures around ambient. The cations and anions are typically multiatom and asymmetric, thereby leading to lower melting temperatures. As a result, modeling these systems is more challenging than simple salts, but perhaps even more important. This is because there are almost an infinite number of cation and anion combinations that can be combined to make an ionic liquid. Figure 7 shows an example of just a few of the classes of cations and anions under investigation. While there

3070


Figure 7. Examples of some of the cations (a) and anions (b) under investigation as ionic liquids. Image courtesy of Robin Rogers.

are a relatively small number of akali halide salts, the composition space for ionic liquids is huge. Simulations offer the possibility of gaining insight into molecular structure property relations that can help direct synthesis efforts toward compound classes having desired properties. In 2001, Hanke, Price, and Lynden-Bell were the first to conduct an atomistic simulation of compounds resembling those in Figure 7.125 They used MD to model the crystalline state of 1,3-dimethylimidazolium chloride ([C1mim][Cl]), 1,3-dimethylimidazolium hexafluorophosphate [C1mim][PF6], 1-ethyl-3methylimidazolium chloride ([C2mim][Cl]), and 1-ethyl-3methylimidazolium hexafluorophosphate ([C2mim][PF6]). They also modeled the liquid state of [C1mim][Cl] and [C1mim][PF6], both of which are relatively high melting substances. They computed molar volumes, average energies, liquid structure in the form of radial distribution functions, and mean square displacements as a function of time. The relatively short simulation times (ca. 100 ps) are now known to be too short to obtain reliable self-diffusivities. Nevertheless, this work opened the door to the possibility of using classical force fields and MD to simulate the properties of this new class of material. Soon after this, several other papers appeared in which imidazolium-based ionic liquids were simulated. Margulis, Stern, and Berne conducted MD simulations of [C4mim][PF6] using an explicit atom model that included fully flexible bond lengths and bond angles, in addition to a dihedral potential.126 The MD simulations utilized 200 ps of equilibration time and 50 ps of production time. Computed liquid phase densities were in very good agreement with experimental values, but again the simulations were too short to obtain reliable dynamics properties. Nevertheless, these authors identified anomalous and complex dynamical behavior, indicative of a supercooled liquid. This anomalous dynamical behavior turns out to have been an extremely important finding as related to the dynamical properties of ionic liquids, and was first appreciated in the modeling community.

Morrow and Maginn also carried out a molecular dynamics study of [C4mim][PF6], using an explicit atom and fully flexible model with a force field similar to that in eq 1.127 Their equilibration times were from 700-1000 ps, and production runs of 4 ns were carried out. Consistent with the work of Margulis, Stern, and Berne, they found that long simulations were required to obtain dynamical properties because the ionic liquids tend to exhibit behavior reminiscent of supercooled liquids at ambient conditions. While densities and other static properties converged within 50-100 ps, rotational time constants were calculated to be much longer than even the 4 ns simulation run time. Since these early studies just eight years ago, there have been over 200 papers published in which MD has been used to examine ionic liquid behaVior. Most of the papers use standard methods and probe liquid structure, dynamics, and solutessolvent interactions in different ionic liquids. It is simply impossible to even list all the papers heresinstead, the interested reader is referred to a recent review for more information.128 This explosion of interest is the result of the happy coincidence that an entire new class of material has been discovered at the same time that the force fields, methods, and computational horsepower described earlier in this article have come together to enable MD to make accurate and efficient property predictions of these materials. Had ionic liquids come of age ten years earlier, these developments would not have been possible. Thus it is that one should expect MD to play an increasingly significant role in discovery efforts of new materials, rather than postpredictive studies of existing materials or even studies aimed at refining materials that are well-established. Brand new classes of materials are the most fertile areas for future growth of MD. 5.4. Physical Properties. To say that MD has evolved into a chemical engineering tool implies the availability of applications ready at present, as well as those that may come in the near future. Concerning physical property predictions at present, molecular simulation (including MC for equilibrium properties) is the best available technology in several cases. Longer term, we can extrapolate by adopting a broader perspective on what is considered a physical property. For example, vapor pressure, density, diffusivity, and viscosity are obvious physical properties in the tradition of “Properties of Gases and Liquids.” More broadly, glass transition temperature, yield strength, and octanol-water partition coefficient could be considered as physical properties. Yet more broadly (and for the future), substrate binding, β-sheet content, and phase order (e.g., smectic, lamellar, bilayer, cylindrical) could be considered as physical properties. In all cases, the procedure is to gather data, make predictions, and refine the model. This procedure is not limited to traditional physical properties, but there is certain logic to doing the obvious ones first. It is not trivial to compete with the accuracy of traditional physical property methods. These have been refined over many years using all available data. Just making the initial predictions for all available data requires a prodigious effort with MD. Refining the MD predictions can require similar effort many times over. Nevertheless, the effort is necessary before it is reasonable to say that MD has an advantage over traditional methods. If accurate predictions have been demonstrated for the property of interest using a test case comparable to that of a traditional physical property method, then this is a property that is available at present. If promising results have been obtained for a few test cases, we classify the property as available in the near future. MD predictions have one saving


grace in that properties can be predicted self-consistently using a force field that was optimized based on other properties. When predictions are accurate without correlating to a specific property and compared to predictions of a traditional method that are outside its training set, we consider that to be a comparable basis even if the number of demonstrations is relatively few. 5.4.1. Common Potential Models (aka. Force Fields). The term atomic force microscope (AFM) is a misnomer. Atomic forces occur on the nanometer scale, not the microscopic level. For example, the diameter of a methane molecule is roughly 0.4 nm. Furthermore, the breadth of an AFM tip incorporates over a hundred atoms interacting with the surface, not individual atoms. Finally, the resolution is not very precise in the range of 0.3-1.5 nm, partially because the tip tends to deform the atoms of the surface and the tip itself. On the other hand, accurately measured physical properties coupled with MD can provide precision of (0.1 K at every position in the range 0.3-1.5 nm. Various physical properties are sensitive to varying details of the potential model. Liquid density is most sensitive to the softness of the potential in the range of 0.3-0.4 nm, especially compressed liquid densities. Vapor pressure, however, is most sensitive to the detail of the potential in the range of 0.5-0.8 nm. Phase equilibrium data can resolve site-site interactions with sensitivity analogous to a differential pressure gage compared to an absolute pressure gage. Transport properties trend toward the same sensitivities as liquid density, but also reflect the attractive potential. Altogether, nothing comes closer to an “atomic force nanoscope” than the combination of MD and physical property measurements. Physical property predictions require an assumed potential model and potential model characterization is based on physical property predictions. It is a problem requiring an iterative solution. The choice of physical properties for characterizing the potential model reflects the intended applications of the developer. The accuracy of physical properties of traditional chemical engineering interest is dominated by the characterization of the nonbonded interaction potential. Once the physical properties of interest have been identified, the procedure seems straightforward but leads to solutions that are neither unique nor universal. For example, the CH3 and CH2 potential models could be characterized with n-alkanes. Then branched alkanes could be included, then alcohols, and so forth. Subjective judgment is required even with the n-alkanes. Should the hydrogen atoms be represented explicitly (EA model) or implicitly in the form of united atoms (UA) that comprise single “CH3” and “CH2” sites? Precisely how should the dihedral angle potentials be evaluated? When polar molecules like alcohols are included, judgment is required to characterize the polarity and the consequential hydrogen bonding. In most cases, polarity is characterized with a distribution of partial charges bonded at specific distances and orientations. There is no universal procedure for assigning this distribution. It is tedious work that has often seemed to be a distraction from the authors’ primary interest. This has led to the development of numerous potential models, each with their own proponents and none that are universally accepted. Early efforts were based on liquid density and solubility parameter near room temperature. These could be evaluated economically with a single simulation, whether MD or MC. As new functionalities were introduced, their potential models could be characterized with limited data. Sometimes data would be estimated if none were available at room temperature. Liquid densities are sensitive to site diameters leading to precise values, but the solubility parameter tends to be relatively insensitive to

3071

details of the potential model. Furthermore, the problem is grossly underspecified when one considers the number of decision variables involved in specifying a distribution of partial charges. To alleviate this problem, radial distribution functions are typically evaluated to display a characteristic peak indicating hydrogen bonding, but the assessment is qualitative. This is the general approach that has been applied in developing the CHARMM129 and COMPASS130 potential models. The limitations of this approach can be appreciated when vapor pressures are predicted. The solubility parameter characterizes the heat of vaporization, but not the vapor pressure. Typical deviations for vapor pressure predictions with these potential models are a factor of 2 or more. Vapor pressure tends to be very sensitive to the characterization of the nonbonded interactions, making it the obvious choice as the physical property for characterizing potential models of interest to chemical engineers. Several general potential models have been developed based on this choice, including the transferable potentials for phase equilibria (TraPPE) model,131 the anisotropic united atom (AUA) model,132 and the SPEADMD model. All three models are united atom models, although EA versions of the TraPPE model have recently emerged.133 The TraPPE model was originally developed to address trends in critical properties more than vapor pressure directly. Therefore, systematic deviations are encountered for vapor pressure, especially for the (hydrocarbon) site types characterized in the early stages. The optimization procedure for the TraPPE model involves multiple simulations with the complete potential model, making it very difficult to refine the model at a later stage. The optimization procedure with the AUA model takes advantage of thermodynamic derivative properties to accelerate refinement of the potential. The sensitivity of the vapor pressure to the potential parameters can be expressed in terms of derivative properties that can be evaluated during the simulation. The optimal potential parameters can then be predicted based on extrapolation from one specific simulation. In principle, the extrapolation can be refined at a later date if further information indicates a need for a change. For example, suppose that the original characterization of the branched CH3 site type was based on isobutane, but later analysis showed substantial deviations for isopentane, isohexane, etc. The AUA approach would permit an estimate of the optimal potential parameters for all compounds without resimulating isobutane. The SPEADMD approach is something in between. The site diameters of the SPEADMD model are fixed for a particular potential model (e.g., a 4-step potential), but the attractive part of the potential is entirely variable and can be fully optimized for the entire database at any time. Naturally, this leads to improved precision for vapor pressure but sacrifices precision for liquid density. These are the potential models that have been applied most often by chemical engineers and are cited often in the discussion below. 5.4.2. Vapor Pressure and Vapor-Liquid Equilibria. Vapor-liquid equilibrium (VLE) is of prime importance in chemical engineering and accurate vapor pressure is essential to accurate VLE. Noting that vapor pressure represents the VLE of the pure fluid, inaccurate end points must lead to inaccurate interpolations. Furthermore, vapor pressure is sensitive to small variations in the intermolecular potential models, changing by 2-3 orders of magnitude in the typical temperature range of interest. The TraPPE model provides accuracy for vapor pressure that is around 30%, with predictions for hydrocarbons being systematically high by about 15%. Roughly 60 compounds are included in the training set for correlating the TraPPE param-

3072


eters. The accuracy of the AUA method is roughly 15% based on a database of about 100 compounds. The SPEADMD (4step, “2580”) model correlates vapor pressure to within roughly 10% accuracy based on a database of about 500 compounds.25 These can be compared to roughly 35% accuracy for traditional prediction methods based on boiling temperature and heat of vaporization.134,135 Predictions of vapor pressure and VLE have been pursued by both MC and MD. We focus here on the MD aspect, but conclusions are similar in both methods. The trend is illustrated by Vahid et al. for a relatively broad database.136 The conclusion is that the accuracy for predictions of VLE approaches 10% when the vapor pressure is predicted by transferable potential models and 3% when the potential model is customized to match the vapor pressure of specific components in the mixture. The latter is comparable to the accuracy of a typical equation of state like the Peng-Robinson model, which is similarly customized to match vapor pressure through its use of acentric factor. Comparing the accuracy for predictions that have not been customized is difficult. The Peng-Robinson model does not include an extension for predicting vapor pressure. Since vapor pressure predictions based on transferable potential models are most accurate, we would suggest that MD or MC should be the recommended method for predicting VLE when neither vapor pressure nor mixture data are available. 5.4.3. Liquid Density. To compute the density in an MD simulation, the temperature and pressure of the system are specified and isothermal-isobaric MD is run. The density is then found by computing the average volume of the system. Liquid density is a property that seems very simple, but influences other properties in a significant way. For example, the accuracy for predictions of viscosity and glass transition temperature can be directly related to accurate density predictions. TraPPE and AUA provide accuracy of roughly 1% for density, while accuracy for the SPEADMD (2580) model is on the order of 2-3%. For comparison, the accuracy of the LeeKesler traditional model is roughly 2-3% while the accuracy of the Hoy group contribution method is roughly 6%. 5.4.4. Self-Diffusivity. Self-diffusivity is a physical property of limited interest by itself. Relatively few experimental data are available and mass transport applications rarely require diffusivity estimates over broad ranges of temperature and density. On the other hand, it is in many ways the simplest transport property, forming an initial basis for simultaneous study of equilibrium and transport properties. Its evaluation for systems such as sorbates in zeolites is of great interest for applications in catalysis and separations. Self-diffusivities can be computed from equilibrium MD using eq 15. A number of correlations of bulk phase self-diffusivity are based loosely on MD simulations. The most widely applied method is the correlation of Liu et al., based on adaptation of theory and simulations of diffusivity in square-well and LennardJones spheres.137 Liu et al. adapt the simulation results as the functional form for correlating properties, however, rather than simulating molecules directly. Liu et al. provide methods for correlating diffusivity in terms of the effective sphere diameter, the attractive energy, and a parameter characterizing the softness or free volume of the molecule.137 Their correlations provide accuracy of about 5% for nonassociating compounds and 20% for associating compounds. The trends in their regressed parameters do not correlate with physical intuition, however, and attempts to use this correlation as the basis for making predictions have led to large deviations.

The SPEADMD method has been evaluated for a database of roughly 50 compounds, comprising the vast majority of known bulk self-diffusivity data.27 The method applies simulations of the reference molecules to estimate the diffusivity in the absence of attractive forces then applies an empirical correction for the influence of attractive forces. Deviations for n-alkanes are roughly 10%, including predictions over the entire range of molecular weight below the entanglement threshold (MW ∼ 1000). Using the n-alkanes as a reference, diffusivity of nonalkanes can be predicted by simulating the reference fluid and extracting an effective value for the n-alkane chain length. Predictions for nonassociating compounds have deviations near 20% and predictions for associating compounds are near 40%. Briefly, the improved physical basis of MD for the specific molecular shape with all its branches and rings permits predictions that approach the accuracy of purely empirical correlations. This is encouraging for the general prediction of transport properties. 5.4.5. Viscosity. The shear viscosity of a bulk fluid is a critical engineering property and can be computed from either an equilibrium or nonequilibrium MD simulation using the methods described in section 3.1. Given its practical importance, it is not surprising that it is a frequently calculated quantity. Interest in the viscosity of solvents, lubricants, and process fluids is high, and MD studies have been especially important in helping understand viscosity under extreme conditions where experimental investigations are difficult. There are too many examples in the literature to provide a comprehensive review, but a few of the significant papers within the chemical engineering field should suffice to give a sense of what is possible. One of the earliest MD calculations of the viscosity of “large” alkane molecules occurred in 1992. Berker and co-workers used nonequilibrium MD to compute the viscosity of n-hexadecane.138 A few years later in 1996, Cummings and co-workers used equilibrium and nonequilibrium MD to compute Newtonian and non-Newtonian shear behavior (including the Newtonian shear viscosity) of a series of long linear alkanes, the largest being tetracosane.139 At about the same time, Allen and Rowley carried out a nonequilibrium MD study to examine the effect of branching and type of force field on computed viscosity.140 Equilibrium MD simulations of alkanes in the gas and liquid phase, up to near the triple point, were carried out by Dysthe, Fuchs, and Rousseau.141 Simulations of realistic lubricant molecules under extreme pressures were also carried out around this time,142,143 and not only were predictions of viscosity made, but a deeper fundamental insight into how viscosity varies with pressure and free volume was also developed. The simulations showed that the commonly used “Barus equation”, which predicts that the logarithm of the shear viscosity varies linearly with pressure, is incorrect. Instead it was shown that the viscosity correlates much better with the available free volume. Systems other than alkanes have been studied. Chemical engineers Wheeler, Fuller, and Rowley showed how to properly implement the Ewald sum in a simulation with the “LeesEdwards” sliding boundary conditions used in the conventional nonequilibrium MD approach.144 This opened the door for a number of other simulations of polar fluids, such as the alcohol simulations described in section 4.5. Finally, we note some recent methodological developments that have been designed to make viscosity calculations easier. The most difficult part of a conventional nonequilibrium MD simulation is that the flow field (shear rate) is imposed on the system and the response (shear stress) is computed. Stresses


are one of the most difficult quantities to compute in a simulation because they fluctuate tremendously. Mu¨ller-Plathe came up with an ingenious method for computing viscosity that essentially “reverses” the conventional method. Originally developed for thermal conductivity calculations, this “reverse nonequilibrium” MD procedure imposes the hard-to-compute stress by periodically exchanging momenta between molecules and the easy-to-compute shear rate is determined by measuring the resulting velocity profile.145 This method does not require alteration of boundary conditions and so can be implemented easily for charged as well as uncharged systems. Although all of the methods outlined here are steady-state approaches, it is also possible to compute viscosity from transient simulations where the relaxation of the system to a perturbation is determined and related to the viscosity.146 These transient methods may be inappropriate for large molecules with long relaxation times, but their advantage is that they can yield viscosity estimates in a fraction of the time required for conventional steady-state methods. 5.4.6. Polymer Properties. The slowness of MD is an obstacle to the study of polymer properties, but not a barrier. This is especially true when the properties of interest are localized, as in the case of glass transition or diffusivity of small molecules through a polymer matrix. At the far end of what can be achieved, MD can also be applied to study polymer dynamics on the time and length scales of polymer entanglement. For current applications, we focus on methods that have been illustrated with several polymers and are on the verge of being competitive with existing empirical methods. Predictions of polymer properties underscore the difficulty of developing MD methods that are competitive with existing empirical methods. Simulations of polymers require greater computational effort. In the case of gas penetrants, for example, the equilibrium concentrations are very small. This necessitates simulating a large number of polymer atoms just to simulate a few gas penetrant molecules in order to get a reasonable average. Therefore, the system size is large. Furthermore, the dynamics of polymer systems are relatively slow, necessitating a long simulation. Long simulations of large systems are extremely expensive. On the other hand, academic motivation is very weak for developing MD methods that would be conclusively superior to empirical methods, as demonstrated with similarly large databases. The computational expense for polymer systems exacerbates this problem. As early proofs of principle, work on polymers focused on determining glass transition temperatures (Tg).48 The glass transition is an important property in polymer processing, and it determines the kinds of applications that may be appropriate for a particular polymer. Early work on polymer simulations also treated diffusion of penetrants like oxygen through polymers. Penetrants are important in the packaging industry. a. Glass Transition Temperature. The basic idea is to simulate the polymer at constant pressure and track the volume (or some surrogate sensitive to the volume) while reducing the temperature. A characteristic break occurs at Tg. This was the approach pioneered by Rigby and Roe.48 The method has been applied to dozens of polymers and is accurate to roughly 50 K. A systematic problem with the method is that cooling rates in MD are typically 20 K/ns, compared to 5 K/min experimentally, a factor of 1011 different. The slower experimental cooling rates allow time for the polymers to optimize their packing at each temperature. As a result, the Tg’s estimated by MD with this approach are systematically lower than the experimental values. Unfortunately, the deviation from experiment varies with

3073

polymer structure, so a simple subtraction of a fixed quantity does not entirely resolve this inaccuracy.147 As an alternative, researchers have computed correlation times of segmental dynamics as a function of temperature. These correlation times show a distinct increase when passing through Tg, as validated against neutron scattering and other experiments. For example, Theodorou and co-workers adapted this approach to atactic polypropylene.148 For comparison, a typical group contribution method for predicting Tg provides accuracy of ∼20 K, as demonstrated by van Krevelen for over 500 polymers.149 Van Krevelen’s analysis also shows that experimentally measured Tg’s vary by ∼10 K. In this context, we must conclude that MD methods provide insights about the relative impacts of molecular structural variations, but are not yet competitive with empirical correlations in an absolute sense. b. Gas Penetrant Permeation. Gases permeate polymers in small amounts but their implications can be large. Particular examples can be cited in food packaging. The effusion of CO2 can limit the shelf life of soda stored in polyester bottles. Similarly, the infusion of O2 can limit the viability of polymer packaging for everything from beer to beef. Controlling the permeation of these penetrants requires formulations that involve a number of barriers and reactants. These considerations help to motivate studies of the molecular structures that can improve formulations. Similar to glass transition temperatures, permeant diffusivity is sensitive to the precision of the polymer’s estimated density. The density controls the amount and size of voids in the polymer, through which the penetrants must diffuse. An early review of the field was performed by Suter and co-workers.150 Boyd also cites several references addressing penetrants.147 A more recent monograph on polymer membranes includes a discussion of the status of this field.151 Overall, a limited number of polymers have been treated so a comprehensive analysis is not possible. Estimated diffusivities are generally accurate to within an order of magnitude. On the other hand, empirical methods of estimating this property are somewhat underdeveloped as well, making MD a competitive methodology by default. Progress in this area might require addressing how composite formulations would impact predictions (e.g., 0.1% nylon in polyester). The length scales of composites would not be accessible by direct MD simulation. 5.5. Industrial Fluid Properties Simulation Challenge. A key milestone in the evolution of molecular dynamics (and simulations in general) as a powerful tool in chemical engineering was the establishment of the Industrial Fluid Properties Simulation Collective (IFPSC) competition in 2002.152 The idea of the competition was to have certain properties measured and withheld from publication for a short period of time. Modelers would then be challenged to predict the properties as best they could over a period of a few months. Competitors submitted their entries to the organizing committee, who then announced the winner at a symposium during the fall AIChE meeting. The goal of the challenge was to “promote the use of molecular simulation methods to predict materials’ properties of industrial relevance”.152 A secondary objective was to establish reliable comparisons between the available methods so that modeling practitioners could choose the best methods for a given problem. A recent publication describes more fully the history of the competition and its origins.153 In the first competition, three “problems” were posed to the group. In the first problem, participants were asked to predict the Px curve for a mixture of dimethyl ether and propylene at

3074


Participation in the competition has been variable and narrow relative to the global simulation community. For the most recent competition, only four groups participated (one from Europe, one from China, and two from the United States). Considering the large number of groups conducting simulations worldwide (see Figure 1), this is a small number. Part of this is due to the fact that these problems are essentially unfunded research and each takes a significant investment in time and resources to complete. While a small monetary award is given to the winners, the main incentive for participating is an altruistic desire to see the field advance and to help in that effort. Despite the relatively low participation levels, the challenge has helped identify weaknesses in the field and has highlighted areas where simulations are mature and able to yield high quality property predictions. For this challenge to have a greater impact, participation levels must increase. One mechanism for doing this is to ask agencies who fund molecular modeling research to help fund (and then encourage) grantees to participate. Figure 8. Comparison of predicted and experimental viscosities at 0.1 MPa for five different polyhydric alcohols. Predictions are from three different groups, all of whom used MD. (Reprinted with permission from ref 154. Copyright 2007 Elsevier.)

a series of temperatures and compositions. They were also asked to predict the pressure and composition of the azeotropic point for a mixture of nitroethane and propylene glycol monomethyl ether, as well as the bubble point pressure for this mixture at two compositions. The second problem asked for densities of a wide range of compounds and mixtures at two different state points. The compounds included alcohols, chlorinated hydrocarbons, pyridine, an alkanol amine, and salt solutions. The third problem asked for the viscosity of nonane, isopropanol, and their mixtures at different state points. The challenge was open to any group, and the first competition had entries from three companies, one national lab, and several university research groups. The first problem was attempted by three teams, two of whom used Monte Carlo and the third used an activity coefficient model based on the COSMO approach. No MD simulations were used. For the second problem (by far the easiest of the three), both MD and Monte Carlo methods were used. For the last problem, three entrants participated and all used MD for their predictions. The results from this first challenge highlighted both the potential of simulations and the areas where more developments were needed. Some of the entries missed the mark badly, and the fact that only one team was able to participate in all three problems showed that methods are still specialized and that it is not easy to compute a wide range of diverse properties. There have been three subsequent challenges, with the results of the fourth IFPSC challenge recently published.153 The results have been mixed. Some properties have been computed with remarkable accuracy, while others have been much more difficult to predict. MD has been the workhorse for predicting transport properties such as viscosity and thermal conductivity, while Monte Carlo has been the dominant method used for phase equilibrium. Figure 8 shows an example from the third challenge in which the viscosity of five different polyhydric alcohols was to be predicted at various state points.154 There is substantial deviation among simulation predictions, and only one entry was able to predict the viscosity at this state point to within the experimental uncertainty. This points out the need for better force fields for transport properties in addition to the development of more effective and easily used methods for computing properties such as viscosity.

6. Conclusions and Future Outlook MD can serve as an engineering tool in applications that range from mundane to inspirational. MD has reached a level of maturity and is technically ready to serve as a tool as evidenced by responses to the IFPSC. From an engineering perspective, MD is near the point of superseding traditional methods of physical property estimation. The assumption of transferable potentials is comparable to a group contribution correlation, but the foundation of the functional forms for the property correlations is as sound as Newton’s laws. Furthermore, a unified set of parameters describes all properties (i.e., the potential parameters). This compares to having one set of parameters for activity coefficients, another for vapor pressure, another for viscosity, and another for thermal conductivity. The unified approach is capable of resolving inconsistencies between trends in viscosity and liquid density, for example. The unified nature of MD is also a two-edged sword because potential models are always approximate. If a detail of the intermolecular potential is lacking and critical for a particular physical property, then predictions might be inaccurate. On the other hand, resolving such discrepancies would suggest ways of systematically refining the potential model to a high degree of resolution. MD then becomes a platform for continuously improVing our insights into molecular interactions. Admittedly, the slowness of MD precludes replacing group contribution methods in a combinatorial approach to molecular design.155 Nevertheless, MD should guide the functional forms used to train the correlations in order to enhance their capacity for extrapolation. Going forward, physical property correlations (including mixture VLE) should be informed to an increasing extent by molecular simulations, not just for general functional forms but for direct computation of the property of interest based on a well-calibrated potential model. Direct simulation is already the superior method for predicting vapor pressure, liquid density, diffusivity (bulk and nanoporous), and viscosity. MD can guide mesoscale research in a similar fashion. At its simplest level, fundamental relations can be tested and extended while retaining much of the traditional formalism at the continuum level. Characterizing the transitions from noslip, to partial slip, to molecular boundary conditions provides an example of how this connection can often be made. For example, it might not be useful to apply MD to study protein flow entering a microtube if the phenomena could be well characterized by a continuum approach with a viscosity from


MD and a slightly modified boundary condition. It is important to apply MD carefully in such situations, thinking critically about what distinctive information is to be gleaned from time scales longer than a microsecond and sizes larger than 10 000 interaction sites. System scales like these are necessary for many biological, polymeric, and self-assembling systems. Adapting to specific applications may require specialized techniques for each case, as illustrated with a few examples. Ignoring solvent-solvent collisions can cut simulation time by a large percentage, while avoiding the pitfalls of an implicit solvent. Langevin equations or similar models can be used to describe the long-term behavior after the phenomenological coefficients have been characterized using relatively short-term MD simulations.99 Pretabulation can accelerate simulations of fluids in solid pores. Reversible lumping of multiple interaction sites may be viable in some cases. Parallel programming has had a major impact on MD simulations, and new simulation paradigms such as the use of graphics processing units (GPUs) will continue to expand the scope of MD. Approaching the mesoscale will require at least some degree of coarse graining. Systematic progress requires that the coarse graining be verifiably consistent with the rigorous molecular dynamics. Achieving that consistency will require combinations of theory and simulation that strip the simulation component to its minimal form. In the National Academy’s report “Beyond the Molecular Frontier,” the following target is set for molecular simulation to become a ubiquitous tool: simulations should be achievable by a nonexpert in less than a day and for less than $1000 per simulation.156 This has already been achieved for thermochemical quantum calculations, where the cost and accuracy of a calculated heat of formation has achieved such a favorable balance when compared to experimental measurement that calculation is now the preferred method. For example, it has been estimated that in 1996 a single heat of formation could be measured calorimetrically for about $70,000, but it could be computed with comparable accuracy for only $20,000.157 By 2000, the experimental cost was $100,000 or more, versus $2,000 for a calculation. Today, it is hard to imagine a company who needs heats of formation not using calculations. There are still many barriers until this type of adoption takes place for MD simulations. The first problem is that, unlike quantum chemical calculations, MD is a “dispersed” method. There are many factors that go into an MD calculation: What is the level of detail used to describe the molecule? What force field and parameter set will be used? What properties will be calculated? What methods will be used to calculate these properties? Unlike quantum calculations where there are mainly varying levels of approximation for the same basic calculation, the “input complexity” of an MD simulation makes it more difficult to set up and run a simulation. Whereas quantum calculations are mainly concerned with obtaining minimum energy structures and thermochemistry, MD simulations are expected to compute a huge range of properties from the condensed phase to the gas phase. Making codes and methods general enough to do this for “nonexperts” is a daunting task, and there has been no real driving force to make this happen outside the commercial software companies. The open source or nominal cost software tools available are getting to the point where ease of use is not the main problem, however, but instead the problem is expertise of use. Experts are still required to set up calculations, to make sure they are reliably used, and to know what questions can and cannot be answered. It should be possible to develop standard MD methods and algorithms that can be used by nonexperts to perform routine calculations for

3075

some of the properties described in this article. For other properties, it is likely that experts will still be required to develop and apply the methods for some time to come. The future of MD is bright. As illustrated in Figure 1, the usage base of MD is growing rapidly, and there are no signs of this slowing down. Since its rather humble beginnings in the 1950s, MD has grown up to become a major component of the modern scientific and engineering research enterprise. Chemical engineers have played a significant role in getting MD to this point, and they will continue to be major users and developers of MD long into the future. Acknowledgment E.J.M. acknowledges support from the National Science Foundation (CBET 0651726) and the Air Force Office of Scientific Research (FA9550-07-1-0443). J.R.E. acknowledges support from Chemstations, Inc. Literature Cited (1) Theodorou, D. N. Progress and Outlook in Monte Carlo Simulations. Ind. Eng. Chem. Res. 2010, 49, doi: 10.1021/ie9019006. (2) Haile, J. M. Molecular Dynamics Simulation: Elementary Methods; Wiley: New York, 1992. (3) Detcheverry, F. A.; Pike, D. Q.; Nealey, P. F.; Muller, M.; de Pablo, J. J. Monte Carlo Simulation of Coarse Grain Polymeric Systems. Phys. ReV. Lett. 2009, 102, 197801. (4) Dauxios, T. Fermi, Pasta, Ulan, and a Mysterious Lady. Phys. Today 2008, 61, 55–87. (5) Fermi, E.; Pasta, J.; Ulam, S. Los Alamos Scientific Laboratory report LA-1940. In Collected Papers of Enrico Fermi, Segre, E., Ed.; University of Chicago Press: Chicago, 1965; Vol. 2. (6) Alder, B. J.; Wainwright, T. E. Phase Transition for a Hard Sphere System. J. Chem. Phys. 1957, 27, 1208. (7) Michael, G. An Interview with Berni Alder. http://www.computerhistory.info/Page1.dir/pages/Alder.html. (8) Leibson,S.StanFrankel,UnrecognizedGenius.http://www.hp9825.com/ html/stan_frankel.html (accessed March 10). (9) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculation by Fast Computing Machine. J. Chem. Phys. 1953, 21, 1087–1092. (10) Rapaport, D. C. Molecular Dynamics Simulation of Polymer Chains With Excluded Volume. J. Phys. A, Math. Gen. 1978, 11, L213–L217. (11) Denlinger, M. A.; Hall, C. K. Molecular Dynamics Simulation Results for the Pressure of Hard-Chain Fluids. Mol. Phys. 1990, 71, 541. (12) Chapela, G. A.; Martinez-Casas, S. E.; Alejandre, J. Molecular Dynamics for Discontinuous Potentials. I. General Method and Simulation of Hard Polyatomic Molecules. Mol. Phys. 1984, 53, 139. (13) Allen, P.; Tildesley, D. J. Computer Simulation of Liquid; Oxford Science Publication: Oxford, U. K., 1987. (14) Frenkel, D.; Smit, B. Understanding Molecular Simulation: from Algorithms to Applications; Academic Press: San Diego, 1996. (15) Giancoli, D. C. Physics for Scientists and Engineers, 3rd ed.; Prentice-Hall: Englewood Cliffs, N.J., 2000. (16) Kofke, D. A.; Mihalick, B.; Schultz, A. Etomica As A Java API And Development Environment For Construction And Implementation Of Molecular Simulations. www.etomica.org (accessed August 18, 2009). (17) Rapaport, D. C. Molecular Dynamics Study of a Polymer Chain in Solution. J. Chem. Phys. 1979, 71, 3299. (18) Rapaport, D. C. The Event Scheduling Problem in Molecular Dynamic Simulation. J. Comput. Phys. 1980, 34, 184. (19) Smith, S. W.; Hall, C. K.; Freeman, B. D. Molecular Dynamics Study of Transport Coefficients for Hard-chain Fluids. J. Chem. Phys. 1995, 102, 1057. (20) Smith, S. W.; Hall, C. K.; Freeman, B. D. Molecular Dynamics for Polymeric Fluids Using Discontinuous Potentials. J. Comp. Phys. 1997, 134, 16–30. (21) Chapela, G. A.; Scriven, L. E.; Davis, H. T. Molecular Dynamics for Discontinuous Potential. IV. Lennard-Jonesium. J. Chem. Phys. 1989, 91 (7), 4307. (22) Cui, J.; Elliott, J. R. Phase Envelopes for Variable Well-Width Square Well Chain Fluids. J. Chem. Phys. 2001, 114, 7283.

3076


(23) Unlu, O.; Gray, N. H.; Gerek, Z. N.; Elliott, J. R. Transferable Step Potentials for the Straight Chain Alkanes, Alkenes, Alkynes, Ethers, and Alcohols. Ind. Eng. Chem. Res. 2004, 43, 1788–1793. (24) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria.1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569–2577. (25) Ucyigitler, S.; Camurdan, M. C.; Turkay, M.; Elliott, J. R. Inferring Transferable Intermolecular Potential Models. Mol. Simul. 2008, 34, 147– 154. (26) Liu, J.-X.; Elliott, J. R. Screening Effects on Hydrogen Bonding in Chain Molecular Fluids: Thermodynamics and Kinetics. Ind. Eng. Chem. Res. 1996, 35, 2369–2377. (27) Gerek, Z. N.; Elliott, J. R. Self-Diffusivity Estimation By Molecular Dynamics. Ind. Eng. Chem. Res., in press. (28) Gibson, J. B.; Goland, A. N.; Milgram, M.; Vineyard, G. H. Dynamics of Radiation Damage. Phys. ReV. 1960, 120, 1229–1253. (29) Rahman, A. Correlation in the Motion of Atoms in Liquid Argon. Phys. ReV. 1964, 136A, 405. (30) Verlet, L. Computer “Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-jones Molecules. Phys. ReV. 1967, 159, 98. (31) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. J. Chem. Phys. 1982, 76, 637. (32) Hansen, J.-P.; Verlet, L. Phase Transitions of the Lennard-Jones System. Phys. ReV. 1969, 184, 151–161. (33) Woodcock, L. V. Isothermal Molecular Dynamics Calculations for Liquid Salts. Chem. Phys. Lett. 1971, 10, 257–261. (34) Anderson, H. C. Molecular Dynamics Simulation at Constant Pressure and/or Temperature. J. Chem. Phys. 1980, 72, 2384–2393. (35) Haile, J. M.; Graben, H. W. On the isoenthalpic-isobaric Ensemble in Classical Statistical Mechanics. Mol. Phys. 1980, 40, 1433–1439. (36) Nose, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255–268. (37) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. ReV. A 1985, 31, 1695–1697. (38) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Simulation of Electrostatic Systems in Periodic Boundary Condition. 1. Lattice Sums and Dielectric Constants. Proc. R. Soc. A 1980, 372, 27–56. (39) Darden, D.; York, D.; Pedersen, L. Particle Mesh EwaldsAn N log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089–10092. (40) Rahman, A.; Stillinger, F. E. Molecular Dynamics Study of Liquid Water. J. Chem. Phys. 1971, 55, 3336–3359. (41) Ryckaert, J.-P.; Bellemans, A. Molecular Dynamics of Liquid n-butane Near Its Boiling Point. Chem. Phys. Lett. 1975, 30 (1), 123–125. (42) Murad, S.; Evans, D. J.; Gubbins, K. E.; Streett, W. B.; Tildesley, D. J. Molecular Dynamics Simulation of Dense Fluid Methane. Mol. Phys. 1979, 37, 725–736. (43) Evans, D. J. On the Representation of Orientation Space. Mol. Phys. 1977, 34, 317–325. (44) Evans, D. J.; Murad, S. Singularity Free Algorithm for Molecular Dynamics Simulation of Rigid Polyatomics. Mol. Phys. 1977, 34, 327– 331. (45) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comp. Phys. 1977, 23, 327–341. (46) Postma, J. P. M.; Berendsen, H. J. C.; Haak, J. R. Thermodynamics of Cavity Formation in WatersA Molecular Dynamics Study. Faraday Symp. Chem. Soc. 1982, 17, 55–67. (47) Singh, U. C.; Brown, F. K.; Bash, P. A.; Kollman, P. A. An Approach to the Application of Free Energy Perturbation Methods Using Molecular Dynamics - Applications to the Transformations of CH3OH f CH3CH3, H3O+ f NH4+, Glycine f Alanine and Alanine f Phenylalanine in Aqueous Solution and to H3O+(H2O) 3 f NH4+(H2O)3 in the Gas Phase. J. Am. Chem. Soc. 1987, 109, 1607–1614. (48) Rigby, D.; Roe, R. J. Molecular-Dynamics Simulation Of Polymer Liquid And Glass 0.1. Glass-Transition. J. Chem. Phys. 1987, 87, 7285– 7292. (49) Kremer, K.; Grest, G. S. Dynamics of Entangled Linear Polymer Melts: A Molecular Dynamics Simulation. J. Chem. Phys. 1990, 92, 5057. (50) Mansfield, K. F.; Theodorou, D. N. Molecular Dynamics Simulation of a Glassy Polymer Surface. Macromolecules 1991, 24, 6283–6294. (51) Plimpton, S. J. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys. 1995, 117, 1–19.

(52) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Klaus Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802. (53) Smith, W.; Forester, T. R. DL_POLY_2.0: a general-purpose parallel molecular dynamics simulation package. J. Mol. Graphics 1996, 14, 136–141. (54) Ponder, J. W. TINKER Home Page. http://dasher.wustl.edu/tinker/ (accessed November 27). (55) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187–217. (56) Pearlman, D. A.; Case, D. A.; Caldwell, J. W.; Ross, W. S.; Cheatham, T. E., III.; DeBolt, S.; Ferguson, D.; G., S.; Kollman, P. AMBER, a Package of Computer Programs for Applying Molecular Mechanics, Normal Mode Analysis, Molecular Dynamics and Free Energy Calculations to Simulate the Structural and Energetic Properties of Molecules. Comput. Phys. Commun. 1995, 91, 1–41. (57) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: A Package for Molecular Simulation and Trajectory Analysis. J. Mol. Model. 2001, 7, 306–317. (58) Accelrys Accelrys Software, Inc. http://accelrys.com/ (accessed November 27). (59) HyperCube. HyperChem Professional Overview. http://www.hyper. com (accessed November 27). (60) Rigby, D. Practical R&D Applications: Current, Intermediate and Long Term Perspectives. In Proceedings of the AIChE Annual Meeting, Philadelphia, PA, 2008. (61) Ashurst, W. T.; Hoover, W. G. Dense Fluid Shear Viscosity Via Nonequilibrium Molecular Dynamics. Phys. ReV. A 1975, 11, 658–678. (62) Evans, D. J.; Morris, G. P. Non-Newtonian Molecular Dynamics. Comp. Phys. Rep. 1984, 1, 297–343. (63) Simmons, A. D.; Cummings, P. T. Nonequilibrium Molecular Dynamics Simulation of Dense Fluid Methane. Chem. Phys. Lett. 1986, 129, 92–98. (64) Cummings, P. T.; Evans, D. J. Nonequilibrium Molecular Dynamics Approaches to Transort Properties and Non-Newtonian Fluid Rehology. Ind. Eng. Chem. Res. 1992, 31, 1237–1252. (65) Edberg, R.; Evans, D. J.; Morris, G. P. Rheoogy of n-Alkanes by Nonequilibrium Molecular Dynamics. J. Chem. Phys. 1987, 86, 4555–4570. (66) Rowley, R. L.; Ely, J. F. Nonequilibrium Molecular Dynamics Dimulations of Structured Molecules. Part I. Isomeric Effects on the Viscosity of Butanes. Mol. Phys. 1991, 72 (4), 831–846. (67) Cummings, P. T.; Varner, T. L. Nonequilibrium Molecular Dynamics Calculation of the Shear Viscosity of Liquid Water. J. Chem. Phys. 1988, 89, 6391–6398. (68) Guevara-Carrion, G.; Nieto-Draghi, C.; Vrabec, J.; Hasse, H. Prediction of Transport Properties by Molecular Simulation: Methanol and Ethanol and their Mixture. J. Phys. Chem. B 2008, 112, 16664–16674. (69) Desgranges, C.; Delhomelle, J. Rheology of Liquids fcc Metals: Equilibrium and Transient-Time Correlaton Function Nonequilibrium Molecular Dynamics Simulations. Phys. ReV. B 2008, 78, 184202. (70) Maginn, E. J.; Bell, A. T.; Theodorou, D. N. Transport Diffusivity of Methane in Silicalite From Equilibrium and Nonequilibrium Simulations. J. Phys. Chem. 1993, 97, 4173. (71) Ahunbay, M. G.; Elliott, J. R.; Talu, O. Effect of Surface Resistances on the Diffusion of Binary Mixtures in the Silicalite Single Crystal Membrane. J. Phys. Chem. B 2005, 109, 923–929. (72) Kioupis, L. I.; Maginn, E. J. Molecular Simulation of Poly-AlphaOlefin Synthetic Lubricants: Impact of Molecular Architecture on Performance Properties. J. Chem. Phys. B 1999, 103, 10781–10790. (73) Moore, J. D.; Cui, S. T.; Cochran, H. T.; Cummings, P. T. Rheology Of Lubricant Basestocks: A Molecular Dynamics Study Of C30 Isomers. J. Chem. Phys. 2000, 113, 8833–8840. (74) Yang, Y.; Pakkanen, T. A.; Rowley, R. L. NEMD Simulations of Viscosity and Viscosity Index for Lubricant-Size Model Molecules. Int. J. Thermophys. 2002, 23, 1441–1454. (75) Gordon, P. A. Development of Intermolecular Potentials for Predicting Transport Properties of Hydrocarbons. J. Chem. Phys. 2006, 125, 014504. (76) Kelkar, M. S.; Shi, W.; Maginn, E. J. Determining the Accuracy of Classical Force Fields for Ionic Liquids: Atomistic Simulation of the Thermodynamic and Transport Properties of 1 Ethyl-3-methylimidazolium Ethylsulfate ([emim][EtSO4]) and Its Mixtures with Water. Ind. Eng. Chem. Res. 2008, 47, 9115–9126. (77) Calambra, N.; Nieto-de Castro, C. A.; Ely, J. F. Shear Viscosity of Molten Alkali Halides from Equilibrium and Nonequilibrium Molecular Dynamics Simulations. J. Chem. Phys. 2005, 122, 224501.

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010 (78) Hulse, R. J.; Ely, J. F.; Wilding, N. V. Transient Nonequilibrium Molecular Dynamic Simulations of Thermal Conductivity. 1. Simple Fluids. Int. J. Thermophys. 2005, 26, 112. (79) Kelkar, M. S.; Maginn, E. J. Rapid Shear Viscosity Calculation by Momentum Impulse Relaxation Molecular Dynamics. J. Chem. Phys. 2005, 123, 224904. (80) Heffelfinger, G. S.; van Swol, F. Diffusion in Lennard-Jones Fluids using Dual Control Volume Grand Canonical Dynamics Simulation(DCVGCMD). J. Chem. Phys. 1994, 100 (10), 7548–7552. (81) Heffelfinger, G. S.; Ford, D. M. Massively Parallel Dual Control Volume Grand Canonical Molecular Dynamics with LADERA I. Gradient Driven Diffusion in Lennard-Jones Fluids. Mol. Phys. 1998, 94, 659–671. (82) Arya, G.; Chang, H.-C.; Maginn, E. J. A Critical Comparison of Equilibrium, Nonequilibrium and Boundary-Driven Molecular Dynamics Techniques for Studying Transport in Microporous Materials. J. Chem. Phys. 2001, 115, 8112–8124. (83) Weeks, J. D.; Chandler, D.; Anderson, H. C. Role of Repulsive Forces in Determining the Equilibrium Structure of Simple Liquids. J. Chem. Phys. 1971, 54, 5237. (84) Peter, C.; Delle Site, L.; Kremer, K. Multiscale simulation of soft matter: From scale bridging to adaptive resolution. Annu. ReV. Phys. Chem. 2008, 59, 545–571. (85) Alley, W. E.; Alder, B. J. Studies in Molecular Dynamics XV. High Temperature Description of the Transport Coefficients. J. Chem. Phys. 1975, 63, 3764–3768. (86) Nguyen, H. D.; Hall, C. K. Molecular dynamics simulations of spontaneous fibril formation by random-coil peptides. PNASsUSA 2004, 101, 16180–16185. (87) Haile, J. M.; O’Connell, J. P. Internal Structure Of A Model Micelle Via Computer-Simulation. J. Phys. Chem. 1984, 88, 6363–6366. (88) Nelson, P. H.; Hatton, T. A.; Rutledge, G. C. Asymmetric growth in micelles containing oil. J. Chem. Phys. 1999, 110, 9673–9680. (89) Abu-Sharkh, B. F.; Hamad, E. Z. Investigation of the microstructure of micelles formed by hard-sphere chains interacting via size nonadditivity by discontinuous molecular dynamics simulation. Langmuir 2004, 20, 254– 259. (90) Reddy, G.; Yethiraj, A. Implicit and Explicit Solvent Models for the Simulation of Dilute Polymer Solutions. Macromolecules 2006, 39, 8536–8542. (91) Dijkstra, M.; R., v. R.; Evans, R. Phase diagram of highly asymmetric binary hard-sphere mixtures. Phys. ReV. E. 1999, 59, 5744– 5771. (92) Noid, W. G.; Liu, P.; Wang, Y.; Chu, J. W.; Ayton, G. S.; Izvekov, S.; Anderson, H. C.; Voth, G. A. The multiscale coarse-graining method. II. Numerical implementation for coarse-grained molecular models. J. Chem. Phys. 2008, 128, 244115. (93) Silbermann, J. R.; Schoen, M.; Klapp, S. H. L. Coarse-grained single-particle dynamics in two-dimensional solids and liquids. Phys. ReV. E. 2008, 78, 011201. (94) Gubbins, K. E.; Moore, J. D. Molecular Modeling of Matter: Impact and Prospects in Engineering. Ind. Eng. Chem. Res. 2010, 49, doi: 10.1021/ ie901909c. (95) Malevanets, A.; Kapral, R. Solute Molecular Dynamics In A Mesoscale Solvent. J. Chem. Phys. 2000, 112, 7260. (96) Padding, J. T.; Louis, A. A. Interplay between hydrodynamic and Brownian fluctuations in sedimenting colloidal suspensions. Phys. ReV. E. 2008, 77, 011402. (97) Donev, A.; Garcia, A. L.; Alder, B. J. Stochastic Event-Driven Molecular Dynamics. J. Comp. Phys. 2008, 227, 2644–2665. (98) Bird, G. A. Recent Advances and Current Challenges for DSMC. Computers Math. App. 1998, 35 (1/2), 1–14. (99) Virnau, P.; Binder, K.; Heinz, H.; Kreer, T.; Muller, M. M. Molecular Simulation and Thermodynamics of Polymer Blends and Melts. In Encyclopedia of Polymer Blends, Isayev, A. I., Ed.; Wiley-VCH: New York, 2009; Vol. 1. Fundamentals. (100) Shaw, D. E., Millisecond-Long Molecular Dynamics Simulations of Proteins on the Anton Machine. In Proceedings of FOMMS 2009: Foundations for InnoVation, Blaine, WA, 2009. (101) June, R. L.; Bell, A. T.; Theodorou, D. N. Transition State Studies of Xenon and SF6 Diffusion in Silicalite. J. Phys. Chem. 1991, 95, 8866– 8878. (102) Tunca, C.; Ford, D. M. A Transition State Theory Approach to Adsorbate Dynamics at Arbitrary Loadings. J. Chem. Phys. 1999, 111, 2751–2760. (103) Dubbledam, D.; Beerdsen, E.; Calero, S.; Smit, B. Dynamically Corrected Transition State Theory Calculations of Self-Diffusion in Anisotropic Nanoporous Materials. J. Phys. Chem. B 2006, 110, 3164–3172.

3077

(104) Maginn, E. J.; Bell, A. T.; Theodorou, D. N. Dynamics of Long n-Alkanes in Silicalite: A Hierarchical Simulation Approach. J. Phys. Chem. 1996, 100, 7155–7173. (105) Fichthorn, K. A.; Miron, R. A. Thermal Desorption of Large Molecules from Solid Surfaces. Phys. ReV. Lett. 2002, 89, 196103. (106) Santiso, E. E.; Gubbins, K. E. Multi-Scale Molecular Modeling of Chemical Reactivity”, Molecular Simulation. Mol. Simul. 2004, 30, 699– 748. (107) Becker, K. E.; Fichthorn, K. A. Accelerated Molecular Dynamics Simulation of the Thermal Desorption of n-Alkanes from the Basal Plane of Graphite. J. Chem. Phys. 2006, 125, 184706. (108) Tait, S. L.; Dohna`lek, Z.; Campbell, C. T.; Kay, B. D. Desorption Energies of n-Alkanes. 2009, unpublished. (109) Paserba, K. R.; Gellman, A. J. Effects of Conformational Isomerism on the Desorption Kinetics of n-Alkanes from Graphite. J. Chem. Phys. 2001, 115, 6737–6751. (110) Zhang, L.; Greenfield, M. L. Relaxation Time, Diffusion and Viscosity Analysis of Model Asphalt Systems Using Molecular Simulation. J. Chem. Phys. 2007, 127, 194502. (111) Magda, J. J.; Tirrell, M.; Davis, H. T. Molecular Dynamics of Narrow, Liquid-Filled Pores. J. Chem. Phys. 1986, 83, 1888–1901. (112) Maginn, E. J.; Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. Simulation of Hydrocarbon Diffusion in Zeolites. Stud. Surf. Sci. Catal. 1997, 105, 1851–1858. (113) Travis, K. P.; Gubbins, K. E. Combined Diffusive and Viscous Transport in a Carbon Slit Pore. Mol. Simul. 2000, 25, 209–227. (114) Snurr, R. Q.; Karger, J. Molecular Simulatons and NMR Measurements of Binary Diffusion in Zeoites. J. Phys. Chem. B 1997, 101, 6469– 6473. (115) Papadopoulos, G. K.; Jobic, H.; Theodorou, D. N. Transport Diffusivity of N2 and CO2 in Silicalite: Coherent Quasielastic Neutron Scattering Measurements and Molecular Dynamics Simulations. J. Phys. Chem. B 2004, 108, 12748–12756. (116) Karger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; Wiley: New York, 1992. (117) Sholl, D. S. Understanding Macroscopic Diffusion of Adsorbed Molecules in Crystaline Nanoporous Materials via Atomistic Simulations. Acc. Chem. Res. 2006, 39, 403–411. (118) Chen, H.; Sholl, D. S. Rapid Diffusion of CH4/H2 Mixtures in Single-Walled Nanotubes. J. Am. Chem. Soc. 2004, 126, 7778–7779. (119) Schuring, A.; Vasenkov, S.; Fritzsche, S. Influence of Boundaries of Nanoporous Crystals on Moelcular Exchange Under the Conditions of Single-File Diffusion. J. Phys. Chem. B 2005, 109, 16711–16717. (120) Huggins, M. L.; Mayer, J. E. Interatomic Distances in Crystals of the Alkali Halides. J. Chem. Phys. 1933, 1, 643. (121) Mayer, J. E. Dispersion and Polarizability and the van der Waals Potential in the Alkali Halides. J. Chem. Phys. 1933, 1, 270. (122) Tosi, M. P.; Fumi, F. G. Ionic Sizes and Born Repulsive Parameters in the NaCl-type Alkali Halides. I. J. Phys. Chem. Solids 1964, 25, 31. (123) Sangster, M. J.; Dixon, M. Interionic Potentials in Alkali Halides and Their Use in Simulations of the Molten Salts. AdV. Phys. 1976, 25, 247. (124) Galambra, N.; Nieto-de Castro, C. A.; Ely, J. F. Shear Viscosity of Molten Alkali Halides from Equilibrium and Nonequilibrium Molecular Dynamics Simulations. J. Chem. Phys. 2005, 122, 224501. (125) Hanke, C. G.; Price, S. L.; Lynden-Bell, R. M. Intermolecular Potentials for Simulations of Liquid Imidazolium Salts. Mol. Phys. 2001, 99, 801. (126) Margulis, C. J.; Stern, H. A.; Berne, B. J. Computer Simulation of a “Green Chemistry” Room-Temperature Ionic Solvent. J. Phys. Chem. B 2002, 106, 12017. (127) Morrow, T. I.; Maginn, E. J. Molecular Dynamics Study of the Ionic Liquid 1-n-Butyl-3-methylimidazolium Hexafluorophosphate. J. Phys. Chem. B 2003, 106, 12807. (128) Maginn, E. J., Atomistic Simulation of Ionic Liquids. In ReViews in Computational Chemistry; Lipkowitz, K. B., Larter, R., Cundari, T., Eds.; Wiley: Hoboken, NJ, 2009; Vol 24. (129) MacKerell, A. D., Jr. Empirical Force Fields for Biological. Macromol.: OVerView Issues J. Comp. Chem. 2004, 25, 1584–1604. (130) McQuaid, M. I.; Sun, H.; Rigby, D. Development and Validation of COMPASS Force Field Parameters for Molecules with Aliphatic Azide Chains. J. Comput. Chem. 2004, 25, 61–71. (131) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 6. United-Atom Description for Ethers, Glycols, Ketones, and Aldehydes. J. Phys. Chem. B 2004, 108, 17596–17605. (132) Boutard, Y.; Ungerer, P.; Teuler, J. M.; Ahunbay, M. G.; Sabater, S. F.; P’erez-Pellitero, J.; Mackie, A. D.; Bourasseau, E. Extension of the

3078


anisotropic united atoms intermolecular potential to amines, amides and alkanols. Fluid Phase Equilib. 2005, 236, 25–41. (133) Chen, B.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 3. Explicit-Hydrogen Description of Normal Alkanes. J. Phys. Chem. B 1999, 103, 5370–5379. (134) Emami, F. S.; Vahid, A.; Elliott, J. R. Finitely limited group contribution correlations for boiling temperatures. J. Chem. Thermodyn. 2009, 41, 530–537. (135) Emami, F. S.; Vahid, A.; Elliott, J. R.; Feyzi, F. Group Contribution Prediction of Vapor Pressure with Statistical Associating Fluid Theory, Perturbed-Chain Statistical Associating Fluid Theory, and Elliott-SureshDonohue Equations of State. Ind. Eng. Chem. Res. 2008, 47, 8401–8411. (136) Vahid, A.; Sans, A. D.; Elliott, J. R. Correlation of Mixture VLE with the SPEADMD model. Ind. Eng. Chem. Res. 2008, 47, 7955–7964. (137) Liu, H.; Silva, C. M.; Macedo, E. A. Unified approach to the selfdiffusion coefficients of dense fluids over wide ranges of temperature and pressure: hard-sphere, square-well, Lennard-Jones and real substances. Chem. Eng. Sci. 1998, 53, 2403. (138) Berker, A.; Chynoweth, S.; Klomp, U. C.; Michopoulos, Y. Nonequilibrium Molecular Dynamics (NEMD) Simulations and the Rheological Properties of Liquid Normal Hexadecane. J. Chem. Soc., Faraday Trans. 1992, 88, 1719–1725. (139) Cui, S. T.; Gupta, S. A.; Cummings, P. T.; Cochran, H. D. Molecular Dynamics Simulation of the Rheology of Normal Decane, Hexadecane, and Tetracosane. J. Chem. Phys. 1996, 105, 1214–1220. (140) Allen, W.; Rowley, R. L. Predicting the Viscosity of Alkanes Using Nonequilibrium Molecular Dynamics: Evaluation of Intermolecular Potential Models. J. Chem. Phys. 1997, 106, 10273–10281. (141) Dysthe, D. K.; Fuchs, A. H.; Rousseau, B. Fluid Transport Properties by Equilibrium Molecular Dynamics. I. Methodology at Extreme Fluid States. J. Chem. Phys. 1999, 110 (8), 4047–4059. (142) Kioupis, L. I.; Maginn, E. J. Impact of Molecular Architecture on the High Pressure Rheology of Hydrocarbon Fluids. J. Phys. Chem. B 2000, 104, 7774–7783. (143) McCabe, C.; Cui, S. T.; Cummings, P. T.; Gordon, P. A.; Saeger, R. B. Examining the Rheology of 9-Octylheptadecane to Giga-Pascal Pressures. J. Chem. Phys. 2001, 114, 1887–1891. (144) Wheeler, D. R.; Fuller, N. G.; Rowley, R. L. Nonequilibrium Molecular Dynamics Simulations of the Shear Viscosity of Liquid Methanol: Adaptation of the Ewald Sum to the Lees-Edwards Boundary Condition. Mol. Phys. 1997, 92, 55–62. (145) Muller-Plathe. Reversing the Perturbation in Nonequilibrium Molecular Dynamics: An Easy Way to Calculate the Shear Viscosity of Fluids. Phys. ReV. E. 1999, 59, 4894–4898. (146) Arya, G.; Maginn, E. J.; Chang, H.-C. Efficient Viscosity Estimation from Molecular Dynamics Simulation Via Momentum Impulse Relaxation. J. Chem. Phys. 2000, 113, 2079–2087. (147) Boyd, R. H. Glass Transition Temperatures from Molecular Dynamics Simulations. Trends Polym. Sci. 1996, (4).

(148) Antoniadis, S. J.; Samara, C. T.; Theodorou, D. N. Molecular dynamics of atactic polypropylene melts. Macromolecules 1998, 31, 7944– 7952. (149) van Krevelen, D. W. Properties of Polymers: Their Correlation with Chemical Structure, Their Numerical Estimation and Prediction from AdditiVe Group Contributions, 3rd ed.; Elsevier: Amsterdam, 1990. (150) Gusev, A. A.; Muller-Plathe, F.; van Gunsteren, W. R.; Suter, U. W. Dynamics of small molecules in bulk polymers. AdV. Polym. Sci. 1994, 116, 207–247. (151) Yampolskii, Y.; Pinnau, I.; Freeman, B. D. Materials Science of Membranes for Gas and Vapor Separation; John Wiley: Hoboken, NJ, 2006. (152) Case, F.; Chaka, A. M.; Friend, D. G.; Frurip, D.; Golab, J.; Johnson, R.; Moore, J. D.; Mountain, R. D.; Olson, J.; Schiller, M.; Storer, J. The First Industrial Fluid Properties Simulation Challenge. Fluid Phase Equilib. 2004, 217, 1–10. (153) Case, F. H.; Brennan, J.; Chaka, A.; Dobbs, K. D.; Friend, D. G.; Gordon, P. A.; Moore, J. D.; Mountain, R. D.; Olson, J. D.; Ross, R. B.; Schiller, M.; Shen, V. K.; Stahlberg, E. A. The Fourth Industrial Fluid Properties Simulation Challenge. Fluid Phase Equilib. 2008, 274, 2–9. (154) Case, F. H.; Brennan, J.; Chaka, A.; Dobbs, K. D.; Friend, D. G.; Frurip, D.; Gordon, P. A.; Moore, J. D.; Mountain, R. D.; Olson, J. D.; Ross, R. B.; Schiller, M.; Shen, V. K.; Stahlberg, E. A. The Third Industrial Fluid Properties Simulation Challenge. Fluid Phase Equilib. 2007, 260, 153– 163. (155) Gani, R. Chemical product design: challenges and opportunities. Comput. Chem. Eng. 2004, 28, 2441–2457. (156) Breslow, R.; Tirrell, M. V. Beyond the Molecular Frontier Challenges for Chemistry and Chemical Engineering; The National Academies Press: Washington, D.C., 2003; p 224. (157) Westmoreland, P. R. Kollman, P. A. Chaka, A. M. Cummings, P. T. Morokuma, K. Neurock, M. Stechel, E. B. Vashishta, P. Applying of Molecular and Materials Modeling; Kluwer: Norwell, MA, 2002. (158) Tait, S. L.; Dohna`lek, Z.; Campbell, C. T.; Kay, B. D., unpublished. (159) Paserba, K. R.; Gellman, A. J. Effects of Conformational Isomerism on the Desorption Kinetics of n-Alkanes from Graphite. J. Chem. Phys. 2001, 115, 6737–6751.

E. J. Maginn† and J. R. Elliott*,‡ Department of Chemical and Biomolecular Engineering, UniVersity of Notre Dame, Notre Dame, Indiana 46556, Department of Chemical and Biomolecular Engineering, UniVersity of Akron, Akron, Ohio 44325-3906 ReceiVed for reView December 2, 2009 IE901898K

Historical Perspective and Current Outlook for Molecular Dynamics As

Recommend Documents