import numpy as np
import matplotlib.pyplot as plt
import sympy as sy
from openmm import unit

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/tmp/ipykernel_3426/1768193541.py in <module>
      1 import numpy as np
----> 2 import matplotlib.pyplot as plt
      3 import sympy as sy
      4 from openmm import unit

ModuleNotFoundError: No module named 'matplotlib'

The Lennard-Jones potential

The Lennard-Jones 12-6 potential between two particles is defined by the following equation, where \(x\) is the distance between the particles, and \(\sigma\) and \(\epsilon\) are two parameters of the potential:

(1)\[ V(x) = 4 \epsilon \left[ \left( \frac{\sigma}{x} \right)^{12} - \left( \frac{\sigma}{x} \right)^6 \right] \]

Lets see the shape of this function:

def LJ (x, sigma, epsilon):
    
    t = sigma/x
    t6 = t**6
    t12 = t6**2
    
    return 4.0*epsilon*(t12-t6)

sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

xlim_figure = [0.01, 6.0]
ylim_figure = [-2.0, 10.0]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstrom
plt.plot(x, LJ(x, sigma, epsilon))
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
../../../_images/LJ_Potential_5_0.png

The way the LJ potential is built, the \(\sigma\) and \(\epsilon\) parameters have a straightforward interpretation. The cut with \(y=0\) is located in \(x=\sigma\):

sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

xlim_figure = [0.01, 6.0]
ylim_figure = [-2.0, 10.0]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstrom
plt.plot(x, LJ(x, sigma, epsilon))
plt.hlines(0, xlim_figure[0], xlim_figure[1], linestyles='dotted', color='gray')
plt.vlines(sigma._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='red')
plt.text(sigma._value+0.02*xlim_figure[1], 0.7*ylim_figure[1], '$\sigma$', fontsize=14)
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
../../../_images/LJ_Potential_7_0.png

And \(\epsilon\) is the depth of the minimum measured from \(y=0\):

sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

xlim_figure = [0.01, 6.0]
ylim_figure = [-2.0, 10.0]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstrom
plt.plot(x, LJ(x, sigma, epsilon))
plt.hlines(0, xlim_figure[0], xlim_figure[1], linestyles='dotted', color='gray')
plt.hlines(-epsilon._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='red')
plt.annotate(text='', xy=(1.0,0.0), xytext=(1.0,-epsilon._value), arrowprops=dict(arrowstyle='<->'))
plt.text(1.0+0.02*xlim_figure[1], -0.7*epsilon._value, '$\epsilon$', fontsize=14)
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
../../../_images/LJ_Potential_9_0.png

Notice that the LJ potential has physical meaning when \(\epsilon>0\) and \(\sigma>0\) only. Actually, the potential vanishes whether \(\epsilon=0\) or \(\sigma=0\).

The Lennard Jones minimum and the size of the particles

The LJ potential has a single minimum located in \(x_{min}\). Lets equal to \(0\) the first derivative of the potential to find the value of \(x_{min}\):

x, sigma, epsilon = sy.symbols('x sigma epsilon', real=True, positive=True)
V = 4.0*epsilon*((sigma/x)**12-(sigma/x)**6)
gradV = sy.diff(V,x)
roots=sy.solve(gradV, x)
x_min = roots[0]

x_min
\[\displaystyle 1.12246204830937 \sigma\]

The minimum is then located in:

(2)\[ x_{min} = 2^{1/6} \sigma \]

where the potential takes the value:

(3)\[ V(x_{min}) = -\epsilon \]
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

x_min = 2**(1/6)*sigma
y_min = -epsilon

xlim_figure = [x_min._value-0.4, x_min._value+0.4]
ylim_figure = [y_min._value-0.1, y_min._value+0.5]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstroms
plt.plot(x, LJ(x, sigma, epsilon))
plt.hlines(y_min._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='gray')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
../../../_images/LJ_Potential_16_0.png

