made non-bonded potential shifts mdp options

Change-Id: I03e500c28ebd677917571cca7a7551523f9eb3af

diff --git a/include/names.h b/include/names.h
index 934ee29f6abedaf3c00abb76639fd2a3a76cc004..a654b3838824288238ba78fec0d1d272417d3a8a 100644 (file)
--- a/include/names.h
+++ b/include/names.h
@@ -57,6 +57,7 @@ extern const char *ens_names[ensNR+1];
 extern const char *ei_names[eiNR+1];
 extern const char *yesno_names[BOOL_NR+1];
 extern const char *bool_names[BOOL_NR+1];
+extern const char *eintmod_names[eintmodNR+1];
 extern const char *eel_names[eelNR+1];
 extern const char *eewg_names[eewgNR+1];
 extern const char *evdw_names[evdwNR+1];
@@ -112,6 +113,7 @@ extern const char *eAdressSITEtype_names[eAdressSITENR+1];
 #define EREFSCALINGTYPE(e) ENUM_NAME(e,erscNR,erefscaling_names)
 #define EBLOCKS(e)     ENUM_NAME(e,ebNR,eblock_names)
 #define EPARAM(e)      ENUM_NAME(e,epNR,eparam_names)
+#define INTMODIFIER(e) ENUM_NAME(e,eintmodNR,eintmod_names)
 #define EELTYPE(e)     ENUM_NAME(e,eelNR,eel_names)
 #define EVDWTYPE(e)    ENUM_NAME(e,evdwNR,evdw_names)
 #define ECONSTRTYPE(e) ENUM_NAME(e,econtNR,econstr_names)

diff --git a/include/types/enums.h b/include/types/enums.h
index ed01ad6e5f83ed33d117e033d7be5d057358df7d..9bc5948fbcfece19e1e466698e5c268ec872ecd5 100644 (file)
--- a/include/types/enums.h
+++ b/include/types/enums.h
@@ -80,6 +80,13 @@ enum {
   ecutsGROUP, ecutsVERLET, ecutsNR
 };
 
+/* Coulomb / VdW interaction modifiers.
+ * grompp replaces eintmodPOTSHIFT_VERLET by eintmodPOTSHIFT or eintmodNONE.
+ */
+enum {
+    eintmodPOTSHIFT_VERLET, eintmodPOTSHIFT, eintmodNONE, eintmodNR
+};
+
 /*
  * eelNOTUSED1 used to be GB, but to enable generalized born with different
  * forms of electrostatics (RF, switch, etc.) in the future it is now selected

diff --git a/include/types/forcerec.h b/include/types/forcerec.h
index 06210d2dd2a543ddc3f2a344b587b34df42b357f..c94e382a089478cd8e3e249f5b654322d03828f0 100644 (file)
--- a/include/types/forcerec.h
+++ b/include/types/forcerec.h
@@ -224,10 +224,12 @@ typedef struct {
   t_forcetable tab14; /* for 1-4 interactions only */
 
   /* PPPM & Shifting stuff */
+  gmx_bool coul_pot_shift;
   real rcoulomb_switch,rcoulomb;
   real *phi;
 
   /* VdW stuff */
+  gmx_bool vdw_pot_shift;
   double reppow;
   real rvdw_switch,rvdw;
   real bham_b_max;

