Minimisation with scipy.optimize.minimize

Minimisation with scipy.optimize.minimize#

scipy.optimize.minimize() is a powerful tool for numerical optimisation that finds the minimum of any scalar function (a function that returns a single numerical value). Unlike the methods we have implemented so far, scipy.optimize.minimize() uses more sophisticated algorithms that automatically handle many of the challenges we have encountered, such as choosing appropriate step sizes and dealing with difficult potential energy surfaces.

The following code demonstrates how minimize() can be used to find the minimum of a harmonic potential.

from scipy.optimize import minimize

def harmonic_potential(r, k, r0):
    """Calculate the energy in eV for a harmonic potential U(r) = k/2 * (r - r0)^2.
    
    Args:
        r (float): Distance parameter to optimise.
        k (float): Spring constant in eV/Ų.
        r0 (float): Equilibrium distance in Å.
    
    Returns:
        float: Potential energy in eV.
    """
    return 0.5 * k * (r - r0)**2

# Initial guess and parameters
x0 = 1.0  # Starting guess for r
k = 36.0  # Spring constant
r0 = 0.74  # Equilibrium distance

result = minimize(harmonic_potential, x0=x0, args=(k, r0))

result

  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: 1.2580072611001788e-17
        x: [ 7.400e-01]
      nit: 1
      jac: [ 2.983e-07]
 hess_inv: [[1]]
     nfev: 6
     njev: 3

Let us break this down step by step:

First, we import the minimize() function from scipy.optimize.

from scipy.optimize import minimize

Next, we define a function for our harmonic potential:

def harmonic_potential(r, k, r0):
    """Calculate the energy in eV for a harmonic potential U(r) = k/2 * (r - r0)^2.
    
    Args:
        r (float): Distance parameter to optimise.
        k (float): Spring constant in eV/Ų.
        r0 (float): Equilibrium distance in Å.
    
    Returns:
        float: Potential energy in eV.
    """
    return 0.5 * k * (r - r0)**2

Note that we have ordered the three arguments to this function so that the parameter we want to minimise with respect to, r, is the first parameter. The second and third parameters, k and r0, will be fixed for a given potential energy surface.

Then, we define variables to store our initial guess for the equilibrium bond length and the parameters that define our harmonic potential.

# Initial guess and parameters
x0 = 1.0  # Starting guess for r
k = 36.0  # Spring constant
r0 = 0.74  # Equilibrium distance

Now we call the minimize() function:

result = minimize(harmonic_potential, x0=x0, args=(k, r0))

We pass in three arguments:

  • The first argument is the name of the function that we want to minimise. In this case our function is called harmonic_potential.

  • The second argument is our starting guess for the parameter we are minimising with respect to, so in this example, we pass our guess for our equilibrium bond length of 1 Å.

  • The third argument is a tuple of parameters that will also get passed to the function we are minimising, that are fixed during the minimisation. In this case, we need to define the shape of our harmonic potential energy surface by providing values for k and r0. Because our harmonic_potential() function takes two additional parameters, the tuple we pass to minimize() needs to contain two parameters. The order of these parameters also needs to match the order of the corresponding arguments in the function being minimised; so our tuple is structured as (k, r0).

minimize() returns an OptimizeResult object that contains lots of information about the optimisation, such as whether it was deemed to have successfully found a minimum, how many iterations it took, and the estimated minimum value.

We can access the estimated minimum value directly using

result.x

array([0.74])

Here, we are optimising a one-dimensional function (the energy varies with bond-length), and we get back a one-dimensional solution array. If we were minimising a function over more dimensions (e.g., the energy of a linear triatomic molecule, which then depends on two bond lengths) then the solution array would contain the minimisation solution for each dimension.

Exercises#

Exercise 1: Lennard-Jones Potential with scipy#

Use scipy.optimize.minimize() to perform a geometry optimisation for the Lennard-Jones potential:

\[U_\mathrm{LJ} = \frac{A}{r^{12}} - \frac{B}{r^6}\]

with \(A\) = 1 × 105 eV Å12 and \(B\) = 40 eV Å6.

Part A:

Write a function lennard_jones(r, A, B) to calculate this potential.

Part B:

Run minimize() with three different starting values: \(r\) = 3.2 Å, \(r\) = 4.4 Å, and \(r\) = 6.0 Å.

Part C:

For each starting point, print:

  • The starting position

  • The final optimised position

  • The number of iterations required

Part D:

Do all starting points converge to the same minimum? Why or why not?

Exercise 2. Multiple Minima in 1,2-Dichloroethane#

The potential energy associated with changes in the Cl-C-C-Cl dihedral angle, \(\theta\), in 1,2-dichloroethane can be modelled as

\[U(\theta) = \frac{1}{2}\left\{A_1\left[1+\cos\theta\right] + A_2\left[1+\cos2\theta\right]+ A_3\left[1+\cos3\theta\right]\right\}\]

with \(A_1\) = 55.229 kJ mol\(^{-1}\), \(A_2\) = 3.3472 kJ mol\(^{-1}\), and \(A_3\) = -58.576 kJ mol\(^{-1}\).

This exercise explores how the choice of starting geometry affects which minimum the optimiser finds.

Part A: Define and Visualise the Potential Energy Surface

  1. Write a function dihedral_potential(theta, A1, A2, A3) that calculates the potential energy for a given dihedral angle \(\theta\) (in radians).

  2. Create an array of \(\theta\) values from \(-\pi\) to \(\pi\) using np.linspace() with at least 200 points.

  3. Calculate the potential energy at each point and create a plot showing \(U(\theta)\) vs \(\theta\).

  4. By examining your plot, identify approximately how many minima exist and roughly where they are located.

Part B: Test Individual Starting Points

  1. Use scipy.optimize.minimize() to find the minimum, starting with \(\theta_0 = 0\) radians as your initial guess. What minimum does the optimiser find?

  2. Repeat the optimisation starting from \(\theta_0 = \pi\) radians. Does it find the same minimum?

  3. Try one more starting point of your choice. Based on these three tests, what pattern do you observe?

Part C: Systematic Exploration of Starting Points

Now we will systematically test many different starting points to understand the overall behaviour.

  1. Create an array of 20 equally spaced starting angles between \(-\pi\) and \(\pi\) using np.linspace(). Store this in a variable called start_angles.

  2. Create two empty lists called start_positions and final_positions.

  3. Write a for loop that:

    • Iterates through each angle in start_angles

    • Performs an optimisation starting from that angle

    • Appends the starting angle to start_positions

    • Appends the optimised angle (from result.x[0]) to final_positions

Part D: Visualise the Basins of Attraction

  1. Create a scatter plot with:

    • Starting angles on the x-axis (start_positions)

    • Final optimised angles on the y-axis (final_positions)

    • Appropriate axis labels and title

    This plot shows you which starting positions lead to which minima. Each horizontal band of points represents a different minimum (basin of attraction).

Part E: Analysis

  1. How many distinct minima does your scatter plot reveal? (Look for horizontal bands of points.)

Note

You may notice points clustered near both \(\theta \approx +3.14\) (near \(+\pi\)) and \(\theta \approx -3.14\) (near \(-\pi\)). Before counting these as separate minima, remember that dihedral angles are periodic: \(\theta = +\pi\) and \(\theta = -\pi\) represent the same molecular conformation (both are 180°).

Check your potential energy plot from Part A: you should see that \(U(+\pi) = U(-\pi)\). These two sets of points represent convergence to the same minimum approached from different directions, not two different minima.

This is a consequence of how the optimiser works: it finds a nearby minimum by following gradients, and doesn’t “wrap around” the periodic boundary. Starting from \(\theta = +2\) radians naturally leads to the minimum expressed as \(\theta \approx +\pi\), whilst starting from \(\theta = -2\) radians leads to the same minimum expressed as \(\theta \approx -\pi\).

  1. For each minimum, estimate the range of starting angles that lead to that minimum. These ranges are called the “basin of attraction” for each minimum.

  2. Based on your results, if you wanted to find the global minimum (the lowest energy structure), what strategy would you need to use? (Hint: look back at your potential energy plot from Part A to see which minimum has the lowest energy.)

Contents