Program ThermalConductivity> simulates the behavior of particles moving in a three-dimensional box interacting through a Lennard-Jones potential. The motion evolves using the numerical solution of the differential equations in Newton's second law. The velocity Verlet algorithm is used to update the positions and velocities of the particles. Periodic boundary conditions are used so that particles that leave the box enter on the other side, much like many video games.

In the program the length of the system in the $z$-direction is three times longer than the other two directions to allow for a better estimate of a temperature gradient in the $z$-direction. The system is divided into slabs in the $z$-direction. The temperature gradient is created by an energy flux due to swapping the slowest moving particle in a “hot slab” with the fastest moving particle in the “cold slab.” The thermal conductivity is determined from the ratio of the energy flux to the temperature gradient.

To begin the simulation we first equilibrate the system at a user specified temperature. Stopping the equilibration starts the swapping, thus creating an energy flux and temperature gradient.

Program ThermalConductivity computes the thermal conductivity $\kappa$ for a system of particles interacting with a Lennard-Jones potential cut off at $3\,\sigma$. The program uses the algorithm by Müller-Plathe that we have discussed. The total number of particles is $N=3k^3$, where the number of slabs in the $z$ direction is $3k$. The default values are $k=10$, $T=0.7$, and $\rho = 0.849$. Run the program until the system reaches equilibrium at constant $T$. (The program uses velocity rescaling to bring the system to the desired temperature.) When the Stop equilibration button is pressed, velocity rescaling is stopped, and energy transfer from the slab at $z=0$ to the middle slab is started. Press the Reset Data button several times to check that the simulation has reached a steady state. The temperature profiles in the left and right halves of the box are plotted separately to give you an idea of how well the averages have converged.

1. What is the estimated value of $\kappa$ for the default values of $N$, $\rho$, and $T$?
2. How do you expect $\kappa$ to change with increasing $T$ for fixed $\rho$? Test your expectation by doubling the temperature.
3. How do you expect $\kappa$ to change with increasing $\rho$ for fixed $T$? Test your expectation by increasing the density by $10\%$.
4. *Increase $N$ while keeping the density and the temperature the same as in part (a) and determine the importance of finite size effects on the computed value of $\kappa$.

Resources

1. Problem 10.15 in Statistical and Thermal Physics: With Computer Applications, 2nd ed., Harvey Gould and Jan Tobochnik, Princeton University Press (2021).
2. Florian Müller-Plathe, “A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity,” J. Chem. Phys. ${\bf 106}$, 6082–6085 (1997).