   # 4.1 Hard Spheres / Hard disks

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 temperature: • 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 representable number.

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: and all .

Now an elastic collision between and occurs, resulting in the new velocities where 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. EXERCISE: 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'' simulation.   F. J. Vesely / University of Vienna