In the Laboratory
Molecular Dynamics Simulations of Simple Liquids
W
Owen F. Speer, Brian C. Wengerter, and Ramona S. Taylor* Department of Chemistry, College of the Holy Cross, Worcester, MA 01610; *
[email protected] Computer simulations are now used to study complex chemical events ranging from the folding of proteins (1) to the high energy desorption of organic molecules from metal surfaces (2). Simulations not only offer an alternative method for performing costly experiments, but they also provide an added molecular-level insight into the chemistry that is not always available via experiment. In the past, computational techniques were introduced to the undergraduate student as an aid in the visualization of molecular structures or as a method of determining reaction energetics. These applications failed to demonstrate the power of the computational methods. Students were left unaware of which computational method is most appropriate for the problem at hand or of the abundance of problems that could be examined via computational techniques. Even worse, the theory behind these techniques was often deemed too demanding for student comprehension, and left unmentioned. Recently, however, many examples of how the theories behind computer simulations can be showcased in the undergraduate classroom have been presented in this Journal (3–9). Students are introduced to the theories behind these computational techniques in the physical chemistry lecture. Consequently, the physical chemistry laboratory is a natural place for students to use these techniques. This laboratory experiment utilizes the technique of molecular dynamics (MD) computer simulation to investigate the microscopic structure of simple liquids. In the spirit of the inquiry-based approach, students are encouraged to explore how varying simulation parameters affects the structure and physical properties of a specific liquid. Students utilize the Amber suite (10) of programs to simulate the behavior of one or more liquids under varying environmental conditions. Amber is commonly utilized in both academia and industry to simulate condensed phase systems involving small molecules or large biomolecules. Although its interface is not graphical, our students have quickly mastered its use. The advantage of Amber over more graphical-based programs is that many statistical mechanical analysis programs exist within the Amber suite. These analysis programs facilitate the calculation of a variety of statistical mechanical properties, such as radial distribution functions, from the simulated data. Liquids are particularly interesting because, even though they do not exhibit the macroscopic order of solids, they possess a definite molecular-level structure. The long-range organizational structure that emerges from this molecular-level structure is not readily apparent with even today’s most sophisticated molecular visualization packages. However, it can be determined from a liquid’s atom–atom radial distribution function. The atom–atom radial distribution function, gαβ(r) is defined as (11–13), ρ (r ) gαβ ( r ) = β (1) ρ bulk β
1330
Journal of Chemical Education
•
where ρβ(r) is the density of β-type atoms, averaged over the length of the simulation, at a distance r from a central αtype atom and ρβbulk is the average density for β-type atoms in the bulk of the liquid. Thus, gαβ(r) is a measure of the local structure of the β-type atoms in the liquid as you move radially away from the central α-type atom. For a random distribution of molecules, gαβ(r) would equal one. In addition, the individual peaks in the radial distribution function can be integrated to determine the coordination number, or the average number of nearest neighbors within each neighbor shell, nc (11, 13): nc = 4πρβbulk ∫ gαβ ( r ) r 2 dr
(2)
For a molecular liquid, such as water, there are three different atomic radial distribution functions that define the structure of the liquid: gOO(r), gHO(r), and gHH(r). Each function highlights a different aspect of the overall structure of the water and the functions taken together form a complete picture of the local environment experienced by each water molecule. Discussions of the radial distribution function and its relationship to the long-range organizational structure of a liquid can be found in many modern physical chemistry textbooks (14–16). Most physical chemistry students realize that liquids have short-range organization—that is, more molecular-level structure than gases and less than solid crystals—but they have little idea of how to make this observation quantitative. Radial distribution functions provide them with visual proof that what they have been taught is, in fact, the case. Many students are excited to merely explore how changing the temperature affects the radial distribution function. Other students wonder how the different macroscopic properties of two liquids will be reflected in their microscopic structure. Finally, a few students become quite ambitious and attempt to observe specific types of intermolecular structure such as the appearance of the π–π stacking in benzene or the dipole– dipole interactions in acetonitrile. Procedure This laboratory has been performed successfully in our junior-level physical chemistry laboratory course. As preparation, the Lennard-Jones potential is introduced in an earlier physical chemistry lecture. Both the elementary basis of classical dynamics simulations and how one would actually perform a MD simulation are then discussed in the prelaboratory lecture. A thorough overview of this material is provided in the Supplemental Material.W This experiment is conducted over two, four-hour laboratory periods. In the first week, the students are taught how to use Amber (10) to perform a MD simulation for a bulk liquid, such as water, and how to view snapshots of their
Vol. 81 No. 9 September 2004
•
www.JCE.DivCHED.org
In the Laboratory
www.JCE.DivCHED.org
•
6
5
4
gOO(r)
simulations (17). At the end of the first week, the students are asked to design a set of MD simulations to investigate some aspect of how modifying the simulation parameters influences the microscopic structures of their liquids. Students may choose either to model two very different liquids, such as H2O and benzene, or to examine the effect varying an environmental parameter, such as temperature, has on the structure of a liquid. Along with deciding on which chemical question to explore, each group of students is asked to formulate a viable hypothesis for what they expect to happen and to construct a tentative schedule of when the simulations will be performed and who will be responsible for starting the various simulation runs so that all of the data can be collected and analyzed before the end of week two. A typical group consists of two to three students. Each group has one SGI Octane 2 workstation at its disposal for the two-week period allocated to do this experiment. Thus, this experiment would require a two-week commitment of four workstations for a laboratory section of 12 students. Alternatively, a single workstation can be used if this experiment is but one of the experiments that the students rotate through during the course of the semester. This latter roundrobin format has worked well in our physical chemistry laboratory course at Holy Cross. Students are given an initial set of equilibrated coordinates (at 25 ⬚C) and a topology file for each of the liquids available for analysis. Each topology file contains the bonding and potential model information for its specific liquid. The liquids available for analysis include water (756 molecules) (18), benzene (288 molecules) (19), ethanol (267 molecules) (20), and acetonitrile (216 molecules) (21). Students perform an equilibration run of 200 ps and a production run of 50 ps, which is adequate to calculate statistically meaningful radial distribution functions for each of their proposed systems. At 25 ⬚C, the size of the equilibrated periodic boxes range from 25.0 Å × 25.0 Å × 25.0 Å for 756 water molecules to 35.0 Å × 35.0 Å × 35.0 Å for 288 benzene molecules. The two observables that students measure are the density and the average microscopic structure of the liquid. The density of a liquid depends both on the mass of the particles that comprise the liquid and the intermolecular spacing between these particles. Thus, one needs to be careful when comparing densities of dissimilar liquids. To this end, students are re-introduced to the concept of number density as opposed to mass density. In a constant pressure simulation (a NPT ensemble where N is the number of molecules and T is the temperature), density can be calculated by dividing the number of particles by the volume of the simulation box. Amber reports this as mass density in the output file (10) and the students convert it to number density. As described above, the structure of the liquid can be derived from its radial distribution function, gαβ(r). It is relatively straightforward to calculate the atomic radial distribution functions for a bulk liquid using the analysis programs that are available within the Amber 6.0 (or higher) suite of programs (10). Detailed instructions on how to calculate radial distribution functions using Amber 6.0 are included in the Supplemental Material.W
200 K 298 K 313 K 343 K
3
2
1
0 2
3
4
5
6
7
rOO / Å Figure 1. O–O radial distribution function for ethanol as a function of temperature.
Hazards No significant hazards exist. Examples of Student Results One of the most common questions posed by the students is how does temperature affect the microscopic structure of the liquid. The oxygen–oxygen radial distribution functions calculated in one such investigation are shown in Figure 1. Here, the students chose to model ethanol at four different temperatures: 200 K, 298 K, 313 K, and 343 K. At 200 K, the O⫺O distribution function is sharply peaked. As the temperature is increased, the peaks begin to broaden and lessen in intensity. Yet, the appearance of liquid-like structure is retained. At 200 K, the probability of finding an ethanol molecule between the first and second neighbor shell is approximately zero. However, upon increased temperature this changes and the likelihood of finding an ethanol molecule between the first and second coordination shells increases. Beyond 6 Å from the central ethanol molecule, no structure is apparent in the liquid. The C⫺C radial distribution function for the same four systems where the C⫺C distance is calculated from the center-of-mass between the two carbon atoms in ethanol’s tail is shown in Figure 2. Notice that the peaks are much less intense for the C⫺C interactions than for the O⫺O interactions. This signifies that the hydrogen bonding interactions are the driving force in determining the local structure in ethanol at all four temperatures. Besides this decrease in peak intensities, two additional items are worth noticing in Figure 2. First, the C⫺C interactions show a much longer-range structure than the O⫺O interactions. This occurs because the carbon tails of ethanol extend approximately 4 Å beyond the O atoms. Second, the nearest-neighbor C⫺C peak at 200 K is split into two distinct peaks. Upon careful analysis, stu-
Vol. 81 No. 9 September 2004
•
Journal of Chemical Education
1331
In the Laboratory
step simulation for each liquid; a complete laboratory writeup for use by the students; and comments on possible modifications to this laboratory are available in this issue of JCE Online.
2.0
gCC(r )
1.5
Literature Cited 1.0
200 K 298 K 313 K 343 K
0.5
0.0 2
3
4
5
6
7
8
9
10
rCC / Å Figure 2. C–C radial distribution function for ethanol as a function of temperature.
dents can show that the two C⫺C nearest neighbor geometries originate from the hydrogen-donating and hydrogenaccepting abilities of each ethanol molecule involved in a hydrogen bond. In agreement with the experimental results (22), the density of ethanol decreases from 0.88 g兾cm3 at 200 K to 0.73 g兾cm3 at 343 K. Using the calculated density values, the O⫺O radial distribution function and eq 2, one finds that the number of nearest neighbors in the first coordination shell of ethanol decreases from 2.0 at 200 K to 1.9 at 343 K. Summary Students are given the opportunity to perform molecular dynamics simulations on a series of molecular liquids using the Amber suite of programs. They are introduced both to the physical theories underlying classical mechanics simulations and to the atom–atom pair distribution function. Using this distribution function, they can graphically “see” that liquids do in fact posses some sort of long-range microscopic structure and that the extent of this structure differs from liquid to liquid and from one temperature to another. In the spirit of the guided-inquiry approach, each group of students is asked to design a set of simulations to answer their selfdetermined question. This allows the students a chance to develop their own scientific creativity. Acknowledgments We would like to thank the Camille and Henry Dreyfus Foundation (Grant # SU-99-076) and the College of Holy Cross for their support of this project. Also RST would like to thank the students in her spring 2001 and fall 2002 Physical Chemistry Laboratory courses for their help and patience in the development of this laboratory. W
Supplemental Material
Instructor’s notes, which include an equipment list, a detailed description of how to obtain radial distribution functions using Amber 6, and CPU times for completing a 50,000 1332
Journal of Chemical Education
•
1. Tiraldo-Rives, J.; Jorgensen, W. L. Biochemistry 1993, 32, 4175. 2. Garrison, B. J.; Delcorte, A.; Krantzman, K. D. Acc. Chem. Res. 2000, 33, 69. 3. Gasyna, Z.; Rice, S. A. J. Chem. Educ. 1999, 76, 1023. 4. Paselk, R. A.; Zoellner, R. W. J. Chem. Educ. 2002, 79, 1192. 5. Whisnant, D. M.; Howe, J. J.; Lever, L. S. J. Chem. Educ. 2000, 77, 199. 6. Martin, N. H. J. Chem. Educ. 1998, 75, 241. 7. Taylor, A. T. S.; Feller, S. E. J. Chem. Educ. 2002, 79, 1467. 8. Saiz, E.; Tarazona, M. P. J. Chem. Educ. 1997, 74, 1350. 9. Cramer, C. J.; Kormos, B. L.; Winget, P.; Audette, V. M.; Beebe, J. M.; Brauer, C. S.; Burdick, W. R.; Cochran, E. W.; Eklov, B. M.; Giese, T. J.; Jun, Y.; Kesavan, L. S. D.; Kinsinger, C. R.; Minyaev, M. E.; Rajamani, R.; Salsbury, J. S.; Stubbs, J. M.; Surek, J. T.; Thompson, J. D.; Voelz, V. A.; Wick, C. D.; Zhang, L. J. Chem. Educ. 2001, 78, 1202. 10. Case, D. A; Pearlman, D. A.; Caldwell, J. W.; Cheatham, T. E., III; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. AMBER 6.0, University of California, San Francisco, 1999. 11. Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford Science Publications: New York, 1987. 12. Barton, A. F. M. J. Chem. Educ. 1973, 50, 119. 13. Contreras, M.; Valenzuela, J. J. Chem. Educ. 1986, 63, 7. 14. Laidler, K. J.; Meiser, J. H.; Sanctuary, B. C. Physical Chemistry; Houghton Mifflin Company: Boston, MA, 2003; pp 907–908. 15. Atkins, P.; de Paula, J. Physical Chemistry, 7th ed.; W. H. Freeman & Company: New York, 2002; pp 709–711. 16. Berry, R. S.; Rice, S. A.; Ross, J. Matter in Equilibrium—Statistical Mechanics and Thermodynamics, 2nd ed.; Oxford University Press: New York, 2002; pp 623–628. 17. Utilities exist within Amber to transform the xyz-formatted coordinate files into the PDB (protein database) format. Thus, any graphics program that will allow students to depict PDB files as ball-and-stick representations will work. At Holy Cross, we are currently using Spartan ’02 for UNIX (Wavefunction, Inc., Irvine, CA). However, programs such as VMD [http:// www.ks.uiuc.edu/Research/vmd (accessed May 2004)] or CHIME [http://www.mdl.com/chime (accessed May 2004)], which are available free of cost for the Windows environment, will also work. 18. Berendsen, H. J. C.; Gigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. 19. Jorgensen, W. L.; Severance, D. L. J. Am. Chem. Soc. 1990, 112, 4768. 20. Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. 21. Bohm, H. J.; McDonald, I. R.; Madden, P. A. Molec. Phys. 1983, 49, 347. 22. CRC Handbook of Chemistry and Physics; Linde, D. R., Ed.; CRC Press: Boca Raton, FL, 1999; pp 6–136.
Vol. 81 No. 9 September 2004
•
www.JCE.DivCHED.org
