def lennard_jones(rsq, epsilon, sigma):
    """
    Lennard-Jones potential.

    u(r) = 4 * epsilon * [(r/sigma)^12 - (r/sigma)^6]
    """
    sigsq = sigma ** 2
    u = 4 * epsilon * ((sigsq / rsq) ** 6 - (sigsq / rsq) ** 3)
    w = 24 * epsilon * (2 * (sigsq / rsq) ** 6 - (sigsq / rsq) ** 3) / rsq
    h = 96 * epsilon * (7 * (sigsq / rsq) ** 6 - 2 * (sigsq / rsq) ** 3
        ) / rsq ** 2
    return u, w, h
