Fixed shift and switch modifiers, particularly for free-energy

[alexxy/gromacs.git] / src / mdlib / sim_util.c

diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c
index 26c709d6feec7f066328d544fd8577582da303e7..59868db7e85ccb4a993cc5e67ad2b184a22c8188 100644 (file)
--- a/src/mdlib/sim_util.c
+++ b/src/mdlib/sim_util.c
@@ -2139,11 +2139,12 @@ void do_constrain_first(FILE *fplog, gmx_constr_t constr,
 
 void calc_enervirdiff(FILE *fplog, int eDispCorr, t_forcerec *fr)
 {
-    double eners[2], virs[2], enersum, virsum, y0, f, g, h;
-    double r0, r1, r, rc3, rc9, ea, eb, ec, pa, pb, pc, pd;
-    double invscale, invscale2, invscale3;
-    int    ri0, ri1, ri, i, offstart, offset;
-    real   scale, *vdwtab, tabfactor, tmp;
+    double   eners[2], virs[2], enersum, virsum, y0, f, g, h;
+    double   r0, r1, r, rc3, rc9, ea, eb, ec, pa, pb, pc, pd;
+    double   invscale, invscale2, invscale3;
+    int      ri0, ri1, ri, i, offstart, offset;
+    real     scale, *vdwtab, tabfactor, tmp;
+    gmx_bool bSwitch, bShift;
 
     fr->enershiftsix    = 0;
     fr->enershifttwelve = 0;
@@ -2152,6 +2153,9 @@ void calc_enervirdiff(FILE *fplog, int eDispCorr, t_forcerec *fr)
     fr->virdiffsix      = 0;
     fr->virdifftwelve   = 0;
 
+    bSwitch = (fr->vdwtype == evdwSWITCH) || (fr->vdw_modifier == eintmodPOTSWITCH);
+    bShift  = (fr->vdwtype == evdwSHIFT)  || (fr->vdw_modifier == eintmodPOTSHIFT);
+
     if (eDispCorr != edispcNO)
     {
         for (i = 0; i < 2; i++)
@@ -2159,27 +2163,41 @@ void calc_enervirdiff(FILE *fplog, int eDispCorr, t_forcerec *fr)
             eners[i] = 0;
             virs[i]  = 0;
         }
-        if ((fr->vdwtype == evdwSWITCH) || (fr->vdwtype == evdwSHIFT))
+        if (bSwitch || bShift)
         {
-            if (fr->rvdw_switch == 0)
+            if (bSwitch && fr->rvdw_switch == 0)
             {
                 gmx_fatal(FARGS,
                           "With dispersion correction rvdw-switch can not be zero "
                           "for vdw-type = %s", evdw_names[fr->vdwtype]);
             }
 
-            scale  = fr->nblists[0].table_elec_vdw.scale;
+            scale  = fr->nblists[0].table_vdw.scale;
             vdwtab = fr->nblists[0].table_vdw.data;
 
             /* Round the cut-offs to exact table values for precision */
             ri0  = floor(fr->rvdw_switch*scale);
             ri1  = ceil(fr->rvdw*scale);
+
+            /* The code below has some support for handling force-switching, i.e.
+             * when the force (instead of potential) is switched over a limited
+             * region. This leads to a constant shift in the potential inside the
+             * switching region, which we can handle by adding a constant energy
+             * term in the force-switch case just like when we do potential-shift.
+             *
+             * For now this is not enabled, but to keep the functionality in the
+             * code we check separately for switch and shift. When we do force-switch
+             * the shifting point is rvdw_switch, while it is the cutoff when we
+             * have a classical potential-shift.
+             */
+            ri0  = (bShift) ? ri1 : ri0;
+
             r0   = ri0/scale;
             r1   = ri1/scale;
             rc3  = r0*r0*r0;
             rc9  = rc3*rc3*rc3;
 
-            if (fr->vdwtype == evdwSHIFT)
+            if (bShift)
             {
                 /* Determine the constant energy shift below rvdw_switch.
                  * Table has a scale factor since we have scaled it down to compensate
@@ -2188,6 +2206,7 @@ void calc_enervirdiff(FILE *fplog, int eDispCorr, t_forcerec *fr)
                 fr->enershiftsix    = (real)(-1.0/(rc3*rc3)) - 6.0*vdwtab[8*ri0];
                 fr->enershifttwelve = (real)( 1.0/(rc9*rc3)) - 12.0*vdwtab[8*ri0 + 4];
             }
+
             /* Add the constant part from 0 to rvdw_switch.
              * This integration from 0 to rvdw_switch overcounts the number
              * of interactions by 1, as it also counts the self interaction.