Crystal Formation Using Lennard-Jones Potential

Introduction

Before we do the math for crystal formation with many atoms, we should first compute the range of oscillation of two atoms in their mutual attractive/repulsive potential energy wells. After all that will be approximately the final mode of oscillation of the atoms in the crystal. We can actually see the two atom oscillation in the animation prior to the crystal formation. After the crystal is formed, but before all kinetic energy is extracted, we can see the oscillation of the atoms in the multi-atom potential.

Calculation of turnaround points for 2 atoms

$V = 4 ε ({(\frac{σ}{r})}^{12} - {(\frac{σ}{r})}^{6})$

(1.1)

Let

$u = {(\frac{σ}{r})}^{6}$

(1.2)

The total energy is E+V where E is the kinetic energy, $\frac{m v^{2}}{2}$

The excess energy with respect to the bottom of the V potential well is

$(V - V_{\min}) - E = 0 = 4 ε (u^{2} - u) + ε - E$

(1.3)

Setting up to solve for u we get:

$4 ε (u^{2} - u) + (ε - E) = 0$

(1.4)

$u = \frac{4 ε \pm \sqrt{16 ε^{2} - 16 (ε - E)} ε}{8 ε} = 1 \pm \sqrt{\frac{E}{ε}}$

(1.5)

Obviously

$r = u^{- \frac{1}{6}} σ = \frac{σ}{{(1 \pm \sqrt{\frac{E}{ε}})}^{\frac{1}{6}}}$

(1.6)

If we set $E = \frac{m v^{2}}{2} = ε$ and set the initial distance $r = 2^{\frac{1}{6}} σ$ where V is minimum and equal to -ε, then m will have just enough energy to escape the potential well. Since we don't want the atoms to go to infinity we must set the kinetic energy slightly less than ε.

Boundary Potential

It was found that using hard boundaries with soft atom-atom collisions was not conducive to maintaining finite kinetic energy. Therefore soft forces for the boundaries were implemented. If the bounds of the box are at x=0,w and y=0,h then a good choice of soft boundary force distribution is

$\begin{array}{l} F_{x} = ε_{b} ({(\frac{σ_{b}}{x})}^{12} - {(\frac{σ_{b}}{w - x})}^{12}) \\ F_{y} = ε_{b} ({(\frac{σ_{b}}{y})}^{12} - {(\frac{σ_{b}}{w - y})}^{12}) \end{array}$

(1.7)

where ε_b and σ_b are analogous to like-named quantities for the atom-atom potential and F_x and F_y are the force components due to this boundary potential.