Comparing Classical Water Models Using Molecular Dynamics To

Apr 16, 2018 - Optionally, the free molecular visualization program Chimera(24) can be used by students to qualitatively see the difference between th...
0 downloads 8 Views 2MB Size
Laboratory Experiment Cite This: J. Chem. Educ. XXXX, XXX, XXX−XXX

pubs.acs.org/jchemeduc

Comparing Classical Water Models Using Molecular Dynamics To Find Bulk Properties Laura J. Kinnaman,† Rachel M. Roller,‡ and Carrie S. Miller*,‡ †

Department of Mathematical Sciences, Morningside College, 1501 Morningside Avenue, Sioux City, Iowa 51106, United States Department of Biology and Chemistry, Azusa Pacific University, 901 E. Alosta Avenue, Azusa, California 91702, United States



S Supporting Information *

ABSTRACT: A computational chemistry exercise for the undergraduate physical chemistry laboratory is described. In this exercise, students use the molecular dynamics package Amber to generate trajectories of bulk liquid water for 4 different water models (TIP3P, OPC, SPC/E, and TIP4Pew). Students then process the trajectory to calculate structural (radial distribution functions) and dynamic (diffusion coefficients) properties of water that they compare to experimental values and to the other models. On the basis of these comparisons, students are also able to draw conclusions regarding the relative efficacy of these water models at modeling properties of real water. KEYWORDS: Upper-Division Undergraduate, Physical Chemistry, Computer-Based Learning, Computational Chemistry, Molecular Mechanics/Dynamics, Molecular Properties/Structure, Water/Water Chemistry, Laboratory Instruction



BACKGROUND Undergraduate students are often exposed to chemical properties through either individual molecules (e.g., dipoles

bonds and, importantly, the hydrogen-bonding network in water. However, they are highly efficient and can produce long simulations with many molecules in a reasonable time frame. Water is an especially interesting molecule to study, as it has a deceptively simple molecular structure, is vital to many chemical reactions, and demonstrates a hydrogen-bonded network. Because different properties of water are most relevant for different chemical problems, a vast array of computational methods have been developed to achieve different results. Our focus is on comparing four common water models: TIP3P,1 TIP4Pew,2 OPC,3 and SPC/E.4 Choices must be made when a model is developed, depending on what properties are important for the model to reliably match. A model might, for example, correctly describe structural properties such as distances between molecules and orientation of molecules, but it might then be unable to reproduce dynamic properties such as diffusion. The overall accuracy of a model can be characterized by comparing its structural, dynamical, and even spectroscopic properties to experimental data or quantum mechanical calculations. Being able to determine when a model is most effective is vital for scientists deciding which water model to use when doing simulations of complex processes such as protein folding,5−9 and interactions of water with semiconductors10−12 or biochemical surfaces.13,14 Although laboratory experiments and articles have been published in this Journal calculating distribution functions,15 transport properties,16,17 or simulations of water,18,19 this is the first paper to combine these calculations into an examination of

Table 1. Diffusion Coefficients of Water, Found Using the Diaphragm-Cell Technique22 temp (°C)

D (×10−9, m2/s)

1 4 5 15 25 35 45

1.149 1.276 1.313 1.777 2.299 2.919 3.575

and geometry) or bulk properties (e.g., boiling point or density) and may have difficulty seeing how the molecular and bulk properties relate to each other. Molecular dynamics (MD) represents a unique opportunity for students to connect these two ideas and see how the behavior of many molecules together, following simple physics laws, can give rise to the bulk properties that are measured. Because real molecules are quantum mechanical, it is necessary to solve the Schrödinger equation to fully capture their behavior. By using a few approximations, ab initio models do so; however, these calculations are limited to very small systems of just a few molecules. In order to actually model systems with enough molecules to reproduce the bulk state, molecules and their dynamics must be treated classically, using Newtonian mechanics. Classical models cannot capture certain effects which rely on electron behavior, such as the breaking of © XXXX American Chemical Society and Division of Chemical Education, Inc.

Received: June 14, 2017 Revised: March 26, 2018

A

DOI: 10.1021/acs.jchemed.7b00385 J. Chem. Educ. XXXX, XXX, XXX−XXX

Journal of Chemical Education

Laboratory Experiment

have random distributions of molecules and thus no peaks corresponding to the structure of neighboring atoms. Liquids have local structure and global disorder, so the RDF shows a few peaks (corresponding to solvation shells) before flattening.

