Set up the hard spheres/disks on a lattice; then assign
random initial velocities to the particles.
Adjust these velocities such that
the total kinetic energy is consistent with some desired
the total momentum (conserved in the simulation) equals zero.
Note that the adjustment of the temperature is only temporary; it will
have to be repeated several times before thermal equilibrium is reached;
even then, will continue to fluctuate.
Now calculate, for each pair of particles in the
system, the time it would take that pair to meet:
where is the sphere diameter, is the distance between the centers of
and , and
Note that the above formula for derives from a straightforward
solution of the quadratic equation
. As discussed in
the appendix, this solution may give rise to numerical errors. A more
secure alternative is, with the same meaning for and ,
(Since must be negative for a pair to meet in the future, we
may restrict the calculation to the case .)
Set up two arrays that contain, for each particle , the smallest
positive collision time
and the next collision
partner . If a particle has no collision partner at positive times
we simply set and , i.e. the largest
This double loop over all pairs need be performed only once,
at the start of the simulation.
Now find the smallest element in the table of free flight times,
calling it . This is the time until the very next encounter
between two particles in the system. The indices of those two particles
are named and .
During the time all particles perform a free flight, thus:
Now an elastic collision between and occurs,
resulting in the new velocities
Since and have new flight directions and speeds,
all pairwise collision times involving these two must be
recalculated. This means that pairs have to be scanned.
This completes the basic hard sphere/disc MD step. Now all
are once more searched for the smallest element, etc.
Molecular dynamics of hard spheres
For a two-dimensional system of hard disks, write subroutines to
a) set up an initial configuration (simplest, though not best: square
lattice;) b) calculate and ; c) perform a pair collision.
Combine these subroutines into an MD code. To avoid the difficulty
mentioned at the end of the preceding figure, one might simply use
reflecting boundary conditions, doing a ``billiard dynamics''
F. J. Vesely / University of Vienna