210 Molecular dynamics simulations
Table 8.1. Molecular dynamics data for
thermodynamic quantities of the Lennard–Jones liquid.
ρ( 1 /σ^3 ) T 0 (ε/kB) T βP/ρ U(ε)
0.88 1.0 0.990 (2) 2.98 (2) −5.704 (1)
0.80 1.0 1.010 (2) 1.31 (2) −5.271 (1)
0.70 1.0 1.014 (2) 1.06 (4) (5) −4.662 (1)
T 0 is the desired temperature;Tis the temperature as determ-
ined from the simulation;ρis the density:ρ =N/V. All
values are in reduced units.
programming exercise
Write a program that simulates the behaviour of a Lennard–Jones liquid with
the proper argon parameters given above.
Check 1To check the program, you can use small particle numbers, such as 32 or
- Check whether the program is time-reversible by integrating for some time
(without rescaling) and then reversing velocities. The system should then return
to its initial configuration (graphical display of the system might be helpful).
Check 2The definite check is to compare your results for argon with literature.
A good value for the equilibration time is 10.0τand rescalings could take place
after every 10 or 20 time steps. A sufficiently long simulation time to obtain
accurate results is 20.0τ (remember the time step is 0.004τ). InTable 8.1
you can find a few values for the potential energy and pressure for different
temperatures. Note that the average temperature in your simulation will not be
precisely equal to the desired value. InFigure 7.1, the pair correlation function
forρ=N/V=1.06 andT=0.827 is shown.
It is interesting to study the specific heat (Eq. (7.37)) in the solid and in the gas
phase. You may compare the behaviour with that of an ideal gas,cV = 3 kB/T
per particle, and for a harmonic solid,cv= 3 kBTper particle (this is the Dulong–
Petit law).
Note that phase transitions are difficult to locate, as there is a strong hysteresis in
the physical quantities there. It is however interesting to obtain information about
the different phases. ForT=1,ρ=0.8 the argon Lennard–Jones system is found
in the liquid phase, and forρ =1.2 andT =0.5 in the solid phase. The gas
phase is found for example withρ=0.3 and andT=3.0. It is very instructive
to plot the correlation function for the three phases and explain how they look.
Another interesting exercise is to calculate the diffusion constant by plotting the
displacement as a function of time averaged over all particles. For times smaller