the relative ability of different water models to predict and describe real properties. Given the ever-increasing use of computational chemistry in the undergraduate curriculum, it is vital that students learn not only how to use common calculational tools, but also how to evaluate the validity of the models built into those tools. In this experiment, students equilibrate bulk liquid water, using their assigned model, under various conditions, analyzing an example dynamic property (the diffusion coefficient) and structural property (the radial distribution function). By examining these properties under different conditions, students gain insight into the bulk behavior of water as it emerges from the dynamics of many individual molecules. By comparison to experimental results and to the other models, students gain an insight into how well their water model works under different conditions.



STUDENT LEARNING OUTCOMES The following are the learning goals for the exercise, to be assessed by the instructor. By the end of the laboratory students will be able to... 1. Produce MD simulation data of equilibrated water 2. Calculate RDFs of water from MD simulation data under different conditions 3. Calculate diffusion coefficients of water from MD simulation data at different temperatures 4. Interpret the physical meaning of RDFs and diffusion coefficients 5. Evaluate the ability of a classical water model to recreate experimental properties of water under different conditions

Diffusion Constant

Water molecules wander from their original positions because of thermal fluctuations, which can cause them to move both away from and toward that original position. Eventually, however, at a rate characterized by D, the diffusion coefficient, the water molecule will be far from where it started. The diffusion coefficient is a proportionality constant which relates the flux of particles to the concentration gradient. Effectively, particles can move faster through a less dense liquid and in the direction where the liquid is becoming less dense. For long time scales, the diffusion constant can be found using Einstein’s relationship:20 D=



EXPERIMENT This experiment consists of three parts: (I) Students learn how to use Amber 1623 to produce an MD trajectory for their

2 ⃗ |⟩ 1 ⟨| ri (⃗ t ) − ri (0) 6 t

where r is the position of the ith particle at time t. The quantity ri(t) − ri(0) is the difference between the position of molecule i at some later time t and that molecule’s original position at time zero. The angle brackets ⟨ and ⟩ indicate an average over all of the molecules in the simulation at time t. Diffusion coefficients also have a proportional relationship with absolute temperature, as predicted by the Stokes−Einstein equation21

D=

kB T 6πηa

which describes a fluid with dynamic viscosity η, no turbulence, and spherical particles of radius a. Although water is not a sphere, it has no turbulence in MD simulations. Note that the dynamic viscosity also depends on temperature, so the relationship between diffusion coefficient and temperature is not perfectly linear. In addition, the diffusion coefficient of water at various temperatures has been found experimentally (Table 1). Radial Distribution Function

Figure 1. (a) Ice-h polymorph of water, showing the underlying hexagonal structure. (b) Equilibrated box of water at 300 K.

The radial distribution function (RDF) describes the average number of specific atoms in the volume of a spherical shell of thickness dr and radius r from a given atom within a molecule, compared to the liquid’s average density.20 As such, there is an oxygen−oxygen RDF, oxygen−hydrogen RDF, and hydrogen− hydrogen RDF for water; this lab focuses on the oxygen− oxygen RDF. The intensity of peaks in an RDF indicate by how much the local density (for a given solvation shell) exceeds the bulk density. Qualitatively, the RDF can also show how organized the structure is. Solids have regular intervals between molecules and hence repetitive peaks in the RDF while gases

assigned water model and how to calculate the RDF and diffusion coefficient from the trajectory. (II) Students calculate the diffusion coefficients of liquid water over a range of temperatures and determine which model best matches the experimental dependency of D on temperature. (III) Students calculate the oxygen−oxygen RDFs of solid-, liquid-, and gasphase water and use their results to determine the phase of water at 250 K and 100 atm. B

DOI: 10.1021/acs.jchemed.7b00385 J. Chem. Educ. XXXX, XXX, XXX−XXX

Journal of Chemical Education

Laboratory Experiment

Figure 2. Mean displacement squared vs time at 280, 298, 330, and 360 K for (a) OPC, (b) SPC/E, (c) TIP3P, and (d) TIP4Pew water models. Dashes correspond to the linear fit of the data. For TIP4Pew (280 K) and OPC (280 and 300 K), only the first 10 ps are included in the fit. Temperatures are the target temperature for the MD simulation.

Research Questions

