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?
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
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.
