next up previous

4.1 Hard Spheres / Hard disks

Set up the hard spheres/disks on a lattice; then assign random initial velocities to the $N$ particles. Adjust these velocities such that 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, $T$ will continue to fluctuate.

Now calculate, for each pair of particles $(i,j)$ in the system, the time $t_{ij}$ it would take that pair to meet:

t_{ij}=\frac{-b - \sqrt{b^{2}-v^{2}(r^{2}-d^{2})}}{v^{2}}

where $d$ is the sphere diameter, $r$ is the distance between the centers of $i$ and $j$, and
$\displaystyle b$ $\textstyle =$ $\displaystyle (\mbox{$\bf r$}_{j}-\mbox{$\bf r$}_{i}) \cdot (\mbox{$\bf v$}_{j}-\mbox{$\bf v$}_{i})$  
$\displaystyle v$ $\textstyle =$ $\displaystyle \left\vert (\mbox{$\bf v$}_{j}-\mbox{$\bf v$}_{i})\right\vert$  

Note that the above formula for $t_{ij}$ derives from a straightforward solution of the quadratic equation $r_{ij}^{2}(t)-d^{2}=0$. As discussed in the appendix, this solution may give rise to numerical errors. A more secure alternative is, with the same meaning for $b$ and $v$,

q = - \left[ b +sgn(b) \sqrt{b^{2}-v^{2}(r^{2}-d^{2})} \right],
min \left\{ q/v^{2},(r^{2}-d^{2})/q \right\}

(Since $b$ must be negative for a pair to meet in the future, we may restrict the calculation to the case $sgn(b)<0$.)

Set up two arrays that contain, for each particle $i$, the smallest positive collision time $t(i)=min(t_{ij})$ and the next collision partner $j(i)$. If a particle has no collision partner at positive times we simply set $j(i)=0$ and $t(i)=[\infty]$, i.e. the largest representable number.

This double loop over all $N(N-1)/2$ 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 $t(i_{0})$. This is the time until the very next encounter between two particles in the system. The indices of those two particles are named $i_{0}$ and $j_{0}$.

During the time $t(i_{0})$ all particles perform a free flight, thus:

\mbox{$\bf r$}_{i}\longrightarrow \mbox{$\bf r$}_{i}+\mbox{$\bf v$}_{i} \cdot t(i_{0})

and all $t(i) \longrightarrow t(i)-t(i_{0})$.

Now an elastic collision between $i=i_{0}$ and $j=j_{0}$ occurs, resulting in the new velocities

\mbox{$\bf v$}_{i}'=\mbox{$\bf v$}_{i} + \Delta \mbox{$\bf v...
...\mbox{$\bf v$}_{j}'=\mbox{$\bf v$}_{j} - \Delta \mbox{$\bf v$}


\Delta \mbox{$\bf v$}=b  \frac{\mbox{$\bf r$}_{ij}}{d^{2}}

Since $i_{0}$ and $j_{0}$ have new flight directions and speeds, all pairwise collision times $t_{ij}$ involving these two must be recalculated. This means that $2N-3$ pairs have to be scanned.

This completes the basic hard sphere/disc MD step. Now all $t_{i}$ are once more searched for the smallest element, etc.
Figure 4.1: Molecular dynamics of hard spheres
{\bf Molecular dynamics simulation of hard sphere...
$\Delta x_{\alpha} \equiv v_{\alpha}  \Delta t \leq L/2-d$.

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 $t(i)$ and $j(i)$; 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.

next up previous
F. J. Vesely / University of Vienna