Computational Techniques in Theoretical Physics
Section 15
Lennard-Jones Fluid


 Key Point
The Physics of Lennard-Jones fluid:
Lennard-Jones fluid

Lennard-Jones fluid is a system of many point particles of equal mass with an atom-atom pair interaction of the following form:

V(r) = V0 [ (r0/r)12 - 2(r0/r)6 ]
where r is the distance between two particles.


Properties of Lennard-Jones potantial:

r=r0:    V(r)=-V0; maximum attractive force.

r=infinity:   V(r)=0;   zero force.

r< r0/21/6 : V(r) > 0; repulsive.

The attractive part with a -1/r6 dependence is called dispersive force (or van der Waals force) which is due to induced dipole moments of electron clouds of two atoms.

The repulsive part is due to Pauli exclusion principle of electrons, which forbids two electrons to occupy the same state. The atoms behave as if they are solid balls of diameter r=r0/21/6.

The actual form of the repulsive force is not exactly 1/r12, but it is a good approximation and is also computationally more convenient.

Lennard-Jones system is a relatively simple system to study by molecular dynamics. The principle to simulation other systems, e.g., biological molecules, is the same except that the forces between atoms become more complicated.

Forces and equations of motion:

The force acting on the second particle located at r2 due to the presence of the first particle located at r1 is

F12 = - dV(r)/dr r/r,    r= r2-r1

where r = |r| is the distance and r is vector pointing from r1 to r2. The vector notation takes care of the direction of the force.

This form of force satisfies Newton's third law---the force acting on the first particle by the second particle is:

F12 = - F21

The Newton's equation for the i-th particle is then:

m d2ri / dt2 = sum j=1 to N Fij

Lennard-Jones fluid in a box.
In a computer simulation, we usually put the particles in a box of size Ld with periodic boundary condition.

We can view the periodic boundary condition as the box repeated infinite number of times in all directions, so that each particle located at the position ri has its images at:

ri+ lLi + mLj + nLk
where l, m, and n are integers, i, j, and k are mutually perpendicular unit vectors of the Cartesian coordinate axes.

Alternatively, we can put the system on a (three-dimensional) torus, but still use the Euclidean distance. That is, the distance between two particles is:

r = sqrt{ dx2 + dy2 + dz2 }
where dx is the x-component of the smaller one of |x2 - x1| and L - |x2-x1|. The shorter one when going from x1 to x2 (in x-direction) on a ring is used.

Likewise, dy and dz are similarly defined.

Care must be taken when calculating the direction of force or the vector r2 - r1 with periodic boundary condition.

Conservation laws
Total momentum
From Newtonian mechanics, we know that such a system has several conservation laws. The total momentum is conserved if there is no external force. This is the case if we simulation the system in a box with periodic boundary conditions. We should expect:
P = sum i=1 to N mvi
    = sum i=1 to N m dri/dt
to be a constant.

Total angular momentum

The total angular momentum,
L = sum i=1 to N m ri x vi
is also conserved if there is no external torque. Unfortunately, this quantity is not conserved if we consider the system in a torus, why?

Total energy

The most important conserved quantity in the simulation is the total energy:

E = K + V
    =   sum i=1 to N 1/2 m vi2
      + sum i=1 to N-1 sum j=i+1 to N V(|ri - rj| )

This equation can be used to test your algorithm. Deviations from energy conservation will occur in simulation due to discretization and roundoff error.

Comment on simulation algorithms:
Note that both the Verlet algorithm and symplectic algorithms do not conserve energy exactly, even though the original Newton's differential equations respect energy conservation. However, there exists quantity which are exactly conserved, and which differs from the total energy by a higher order term in h.

By the way, the total momentum conservation is obeyed exactly for the algorithms (why?).

Statistical mechanics of Lennard-Jones fluid
Why statistical mechanics?

The objective of simulaiton of fluids is to derive its average macrosopic physical quantities such as the relationship between pressure, volume and temperature.

Statistical mechanics is the theory that bridges the microscopic random motions with equilibrium and nonequilibrium macroscopic behavior.

In the molecular dynamics simulation scheme outlined above, we have implied that the system has a fixed number of particles N in a box of volume V=Ld. Also the dynamics is such that the total energy E is a constant fixed at some value.

In statistical mechanics, systems with fixed N, V, and E form so-called microcanonical ensemble. Our molecular dynamics simulation is thus a microcanonical one.

It is possible to have canonical molecular dynamics, where the temperature T, number of particle N, and the volume are constants. We'll study those methods later on.

Micro-Macro connection
The prefix micro- refers to something that is very small, e.g., microscopic, micro-computer, and microfilm.

The prefix macro- is the opposite, e.g., macroscopic, macrocosm.

The ancient greek natural philosophy is to interpret all macroscopic phenomena by the smallest constitutes of matter---atoms. This really becomes possible only in modern times.

Back to our molecular dynamics. It is not very useful if you just print out the positions of thousand particles at each time step. These `microscopic' quantities have to be somehow analyzed in some way so that we can compare them with real experiments. The theory for this micro-macro connection is statistical mechanics.

One of the most important features of physical system (e.g., our molecular dynamical system) is that it tends to reach a state that the macroscopic quantities do not vary with time if we leave the system undisturbed. This state is called an equilibrium state.

Statistical mechanics studies mostly equilibrium systems. If we start our molecular system with some arbitrary initial condition, it will not be in thermodynamic equilibrium. However, if we let the system develop according to Newton's law of motion, it will get into equilibrium after sometime.

Statistical mechanics has many predictions about the equilibrium properties which we'll just quote. Read your StatMech text book for more information.


Microcanonical distribution (fundamental distribution)


For system at fixed N, V, and E, the probability density of finding the system in the phase space:

X=(p,q)= (p1,p2, ... , q1,q2, ... )
is a constant on the energy surface, and zero outside the energy surface.

In other words, the phase space point X can be anywhere equally likely on the constant energy hypersurface. The probability being in a particular region is proportional to the volume of the region.

We can view the above statement as the basic assumption of statistical mechanics. Its validity depends on two important properties of the system:

Ergodicity and mixing.
Complex systems like our Jennard-Jones particles in a box have these properties, while a simple system like the harmonic oscillator does not.

Ergodicity means that for almost all starting point X0 in phase space, following the dynamics of the system, it is possible to reach a neighborhood of any other point X.

This allows us to replace the trajectory time average by phase space average or vice versa. This forms the basic link between statistical-mechanical theory and molecular dynamics simulation. Mixing means that the phase space points spread uniformly in the long time limit.