constant pressure. The calculations together have about a 20 min run time on a dual processor iMac. Further computational details are available in the Instructor Notes in the Supporting Information. Optionally, the free molecular visualization program Chimera24 can be used by students to qualitatively see the difference between the original ice structure and their equilibrated system (see Figure 1). In part I of the lab, students learn how to use Amber’s builtin postprocessing tool to pull the coordinates of the atoms from the trajectory files and calculate the position data needed to find diffusion coefficients and radial distribution functions. In part II, students determine the diffusion coefficient from the slope of the graph of mean displacement squared versus time for a range of temperatures. They fit the experimental diffusion coefficients versus temperature to a curve, which is then compared to the calculated values of each model to determine which model best predicts the temperature dependency of the diffusion coefficient. In part III, students calculate the RDF for their model under conditions where water is expected to be a solid (10 K, 1 atm), a liquid (298 K, 1 atm), and a gas (500 K, 1 atm) and use these results to assign the phase of their model at 250 K and 100 atm. Each part of the lab takes about 90 min per pair. If lab time is limited, the instructor can have the students perform only part II or part III as desired. Several components of the software used in this laboratory can be replaced with alternates as needed: the simulations and analysis can be performed with any molecular dynamics package (e.g., CHARMM25 or NAMD26).

The research questions, which students answer during the course of this lab, are the following: 1. How well do different water models reproduce experimental results for liquid water’s dynamics and structure at 298 K and 1 atm? 2. How well does each model reproduce the dynamics of liquid water through a range of temperatures? 3. How well does each water model reproduce the structure of water at a given temperature and pressure? Specifically, can the model correctly predict the phase of water at 250 K and 100 atm? Students answer the first and second questions on the basis of their calculated diffusion coefficients (a dynamical property). They answer the first and third questions on the basis of their calculated RDFs (a structural property). The comparison of both results to experiment and to the results of the other water models allows them to evaluate their assigned model in terms of its ability to represent structure and/or dynamics of water accurately. Process

The molecular dynamics package Amber16 is used to construct a box of 432 water molecules (using the OPC, SPC/E, TIP3P, or TIP4Pew model), defrost the system to the target temperature (at constant volume), equilibrate the system, and obtain the trajectory from which the bulk properties can be calculated. Both equilibration and production are run at C

DOI: 10.1021/acs.jchemed.7b00385 J. Chem. Educ. XXXX, XXX, XXX−XXX

Journal of Chemical Education

Laboratory Experiment

Visualization of the PDB files can be done with many different software packages (e.g., Chimera or VMD27). Graphing the results and the RDFs can likewise be done in a variety of software packages (e.g., Grace28 or Excel).

Table 3. Average Densities of Water Models at 1 atm and Various Temperatures Average T (K)

■ ■

Density (g/cm3)

Average T (K)

OPC

HAZARDS There are no hazards involved in this experiment.

9.8 243a 272 285 320 350 490

RESULTS AND DISCUSSION The calculated mean square displacement versus time data for the four water models are given in Figure 2. Linear fits of the

Average Temperature 272 285 320 350

DModel (×10−9, m2/s) 0.953 0.989 3.868 5.628

± ± ± ±

0.003 0.003 0.004 0.003

DExp (×10−9, m2/s) OPC 1.076 1.626 3.750 6.315

± ± ± ±

0.007 0.007 0.007 0.007

9.8 246a 274 287 322 355 485

(DModel − DExp)2 0.014 ± 0.002 0.411 ± 0.010 0.022 ± 0.002 0.320 ± 0.009 sum = 0.767 ± 0.013

274 289 320 349

1.440 2.219 4.355 5.437

± ± ± ±

SPC/E 0.001 1.152 0.001 1.822 0.002 3.750 0.002 6.218

± ± ± ±

0.007 0.007 0.007 0.007

0.084 ± 0.004 0.154 ± 0.006 0.401 ± 0.009 0.440 ± 0.010 sum = 1.080 ± 0.015

274 287 322 355

3.294 3.635 6.976 8.711

± ± ± ±

0.002 0.001 0.002 0.002

TIP3P 1.152 1.723 3.900 6.809

± ± ± ±

0.007 0.007 0.007 0.007

4.599 ± 0.032 3.640 ± 0.027 9.670 ± 0.045 4.183 ± 0.030 sum = 22.092 ± 0.068

274 287 323 351

0.944 2.074 4.139 5.829

± ± ± ±

TIP4Pew 0.003 1.152 0.003 1.723 0.003 3.975 0.004 6.412

± ± ± ±

0.007 0.007 0.007 0.007

0.042 ± 0.003 0.121 ± 0.005 0.039 ± 0.003 0.209 ± 0.008 sum = 0.412 ± 0.010

SPC/E 0.9197 0.8949 0.8902 0.8882 0.9946 0.9750 0.8522

9.8 246a 274 289 320 349 486

0.9846 1.0284 1.0040 0.9972 0.9597 0.9283 0.7225

9.7 240a 274 287 323 351 485

TIP3P

Table 2. Calculated Diffusion Coefficients for Liquid Water at Various Temperatures

a

Density (g/cm3) 0.9752 0.9471 1.0109 1.0113 0.9805 0.9665 0.8080

TIP4Pew 0.9701 0.9404 0.9309 0.9899 0.9867 0.9631 0.8211

100 atm.

only 289 K, so it is possible the model is less accurate than the initial result suggests. Second, the results from the simulation are outputted in terms of timesteps and must be converted to time units. The Amber input file for the production run (available in the Supporting Information) uses a time step of 1 fs. Third, the results in Table 2 are from single-trajectory runs, and thus, the results may vary from those of trajectories starting with a significantly different initial velocity distribution (although the trends will be the same). Further details are available in the Instructor Notes in Supporting Information. In order to determine which model best predicts the temperature dependency of water’s diffusion coefficient with water, the experimental values (Table 1) were fit to a secondorder polynomial. This fit equation, D = 0.0003485T2 − 0.1511T + 16.38, was used to determine what the experimental value would be at the average temperature associated with each simulation (DExp). Finally, in order to compare the models’ results, we find the square of the difference between DModel and DExp. The closer the sum of (DModel − DExp)2 is to zero, the better the water model describes the temperature dependency of diffusion. The best results are from TIP4Pew, while the worst are from TIP3P. SPC/E is significantly better than TIP3P, but worse than OPC and TIP4Pew, even though its prediction of the diffusion coefficient at 289 K is quite close to the experimental value for 298 K. Average temperatures and densities are given in Table 3. As expected, the densities fluctuate around 1 g/cm3, with the lowest values at 485 K. However, these gas-phase values are still 4 orders of magnitude greater than that of real water vapor. Because these models include nonzero interactions at all temperatures, the gas phase is nonideal. Figures 3 and 4 show the oxygen−oxygen radial distribution functions, organized by condition and by model, respectively. A qualitative comparison is sufficient for students to achieve their learning goal of interpreting the physical meaning of each RDF (although quantitative analyses such as the positions of local extrema and the integration of peaks can be done at the instructor’s discretion). As expected, the solid-phase results at 10 K and 1 atm (Figure 3a) display the most distinct and well-spaced peaks for all models. Although TIP4Pew, SPC/E, and OPC all have

data were generally performed for the entire 20 ps run with the exceptions of OPC at 280 and 298 K, and TIP4Pew at 280 K, which exhibited nonlinear behavior after 10 ps. However, extending the trajectory to 200 ps does not improve the results over using the first 10−20 ps. The slope of the mean square displacement versus time can vary in different regions of a long (200 ps) simulation,29 so the overall fit for our models was not improved using a long trajectory. In addition, the trend of best to worst model was the same regardless of trajectory length. The self-diffusion coefficients (DModel), calculated from the slopes of the displacement data, are given in Table 2. As expected, the coefficients increase as the temperature increases. A few details are important to note. First, the target temperature in the MD input file is not always equal to the actual average temperature of the simulation. For example, at the target temperature of 298 K, SPC/E (at 2.219 × 10−9 m2/s) has the closest result to the experimental value of 2.3 × 10−9 m2/s;22 however, the average temperature of the simulation was D

DOI: 10.1021/acs.jchemed.7b00385 J. Chem. Educ. XXXX, XXX, XXX−XXX

Journal of Chemical Education

Laboratory Experiment

Figure 3. Oxygen−oxygen RDFs for the four water models at (a) 10 K and 1 atm, (b) 298 K and 1 atm, (c) 500 K and 1 atm, and (d) 250 K and 100 atm.

