Molecular Dynamics Exercises

Daniel V. Schroeder, September 2014

These exercises make use of the Interactive Molecular Dynamics web application. They are designed for use in undergraduate courses at the beginning and intermediate levels. The exercises illustrate and reinforce the concepts of thermal physics using numerical data for a nontrivial system, and highlight the idealizations ordinarily made in undergraduate courses, often showing when and how these idealizations break down. Some of the exercises could be developed further into advanced projects.

A more generic version of this collection of exercises appears as Section VI of the accompanying article in the American Journal of Physics. Whereas that version is intended for use with any sufficiently flexible molecular dynamics software, the version here includes specific instructions and hints for use with the accompanying web application. The exercises (like the software) use natural units, but could be augmented by asking students to convert results to conventional units. Many of the exercises require that students use a spreadsheet program or other plotting software.

1. Basic phase behavior. Run the simulation under a variety of conditions, clicking the Faster and Slower buttons (while the simulation is running) to add and remove energy, and using the sliders to adjust the number of atoms and the volume (which is actually an area in two dimensions). Be sure to wait for the system to equilibrate at each setting. Keep exploring different settings until you have produced each of the following states of matter: (a) a pure gas, with plenty of space around most atoms and no large clumps; (b) a pure liquid, with little space between atoms but no long-range order; (c) a pure solid, with all the atoms in orderly rows; (d) a liquid droplet surrounded by gas; (e) a solid crystal surrounded by gas. For one example of each of these states, write down the following data: number of particles, volume, total energy, kinetic energy, potential energy, temperature, and pressure. (Click the Data button to show the kinetic and potential energies.) Summarize your results in a paragraph, describing the appearance of each phase and any other interesting behavior that you notice.
2. Comparison to an ideal gas. Set the number of particles ($N$) to 100 and the volume ($V\!$, actually an area in two dimensions) to approximately 5000. With the simulation running, add or remove energy (using the Faster and Slower buttons) until the temperature ($T$) remains stable at approximately 1.0, then note the pressure ($P$) and compute the ratio $PV/Nk_BT$ (remembering that in the natural units used by the simulation, $k_B=1$). What is this ratio for an ideal gas, and how does your result compare? Repeat (using the same $N$ and $V$) for $T\approx0.5$ and $T\approx0.3$, being sure to remove enough energy for the system to equilibrate at these approximate temperatures. Describe the way in which this system’s behavior differs from that of an ideal gas, and explain the reason for this difference, noting the visual appearance of the system at each temperature.
3. Free expansion. Set up an experiment in which a nearly ideal gas (perhaps 100 atoms with an initial volume of about 5000) expands into a vacuum to approximately double its volume. (Pause the simulation before you drag the slider to increase the volume.) If the initial temperature is about 2.0, what is the final temperature, and why? Repeat the experiment with a smaller initial volume (perhaps 1000), and explain the results. Then try it with a smaller initial temperature as well (perhaps 1.0), and again explain the results.
4. Heat capacities and equipartition. Use the simulation to measure the heat capacity (at constant volume) of the Lennard-Jones system when it is (a) a nearly ideal gas, and (b) a single solid crystal, with at least 100 atoms in each case. For each of these measurements you will need to measure the temperatures at two slightly different energies ($E$), then subtract these nearby values to obtain $\Delta E$ and $\Delta T$. Thus, a good procedure is to note the energy and temperature, then click the Faster button to add a bit of energy, and note the energy and temperature again after the system has equilibrated. Compare your results to the predictions of the equipartition theorem, thinking carefully about how many degrees of freedom the system has. (Don’t expect perfect agreement, but don’t expect enormous disagreements either.)
5. Heat capacity at constant pressure. Use the simulation to measure the heat capacity at constant pressure of the Lennard-Jones system when it is (a) a nearly ideal gas, and (b) a single solid crystal. In each case you’ll have to increase the volume by a small percentage, then add energy until the pressure reaches its previous value (use the $+1\%$ and $-1\%$ buttons to make fine adjustments to the energy). In calculating the results, be sure to include the $PdV$ term in $C_P = (dE+PdV)/dT$. For the gas, compare your result to the exact formula derived in textbooks. For the solid, check that $C_P>C_V$.
6. Pressure and energy as functions of temperature. For a simulated system of 100 atoms or more, with the volume fixed in the range of 10 to 20 units per atom, measure the total energy, temperature, and pressure over the full range of temperatures from 0 to 1.0, in intervals of 0.1 or less. Use the Faster and Slower buttons to adjust the system’s energy, being sure to let the system equilibrate at each temperature before recording your data. Note that you can press the Data button to show a panel of data-gathering controls. With the default settings, the “Write data” button will record all the information you need, and you can copy and paste all the data into a spreadsheet when you’re finished. Then plot a graph of pressure vs. temperature and another graph of energy vs. temperature. Comment on the portions of these graphs that can be understood in terms of the ideal gas law and the equipartition theorem, and on the portions that cannot be so simply understood (and why). How does the low-temperature behavior of the heat capacity differ from that of a real-world solid?
7. Heat capacity and entropy. From your data in the previous problem, construct a graph of the heat capacity at constant volume, $C_V$, as a function of temperature. You may have to do some smoothing to reduce the effects of noise in the data. Then construct a table and graph of $C_V/T$ vs. $T$, and numerically integrate this function (to an accuracy of one or two significant figures) to determine the entropy as a function of temperature, relative to the entropy at $T=0.1$. Why can't you determine the absolute entropy, relative to $T=0$? Why doesn't this limitation affect real-world materials?
8. Critical point. Set the number of atoms to at least 1000 (more is better) and the volume, in natural units, to approximately three times the number of atoms. Add and remove energy (using the Faster and Slower buttons) to carefully explore the behavior of the system over the temperature range from about 0.4 to 0.7, and describe how its appearance changes over this range. What is your best estimate of the critical temperature of this system, and what is the corresponding pressure?
9. Phase diagram. Map out the approximate phase diagram of the two-dimensional Lennard-Jones system, by adjusting both the temperature and the volume to find the various phase boundary lines (where two phases coexist in equilibrium). Keep the number of atoms fixed (preferably at 500 or more). It’s easiest to start at a large volume and low temperature, so the system consists of a single solid crystal surrounded by a low-density gas. Use the Faster button to add energy gradually, letting the system equilibrate at various temperatures and recording the temperature and pressure (using the Data panel and “Write data” button) after each equilibration. Be sure to note the approximate triple point, where the solid crystal (with atoms in orderly rows) melts into a liquid (with no long-range order). The critical point is the subject of the previous problem. Finally, reduce the volume (and the energy) to try to locate the solid-liquid phase boundary at pressures somewhat above that of the triple point. Copy and paste all your data into a spreadsheet and then plot all of your pressure-temperature points. Print this plot, sketch in the approximate phase boundary lines, and annotate the plot with descriptions of the system’s appearance under the various conditions.
10. Phase boundary ambiguities. Phase boundary lines are sharp (because the properties of the system across a boundary are discontinuous) only in the limit of an infinitely large system. As a follow-up to the previous problem, explore how the phase boundary locations depend on the volume of the system and on the number of atoms. For example, try plotting the liquid-gas phase boundary for systems with different sizes but similar average densities. Explain the results qualitatively, by considering what fraction of the atoms in the liquid droplet are near the surface.
11. Velocity distribution. Record the instantaneous velocities ($x$ and $y$ components) for 1000 or more atoms in equilibrium at a temperature of about 0.5 in natural units. To record the data, press the Data button to show the data-gathering controls, then choose “All atoms” from the Data type menu, and press the “Write data” button. Copy and paste this data into a spreadsheet, then plot a histogram of the $v_x$ values, using about 20 bins to cover the velocity interval $-2.0$ to $+2.0$. (One way to make a histogram is to first create a column of the bin upper limits and then, in the next column, use the COUNTIF spreadsheet function to count how many entries in the list of $v_x$ values are less than each limit; use a third column to subtract these numbers to obtain the number of $v_x$ values in each bin, and plot this column on the vertical axis.) Do the same for the $v_y$ values. Also plot the expected results according to the Maxwell-Boltzmann velocity distribution, which for either component $v_i$ is $\sqrt{m/(2\pi k_B T)} \exp({-mv_i^2/2k_B T})$. (This is the function that, when multiplied by any small velocity interval $dv_i$, gives the probability of finding a single atom within this interval. To compare it to your simulation results, you’ll have to take into account the number of atoms and the sizes of the histogram bins.) Repeat this whole procedure for a different temperature (adjusting the histogram range if necessary), and discuss the results. Does it matter whether the simulated material is in a solid, liquid, or gas state?
12. Speed distribution. As in the previous problem, record the instantaneous velocities of 1000 or more atoms in equilibrium. Use a spreadsheet to calculate the speed of each atom, plot a histogram of the speeds, and compare to the theoretical prediction (i.e., the two-dimensional Maxwell speed distribution). Why does the speed distribution equal zero at $v = 0$, where the distributions for $v_x$ and $v_y$ have their peaks?
13. Gas density in a gravitational field. Set the size of the container to $100\times100$, the number of atoms to 500, the gravitational constant to 0.02, and the temperature (via the Faster and Slower buttons) to about 1.0. The system is now an “atmosphere” whose density decreases with altitude. How does the typical gravitational energy compare to the typical kinetic energy? Use the “All atoms” setting in the Data panel to record the positions of all the atoms at one instant, and copy this data into a spreadsheet. Then use the spreadsheet to plot a graph of relative density ($\rho$) as a function of altitude, dividing the full height of the container into ten bins. (For better statistics you may want to combine data from several “snapshots” of the system.) Fit an exponential function to this data (perhaps plotting it on a semi-log scale), and compare the result to the standard formula for an isothermal atmosphere, $\rho\propto \exp({-mgy/k_BT})$. Is this system isothermal? Would you expect it to be? Why or why not? What happens if you reduce the temperature to 0.5?
14. Thermal expansion. Devise a numerical experiment to measure the thermal expansion of a two-dimensional Lennard-Jones solid. (Hints: Make an elongated solid and use a small amount of gravity to hold it along the bottom of the container. Anchor an atom near one end and select an atom near the other end, then use the Data panel, with Data type “Selected atom,” to repeatedly record the position of this selected atom. Use a spreadsheet to compute the average distance between the anchored and selected atoms.) Explore the temperature range from 0 up to about 0.10, being sure to acquire enough data to see the signal through the statistical noise. Does the expansion depend linearly on temperature over this range? What is the approximate thermal expansion coefficient? How does this behavior compare to that of real solids at low temperatures?
15. Brownian motion. Set up a simulation of approximately 50 atoms arranged in a stable crystalline shape (not necessarily symmetrical), surrounded by plenty of empty space so the overall volume per atom is 10 or so. (To form a stable crystal, you may need to pause the simulation and drag some atoms into place.) Add or remove energy until the temperature is between 0.06 and 0.08, and run the simulation for a while. You should see the crystal bounce around randomly, with each bounce off of a wall changing its overall velocity as energy is exchanged between its macroscopic and microscopic degrees of freedom. Then use the Data panel, with “More detail” checked, to measure the system’s total momentum ($x$ and $y$ components) repeatedly, at regular intervals that are long enough for at least one or two bounces, on average, to occur between measurements. Make at least 100 such measurements (more is better). Copy this data into a spreadsheet and then calculate the average values of $p_x^2$ and $p_y^2$; compare these averages to the prediction of the equipartition theorem. Also plot histograms of the $p_x$ and $p_y$ values, and compare them to the prediction of the Maxwell-Boltzmann distribution.
16. Brownian bouncing ball. As in the previous problem, set up a simulation of a small solid crystal surrounded by empty space. This time, to minimize the effects of the crystal’s orientation, it’s best to make it as nearly round as possible. Freeze the crystal’s overall motion when it is near the middle of the region, then turn on a downward gravitational force of 0.001 in natural units. Run the simulation to let the crystal bounce off the bottom surface a few times, then adjust the temperature to somewhere in the range between 0.06 and 0.10. Use the Data panel to measure the system’s total gravitational energy at regular intervals (far enough apart in time for the crystal to move a significant distance). After making at least a few hundred such measurements, copy the data into a spreadsheet, plot a histogram of the gravitational energy values, and compare to the prediction of the Boltzmann distribution.
17. Fluctuations. The simulation calculates temperature by taking a time average of the average kinetic energy per particle. The time average is needed because the instantaneous average kinetic energy per particle fluctuates significantly for such a small system. To study these fluctuations, start with about 50 atoms in a volume of about 250, at a temperature of about 0.4, so the system consists of a liquid droplet surrounded by gas. While the simulation runs, use the Data panel to measure the total kinetic energy about a hundred times at regular intervals. Copy the kinetic energy values into a spreadsheet and for each, calculate the average kinetic energy per particle; then calculate the standard deviation of these measurements. Repeat this process for systems of about 100, 200, and 500 particles, keeping the density and temperature approximately fixed. How does the standard deviation vary with $N$? Next, hold $N$ fixed and repeat the process at lower and higher temperatures where the system is entirely solid or entirely gaseous. Describe and interpret your results as completely as you can.
18. Reversibility and chaos. Use the Presets menu to set up the configurations of Fig. 3(a), (b), and (c), and watch how each evolves over time. Reverse the motion after a short time, and check whether the reversed motion restores the initial configuration (at least approximately). In each case, determine the approximate time limit beyond which a reversal will not restore the initial state. (You may need to reduce the number of steps per frame, and/or pause the simulation and use the Step button.) Then set up the configuration of Fig. 4, and determine the approximate number of bounces before the motions of the two molecules are no longer (approximately) synchronized.
19. Thermal conductivity. Fill the simulated space with a solid of several hundred atoms, with the left half at a temperature of about 0.1 and the right half at $T=0$. (Setting up this initial state from scratch is somewhat laborious, but it is already set up for you in the Presets menu.) Starting with this out-of-equilibrium state, step the system forward by small time increments and watch it equilibrate. Save the state periodically during this process, and use a spreadsheet to calculate the average temperature separately for the left and right halves of the system at each time. Plot these temperatures vs. time, and determine the approximate initial value of the slope of each graph, $dT/dt$. Finally, use this result to estimate the thermal conductivity, $k_t$, of the two-dimensional Lennard-Jones solid. This quantity is defined by the two-dimensional version of the Fourier heat conduction law, $Q/\Delta t = -k_t L \,dT/dx$, where $x$ is the coordinate along which the temperature varies, $L$ is the cross-sectional length of the material (measured perpendicular to $x$), and $Q$ is the amount of heat that flows across the boundary in time $\Delta t$. You'll need to know the heat capacity of the material, which you can determine from your data or from exercise 4 or 6 above.
20. Diffusion. Configure the simulation to model a fluid of 1000 atoms, in a volume of 1600, at a temperature of approximately 1.0. Select one atom that is initially near the center of the container and as the simulation runs, record its position at regular intervals of one time unit, for at least 200 time units. (If it reaches the edge of the container in less than 200 time units, discard the data and try again.) Then, using a spreadsheet or other computing environment, compute the squared displacement, $(\Delta x)^2+(\Delta y)^2$, for each of the (200 or so) one-unit time intervals in your data set. Average these values to obtain the mean squared displacement (MSD). Similarly, use the same data set to calculate the MSD for time intervals ($\Delta t$) of 2, 5, 10, and 20 units. Plot the MSD vs. $\Delta t$ and notice that the graph is approximately linear; this is the characteristic behavior of diffusive motion (or a so-called random walk). The slope of the line is closely related to the diffusion constant, $D$; in two dimensions, the MSD is $4D\Delta t$. Estimate the diffusion constant from your data, then repeat the analysis, holding the fluid density fixed, at temperatures of approximately 0.5 and 2.0.