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?

  • 2
    $\begingroup$ 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. $\endgroup$
    – B. Kelly
    Aug 18, 2022 at 2:17

1 Answer 1



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)


To apply this potential, you just need to calculate $$ f = \sum_{i<j} V_{\rm LJ}(R_{ij}) $$ 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

Gradient methods

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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .