Minimisation with scipy.optimize.minimize

Contents

Minimisation with scipy.optimize.minimize#

scipy.optimize.minimize() is a powerful tool for numerical optimization that finds the minimum of any scalar function (a function that returns a single numerical value).

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
    
    Parameters:
        r (float): Distance parameter to optimize
        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
    
    Parameters:
        r (float): Distance parameter to optimize
        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 for 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 optimisation 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#

1.#

Use scipy.optimize.minimize() to perform a geometry optimisation for your Lennard-Jones potential, again with starting values of \(r\) = 3.2 Å, \(r\) = 4.4 Å, and \(r\) = 6.0 Å.

2.#

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

\[U = \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.3.472 kJ mol−1, and \(A_3\) = -58.576 kJ mol−1.

  • Plot this function for \(-\pi \leq \theta \leq \pi\).

  • Using scipy.optimize.minimize() test the optimisation behaviour for different initial guesses of \(\theta\). In what ways does the minimisation result depend on the choice of initial guess?

Contents