Metropolis Monte Carlo (MMC) Algorithm
The MMC algorithm implements sequential steps that allow the chosen
atomic sphere to move or not to a new spatial position consistent with
the atomic bonding structure. The list of the (x, y, z) positions for
the centers of each individual sphere were updated after each MMC step
as these are the values to be used for the desired calculations of
average thickness of a monolayer and the formation of gauche-bonds.
Using classical physics, the algorithm computes the total energy of a
sphere at a particular “initial” position and also at a randomly
chosen new position, a “final” position. The algorithm specifies the
calculations to be done to either accept or reject the possible
new-final position.
The total energy of a sphere is given by: i) the bond stretching energy,\(E_{b}\), between neighboring bonded spheres (on the same molecule)
and, ii) the dispersion energy, \(E_{d}\), between non-bonded spheres on
the same molecule, or spheres on different molecules . Thus, the total
energy \(E,\ \) for a particular sphere is
\(E=\ E_{b}+\ E_{d}\) [1]
Where \(E_{b}\), was represented by
\(E_{b}={3.8\ \ \left(r-3\right)}^{4}\text{\ \ }\) [2]
Here the unit is one-half Å, and r represents the center-to-center
distance of bonded (adjacent) spheres. The energy of 3.8 eV is the
energy of a single C-C bond, . The value of 3.8
(eV/units4) was used to cause the bonding energy to
reach 3.8 eV at a distance of ± 1 unit from the minimum.
\(E_{d}\) was calculated from equation 3
\(E_{d}=\left(\frac{j}{r^{8}}-\frac{k}{r^{6}}+l\right)\frac{5}{l}\ \)[3]
Where, as above, r represents the distance between two non-bonded
spheres. r6 is the well-known attractive London
dispersion energy . The short-range repulsion ofr-8 was used rather than the typicalr-12 (. We chose r-8 in
order to permit the formation of gauche bonds, by making the electron
cloud less repulsive than when using r-12 ─ if
no significant number of gauche bonds arise usingr-8 , then they will not arise usingr-12 . \(j,\ k\ \text{and}\text{\ l}\) are
coefficients to assign an energy of 0 as the minimum, based on the
equilibrium distance of 4.9 Å between parallel hydrocarbon chains .
Notice that the total energy does not include electrostatic energy, as
this was ignored so that our simulation excludes BL.
In this system, the solvent oil is assumed to be a neutral background
which does not play a significant role in the formation of gauche bonds,
hence it was not necessary to consider it.
An MMC algorithm starts by selecting a sphere to move. This simulation
followed a sequential order starting from the “center” sphere in each
molecule. This center of the molecule was held constant throughout the
simulation. Molecules were prevented from moving perpendicular to the
long axis (Fig 1).
PLACE FIGURE 1 HERE
Since each molecule’s center was held fixed throughout the simulation,
its location and the distances between atoms, and therefore, molecules
were known.
An MMC step was started by recording the spatial position of the
selected sphere (initial position) and by computing its initial
energy\(\ E_{i}\) (eq. 1). Note that the \(E_{b}\ \)(eq. 2) component
was computed using the spheres bonded to the selected sphere in this
initial position, while \(E_{d}\) (eq. 3) was computed using the
“non-bonding” neighboring spheres. A random number was then used to
select a possible new position (final position) for the selected sphere.
The final energy \(E_{f}\) was computed for this group at the new
coordinates. The newly computed \(E_{f}\) was subtracted from \(E_{i}\)(eq. 4)
\(\Delta E=E_{f}-E_{i}\) [4]
The final calculation in the MMC step was to determine if the new
position was to be accepted by looking at the energy change \(\Delta E\)(eq. 4) value. If \(\Delta E\) was negative, a lower energy state was
achieved, hence the move to the new position was accepted. If the change
in energy was positive, a new calculation was made to decide if this
higher energy state was acceptable. Equation 5, which takes into account
the entropy of the system, details the calculation made
\(g=\ e^{(\frac{-\Delta E}{k_{B}T})}\ \) [5]
where kB is Boltzmann’s constant (eV) and T = 300
°K for this simulation.
A random number R was then chosen in the range \(0\leq R<1\).
If \(R\leq g\), then the new position was accepted if the move was
allowed. However, if R > g the new position
was not accepted.
After the last calculation, a new sphere was selected and all the
previous steps repeated. Spheres were selected sequentially starting at
the center. Once all spheres in the molecule were treated with one MMC
step, the software routine comes back to the first central-fixed group
and started over. Each group was attempted to be moved sufficient times
for a total of 100,000 times (105 MMC steps). Ten
replicates were run and averaged to obtain the results.