Minimisation Libraries

Part of the power of the Python programming language comes from the vast ecosystem of pre-existing libraries that are available to use. Previously, you have written your own minimisation algorithms in Python, using loops and functions. Numerical function minimisation is a very common scientific problem, and, perhaps unsurprisingly, there is a library that we can use to do this. The library is called scipy (pronounced Sigh-Pie), and includes an enormous number of common scientific computing utilities. In this section, we will focus on the minimize method (note the American spelling).

The minimize method is located within the optimize module of the scipy library, so to access this, we must import it.

from scipy.optimize import minimize

Let us first investigate the documentation of the minimize function, this is available online or it can be printed in the Jupyter Notebook using

print(minimize.__doc__)

or

minimize?

Because this documentation is long we will not reproduce it fully here, but it is useful to know where to look when figuring out how a new function or method behaves.

From the documentation, we can see that the minimize function takes two required arguments:

• fun: this is the function that we will minimize.

• x0: this is our initial guess of the parameter used in the minimisation.

The only optional argument we will look at from the documentation is args, which are the invariant (unchanging) parameters in our model that should be passed to our function. For example, in the case of our Lennard-Jones potential, the precise shape of the function depends on the invariant parameters \(A\) and \(B\), so if we want to ue minimize to find the minimum of our function, we need to specify the relevant values of \(A\) and \(B\).

To minimize the potential energy from the Lennard-Jones potential, we write a function for this, this is the fun discussed above. The parameter we are minimising with respect to (in this case, the interatomic separation, \(r\)) needs to be the first argument in our function, followed by any additional parameters that we are not varying as part of our minimisation (here \(A\) and \(B\)).

def lennard_jones(r, A, B):
    """
    Lennard-Jones potential energy.
    
    Args:
        r (float): the distances to find the energy for (Å).
        A (float): the repulsive prefactor (eV/Å^12).
        B (float): the attractive prefactor (eV/Å^6).
    
    Returns:
        (float): the potential energy at each r.
    """
    return A / r**12 - B / r**6

Note in the documentation, it says that the minimize function should be passed a function with the form

fun(x, *args) -> float

Above, this is mapped out, were r is our x and A, B are the *args (note the * means this can be any number of args). The args are not varied during the minimisation process. The minimization function works as follows.

A = 1e5
B = 40

result = minimize(lennard_jones, [3.5], args=(A, B))

print(result)

  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: -0.003999999054071188
        x: [ 4.135e+00]
      nit: 8
      jac: [-5.648e-06]
 hess_inv: [[ 5.471e+01]]
     nfev: 18
     njev: 9

The object that is returned by the minimize function is a special object which gives a lot of information about the minimisation. The value we are interested in is the minimum energy difference, which is the x attribute of this object.

print(result.x)

[4.13485048]

This agrees well with the value that we found when writing our own minimisation methods.

Exercise

Write a function to model the following potential energy function,

\[ E = \frac{1}{2} \{A_1[1 + \cos{(\phi)}] + A_2[1 + \cos{(2\phi)}] + A_3[1 + \cos{(3\phi)}]\}, \]

this function describes the potential energy of the Cl-C-C-Cl dihedral angle in a molecule of 1, 2-dichloroethane, where \(A_1 = 55.229\) kJ/mol, \(A_2 = 3.3472\) kJ/mol, \(A_3 = -58.576\) kJ/mol and \(\phi\) is the dihedral angle.

1. Plot this function between \(-\pi\) and \(\pi\), remember to label your axes.

2. Using minimize test different initial values of x0 for \(\phi\). Consider why the result is so dependent on the starting position for this function (it may help to plot your result on the function).

Tip

It is advised to use the NumPy trigonmetric functions in this exercise.