Molecular Dynamics of Lennard-Jones fluids.

Recover the Ex. 3 of Lecture XI, but using the integration method indicated in the following.

Use the Verlet Leapfrog algorithm: 
r(t + δt)  =  2r(t) – r(t – δt) + δt2a(t) 
where r(t) and a(t) are the (vector) position and acceleration of the particle at time t, respectively. 
The algorithm thus provides a way to take the positions of the particles at one time and generate 
their new positions at a later time; the velocities of the particles can be computed by: 
v(t)  =  [r(t + δt) – r(t – δt)]/2δt

 Run the code for the microcanonical Lennard-Jones fluid with at least 108 
particles (and more if you like) at the following thermodynamic state points: 

Reduced density  Reduced temperature 
liquid 1:  			0.7                     1.1 
liquid 2:    			0.85                   1.1 
gas:       			0.05                   1.1 
supercritical fluid:   	0.4                     1.4

and calculate each of the following quantities for each thermodynamic state point: 
1. The velocity autocorrelation function, Cvv(t)  =  〈v(0)•v(t)〉/〈v**2〉.  (Compute it for the 
velocity of each particle , and then average over the particles in the simulation).  Be sure to run your simulations for long enough that you have  good statistics at both early and later times (i.e. your trajectory should be roughly 100 times longer than the longest time for which velocity correlations persist). 
2. The mean-square displacement of the particles, R**2(t)  =  〈r(t) – r(0)|**2〉.  As discussed in the notes, your plots should show ballistic behavior at early times, and then linear 
behavior at long times.  Use a least-squares fit of the slopes of your plots at long times to 
calculate the diffusion constant, D. 
3. The radial distribution function, g(r)  =  〈ρ(0)ρ(r)〉/ρ**2.  (See e.g. Chandler p. 199)
4. The force autocorrelation function, Cff(t)  =  〈f(0)•f(t)〉/〈f**2〉.  Physically interpret all of 
your plots; how does Cff compare to the Cvv you calculated above? 

Briefly discuss how long you ran each simulation and how this affected the statistics of all your calculated results. Why did you have to run each simulation as long as you did?  
How are the properties of the various state points different from each other? 
How well does the diffusion constant you get from integrating Cvv 
compare with that obtained from the slope of the mean square displacement? 

Note:  since the simulations are done at constant E, V and N, you cannot run trajectories at 
precise temperatures.  This means that you will have to estimate the temperature for each 
thermodynamic state point using the output of the simulations.  Be sure to show how you 
estimated your temperature for each state point. 

Many checks are possible with existing literature, for example: J. Phys. Chem. B 104 6240; ibid. 6248; ibid., 6258.