lennard jones gradient and energy

jlge.py

def calc_energy_and_gradient(X, Y):
    """ Calculate the Lennard-Jones energy and
        Energy gradient between particles.
        X and Y are lists containing x and y
        coordinates, respectively.
    """

    # Holds the energy and the energy gradient
    energy = 0.0
    gradient = numpy.zeros((2, len(X)))

    # Double loop over all particles
    for i in range(len(X)):
        for j in range(len(X)):
            if i>j:

                # Fortran'esque calculation of inverse distances
                # to 2nd, 6th and 12th power.
                rij2 = (X[i]-X[j])**2 + (Y[i]-Y[j])**2

                irij2 = 1.0 / rij2

                irij6 = irij2*irij2*irij2

                irij12 = irij6*irij6


                # Summation of the Lennard-Jones energy
                energy += irij12 - irij6

                # Summation of the Lennard-Jones energy gradient
                gradient[0][i] += -(X[j]-X[i])*irij2 * (irij12 - 0.5*irij6)
                gradient[0][j] +=  (X[j]-X[i])*irij2 * (irij12 - 0.5*irij6)

                gradient[1][i] += -(Y[j]-Y[i])*irij2 * (irij12 - 0.5*irij6)
                gradient[1][j] +=  (Y[j]-Y[i])*irij2 * (irij12 - 0.5*irij6)


    energy   *=  4.0
    gradient *= 48.0

    return energy, gradient