Program Viscosity simulates the behavior of particles moving in a three-dimensional box interacting through a Lennard-Jones potential. The program uses the velocity Verlet algorithm to update the positions and velocities of the particles according to Newton's second law. Periodic boundary conditions are used.

In this 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 the $x$-velocity gradient in the $z$-direction. The system is divided into slabs in the $z$-direction. The $x$-velocity gradient is created by a momentum flux due to swapping the smallest $x$-velocity in a “high momentum slab” with the largest x-velocity in the “low momentum slab.” The viscosity is determined from the ratio of the momentum flux to the $x$-velocity gradient.

To begin the simulation we first equilibrate the system at a user specified temperature. We then stop the equilibration, which starts the swapping of velocities, creating a momentum flux and $x$-velocity gradient.

Program Viscosity uses $\eta = - \frac{P_{xz}}{(\Delta v_x/\Delta z)}$ to determine the viscosity for a system of particles interacting via the Lennard-Jones potential with a cutoff equal to $3\,\sigma$. The 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 the desired constant value of $T$. (The program uses velocity rescaling.) When the Stop equilibration button is pressed, the process of momentum transfer from slab at $z=0$ to the middle slab is started. Press the Reset Data button several times to ensure that the system has reached a steady state. The mean $x$ momentum in the left and right halves of the box are plotted separately so that you can determine how well they agree.

1. What is the value of the viscosity for the default values of $N$, $\rho$, and $T$?
2. How do you expect $\eta$ to change with increasing $T$ for fixed $\rho$? Test your expectation using Program Viscosity by doubling the temperature.
3. How do you expect $\eta$ to change with increasing $\rho$ for fixed $T$? Test your expectation using Program Viscosity 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 $\eta$.

Resources

1. Problem 10.17 in Statistical and Thermal Physics: With Computer Applications, 2nd ed., Harvey Gould and Jan Tobochnik, Princeton University Press (2021).
2. Florian Müller-Plathe, “Reversing the perturbation in nonequilibrium molecular dynamics: An easy way to calculate the shear viscosity of fluids,” Phys. Rev. E ${\bf 59}$, 4894–4898 (1999).