Numerical Integration

In a computer experiment we will not be able to generate the true trajectory of a system with a given set of initial positions and velocities. For all potentials $U$ used in real applications only numerical integration techniques can be applied. These techniques are based on a discretization of time and a repeated calculation of the forces on the particles. Many such methods have been devised [23] (look for "Integration of Ordinary Differential Equations"). However, what we are looking for is a method with special properties: long time energy conservation and short time reversibility. It turns out that symplectic methods (they conserve the phase space measure) do have these properties. Long time energy conservation ensures that we stay on (in fact close) to the constant energy hypersurface and the short time reversibility means that the discretized equations still exhibit the time reversible symmetry of the original differential equations. Using these methods the numerical trajectory will immediately diverge from the true trajectory (the divergence is exponential) but as they stay on the correct hypersurface they still sample the same microcanonical ensemble. On the other hand, a short time accurate method will manage to stay close to the true trajectory for a longer time and ultimately will also exponentially diverge but will not stay close to the correct energy hypersurface and therefore will not give the correct ensemble averages.

Our method of choice is the velocity Verlet algorithm [24,25]. It has the advantage that it uses as basic variables positions and velocities at the same time instant t. The velocity Verlet algorithm looks like a Taylor expansion for the coordinates:

$${\bf R}(t+\Delta t) = {\bf R}(t) + {\bf V}(t) \Delta t + {{\bf F}(t) \over 2 M } \Delta t^2 \enspace .$$ (9)

This equation is combined with the update for the velocities

$${\bf V}(t+\Delta t) = {\bf V}(t) + {{\bf F}(t + \Delta t) + {\bf F}(t) \over 2 M } \Delta t^2 \enspace .$$ (10)

The velocity Verlet algorithm can easily be cast into a symmetric update procedure that looks in pseudo code

       V(:) := V(:) + dt/(2*M(:))*F(:)
       R(:) := R(:) + dt*V(:)
       Calculate new forces F(:)
       V(:) := V(:) + dt/(2*M(:))*F(:)

To perform a computer experiment the initial values for positions and velocities have to be chosen together with an appropriate time step (discretization length) $\Delta t$ . The choice of $\Delta t$ will be discussed in more detail in a later chapter about ab initio molecular dynamics. The first part of the simulation is the equilibration phase in which strong fluctuation may occur. Once all important quantities are sufficiently equilibrated, the actual simulation is performed. Finally, observables are calculated from the trajectory. Some quantities that can easily be calculated are (for these and other quantities see the books by Frenkel and Smit [13] and Allen and Tildesley [12] and references therein)

• The average temperature

  $$\left\langle {1 \over 2} M {\bf V}^2 \right\rangle = {1 \over 2 } k_{\rm B} T$$ (11)

  
  

• The diffusion constant

  $$D = {1 \over 6 \tau} \Delta{\bf R}^2(\tau) = {1 \over 6 \tau} \langle ({\bf R}(\tau) - {\bf R}(0) )^2 \rangle$$ (12)

  
  
  for large times $\tau$ .
• The pair correlation function

  $$g(r) = {V \over N} \left\langle \sum_i \sum_{j \neq i} \delta({\bf r} - {\bf R}_{ij}) \right\rangle$$ (13)

  
  

• The temporal Fourier transform of the velocity autocorrelation function $\langle {\bf V}(t) \cdot {\bf V}(0) \rangle$ is proportional to the density of normal modes (in a purely harmonic system).