Description
1. Simulating Argon as a Lennard–Jones fluid (42 points)
Download the original paper on using the Verlet algorithm to simulate the molecular dynamics of Argon from Stud.IP ▶ Ubung: Methods of Computational Physic¨ s ▶ Files ▶ Additional Material. Read it and implement the described simulation with the setup, parameters and reduced units given. The starting temperature and the density ρ∗ should be adjustable parameters.
For now, ignore the optimisation of keeping neighbourhood tables suggested by Verlet. This will be the task in part b).
Instead of the original Verlet algorithm that is described in the publication, we suggest to use the Velocity Verlet algorithm in the following form:
v˜ , (1)
r(t + h) = r(t) + hv˜(t), (2) v . (3)
The new force F(t + h) is calculated between step (2) and (3). You can validate your integration by running the simulation for a number of integration steps then reversing the velocities, and running the simulation again for the same number of steps. Your final positions should be the same as the initial positions, except for rounding errors.
a) Run a simulation with 38 and ρ∗ = 0.55 for a total time of t∗ = 1, and plot the melting factor ρk against time. For the initial fcc lattice configuration, you should get ρk = 3N. Furthermore, plot the distribution of each velocity component at the end of the simulation and discuss your result. (6 points)
1
b) Implement the neighbourhood tables as suggested by Verlet and report the speed-up you can achieve. Use the suggested parameters, i.e. 3 and n = 16. (4 points) c) Run a simulation with 38 and ρ∗ = 0.55 for a total time of t∗ = 48, and plot the temperature T∗ against time. When is the system equilibrated? Calculate the time average ⟨T∗⟩, ignoring measurements before equilibration. Discuss your results, in particular comparing . (6 points)
d) Implement a momentum rescaling, which kicks in occasionally during the equilibration (e.g. every 20 steps) to drive the temperature towards . Then repeat the simulation and analysis made in part c), again comparing ⟨T∗⟩ with . (6 points)
e) Now implement additional observables, such that you can characterise the phase of your simulated Argon. First, we need a histogram of the (reduced) particle pair distances at each time step (left edge: 0, right edge: 5, bin width ∆r∗ = 0.05). You can update it whenever you need to iterate over all particle pairs (e.g. to calculate the forces). At the end, divide your histogram by the number of updates, the number of particles, by the volume of the spherical shell (4πr∗2∆r∗), and by ρ∗. The result is a discretised version of g(r∗), the radial pair-correlation function, which is our first observable to characterise the system state. After the simulation, you can plot g(r∗) versus r∗.
The second is the mean-squared displacement ⟨r∗(t∗)2⟩N, where ri∗ is the distance travelled by the ith particle, and ⟨·⟩N denotes the average over all particles. The third observable is the melting factor ρk(t∗). The displacement and the melting factor should be plotted against time.
Argon should be solid for T∗ = 0.5 and ρ∗ = 1.2, liquid for T∗ = 1.0 and ρ∗ = 0.8, and in its gas phase for T∗ = 3.0 and ρ∗ = 0.3.
Form hypotheses how the three observables might look like in the three phases, then run the corresponding simulations and discuss your findings.
(16 points)
f) Lastly, let us measure the pressure of the system, which can be calculated using
. (4)
The last term is a correction factor that takes into account that we cut off the potential at rν in the simulation. Since we did not measure g(r) beyond rν, you can set g(r) = 1 (which is its asymptotic value when the system is only locally ordered). Can you reproduce the reference value p/(ρkBT) = 1.31(2) for ρ∗ = 0.8 and T∗ = 1.0?
(4 points)
2
