Lecture 2

Hard Spheres

27-534/27-734 ∙ Methods of Computational Materials Science

 Hard sphere model

‣ Our rst goal is to study a system of hard spheres.

Length L

These are the hard

spheres.
Hard spheres collide
against one another
Diameter 2σ

John Dalton (1766-1844)

These are the walls of

Hard spheres the box.
collide against
the wall

‣ Scienti cally, this model is important because it provides a simple

description of an ideal gas (as de ned in thermodynamics).
2 | Dabo
fi
fi
fi
e seen, the fluctuations
Hardthan
ansition region sphere
in model
and pure fiuid (A/Ae
larger in the middle of
n near the ‣ Historically,
ends. The this model is also important because it is one of the
rst models
dicate the middle range ever simulated on a computer
of the intervals show
ertical line and 4 show
Temperature over pressure T/P

he line. For the shorter Crystal

y this medium vertical
Progress of Theoretical Physics Supplement No. 178, 2009 1

In Memoriam:

ig. 1 indicate the esti-Fluid branch Thomas E. Wainwright

September 22, 1927 – November 27, 2007

Downloaded from [Link] by guest on 27 August 2020

Berni J. Alder

pressure determination
Lawrence Livermore National Laboratory, P.O. Box 808, Livermore,

Berni J. Alder (1925-2020)

CA 94551-9900, USA

are made by comparing Crystal Fluid

+Fluid
ons each with various
arison was typically
ities. The fact that the
Crystal branch
g starting conditions System of hard spheres undergoing
indication that phase a phase region, showing crystalline
mpled. Four diferent Volume V and fluid regions (the curves record
Thomas E. Wainwright
the trajectories of the hards spheres)
ed: (1) all the particles (1927-2007)
Thomas E. Wainwright (May 1972)

ith only one particle in ___

attice positions andAlder and Wainwright, ‘Phase Transition in Elastic Disks’, Physical Review 123, 359 (1962)
3 | Dabo
in Fic. 2. The traces of the centers of particles in the phase-
fi
 Kinetic energy

‣ The starting point of most models is to write the total energy

‣ There are two contributions to the total energy of a system of

classical particles: the kinetic energy K and potential energy U

‣ K is always de ned in the same way:

De nition | Kinetic energy

m1
v1 m2
Consider a system of N atoms with
velocities vi and masses mi. The kinetic
m3
energy K of this system is v2
N v3
X 1 vN
K= mi vi2
i=1
2
mN
Note that the K is always positive.

4 | Dabo
fi
fi
 Potential energy

‣ Takeaway. Expressing the kinetic energy K is relatively simple.

‣ But expressing the potential energy U is in general much more

complicated

‣ This is where it is very important to choose an appropriate model

‣ Fortunately, in the case of the hard sphere model (and of some

other models that we will see), U has relatively simple form

N X
X N This prime means that we
1 0 exclude i = j (a particle cannot
U= (rij ) interact with itself)
2 i=1 j=1
This one half is here because
each pair has been counted
twice (e.g., 1–7 is the same as Distance between i and j
7–1)
The function ϕ is called the pair interaction
5 | Dabo
 Pair interaction function

2σ rij
Pair interaction
ϕ(rij) Sphere i Sphere j

Question. Can you guess the shape

of the pair interaction ϕ as a function
of the pair distance rij?

rij
Pair distance
6 | Dabo rij<2σ rij=2σ rij>2σ
 Pair interaction function

2σ rij
Pair interaction Two spheres cannot have a
distance lower than 2σ.
ϕ(rij) Sphere i Sphere j
The energy is in nite.

+∞

When the spheres are far

apart they don’t see one
another

rij
Pair distance
7 | Dabo rij<2σ rij=2σ rij>2σ
fi
 Total‣ energy
There are two contributions to the total energy of a system of
classical particles: the kinetic energy K and potential energy U

‣ Summary. Fordefined
‣ K is always the hard sphere
in the model,
same way for anythe kinetic energy K,
system
potential energy U, and their sum, the total energy E, can be
Definition | Kinetic energy
expressed as m1
v 1 m2
Consider a system of N atoms with
velocities vi N
and masses mi. The kinetic
X1
energy K of this system is
m3
v2
2
K= Nmi v
X 1 i2
v3
vN
K =2 mi v i
i=1 2
i=1
mN
Note that the K is always positive. ϕ(rij)
N X
X N
7 | Dabo
1 0
U= (rij ) = 0
2 i=1 j=1
rij
___
E=K+U=K Note. The potential energy U is 0 except when
two particles collide. But collisions are so brief
that we can consider U to be always 0.
8 | Dabo We will discuss collisions in a few slides.
 Molecular dynamics

‣ At this point, we know how to write the energies for the hard
sphere model.
Length L

These are the hard

spheres.
Hard spheres collide
against one another
Diameter 2σ

These are the walls of

Hard spheres the box.
collide against
the wall

‣ Question. What do we need to do to simulate the dynamical

evolution of this system?
9 | Dabo
 Molecular dynamics

‣ Answer. We need to write and integrate the equations of motion,

i.e., Newton's laws.

dvi

{
m = Fi
dt
dri
= vi
dt Question. What are the forces in the
hard sphere model?

10 | Dabo
 Molecular dynamics

‣ Answer. We need to write and integrate the equations of motion,

i.e., Newton's laws.

dvi

{
m = 0Fi
dt
dri
= vi
dt
Answer. We have seen that the potential U
is zero (expect during the very brief
collisions, which will be discussed shortly).
Therefore, the forces which are derivatives
of the energy are also zero.

‣ Now, the last step before coding is to rewrite the problem into
vector form using a process called time discretization.
11 | Dabo
 Time discretization

‣ Time discretization is the process of converting the

continuous time variable into a discrete time variable so
that the time equations can be solved numerically

Continuous time t

Discrete time
t1 t2 … tn …
discrete time point Δt
time step

‣ To perform this discretization, we introduce the notations

ri,n = ri(tn) vi,n = vi(tn)

particle index time index
12 | Dabo
 Time discretization (cont’d)

‣ Using these notations, we can discretize the equations of

motions use the nite-difference approach
This dot means
time derivative . ri,n+1 – ri,n

{ {
ri = vi = vi,n
∆t
. vi,n+1 – vi,n
mvi = 0
time discretization m =0
∆t
to obtain

ri,n+1 = ri,n + ∆t vi,n

vi,n+1 = vi,n

‣ We are almost done. Let us plug these equations

into Matlab (or Python) and visualize the simulation.
13 | Dabo
fi
 Simulation

ri,n+1 = ri,n + ∆t vi,n

vi,n+1 = vi,n

‣ Download the code hard_spheres.m from CANVAS

‣ Read the code and understand how it works (ask questions or

use the help command or the dialog box, if needed)

‣ Complete the code with the above equations

‣ Question. What are the two missing items in this simulation?

14 | Dabo
 Simulation

ri,n+1 = ri,n + ∆t vi,n

vi,n+1 = vi,n

‣ Download the code hard_spheres.m from CANVAS

‣ Read the code and understand how it works (ask questions or

use the help command or the dialog box, if needed)

‣ Complete the code with the above equations

‣ Question.
Answer. First,
Whatweare
needthetotwo
code
missing
boundary
itemsconditions,
in this simulation?
i.e., what
happens when a sphere hits a wall. Second, we need to code
collisions, i.e., what happens when a sphere hits another sphere.

14 | Dabo
 Hard-wall boundary conditions

‣ Boundary conditions. Particle-wall collisions happen in 4 cases.

?
2σ
L
?
?

Question 1. What are

these 4 cases?
0

?
15 | Dabo
0 L
 Hard-wall boundary conditions

‣ Boundary conditions. Particle-wall collisions happen in 4 cases.

If y > L – σ
2σ
L
If x < σ
If x > L – σ

Question 2. What should

we do in each case?
0

If y < σ

16 | Dabo
0 L
 Hard-wall boundary conditions

‣ Boundary conditions. Particle-wall collisions happen in 4 cases.

If y > L – σ
vy → – v y 2σ
L
If x < σ
vx → – v x If x > L – σ
vx → – v x

Question 3. Can we
reduce the number of “If”
statements?
0

If y < σ
vy → – v y
17 | Dabo
0 L
 Hard-wall boundary conditions

‣ Boundary conditions. Particle-wall collisions happen in 4 cases.

2σ
L
If x < σ or x > L – σ
vx → – v x
If y < σ or y > L – σ
vy → – v y

Question 4. Can you

write these two
statements in Matlab?
0

18 | Dabo
0 L
 Hard-wall boundary conditions

‣ Examine this piece of code carefully

If x < σ or x > L – σ if r(i,1,n+1) < sigma | r(i,1,n+1) > L – sigma

v(i,1,n) = – v(i,1,n)
vx → – v x end

If y < σ or y > L – σ if r(i,2,n+1) < sigma | r(i,2,n+1) > L – sigma

vy → – v y v(i,2,n) = – v(i,2,n)
end

___
Note. The logical operation ‘or’ is denoted as a single pipe (|) or double pipe (||) in
MATLAB. The double pipe is the short-circuited version of the logical operation;
the second term is not evaluated if the rst term is true.
19 | Dabo 19
fi
 Hard-wall boundary conditions

‣ Examine this piece of code carefully

1 is the rst space dimension
i is the index of the particle
n +1 is the next time index

If x < σ or x > L – σ if r(i,1,n+1) < sigma | r(i,1,n+1) > L – sigma

v(i,1,n) = – v(i,1,n)
vx → – v x end
2 is the second
n is the current time index
space dimension
If y < σ or y > L – σ if r(i,2,n+1) < sigma | r(i,2,n+1) > L – sigma
vy → – v y v(i,2,n) = – v(i,2,n)
end

___
Note. The logical operation ‘or’ is denoted as a single pipe (|) or double pipe (||) in
MATLAB. The double pipe is the short-circuited version of the logical operation;
the second term is not evaluated if the rst term is true.
19 | Dabo 19
fi
fi
 Hard-wall boundary conditions

‣ Examine this piece of code carefully

1 is the rst space dimension
i is the index of the particle
n +1 is the next time index

If x < σ or x > L – σ if r(i,1,n+1) < sigma | r(i,1,n+1) > L – sigma

v(i,1,n) = – v(i,1,n)
vx → – v x end
2 is the second
n is the current time index
space dimension
If y < σ or y > L – σ if r(i,2,n+1) < sigma | r(i,2,n+1) > L – sigma
vy → – v y v(i,2,n) = – v(i,2,n)
end
‣ Coding. Add these lines to the code hard_spheres.m (making
sure to update the position and velocity at n+1).
___
Note. The logical operation ‘or’ is denoted as a single pipe (|) or double pipe (||) in
‣ Challenge coding. How would you code periodic boundary
MATLAB. The double pipe is the short-circuited version of the logical operation;
conditions instead of hard walls thein Matlab/Python?
second term is not evaluated if the rst term is true.
19 | Dabo 19
fi
fi
 Summary

‣ Today, we have introduced the hard sphere model.

‣ It is a simpli ed description for an ideal gas and was one of the

rst models ever simulated on a computer.

‣ We know how to write the kinetic energy K, potential energy U,

and the total energy E for the hard sphere model.

‣ The potential energy U of the hard spheres involves a pair

interaction function ϕ that ‘goes to zero at in nity and goes to
in nity at zero.’

‣ We know how to handle particle-wall collisions using appropriate

if statements in Matlab/Python.

20 | Dabo
fi
fi
fi
fi