Added two new bonded functions.

The first one implements a potential for linear angles without any
discontinuities at 180 degrees. A paper is being written about this. The
potential has [ angles ] type 9 or F_LINEAR_ANGLES internally.
A free energy option has been implemented as well for the linear angles.
The second is an anharmonic polarization function, based on the normal harmonic
term plus a forth order term, as introduced by the Roux group in
J. Chem. Theor. Comput. 6 (2010) p. 774-786. The potential has [ polarization ]
type 2 or F_ANHARM_POL internally.

Change-Id: I0497ed20baf1c1c0b2a6c2354eed77bc9ed75c6d

diff --git a/include/bondf.h b/include/bondf.h
index 9d40a3fa921de8409fb8c5a8e36369847848d955..db28fd376f2490c95ebbd4726ea319af2bf225ed 100644 (file)
--- a/include/bondf.h
+++ b/include/bondf.h
@@ -138,10 +138,10 @@ void do_dih_fup(int i,int j,int k,int l,real ddphi,
  *
  *************************************************************************/
   t_ifunc bonds,g96bonds,morse_bonds,cubic_bonds,FENE_bonds,restraint_bonds;
-  t_ifunc angles,g96angles,cross_bond_bond,cross_bond_angle,urey_bradley,quartic_angles;
+  t_ifunc angles,g96angles,cross_bond_bond,cross_bond_angle,urey_bradley,quartic_angles,linear_angles;
   t_ifunc pdihs,idihs,rbdihs;
   t_ifunc tab_bonds,tab_angles,tab_dihs;
-  t_ifunc polarize,water_pol,thole_pol,angres,angresz,unimplemented;
+  t_ifunc polarize,anharm_polarize,water_pol,thole_pol,angres,angresz,unimplemented;
 
 #ifdef __cplusplus
 }

