Calculating the Thermal Conductivity

Here is one possible answer to our question. Since an MD simulation shows how a system will evolve in time, we might consider calculating the thermal conductivity by applying a temperature gradient, dT/dx, to the system and then seeing if there is a flux of heat, J. You might imagine this as a molecular-level equivalent of the thought experiment in the page discussing macroscopic basics of heat transfer, where we imagined heat flux through a metal sheet. There we imagined the metal sheet with its two large faces in contact with thermal baths at two different temperatures. Now at a much smaller length scale, we could imagine placing walls at the left and right sides of the simulation box and maintaining them at different temperatures. Then when a molecule bounced off of one of the walls it could transfer energy to or from the wall, if the molecule were hotter than or colder than the wall, respectively. We could watch the flow of heat through the fluid of molecules between the two walls. The heat flux could be obtained by creating an invisible vertical plane inside the simulation, where we would keep track of the velocities of molecules moving across the plane. The vector sum of these velocities would be simply related to the net flux of heat moving from the hot side of the simulation to the cold side (remember the higher the temperature, the faster the atoms move). Of course, many molecules would be moving across the plane in both directions (toward the left and right), but overall we would see more energy moving in one direction (toward the cold wall) than the other.

There is, however, a fairly big problem with this approach. Remember that we're dealing with a small system here of maybe a few thousand atoms. Introducing a wall into such a system will be a rather large perturbation to the system. Molecules near surfaces experience a different environment than those in the “bulk” surrounded only by other fluid molecules. For example, molecules near a surface may be less mobile or may pack more densely. In a macroscopic system of order 10²³ molecules, only a negligible fraction are at the boundaries, but we cannot simulate 10²³ molecules. With a small system, we would have a large fraction of our molecules in this surface environment rather than in the bulk, so we wouldn't really calculate the properties of the bulk fluid that we're interested in. 

Periodic Boundary Conditions

To avoid this problem, we would like to use a method that does not introduce artificial walls into the system. The simulations in this module, and most molecular simulations in general, use a trick called “periodic boundary conditions” to avoid these spurious edge effects. If you’ve looked at the simulation applet here, you may have noticed that when a molecules leaves the right side of the box, another molecule enters at the left. (This is also used in some video games.) In effect, the simulation box is surrounded by identical images of itself, so that the box feels like it is in the middle of an infinite phase.

Algorithm to Calculate the Thermal Conductivity

Now, with periodic boundary conditions, how do we calculate the thermal conductivity? We’ve gotten rid of any “walls” surrounding the system, which is now considered to be effectively infinite in size, so the simple scheme suggested above will not work. Over the years, researchers have suggested several methods, as described by Allen & Tildesley (1987). [Click for more details.] The method we use here was suggested in 1997 by Mueller-Plathe. Our discussion closely follows his paper. We again rely on the basic equation:

To impose a heat flux on the system, the simulation box is divided in N “slabs” as shown in the figure below. The instantaneous temperature in each slab can be calculated from the kinetic energy of the molecules in that region. In three dimensions we have

We will define slab 0 as the “cold” slab and slab N/2 as the “hot” one. The heat flux is created by perodically exchanging the velocity vectors of a molecule in the cold slab and one in the hot slab in such a way that the temperature increases in the hot slab and decreases in the cold slab. This can usually be done by choosing the hottest molecule in the cold slab and the coldest molecule in the hot slab and swapping their velocities. This procedure transfers kinetic energy from one region to another; in other words it produces a flux of energy or heat. This leads to a temperature difference between the hot and cold slabs. Now what happens?

Having one end of the simulation hot and one cold is clearly not the equilibrium state. So, just like in a real system, the system tries to reach equilibrium. In this case, that means energy will be transferred to try to bring the entire system to the same temperature. But if we keep doing the velocity exchanges, we will always have hot and cold regions. The system will eventually reach a steady state where the energy transfer from the unphysical velocity exchanges is exactly balanced by a heat flux in the opposite direction. It is this heat flux that depends on the thermal conductivity.

Note that the amount of heat flux created each time the two velocity vectors are swaped is most likely to be different because the hottest atom in the cold slab and the coldest atom in the hot slab are different each time. However, the heat flux approaches a steady state value on average.

For the applet here, N is set to 8. You will see 8 boxes at the bottom of the applet that report the temperature or density in the slab directly above. You can toggle the buttons to report instantaneous temperatures or densities. To avoid bulk motion as disscussed before, it is desirable to keep the temperature gradient low (i.e. less swaps per unit time) so that there will not be too much density difference across the slabs (density differences causes bulk motion). To calculate the thermal conductivity, we divide the heat flux (imposed by the velocity exchanges) by the temperature gradient (that arises in the simulation). Basically,