CL 654: Multiscale Modelling and Simulation

SOLUTION TO ASSIGNMENT
TOTAL: 40 Marks

There are SEVEN questions on this assignment. Answer all questions, clearly showing the
steps.

1. [5 marks] For a collection of N ideal monoatomic gas atoms (where N is very large) in volume V , the
canonical ensemble partition function is given by
 3N/2
2πm
h2 β
Q(N, V, β) = VN
N!

3N
where m is the mass of an atom and h is the Planck’s constant. Show that U = where β = 1/(kB T )
  2β
∂ ln Q
(kB : Boltzmann constant). You may use the result U = kB T 2 .
∂T V,N

Solution: We have
 
2 ∂ ln Q
U = kB T
∂T V,N

Also,

 3N/2 
2πm
 3N/2  
 h2 β 
N 2πm N 3N 2πm
ln Q = ln  V  = ln + ln V − ln N ! = ln + N ln V − ln N !

 N!  h2 β 2 h2 β

Now,
            
∂ ln Q ∂ 3N 2πm 3N ∂ 2πm 3N ∂ ln β
= ln = ln − ln β = −
∂T V,N ∂T 2 h2 β V,N 2 ∂T h2 V,N 2 ∂T V,N

Using β = 1/(kB T ), we get

  
1

∂ ln Q

3N  ∂ ln kB T 3N

∂ ln(kB T )

3N
=−  = =
∂T V,N 2 ∂T 2 ∂T V,N 2T
V,N

Therefore,  
∂ ln Q 3N 3N kB T 3N
U = kB T 2 = kB T 2 × = =
∂T V,N 2T 2 2β

1
2. [5 marks] The entropy in terms of the canonical ensemble partition function is given by
 
∂ ln Q
S = kB ln Q + kB T
∂T V,N

Using this result, for the N -particle ideal monoatomic gas in problem 1 occupying volume V , derive the
Sackur-Tetrode equation: "  3/2 #
S V 4πmU 5
= ln +
kB N N 3h2 N 2

You may use Stirling’s approximation (valid for large N ): ln N ! = N ln N − N and the result obtained in
problem 1 (i.e., U = 3N/2β).

Solution: We have, from the result in problem 1,

   
2 ∂ ln Q 3N ∂ ln Q 3N 3
U = kB T = ⇒ kB T = = N kB
∂T V,N 2β ∂T V,N 2T β 2

Also from problem 1,  

3N 2πm
ln Q = ln + N ln V − ln N !
2 h2 β

Now,
   
∂ ln Q 3N kB 2πm 3
S = kB ln Q + kB T = ln + N kB ln V − kB ln N ! + N kB
∂T V,N 2 h2 β 2

Using Stirling’s approximation, ln N ! = N ln N − N , we get

 
3N kB 2πm 3
S= ln 2
+ N kB ln V − N kB ln N + N kB + N kB
2 h β 2
  " 3/2 #  
S 3 2πm 3 2πm V 5
⇒ = ln + ln V − ln N + 1 + = ln + ln +
kB N 2 h2 β 2 h2 β N 2
"  3/2 #
S V 2πm 5
⇒ = ln +
kB N N h2 β 2

Finally, as U = 3N/(2β), we can replace β with 3N/(2U ) to get

"  3/2 #
S V 4πmU 5
= ln 2
+
kB N N 3h N 2

2
3. [5 marks] Consider an argon atom interacting with an oxygen molecule (O2 ). The non-bonded van der
Waals interaction between argon and oxygen can be described by the Lennard-Jones (LJ) potential:
" 12  6 #
LJ σij σij
Uij (r) = 4ϵij −
rij rij

where i and j represent the species (argon and oxygen in this case), rij is the interatomic distance, and
σij and ϵij are the LJ parameters. The oxygen atoms are at positions (3.2Å, 0, 0) and (3.2Å, 0, 1.2Å), and
the argon atom is at the origin (0, 0, 0). If σO = 3.0 Å, ϵO /kB = 52 K, σAr = 3.4 Å, ϵAr /kB = 120 K,
determine the total non-bonded van der Waals interaction energy (in J/mol) between the argon atom and
the oxygen molecule (O2 ). The Ar − O LJ parameters may be calculated using the Lorentz-Berthelot
mixing rules:
σii + σjj √
σij = and ϵij = ϵii ϵjj
2
Note: The ϵO /kB values in K can be converted to ϵ values in J/mol by multiplying with R = 8.314
J/(mol·K). For example, ϵAr /kB = 120 K ⇒ ϵAr = 120 × 8.314 = 997.68 J/mol.

Solution: We have two pairwise non-bonded interactions, between an argon atom (numbered as 1) and
two oxygen atoms (numbered 2 and 3). The distance between the argon atom and the first oxygen atom
is given by
p p
r12 = (x2 − x1 )2 + (y2 − y1 )2 + (z2 − z1 )2 = (3.2 − 0)2 + (0 − 0)2 + (0 − 0)2 = 3.2 Å

The distance between the argon atom and the second oxygen atom is given by
p p
r13 = (x3 − x1 )2 + (y3 − y1 )2 + (z3 − z1 )2 = (3.2 − 0)2 + (0 − 0)2 + (1.2 − 0)2 = 3.4176 Å

Given: σO = 3.0 Å, σAr = 3.4 Å, ϵO /kB = 52 K, ϵAr /kB = 120. So, ϵO = 52 × 8.314 J/mol = 432.328
J/mol and ϵAr = 120 × 8.314 J/mol = 997.68 J/mol. The Ar − O cross interaction parameters are given
by
σO + σAr 3 + 3.4 √ √
σAr−O = = = 3.2 Å and ϵAr−O = ϵAr ϵO = 997.68 × 432.328 = 656.753 J/mol
2 2

The total non-bonded van der Waals interaction energy (in J/mol) between the argon atom and the oxygen
molecule is given by the sum of the pairwise interactions between the argon atom and the two oxygen
atoms constituting the O2 molecule. So,
" 12  6 # " 12  6 #
LJ LJ σ12 σ12 σ13 σ13
UAr−O2 = U12 + U13 = 4ϵ12 − + 4ϵ13 −
r12 r12 r13 r13

We have σ12 = σ13 = σAr−O = 3.2 Å and ϵ12 = ϵ13 = ϵAr−O = 656.753 J/mol. Substituting the values,
we get

" 12  6 # " 12  6 #

LJ LJ 3.2 3.2 3.2 3.2
U12 + U13 = 4 × 656.753 × − + 4 × 656.753 × −
3.2 3.2 3.4176 3.4176

LJ LJ
⇒ U12 + U13 = 2627.012 × (0 − 0.21977) ≈ −577.34 J/mol

3
4. [5 marks] Consider two atoms i and j directly bonded to each other. The bond stretching energy is given
by the following harmonic potential: Ub = Kb (rij − r0 )2 where Kb , r0 are the potential parameters, and
rij is the centre-to-centre distance between atoms i and j. Derive an expression for the force acting on
atom i. If Kb = 270 kcal/(mol.Å2 ) and r0 = 1.5 Å, what is the value (magnitude) of the force acting on
atom i?
Note: The force on atom i in the x direction is given by
∂Ub ∂Ub ∂rij
fxi = (−∇Ui )x = − =−
∂xi ∂rij ∂xi
p
where rij = |⃗rij | = |(⃗rj − ⃗ri )| = (xj − xi )2 + (yj − yi )2 + (zj − zi )2 .
Similarly, the force components on atom i in y and z directions are fyi = −∂Ub /∂yi and fzi = −∂Ub /∂zi
respectively. The total force is given by f⃗i = fxi î + fyi ĵ + fzi k̂

Note: The force on atom i in the x direction is given by

∂Ub ∂Ub ∂rij
fxi = (−∇Ui )x = − =−
∂xi ∂rij ∂xi
p
where rij = |⃗rij | = |(⃗rj − ⃗ri )| = (xj − xi )2 + (yj − yi )2 + (zj − zi )2 .
Similarly, the force components on atom i in y and z directions are fyi = −∂Ub /∂yi and fzi = −∂Ub /∂zi
respectively. The total force is given by f⃗i = fxi î + fyi ĵ + fzi k̂

Solution: The force on atom i in the x direction is given by

!−1
2 2 2
∂Ub ∂Ub ∂rij ∂Ub ∂rij ∂rij 1 ∂Ub ∂rij
fxi = (−∇Ui )x = − =− =− =−
∂xi ∂rij ∂xi ∂rij ∂xi ∂rij 2rij ∂rij ∂xi

We further have
2
∂rij ∂
= {(xj − xi )2 + (yj − yi )2 + (zj − zi )2 } = −2(xj − xi )
∂xi ∂xi
∂Ub
and = 2Kb (rij − r0 ). So,
∂rij
2
1 ∂Ub ∂rij 2Kb
fxi = − = (rij − r0 )(xj − xi )
2rij ∂rij ∂xi rij

Similarly,
2Kb 2Kb
fyi = (rij − r0 )(yj − yi ) and fzi = (rij − r0 )(zj − zi )
rij rij

The total force of atom i due to its bonded interaction with atom j is given by
2Kb h i
f⃗ i = fxi î + fxj ĵ + fxk k̂ = (rij − r0 ) (xj − xi ) î + (yj − yi ) ĵ + (zj − zi ) k̂
rij

2Kb ⃗rij
⇒ f⃗ i = (rij − r0 ) ⃗rij = 2Kb (rij − r0 ) = 2Kb (rij − r0 ) r̂ij
rij rij

where ⃗rij = ⃗rj − ⃗ri is the bond vector from atom i to atom j and r̂ij is the unit vector along ⃗rij . The
magnitude of the force acting on atom i is

f i = |f⃗ i | = 2Kb (rij − r0 ) |r̂ij | = 2Kb (rij − r0 ) = 2 × 270 (rij − 1.5) = 540 (rij − 1.5) kcal/(mol · Å)

where rij is in Å.

4
5. [5 marks] Consider the cubic simulation box shown
on the right. Periodic boundary conditions apply in
all three directions. The Lennard-Jones potential pa-
rameters of the particles in the system are σ = 3.4
Å and ϵ = 500 J/mol. The length of the simula-
tion box is 25 Å. The positions of particles i and j
in the simulation box are (−12, 5, 4) and (10, 5, 6)
respectively, where the coordinates are in Å. If the
Lennard-Jones interactions are truncated at a cutoff
distance of 10 Å and the minimum image convention
is used, determine the Lennard-Jones energy of in-
teraction between i and j (or its appropriate image).

Solution: Let us first calculate the distance between particle i and the minimum image of particle j. In
the x direction, we have xij = xj − xi = 10 − (−12) = 22. As xij is greater than the Lennard-Jones cutoff
distance, no LJ interaction between i and j in the central simulation box will be considered. Further, as
xij > (half the box length), the distance between particle i and the minimum image of particle j in the x
direction is xij,min = xij − (box length) = 22 − 25 = −3Å.
In the y direction, we have yij = yj − yi = 5 − 5 = 0. As yij < (half the box length), yij,min = yij = 0 Å.
In the z direction, we have zij = zj − zi = 6 − 4 = 2. As zij < (half the box length), zij,min = zij = 2 Å.
The overall distance between atom i and the minimum image of atom j is
q p
rij,min = x2ij,min + yij,min
2 2
+ zij,min = (−3)2 + (0)2 + (2)2 = 3.60555 Å

The corresponding LJ energy is

" 12  6 # " 12  6 #
LJ σij σij 3.4 3.4
Uij = 4ϵij − = 4 × 500 × − = −417.465 J/mol
rij rij 3.60555 3.60555

5
6. [5 marks] In a grand canonical Monte Carlo (GCMC) simulation, a particle is inserted at a randomly
chosen position in a cubic simulation box already containing 100 particles. The length of the box is 70 Å.
The non-bonded particle interactions are described by the Lennard-Jones (LJ) potential with σ = 3.4 Å
and ϵ = 0.24 kcal/mol. There are four other particles in the system within the LJ cutoff distance from the
inserted particle; these particles are at a distance of 3.3 Å, 3.5Å, 3.8 Å and 4 Å from the inserted particle.
The chemical potential is µ = −4 kcal/mol. The molar mass of the particles is 40 g/mol. The values
of the Avogadro number, Boltzmann’s constant and Planck’s constant are NAv = 6.022 × 1023 mol−1 ,
kB = 1.38 × 10−23 J/K and h = 6.626 × 10−34 J· s respectively. Determine the acceptance probability of
the move at T = 200 K.
√
Note: 1 cal = 4.184 J. The thermal de Broglie wavelength is defined as Λ = h/ 2πmkB T where m is
the mass of 1 molecule/particle. Convert molar mass to m by dividing with Avogadro number. As the
energies are given in per mole units, use β = 1/(RT ).

Solution: The mass of one particle is

M 40
m= = = 6.64 × 10−23 g = 6.64 × 10−26 kg
NAv 6.022 × 1023

The thermal de Broglie wavelength at T = 200 K is

h 6.626 × 10−34
Λ= √ = = 1.953 × 10−11 m = 0.1953 Å
2πmkB T 2π × 6.64 × 10−26 × 1.38 × 10−23 × 200

The change is potential energy, ∆U , is due to the interactions between the inserted particle and the other
particles in the system within the LJ cutoff distance of the inserted particle. Let the inserted atom be i
and the four particles within the cutoff distance be numbered 1, 2, 3 and 4. So,
" 12  6 # " 12  6 #
σi1 σi1 3.4 3.4
Ui1 = 4ϵi1 − = 4 × 0.24 × − = 0.2253 kcal/mol
ri1 ri1 3.3 3.3
" 12  6 # " 12  6 #
σi1 σi1 3.4 3.4
Ui2 = 4ϵi1 − = 4 × 0.24 × − = −0.1288 kcal/mol
ri1 ri1 3.5 3.5
" 12  6 # " 12  6 #
σi3 σi3 3.4 3.4
Ui3 = 4ϵi3 − = 4 × 0.24 × − = −0.2398 kcal/mol
ri3 ri3 3.8 3.8
" 12  6 # " 12  6 #
σi4 σi4 3.4 3.4
Ui4 = 4ϵi4 − = 4 × 0.24 × − = −0.2255 kcal/mol
ri4 ri4 4.0 4.0

So, ∆U = Ui1 + Ui2 + Ui3 + Ui4 = 0.2253 − 0.1288 − 0.2398 − 0.2255 = −0.3688 kcal/mol.
Given: µ = −4 kcal/mol, N = 100 and V = 703 = 343×103 Å. The acceptance probability of the insertion
move at T = 200 K is given by

343 × 103 (−4 + 0.3688) × 4.184 × 103

 
V
exp [β(µ − ∆U )] = exp
Λ3 (N + 1) 0.19533 × (100 + 1) kB T

 
15192.941
= 455896.2463 × exp − = 455896.2463 × exp(−9.13696) = 49.061
(8.314 × 200)

As this value is greater than 1, the acceptance probability is 1.

Note that while evaluating the term exp [β(µ − ∆U )], we have changed the units from kcal/mol to J/mol
by multiplying with a factor of 4.184 × 103 .

6
7. [10 marks (= 3 + 5 + 2)] Consider four consecutively bonded atoms, i, j, k and l, in a linear molecule
where i is bonded to j, j is bonded to k and k is bonded to l. The position of i, j, k and l (in Å) are (0,
0, 0), (1.5, 0, 0), (2.6, 0.9 , 0) and (3.6, 1.9, 0.6) respectively. The angle bending and dihedral potentials
are given by:
Uangle = Kθ (θijk − θ0 )2
1 1 1
Udihed = K1 [1 + cos(ϕijkl )] + K2 [1 − cos(2ϕijkl )] + K3 [1 + cos(3ϕijkl )]
2 2 2
where K = 150 kcal/(mol.rad2 ), θ0 = 120◦ = 2π/3 rad, K1 = 1.740 kcal/mol, K2 = −0.157 kcal/mol,
K3 = 0.279 kcal/mol. Calculate the following:

(a) [3 marks] The bond angle formed by atoms i, j and k (i.e., θijk ).
(b) [5 marks] The dihedral angle formed by the atoms i, j, k and l (i.e., ϕijkl ). The dihedral angle can
be obtained as the angle between ⃗rij × ⃗rjk and ⃗rjk × ⃗rkl .
(c) [2 marks] The dihedral potential energy for atoms i, j, k and l.

Note: The bond angle formed by atoms i, j and k, where j is the central atom, is given by the angle
between vectors r⃗ji and ⃗rjk , where ⃗rji = ⃗ri − ⃗rj , ⃗rjk = ⃗rk − ⃗rj , with ⃗ri , ⃗rj , ⃗rk being the position vectors
of atoms i and j respectively. If the coordinates of atoms i are xi , yi , zi , its position vector is defined as
⃗ri = xi î + yi ĵ + zi k̂, where î, ĵ, k̂ are the unit vectors along the x, y and z axes respectively.

Solution:
(a) The vectors ⃗rji and ⃗rjk are given by

⃗rji = ⃗ri − ⃗rj = (xi − xj ) î + (yi − yj ) ĵ + (zi − zj ) k̂ = (0 − 1.5) î + (0 − 0) ĵ + (0 − 0) k̂ = −1.5 î

⃗rjk = ⃗rk − ⃗rj = (xk − xj ) î + (yk − yj ) ĵ + (zk − zj ) k̂ = (2.6 − 1.5) î + (0.9 − 0) ĵ + (0 − 0) k̂ = 1.1 î + 0.9 ĵ

The magnitudes of vectors r⃗ji and ⃗rjk are

p p
|⃗rji | = (−1.5)2 = 1.5 Å and |⃗rjk | = (1.1)2 + (0.9)2 = 1.42 Å

The bond angle formed by the atoms i, j and k is given by

  !  
−1 ⃗rji · ⃗rjk −1 (−1.5 î) · (1.1 î + 0.9 ĵ) −1 −1.5 × 1.1
θijk = cos = cos = cos = cos−1 (−0.77465) = 140.77◦
|⃗rji ||⃗rjk | 1.5 × 1.42 2.13

(b) From part (a), the vectors ⃗rij and ⃗rjk are

⃗rij = −⃗rji = 1.5 î and ⃗rjk = 1.1 î + 0.9 ĵ

The vector ⃗rkl is given by

⃗rkl = ⃗rl − ⃗rk = (xl − xk ) î + (yl − yk ) ĵ + (zl − zk ) k̂ = (3.6 − 2.6) î + (1.9 − 0.9) ĵ + (0.6 − 0) k̂ = î + ĵ + 0.6 k̂

The cross product ⃗rij × ⃗rjk is given by

î ĵ k̂
⃗rij × ⃗rjk = 1.5 0 0 = (1.5 × 0.9 − 0 × 1.1) k̂ = 1.35 k̂
1.1 0.9 0

The cross product ⃗rjk × ⃗rkl is given by

î ĵ k̂
⃗rjk × ⃗rkl = 1.1 0.9 0 = (0.9 × 0.6) î − (1.1 × 0.6) ĵ + (1.1 × 1 − 0.9 × 1) k̂ = 0.54 î − 0.66 ĵ + 0.2 k̂
1.0 1.0 0.6

The magnitudes of the cross products, i.e., |⃗rij × ⃗rjk | and |⃗rjk × ⃗rkl |, are

7
 p 2 p 2
|⃗rij × ⃗rjk | = (1.35)2 = 1.35 Å and |⃗rjj × ⃗rkl | = (0.54)2 + (−0.66)2 + (0.2)2 = 0.8759 Å

So, the dihedral angle (ϕijkl ) is

   
(⃗rij × ⃗rjk ) · (⃗rjk × ⃗rkl ) 1.35 × 0.2
ϕijkl = cos−1 = cos−1 = cos−1 (0.2283) = 76.803◦
|⃗rij × ⃗rjk ||⃗rjk × ⃗rkl | 1.35 × 0.8759

(c) The dihedral potential energy is given by

1 1 1
Udihed = K1 [1 + cos(ϕijkl )] + K2 [1 − cos(2ϕijkl )] + K3 [1 + cos(3ϕijkl )]
2 2 2
1.74 0.157 0.279
⇒ Udihed = [1 + cos(76.803◦ )] − [1 − cos(2 × 76.803◦ )] + [1 + cos(3 × 76.803◦ )]
2 2 2

⇒ Udihed = 1.0686 − 0.1488 + 0.0506 = 0.9704 kcal/mol