[gmx-developers] Buckingham potential force calculation correct ?

[Mathias PUETZ]

Hi,

while tuning the force kernels for BlueGene I stumbled across the 
following.

The C (also Fortran) code generated by mknb for calculating the Buckingham 
potential is the following:

  rinvsq = rinv * rinv
  br = cexp2 * rsq * rinv
  Vvdwexp = cexp1 * exp (-br) 
  Vvdw6 = C6 / (rinvsq*rinvsq*rinvsq)
  fscal = (br * Vvdwexp - 6.0 * Vvdw6) * rinvsq;

  fx = dx * fscal
  fy = dy * fscal
  fz = dz * fscal

This does not look right to me, because the br * Vvdwexp part in the 
calculation of
fscal should only be scaled by rinv and not rinvsq as deduced by following 
calculation:

The Buckingham potential is defined by the following form:

  V(r) =   C1 * exp ( - C2 * r) -  C6 / r**6

The force on particle 1 is calculated by the formula:

  Fx = -  dx / r  *  d/dr V(r)
       =  -dx / r  * { - C1*C2*r * exp( -C2 * r)  + 6.0 * C6 / r**7 }
       = -dx *  { 1/r**2  * 6.0 * C6 / r**6  + 1/r * br * C1 * exp (-C2 * 
r) }
       = -dx * { rinvsq * 6.0 * Vvdw6  + rinv * br * Vvdwexp }

      != -dx * rinvsq * { 6.0 * Vvdw6  + br * Vvdwexp }

Am I making a mistake here or is the GROMACS code wrong ?

Viele Grüsse / Best regards,
Dr. Mathias Pütz

IT Specialist for Application Perfomance

Deep Computing - Strategic Growth Business
IBM Systems & Technology Group

e-mail:  mpuetz at de.ibm.com
mobile: + 49-(0)160-7120602
fax:         + 49-(0)6131-84-6660

snailmail:
  IBM Deutschland GmbH
  Department B458
  Hechtsheimer Str. 2 / Building 12
  55131 Mainz
  Germany
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20061002/a4e91ac9/attachment.html>