This way two particles in the equilibrium position will be placed at a \(2^{1/6} \sigma\) distance. The potential is thereby modeling two “soft spheres” atracting each other very lightly. Their radii, given that both particles are equal, are equal to \(r\):

(4)\[ r = \frac{1}{2} x_{min} = 2^{-5/6} \sigma \]

And we say these spheres are “soft” because their volume is not limited by a hard-wall potential, they can penetrate each other suffering a not infinite repulsive force.

Time period of the small harmonic oscillations around the minimum

If we want to perform a molecular simulation of this two particles we should wonder how big the integrator timestep must be. To answer this question we can study the harmonic approximation around the minimum. Lets calculate the time period, \(\tau\), of a small harmonic oscillation around the minimum:

x, sigma, epsilon = sy.symbols('x sigma epsilon', real=True, positive=True)
V = 4.0*epsilon*((sigma/x)**12-(sigma/x)**6)
gradV = sy.diff(V,x)
grad2V = sy.diff(V,x,x)

x_min = sy.solve(gradV,x)[0]
k_harm = grad2V.subs(x, x_min)

k_harm
\[\displaystyle \frac{57.1464378708551 \epsilon}{\sigma^{2}}\]

The harmonic constant of the second degree Taylor polynomial of the LJ potential at \(x=x_{min}\) is then:

(5)\[ k_{harm} = 36·2^{2/3} \frac{\epsilon}{\sigma^2} \]

The oscillation period of a particle with \(m\) mass in an harmonic potential defined by \(\frac{1}{2} k x²\) is:

(6)\[ \tau = 2 \pi \sqrt{ \frac{m}{k}} \]

As such, the period of the small harmonic oscillations around the LJ minimum of particle with \(m\) mass is:

(7)\[ \tau = 2 \pi \sqrt{ \frac{m}{k_{harm}}} = \frac{\pi}{3·2^{1/3}} \sqrt{\frac{m\sigma^2}{\epsilon}} \]

With the mass and parameters taking values of amus, angstroms and kilocalories per mole, the time period is in the order of:

mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

tau = 2*np.pi * np.sqrt(mass/k)

print(tau)

0.5746513694274475 ps

But, is this characteristic time a good threshold for a LJ potential? If the oscillations around the minimum are not small enough, the harmonic potential of the second degree term of the taylor expansion is easily overcome by the sharp left branch of the LJ potential:

sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

x_min = 2**(1/6)*sigma
y_min = -epsilon

xlim_figure = [x_min._value-0.2, x_min._value+0.2]
ylim_figure = [y_min._value-0.1, y_min._value+0.6]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstroms
plt.plot(x, LJ(x, sigma, epsilon))
plt.plot(x, 0.5*k*(x-x_min)**2+y_min)
plt.hlines(y_min._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='gray')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
../../../_images/LJ_Potential_28_0.png

Let’s imagine the following situation. Let a particle be in the harmonic potential at temperature of 300K. Will the particle be more constrained in space than in the well of the LJ potential? Will the particle feel the harmonic potential softer or sharper than the LJ? Lets make some numbers to evaluate if the oscillation time period of the harmonic approximation can be a good time threshold for the integration timestep of a molecular dynamics of the LJ potential.

The standard deviation of an harmonic oscillation with the shape \(\frac{1}{2}k x^2\) in contact with a stochastic thermal bath can be computed as:

(8)\[ \beta = \frac{1}{k_{\rm B} T} \]
(9)\[ Z_x = \int_{-\infty}^{\infty} {\rm e}^{- \beta \frac{1}{2}k x^2} = \sqrt{\frac{2 \pi}{\beta k}} \]
(10)\[ \left< x \right> = \frac{1}{Z_x} \int_{-\infty}^{\infty} x {\rm e}^{-\beta \frac{1}{2}k x^2} = 0 \]
(11)\[ \left< x^2 \right> = \frac{1}{Z_x} \int_{-\infty}^{\infty} x^{2} {\rm e}^{-\beta \frac{1}{2}k x^2} = \frac{1}{Z_x} \sqrt{\frac{2 \pi}{\beta³ k^3}} = \frac{1}{\beta k} \]
(12)\[ {\rm std} = \left( \left< x^2 \right> -\left< x \right>^2 \right)^{1/2} = \sqrt{ \frac{k_{\rm B}T}{k} } \]

This way, in the case of the harmonic potential obtained as the second degree term of the Taylor expansion around the LJ minimum:

mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole
temperature = 300 * unit.kelvin
kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA

k = 36 * 2**(2/3) * epsilon/sigma**2
std = np.sqrt(kB*temperature/k)

x_min = 2**(1/6)*sigma
y_min = -epsilon

xlim_figure = [x_min._value-0.4, x_min._value+0.4]
ylim_figure = [y_min._value-0.1, y_min._value+0.6]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstroms
plt.plot(x, LJ(x, sigma, epsilon))
plt.plot(x, 0.5*k*(x-x_min)**2+y_min)
plt.hlines(y_min._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='gray')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.axvspan(x_min._value - std._value, x_min._value + std._value, alpha=0.2, color='red')
plt.annotate(text='', xy=(x_min._value, y_min._value - 0.5*(y_min._value-ylim_figure[0])),
             xytext=(x_min._value-std._value, y_min._value - 0.5*(y_min._value-ylim_figure[0])),
             arrowprops=dict(arrowstyle='<->'))
plt.text(x_min._value-0.6*std._value, y_min._value - 0.4*(y_min._value-ylim_figure[0]), '$std$', fontsize=14)
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
../../../_images/LJ_Potential_30_0.png

The harmonic potential is too soft as approximation. Its oscillation time used as threshold to choose the integration timestep can yield to numeric problems. Let’s try with a stiffer potential, let’s double the harmonic constant:

mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole
temperature = 300 * unit.kelvin
kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA

k = 36 * 2**(2/3) * epsilon/sigma**2
std = np.sqrt(kB*temperature/k)

x_min = 2**(1/6)*sigma
y_min = -epsilon

xlim_figure = [x_min._value-0.4, x_min._value+0.4]
ylim_figure = [y_min._value-0.1, y_min._value+0.6]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstroms
plt.plot(x, LJ(x, sigma, epsilon))
plt.plot(x, 0.5*k*(x-x_min)**2+y_min)
plt.plot(x, k*(x-x_min)**2+y_min, label='2k_{harm}')
plt.hlines(y_min._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='gray')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.axvspan(x_min._value - std._value, x_min._value + std._value, alpha=0.2, color='red')
plt.annotate(text='', xy=(x_min._value, y_min._value - 0.5*(y_min._value-ylim_figure[0])),
             xytext=(x_min._value-std._value, y_min._value - 0.5*(y_min._value-ylim_figure[0])),
             arrowprops=dict(arrowstyle='<->'))
plt.text(x_min._value-0.6*std._value, y_min._value - 0.4*(y_min._value-ylim_figure[0]), '$std$', fontsize=14)
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
../../../_images/LJ_Potential_32_0.png

Lets take then, as reference, an harmonic potential with constant equal to \(2k_{harm}\) could be a better idea. Lets compute then the new time threshold to choose the integration timestep:

(13)\[ \tau' = 2 \pi \sqrt{ \frac{m}{2k_{harm}}} = \frac{\pi}{3·2^{5/6}} \sqrt{\frac{m\sigma^2}{\epsilon}} = \frac{1}{\sqrt{2}} \tau \]
mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

tau = 2*np.pi * np.sqrt(mass/(2*k))

print(tau)

0.4063398801402841 ps

It is an accepted rule of thumb that the integration timestep must be as large as \(\tau / 10\), being \(\tau\) the oscillation time period of the fastest possible vibration mode. So finally, in this case the integration time step should not be longer than:

mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

tau = 2*np.pi * np.sqrt(mass/(2*k))

print(tau/10.0)

0.04063398801402841 ps