# FAO developers: vdW energies, forces and derivatives

From: Floris Buelens (floris_buelens_at_yahoo.com)
Date: Mon Jan 29 2007 - 05:23:43 CST

Hello NAMD users and developers,

I have a question about the relationship between Lennard-Jones energy and force in NAMD. From an old NAMD programming guide (v1.5, 1998), I read:

vdW energy = (A / r^12) - (B / r^6) - the basic Lennard-Jones potential function;

It's then stated that the force is obtained by differentiating that function:

vdW force = (12A / r^13) - (6B / r^7)

However I can't seem to make this agree with the implementation I find in the 2.6 source... ComputeNonbondedUtil.C pre-calculates tables of r^2 versus various nonbonded energy parameters at startup. For the vdW forces, the calculation and differentiation all appear to be based on r^2, not r - I interpret:

r_2 = 1 / r^2
r_6 = r^(-6) = (r_2) ^ 3;
r_12 = (r_6) ^ 2;
van der Waals energy = A * r_12 - B * r_6

I'm clear on this so far, but then the routine stores vdwa_gradient and vdwb_gradient as:

vdwa_gradient = -6 * r_2 * r_12 = -6 * r^(-14) [this corresponds to dE / dr^2 of E = (r^2)^-6 -> -6*(r^2)^-7]
vdwb_gradient = -3 * r_2 * r_6 = -3 * r^(-8) [this corresponds to dE / dr^2 of E = (r^2)^-3 -> -3(r^2)^-4]

(assuming for convenience that switching function is not active.)

These gradient values are then later (in ComputeNonbondedBase2) multiplied by A and B respectively to give the vdW force.

My understanding and the NAMD programming guide suggest we need to calculate dE / dr to get the force, while as I understand it, NAMD 2.6 calculates dE / dr^2. (12A / r^13) is not the same as (6A / r^14), and (6B / r^7) is not the same as (3B / r^8) - I don't understand the discrepancy. I need to get this straight in order to implement a modified vdW potential for free energy calculations. Can anybody help?

Floris Buelens
Department of Crystallography, Birkbeck College, London

p.s. - for comparison, the equivalent formula for electrostatics follows my expectations:
energy = Cqiqj / r = Cqiqj * r^-1
force = -1 * (Cqiqj) ^ -2 = -Cqiqj / r^2
... from the programming guide, corresponds to:

fast_energy = 1.0/r