Comments on PHY480 HW#5
Aleksandar Donev and David Tomanek

Theoretical background

The atoms in this cluster interact via a central potential:

The force is then radial and is given via the gradient of the potential:

In the cluster of atoms though, we have to sum the force on each atom over all the other atoms:

where now the radius vector is in fact the radius vector between the pair of interacting atoms:

The equation of motion for each atom is then Newtons second law:

This is however a second-order ODE, and Runge-Kutta only solves first order ODE's. So we transform this into a system of 2 first order equations via using the velocity as an intermediate result. This way we get a system that resembles the form we studied in class:

but now the vector Y is composed from 6N components, where N is the number of atoms, and we have 3 coordinates and 3 velocities per atom:

These are the equations of motion ready to be used in the Fortran program.

In this molecular dynamics simulation we are studying the time evolution of a cluster of Lenard-Jones atoms which are initially at rest and have given positions (in some units). Therefore, the only important input parameter is the size of the atoms sigma (the other parameters just change the time scale).

There are two major types of clusters depending on the initial conditions:

-Bound clusters result when the total initial energy is negative. In this case, the atoms can never separate to infinity, but must stay bound at all times. The r.m.s. radius of the cluster then oscillates between some finite sizes.

-Unbound clusters result when the total initial energy is positive (the atoms are compressed too much). In this case, the initial potential energy gets quickly converted into kinetic energy and the atoms of the cluster fly apart. The r.m.s. radius of the cluster then increases without bound.

So everything depends upon the initial potential (since the kinetic is zero) energy of the cluster. Here is a plot that shows how this energy depends on the value of sigma for the input file Na3.xyz. The plot shows that for values of sigma up to about 1.4 the cluster is bound, and for larger sizes it is unbound.

Here are the results of the ljmd.f program for a bound cluster, demonstrating all of the above claims. This is a plot that shows the energy conversion in time, this is a plot that shows the r.m.s. radius in time, and this is a plot that shows the orbits of the particles in time.

Here are the results of the program for an unbound cluster demonstrating all of the above claims. This is a plot that shows the energy conversion in time, this is a plot that shows the r.m.s. radius in time, and this is a plot that shows the orbits of the particles in time.

Looking at the output of the ljmd program, we see that the mechanically conserved constants, namely the total energy, linear and angular momentum, the center of mass, etc., are all constant within some tolerance. If we increase the time step size h, we see that the variation in the conserved quantities increases, and for some step sizes the numerical integration becomes unstable.

The way the step size is usually determined in practice is by first deciding what kind of relative accuracy over the entire run (say 5%) we wish to achieve, and then plotting some conserved quantity, usually the total energy, versus time. We then increase the step-size slowly until it becomes too large and the energy deviations become larger than the desired accuracy. Then, we freeze the time step at that value.