diff --git a/include/types/idef.h b/include/types/idef.h
index efa74f8059fd1c06e38955204ea18c0b28cccfe0..20180956454e2c7c6be60d5f4077005e85e9059f 100644 (file)
--- a/include/types/idef.h
+++ b/include/types/idef.h
@@ -67,6 +67,7 @@ enum {
   F_RESTRBONDS,
   F_ANGLES, 
   F_G96ANGLES,
+  F_LINEAR_ANGLES,
   F_CROSS_BOND_BONDS,
   F_CROSS_BOND_ANGLES,
   F_UREY_BRADLEY,
@@ -101,6 +102,7 @@ enum {
   F_POLARIZATION,
   F_WATER_POL,
   F_THOLE_POL,
+  F_ANHARM_POL,
   F_POSRES,
   F_DISRES,
   F_DISRESVIOL,
@@ -147,6 +149,7 @@ typedef union
    */
   struct {real a,b,c;                                     } bham;
   struct {real rA,krA,rB,krB;                             } harmonic;
+  struct {real klinA,aA,klinB,aB;                          } linangle;
   struct {real lowA,up1A,up2A,kA,lowB,up1B,up2B,kB;        } restraint;
   /* No free energy supported for cubic bonds, FENE, WPOL or cross terms */ 
   struct {real b0,kb,kcub;                                 } cubic;
@@ -156,6 +159,7 @@ typedef union
   struct {real theta,ktheta,r13,kUB;                       } u_b;
   struct {real theta,c[5];                                 } qangle; 
   struct {real alpha;                                      } polarize;
+  struct {real alpha,drcut,khyp;                           } anharm_polarize;
   struct {real al_x,al_y,al_z,rOH,rHH,rOD;                 } wpol;
   struct {real a,alpha1,alpha2,rfac;                       } thole;
   struct {real c6,c12;                                    } lj;

diff --git a/include/types/nrnb.h b/include/types/nrnb.h
index 14de6c9ef4f9fe7a2bf9bbbb4cfa4c594a6554ca..8f18e18b5d545c69b6a28090c4da6e04a8f835be 100644 (file)
--- a/include/types/nrnb.h
+++ b/include/types/nrnb.h
@@ -100,13 +100,14 @@ enum
     eNR_CONV,                 eNR_SOLVEPME,eNR_NS,      eNR_RESETX,
     eNR_SHIFTX,               eNR_CGCM,                 eNR_FSUM,
     eNR_BONDS,                eNR_G96BONDS,             eNR_FENEBONDS,
-    eNR_TABBONDS,             eNR_RESTRBONDS,
+    eNR_TABBONDS,             eNR_RESTRBONDS,           eNR_LINEAR_ANGLES,
     eNR_ANGLES,               eNR_G96ANGLES,            eNR_QANGLES,
     eNR_TABANGLES,            eNR_PROPER,               eNR_IMPROPER,
     eNR_RB,                   eNR_FOURDIH,              eNR_TABDIHS,
     eNR_DISRES,               eNR_ORIRES,               eNR_DIHRES,
     eNR_POSRES,               eNR_ANGRES,               eNR_ANGRESZ,
     eNR_MORSE,                eNR_CUBICBONDS,           eNR_WALLS,
+    eNR_POLARIZE,             eNR_ANHARM_POL,
     eNR_WPOL,                 eNR_THOLE,                eNR_VIRIAL,
     eNR_UPDATE,               eNR_EXTUPDATE,            eNR_STOPCM,
     eNR_PCOUPL,               eNR_EKIN,                 eNR_LINCS,

diff --git a/src/gmxlib/bondfree.c b/src/gmxlib/bondfree.c
index 273804f7010efec7e8c0c56245f693dae4fddbea..f0e5f181e3095177bfca4ad21f5154a41d1ed99b 100644 (file)
--- a/src/gmxlib/bondfree.c
+++ b/src/gmxlib/bondfree.c
@@ -517,6 +517,63 @@ real polarize(int nbonds,
   return vtot;
 }
 
+real anharm_polarize(int nbonds,
+                     const t_iatom forceatoms[],const t_iparams forceparams[],
+                     const rvec x[],rvec f[],rvec fshift[],
+                     const t_pbc *pbc,const t_graph *g,
+                     real lambda,real *dvdlambda,
+                     const t_mdatoms *md,t_fcdata *fcd,
+                     int *global_atom_index)
+{
+  int  i,m,ki,ai,aj,type;
+  real dr,dr2,fbond,vbond,fij,vtot,ksh,khyp,drcut,ddr,ddr3;
+  rvec dx;
+  ivec dt;
+
+  vtot = 0.0;
+  for(i=0; (i<nbonds); ) {
+    type  = forceatoms[i++];
+    ai    = forceatoms[i++];
+    aj    = forceatoms[i++];
+    ksh   = sqr(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].anharm_polarize.alpha; /* 7*/
+    khyp  = forceparams[type].anharm_polarize.khyp;
+    drcut = forceparams[type].anharm_polarize.drcut;
+    if (debug)
+      fprintf(debug,"POL: local ai = %d aj = %d ksh = %.3f\n",ai,aj,ksh);
+  
+    ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);   /*   3          */
+    dr2  = iprod(dx,dx);                       /*   5          */
+    dr   = dr2*gmx_invsqrt(dr2);                       /*  10          */
+
+    *dvdlambda += harmonic(ksh,ksh,0,0,dr,lambda,&vbond,&fbond);  /*  19  */
+
+    if (dr2 == 0.0)
+      continue;
+    
+    if (dr > drcut) {
+        ddr    = dr-drcut;
+        ddr3   = ddr*ddr*ddr;
+        vbond += khyp*ddr*ddr3;
+        fbond -= 4*khyp*ddr3;
+    }
+    fbond *= gmx_invsqrt(dr2);                 /*   6          */
+    vtot  += vbond;/* 1*/
+
+    if (g) {
+      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;
+      fshift[ki][m]+=fij;
+      fshift[CENTRAL][m]-=fij;
+    }
+  }                                    /* 72 TOTAL     */
+  return vtot;
+}
+
 real water_pol(int nbonds,
               const t_iatom forceatoms[],const t_iparams forceparams[],
               const rvec x[],rvec f[],rvec fshift[],
@@ -804,6 +861,76 @@ real angles(int nbonds,
   return vtot;
 }
 
+real linear_angles(int nbonds,
+                   const t_iatom forceatoms[],const t_iparams forceparams[],
+                   const rvec x[],rvec f[],rvec fshift[],
+                   const t_pbc *pbc,const t_graph *g,
+                   real lambda,real *dvdlambda,
+                   const t_mdatoms *md,t_fcdata *fcd,
+                   int *global_atom_index)
+{
+  int  i,m,ai,aj,ak,t1,t2,type;
+  rvec f_i,f_j,f_k;
+  real L1,kA,kB,aA,aB,dr,dr2,va,vtot,a,b,klin,dvdl;
+  ivec jt,dt_ij,dt_kj;
+  rvec r_ij,r_kj,r_ik,dx;
+    
+  L1   = 1-lambda;
+  vtot = 0.0;
+  dvdl = 0.0;
+  for(i=0; (i<nbonds); ) {
+    type = forceatoms[i++];
+    ai   = forceatoms[i++];
+    aj   = forceatoms[i++];
+    ak   = forceatoms[i++];
+    
+    kA = forceparams[type].linangle.klinA;
+    kB = forceparams[type].linangle.klinB;
+    klin = L1*kA + lambda*kB;
+    
+    aA   = forceparams[type].linangle.aA;
+    aB   = forceparams[type].linangle.aB;
+    a    = L1*aA+lambda*aB;
+    b    = 1-a;
+    
+    t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
+    t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
+    rvec_sub(r_ij,r_kj,r_ik);
+    
+    dr2 = 0;
+    for(m=0; (m<DIM); m++) 
+    {
+        dr     = - a * r_ij[m] - b * r_kj[m];
+        dr2   += dr*dr;
+        dx[m]  = dr;
+        f_i[m] = a*klin*dr;
+        f_k[m] = b*klin*dr;
+        f_j[m] = -(f_i[m]+f_k[m]);
+        f[ai][m] += f_i[m];
+        f[aj][m] += f_j[m];
+        f[ak][m] += f_k[m];
+    }
+    va    = 0.5*klin*dr2;
+    dvdl += 0.5*(kB-kA)*dr2 + klin*(aB-aA)*iprod(dx,r_ik); 
+            
+    vtot += va;
+    
+    if (g) {
+        copy_ivec(SHIFT_IVEC(g,aj),jt);
+      
+        ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
+        ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
+        t1=IVEC2IS(dt_ij);
+        t2=IVEC2IS(dt_kj);
+    }
+    rvec_inc(fshift[t1],f_i);
+    rvec_inc(fshift[CENTRAL],f_j);
+    rvec_inc(fshift[t2],f_k);
+  }                                           /* 57 TOTAL      */
+  *dvdlambda = dvdl;
+  return vtot;
+}
+
 real urey_bradley(int nbonds,
                  const t_iatom forceatoms[],const t_iparams forceparams[],
                  const rvec x[],rvec f[],rvec fshift[],

diff --git a/src/gmxlib/ifunc.c b/src/gmxlib/ifunc.c
index d0444c565be5bdf7e869df92cbfbad25af5bc675..13f0cf166811344cccaba9e5a3644a54c9c5fe2c 100644 (file)
--- a/src/gmxlib/ifunc.c
+++ b/src/gmxlib/ifunc.c
@@ -100,6 +100,7 @@ const t_interaction_function interaction_function[F_NRE]=
   def_bonded  ("RESTRAINTPOT", "Restraint Pot.", 2, 4, 4,  eNR_RESTRBONDS,  restraint_bonds ),
   def_angle   ("ANGLES",   "Angle",           3, 2, 2,  eNR_ANGLES, angles        ),
   def_angle   ("G96ANGLES","G96Angle",        3, 2, 2,  eNR_ANGLES, g96angles     ),
+  def_angle   ("LINEAR_ANGLES", "Lin. Angle", 3, 2, 2,  eNR_LINEAR_ANGLES, linear_angles ),
   def_bonded  ("CROSS_BOND_BOND", "Bond-Cross", 3, 3, 0,0,          cross_bond_bond ),
   def_bonded  ("CROSS_BOND_ANGLE","BA-Cross",   3, 4, 0,0,          cross_bond_angle ),
   def_angle   ("UREY_BRADLEY","U-B",          3, 4, 0,  0,          urey_bradley ),
@@ -134,6 +135,7 @@ const t_interaction_function interaction_function[F_NRE]=
   def_bondnb  ("POLARIZATION", "Polarization",2, 1, 0,  0,          polarize      ),
   def_bonded  ("WATERPOL", "Water Pol.",      5, 6, 0,  eNR_WPOL,   water_pol     ),
   def_bonded  ("THOLE",    "Thole Pol.",      4, 3, 0,  eNR_THOLE,  thole_pol     ),
+  def_bondnb  ("ANHARM_POL", "Anharm. Pol.",2, 3, 0, 0,          anharm_polarize      ),
   def_bonded  ("POSRES",   "Position Rest.",  1, 3, 3,  eNR_POSRES, unimplemented ),
   def_bonded  ("DISRES",   "Dis. Rest.",      2, 6, 0,  eNR_DISRES, ta_disres     ),
   def_nofc    ("DRVIOL",   "D.R.Viol. (nm)"                                       ),    

diff --git a/src/gmxlib/nrnb.c b/src/gmxlib/nrnb.c
index e99331eda7a85280be868ef0f6e13da1141103e1..85567cde38a4884b14c96f661092ffc71adb47f8 100644 (file)
--- a/src/gmxlib/nrnb.c
+++ b/src/gmxlib/nrnb.c
@@ -214,6 +214,7 @@ static const t_nrnb_data nbdata[eNRNB] = {
     { "FENE Bonds",                    58  },
     { "Tab. Bonds",                    62  },
     { "Restraint Potential",           86  },
+    { "Linear Angles",                 57  },
     { "Angles",                        168 },
     { "G96Angles",                     150 },
     { "Quartic Angles",                160 },
@@ -232,6 +233,8 @@ static const t_nrnb_data nbdata[eNRNB] = {
     { "Morse Potent.",                 58  },
     { "Cubic Bonds",                   54  },
     { "Walls",                         31  },
+    { "Polarization",                  59  },
+    { "Anharmonic Polarization",       72  },
     { "Water Pol.",                    62  },
     { "Thole Pol.",                    296 },
     { "Virial",                        18  },

diff --git a/src/gmxlib/tpxio.c b/src/gmxlib/tpxio.c
index ec6863aaaa84ab2d89c534e7eb1f3c9ddc4e616f..2edea91884c4f1d36c6f34e01d4983a782fbd926 100644 (file)
--- a/src/gmxlib/tpxio.c
+++ b/src/gmxlib/tpxio.c
@@ -64,7 +64,7 @@
 #include "mtop_util.h"
 
 /* This number should be increased whenever the file format changes! */
-static const int tpx_version = 75;
+static const int tpx_version = 76;
 
 /* This number should only be increased when you edit the TOPOLOGY section
  * of the tpx format. This way we can maintain forward compatibility too
@@ -126,6 +126,7 @@ static const t_ftupd ftupd[] = {
   { 43, F_TABBONDS          },
   { 43, F_TABBONDSNC        },
   { 70, F_RESTRBONDS        },
+  { 76, F_LINEAR_ANGLES     },
   { 30, F_CROSS_BOND_BONDS  },
   { 30, F_CROSS_BOND_ANGLES },
   { 30, F_UREY_BRADLEY      },
@@ -161,6 +162,7 @@ static const t_ftupd ftupd[] = {
   { 69, F_VTEMP             },
   { 66, F_PDISPCORR         },
   { 54, F_DHDL_CON          },
+  { 76, F_ANHARM_POL        }
 };
 #define NFTUPD asize(ftupd)
 
@@ -1031,6 +1033,12 @@ void do_iparams(t_fileio *fio, t_functype ftype,t_iparams *iparams,
       iparams->pdihs.cpB  = iparams->pdihs.cpA;
     }
     break;
+  case F_LINEAR_ANGLES:
+    gmx_fio_do_real(fio,iparams->linangle.klinA);
+    gmx_fio_do_real(fio,iparams->linangle.aA);
+    gmx_fio_do_real(fio,iparams->linangle.klinB);
+    gmx_fio_do_real(fio,iparams->linangle.aB);
+    break;
   case F_FENEBONDS:
     gmx_fio_do_real(fio,iparams->fene.bm);
     gmx_fio_do_real(fio,iparams->fene.kb);
@@ -1094,6 +1102,11 @@ void do_iparams(t_fileio *fio, t_functype ftype,t_iparams *iparams,
   case F_POLARIZATION:
     gmx_fio_do_real(fio,iparams->polarize.alpha);
     break;
+  case F_ANHARM_POL:
+    gmx_fio_do_real(fio,iparams->anharm_polarize.alpha);
+    gmx_fio_do_real(fio,iparams->anharm_polarize.drcut);
+    gmx_fio_do_real(fio,iparams->anharm_polarize.khyp);
+    break;
   case F_WATER_POL:
     if (file_version < 31) 
       gmx_fatal(FARGS,"Old tpr files with water_polarization not supported. Make a new.");

diff --git a/src/gmxlib/txtdump.c b/src/gmxlib/txtdump.c
index 4a362578cce8b959dc28c58f1eddb920e0e501f7..ba1ce00e95b67f0b487dab56a160f3ef8a6aef0b 100644 (file)
--- a/src/gmxlib/txtdump.c
+++ b/src/gmxlib/txtdump.c
@@ -770,6 +770,11 @@ void pr_iparams(FILE *fp,t_functype ftype,t_iparams *iparams)
            iparams->cross_ba.r1e,iparams->cross_ba.r2e,
            iparams->cross_ba.r3e,iparams->cross_ba.krt);
     break;
+  case F_LINEAR_ANGLES:
+    fprintf(fp,"klinA=%15.8e, aA=%15.8e, klinB=%15.8e, aB=%15.8e\n",
+            iparams->linangle.klinA,iparams->linangle.aA,
+            iparams->linangle.klinB,iparams->linangle.aB);
+    break;
   case F_UREY_BRADLEY:
     fprintf(fp,"theta=%15.8e, ktheta=%15.8e, r13=%15.8e, kUB=%15.8e\n",
            iparams->u_b.theta,iparams->u_b.ktheta,iparams->u_b.r13,iparams->u_b.kUB);
@@ -823,6 +828,12 @@ void pr_iparams(FILE *fp,t_functype ftype,t_iparams *iparams)
   case F_POLARIZATION:
     fprintf(fp,"alpha=%15.8e\n",iparams->polarize.alpha);
     break;
+  case F_ANHARM_POL:
+    fprintf(fp,"alpha=%15.8e drcut=%15.8e khyp=%15.8e\n",
+            iparams->anharm_polarize.alpha,
+            iparams->anharm_polarize.drcut,
+            iparams->anharm_polarize.khyp);
+    break;
   case F_THOLE_POL:
     fprintf(fp,"a=%15.8e, alpha1=%15.8e, alpha2=%15.8e, rfac=%15.8e\n",
            iparams->thole.a,iparams->thole.alpha1,iparams->thole.alpha2,

diff --git a/src/kernel/convparm.c b/src/kernel/convparm.c
index fdef7ea8155732c56a97cc607e1453081ac8feb1..1898c24945d3b9f63fe7694fd982e2554d4756a5 100644 (file)
--- a/src/kernel/convparm.c
+++ b/src/kernel/convparm.c
@@ -159,6 +159,12 @@ static void assign_param(t_functype ftype,t_iparams *newparam,
     for(i=0; i<5; i++)
       newparam->qangle.c[i]=old[i+1];
     break;
+  case F_LINEAR_ANGLES:
+    newparam->linangle.klinA = old[0];
+    newparam->linangle.aA    = old[1];
+    newparam->linangle.klinB = old[2];
+    newparam->linangle.aB    = old[3];
+    break;
   case F_ANGLES:
   case F_BONDS:
   case F_HARMONIC:
@@ -183,6 +189,11 @@ static void assign_param(t_functype ftype,t_iparams *newparam,
   case F_POLARIZATION:
     newparam->polarize.alpha = old[0];
     break;
+  case F_ANHARM_POL:
+    newparam->anharm_polarize.alpha = old[0];
+    newparam->anharm_polarize.drcut = old[1];
+    newparam->anharm_polarize.khyp  = old[2];
+    break;
   case F_WATER_POL:
     newparam->wpol.al_x   =old[0];
     newparam->wpol.al_y   =old[1];

diff --git a/src/kernel/topdirs.c b/src/kernel/topdirs.c
index 14dd59fc00a7388bf45795a81912836e442d6509..fe5e9acf508a979aba5176ea11b09f35f4871317 100644 (file)
--- a/src/kernel/topdirs.c
+++ b/src/kernel/topdirs.c
@@ -94,6 +94,8 @@ int ifunc_index(directive d,int type)
       return F_QUARTIC_ANGLES;
     case 8:
       return F_TABANGLES;
+    case 9:
+      return F_LINEAR_ANGLES;
     default:
       gmx_fatal(FARGS,"Invalid angle type %d",type);
       break;
@@ -193,7 +195,15 @@ int ifunc_index(directive d,int type)
     }
     break;
   case d_polarization:
-    return F_POLARIZATION;
+    switch (type) {
+    case 1:
+      return F_POLARIZATION;
+    case 2:
+      return F_ANHARM_POL;
+    default:
+      gmx_fatal(FARGS,"Invalid polarization type %d",type);
+    }
+    break;
   case d_thole_polarization:
     return F_THOLE_POL;
   case d_water_polarization:

diff --git a/src/mdlib/shellfc.c b/src/mdlib/shellfc.c
index 91d11ecd79b8a44e63d3b424cd5e2437c95ca450..3295cffbf0cc555a9054c72a14bae7089d72610b 100644 (file)
--- a/src/mdlib/shellfc.c
+++ b/src/mdlib/shellfc.c
@@ -200,7 +200,7 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog,
   int         i,j,nmol,type,mb,mt,a_offset,cg,mol,ftype,nra;
   real        qS,alpha;
   int         aS,aN=0; /* Shell and nucleus */
-  int         bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_WATER_POL };
+  int         bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_ANHARM_POL, F_WATER_POL };
 #define NBT asize(bondtypes)
   t_iatom     *ia;
   gmx_mtop_atomloop_block_t aloopb;
@@ -304,6 +304,7 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog,
          case F_HARMONIC:
          case F_CUBICBONDS:
          case F_POLARIZATION:
+         case F_ANHARM_POL:
            if (atom[ia[1]].ptype == eptShell) {
              aS = ia[1];
              aN = ia[2];
@@ -361,6 +362,7 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog,
              shell[nsi].k    += ffparams->iparams[type].cubic.kb;
              break;
            case F_POLARIZATION:
+           case F_ANHARM_POL:
              if (qS != atom[aS].qB)
                gmx_fatal(FARGS,"polarize can not be used with qA != qB");
              shell[nsi].k    += sqr(qS)*ONE_4PI_EPS0/