[gmx-developers] introduce effective potential

[Nguyen Hoang Phuong]

Dear All,

I would like to replace the dihedral potential energy function (formed
by 4 atoms i,j,k,l) by the effective potential

U_eff(\theta_{ijkl}) = KT*q/(q-1)*log(1 + KT*(q-1)*[U(\theta_{ijkl}) + epsilon])

where q and epsilon are the parameters, U(\theta_{ijkl}) is the original
potential. This effective potential is smoother (depends on value of q) than the originial
one and therefore the sampling is enhanced.

This is equivalent to introduce the scale factor:

q/(1 + KT*(q-1)*[U(\theta_{ijkl}) + epsilon])

into the forces which acting on the atoms i,j,k,l

To this end, I have modified to code as follow (in bondfree.c)

/***************************************************
real pdihs(int nbonds,
           t_iatom forceatoms[],t_iparams forceparams[],
           rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
           matrix box,real lambda,real *dvdlambda,
           t_mdatoms *md,int ngrp,real egnb[],real egcoul[],
           t_fcdata *fcd)
{
  real scale; //PHUONG

  int  i,type,ai,aj,ak,al;
  rvec r_ij,r_kj,r_kl,m,n;
  real phi,cos_phi,sign,ddphi,vpd,vtot;

  vtot = 0.0;

  for(i=0; (i<nbonds); ) {
    type = forceatoms[i++];
    ai   = forceatoms[i++];
    aj   = forceatoms[i++];
    ak   = forceatoms[i++];

phi=dih_angle(box,x[ai],x[aj],x[ak],x[al],r_ij,r_kj,r_kl,m,n,
                  &cos_phi,&sign);                      /*  84          */

    *dvdlambda += dopdihs(forceparams[type].pdihs.cpA,
                          forceparams[type].pdihs.cpB,
                          forceparams[type].pdihs.phiA,
                          forceparams[type].pdihs.phiB,
                          forceparams[type].pdihs.mult,
                          phi,lambda,&vpd,&ddphi);

    vtot += vpd;
    do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
               f,fr,g,x);                               /* 112          */

/***PHUONG***********************/

   scale =  q/(1 + KT*(q-1)*(vpd + epsilon));

   svmul(scale,f[ai],f[ai]);
   svmul(scale,f[aj],f[aj]);
   svmul(scale,f[ak],f[ak]);
   svmul(scale,f[al],f[al]);

/********************************/

#ifdef DEBUG
    fprintf(debug,"pdih: (%d,%d,%d,%d) cp=%g, phi=%g\n",
            ai,aj,ak,al,cos_phi,phi);
#endif
  } /* 223 TOTAL        */

  return vtot;
}
/******************************************/

My questions are:

1) Is this modification correct?

2) I would like to introduce the parameters q and epsilon into the mdp
file, e.g.

q       = 2
epsilon = 10

What I have to do such that gromacs can recognizes these parameters and
put into the right place in the above code?

I appreciate for any suggestion.

Phuong








implement the method which is use to smooth potential
energy function introfuced by Straub (



On Mon, 5 Apr 2004, Nguyen Hoang Phuong wrote:

>
> Dear Developers,
>
> let consider a dihedral angle formed by 4 atoms: i,j,k,l. The
> dihedral potential energy is V(\theta_{ijkl}). I would like to modify the
> force acting on the atoms i,j,k,l as the following
>
> f'_i = (V/KT)*f_i
> f'_j = (V/KT)*f_j
> f'_k = (V/KT)*f_k
> f'_l = (V/KT)*f_l
>
> where f_i is the standard force (= dV/dr_i).
>
> I had a look at the subroutine do_dih_fup.c of the bondfree.c but I could
> not find the expressions of the forces where I can change. Can anyone
> help me? Thank you very much in advance.
>
> Phuong
>
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.
>