Computational Physics

8.3 A molecular dynamics simulation program for argon 209

each of the momenta in order to make the total momentum zero. Now the kinetic
energy is calculated and then all momenta are rescaled to arrive at the desired kinetic
energy.
When calculating the forces, the minimum image convention should be adopted.
It is advisable to start without using a neighbour list. For the minimum image
convention it should be checked for each pair(i,j)whether the difference of thex-
componentsxi−xjis larger or smaller thanL/2 in absolute value. If it is larger, then
an amountLshould be added to or subtracted from this difference to transform it to
a value which is smaller thanL/2 (in absolute value). In many codes, this translation
is implemented as follows:

x→x−[x/L]∗L, (8.21)

where[]denotes the integer part. This procedure is then repeated for they- and
z-components. Potential and force may be adjusted according toEqs (8.12)and
(8.13).
The equations of motion are solved using the leap-frog or the velocity form of
the Verlet algorithm. A good value for the time step is 10−^14 s which in units of
(mσ^2 /ε)^1 /^2 is equal to about 0.004. Using the argon mass as the unit of mass,σas
the unit of distance andτ=(mσ^2 /ε)^1 /^2 as the unit of time, thex-component of the
force acting on particleiresulting from the interaction with particlejis given by

Fijx=(xi−xj)( 48 rij−^14 − 24 r−ij^8 ) (8.22)

with similar expressions for they- andz-components.
After each step in the Verlet/leap-frog algorithm, each particle should be checked
to see whether it has left the volume. If this is the case, it should be translated over
a distance±Lalong one or more of the Cartesian axes in order to bring it back into
the system in accordance with the periodic boundary conditions.
During equilibration, the velocities (momenta) should be rescaled at regular
intervals. The user might specify the duration of this phase and the interval between
momentum rescalings.
During the production phase, the following quantities should be stored in a file
at each time step: the kinetic energy, potential energy, and the virial
∑

ij

rijF(rij). (8.23)

Furthermore, the program should keep a histogram-array containing the numbers
of pairs found with a separation betweenrandr+for, say,=L/200 from
which in the end the correlation function can be read off.