Gravity Leveling of a 2D Liquid

Here I model the leveling of a 2D liquid interface. The liquid discs interact via a Lennard Jones force the potential of which has the form
V(r)=4*epsilon*((sigma/r)^12-(sigma/r)^6)
where r is the separation between disc centers, sigma has units of length, and epsilon is the depth of the potential well. The first term provides a very strong repelling force and the second term provides a weaker attractive force. Because the repelling force is so strong, the liquid is almost in-compressible as one expect for a liquid. The borders of the container also have an exponential potential to keep the discs from escaping. The simulation takes several minutes to finish while storing some image and plot videos. I have therefore provided a movie so that the animation can be run more rapidly. I provide two plot windows: The x and y centers of mass of the discs and a histogram of the distribution of discs as a function of vertical height. I give a choice of either stepped or sloped starting distributions of the discs. In both cases the forces due to gravity tend to level the discs although the top layer(s) remain a bit ragged since then there is nothing to cause the top layer discs sideways except their thermal energy. Of course, the gravity fall of the initially top discs gives them a lot of energy, in fact too much energy for stability of the liquid. To avoid this problem I chose to provide a drag coefficient that is a force proportional to disc speed. This provides stability but also reduces the thermal speed of the discs which results in a more ragged top surface of the liquid. To improve the leveling of the liquid I provide an annealing period which is an increase of thermal speed near the end of the simulation. An important result that the learner should note is that the final geometry is that all but the boundary discs exhibit 6-fold symmetry where all the 6 nearest neighbors are the same distance and at 60 degree angles from each other. This is the most stable geometry and it roughly prevails even during the most erratic parts of the leveling process. The distances are just a little less than the zero force discance of Vsolid(r) because the liquid is slightly compressed by the continous downward gravity force. One should expect that the density of the lower part of the liquid will be higher than the upper parts because the lower parts have all the gravity forces of the column of discs above them. I attempted to show that the this was the case but result was inconclusive for the number of discs that I was able to model.a When playing the movie, note that some discs are vaporized when they hit the right side of the container. This is actually very realistic for the rather violent liquid movement that results initially. If you choose parameters that are too aggressive, then the container might lose some of the discs. If that happens, it will be reported but the results for the centroid position plot will not be affected by these losses.