Computational Techniques in Theoretical Physics

Section 17
Lennard-Jones Fluid, Part II

 Key Point

Statistical Mechanics of Lennard-Jones fluid:
Maxwell velocity distribution
If we look at a particular particle, its velocity  changes at random.

In fact, the values of velocity obey some probability distribution.

Equivalently, if we look at the velocities of all the particles at a given time, they should form the same distribution patern.

According to laws of physics, the velocities obeys the Maxwell distribution:

 P(v) ~ exp{ - 1/2  m v2 / (kB T) }
where kB is Boltzmann constant, and T is temperature.

The proportionality constant is determined by the normalization condition:

int P(v) d3 = 1
We can also consider the distribution of the speed v = |v| = sqrt{ vx2 +  vy + vz2 }, which is (in three dimensions):
P(v) dv  = exp{-m v2/2kBT} dv
Equipartition theorem

Based on the Maxwell distribution, the average kinetic energy in each degree of freedom obeys the equipartition theorem:

< 1/2  mi vix2 > =
< 1/2 mi viy2 >  =
< 1/2 mi viz2 >  =
1/3 < 1/2 mi vi2 > = 1/2 kB T
where  vix is the x-component of the velocity of the i-th particle. The other two are the y- and z- component respectively.

This equation can be conveniently used to calculate the temperature of the system.

Equation of state

The equation of state is a relation between volume, temperature, and pressure.

Determination of the equation:

We can calculate the pressure from:

PV = N kB T +
         1/3 < sum i<j ( ri - rj ) * Fij >

where the summation runs over each pair of interactions exactly once [N(N-1)/2 terms]; each term is the scalar (dot) product of the pair interaction force and the distance vector. The above relationship is known as virial theorem.

Pair correlation function g(r)

This is an important function which characterizes how two particles are related (the magnitude of interaction between each other) in the system.

Given a particle at some position, g(r) d3r  is the average number of particles in the
volume element dV=d3r, located at a distance smaller than |r| away from the given particle.

This function can be easily calculated
from molecular dynamics, by taken a statistics of number of particles within a certain distance |r| from any given particle:

g(r) = < Nr > / d3r

If we know the pair correlation function, the
average total potential energy and equation of state can be calculated:

< V > = < sum i<j V( | ri - rj | ) >
= 2pi N2/V  int 0 to infinity V(r) g(r) r2 dr.
 PV = NkB T
 - 2piN2/3 V   int 0 to infinity  dV(r)/dr  g(r)r3dr.
Computer implementation  of molecular dynamics.
Basic strategy:
In any numerical calculation, it is best to reduce quantities to dimensionless numbers.
This way we avoid round-off errors etc.

For the Lennard-Jones system, we have
natural scales from the potential:

V(r) = V0 [ (r0/r)12 - 2(r0/r)6 ]

The length can be measured in terms of the parameter r0.  The scale for energy is V0.

We measure temperature in unit of  V0 /kB so that the combination kBT/V0 is a dimensionless number.

We use the combination sqrt{ m r02/ 48 V0} as the unit of time, which is proportional to
the vibrational period of a two-body Lennard-Jones particle.

In terms of the dimensionless variables, the Hamilton's equation for the Lennard-Jones system is:

dp'i/dt' = sum j =/= i (r'i-r'j) [r'ij-14 - 1/2 r'ij-8]

dr'i/dt' = p'i = v'i

where the reduced coordinates and time are defined as:

r'i ri / r0

t' = t * sqrt{ m r02/ 48 V0}

r'ij = | r'i-r'j |

For each particle we need to calculate the total force acting on it by other particles.  In the following particular implementation,
we will not use periodic boundary condition. Rather, we put the particles inside a box with repulsive hard wall.  Each particle, besides interacting with every other particles, also interacts with the walls, by the form of wall potential:

Vwall( x ) = a/x9

where x is the perpendicular distance to the wall.

Clearly, the force quickly decreases to zero away from the wall.

In two dimensions, we have four terms for each of the four walls.

This is purely for a better, more realistic graphic display.  It is not done in usual molecular dynamics simulation which typically use periodic boundary conditions.  It is not difficult to modify the program to adopt a periodic boundary condition.

After we calculated the forces due to walls for each particle, we go through each pair of particles (i,j) exactly once.  There are
two forces associated with each pair of particles:

They are equal in magnitude and opposite in sign.  We accumulate this force for particle i (due to particle j) and for particle j (due to particle i).

Computer algorithm:
 To make the program very concise, we use a one-dimensional arrays for positions, momenta (same as velocity when mass is taken to be 1), and forces.  Thus:

The computation can be conducted in two steps:


The C-code for computing the force acting on a particle in Lennard-Jones fluid:

The following function calculates the force
f[ ]  on each particle when the particles are at locations specified by q[ ].

int L;

void force( double q[], double f[], int n)
   double xl, xr, r2, LJforce, dx, dy;
   int  i, j;
   for (i = 0; i < n; ++i) {
      xl = q[i];
      xr = L - q[i];
      f[i] = pow(xl, -9) - pow(xr, -9);

   for (i = 1; i < n/2; ++i)
      for (j = 0; j < i; ++j) {
         dx = q[2*j  ] - q[2*i  ];
         dy = q[2*j+1] - q[2*i+1];
         r2 = dx*dx + dy*dy;
         LJforce =
        ( pow(r2,-7) - 0.5*pow(r2,-4) );
         f[2*i  ] = f[2*i  ] - LJforce*dx;
         f[2*i+1] = f[2*i+1] - LJforce*dy;
         f[2*j  ] = f[2*j  ] + LJforce*dx;
         f[2*j+1] = f[2*j+1] + LJforce*dy;

The use of symplectic algorithm:

The application of the symplectic algorithm is straightforward. The function takes the current values of positions, momenta, and forces, then derive the values at time elapsed by h or h/2.

Note that we do not need to calculate the forces twice.  We rewrite the algorithm C in the following equivalent way (assuming mass mi = 1):

pn+1/2 = pn + h/2 f(qn)

qn+1 = qn + h pn+1/2

pn+1 = pn+1/2 + h/2 f(qn+1)

In the program, before this function is called, one must call the  force() function just once to initialized the force f[ ].

The symplectic algorithm C-code

#define nmax 200

void symplectic(double p[], double q[], double f[], int n)
   const double h = 0.06;
   double t[nmax];
   int i;

   for (i = 0; i < n; ++i) {
      t[i] = p[i] + 0.5 * h * f[i];
      q[i] = q[i] + h * t[i];
   for (i = 0; i < n; ++i)
      p[i] = t[i] + 0.5 * h * f[i];

We have not discussed how to choose the step size h.  Clearly, h can not be too large, otherwise, our numerical trajectory will
deviate greatly from the true path.  But more importantly, the algorithm becomes unstable for large value of h.  As a rule of thumb, the value h should be about 1/10 or 1/20 of the smallest vibrational period of the system.  For the Lennard-Jones particle system, this means that the dimensionless h should be about 0.1
or smaller.