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;
fr->enershiftsix = 0;
fr->enershifttwelve = 0;
eners[i] = 0;
virs[i] = 0;
}
- if (fr->vdwtype == evdwSWITCH || fr->vdwtype == evdwSHIFT ||
- fr->vdw_modifier == eintmodPOTSWITCH ||
- fr->vdw_modifier == eintmodFORCESWITCH)
+ if ((fr->vdw_modifier == eintmodPOTSHIFT) ||
+ (fr->vdw_modifier == eintmodPOTSWITCH) ||
+ (fr->vdw_modifier == eintmodFORCESWITCH) ||
+ (fr->vdwtype == evdwSHIFT) ||
+ (fr->vdwtype == evdwSWITCH))
{
- if (fr->rvdw_switch == 0)
+ if (((fr->vdw_modifier == eintmodPOTSWITCH) ||
+ (fr->vdw_modifier == eintmodFORCESWITCH) ||
+ (fr->vdwtype == evdwSWITCH)) && 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.
+ *
+ * For a pure potential-shift the potential has a constant shift
+ * all the way out to the cutoff, and that is it. For other forms
+ * we need to calculate the constant shift up to the point where we
+ * start modifying the potential.
+ */
+ ri0 = (fr->vdw_modifier == eintmodPOTSHIFT) ? ri1 : ri0;
+
r0 = ri0/scale;
r1 = ri1/scale;
rc3 = r0*r0*r0;
rc9 = rc3*rc3*rc3;
- if (fr->vdwtype == evdwSHIFT ||
- fr->vdw_modifier == eintmodFORCESWITCH)
+ if ((fr->vdw_modifier == eintmodFORCESWITCH) ||
+ (fr->vdwtype == evdwSHIFT))
{
/* Determine the constant energy shift below rvdw_switch.
* Table has a scale factor since we have scaled it down to compensate
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];
}
+ else if (fr->vdw_modifier == eintmodPOTSHIFT)
+ {
+ fr->enershiftsix = (real)(-1.0/(rc3*rc3));
+ fr->enershifttwelve = (real)( 1.0/(rc9*rc3));
+ }
+
/* 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.
*/
eners[0] += 4.0*M_PI*fr->enershiftsix*rc3/3.0;
eners[1] += 4.0*M_PI*fr->enershifttwelve*rc3/3.0;
+
+ /* Calculate the contribution in the range [r0,r1] where we
+ * modify the potential. For a pure potential-shift modifier we will
+ * have ri0==ri1, and there will not be any contribution here.
+ */
for (i = 0; i < 2; i++)
{
enersum = 0;
virs[i] -= virsum;
}
- /* now add the correction for rvdw_switch to infinity */
+ /* Alright: Above we compensated by REMOVING the parts outside r0
+ * corresponding to the ideal VdW 1/r6 and /r12 potentials.
+ *
+ * Regardless of whether r0 is the point where we start switching,
+ * or the cutoff where we calculated the constant shift, we include
+ * all the parts we are missing out to infinity from r0 by
+ * calculating the analytical dispersion correction.
+ */
eners[0] += -4.0*M_PI/(3.0*rc3);
eners[1] += 4.0*M_PI/(9.0*rc9);
virs[0] += 8.0*M_PI/rc3;
evdw_names[fr->vdwtype]);
}
- /* TODO: remove this code once we have group LJ-PME kernels
- * that calculate the exact, full LJ param C6/r^6 within the cut-off,
- * as the current nbnxn kernels do.
- */
+ /* When we deprecate the group kernels the code below can go too */
if (fr->vdwtype == evdwPME && fr->cutoff_scheme == ecutsGROUP)
{
/* Calculate self-interaction coefficient (assuming that
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff