# Lennard Jones Clusters in Python

I would like to explore the global optimisation of Lennard-Jones clusters in Python. I do not come from a physical chemistry background. I would like to be able to calculate the cluster energy pointwise. That is, I would like a python function $$f : X \rightarrow \mathbb{R}$$ which I can use to test my methods on.

Do there exist simple python packages that have this functionality, or how can the clusters be implemented in python?

• Check out the codes/software that uses Jax. github.com/google/jax-md a year or two ago when I first looked at it, it only did Lennard-Jones, but it looks like the author has really put alot of work into it. Commented Aug 18, 2022 at 2:17

# Potential

The Lennard-Jones potential has a simple form: $$V_{\rm LJ} (r) = 4\epsilon \left[ \left( \frac \sigma r \right)^{12} - \left( \frac \sigma r \right)^{6} \right].$$ In python, this could be

def V_LJ(r, epsilon, sigma):
'''Computes the Lennard-Jones potential'''
return 4.0*epsilon*((sigma/r)**12 - (sigma/r)**6)


# Application

To apply this potential, you just need to calculate $$f = \sum_{i where $$R_{ij}$$ is the distance between atoms $$i$$ and $$j$$. Assume that your coordinates are in a numpy array coords = numpy.zeros((3,Natoms)). Then, you can compute $$R_{ij}$$ as Rij = numpy.linalg.norm(coords[i]-coords[j]).

Evaluation of the full potential energy would then be

def LJ_potential(atoms):
'''Computes Lennard-Jones potential for atoms'''
Natoms = atoms.shape[1]
f = 0.0
for i in range(Natoms):
for j in range(i):
Rij = numpy.linalg.norm(coords[i]-coords[j])
f += V_LJ(Rij, epsilon, sigma)
return f


Note that the function also has a simple gradient $$\frac {dV_{\rm LJ}} {dr} (r) = \frac {4\epsilon} {\sigma} \left[ -12 \left( \frac \sigma r \right)^{13} + 6 \left( \frac \sigma r \right)^{7} \right],$$ so you can also use gradient descent methods for optimization, or molecular dynamics methods.