Minimisation

We have created a structure of atoms, and we're able to calculate the potential between those atoms, and the force on each atom. From there we are going to minimise this structure. When we minimise the structure, we will be moving the atoms so that the structure keeps its symmetry, but not necessarily its shape, such that the potential of that structure is at a minimum. At that configuration, the atoms should be completely stable, that is, the force on them should be zero and they should not move.

The method that we will use to minimise the structure is called the steepest descent method. The steepest descent method initially takes the negative gradient of the function (in our case, the force), and then determines a value $\alpha$, to solve the following one-dimensional minimisation problem:

$$\min_{\alpha} f(x_k - \alpha\nabla f(x_k))$$

Note that $f(x)$ in our case is the potential. Once we have determined $\alpha$, we then update the position of the atoms, using the following formula:

$x_{x + 1} = x_k - \alpha_k\nabla f(x_k)$

We then find the new negative gradient, and then repeat the process, until we have reached the point where the total force between all the atoms is sufficiently low to say that the structure is stable.

There are several methods of solving the one-dimensional minimisation problem. The one that we will use is called the golden section search. The golden section search takes an interval, $\lbrack a, b\rbrack$, on a function, and finds the minimum point on that interval. It does this by choosing two points, $x_1$ and $x_2$, and then calculates the functions at those points, $f(x_1)$ and $f(x_2)$. If $f(x_1) > f(x_2)$, it repeats the process on the interval $\lbrack f(x_1), b\rbrack$, and then chooses a new $x_1$ and $x_2$.

Now if we choose $x_1$ and $x_2$ to have the relative positions $\tau$ and $1 - \tau$, where $\tau^2 = 1 - \tau$, then we do not need to calculate the function at one of the new points, as one of the points will have evaluated to be the previous lowest point value.

One problem with the golden section search, is that if the function is not unimodal on the starting interval, it does not guarantee that the global minimum will be found. This means that it could find a minimum that is greater than the current potential. The nature of our potential function is that it is next to impossible for us to give it an interval that is unimodal, however, we can ensure that the result will be less or equal to the current potential. We do this by setting the starting interval such that the first $x_1$ chosen is zero. This is achieved by using the following formula:

$a = - x_2(\tau + 1)$ $b = x_2(\tau + 2)$

where $x_2$ is any number we choose. For our program, it is sufficient to choose $x_2$ to be one, as it is equivalent to moving the atoms in an amount equal to the force on the atoms.

The following is our implementation of the golden section search. It takes three arguments, the structure of atoms, the search direction, which is the gradient, and the tolerance, which is what the difference between $f_1$ and $f_2$ must be below in order to be considered minimised.

<Golden Section Search>=
def goldenSectionSearch(atoms, searchDirection, tol):
    t = (math.sqrt(5) - 1) / 2
    a = -t - 2
    b = t + 1
    x1 = a + (1 - t) * (b - a)
    f1 = evalPotential(nextAtoms(atoms, x1, searchDirection))
    x2 = a + (t * (b - a))
    f2 = evalPotential(nextAtoms(atoms, x2, searchDirection))
    while (abs(f1 - f2) > tol):
        if f1 > f2:
            a = x1
            x1 = x2
            f1 = f2
            x2 = a + (t * (b - a))
            f2 = evalPotential(nextAtoms(atoms, x2, searchDirection))
        else:
            b = x2
            x2 = x1
            f2 = f1
            x1 = a + (1 - t) * (b - a)
            f1 = evalPotential(nextAtoms(atoms, x1, searchDirection))
    return x1

Used below.

The nextAtoms routine evaluates the position of the new structure given a search direction and $\alpha$.

<Next Atoms>=
def nextAtoms(atoms, alpha, searchDirection):
    newAtoms = []
    for i in range(len(atoms)):
        newAtoms.append((atoms[i][0] + alpha * - searchDirection[i][0],
                         atoms[i][1] + alpha * - searchDirection[i][1],
                         atoms[i][2] + alpha * - searchDirection[i][2]))
    return newAtoms

Used below.

Now that we have the golden section search and next atoms routines, we can make our minimisation function. The convergence condition depends on the total force in the system compared to a function of the number of atoms - the more atoms the higher the total force will be. At each iteration we calculate a new tolerance for the golden section search, this is based on the proportional change of the last iteration in the potential energy.

<Minimise Initialisation and Loop>=
tol = 0.01
atoms = createSimpleCube(rside, nside)
force = evalForce(atoms)
totalPotential = evalPotential(atoms)
forcetol = 0.01 * len(atoms)
totalForce = forcetol + 1
<Create Balls>
while totalForce > forcetol:
    alpha = goldenSectionSearch(atoms, force, tol)
    atoms = nextAtoms(atoms, alpha, force)
    oldPotential = totalPotential
    force = evalForce(atoms)
    totalPotential = evalPotential(atoms)
    tol = abs((totalPotential - oldPotential) / totalPotential) * 10
    <Move Balls>
    totalForce = 0
    for f in force:
        totalForce = totalForce + f[0] ** 2 + f[1] ** 2 + f[2] ** 2

Used below.

*