Exercise

The first step was to write a function for the first derivative of the potential energy.

def first_derivative(r, A, B):
    """
    The first derivative of the Lennard-Jones potential model. 
    
    Args:
        r (float): Atom-atom distance (Å).
        A (float): Interaction parameter (eVÅ^12).
        B (float): Interaction parameter (eVÅ^6).
        
    Returns:
        (float): Potential energy.
    """
    return -12. * A / r ** 13 + 6 * B / r ** 7

We when implement a gradient descent method. alpha can have different values, but alpha=100 finds the minimum within 30 steps from r=6.0 initially.

A = 1e5
B = 40

alpha = 100
r = 6.0
r_list = [r]
for i in range(30):
    E_dash = first_derivative(r, A, B)
    r = r - alpha * E_dash
    r_list.append(r)

It is easier to show the convergence with a plot.

import matplotlib.pyplot as plt

plt.plot(range(31), r_list, 'o--')
plt.xlabel('Iteration')
plt.ylabel('r')
plt.show()