[gmx-users] GROMACS and FENE
Michael Patra
patra at lorentz.leidenuniv.nl
Mon Sep 1 20:55:01 CEST 2003
> > FENE is primarily for the bead-spring-type polymer model. It is a bonded
> > interaction. Explicitly, its standard form is:
> >
> > E = -0.5 K R0^2 * ln[1 - (r/R0)^2]
> >
> > r = distance (btw two connected beads)
> >
> > coeff1 = K (energy/distance^2)
> > coeff2 = R0 (distance)
> Ok, that is easy enough. I suggest you try to implement that e.g.
> instead of the Morse potential. You will then only ave to edit a single
> routine (morse in bondfree.c) and (ab)use the Morse parameters. If you
> have a working function (with energy conservation that is) then you are
> welcome to send me the source code for inclusion in a future release of
> gromacs.
Actually, it is not that easy at all. At some point in time I was implementing
the FENE spring into Gromacs. You can find the diff file at the end of this mail
(without warranty since I introduced also quite a few other changes into the
source code, and now had to extract only the changes relating to the FENE
spring).
As can be seen from the equation above, the force becomes infinite as r->R0.
From time to time (in my experiments: once every 100000 or so time steps) two
bonded particles would become so separated that the resulting force would make
the simulation explode. And even when the simulation runs normally, I guess I
don't even need to comment on the quality of energy conservation. In my
personal opinion, the FENE spring is useless for efficient simulations.
The only work around that I could find, was to cutoff the forces at some maximum
value. This is the reason why you will find two different FENE springs in the
code. Using the cutoff-ed version allows to use a similar time step as when
using harmonic springs. Unfortunately I have to say that the simulation results
were surprisingly much dependent on the precise value of the cutoff. In my
personal opinion, using the cutoff FENE spring therefore is somewhat
questionable.
-----------------------------------------------------------------------------
Dr. Michael Patra
Biophysics and Statistical Mechanics Group
Laboratory of Computational Engineering
Helsinki University of Technology
P.O.Box 9203 (Tekniikantie 14)
FIN-02015 HUT
Finland
Phone: +358-9-451 5724
Fax: +358-9-451 4830
-----------------------------------------------------------------------------
diff -ur orig/gromacs-3.0.5/include/bondf.h gromacs-3.0.5/include/bondf.h
--- orig/gromacs-3.0.5/include/bondf.h Fri Oct 26 15:22:00 2001
+++ gromacs-3.0.5/include/bondf.h Fri Jan 11 22:46:01 2002
@@ -102,7 +102,7 @@
* Bonded force functions
*
*************************************************************************/
- extern t_ifunc bonds,g96bonds,morsebonds,cubicbonds;
+ extern t_ifunc bonds,g96bonds,morsebonds,cubicbonds,fenebonds,fene2bonds;
extern t_ifunc angles,g96angles;
extern t_ifunc pdihs,idihs,rbdihs;
extern t_ifunc water_pol,posres,angres,angresz,do_14,unimplemented;
diff -ur orig/gromacs-3.0.5/include/types/idef.h gromacs-3.0.5/include/types/idef.h
--- orig/gromacs-3.0.5/include/types/idef.h Fri Oct 26 15:22:00 2001
+++ gromacs-3.0.5/include/types/idef.h Fri Jan 11 22:47:21 2002
@@ -57,6 +57,8 @@
F_CUBICBONDS,
F_CONNBONDS,
F_HARMONIC,
+ F_FENE,
+ F_FENE2,
F_ANGLES,
F_G96ANGLES,
F_PDIHS,
@@ -123,6 +125,7 @@
*/
struct {real doh,dhh; } settle;
struct {real b0,cb,beta; } morse;
+ struct {real k,R; } fene;
struct {real pos0[DIM],fc[DIM]; } posres;
struct {real rbc[NR_RBDIHS]; } rbdihs;
struct {real a,b,c,d,e,f; } dummy;
diff -ur orig/gromacs-3.0.5/include/types/nrnb.h gromacs-3.0.5/include/types/nrnb.h
--- orig/gromacs-3.0.5/include/types/nrnb.h Fri Oct 26 15:22:00 2001
+++ gromacs-3.0.5/include/types/nrnb.h Fri Jan 11 22:48:02 2002
@@ -88,6 +88,7 @@
eNR_PSHAKEINITLD, eNR_PSHAKEINITMD, eNR_PSHAKE,
eNR_DUM2, eNR_DUM3, eNR_DUM3FD,
eNR_DUM3FAD, eNR_DUM3OUT, eNR_DUM4FD,
+ eNR_FENE, eNR_FENE2,
eNRNB
};
diff -ur orig/gromacs-3.0.5/src/gmxlib/bondfree.c gromacs-3.0.5/src/gmxlib/bondfree.c
--- orig/gromacs-3.0.5/src/gmxlib/bondfree.c Fri Oct 26 15:21:56 2001
+++ gromacs-3.0.5/src/gmxlib/bondfree.c Fri Jan 11 23:20:02 2002
@@ -253,6 +253,118 @@
}
+
+/* FENE potential (finite extensible spring)
+ Implementation by Michael Patra <patra at lorentz.leidenuniv.nl>
+
+ This describes a spring with stiffness k that cannot be extended to
+ a length r of more than R. The potential are force are hence
+ V(r) = - 1/2 k R^2 ln( 1 - r^2/R^2 )
+ F(r) = - k r R^2 / ( R^2 - r^2 )
+ For fene2bonds everything is cutoff for r>0.99*R, allowing a significantly
+ larger integration time step.
+*/
+
+real fenebonds(int nbonds,
+ t_iatom forceatoms[],t_iparams forceparams[],
+ rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
+ matrix box,real lambda,real *dvdl,
+ t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
+{
+ const real half=0.5;
+ const real one=1.0;
+ real dr,dr2,fbond,vbond,fij,k,R,Rsq,vtot;
+ rvec dx;
+ int i,m,ki,type,ai,aj;
+ ivec dt;
+
+ vtot = 0.0;
+ for(i=0; (i<nbonds); ) {
+ type = forceatoms[i++];
+ ai = forceatoms[i++];
+ aj = forceatoms[i++];
+
+ k = forceparams[type].fene.k;
+ R = forceparams[type].fene.R;
+ Rsq = R*R;
+
+ pbc_rvec_sub(x[ai],x[aj],dx); /* 3 */
+ dr2 = iprod(dx,dx); /* 5 */
+ dr = sqrt(dr2); /* 10 */
+
+ vbond = - half*k*Rsq*log( one - dr2/Rsq ); /* 15 */
+ fbond = -k*dr*Rsq / ( Rsq - dr2 ); /* 14 */
+
+ vtot += vbond; /* 1 */
+
+ ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
+ ki=IVEC2IS(dt);
+
+ for (m=0; (m<DIM); m++) { /* 15 */
+ fij=fbond*dx[m];
+ f[ai][m]+=fij;
+ f[aj][m]-=fij;
+ fr->fshift[ki][m]+=fij;
+ fr->fshift[CENTRAL][m]-=fij;
+ }
+ }
+ return vtot; /* Total: 63 */
+}
+
+real fene2bonds(int nbonds,
+ t_iatom forceatoms[],t_iparams forceparams[],
+ rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
+ matrix box,real lambda,real *dvdl,
+ t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
+{
+ const real half=0.5;
+ const real one=1.0;
+ const real ninenine=0.99; /* cutoff for the r->R */
+ real dr,dr2,fbond,vbond,fij,k,R,Rsq,vtot;
+ rvec dx;
+ int i,m,ki,type,ai,aj;
+ ivec dt;
+
+ vtot = 0.0;
+ for(i=0; (i<nbonds); ) {
+ type = forceatoms[i++];
+ ai = forceatoms[i++];
+ aj = forceatoms[i++];
+
+ k = forceparams[type].fene.k;
+ R = forceparams[type].fene.R;
+ Rsq = R*R;
+
+ pbc_rvec_sub(x[ai],x[aj],dx); /* 3 */
+ dr2 = iprod(dx,dx); /* 5 */
+ dr = sqrt(dr2); /* 10 */
+ if( dr > ninenine*R ) { /* 1 */
+ dr = ninenine*R;
+ dr2 = dr * dr;
+ }
+
+ vbond = - half*k*Rsq*log( one - dr2/Rsq ); /* 15 */
+ fbond = -k*dr*Rsq / ( Rsq - dr2 ); /* 14 */
+
+ vtot += vbond; /* 1 */
+
+ ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
+ ki=IVEC2IS(dt);
+
+ for (m=0; (m<DIM); m++) { /* 15 */
+ fij=fbond*dx[m];
+ f[ai][m]+=fij;
+ f[aj][m]-=fij;
+ fr->fshift[ki][m]+=fij;
+ fr->fshift[CENTRAL][m]-=fij;
+ }
+ }
+ return vtot; /* Total: 64 */
+}
+
+
+
+
real bonds(int nbonds,
t_iatom forceatoms[],t_iparams forceparams[],
rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
diff -ur orig/gromacs-3.0.5/src/gmxlib/ifunc.c gromacs-3.0.5/src/gmxlib/ifunc.c
--- orig/gromacs-3.0.5/src/gmxlib/ifunc.c Fri Oct 26 15:21:57 2001
+++ gromacs-3.0.5/src/gmxlib/ifunc.c Fri Jan 11 23:41:16 2002
@@ -73,6 +73,8 @@
def_bond ("MORSE", "Morse", 2, 3, 0, eNR_MORSE, morsebonds),
def_bond ("CUBICBONDS","Cubic Bonds", 2, 3, 0, eNR_CUBICBONDS, cubicbonds),
def_bondnb ("CONNBONDS","Connect Bonds", 2, 0, 0, 0, unimplemented),
+ def_bond ("FENE", "FENE", 2, 2, 2, eNR_FENE, fenebonds),
+ def_bond ("FENE", "FENE", 2, 2, 2, eNR_FENE2, fene2bonds),
def_bonded ("HARMONIC", "Harmonic Pot.", 2, 2, 2, eNR_BONDS, bonds ),
def_angle ("ANGLES", "Angle", 3, 2, 2, eNR_ANGLES, angles ),
def_angle ("G96ANGLES","G96Angle", 3, 2, 2, eNR_ANGLES, g96angles),
diff -ur orig/gromacs-3.0.5/src/gmxlib/mshift.c gromacs-3.0.5/src/gmxlib/mshift.c
--- orig/gromacs-3.0.5/src/gmxlib/mshift.c Fri Oct 26 15:21:57 2001
+++ gromacs-3.0.5/src/gmxlib/mshift.c Fri Jan 11 23:05:17 2002
@@ -196,6 +196,7 @@
if ((tp == F_BONDS) || (tp == F_G96BONDS) ||
(tp == F_MORSE) || (tp == F_SHAKE) ||
(tp == F_CUBICBONDS) || (tp == F_CONNBONDS) ||
+ (tp == F_FENE) || (tp == F_FENE2) ||
(tp == F_HARMONIC) ||
(interaction_function[tp].flags & IF_DUMMY))
nbond[iaa]++;
diff -ur orig/gromacs-3.0.5/src/gmxlib/nrnb.c gromacs-3.0.5/src/gmxlib/nrnb.c
--- orig/gromacs-3.0.5/src/gmxlib/nrnb.c Fri Oct 26 15:21:56 2001
+++ gromacs-3.0.5/src/gmxlib/nrnb.c Fri Jan 11 23:09:05 2002
@@ -179,7 +179,9 @@
{ "Dummy3fd", 95 },
{ "Dummy3fad", 176 },
{ "Dummy3out", 87 },
- { "Dummy4fd", 110 }
+ { "Dummy4fd", 110 },
+ { "FENE", 64 },
+ { "FENE", 63 }
};
void init_nrnb(t_nrnb *nrnb)
@@ -363,6 +365,11 @@
do_real(iparams->cubic.kcub);
break;
case F_CONNBONDS:
+ break;
+ case F_FENE:
+ case F_FENE2:
+ do_real(iparams->fene.k);
+ do_real(iparams->fene.R);
break;
case F_WPOL:
do_real(iparams->wpol.kx);
diff -ur orig/gromacs-3.0.5/src/gmxlib/txtdump.c gromacs-3.0.5/src/gmxlib/txtdump.c
--- orig/gromacs-3.0.5/src/gmxlib/txtdump.c Fri Oct 26 15:21:56 2001
+++ gromacs-3.0.5/src/gmxlib/txtdump.c Sun Jan 13 23:46:21 2002
@@ -435,6 +436,14 @@
break;
case F_CONNBONDS:
fprintf(fp,"\n");
+ break;
+ case F_FENE:
+ fprintf(fp,"k=%15.8e, R=%15.8e\n",
+ iparams->fene.k,iparams->fene.R);
+ break;
+ case F_FENE2:
+ fprintf(fp,"k=%15.8e, R=%15.8e (with cutoff)\n",
+ iparams->fene.k,iparams->fene.R);
break;
case F_WPOL:
fprintf(fp,"kx=%15.8e, ky=%15.8e, kz=%15.8e, rOH=%9.6f, rHH=%9.6f, rOD=%9.6f\n",
diff -ur orig/gromacs-3.0.5/src/kernel/convparm.c gromacs-3.0.5/src/kernel/convparm.c
--- orig/gromacs-3.0.5/src/kernel/convparm.c Fri Oct 26 15:21:56 2001
+++ gromacs-3.0.5/src/kernel/convparm.c Fri Jan 11 23:51:09 2002
@@ -92,6 +92,11 @@
break;
case F_CONNBONDS:
break;
+ case F_FENE:
+ case F_FENE2:
+ new->fene.k =old[0];
+ new->fene.R =old[1];
+ break;
case F_WPOL:
new->wpol.kx =old[0];
new->wpol.ky =old[1];
diff -ur orig/gromacs-3.0.5/src/kernel/topdirs.c gromacs-3.0.5/src/kernel/topdirs.c
--- orig/gromacs-3.0.5/src/kernel/topdirs.c Fri Oct 26 15:21:56 2001
+++ gromacs-3.0.5/src/kernel/topdirs.c Fri Jan 11 23:13:26 2002
@@ -60,6 +60,10 @@
return F_CONNBONDS;
else if (type == 6)
return F_HARMONIC;
+ else if (type == 7)
+ return F_FENE;
+ else if (type == 8)
+ return F_FENE2;
else
fatal_error(0,"Invalid bond type %d",type);
case d_angles:
diff -ur orig/gromacs-3.0.5/src/kernel/toputil.c gromacs-3.0.5/src/kernel/toputil.c
--- orig/gromacs-3.0.5/src/kernel/toputil.c Fri Oct 26 15:21:56 2001
+++ gromacs-3.0.5/src/kernel/toputil.c Fri Jan 11 23:14:04 2002
@@ -248,6 +248,12 @@
case F_HARMONIC:
f = 5;
break;
+ case F_FENE:
+ f = 7;
+ break;
+ case F_FENE2:
+ f = 8;
+ break;
case F_PDIHS:
case F_RBDIHS:
bDih=TRUE;
diff -ur orig/gromacs-3.0.5/src/ngmx/manager.c gromacs-3.0.5/src/ngmx/manager.c
--- orig/gromacs-3.0.5/src/ngmx/manager.c Fri Oct 26 15:21:56 2001
+++ gromacs-3.0.5/src/ngmx/manager.c Fri Jan 11 23:15:10 2002
@@ -117,6 +117,8 @@
add_bonds(man,idef->functype,&idef->il[F_MORSE],bB);
add_bonds(man,idef->functype,&idef->il[F_CUBICBONDS],bB);
add_bonds(man,idef->functype,&idef->il[F_CONNBONDS],bB);
+ add_bonds(man,idef->functype,&idef->il[F_FENE],bB);
+ add_bonds(man,idef->functype,&idef->il[F_FENE2],bB);
add_bonds(man,idef->functype,&idef->il[F_SHAKE],bB);
add_bonds(man,idef->functype,&idef->il[F_SETTLE],bB);
}
More information about the gromacs.org_gmx-users
mailing list