From: Jim Phillips (jim_at_ks.uiuc.edu)
Date: Mon Jan 29 2007 - 16:55:40 CST
Hi,
The interpolation tables of U(r^2) are used to calculate dU/dr^2. This is
then used to calculate forces without calling sqrt() with the formula:
force = -dU/dr = -2 r dU/dr^2
The force vector is calculated simply by multiplying the vector between
the two interacting atoms by -2 dU/dr^2. Note that to save operations the
factor of 2 is already included in some of the table coefficients.
If the actual force was interpolated from the table it would be necessary
to call sqrt() to calculate the unit vector between atoms, which is very
expensive compared to multiplication and addition.
-Jim
On Mon, 29 Jan 2007, Floris Buelens wrote:
> 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?
>
> Thanks in advance!
>
> 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
> fast_gradient = -1.0/r2
>
>
>
>
>
>
>
>
> ___________________________________________________________
> New Yahoo! Mail is the ultimate force in competitive emailing. Find out more at the Yahoo! Mail Championships. Plus: play games and win prizes.
> http://uk.rd.yahoo.com/evt=44106/*http://mail.yahoo.net/uk
>
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:44:20 CST