The particles interact by the Lennard-Jones potential which looks like V(r)=V₀(r^-12-r^-6) where r is dimensionless. (In particular V(inf)=0, V(0)=+inf and there is a minimum at V(r=2^1/6)=-V₀/4) Hence, the Lennard-Jones potential is very repulsive when the particles approach each other at a distance less than r=0.9 but attractive in the range r=1..+inf with an attraction max at r=2^1/6.

In the simulation, the potential is set to 0 for distances bigger than a certain cutoff raduis (normally 2.5 to 3) to limit the "force range". Note that the periodic boundary conditions make force calculations a bit tricky since particles can interact across the walls of the box.

The graph on the left shows V(r)/V₀.

First, there are several pages of simulation dumps which are mainly interesting for debugging and seeing what is going on and represent the simulation in action and no final results:

• Red crosses display an x-y dump (z axis projected) of all particles.
• The curve on the top is the current PCF.
• The blue curve is the averaged PCF with error bars. This is not displayed for the first pcf_skip cycles.
• The curve on the bottom is normalized particle velocity versus particle index and can be used to detect errors in the simulation like excessively fast particles due to too large time steps.

Then, there is one page with the computed PCF which is the most important final result of the simulation. The title displays the number of particles (N), the density (rho), the number of iterations performed (iter), the number of PCF samples taken (pcf_samps), the averaged kinetic energy (<E_kin>), the (averaged) temperature (T) and finally how many seconds the complete simulation took.

Note that the PCF should always start at zero and converge to 1 for large radi.

Finally, there are 3 pages which display the pressure, the total energy per particle and the temperature as a function of simulation time. These 3 pages display essentially the same functions but the first page only displays the first 1% of the sim, the second one the first 10% and the last one displays the complete sim.

The dotted lines are the averages of the three curves (real values on the right top); note that collecting values for the average is started after pcf_skip iterations.

The most important thing is to ensure that pcf_skip is so large that it skips the initial iterations where the equilibrium state is not yet reached. Also, the total energy per particle (green curve) should stay constant throughout the sim (energy conservation) unless COOL_FACT is switched on, of course.

Solid

Liquid

Gaseous

First, the fluid in solid state: One can see the high peak which is near the minimum of the Lennard-Jones potential (r=2^1/6=1.1225). There are no particles less distant than 1 since the potential is repulsive in this area. One can see a double peak near r=2 which are the second- and third-next particles and a further peak near 3 for the fourth-next neighbours.

The liquid state with an inital speed around 1 (0.8; any small value will do, even 0 as long as active cooling using COOL_FACT is switched off).

Note the characteristic smooth "damped wave" (sort of). It represents the transition between the highly ordered solid state and the unordered gas state.

Finally, the gas state. Particles can approach each other more closely (like r=0.4) since their kinetic energy is very large. Note that there is very little distance-related correlation: The PCF changes pretty sharply from 0 to 1 which is a sign for the lack of order in that phase.

Due to cooling, the energy is not conserved and drops to negative (binding) energy. Temperature drops to nearly zero, and pressure settles down to a minimum given by the Lennard-Jones potential and the density. (One can have negative pressure when the particle density rho is so small that the particle distance is smaller than the radius of the minimum in the LJ potential.)

In the T,E,p versus time plot, one can see that equilibrium is reached after about 400 simulation cycles. Then, all parameters stay fairly constant.

Due to the much larger speed of the particles (speed 1000 times larger than for liquid but time steps only 50 times smaller time step), the equilibrium is reached much faster (about 50 cycles). However, thermal fluctuations are much larger.

Parameters:

int pcf_samples=8;
int pcf_skip=10000;
int dumpfreq=100;
int currfreq1=10;
int currfreq0=2;
double pcf_sigma=0.01;
sim.Initialize(/*N=*/1000,
    /*dt=*/1e-3,
    /*rho=*/1.05,
    /*cutoff_r=*/3.5,
    /*initial_v=*/10);
sim.InitializePCF(200,3.5);

The initial velocity was set to 10 and the system was cooled down during the first 10000 iterations (dt=1e-3) using COOL_FACT *0.997.

Parameters:

int pcf_samples=8;
int pcf_skip=400;
int dumpfreq=100;
int currfreq1=10;
int currfreq0=2;
double pcf_sigma=0.01/1.2;
sim.Initialize(/*N=*/1000,
    /*dt=*/0.5e-3,
    /*rho=*/0.95,
    /*cutoff_r=*/3.5,
    /*initial_v=*/0.8);
sim.InitializePCF(300,3.5);

Parameters:

int pcf_samples=8;
int pcf_skip=400;
int dumpfreq=100;
int currfreq1=10;
int currfreq0=2;
double pcf_sigma=0.01;
sim.Initialize(/*N=*/1000,
    /*dt=*/1e-5,
    /*rho=*/0.95,
    /*cutoff_r=*/3.5,
    /*initial_v=*/1000);
sim.InitializePCF(200,3.5);