GMX_LIBGMX_EXPORT
void
gmx_nonbonded_set_kernel_pointers(FILE * fplog,
- t_nblist * nl);
+ t_nblist * nl,
+ gmx_bool bElecAndVdwSwitchDiffers);
gmx_bool bDoForces;
real rcoulomb, rvdw, sh_invrc6;
gmx_bool bExactElecCutoff, bExactVdwCutoff;
- real rcutoff, rcutoff2, rswitch, d, d2, swV3, swV4, swV5, swF2, swF3, swF4, sw, dsw, rinvcorr;
- const real * tab_ewald_F;
- const real * tab_ewald_V;
- real tab_ewald_scale, tab_ewald_halfsp;
+ real d, d2, sw, dsw, rinvcorr;
+ real elec_swV3, elec_swV4, elec_swV5, elec_swF2, elec_swF3, elec_swF4;
+ real vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4;
+ gmx_bool bConvertEwaldToCoulomb, bComputeElecInteraction;
+ const real * ewtab;
+ int ewitab;
+ real ewrt, eweps, ewtabscale, ewtabhalfspace, sh_ewald;
+
+ sh_ewald = fr->ic->sh_ewald;
+ ewtab = fr->ic->tabq_coul_FDV0;
+ ewtabscale = fr->ic->tabq_scale;
+ ewtabhalfspace = 0.5/ewtabscale;
x = xx[0];
f = ff[0];
rvdw = fr->rvdw;
sh_invrc6 = fr->ic->sh_invrc6;
- /* Ewald (PME) reciprocal force and energy quadratic spline tables */
- tab_ewald_F = fr->ic->tabq_coul_F;
- tab_ewald_V = fr->ic->tabq_coul_V;
- tab_ewald_scale = fr->ic->tabq_scale;
- tab_ewald_halfsp = 0.5/tab_ewald_scale;
+ if (fr->coulomb_modifier == eintmodPOTSWITCH)
+ {
+ d = fr->rcoulomb-fr->rcoulomb_switch;
+ elec_swV3 = -10.0/(d*d*d);
+ elec_swV4 = 15.0/(d*d*d*d);
+ elec_swV5 = -6.0/(d*d*d*d*d);
+ elec_swF2 = -30.0/(d*d*d);
+ elec_swF3 = 60.0/(d*d*d*d);
+ elec_swF4 = -30.0/(d*d*d*d*d);
+ }
+ else
+ {
+ /* Avoid warnings from stupid compilers (looking at you, Clang!) */
+ elec_swV3 = elec_swV4 = elec_swV5 = elec_swF2 = elec_swF3 = elec_swF4 = 0.0;
+ }
- if (fr->coulomb_modifier == eintmodPOTSWITCH || fr->vdw_modifier == eintmodPOTSWITCH)
+ if (fr->vdw_modifier == eintmodPOTSWITCH)
{
- rcutoff = (fr->coulomb_modifier == eintmodPOTSWITCH) ? fr->rcoulomb : fr->rvdw;
- rcutoff2 = rcutoff*rcutoff;
- rswitch = (fr->coulomb_modifier == eintmodPOTSWITCH) ? fr->rcoulomb_switch : fr->rvdw_switch;
- d = rcutoff-rswitch;
- swV3 = -10.0/(d*d*d);
- swV4 = 15.0/(d*d*d*d);
- swV5 = -6.0/(d*d*d*d*d);
- swF2 = -30.0/(d*d*d);
- swF3 = 60.0/(d*d*d*d);
- swF4 = -30.0/(d*d*d*d*d);
+ d = fr->rvdw-fr->rvdw_switch;
+ vdw_swV3 = -10.0/(d*d*d);
+ vdw_swV4 = 15.0/(d*d*d*d);
+ vdw_swV5 = -6.0/(d*d*d*d*d);
+ vdw_swF2 = -30.0/(d*d*d);
+ vdw_swF3 = 60.0/(d*d*d*d);
+ vdw_swF4 = -30.0/(d*d*d*d*d);
}
else
{
- /* Stupid compilers dont realize these variables will not be used */
- rswitch = 0.0;
- swV3 = 0.0;
- swV4 = 0.0;
- swV5 = 0.0;
- swF2 = 0.0;
- swF3 = 0.0;
- swF4 = 0.0;
+ /* Avoid warnings from stupid compilers (looking at you, Clang!) */
+ vdw_swV3 = vdw_swV4 = vdw_swV5 = vdw_swF2 = vdw_swF3 = vdw_swF4 = 0.0;
}
bExactElecCutoff = (fr->coulomb_modifier != eintmodNONE) || fr->eeltype == eelRF_ZERO;
bExactVdwCutoff = (fr->vdw_modifier != eintmodNONE);
+ /* For Ewald/PME interactions we cannot easily apply the soft-core component to
+ * reciprocal space. When we use vanilla (not switch/shift) Ewald interactions, we
+ * can apply the small trick of subtracting the _reciprocal_ space contribution
+ * in this kernel, and instead apply the free energy interaction to the 1/r
+ * (standard coulomb) interaction.
+ *
+ * However, we cannot use this approach for switch-modified since we would then
+ * effectively end up evaluating a significantly different interaction here compared to the
+ * normal (non-free-energy) kernels, either by applying a cutoff at a different
+ * position than what the user requested, or by switching different
+ * things (1/r rather than short-range Ewald). For these settings, we just
+ * use the traditional short-range Ewald interaction in that case.
+ */
+ bConvertEwaldToCoulomb = ((icoul == GMX_NBKERNEL_ELEC_EWALD) && (fr->coulomb_modifier != eintmodPOTSWITCH));
+
/* fix compiler warnings */
nj1 = 0;
n1C = n1V = 0;
n1V = tab_elemsize*n0;
}
- /* With Ewald and soft-core we should put the cut-off on r,
- * not on the soft-cored rC, as the real-space and
- * reciprocal space contributions should (almost) cancel.
+
+ /* Only process the coulomb interactions if we have charges,
+ * and if we either include all entries in the list (no cutoff
+ * used in the kernel), or if we are within the cutoff.
*/
- if (qq[i] != 0 &&
- !(bExactElecCutoff &&
- ((icoul != GMX_NBKERNEL_ELEC_EWALD && rC >= rcoulomb) ||
- (icoul == GMX_NBKERNEL_ELEC_EWALD && r >= rcoulomb))))
+ bComputeElecInteraction = !bExactElecCutoff ||
+ ( bConvertEwaldToCoulomb && r < rcoulomb) ||
+ (!bConvertEwaldToCoulomb && rC < rcoulomb);
+
+ if ( (qq[i] != 0) && bComputeElecInteraction)
{
switch (icoul)
{
case GMX_NBKERNEL_ELEC_COULOMB:
- case GMX_NBKERNEL_ELEC_EWALD:
- /* simple cutoff (yes, ewald is done all on direct space for free energy) */
Vcoul[i] = qq[i]*rinvC;
FscalC[i] = Vcoul[i];
+ /* The shift for the Coulomb potential is stored in
+ * the RF parameter c_rf, which is 0 without shift
+ */
+ Vcoul[i] -= qq[i]*fr->ic->c_rf;
break;
case GMX_NBKERNEL_ELEC_REACTIONFIELD:
gmx_fatal(FARGS, "Free energy and GB not implemented.\n");
break;
+ case GMX_NBKERNEL_ELEC_EWALD:
+ if (bConvertEwaldToCoulomb)
+ {
+ Vcoul[i] = qq[i]*(rinvC-sh_ewald);
+ FscalC[i] = qq[i]*rinvC;
+ }
+ else
+ {
+ ewrt = rC*ewtabscale;
+ ewitab = (int) ewrt;
+ eweps = ewrt-ewitab;
+ ewitab = 4*ewitab;
+ FscalC[i] = ewtab[ewitab]+eweps*ewtab[ewitab+1];
+ rinvcorr = rinvC-sh_ewald;
+ Vcoul[i] = qq[i]*(rinvcorr-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+FscalC[i])));
+ FscalC[i] = qq[i]*(rinvC-rC*FscalC[i]);
+ }
+ break;
+
case GMX_NBKERNEL_ELEC_NONE:
FscalC[i] = 0.0;
Vcoul[i] = 0.0;
if (fr->coulomb_modifier == eintmodPOTSWITCH)
{
- d = rC-rswitch;
+ d = rC-fr->rcoulomb_switch;
d = (d > 0.0) ? d : 0.0;
d2 = d*d;
- sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
- dsw = d2*(swF2+d*(swF3+d*swF4));
+ sw = 1.0+d2*d*(elec_swV3+d*(elec_swV4+d*elec_swV5));
+ dsw = d2*(elec_swF2+d*(elec_swF3+d*elec_swF4));
+ FscalC[i] = FscalC[i]*sw - rC*Vcoul[i]*dsw;
Vcoul[i] *= sw;
- FscalC[i] = FscalC[i]*sw + Vcoul[i]*dsw;
+
+ FscalC[i] = (rC < rcoulomb) ? FscalC[i] : 0.0;
+ Vcoul[i] = (rC < rcoulomb) ? Vcoul[i] : 0.0;
}
}
- if ((c6[i] != 0 || c12[i] != 0) &&
- !(bExactVdwCutoff && rV >= rvdw))
+ if ((c6[i] != 0 || c12[i] != 0) && ( !bExactVdwCutoff || rV < rvdw ) )
{
switch (ivdw)
{
if (fr->vdw_modifier == eintmodPOTSWITCH)
{
- d = rV-rswitch;
+ d = rV-fr->rvdw_switch;
d = (d > 0.0) ? d : 0.0;
d2 = d*d;
- sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
- dsw = d2*(swF2+d*(swF3+d*swF4));
+ sw = 1.0+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5));
+ dsw = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4));
+ FscalV[i] = FscalV[i]*sw - rV*Vvdw[i]*dsw;
Vvdw[i] *= sw;
- FscalV[i] = FscalV[i]*sw + Vvdw[i]*dsw;
FscalV[i] = (rV < rvdw) ? FscalV[i] : 0.0;
Vvdw[i] = (rV < rvdw) ? Vvdw[i] : 0.0;
Fscal = 0;
- if (icoul == GMX_NBKERNEL_ELEC_EWALD &&
- !(bExactElecCutoff && r >= rcoulomb))
+ if (bConvertEwaldToCoulomb && ( !bExactElecCutoff || r < rcoulomb ) )
{
- /* Because we compute the soft-core normally,
- * we have to remove the Ewald short range portion.
- * Done outside of the states loop because this part
- * doesn't depend on the scaled R.
+ /* See comment in the preamble. When using Ewald interactions
+ * (unless we use a switch modifier) we subtract the reciprocal-space
+ * Ewald component here which made it possible to apply the free
+ * energy interaction to 1/r (vanilla coulomb short-range part)
+ * above. This gets us closer to the ideal case of applying
+ * the softcore to the entire electrostatic interaction,
+ * including the reciprocal-space component.
*/
- real rs, frac, f_lr;
- int ri;
+ real v_lr, f_lr;
- rs = rsq*rinv*tab_ewald_scale;
- ri = (int)rs;
- frac = rs - ri;
- f_lr = (1 - frac)*tab_ewald_F[ri] + frac*tab_ewald_F[ri+1];
- FF = f_lr*rinv;
- VV = tab_ewald_V[ri] - tab_ewald_halfsp*frac*(tab_ewald_F[ri] + f_lr);
+ ewrt = r*ewtabscale;
+ ewitab = (int) ewrt;
+ eweps = ewrt-ewitab;
+ ewitab = 4*ewitab;
+ f_lr = ewtab[ewitab]+eweps*ewtab[ewitab+1];
+ v_lr = (ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+f_lr));
+ f_lr *= rinv;
for (i = 0; i < NSTATES; i++)
{
- vctot -= LFC[i]*qq[i]*VV;
- Fscal -= LFC[i]*qq[i]*FF;
- dvdl_coul -= (DLF[i]*qq[i])*VV;
+ vctot -= LFC[i]*qq[i]*v_lr;
+ Fscal -= LFC[i]*qq[i]*f_lr;
+ dvdl_coul -= (DLF[i]*qq[i])*v_lr;
}
}
bExactElecCutoff = (fr->coulomb_modifier != eintmodNONE) || fr->eeltype == eelRF_ZERO;
bExactVdwCutoff = (fr->vdw_modifier != eintmodNONE);
- bExactCutoff = bExactElecCutoff || bExactVdwCutoff;
+ bExactCutoff = bExactElecCutoff && bExactVdwCutoff;
if (bExactCutoff)
{
velec = 0;
vvdw = 0;
- if (bExactCutoff && rsq > rcutoff2)
+ if (bExactCutoff && rsq >= rcutoff2)
{
continue;
}
/* Vanilla cutoff coulomb */
velec = qq*rinv;
felec = velec*rinvsq;
+ /* The shift for the Coulomb potential is stored in
+ * the RF parameter c_rf, which is 0 without shift
+ */
+ velec -= qq*fr->ic->c_rf;
break;
case GMX_NBKERNEL_ELEC_REACTIONFIELD:
}
if (bExactElecCutoff)
{
- felec = (rsq <= rcoulomb2) ? felec : 0.0;
- velec = (rsq <= rcoulomb2) ? velec : 0.0;
+ felec = (rsq < rcoulomb2) ? felec : 0.0;
+ velec = (rsq < rcoulomb2) ? velec : 0.0;
}
vctot += velec;
} /* End of coulomb interactions */
}
if (bExactVdwCutoff)
{
- fvdw = (rsq <= rvdw2) ? fvdw : 0.0;
- vvdw = (rsq <= rvdw2) ? vvdw : 0.0;
+ fvdw = (rsq < rvdw2) ? fvdw : 0.0;
+ vvdw = (rsq < rvdw2) ? vvdw : 0.0;
}
vvdwtot += vvdw;
} /* end VdW interactions */
nb_kernel_list_add_kernels(kernellist_avx_256_double, kernellist_avx_256_double_size);
#endif
#if (defined GMX_CPU_ACCELERATION_SPARC64_HPC_ACE && defined GMX_DOUBLE)
- nb_kernel_list_add_kernels(kernellist_sparc64_hpc_ace_double,kernellist_sparc64_hpc_ace_double_size);
+ nb_kernel_list_add_kernels(kernellist_sparc64_hpc_ace_double, kernellist_sparc64_hpc_ace_double_size);
#endif
; /* empty statement to avoid a completely empty block */
}
void
-gmx_nonbonded_set_kernel_pointers(FILE *log, t_nblist *nl)
+gmx_nonbonded_set_kernel_pointers(FILE *log, t_nblist *nl, gmx_bool bElecAndVdwSwitchDiffers)
{
const char * elec;
const char * elec_mod;
}
}
- /* Give up. If this was a water kernel, leave the pointer as NULL, which
- * will disable water optimization in NS. If it is a particle kernel, set
- * the pointer to the generic NB kernel.
+ /* For now, the accelerated kernels cannot handle the combination of switch functions for both
+ * electrostatics and VdW that use different switch radius or switch cutoff distances
+ * (both of them enter in the switch function calculation). This would require
+ * us to evaluate two completely separate switch functions for every interaction.
+ * Instead, we disable such kernels by setting the pointer to NULL.
+ * This will cause the generic kernel (which can handle it) to be called instead.
+ *
+ * Note that we typically already enable tabulated coulomb interactions for this case,
+ * so this is mostly a safe-guard to make sure we call the generic kernel if the
+ * tables are disabled.
+ */
+ if ((nl->ielec != GMX_NBKERNEL_ELEC_NONE) && (nl->ielecmod == eintmodPOTSWITCH) &&
+ (nl->ivdw != GMX_NBKERNEL_VDW_NONE) && (nl->ivdwmod == eintmodPOTSWITCH) &&
+ bElecAndVdwSwitchDiffers)
+ {
+ nl->kernelptr_vf = NULL;
+ nl->kernelptr_f = NULL;
+ }
+
+ /* Give up, pick a generic one instead.
+ * We only do this for particle-particle kernels; by leaving the water-optimized kernel
+ * pointers to NULL, the water optimization will automatically be disabled for this interaction.
*/
if (nl->kernelptr_vf == NULL && !gmx_strcasecmp_min(geom,"Particle-Particle"))
{
fprintf(debug,
"WARNING - Slow generic NB kernel used for neighborlist with\n"
" Elec: '%s', Modifier: '%s'\n"
- " Vdw: '%s', Modifier: '%s'\n"
- " Geom: '%s', Other: '%s'\n\n",
- elec, elec_mod, vdw, vdw_mod, geom, other);
+ " Vdw: '%s', Modifier: '%s'\n",
+ elec, elec_mod, vdw, vdw_mod);
}
}
}
-
return;
}
/* We don't need the non-perturbed interactions */
continue;
}
- /* Neighborlists whose kernelptr==NULL will always be empty */
+
if(kernelptr != NULL)
{
(*kernelptr)(&(nlist[i]), x, f, fr, mdatoms, &kernel_data, nrnb);
}
+ else
+ {
+ gmx_fatal(FARGS, "Non-empty neighborlist does not have any kernel pointer assigned.");
+ }
}
}
}
if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
{
real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
- sprintf(warn_buf, "The switching range for should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
+ sprintf(warn_buf, "The switching range should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
percentage,ir->rcoulomb_switch,ir->rcoulomb,ir->ewald_rtol);
warning(wi, warn_buf);
}
fr->nbkernel_elec_modifier = fr->coulomb_modifier;
fr->nbkernel_vdw_modifier = fr->vdw_modifier;
+ fr->rvdw = cutoff_inf(ir->rvdw);
+ fr->rvdw_switch = ir->rvdw_switch;
+ fr->rcoulomb = cutoff_inf(ir->rcoulomb);
+ fr->rcoulomb_switch = ir->rcoulomb_switch;
+
fr->bTwinRange = fr->rlistlong > fr->rlist;
fr->bEwald = (EEL_PME(fr->eeltype) || fr->eeltype == eelEWALD);
/* If the user absolutely wants different switch/shift settings for coul/vdw, it is likely
* going to be faster to tabulate the interaction than calling the generic kernel.
+ * However, if generic kernels have been requested we keep things analytically.
*/
- if (fr->nbkernel_elec_modifier == eintmodPOTSWITCH && fr->nbkernel_vdw_modifier == eintmodPOTSWITCH)
+ if (fr->nbkernel_elec_modifier == eintmodPOTSWITCH &&
+ fr->nbkernel_vdw_modifier == eintmodPOTSWITCH &&
+ bGenericKernelOnly == FALSE)
{
if ((fr->rcoulomb_switch != fr->rvdw_switch) || (fr->rcoulomb != fr->rvdw))
{
fr->bcoultab = TRUE;
+ /* Once we tabulate electrostatics, we can use the switch function for LJ,
+ * which would otherwise need two tables.
+ */
}
}
else if ((fr->nbkernel_elec_modifier == eintmodPOTSHIFT && fr->nbkernel_vdw_modifier == eintmodPOTSHIFT) ||
fr->nbkernel_elec_modifier == eintmodEXACTCUTOFF &&
(fr->nbkernel_vdw_modifier == eintmodPOTSWITCH || fr->nbkernel_vdw_modifier == eintmodPOTSHIFT))))
{
- if (fr->rcoulomb != fr->rvdw)
+ if ((fr->rcoulomb != fr->rvdw) && (bGenericKernelOnly == FALSE))
{
fr->bcoultab = TRUE;
}
fr->epsilon_r = ir->epsilon_r;
fr->epsilon_rf = ir->epsilon_rf;
fr->fudgeQQ = mtop->ffparams.fudgeQQ;
- fr->rcoulomb_switch = ir->rcoulomb_switch;
- fr->rcoulomb = cutoff_inf(ir->rcoulomb);
/* Parameters for generalized RF */
fr->zsquare = 0.0;
fr->egp_flags = ir->opts.egp_flags;
/* Van der Waals stuff */
- fr->rvdw = cutoff_inf(ir->rvdw);
- fr->rvdw_switch = ir->rvdw_switch;
if ((fr->vdwtype != evdwCUT) && (fr->vdwtype != evdwUSER) && !fr->bBHAM)
{
if (fr->rvdw_switch >= fr->rvdw)
bTab = fr->bcoultab || fr->bvdwtab || fr->bEwald;
bSep14tab = ((!bTab || fr->eeltype != eelCUT || fr->vdwtype != evdwCUT ||
+ fr->coulomb_modifier != eintmodNONE ||
+ fr->vdw_modifier != eintmodNONE ||
fr->bBHAM || fr->bEwald) &&
(gmx_mtop_ftype_count(mtop, F_LJ14) > 0 ||
gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0 ||
}
}
}
+ else if ((fr->eDispCorr != edispcNO) &&
+ ((fr->vdw_modifier == eintmodPOTSWITCH) || (fr->vdw_modifier == eintmodPOTSHIFT)))
+ {
+ /* Tables might not be used for the potential modifier interactions per se, but
+ * we still need them to evaluate switch/shift dispersion corrections in this case.
+ */
+ make_nbf_tables(fp, oenv, fr, rtab, cr, tabfn, NULL, NULL, &fr->nblists[0]);
+ }
+
if (bSep14tab)
{
/* generate extra tables with plain Coulomb for 1-4 interactions only */
int maxsr, int maxlr,
int ivdw, int ivdwmod,
int ielec, int ielecmod,
- int igeometry, int type)
+ int igeometry, int type,
+ gmx_bool bElecAndVdwSwitchDiffers)
{
t_nblist *nl;
int homenr;
}
/* This will also set the simd_padding_width field */
- gmx_nonbonded_set_kernel_pointers( (i == 0) ? log : NULL, nl);
+ gmx_nonbonded_set_kernel_pointers( (i == 0) ? log : NULL, nl, bElecAndVdwSwitchDiffers);
/* maxnri is influenced by the number of shifts (maximum is 8)
* and the number of energy groups.
* cache trashing.
*/
int maxsr, maxsr_wat, maxlr, maxlr_wat;
- int ielec, ielecf, ivdw, ielecmod, ielecmodf, ivdwmod, type;
+ int ielec, ivdw, ielecmod, ivdwmod, type;
int solvent;
int igeometry_def, igeometry_w, igeometry_ww;
int i;
+ gmx_bool bElecAndVdwSwitchDiffers;
t_nblists *nbl;
/* maxsr = homenr-fr->nWatMol*3; */
}
/* Determine the values for ielec/ivdw. */
- ielec = fr->nbkernel_elec_interaction;
- ivdw = fr->nbkernel_vdw_interaction;
- ielecmod = fr->nbkernel_elec_modifier;
- ivdwmod = fr->nbkernel_vdw_modifier;
- type = GMX_NBLIST_INTERACTION_STANDARD;
+ ielec = fr->nbkernel_elec_interaction;
+ ivdw = fr->nbkernel_vdw_interaction;
+ ielecmod = fr->nbkernel_elec_modifier;
+ ivdwmod = fr->nbkernel_vdw_modifier;
+ type = GMX_NBLIST_INTERACTION_STANDARD;
+ bElecAndVdwSwitchDiffers = ( (fr->rcoulomb_switch != fr->rvdw_switch) || (fr->rcoulomb != fr->rvdw));
fr->ns.bCGlist = (getenv("GMX_NBLISTCG") != 0);
if (!fr->ns.bCGlist)
type = GMX_NBLIST_INTERACTION_ADRESS;
}
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ], &nbl->nlist_lr[eNL_VDWQQ],
- maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type);
+ maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_VDW], &nbl->nlist_lr[eNL_VDW],
- maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type);
+ maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_QQ], &nbl->nlist_lr[eNL_QQ],
- maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type);
+ maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATER], &nbl->nlist_lr[eNL_VDWQQ_WATER],
- maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type);
+ maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_QQ_WATER], &nbl->nlist_lr[eNL_QQ_WATER],
- maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_w, type);
+ maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATERWATER], &nbl->nlist_lr[eNL_VDWQQ_WATERWATER],
- maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type);
+ maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_QQ_WATERWATER], &nbl->nlist_lr[eNL_QQ_WATERWATER],
- maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_ww, type);
+ maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
/* Did we get the solvent loops so we can use optimized water kernels? */
if (nbl->nlist_sr[eNL_VDWQQ_WATER].kernelptr_vf == NULL
if (fr->efep != efepNO)
{
- if ((fr->bEwald) && (fr->sc_alphacoul > 0)) /* need to handle long range differently if using softcore */
- {
- ielecf = GMX_NBKERNEL_ELEC_EWALD;
- ielecmodf = eintmodNONE;
- }
- else
- {
- ielecf = ielec;
- ielecmodf = ielecmod;
- }
-
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_FREE], &nbl->nlist_lr[eNL_VDWQQ_FREE],
- maxsr, maxlr, ivdw, ivdwmod, ielecf, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+ maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_VDW_FREE], &nbl->nlist_lr[eNL_VDW_FREE],
- maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+ maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_QQ_FREE], &nbl->nlist_lr[eNL_QQ_FREE],
- maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielecf, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+ maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
}
}
/* QMMM MM list */
if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
{
init_nblist(log, &fr->QMMMlist, NULL,
- maxsr, maxlr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD);
+ maxsr, maxlr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD, bElecAndVdwSwitchDiffers);
}
if (log != NULL)
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;
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++)
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
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.
sfree(td->f);
}
-static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr)
+static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
+ gmx_bool b14only)
{
/* Fill the table according to the formulas in the manual.
* In principle, we only need the potential and the second
int i;
double reppow, p;
double r1, rc, r12, r13;
- double r, r2, r6, rc6;
+ double r, r2, r6, rc2, rc6, rc12;
double expr, Vtab, Ftab;
/* Parameters for David's function */
double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
/* Parameters for the switching function */
double ksw, swi, swi1;
/* Temporary parameters */
- gmx_bool bSwitch, bShift;
+ gmx_bool bSwitch, bForceShift, bPotentialShift;
double ewc = fr->ewaldcoeff;
+ double Vcut = 0;
- bSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
- (tp == etabCOULSwitch) ||
- (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch));
- bShift = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
- (tp == etabShift));
+ if(b14only)
+ {
+ bSwitch = FALSE;
+ bForceShift = FALSE;
+ bPotentialShift = FALSE;
+ }
+ else
+ {
+ bSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
+ (tp == etabCOULSwitch) ||
+ (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSWITCH)));
+ bForceShift = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
+ (tp == etabShift));
+ bPotentialShift = ((tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSHIFT)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSHIFT)));
+ }
reppow = fr->reppow;
{
ksw = 0.0;
}
- if (bShift)
+ if (bForceShift)
{
if (tp == etabShift)
{
fp = xvgropen("switch.xvg", "switch", "r", "s");
#endif
+ if(bPotentialShift)
+ {
+ rc2 = rc*rc;
+ rc6 = 1.0/(rc2*rc2*rc2);
+ if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
+ {
+ rc12 = rc6*rc6;
+ }
+ else
+ {
+ rc12 = pow(rc, -reppow);
+ }
+
+ switch (tp)
+ {
+ case etabLJ6:
+ /* Dispersion */
+ Vcut = -rc6;
+ break;
+ case etabLJ12:
+ /* Repulsion */
+ Vcut = rc12;
+ break;
+ case etabCOUL:
+ Vcut = 1.0/rc;
+ break;
+ case etabEwald:
+ case etabEwaldSwitch:
+ Vtab = gmx_erfc(ewc*rc)/rc;
+ break;
+ case etabEwaldUser:
+ /* Only calculate minus the reciprocal space contribution */
+ Vtab = -gmx_erf(ewc*rc)/rc;
+ break;
+ case etabRF:
+ case etabRF_ZERO:
+ /* No need for preventing the usage of modifiers with RF */
+ Vcut = 0.0;
+ break;
+ case etabEXPMIN:
+ Vcut = exp(-rc);
+ break;
+ default:
+ gmx_fatal(FARGS, "Cannot apply new potential-shift modifier to interaction type '%s' yet. (%s,%d)",
+ tprops[tp].name, __FILE__, __LINE__);
+ }
+ }
+
for (i = td->nx0; (i < td->nx); i++)
{
r = td->x[i];
gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)",
tp, __FILE__, __LINE__);
}
- if (bShift)
+ if (bForceShift)
{
/* Normal coulomb with cut-off correction for potential */
if (r < rc)
}
}
}
+ if (bPotentialShift)
+ {
+ Vtab -= Vcut;
+ }
if (ETAB_USER(tp))
{
init_table(out, nx, nx0,
(tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale,
&(td[k]), !bReadTab);
- fill_table(&(td[k]), tabsel[k], fr);
+ fill_table(&(td[k]), tabsel[k], fr, b14only);
if (out)
{
fprintf(out, "%s table with %d data points for %s%s.\n"
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff