Page 379 - Numerical Methods for Chemical Engineering

368 7 Probability theory and stochastic simulation

the position of lowest F(x) found by any particle to date. The update in the velocity and
position of each particle is then
[α] [α] [α] [α] [α]
v ← v + c 1 u k p − x + c 2 u g k − x k ∈ [1, dim(x)]
k k k k k k
x [α] ← x [α] + v [α] (7.280)
u, u are vectors of independent random numbers, uniformly distributed on [0, 1] and c 1 ,

c 2 are tunable learning parameters. A maximum allowable velocity may be specifed as
well. Write a program that uses PSO with approximately 10–50 particles, and compare its
performance to simulated annealing and genetic algorithms for the cost function (7.265).
7.C.1. Atoms of an ideal gas such as argon interact only through dispersion forces that are
modeled by the Lennard–Jones pairwise interaction:

12 6
σ σ
U LJ (R αβ ) = 4ε − r αβ = R α − R β (7.281)

R αβ R αβ
This interaction includes a strong short-range repulsion when R αβ « σ and a weaker long-
range attraction when R αβ » σ.As U LJ (R αβ ) → 0as R αβ →∞, it is common to set the
interaction energy to zero with R αβ > r cut , where a common choice of cutoff radius is
r cut = 3.5σ.
For a system of N Lennard–Jones atoms, the total potential energy of the system is the
sum of the pairwise interactions from each unique pair of atoms:
N N

U tot (q) = U LJ (R αβ ) (7.282)
α=1 β=α+1
q is the state vector of all atomic positions. Using a cutoff radius, when computing U tot we
need consider only those pairs of atoms with R αβ ≤ r cut .
If the molecular weight of each atom is m a , we simulate the system at density ρ by
using periodic boundary conditions on a cubic box of dimension L × L × L, such that
3
ρL = m a N. We assume that this simulation cell is surrounded by identical copies, with
images of atom α at
R [m 1 ,m 2 ,m 3 ] = R α + m 1 Le [1] + m 2 Le [2] + m 3 Le [3] m j = 0, ±1, ±2,... (7.283)
α
If L > r cut , then for each pair (α, β) of atoms in (7.282), we need only consider the images
of each atom that are closest to each other.
3
At constant N, constant volume V = L , and constant temperature T, we sample the
system at equilibrium using the Metropolis Monte Carlo method for the Boltzmann distri-
bution P(q) ∝ exp[−U tot (q)/k b T ]. The state is updated by selecting at random one or more
atoms, and then randomly translating them. If only one atom is selected, then only a small
fraction of the pairwise interactions in (7.282) change value and need to be recomputed.
In practice, the CPU time necessary to compute the pairwise interactions is reduced by a
strategy such as using a cell list, in which the simulation box is divided into nonoverlapping
subcells of length just greater than r cut . Then, each atom is assigned to one of these subcells,
and when computing the total energy, only atoms in the same or neighboring subcells need
to be considered.