Key Point

• The statistical mechanical description of Lennard-Jones fluid.
• Computer algorithm.

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 v² / (k_B T) }

where k_B is Boltzmann constant, and T is temperature.

The proportionality constant is determined by the normalization condition:
 

int P(v) d³ v  = 1

We can also consider the distribution of the speed v = |v| = sqrt{ v_x² +  v_y²  + v_z² }, which is (in three dimensions):
 

P(v) dv  = exp{-m v²/2k_BT} dv
 
 

Equipartition theorem

 

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

< 1/2  m_i v_ix² > =
< 1/2 m_i v_iy² >  =
< 1/2 m_i v_iz² >  =
1/3 < 1/2 m_i v_i² > = 1/2 k_B T

where  v_ix 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

 
Definition:

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

Determination of the equation:

We can calculate the pressure from:
 

PV = N k_B T +
         1/3 < sum _i<j ( r_i - r_j ) _* F_ij >

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) d³r  is the average number of particles in the
volume element dV=d³r, 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) = < N_r > / d³r

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 N²/V  int _{0 to infinity} V(r) g(r) r² dr.
 
 
 

 PV = NkB T

 - 2piN²/3 V   int _{0 to infinity}  dV(r)/dr  g(r)r³dr.
 
 
 

 

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) = V₀ [ (r₀/r)¹² - 2(r₀/r)⁶ ]

The length can be measured in terms of the parameter r₀.  The scale for energy is V₀.

We measure temperature in unit of  V₀ /k_B so that the combination k_BT/V₀ is a dimensionless number.
 

We use the combination sqrt{ m r₀²/ 48 V₀} 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 =  r_i / r₀

t' = t _* sqrt{ m r₀²/ 48 V₀}

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:

V_wall( x ) = a/x⁹
 

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:
 

• the force acting on particle i from particle  j
• the other way around.

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:
 

• q[0] is the x-component of the 1st particle;
• q[1] is the y-component of the 1st particle;
• q[2] is the x-component of the 2nd particle,
• q[i]  similar for other quantities.

The computation can be conducted in two steps:
 

• Caculation of forces
• Numerical integration of Newton's equations.

 
 
 

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 m_i = 1):
 
 

p_n+1/2 = p_n + h/2 f(q_n)

q_n+1 = q_n + h p_n+1/2

p_n+1 = p_n+1/2 + h/2 f(q_n+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];
   }
   force(q,f,n);
   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.
 
 
 
 

Homework