similar first four or five peaks, TIP3P not only matches those but also clearly has additional peaks in the intermediate regions. TIP3P’s behavior is consistent with being understructured: the probability of finding another oxygen atom is less localized than for the other models. The liquid-phase results at 298 K and 1 atm (Figure 3b) show that SPC/E and TIP4Pew are similarly structured in terms of the height of their first peak and depth of their first minimum, while TIP3P is relatively understructured with both a smaller first peak and a nearly nonexistent first minimum. In contrast, OPC is overstructured, with a much larger first peak and zero density of oxygen atoms between the first and second peaks. In the gas phase, at 500 K and 1 atm (Figure 3c), all four models show a similar lack of structure as they flatten almost immediately after a first peak. This first peak is a result of the relatively high density (roughly 80% of the liquid phase, as shown in Table 3) of the gas. Similar to the liquid results, OPC retains the most structure, as it has a recognizable dip between its first peak and a second shoulder. At 250 K and 100 atm, water is a solid.30 However, not all of these four water models predict that phase under those conditions. While analysis of Figure 3d appears to indicate that the predicted phases are gas (TIP3P) and solid (TIP4Pew, SPC/E, and OPC), comparisons of each model to its own RDFs at different temperatures and pressures, as in Figure 4, give more accurate predictions. TIP3P is consistently understructured; its RDF at 250 K and 100 atm is most similar to its RDF at 298 K and 1 atm. Therefore, TIP3P predicts a liquid phase (because its liquid phase is gas-like). TIP4Pew and SPC/

E both predict a solid phase: In contrast to their liquid RDFs, they have regions of zero density between their first and second peak, as well as distinct third and fourth peaks. OPC is consistently overstructured; its RDF at 250 K and 100 atm is most similar to its 298 K and 1 atm result. Therefore, OPC predicts a liquid phase (because its liquid phase is solid-like). The result for TIP4Pew is consistent with that of Aragones and Vega.30



STUDENT RESULTS This laboratory experiment has been performed in a physical chemistry lab for three years, albeit with revisions for the last year. Students worked in pairs to perform the MD simulations for their assigned water model. The groups shared their data with their peers, and then wrote full lab reports individually. Students found learning and using Terminal commands to be challenging (it is recommended that instructors familiarize themselves with these commands beforehand). It is also important that students have time to read and ask questions about the background information before they make their hypotheses. Students were able to successfully answer the research questions for parts II and III on their own, but needed guidance to answer the research question for part I. Some students struggled with the final conclusions, but most were able to put forth some sort of suggested improvement for their model.



CONCLUSIONS We have developed a computational chemistry laboratory that is appropriate for a physical chemistry course. Students prepare, E

DOI: 10.1021/acs.jchemed.7b00385 J. Chem. Educ. XXXX, XXX, XXX−XXX

Journal of Chemical Education

Laboratory Experiment

Figure 4. Oxygen−oxygen RDFs at 1 atm and 10, 298, and 500 K and 100 atm and 250 K for (a) OPC, (b) SPC/E, (c) TIP3P, and (d) TIP4Pew water models.



equilibrate, and run production simulations of one of four water models using Amber, then analyze a structural quantity (the oxygen−oxygen radial distribution function) and a dynamical quantity (the diffusion coefficient) over a range of conditions. These quantities can be compared among each model and compared to experimental values in order to evaluate the relative accuracy of each water model.



(1) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79 (2), 926. (2) Horn, H.; Swope, W.; Pitera, J.; Madura, J.; Dick, T.; Hura, G.; Head-Gordon, T. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120 (20), 9665. (3) Izadi, S.; Anandakrishnan, R.; Onufriev, A. Building water models: a different approach. J. Phys. Chem. Lett. 2014, 5 (21), 3863. (4) Berendsen, H.; Grigera, J.; Straatsma, T. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91 (24), 6269. (5) Kauzmann, W. Denaturation of Proteins and Enzymes. In The Mechanism of Enzyme Action; McElroy, W. D., Glass, B., Eds.; The Johns Hopkins Press: Baltimore, MD, 1954; pp 70−110. (6) Jackson, S. How do small single-domain proteins fold? Folding Des. 1998, 3 (4), R81−91. (7) Southall, N.; Dill, K.; Haymet, A. A View of the Hydrophobic Effect. J. Phys. Chem. B 2002, 106 (3), 521−533. (8) Prabhu, N.; Sharp, K. Protein-Solvent Interactions. Chem. Rev. 2006, 106 (5), 1616−1623. (9) Sklenar, H.; Eisenhaber, F.; Poncin, M.; Lavery, R. Theoretical Biochemistry and Molecular Biophysics; Adenine Press: New York, 1990. (10) Fujishima, A.; Honda, K. Electrochemical Photolysis of Water at a Semiconductor Electrode. Nature 1972, 238, 37−38. (11) Fujishima, A.; Rao, T.; Tryk, D. Titanium dioxide photocatalysis. J. Photochem. Photobiol., C 2000, 1, 1−21. (12) Tryk, D.; Fujishima, A.; Honda, K. Recent topics in photoelectrochemistry: achievements and future prospects. Electrochim. Acta 2000, 45 (15), 2363−2376.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available on the ACS Publications website at DOI: 10.1021/acs.jchemed.7b00385. Lab handout, instructor notes, input files, perl script, PDB for initial ice structure, and list of commands (ZIP)



REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Carrie S. Miller: 0000-0003-0990-4383 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors thank Melanie Luallen for her help confirming the results. F

DOI: 10.1021/acs.jchemed.7b00385 J. Chem. Educ. XXXX, XXX, XXX−XXX

Journal of Chemical Education

Laboratory Experiment

(13) Furse, K.; Corcelli, S. The dynamics of water at DNA interfaces: computational studies of Hoechst 33258 bound to DNA. J. Am. Chem. Soc. 2008, 130 (39), 13103−13109. (14) Nadassy, K.; Wodak, S.; Janin, J. Structural features of proteinnucleic acid recognition sites. Biochemistry 1999, 38 (7), 1999−2017. (15) Speer, O. F.; Wengerter, B. C.; Taylor, R. S. Molecular Dynamics Simulations of Simple Liquids. J. Chem. Educ. 2004, 81 (9), 1330−1332. (16) Eckler, L. H.; Nee, M. J. A Simple Molecular Dynamics Lab to Calculate Viscosity as a Function of Temperature. J. Chem. Educ. 2016, 93 (5), 927−931. (17) Fuson, M. M. Anisotropic Rotational Diffusion Studied by Nuclear Spin Relaxation and Molecular Dynamics Simulation: An Undergraduate Physical Chemistry Laboratory. J. Chem. Educ. 2017, 94 (4), 521−525. (18) Saiz, E.; Tarazona, M. P. Molecular Dynamics and the Water Molecule: An Introduction to Molecular Dynamics for Physical Chemistry Students. J. Chem. Educ. 1997, 74 (11), 1350−1354. (19) Lykos, P. Modeling the Hydrogen Bond within Molecular Dynamics. J. Chem. Educ. 2004, 81 (1), 147−149. (20) Allen, M.; Tildesley, D. Computer Simulations of Liquids; Oxford University Press, Inc.: New York, 1987. (21) Einstein, A. Ü ber die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flü ssigkeiten suspendierten Teilchen. Ann. Phys. 1905, 322, 549. (22) Mills, R. Self-diffusion in normal and heavy water in the range 1−45.deg. J. Phys. Chem. 1973, 77, 685−688. (23) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, W.; Zhang, K. M.; Merz, B.; Roberts, S.; Hayik, A.; Roitberg, G.; Seabra, J.; Swails, A. W.; Götz, I.; Kolossváry, R. C.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, X.; Ye, J.; Wang, M.-J.; Hsieh, G.; Cui, D. R.; Roe, D. H.; Mathews, M. G.; Seetin, R.; Salomon-Ferrer, C.; Sagui, Q.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12; University of California: San Francisco, 2012. (24) Pettersen, E.; Goddard, T.; Huang, C.; Couch, G.; Greenblatt, D.; Meng, E.; Ferrin, T. UCSF ChimeraA Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25 (13), 1605−1612. (25) Brooks, B.; Brooks, C. L., III; Mackerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: the biomolecular simulation program. J. Comput. Chem. 2009, 30 (10), 1545−1614. (26) Phillips, J.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781. (27) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (28) Grace, a WYSIWYG 2D plotting tool for numerical data. http:// plasma-gate.weizmann.ac.il/Grace/(accessed Feb 2018). (29) Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105, 9954−9960. (30) Aragones, J.; Vega, C. Plastic crystal phases of simple water models. J. Chem. Phys. 2009, 130, 244504.

G

DOI: 10.1021/acs.jchemed.7b00385 J. Chem. Educ. XXXX, XXX, XXX−XXX