diff --git a/include/types/inputrec.h b/include/types/inputrec.h
index 008b6a987b5aef33b33faf7a53a587c73dd0e365..b62fc949759cc3b1a22a9ae5e5b471e5c3f5dc07 100644 (file)
--- a/include/types/inputrec.h
+++ b/include/types/inputrec.h
@@ -317,6 +317,7 @@ typedef struct {
   real rlistlong;      /* long range pairlist cut-off (nm)             */
   real rtpi;            /* Radius for test particle insertion           */
   int  coulombtype;    /* Type of electrostatics treatment             */
+  int  coulomb_modifier; /* Modify the Coulomb interaction              */
   real rcoulomb_switch; /* Coulomb switch range start (nm)             */
   real rcoulomb;        /* Coulomb cutoff (nm)                         */
   real epsilon_r;       /* relative dielectric constant                 */ 
@@ -334,6 +335,7 @@ typedef struct {
   int  sa_algorithm;    /* Algorithm for SA part of GBSA                */
   real sa_surface_tension; /* Energy factor for SA part of GBSA */
   int  vdwtype;         /* Type of Van der Waals treatment              */
+  int  vdw_modifier;    /* Modify the VdW interaction                   */
   real rvdw_switch;     /* Van der Waals switch range start (nm)        */
   real rvdw;               /* Van der Waals cutoff (nm)                */
   int  eDispCorr;       /* Perform Long range dispersion corrections    */

diff --git a/share/html/online/mdp_opt.html b/share/html/online/mdp_opt.html
index e528c80b83f011c429531a8851fd1a75219a9b1f..ef2dd77c68e549048414ecf8b8cd6ee53aa84653 100644 (file)
--- a/share/html/online/mdp_opt.html
+++ b/share/html/online/mdp_opt.html
@@ -36,8 +36,8 @@ IF YOU'RE NOT SURE ABOUT WHAT YOU'RE DOING, DON'T DO IT!
 <li><a HREF="#tpi"><b>test particle insertion</b></a>(rtpi)
 <li><A HREF="#out"><b>output control</b></A> (nstxout, nstvout, nstfout, nstlog, nstcalcenergy, nstenergy, nstxtcout, xtc-precision, xtc-grps, energygrps)
 <li><A HREF="#nl"><b>neighbor searching</b></A> (cutoff-scheme, nstlist, ns-type, pbc, periodic-molecules, verlet-buffer-drift, rlist, rlistlong)
-<li><A HREF="#el"><b>electrostatics</b></A> (coulombtype, rcoulomb-switch, rcoulomb, epsilon-r, epsilon-rf)
-<li><A HREF="#vdw"><b>VdW</b></A> (vdwtype, rvdw-switch, rvdw, DispCorr)
+<li><A HREF="#el"><b>electrostatics</b></A> (coulombtype, coulomb-modifier, rcoulomb-switch, rcoulomb, epsilon-r, epsilon-rf)
+<li><A HREF="#vdw"><b>VdW</b></A> (vdwtype, vdw-modifier, rvdw-switch, rvdw, DispCorr)
 <li><A HREF="#table"><b>tables</b></A> (table-extension, energygrp-table)
 <li><A HREF="#ewald"><b>Ewald</b></A> (fourierspacing, fourier-nx, fourier-ny, fourier-nz, pme-order, ewald-rtol, ewald-geometry, epsilon-surface, optimize-fft)
 <li><A HREF="#tc"><b>Temperature coupling</b></A> (tcoupl, nsttcouple, tc-grps, tau-t, ref-t)
@@ -405,7 +405,7 @@ Native GPU acceleration is only supported with <b>Verlet</b>. With GPU-accelerat
 scaling <b>rcoulomb</b> and the grid spacing. This can be turned off with 
 <tt>-notunepme</tt>.
 
-<b>Verlet<\b> is somewhat faster than <b>group</b> when there is no water, or if <b>group</b> would use a pair-list buffer to conserve energy.
+<b>Verlet</b> is somewhat faster than <b>group</b> when there is no water, or if <b>group</b> would use a pair-list buffer to conserve energy.
 </dd>
 </dl></dd>
 
@@ -660,6 +660,20 @@ i.e. both to the user supplied function and the PME Mesh correction part.</dd>
 
 </dl></dd>
 
+<dt><b>coulomb-modifier:</b></dt>
+<dd><dl compact>
+<dt><b>Potential-shift-Verlet</b></dt>
+<dd>Selects <b>Potential-shift</b> with the Verlet cutoff-scheme,
+as it is (nearly) free; selects <b>None</b> with the group cutoff-scheme.</dd>
+<dt><b>Potential-shift</b></dt>
+<dd>Shift the Coulomb potential by a constant such that it is zero at the cut-off.
+This makes the potential the integral of the force. Note that this does not
+affect the forces or the sampling.</dd>
+<dt><b>None</b></dt>
+<dd>Use an unmodified Coulomb potential.</dd>
+</dl></dd>
+
+
 <A NAME="el2">
 <dt><b>rcoulomb-switch: (0) [nm]</b></dt>
 <dd>where to start switching the Coulomb potential</dd>
@@ -719,6 +733,19 @@ When <b>coulombtype</b> is not set to <b>User</b> the values
 for <tt>f</tt> and <tt>-f'</tt> are ignored.</dd>
 </dl></dd>
 
+<dt><b>vdw-modifier:</b></dt>
+<dd><dl compact>
+<dt><b>Potential-shift-Verlet</b></dt>
+<dd>Selects <b>Potential-shift</b> with the Verlet cutoff-scheme,
+as it is (nearly) free; selects <b>None</b> with the group cutoff-scheme.</dd>
+<dt><b>Potential-shift</b></dt>
+<dd>Shift the Van der Waals potential by a constant such that it is zero at the cut-off.
+This makes the potential the integral of the force. Note that this does not
+affect the forces or the sampling.</dd>
+<dt><b>None</b></dt>
+<dd>Use an unmodified Van der Waals potential.</dd>
+</dl></dd>
+
 <dt><b>rvdw-switch: (0) [nm]</b></dt>
 <dd>where to start switching the LJ potential</dd>
 
@@ -1906,6 +1933,7 @@ reals to your subroutine. Check the inputrec definition in
 <A HREF="#bond">constraints</A><br>
 <A HREF="#neq">cos-acceleration</A><br>
 <A HREF="#el">coulombtype</A><br>
+<A HREF="#el">coulomb-modifier</A><br>
 <A HREF="#free">couple-intramol</A><br>
 <A HREF="#free">couple-lambda0</A><br>
 <A HREF="#free">couple-lambda1</A><br>
@@ -2018,7 +2046,8 @@ reals to your subroutine. Check the inputrec definition in
 <A HREF="#user">userreal2</A><br>
 <A HREF="#user">userreal3</A><br>
 <A HREF="#user">userreal4</A><br>
-<A HREF="#el">vdwtype</A><br>
+<A HREF="#vdw">vdwtype</A><br>
+<A HREF="#vdw">vdw-modifier</A><br>
 <A HREF="#nl">verlet-buffer-drift</A><br>
 <A HREF="#out">xtc-grps</A><br>
 <A HREF="#out">xtc-precision</A><br>

diff --git a/src/gmxlib/names.c b/src/gmxlib/names.c
index f759537ba230d4c42ebee0e7d856d1f14207cd34..0d428e8f470564c1986806e720f83c76db0e15d9 100644 (file)
--- a/src/gmxlib/names.c
+++ b/src/gmxlib/names.c
@@ -94,6 +94,10 @@ const char *econstr_names[econtNR+1] = {
   "Lincs", "Shake", NULL
 };
 
+const char *eintmod_names[eintmodNR+1] = { 
+  "Potential-shift-Verlet","Potential-shift","None", NULL
+};
+
 const char *egrp_nm[egNR+1] = { 
   "Coul-SR","LJ-SR","Buck-SR", "Coul-LR", "LJ-LR", "Buck-LR",
   "Coul-14", "LJ-14", NULL

diff --git a/src/gmxlib/tpxio.c b/src/gmxlib/tpxio.c
index 4c5bc2c7ba5a7e8743ee621930ba8a96740de876..da659df491549c0dea599db5a482c2ba2d4668fe 100644 (file)
--- a/src/gmxlib/tpxio.c
+++ b/src/gmxlib/tpxio.c
@@ -73,7 +73,7 @@
 static const char *tpx_tag = TPX_TAG_RELEASE;
 
 /* This number should be increased whenever the file format changes! */
-static const int tpx_version = 80;
+static const int tpx_version = 81;
 
 /* This number should only be increased when you edit the TOPOLOGY section
  * or the HEADER of the tpx format.
@@ -690,9 +690,25 @@ static void do_inputrec(t_fileio *fio, t_inputrec *ir,gmx_bool bRead,
     gmx_fio_do_int(fio,ir->coulombtype); 
     if (file_version < 32 && ir->coulombtype == eelRF)
       ir->coulombtype = eelRF_NEC;      
+    if (file_version >= 81)
+    {
+        gmx_fio_do_int(fio,ir->coulomb_modifier); 
+    }
+    else
+    {
+        ir->coulomb_modifier = (ir->cutoff_scheme == ecutsVERLET ? eintmodPOTSHIFT : eintmodNONE);
+    }
     gmx_fio_do_real(fio,ir->rcoulomb_switch); 
     gmx_fio_do_real(fio,ir->rcoulomb); 
     gmx_fio_do_int(fio,ir->vdwtype);
+    if (file_version >= 81)
+    {
+        gmx_fio_do_int(fio,ir->vdw_modifier); 
+    }
+    else
+    {
+        ir->vdw_modifier = (ir->cutoff_scheme == ecutsVERLET ? eintmodPOTSHIFT : eintmodNONE);
+    }
     gmx_fio_do_real(fio,ir->rvdw_switch); 
     gmx_fio_do_real(fio,ir->rvdw); 
     if (file_version < 67) {

diff --git a/src/gmxlib/txtdump.c b/src/gmxlib/txtdump.c
index 25da6e1069c722c4164c3340343c45b42196b664..c5cf5bc32a393295713cc084c80399e73541052c 100644 (file)
--- a/src/gmxlib/txtdump.c
+++ b/src/gmxlib/txtdump.c
@@ -700,9 +700,11 @@ void pr_inputrec(FILE *fp,int indent,const char *title,t_inputrec *ir,
     PR("rlistlong",ir->rlistlong);
     PR("rtpi",ir->rtpi);
     PS("coulombtype",EELTYPE(ir->coulombtype));
+    PS("coulomb-modifier",INTMODIFIER(ir->coulomb_modifier));
     PR("rcoulomb-switch",ir->rcoulomb_switch);
     PR("rcoulomb",ir->rcoulomb);
     PS("vdwtype",EVDWTYPE(ir->vdwtype));
+    PS("vdw-modifier",INTMODIFIER(ir->vdw_modifier));
     PR("rvdw-switch",ir->rvdw_switch);
     PR("rvdw",ir->rvdw);
     if (ir->epsilon_r != 0)

diff --git a/src/kernel/readir.c b/src/kernel/readir.c
index 81d428bf982401d724589b5cfc174c5e13c98319..a6adbfa9ae6acd95db43675821e045a2547e95f4 100644 (file)
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@ -186,6 +186,21 @@ static int lcd(int n1,int n2)
   return d;
 }
 
+static void process_interaction_modifier(const t_inputrec *ir,int *eintmod)
+{
+    if (*eintmod == eintmodPOTSHIFT_VERLET)
+    {
+        if (ir->cutoff_scheme == ecutsVERLET)
+        {
+            *eintmod = eintmodPOTSHIFT;
+        }
+        else
+        {
+            *eintmod = eintmodNONE;
+        }
+    }
+}
+
 void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts,
               warninp_t wi)
 /* Check internal consistency */
@@ -221,8 +236,17 @@ void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts,
         warning_error(wi,"rlist should be >= 0");
     }
 
+    process_interaction_modifier(ir,&ir->coulomb_modifier);
+    process_interaction_modifier(ir,&ir->vdw_modifier);
+
     if (ir->cutoff_scheme == ecutsGROUP)
     {
+        if (ir->coulomb_modifier != eintmodNONE ||
+            ir->vdw_modifier != eintmodNONE)
+        {
+            warning_error(wi,"potential modifiers are not supported (yet) with the group cut-off scheme");
+        }
+
         /* BASIC CUT-OFF STUFF */
         if (ir->rlist == 0 ||
             !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
@@ -1562,6 +1586,7 @@ void get_ir(const char *mdparin,const char *mdparout,
   CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
   CTYPE ("Method for doing electrostatics");
   EETYPE("coulombtype",        ir->coulombtype,    eel_names);
+  EETYPE("coulomb-modifier",   ir->coulomb_modifier,    eintmod_names);
   CTYPE ("cut-off lengths");
   RTYPE ("rcoulomb-switch",    ir->rcoulomb_switch,    0.0);
   RTYPE ("rcoulomb",   ir->rcoulomb,   -1);
@@ -1570,6 +1595,7 @@ void get_ir(const char *mdparin,const char *mdparout,
   RTYPE ("epsilon-rf",  ir->epsilon_rf, 0.0);
   CTYPE ("Method for doing Van der Waals");
   EETYPE("vdw-type",   ir->vdwtype,    evdw_names);
+  EETYPE("vdw-modifier",       ir->vdw_modifier,    eintmod_names);
   CTYPE ("cut-off lengths");
   RTYPE ("rvdw-switch",        ir->rvdw_switch,        0.0);
   RTYPE ("rvdw",       ir->rvdw,       -1);

diff --git a/src/kernel/tpbcmp.c b/src/kernel/tpbcmp.c
index 6677703f0ff3d4acc214150c1bcc31c68a829428..4d6fca8cb757256c9bfe4a5af5d2f6d5da2cd764 100644 (file)
--- a/src/kernel/tpbcmp.c
+++ b/src/kernel/tpbcmp.c
@@ -634,10 +634,11 @@ static void cmp_inputrec(FILE *fp,t_inputrec *ir1,t_inputrec *ir2,real ftol, rea
   cmp_real(fp,"inputrec->rlistlong",-1,ir1->rlistlong,ir2->rlistlong,ftol,abstol);
   cmp_real(fp,"inputrec->rtpi",-1,ir1->rtpi,ir2->rtpi,ftol,abstol);
   cmp_int(fp,"inputrec->coulombtype",-1,ir1->coulombtype,ir2->coulombtype);
+  cmp_int(fp,"inputrec->coulomb_modifier",-1,ir1->coulomb_modifier,ir2->coulomb_modifier);
   cmp_real(fp,"inputrec->rcoulomb_switch",-1,ir1->rcoulomb_switch,ir2->rcoulomb_switch,ftol,abstol);
   cmp_real(fp,"inputrec->rcoulomb",-1,ir1->rcoulomb,ir2->rcoulomb,ftol,abstol);
   cmp_int(fp,"inputrec->vdwtype",-1,ir1->vdwtype,ir2->vdwtype);
-  cmp_real(fp,"inputrec->rvdw_switch",-1,ir1->rvdw_switch,ir2->rvdw_switch,ftol,abstol);
+  cmp_int(fp,"inputrec->vdw_modifier",-1,ir1->vdw_modifier,ir2->vdw_modifier);  cmp_real(fp,"inputrec->rvdw_switch",-1,ir1->rvdw_switch,ir2->rvdw_switch,ftol,abstol);
   cmp_real(fp,"inputrec->rvdw",-1,ir1->rvdw,ir2->rvdw,ftol,abstol);
   cmp_real(fp,"inputrec->epsilon_r",-1,ir1->epsilon_r,ir2->epsilon_r,ftol,abstol);
   cmp_real(fp,"inputrec->epsilon_rf",-1,ir1->epsilon_rf,ir2->epsilon_rf,ftol,abstol);

diff --git a/src/mdlib/forcerec.c b/src/mdlib/forcerec.c
index 82db34d70014688764f744d4181d2a2744b5e077..deee5d27927ca15fe15d2a4c148dbd0b42fbbd44 100644 (file)
--- a/src/mdlib/forcerec.c
+++ b/src/mdlib/forcerec.c
@@ -1574,10 +1574,6 @@ void init_interaction_const(FILE *fp,
                             const t_forcerec *fr)
 {
     interaction_const_t *ic;
-    gmx_bool shLJ,shCoul;
-
-    shLJ   = (getenv("GMX_NO_SHIFT_LJ") == NULL);
-    shCoul = (getenv("GMX_NO_SHIFT_COUL") == NULL);
 
     snew(ic, 1);
 
@@ -1585,7 +1581,7 @@ void init_interaction_const(FILE *fp,
 
     /* Lennard-Jones */
     ic->rvdw        = fr->rvdw;
-    if (shLJ)
+    if (fr->vdw_pot_shift)
     {
         ic->sh_invrc6 = pow(ic->rvdw,-6.0);
     }
@@ -1602,7 +1598,7 @@ void init_interaction_const(FILE *fp,
 
     /* Ewald */
     ic->ewaldcoeff  = fr->ewaldcoeff;
-    if (shCoul)
+    if (fr->coul_pot_shift)
     {
         ic->sh_ewald = gmx_erfc(ic->ewaldcoeff*ic->rcoulomb);
     }
@@ -1623,7 +1619,7 @@ void init_interaction_const(FILE *fp,
         /* For plain cut-off we might use the reaction-field kernels */
         ic->epsilon_rf = ic->epsilon_r;
         ic->k_rf       = 0;
-        if (shCoul)
+        if (fr->coul_pot_shift)
         {
             ic->c_rf   = 1/ic->rcoulomb;
         }
@@ -2011,6 +2007,9 @@ void init_forcerec(FILE *fp,
     fr->rlistlong  = cutoff_inf(ir->rlistlong);
     fr->eeltype    = ir->coulombtype;
     fr->vdwtype    = ir->vdwtype;
+
+    fr->coul_pot_shift = (ir->coulomb_modifier == eintmodPOTSHIFT);
+    fr->vdw_pot_shift  = (ir->vdw_modifier     == eintmodPOTSHIFT);
     
     fr->bTwinRange = fr->rlistlong > fr->rlist;
     fr->bEwald     = (EEL_PME(fr->eeltype) || fr->eeltype==eelEWALD);

diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c
index 92486e389b13ed068fee22f2d446adf332b36e60..a1c5473c0bf602eec0a35c19100c96d4ecddae43 100644 (file)
--- a/src/mdlib/sim_util.c
+++ b/src/mdlib/sim_util.c
@@ -2148,7 +2148,7 @@ void calc_enervirdiff(FILE *fplog,int eDispCorr,t_forcerec *fr)
       /* Contribution beyond the cut-off */
       eners[0] += -4.0*M_PI/(3.0*rc3);
       eners[1] +=  4.0*M_PI/(9.0*rc9);
-      if (fr->cutoff_scheme == ecutsVERLET && fr->ic->sh_invrc6 != 0) {
+      if (fr->vdw_pot_shift) {
           /* Contribution within the cut-off */
           eners[0] += -4.0*M_PI/(3.0*rc3);
           eners[1] +=  4.0*M_PI/(3.0*rc9);