Franz J. Vesely > MolSim Tutorial > Molecular Dynamics II  
 
 





 
< >
Vesely MolSim Tutorial







3. Molecular Dynamics Simulation: Standard (NVE)

Fundamentally different methods for () Hard bodies with impulsive collisons and () Particles interacting via smooth potential forces.

Sections

3.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
  • the total kinetic energy is consistent with some desired temperature: $E_{k}=3NkT/2$
  • 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, $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{\textstyle -b - \sqrt{b^{2}-v^{2}(\textstyle r^{2}-d^{2})}}{v^{2}} $

where $d$ is the sphere diameter, $r$ is the distance between the centers of $i$ and $j$, and

$ \begin{eqnarray} b &=& ( r_{j}- r_{i}) \cdot ( v_{j}- v_{i}) \\ v & = & \left| ( v_{j}- v_{i})\right| \end{eqnarray} $


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], \;\;\;\;t_{ij}= 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:

$ r_{i}\longrightarrow r_{i}+ 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

$ v_{i}'= v_{i} + \Delta v , \;\;\; v_{j}'= v_{j} - \Delta v $

where

$ \Delta v=b \frac{\textstyle r_{ij}}{\textstyle 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 3.1: Molecular dynamics of hard spheres

Hard Spheres MD:

Immediately after a collision, for each particle $i$ in the system the time $t(i)$ to its next collision and the partner $j(i)$ at that collision are assumed to be known.
  • Determine the smallest positive element $t(i_{0})$ among the $t(i)$, identify the corresponding particle $i_{0}$ and its collision partner $j_{0}\equiv j(i_{0})$.

  • Let all particles follow their free flight paths for a period $\Delta t \equiv t(i_{0})$; subtract $\Delta t$ from each $t(i)$.

  • Perform the elastic collision between $i_{0}$ and $j_{0}$; after the collision these spheres have the new velocities $ v' = v \pm \Delta v , \;\;\;{\rm with} \;\; \Delta v =b \frac{\textstyle r_{ij}}{\textstyle d^{2}} $

  • Recalculate all times $t(i)$ that involve either $i_{0}$ or $j_{0}$, i.e. for $i=i_{0}$, $i=j(i_{0})$, $i=j_{0}$, and $i=j(j_{0})$.

  • Go to (1).
At low densities the large free path may create problems with the periodic boundary conditions, some particle suddenly appearing where it overlaps another. One therefore limits the time allowed for free flight such that for each particle and each coordinate $\alpha$ the free flight displacement fulfills $\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.



3.2 Continuous Potentials

For continously varying pair potentials we have for a particle $i$ at any time $t$

$ \ddot{ r}_{i}(t)=\frac{\textstyle 1}{\textstyle m} \sum_{j\neq i} K_{ij}(t)\;\;\;\;{\rm with} \;\;\;\; K_{ij}\equiv -\nabla_{i}u(r_{ij}) $


Considering the Lennard-Jones interaction, we find for the pair force

$ K_{ij}=-24 \frac{\epsilon}{\sigma^{2}} \left[ 2 \left(\frac{\textstyle r_{ij}}{\textstyle \sigma}\right)^{-14}\right. -\left. \left(\frac{\textstyle r_{ij}}{\textstyle \sigma}\right)^{-8}\right] r_{ij} $

where $ r_{ij} \equiv r_{j}- r_{i}$.

The above-mentioned nearest image convention (NIC) is used in the evaluation of the force acting on a particle.

Having determined this total force, the equation of motion for particle $i$ may be numerically integrated. A widely used technique is Verlet's algorithm

$ r_{i}(t_{n+1})= 2 r_{i}(t_{n})- r_{i}(t_{n-1}) + b_{i}(t_{n})(\Delta t)^{2} $

(with $ b_{i}\equiv \sum_{j\neq i} K_{ij}/m$).


PROJECT MD (LENNARD-JONES): Augment the subroutine module ENERGY such that it computes, for each Lennard-Jones particle $i$ in the system, the total force exerted on it by all other particles $j$: $ K_{i} \equiv \sum_{j \neq i} K_{ij}$, with $ K_{ij}$ as given above; remember to apply the nearest image convention.

Write a subroutine MOVE to integrate the equations of motion by a suitable algorithm such as Verlet's. Having advanced each particle for one time step, do not forget to apply periodic boundary conditions to retain them all in the simulation box.

Write a main routine that puts the subroutines STARTCONF, ENERGY and MOVE to work. Test your first MD code by monitoring the mechanically conserved quantities.

Do a number of MD steps - say, $50$-$100$ - and average the quantity $| v^{*}|^{2}$ to estimate the actual temperature. To adjust the temperature to a desired value, scale all velocity components of all particles in a suitable way. Repeat this procedure up to 10 times. After $500$- $100$ steps the fluid will normally be well randomized in space, and the temperature will be steady - though fluctuating slightly.



vesely may-06

< >