[gmx-users] BUG!!!!: FENE bond stretching potential

[Berk Hess]

>From: Jiang ran <ran_jiong_lang2 at hotmail.com>
>Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>To: <gmx-users at gromacs.org>
>Subject: [gmx-users] BUG!!!!: FENE bond stretching potential
>Date: Thu, 28 Jun 2007 23:07:55 +0800
>
>
>Hello people,My professor and I tried to use the FENE potentials to relax a 
>protein from a unstable configuration for publishing paper. However,  these 
>potentials didn't seem to work. Checking both the manual and the code of 
>the GROMACS software package for moleculair dynamics we found some things 
>we do not understand both me as even my professor.First one thing: In file 
>bondfree.c in gmxlib  the  force constant  kb (parameter: real kb) for the 
>FENE potential is not included anywhere in the force neither potential 
>calculation .      bm   = forceparams[type].fene.bm;     kb   = 
>forceparams[type].fene.kb;    bm2 = bm*bm;    omdr2obm2  = one - dr2/bm2;   
>  vbond      = -half*bm2*log(omdr2obm2);    fbond      = 
>-bm/omdr2obm2;Second one next thing:In the following lines:The scalar force 
>is composed over its dimension by dx[m]:fij=bond*dx[m]However in stead of 
>taking a unit normalized vector the original connecting vector is taken. 
>The manual isn't clear about this either, but looking at the equations it 
>doens't make a lot of sense to me and my professor to use the original 
>connecting vector.Greetings from China,Jiang

Indeed the code is wrong (apparently nobody ever tried to use it,
including the person who coded it (myself)).

The potential and force should be multiplied with kb and bm should
not be in the force:
    vbond      = -half*kb*bm2*log(omdr2obm2);
    fbond      = -kb/omdr2obm2;

I also noticed that in the manual there is a minus sign missing
in front of the force.

Thanks for reporting this bug,

Berk.

_________________________________________________________________
Play online games with your friends with Messenger 
http://www.join.msn.com/messenger/overview